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

    to compile on the CAAM server type

         gcc iaf2net.c -o iaf2netexe -lm

    from the Unix prompt. To execute type

         iaf2netexe

    this will write the solution to iaf2netout

    you may read and view this in matlab by typing

   >> raster

    once you have copied raster.m to your directory */

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

main()
{
        int i,j,k;

        int fire[100];

        double EL, Vreset, taum, rm, A, Rm, Ie, EK, tausra, ginc, T,
               Vth, dt, s, st, t, Eex, Einh, xx;

        double v[100],          /* the right hand side */
               gsra[100],           /* the solution */
               S[100], tim[100],
               gsyn[100], gsynE[100], I[100];

        double W[100][100],     /* the weight matrix */
               WE[100][100];    /* the weight.*Esyn matrix */


         FILE *fp;

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

EL = -65;
Vreset = -65;
taum = 10;
rm = 1; 
A = .1;
Rm = rm/A; 
Ie = 10;  
EK = -70;
tausra = 10;
ginc = 50*0.06/rm;
T = 200;
Vth = -50;
dt = 0.04;

s = taum/dt;
st = tausra/dt;

Eex = 0;
Einh = -70;

for (i=0;i<=6;i=i+1){
    S[i] = 0;
    tim[i] = -100;
    gsra[i] = 0;
    I[i] = 0;
    v[i] = EL;
    for (j=0;j<=6;j=j+1){
        W[i][j] = 0;
        WE[i][j] = 0;
    }
}

  /* W(3:5,1:2) = 1;  Esyn(3:5,1:2) = Eex; */

W[2][0] = 10; W[2][1] = 10; WE[2][0] = Eex*W[2][0]; WE[2][1] = Eex*W[2][1];
W[3][0] = 10; W[3][1] = 10; WE[3][0] = Eex*W[3][0]; WE[3][1] = Eex*W[2][1];
W[4][0] = 10; W[4][1] = 10; WE[4][0] = Eex*W[4][0]; WE[4][1] = Eex*W[2][1];

   /*  W(6:7,3:5) = 1;  Esyn(6:7,3:4) = Eex;  */

W[5][2] = 20; W[5][3] = 20; WE[5][2] = Eex*W[5][2]; WE[5][3] = Eex*W[5][3];
W[6][2] = 20; W[6][3] = 20; WE[6][2] = Eex*W[6][2]; WE[6][3] = Eex*W[6][3];

  /*  W(6:7,5) = 0.1;  Esyn(6:7,5) = Einh;  */

W[5][4] = 1; WE[5][4] = Einh*W[5][4];
W[6][4] = 1; WE[6][4] = Einh*W[6][4];

t = 0;
j = 0;

while ( j*dt < T ){
   
    j = j + 1;
    t = (j-1)*dt;

    I[0] = 0;
    I[1] = 0;

    if (t > 1) {
       I[0] = Ie;
       I[1] = Ie;
    }

   for (i=0;i<=6;i=i+1)
       gsra[i] = st*gsra[i]/(1+st);

   /* form gsyn = rm*W*S   and  gsynE = rm*WE*S  */

   for (i=0;i<=6;i=i+1){
       gsyn[i] = 0;
       gsynE[i] = 0;
       for (k=0;k<=6;k=k+1){
           gsyn[i] = gsyn[i] + rm*W[i][k]*S[k];
           gsynE[i] = gsynE[i] + rm*WE[i][k]*S[k];
       }
   }

   for (i=0;i<=6;i=i+1)
       v[i]=(s*v[i]+EL+Rm*I[i]+rm*gsra[i]*EK+gsynE[i])/(s+1+rm*gsra[i]+gsyn[i]);

   xx = step(2);

   for (i=0;i<=6;i=i+1){

       fire[i] = 0;

       if (v[i] > Vth) {

          fire[i] = 1;

          tim[i] = t;

          v[i] = Vreset;

          gsra[i] = (1+st)*gsra[i]/st + ginc;

        }

   
        S[i] = (t-tim[i])*exp(3*(tim[i]-t));  

        fprintf(fp, "%d\n", fire[i]);

      /*   fprintf(fp,"%lf %lf %lf\n",v[i], gsra[i], S[i]);  */
    }

}   /* end while */

fclose(fp);

}  /* end main */