double am(double v)
{
	return(     ((127.0/105.0*v)+(201.0/7.0))  / ( 10.0-(10.0*exp((-201.0/70.0)-(127.0/1050.0*v))  )    )         ) ;
}

double bm(double v)
{
	return(4.0*  exp( (-188.0/63.0)-(127.0/1890.0*v) )         );
}

double mbar(double v)
{
	return( am(v)/(am(v)+bm(v) ) );
}


double ah(double v)
{
        return(7.0/100.0 *  exp( (-94.0/35.0)-(127.0/2100.0*v) )         );
}

double bh(double v)
{
        return(1.0/ (1.0+ exp( (-83.0/35.0)-(127.0/1050.0*v))   )         );
}
double hbar(double v)
{
	return( ah(v)/(ah(v)+bh(v) ) );
}


double an(double v)
{
        return(     ((127.0/105.0*v)+(166.0/7.0))  / ( 100.0-(100.0*exp((-83.0/35.0)-(127.0/1050.0*v))  )    )         ) ;
}

double bn(double v)
{
        return(1.0/8.0 *  exp( (-59.0/140.0)-(127.0/8400.0*v) )         );
}

double nbar(double v)
{
	return( an(v)/(an(v)+bn(v) ) );
}


double ma(double v)
{
        return(1.0/ (1.0+ exp( (v-va)/sa  )  )         );
}


double hai(double v)
{
        return(1.0/ (1.0+ exp( (v-vb)/sb  )  )         );
}



double zv(double v)
{
        return(1.0/ (1.0+ exp(-slope_ca_act/100.0* (v-zb) )  )         );
}



double psbar(double v)
{
        return ( 1 / ( exp(-2.0*(v+45.0)) + 1 ) ); //note its v not vs//0.15//-50
} 
double pstau(double v)
{
	return ((tps/(1+exp (-(v+50)/psst) )  ) + tpsmin ) ;
}

double bbar(double v)
{
        return ( 1 / ( exp(-2.0*(v+42.0)) + 1 ) ); //note its v not vs//0.15//-50
} 

double d_c (double v, double c, double z)
{
	//fprintf(rate_ca,"%lf\t%lf\n",c,1000*(rho*  (      ( (kca*z*(vca-v))/(1.0+2.0*c) ) - c      )    ));
	return(rho*  (      ( (kca*z*(vca-v))/(1.0+2.0*c) ) - c      )    );
}


double d_n (double v, double n )
{
	return(lambdan * (an(v) *(1.0-n) - bn(v) * n)       );
}

double d_h (double v, double h )
{
        return(lambdah * (ah(v) *(1.0-h) - bh(v) * h)       );
}


double d_z (double v, double z)
{
	return(  (zv(v)-z )/tauz );
}


double d_ha (double v, double ha)
{
        return(  (hai(v)-ha )*ka );
}

double f_ps(double v, double ps)
{
        return ( (psbar(v) - ps)/pstau(v));
}
double f_kf(double v, double b)
{
        return ( (bbar(v) - b)/btau);
}


double ina(double v, double h)
{
	return(gna * pow(mbar(v),3)*h*(v-vna) );
}


double ica(double v, double c,double z )
{
	tm_ica=gca * z/(ca_inact_var + c) *(v-vca);
        return(gca * z/(ca_inact_var + c) *(v-vca) ); //0.5
}

double ik(double v, double n )
{
	return(gk * pow(n,4)*(v-vk) );
}



double ikca(double v, double c)
{
	 tm_ikca=gkca * c/(0.5+c) * (v-vk) ;
	return (    gkca * c/(0.5+c) * (v-vk)     );
}

double ia(double v, double ha )
{
        return(ga * pow(ma(v),3)*ha * (v-vk) );
}

double i_ps(v,ps)
double v,ps;
{
	//t_ips=fopen("tm_ips.xls","w");
	//ips= G_PS * ps * (v-vk);
	//fprintf(t_ips,"%lf\t%lf\n",time_,ips);
        return (G_PS * ps * (v-vk));
}

double i_kf(v,b)
double v,b;
{
	return (G_kf * b * (v-vk));
}

double isyn(double v )
{
        return (gsyn * (v-vsyn) );
}


