/*
   PRrup.c = filename for this program  (09-21-90)

   This program computes the transport of molecules across a lipid bilayer 
   membrane during an electroporation experiment. The model of the membrane is 
   described in the paper "Electroporation: a Unified Theory". The molecular 
   transport calculation assumes that drift in the electric field is the
   dominant transport mechanism.
*/

#include <stdio.h>
#include <sys/types.h>
#include <time.h>
#include "math.h"
#include "PRcon.h"
  
/*  This program runs the electroporation simulation.  */
/*  This version is a continuous conduction model  */
/*  This version models birth and death of pores and uses cdifc.c */
/*  This version models membrane rupture */

main(argc,argv)
int argc;
char *argv[];
{
/* DEFINE LOCAL VARIABLES */
  FILE *fopen(); /* subroutine in stdio */
/*  FILE  *fp; /* dummy file pointer */
  char *ifile; /* input file name */
  char *ofile1; /* output1 file name */
  char *ofile2; /* output2 file name */
  char *sfile; /* status file name */
  char *ctime(); /* character string used ofr writing time in status file */
  time_t clock, time(); /* used for reading time from system clock */
  double Vcalc(); /* Computes applied voltage */
  double eqct(); /* Steps mambrane voltage forward in time */
  double N;/* number of pres */
  double alp; /* creation rate */
  double Fcalc(); /* calculates creation rate */
  double Acalc(); /* calculates area of membrane covered by pores */
  double area;   /* aqueous area calculated by Acalc() */
  double Rcalc(); /* calculates mean pore radius */
  double R_av; /* mean pore radius calculated by Rcalc() */
  double xb; /* size of largest pore in units of l */
  int err; /* number of arguments read by input routine */
  int lb; /* number of times potential etc has been stored */
  int Nold = 0;
  int cdifc();
  register int i, j, k, ll, ii;
  long offset;
  int michele_count ;

/* READ ARGUMENTS OF PROGRAM CALL */
  if(argc != 5)stop("Error in reading arguments\n");
  ifile = argv[1];
  ofile1 = argv[2];
  ofile2 = argv[3];
  sfile = argv[4];
iprintf("echo: %s %s %s %s\n", argv[1], argv[2], argv[3], argv[4]);

/*  READ INPUT FILE  */
  err = inrd(ifile);  
  if (err != NIN) stop("Error in reading input file\n");

iprintf("/* ALLOCATE STORAGE FOR EXTERNAL VARIABLES */\n") ;
  f = (double *) calloc((6 + nIs) * (n + 2), sizeof(double));
  g = f + n + 2;
  alpha = g + n + 2;
  Gs = alpha + n + 2;
  Gp = Gs + n + 2;
  dVdx = Gp + n + 2;
  if(nIs > 0) gPI = dVdx + n + 2;

iprintf("/* ALLOCATE STORAGE FOR OUTPUT ARRAYS */\n") ;
  fout = (double *) calloc((nout + 1) * (n + 1), sizeof(double));
  nnout = (int *) calloc(nout + 1,sizeof(int));

iprintf("/* INITIALIZE THINGS */\n") ;
  init(); 
  lb = 0;
  michele_count = 0 ;

iprintf("/* LOOP OVER TIME STEPS */\n") ;
  k = 0;
  for(j = 0; j < nt; j++){
    /*time(&clock); */

/* WRITE IN STATUS FILE */
#if 0 
      if(!j){
      fp = fopen(sfile,"w");
      fprintf(fp,"Loop started: nt= %d, j = %d, nout[%d] = %d "
	    ,nt,j,k,nnout[k]);
      fprintf(fp,ctime(&clock));
      offset = ftell(fp);
      fclose(fp);
      fflush(fp);
    } else {
      fp = fopen(sfile,"r+");
      fseek(fp,offset,0);
      fprintf(fp,"Loop started: nt= %d, j = %d, nout[%d] = %d "
	    ,nt,j,k,nnout[k]);
      fprintf(fp,ctime(&clock));
      fflush(fp);
    }
    fclose(fp);
#endif
michele_count = michele_count + 1 ;

if(j % 200 == 0) printf("iteration %d\n", j) ;
/* WRITE DENSITY FUNCTION IN OUTPUT ARRAY */
    if(j == nnout[k]){
      k = k + 1;
      /* CHANGE loop index from k*(n+1) to fix zeroes in last grid space */
      for(i = 0; i < n + 1; i++) fout[(k * (n + 1)) + i] = f[i + 1];
    }

iprintf("/* EVALUATE APPLIED VOLTAGE */\n") ;
    V = Vcalc(j);

iprintf("/* CALCULATE NEW MEMBRANE POTENTIAL */\n") ;
    U = eqct(U,V,tau0,tau1,tau2,dt);
    if((j + 1) % nn == 0){
      UU[(j + 1) / nn] = ((double) (j + 1)) * dt * t0;
      UU[n0 + (j + 1) / nn] = U;
    }

iprintf("/* COMPUTE FORCE AND BIRTH AND DEATH RATES */\n") ;
    alp = Fcalc();

iprintf("/* STEP DENSITY FUNCTION IN TIME */\n") ;
    err = cdifc(f,dVdx,alp,beta,dt,1,n,flag);/* shouldn't 1 be nt? */
    if (err) stop("Error in inverting matrix\n"); /* stop if error */ 

iprintf("/* COMPUTE AREA COVERED BY PORES */\n") ;
    area = Acalc();

iprintf("/* COMPUTE AVERAGE PORE RADIUS */\n") ;
    R_av = Rcalc();

iprintf("/* UPDATE BIG PORE POPULATION */\n") ;
    Nl = Nl + (2. * f[n] - .5 * f[n - 1]) * dt;


iprintf("/* UPDATE MEMBRANE CONDUCTANCE */\n") ;
/*************************************************************************/
/* 10-06-90 Change initial G from 5.e-7 to 8.85e-7; sim ext cell pathway */
/* CHANGE TO G=0.0 TO SIMULATE PLANAR MEMBRANE in BarnettWeaver article  */
/*************************************************************************/
    G = 0.0 ;
    for(ii = 0; ii < nIs; ii++) GPI[ii] = 0.;
    for (i = n; i > 0; i--){  /* Loop through pores */
      G = G + f[i] * g[i];  /* Contribution of small pores  */
      for(ii = 0; ii < nIs; ii++) 
	GPI[ii] = GPI[ii] + gPI[i + (n + 2) * ii] * f[i] * U;
    }
    for(ll = lb; ll >= 0; ll--){  /* Contribution from big pores */
      xbig = xp + (1. - xp) * exp((j - nn * ll) * 2. * Gama * q * q * dt);
      G = G + (pi * sigma * l * xbig * xbig / hh) * (double) bp[ll];
      if(bp[ll] != 0) xb = xbig;
    }
    tau1 = tau0 / (1. + Ri * G);
    tau2 = tau0 / (Ri * G);
    for(ii = 0; ii < nIs; ii++) 
      fPI[ii] = fPI[ii] + GPI[ii];

iprintf("/* STORE DATA FOR POTENTIAL VERSUS TIME PLOT */\n") ;
    if((j + 1)%nn == 0) {
      lb = lb + 1;
      if((bp[(j + 1) / nn] = (Nl - Nold)) >= 1.) Nold = Nl;
      else bp[(j + 1) / nn] = 0;
      UU[2 * n0 + (j + 1) / nn] = G;
      N = 0.;
      for(i = n; i > 0; i--) N = N + f[i];
      UU[3 * n0 + (j + 1) / nn] = N;
      UU[4 * n0 + (j + 1) / nn] = area; /* total aqueous area */
      UU[5 * n0 + (j + 1) / nn] = R_av; /* mean pore radius */
      UU[6 * n0 + (j + 1) / nn] = Nl;
      for(ii = 0; ii < nIs; ii++) 
	UU[(7 + ii) * n0 + (j + 1) / nn] = fPI[ii] * t0 * dt;
      /*wrout(ofile1, ofile2);*/
    }
  /*  if(xb >= xrup) break;*/
  }
   wrout(ofile1,ofile2);
}


