Hello, I need help in solving 3 errors in my C++ code
The first error is** line 95 'output_minus_10' was not declared in the scope**
The Second error is line 96 'output_plus_10' was not declared in the scope
The warning is line 104 unused variable 'R_junction'
#include <cstdio>
#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#define NAME 2000 //define size of output file
#define DT 0.01 //units is mS
using namespace std;
//Define global variables
double V_H ;
double V_L ;
double gamma_H ;
double gamma_L ;
double V_H2 ;
double V_L2;
double gamma_H2;
double gamma_L2;
double alpha_coef;
double beta_coef;
double V_alpha;
double V_beta ;
double alpha_coef2;
double beta_coef2;
double V_alpha2 ;
double V_beta2 ;
int Nchans;
double n_LL0 ;
double n_LH0;
double n_HL0;
double n_HH0;
//Filename to select model parameters
const char *params_file = "cx43cx45_model.dat";
//Function to read files and retune a vector of content
vector<double> fileToVector(const char *name);
int main()
{
//Read the parameters from file and store in vector modelparams
vector<double> modelparams = fileToVector(params_file);
vector<double> modelparams_plus_10 = fileToVector(params_file);
vector<double> modelparams_minus_10 = fileToVector(params_file);
//Assign parameter values from file and assign to global model variables
//Increase at 10% increments
//Exclude (end 5)
V_H = modelparams[0];
V_L = modelparams[1];
gamma_H = modelparams[2];
gamma_L = modelparams[3];
alpha_coef = modelparams[4]*1e-3;//units adjusted to ms
beta_coef = modelparams[5]*1e-3;//units adjusted to ms
V_alpha = modelparams[6];
V_beta = modelparams[7];
V_H2 = modelparams[8];
V_L2 = modelparams[9];
gamma_H2 = modelparams[10];
gamma_L2 = modelparams[11];
alpha_coef2 = modelparams[12]*1e-3;//units adjusted to ms
beta_coef2 = modelparams[13]*1e-3;//units adjusted to ms
V_alpha2 = modelparams[14];
V_beta2 = modelparams[15];
Nchans = (int) (modelparams[16]);
n_LL0 = modelparams[17];
n_LH0 = modelparams[18];
n_HL0 = modelparams[19];
n_HH0 = 1-(n_LL0+n_LH0+n_HL0+n_HH0);
for (int i = 0; i < 16; i++)
{
modelparams_plus_10[i] = modelparams[i]*1.1;
modelparams_minus_10[i] = modelparams[i]*0.9;
}
//Define output file (C style code) to store model data
FILE *output;
char output_file_name[ NAME ];
printf( "\nEnter an output file name: " );
scanf( "%s", output_file_name );
bool print_header = 0;//Flag to print file header
output = fopen( output_file_name, "w" );//Open file
output_plus_10 = fopen( "cx45_model_VH_2.dat", "w" );//Open file
output_minus_10 = fopen( "cx45_model_VH_1.dat", "w" );//Open file
//Declare variables for calculating state conductance
double alpha_1, alpha_2, alpha_3, alpha_4, beta_1, beta_2, beta_3, beta_4;
double V_LL1, V_LL2, V_LH1, V_LH2, V_HL1, V_HL2, V_HH1, V_HH2, Vj;
double g_1, g_2, delta_g_1, delta_g_2, R_1, R_2, A, B, C, D, determinant, V_1, V_2, gamma_1, gamma_2;
double g_HH, g_LH, g_HL, g_LL;
double n_HH,n_LH,n_HL,n_LL;
double g_junction, R_junction, I_junction;
double tolerance = 0.001;
double VJ_array[] = {-130,-100,-70,-40,-10,0,10,40,70,100,130};//Array containing voltage step protocol in units of mV
//Voltage clamp loop (steps through voltages in VJ_array
for (int k=0;k< 11; k++)
{
double t=0;
int count = 0;
Vj = VJ_array[k];
// Procedure for calculating g_LL
g_1 = 10.0; // Initial guess for g_1
g_2 = 10.0; // Initial guess for g_2
V_1 = V_L; // State of hemichannel 1
V_2 = V_L2; // State of hemichannel 2
gamma_1 = gamma_L;
gamma_2 = gamma_L2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1; // Residual equation 1
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2; // Residual equation 2
while ( fabs(R_1) > tolerance || fabs(R_2) > tolerance )
{
A = -((gamma_1*Vj*g_2/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2))))-1.0; // Partial derivative of equation 1 wrt g_1
B = (gamma_1*Vj*g_1/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2))); // Partial derivative of equation 1 wrt g_2
C = -(gamma_2*Vj*g_2/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2))); // Partial derivative of equation 2 wrt g_1
D = ((gamma_2*Vj*g_1/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2))))-1.0; // Partial derivative of equation 2 wrt g_2
determinant = A*D-B*C; // Determinant
delta_g_1 = (1.0/determinant)*(D*(-R_1)-B*(-R_2)); // Change in g_1
delta_g_2 = (1.0/determinant)*(-C*(-R_1)+A*(-R_2)); // Change in g_2
g_1 = g_1+delta_g_1; // Updated value of g_1
g_2 = g_2+delta_g_2; // Updated value of g_2
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1; // Residual for equation 1
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2; // Residual for equation 2
}
g_LL = (g_1*g_2)/(g_1+g_2); // Junction conductance in the LL state
V_LL1 = -Vj*(g_2/(g_1+g_2)); // Voltage drop across the first hemichannel
V_LL2 = Vj*(g_1/(g_1+g_2)); // Voltage drop across the second hemichannel
// Procedure for calculating g_LH
g_1 = 10.0;
g_2 = 10.0;
V_1 = V_L;
V_2 = V_H2;
gamma_1 = gamma_L;
gamma_2 = gamma_H2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
while ( fabs(R_1) > tolerance || fabs(R_2) > tolerance )
{
A = -((gamma_1*Vj*g_2/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2))))-1.0;
B = (gamma_1*Vj*g_1/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2)));
C = -(gamma_2*Vj*g_2/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2)));
D = ((gamma_2*Vj*g_1/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2))))-1.0;
determinant = A*D-B*C;
delta_g_1 = (1.0/determinant)*(D*(-R_1)-B*(-R_2));
delta_g_2 = (1.0/determinant)*(-C*(-R_1)+A*(-R_2));
g_1 = g_1+delta_g_1;
g_2 = g_2+delta_g_2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
}
g_LH = (g_1*g_2)/(g_1+g_2);
V_LH1 = -Vj*(g_2/(g_1+g_2));
V_LH2 = Vj*(g_1/(g_1+g_2));
// Procedure for calculating g_HL
g_1 = 10.0;
g_2 = 10.0;
V_1 = V_H;
V_2 = V_L2;
gamma_1 = gamma_H;
gamma_2 = gamma_L2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
while ( fabs(R_1) > tolerance || fabs(R_2) > tolerance )
{
A = -((gamma_1*Vj*g_2/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2))))-1.0;
B = (gamma_1*Vj*g_1/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2)));
C = -(gamma_2*Vj*g_2/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2)));
D = ((gamma_2*Vj*g_1/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2))))-1.0;
determinant = A*D-B*C;
delta_g_1 = (1.0/determinant)*(D*(-R_1)-B*(-R_2));
delta_g_2 = (1.0/determinant)*(-C*(-R_1)+A*(-R_2));
g_1 = g_1+delta_g_1;
g_2 = g_2+delta_g_2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
}
g_HL = (g_1*g_2)/(g_1+g_2);
V_HL1 = -Vj*(g_2/(g_1+g_2));
V_HL2 = Vj*(g_1/(g_1+g_2));
// Procedure for calculating g_HH
g_1 = 10.0;
g_2 = 10.0;
V_1 = V_H;
V_2 = V_H2;
gamma_1 = gamma_H;
gamma_2 = gamma_H2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
while ( fabs(R_1) > tolerance || fabs(R_2) > tolerance )
{
A = -((gamma_1*Vj*g_2/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2))))-1.0;
B = (gamma_1*Vj*g_1/(V_1*(g_1+g_2)*(g_1+g_2)))*exp((Vj/V_1)*(g_2/(g_1+g_2)));
C = -(gamma_2*Vj*g_2/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2)));
D = ((gamma_2*Vj*g_1/(V_2*(g_1+g_2)*(g_1+g_2)))*exp(-(Vj/V_2)*(g_1/(g_1+g_2))))-1.0;
determinant = A*D-B*C;
delta_g_1 = (1.0/determinant)*(D*(-R_1)-B*(-R_2));
delta_g_2 = (1.0/determinant)*(-C*(-R_1)+A*(-R_2));
g_1 = g_1+delta_g_1;
g_2 = g_2+delta_g_2;
R_1 = gamma_1*exp((Vj/V_1)*(g_2/(g_1+g_2)))-g_1;
R_2 = gamma_2*exp(-(Vj/V_2)*(g_1/(g_1+g_2)))-g_2;
}
g_HH = (g_1*g_2)/(g_1+g_2);
V_HH1 = -Vj*(g_2/(g_1+g_2));
V_HH2 = Vj*(g_1/(g_1+g_2));
alpha_1 = 2.0*alpha_coef/(1.0+exp(-V_LH1/V_alpha));
alpha_2 = 2.0*alpha_coef2/(1.0+exp(-V_HL2/V_alpha2));
alpha_3 = 2.0*alpha_coef/(1.0+exp(-V_LL1/V_alpha));
alpha_4 = 2.0*alpha_coef2/(1.0+exp(-V_LL2/V_alpha2));
beta_1 = beta_coef*exp(-V_HH1/V_beta);
beta_2 = beta_coef2*exp(-V_HH2/V_beta2);
beta_3 = beta_coef*exp(-V_HL1/V_beta);
beta_4 = beta_coef2*exp(-V_LH2/V_beta2);
n_HH=n_HH0;
n_LH=n_LH0;
n_HL=n_HL0;
n_LL=n_LL0;
//Solve differential equations using the forward euler method
while (t < 4000)//Total of 4 second simulation
{
g_junction = Nchans*(n_LL*g_LL+n_LH*g_LH+n_HL*g_HL+n_HH*g_HH)*pow(10.0,-3); //units nanoSiemens
I_junction = g_junction*Vj; //units of picoAmps
n_LL = n_LL + DT*(beta_4*n_LH+beta_3*n_HL-(alpha_3+alpha_4)*n_LL);
n_LH = n_LH + DT*(beta_1*n_HH-(alpha_1+beta_4)*n_LH+alpha_4*n_LL);
n_HL = n_HL + DT*(beta_2*n_HH-(alpha_2+beta_3)*n_HL+alpha_3*n_LL);
n_HH = n_HH + DT*(-(beta_1+beta_2)*n_HH+alpha_1*n_LH+alpha_2*n_HL);
////////////File output code ////////////////
//fprintf(output, "%4.10f\n" ,I_junction);
fprintf(output, "%4.10f\n" ,I_junction);
fprintf(output_minus_10, "%4.10f\n" ,I_junction);
fprintf(output_plus_10, "%4.10f\n" ,I_junction);
if( !(count % 1000) )//1 is true, 0 is false, so to print when modulo is 0, need !
{
if(print_header==0){
fprintf(output,"Time(ms)\tV_j(mV)\tConductance(nS)\tn1H2H\tn1L2H\tn1H2L\tn1L2L\tg1H2H\tg1L2H\tg1H2L\tg1L2L\talpha1\talpha2\talpha3\talpha4\tbeta1\tbeta2\tbeta3\tbeta4\n");
print_header++;
}
fprintf(output, "%f \t%f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f\n"
,t,Vj,g_junction,n_HH,n_LH,n_HL,n_LL,g_HH,g_LH,g_HL,g_LL,alpha_1,alpha_2,alpha_3,alpha_4,beta_1,beta_2,beta_3,beta_4,I_junction);
fprintf(output_minus_10, "%f \t%f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f\n"
,t,modelparams_minus_10[0],modelparams_minus_10[1],modelparams_minus_10[2],modelparams_minus_10[3],modelparams_minus_10[4],modelparams_minus_10[5],modelparams_minus_10[6],modelparams_minus_10[7],modelparams_minus_10[8],modelparams_minus_10[9],modelparams_minus_10[10],modelparams_minus_10[11],modelparams_minus_10[12],modelparams_minus_10[13],modelparams_minus_10[14],modelparams_minus_10[15],beta_3,beta_4,I_junction);
fprintf(output_plus_10, "%f \t%f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f \t%4.10f\n"
,t,modelparams_plus_10[0],modelparams_plus_10[1],modelparams_plus_10[2],modelparams_plus_10[3],modelparams_plus_10[4],modelparams_plus_10[5],modelparams_plus_10[6],modelparams_plus_10[7],modelparams_plus_10[8],modelparams_plus_10[9],modelparams_plus_10[10],modelparams_plus_10[11],modelparams_plus_10[12],modelparams_plus_10[13],modelparams_plus_10[14],modelparams_plus_10[15],beta_3,beta_4,I_junction);
}
t = t+DT;
count = count + 1;
}
}
fclose( output );
return 0;
}
vector<double> fileToVector(const char *name){
vector<double> result;
ifstream input (name);
string lineData;
while(getline(input, lineData))
{
if(lineData.empty()){
cout<<"File empty"<<endl;
}
else{
double d;
double row;
stringstream lineStream(lineData);
while (lineStream >> d){
row = d;
}
result.push_back(row);
}
}
return result;
}