double il(double v )
{
        return (gl * (v-vl) );
}
double il_s(double v )
{
        return (gl_s * (v-vl) );
}
double il_pd(double v )
{
        return (gl_pd * (v-vl) );
}
double il_a(double v )
{
        return (gl_a * (v-vl) );
}

double f_vs(double vs, double vpd, double Va)
{
	
	return (iext_s - (isyn(vs) +  il_s(vs) + Gspd * (vs-vpd) + Gsa * (vs-Va)) /CMs);
}

double f_vpd(double vpd, double vs,double vd,double Va,double n, double h,double ps,double mk)
{
	
	return ( - ( il_pd(vpd) + i_ps(vpd,ps) - Gdpd * (vd - vpd) - Gapd* (Va-vpd) - Gspd * (vs-vpd)) /CMd);
}

double f_vd(double vpd, double vd,double c, double z,double ha,double Va, double b)
{
	//imk= i_mk(vd,mk);
	return( iext_d -( i_kf(vd,b) + ica(vd,c,z) + ikca(vd,c) + ia(vd,ha) + il(vd) + Gdpd * (vd - vpd) + Gda*(vd-Va) ) /CMd);
}

double f_va(double vpd, double Va, double n, double h, double vd, double vs)
{
	return( iext - ( ina(Va,h) + ik(Va,n) + il_a(Va) + Gapd* (Va-vpd) - Gda*(vd-Va)-Gsa * (vs-Va) ) /CMa);
}

void deriv_()
{
        // first neuron
        deriv[0] = f_vs(arg[0],arg[9],arg[7]);
	deriv[9] = f_vpd(arg[9],arg[0],arg[6],arg[7],arg[2],arg[3],arg[8],arg[10]);
	deriv[6] = f_vd(arg[9],arg[6],arg[1],arg[4],arg[5],arg[7], arg[10]);
	deriv[7] = f_va(arg[9],arg[7],arg[2],arg[3],arg[6],arg[0]);
	
	deriv[2] = d_n(arg[7],arg[2]);
	deriv[3] = d_h(arg[7],arg[3]);
        
	deriv[1] = d_c(arg[6],arg[1],arg[4]);
        deriv[4] = d_z(arg[6],arg[4]);
	deriv[5] = d_ha(arg[6],arg[5]);
	
	deriv[8] = f_ps(arg[9],arg[8]);
	deriv[10] = f_kf(arg[6],arg[10]);
}

void rkm()
{
        double old[N], epsilon[N];
        double k1[N], k2[N], k3[N], k4[N],k5[N];
        int i;

        start:
        for (i=0; i < N; i++) {
                arg[i]=state[i];
                //printf("%f\t",state[i]);
                }
        deriv_();
        for (i=0 ;i < N; i++)
        {
                old[i] = state[i];
                k1[i] = step*deriv[i]/3;
                arg[i] = state[i] + k1[i];
        }

        deriv_();
        for (i=0 ;i < N; i++)
        {
                k2[i] = step*deriv[i]/3;
                arg[i] = state[i] + k1[i]/2 +k2[i]/2;
        }

        deriv_();
        for (i=0; i < N; i++)
        {
                k3[i] = step*deriv[i]/3;
                arg[i] = state[i] + 3*k1[i]/8 +9*k3[i]/8;
        }

        deriv_();
        for (i=0; i < N; i++)
        {
                k4[i] = step*deriv[i]/3;
                arg[i] = state[i] + 3*k1[i]/2 - 9*k3[i]/2 + 6*k4[i];
        }

        deriv_();
        for (i=0; i < N; i++)
        {
                k5[i] = step*deriv[i]/3;
                state[i] = state[i] + (k1[i] +4*k4[i] +k5[i])/2 ;
                epsilon[i] = (k1[i] - 9*k3[i]/2 + 4*k4[i] - k5[i]/2)/5;
        }
 if(fabs(epsilon[0]) >= E_MAX || fabs(epsilon[4]) >= E_MAX)
        {
                for(i=0; i < N; i++)
                        state[i] = old[i];
                step = step/2;
                goto start;
        }
        time_ = time_ + step;
        if (fabs(epsilon[0]) <= E_MIN || fabs(epsilon[4]) <= E_MIN)
                step = step*2;
}