stop(st)
char *st;
{
  FILE *fp,*fopen();
  fp = fopen("errfile","w");
  fprintf(fp,st);
  exit();
}
 
double Vcalc(j)
     int j;
{
  double V;

  if(! strcmp(pulse_type, "exponential"))
    V = V0 * exp( - (double) j * t0 * dt / tau_pulse ) ;
  else{
    if( j <= (tpls / dt)) {
      if(! strcmp(pulse_type, "square_wave")) V = V0;
      if(! strcmp(pulse_type, "bipolar")){
	if(j < (tpls / (2*dt))) V = V0 ;
	else V = -V0 ;
      }
    }
    else V = 0.0 ;
  }
  return(V);
}
    
double Fcalc()
{
  int i;
  double x, alp;
  for(i = 0; i < n + 2; i++){
    x = dmin + ((double) i - .5) * q;
    dVdx[i] = (2. * (gama - (Gama + a * alpha[i] * alpha[i] * U * U) * x)) * q;
     }
  alp = crt * exp(a * qp * dmin * dmin * U * U);
  return(alp);
}

double Acalc()
{
  int i;
  double x, area, sum_area = 0;
  for(i = 1; i < n + 2; i++){
    x = dmin + (q * ((double) i - .5));/* x is dimensionless radius, x*l = r */
    area = f[i] * pi * x*x * l*l;
    sum_area += area;
  }
  return(sum_area) ;
}

