/*  traub1.c	This is a rough c translation of traub1.m

    to compile on the CAAM server type

         gcc traub1.c -o traub1exe -lm

    from the Unix prompt. To execute type

         traub1exe

    this will write the solution to traub1out

    you may read and view this in matlab by typing

   >> load traub1out -ascii

   >> plot(traub1out(:,1))

    similarly for the other columns */

#include <stdio.h>
#include <math.h>

double am(double);
double bm(double);
double ah(double);
double bh(double);
double an(double);
double bn(double);
double ay(double);
double by(double);
double as(double);
double bs(double);
double ar(double);
double br(double);
double aq(double,double);
double bq(double,double);

main()
{
        int j;

        double v,xi,ina,ik,ica,ikca,dt,T,m,n,h,top,bot,Is,tmp,
         vL,vna,vk,vca,gL,gna,gca,gk,gkca,Cm,c,A,d,D,xib,
         y,s,r,q,gpk,gpca,gpna,gpkca;


         FILE *fp;

fp = fopen("traub1out","w");

vL = -15;
vna = 115;
vk = -15;
vca = 140;

gL = .01;
gna = 3.32;
gca = 6.64;
gk = 3.98;
gkca = 0.1;  

Cm = 3*3320*1e-5;	

c = 5200; 
A = 3320;
d = 5e-4;
D = c/(A*d);

xib = 0.1; 

j = 1; dt = 0.025;  T=60;

v = 0;
xi = v;
ina = v;
ik = v;
ica = v;
ikca = v;

v = vL;

m = am(v)/(am(v)+bm(v));
h = ah(v)/(ah(v)+bh(v));
n = an(v)/(an(v)+bn(v));
y = ay(v)/(ay(v)+by(v));
s = as(v)/(as(v)+bs(v));
r = ar(v)/(ar(v)+br(v));
q = aq(v,xi)/(aq(v,xi)+bq(v,xi)); 

while ( j*dt < T ){
   
    Is = 0;
    if (j*dt > 5 && j*dt < 25) {
        Is = 0.3; 	
    }
    
    tmp = am(v);
    m = (dt*tmp + m)/(1+dt*(tmp+bm(v)));
    
    tmp = ah(v);
    h = (dt*tmp+h)/(1+dt*(tmp+bh(v)));
    
    tmp = an(v);
    n = (dt*tmp+n)/(1+dt*(tmp+bn(v)));
    
    tmp = ay(v);
    y = (dt*tmp+y)/(1+dt*(tmp+by(v)));
    
    tmp = as(v);
    s = (dt*tmp+s)/(1+dt*(tmp+bs(v)));
    
    tmp = ar(v);
    r = (dt*tmp+r)/(1+dt*(tmp+br(xi)));
    
    tmp = aq(v,xi);
    q = (dt*tmp+q)/(1+dt*(tmp+bq(v,xi)));
    
    gpna = gna*pow(m,3)*h;
    gpk = gk*pow(n,4)*y;
    gpca = gca*pow(s,5)*r;
    gpkca = gkca*q;
    
    xi = (xi - dt*D*(1e-3)*gpca*(v-vca))/(1+dt*xib);
        
    top =  v + (dt/Cm)*(Is+gL*vL+gpna*vna+gpk*vk+gpca*vca+gpkca*vk);
    
    bot = 1 + (dt/Cm)*(gL + gpna + gpk + gpca + gpkca);
    
    v = top/bot;    

    ina = gpna*(v-vna);
    ik = gpk*(v-vk);
    ica = gpca*(v-vca);
    ikca = gpkca*(v-vk);

    fprintf(fp,"%lf %lf %lf %lf %lf \n",v,xi,ina,ik,ica);
    
    j = j + 1;
    
}   /* end while */

fclose(fp);

}  /* end main */

double am(double v)
{
double val;
val = 0.32*(13-v)/(exp((13-v)/4)-1);
return val;
}

double bm(double v)
{
double val;
val = 0.28*(v-40)/(exp((v-40)/5)-1);
return val;
}

double ah(double v)
{
double val;
val = 0.128*exp((17-v)/18);
return val;
}

double bh(double v)
{
double val;
val = 4/(exp((40-v)/5)+1);
return val;
}

double an(double v)
{
double val;
val = 0.032*(15-v)/(exp((15-v)/5)-1);
return val;
}

double bn(double v)
{
double val;
val = 0.5*exp((10-v)/40);
return val;
}

double ay(double v)
{
double val;
val = 0.028*exp((15-v)/15) + 2/(exp((85-v)/10)+1);
return val;
}

double by(double v)
{
double val;
val = 0.4/(exp((40-v)/10)+1);
return val;
}

double as(double v)
{
double val;
val = 0.04*(60-v)/(exp((60-v)/10)-1);
return val;
}

double bs(double v)
{
double val;
val = 0.005*(v-45)/(exp((v-45)/10)-1);
return val;
}

double ar(double v)
{
double val;
val = 0.005;
return val;
}

double br(double xi)
{
double val;
val = 0.025*(200-xi)/(exp((200-xi)/20)-1);
return val;
}

double aq(double v,double xi)
{
double val;
val = exp(v/27)*0.005*(200-xi)/(exp((200-xi)/20)-1);
return val;
}

double bq(double v,double xi)
{
double val;
val = 0.002;
return val;
}