void rkm_()
{
        double k1[N],k2[N],k3[N],k4[N];
        int i;

        for(i=0;i<N;i++)
        {
        arg[i] = state[i];
        //printf("%d\t%lf\n",i,state[i]);
        }

deriv_();

        for(i=0;i<N;i++)
        {
        k1[i] = step *(deriv[i]);
        arg[i] = state[i]+k1[i]/2;
        }

deriv_();

        for(i=0;i<N;i++)
        {
        k2[i] = step * deriv[i];
        arg[i] = state[i]+ k2[i]/2;
        }


deriv_();

        for(i=0;i<N;i++)
        {
        k3[i] = step * deriv[i];
        arg[i] = state[i]+ k3[i];
        }

deriv_();

        for(i=0;i<N;i++)
        {
        k4[i] = step * deriv[i];
        state[i] = state[i] + (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
        }

        time_=time_+step;


}


cur()
{
FILE *psinf,*ps_tau,*ps, *vnull,*vaxonnull, *it, *ip, *itot,*zvbar,*ca,*var , *hinf , *nhrel; 
double v,p,i,q,m,k=0.43,a=0.5,s,t,u,iapd=0;
double n,h;
//vnull = fopen("vnull.xls","w");
//vaxonnull = fopen("vaxonnull.xls","w");
//
//ca = fopen("canull.xls","w"); 
//var = fopen("var.xls","w");
//zvbar = fopen("zvbar.xls","w");   
//
//psinf = fopen("psbar.xls","w");
//ps_tau = fopen("pstau.xls","w");
//
//hinf = fopen("hbar.xls","w");
//
//nhrel = fopen("nhrelation.xls","w");
	for(v=-80.0;v<=60;v=v+0.1)
	{
		
		//fprintf(psinf,"%lf\t%lf\n",v,  psbar(v) );
		//fprintf(ps_tau,"%lf\t%lf\n",v,  pstau(v) );
	
		m=kca*zv(v)*(vca-v) ;		
		//fprintf(ca,"%lf\t%lf\n",  (-1 + sqrt(1+8*m) )/4 ,v );
		//fprintf(ca,"%lf\t%lf\n",v,  m );
		//fprintf(zvbar,"%lf\t%lf\n",v,  zv(v) );
		//i = (ina(v,hbar(v))+ ik(v,nbar(v)) + ia(v,hai(v))+ il(v) );
		i = ( ia(v,hai(v))+ i_kf(v,bbar(v) ) + il(v) );
		q = gca * zv(v) * (v-vca);
		p = gkca * (v-vk);
		
		s = p+i;
		t = p*k + q + i*(a+k);
		u = a*(q+i*k);
		
		//fprintf(vnull,"%lf\t%lf\n",v,(-i*a-q)/(i+p));
		//fprintf(vnull,"%lf\t%lf\n",v,-(a*(i+q-w))/(-w+i+q+p));
		//fprintf(vnull, "%lf\t%lf\n",(-t + sqrt( pow(t,2) -4*s*u )) / (2*s) ,v      ); //(2*s) very important 5/2*6 is not the same as 5/(2*6)
		//fprintf(var,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",v,i,q,p,s,t,u);
	
		//fprintf(vaxonnull, "%lf\t%lf\n",v, (iext - gk* pow(nbar(v),4)*(v-vk) - il_a(v)) /(gna* pow(mbar(v),3)*(v-vna)));
		//fprintf(hinf,"%lf\t%lf\n",v,hbar(v));
		
	}

		for(v=-74.0;v<=60;v=v+0.001)
		{
			iapd = 0.227924;
			//fprintf(vaxonnull, "%lf\t%lf\n",v, (iext  - iapd - gk* pow(nbar(v),4)*(v-vk) - il_a(v)) /(gna* pow(mbar(v),3)*(v-vna)));//-0.274053
		}


	for(n=0;n<0.6;n=n+0.01)
	{
		h=-1.7253*n+0.9828; //n=(0.9822-h)/1.7253
		//fprintf(nhrel,"%lf\t%lf\n",n,h);
	}

}