double Rcalc()
{
  int i;
  double x, r_i_tot, sum_r_i_tot = 0., sum_pores = 0., R_av;
  for(i = n+1; i > 0; i--) {
    x = dmin + (q * ((double) i - .5));
    r_i_tot = f[i] * x * l ;
    sum_r_i_tot += r_i_tot ;
    sum_pores += f[i];
  }
  R_av = sum_r_i_tot / sum_pores;
  return(R_av);
}

/* This version is for use with cdifc.c  */

init()
{
  register int i, ii;
  double parsegian ;
       /* value of Parsegian function determined by Parseg_calc() */
  double Parseg_calc() ;
       /* Parsegian function which enters into chemical potential of an    */
       /*   ion in a pore modeled as a point charge on the axis of an      */
       /*   infinite cylinder, cited in Barnett article and in Pastushenko */
       /*           and Chizmadzhev                                        */
       /*   Parsegian, Nature, vol 221, p 844-846, 1969                    */
       /*   Parsegian, Ann.N.Y.Acad.Sci., vol 264, p 161-174, 1975         */
       /*   Smythe, W. R., Static and Dynamic Electricity, MacGraw-Hill    */
       /*           1968                                                   */
  double x, EkT, sum, xx, H, aa = -2.1, bb = 2.09, cc = -.95, ccon, ff, alp;
  double dd, nm = 1.3e16, Gpp, Gps;
  double fff, sc;
  double Area0, Area_tot, Areai = 0.;
  double radius_tot, sum_radius_i = 0.;
  double dummy;  /* dummy variable used for calculations */
 
  parsegian = Parseg_calc(el) ;
  printf("parsegian = %lf\n", parsegian) ;
/******* change ccon to calculation in the article p. 167  - 910730 ********/
/*   ccon enters calculation of chemical potential of ion in a pore        */
/*   ccon = 0.08 * e * e / (l * kT) ;                                      */
  ccon = e * e / ( el * l * kT);
/*   removed factor of 1/4pi, probably cited with a units error in article */
/***************************************************************************/
  dmin = rmin / l;
  q = (1. - dmin) / ((double) n + .5);
  t0 = l * l * q * q / D;
  gape = gap / kT;
  beta = drat * t0 * exp(- gape) / (q * l);
  crt = crat * t0 * exp(- gape);
  alp = crt;
  hh = h / l;
  gama = pi * l * gmma / kT;
  Gama = pi * l * l * Gamma / kT;
  xp = gama / Gama;
  xrup = rrup / l;
  RR = R / l;
  Ze = Z * e / l;
  tpls = tplse / t0;
  nt = t / (t0 * dt) + 1;
  if((nn = (nt - 1) / (n0 - 1)) == 0){
    nn = 1;
    n0 = nt + 1;
  }
  else n0 = nt / nn + 1;
  if(nIs){
    GPI = (double *) calloc(nIs, sizeof(double));
    fPI = (double *) calloc(nIs, sizeof(double));
  }
  UU = (double *) calloc((7 + nIs) * n0, sizeof(double));
  bp = (int *) calloc(n0, sizeof(int));
  UU[0] = 0.;
  UU[n0] = 0.;
  G = 8.85e-7;
  a = l * (ew - el) / (KK * hh * kT);
  Nl = 0.;
  U = 0.;
  sum = 0.;
  Area_tot = 0.;
  radius_tot = 0.;
  for(i = 0; i < nout; i++) nnout[i] = tout[i] / (t0 * dt);
  if(nout == 0) nnout[0] = -1;
/*  fout[n + 1] = l; */
  f[0] = (alp / beta) * (1. + ((gama - Gama * dmin) * q));
  f[1] = (alp / beta) * (1. - ((gama - Gama * dmin) * q));
/*Area0 = pi * pow((rmin + 0.5 * (l - rmin) / (double)n), 2.0) ;*/
  for(i = n + 1; i > 0; i--){
    x = dmin + (q * ((double) i - .5));
    fout[i - 1] = x * l;
    dd = dmin + q / 2.;
    EkT = (2. * gama * (x - dd) - Gama * ((x * x) - (dd * dd))); 
    /***** changed single ampersand to double ampersand *****/
    if(fabs(EkT) < UFLOW && x < xp) f[i] = f[1] * exp( - EkT);
    else f[i] = 0.;
    sum = sum + f[i];
    Areai = f[i] * x * x * l * l;
    Area_tot += Areai;
    sum_radius_i = f[i] * x * l ;
    radius_tot += sum_radius_i ;
    Gs[i] = 2 * x * l * sigma;
    xx = RR / x;
    if(xx > 1.) H = 0.;
    else H = (x - RR) * (x - RR) *
      (1. + (aa * xx) + (bb * xx * xx * xx) + (cc * xx * xx * xx * xx * xx));
/*****************************************************************/
/*  ff = .02 * ccon * hh / (x * x);                              */
/*  from Jordan (see book chapter)                               */
/*  replacing ff (mu_0 in article) with article's expression for */
/*  Born energy, trying to duplicate mechanical rupture shown    */
/*  maximum value of Parsegian function is 0.25                  */
/*  but value read from graph looks like P(1/40) = 0.1625 (ruler)*/
/*  BUT from Parsegian Ann.N.Y.Acad.Sci. Vol 264, 172            */
/*       P(0.025) = 0.170                                        */

    ff = ccon * parsegian / x ;

/*  change this too... sigma in units of conductivity per length? */
/*   if (fabs(ff) < 30.) Gp[i] = (pi * sigma * l / hh) * H * exp(-ff); */
/*   where did the 30 come from anyway?                           */
/*   x^2 (or r^2) dependence built into H calculations-- ok as is */
/******************************************************************/
    if (fabs(ff) < 30.) Gp[i] = (pi / hh) * sigma * l * H * exp(-ff) ; 
    else Gp[i] = 0;
    alpha[i] = 1. / (1. + Gp[i] / Gs[i]);
    g[i] = (Gs[i] * Gp[i]) / (Gs[i] + Gp[i]);
    G = G + f[i] * g[i];
    rp = x * l;
    Gps = 2 * x * l;
    for(ii = 0; ii < nIs; ii++) {
      if(rp > RPI[ii]){
	if(shape[ii]){
	  A0 = (RPI[ii] * sqrt(rp * rp - RPI[ii] * RPI[ii])    /* Disk-shaped molecules */
		      + 2 * rp * rp * (pi / 2. - asin(RPI[ii] / rp)));
	  sc = 1.;
	} else {
	  xx = RPI[ii] / rp ;
	  if(xx > 1.) H = 0.;
	  else H = (1. + (aa * xx) + (bb * xx * xx * xx) + (cc * xx * xx * xx * xx * xx));
	  A0 = H * pi * (rp * rp - 2 * rp * RPI[ii] + RPI[ii] * RPI[ii]); /* spherical molecules */
	  sc = 1./1.77; /* 1.77 corrects for different coefficient in Stokes's law */
	}
	fff = Zeff[ii] * ff;
	fPI[ii] = 0;
	if (fabs(fff) < 30.) {
	  Gpp = A0 / h * exp(-fff);
	} else {
	  Gpp = 0;
	}
	/*** in gPI, n changed to n+2 to be consistent with allocation ***/
	gPI[i + (n+2) * ii] = ZPI[ii] * sc * nm * (Gpp * Gps / (Gpp + Gps));
      } else {
	gPI[i + (n+2) * ii] = 0;
      }
    }
  }
/*Area_tot += Area0;*/
  Gs[0] = - Gs[1];
  alpha[0] = 1.;
  UU[2 * n0] = G;
  UU[3 * n0] = sum;
  UU[4 * n0] = Area_tot;
  UU[5 * n0] = radius_tot / sum ;
  UU[6 * n0] = 0.;
  tau0 = Ri * C / t0;
  tau1 = tau0 / (1. + (Ri * G));
  tau2 = tau0 / (Ri * G);
  return;
}

/***********************************************************  910730  ***/
/* function Parseg_calc(el) reads the value of the Parsegian function   */
/* from a Table listed in Parsegian, Ann.N.Y.Sci.Acad., vol 264, p 172  */
/* it would be better to modify this to calculate the function itself,  */
/* instead of approximating by connect-the-dots                         */
/************************************************************************/

double Parseg_calc(el)
     double el ;
{
  double parsegian ;
  double R ;

  R = el / ew ;

  if (R < 0.005 || R == 0.005)
    parsegian = (0.095/0.005) * R ;
  else if ( (R > 0.005 && R < 0.01) || R == 0.01)
    parsegian = 0.095 + ((0.123 - 0.095)/(0.01 - 0.005) * (R - 0.005)) ;
  else if ( (R > 0.01 && R < 0.015) || R == 0.015)
    parsegian = 0.123 + ((0.143 - 0.123)/(0.015 - 0.01) * (R - 0.01)) ;
  else if ( (R > 0.015 && R < 0.02) || R == 0.02)
    parsegian = 0.143 + ((0.158 - 0.143)/(0.02 - 0.015) * (R - 0.015)) ;
  else if ( (R > 0.02 && R < 0.025) || R == 0.025)
    parsegian = 0.158 + ((0.170 - 0.158)/(0.025 - 0.02) * (R - 0.02)) ;
  else if ( (R > 0.025 && R < 0.03) || R == 0.03)
    parsegian = 0.170 + ((0.1795 - 0.17)/(0.03 - 0.025) * (R - 0.025)) ;
  else if ( (R > 0.03 && R < 0.035) || R == 0.035)
    parsegian = 0.1795 + ((0.188 - 0.1795)/(0.035 - 0.03) * (R - 0.03)) ;
  else if ( (R > 0.035 && R < 0.04) || R == 0.04)
    parsegian = 0.188 + ((0.195 - 0.188)/(0.04 - 0.035) * (R - 0.035)) ;
  else if ( (R > 0.04 && R < 0.045) || R == 0.045)
    parsegian = 0.195 + ((0.201 - 0.195)/(0.045 - 0.04) * (R - 0.04)) ;
  else if ( (R > 0.045 && R < 0.05) || R == 0.05)
    parsegian = 0.201 + ((0.207 - 0.201)/(0.05 - 0.045) * (R - 0.045)) ;
  else if ( (R > 0.05 && R < 0.055) || R == 0.055)
    parsegian = 0.207 + ((0.212 - 0.207)/(0.055 - 0.05) * (R - 0.05)) ;
  else if ( (R > 0.055 && R < 0.06) || R == 0.06)
    parsegian = 0.212 + ((0.216 - 0.212)/(0.06 - 0.055) * (R - 0.055)) ;
  else if ( (R > 0.06 && R < 0.065) || R == 0.065)
    parsegian = 0.216 + ((0.220 - 0.216)/(0.065 - 0.06) * (R - 0.06)) ;
  else if ( (R > 0.065 && R < 0.07) || R == 0.07)
    parsegian = 0.220 + ((0.223 - 0.220)/(0.07 - 0.065) * (R - 0.065)) ;
  else if ( (R > 0.07 && R < 0.075) || R == 0.075)
    parsegian = 0.223 + ((0.227 - 0.223)/(0.075 - 0.07) * (R - 0.07)) ;
  else if ( (R > 0.075 && R < 0.08) || R == 0.08)
    parsegian = 0.227 + ((0.230 - 0.227)/(0.08 - 0.075) * (R - 0.075)) ;
  else if ( (R > 0.08 && R < 0.085) || R == 0.085)
    parsegian = 0.230 + ((0.232 - 0.230)/(0.085 - 0.08) * (R - 0.08)) ;
  else if ( (R > 0.085 && R < 0.09) || R == 0.09)
    parsegian = 0.232 + ((0.235 - 0.232)/(0.09 - 0.085) * (R - 0.085)) ;
  else if ( (R > 0.09 && R < 0.095) || R == 0.095)
    parsegian = 0.235 + ((0.237 - 0.235)/(0.095 - 0.09) * (R - 0.09)) ;
  else if ( (R > 0.095 && R < 0.1) || R == 0.1)
    parsegian = 0.237 + ((0.239 - 0.237)/(0.1 - 0.095) * (R - 0.095)) ;
  else if ( (R > 0.1 && R < 0.15) || R == 0.15)
    parsegian = 0.239 + ((0.249 - 0.239)/(0.15 - 0.1) * (R - 0.1)) ;
  else if ( (R > 0.15 && R < 0.2) || R == 0.2)
    parsegian = 0.249 ;
  else if ( (R > 0.2 && R < 0.25) || R == 0.25)
    parsegian = 0.249 + ((0.244 - 0.249)/(0.25 - 0.2) * (R - 0.2)) ;
  else if ( (R > 0.25 && R < 0.3) || R == 0.3)
    parsegian = 0.244 + ((0.235 - 0.244)/(0.3 - 0.25) * (R - 0.25)) ;
  else if ( (R > 0.3 && R < 0.35) || R == 0.35)
    parsegian = 0.235 + ((0.224 - 0.235)/(0.35 - 0.3) * (R - 0.3)) ;
  else if ( (R > 0.35 && R < 0.4) || R == 0.4)
    parsegian = 0.224 + ((0.212 - 0.224)/(0.4 - 0.35) * (R - 0.35)) ;
  else if ( (R > 0.4 && R < 0.45) || R == 0.45)
    parsegian = 0.212 + ((0.197 - 0.212)/(0.45 - 0.4) * (R - 0.4)) ;
  else if ( (R > 0.45 && R < 0.5) || R == 0.5)
    parsegian = 0.197 + ((0.182 - 0.197)/(0.5 - 0.45) * (R - 0.45)) ;
  else if ( (R > 0.5 && R < 0.55) || R == 0.55)
    parsegian = 0.182 + ((0.166 - 0.182)/(0.55 - 0.5) * (R - 0.5)) ;
  else if ( (R > 0.55 && R < 0.6) || R == 0.6)
    parsegian = 0.166 + ((0.149 - 0.166)/(0.6 - 0.55) * (R - 0.55)) ;
  else if ( (R > 0.6 && R < 0.65) || R == 0.65)
    parsegian = 0.149 + ((0.131 - 0.149)/(0.65 - 0.6) * (R - 0.6)) ;
  else if ( (R > 0.65 && R < 0.7) || R == 0.7)
    parsegian = 0.131 + ((0.114 - 0.131)/(0.7 - 0.65) * (R - 0.65)) ;
  else if ( (R > 0.7 && R < 0.75) || R == 0.75)
    parsegian = 0.114 + ((0.0969 - 0.114)/(0.75 - 0.7) * (R - 0.7)) ;
  else if ( (R > 0.75 && R < 0.8) || R == 0.8)
    parsegian = 0.0969 + ((0.0769 - 0.0969)/(0.8 - 0.75) * (R - 0.75)) ;
  else if ( (R > 0.8 && R < 0.85) || R == 0.85)
    parsegian = 0.0769 + ((0.0579 - 0.0769)/(0.85 - 0.8) * (R - 0.8)) ;
  else if ( (R > 0.85 && R < 0.9) || R == 0.9)
    parsegian = 0.0579 + ((0.0388 - 0.0579)/(0.9 - 0.85) * (R - 0.85)) ;
  else if ( (R > 0.9 && R < 0.95) || R == 0.95)
    parsegian = 0.0388 + ((0.0197 - 0.0388)/(0.95 - 0.9) * (R - 0.9)) ;
  else if ( (R > 0.95 && R < 1.0) || R == 1.0)
    parsegian = 0.0197 + ((0.0 - 0.197)/(1.0 - 0.95) * (R - 0.95)) ;
  else{
    printf("Invalid input parameter el = %lf, Ratio el/ew = %lf\nAbort...\n", el
, R) ;
    exit(0) ;
  }

  return(parsegian) ;
}