/* afvp.c                                                                     */
/* Copyright (c)2006 by  Jie Lian, All Rights Reserved                        */
/* This program simulates the interaction between AF and Vp.                  */
/* See Lian J., Mussig D., Lang V.  Computer Modeling of Ventricular          */
/* Rhythm  During Atrial Fibrillation and Ventricular Pacing.                 */
/* IEEE Transactions on Biomedical Engineering 53(8): 1512-1520, 2006.        */ 
/* Contact J. Lian  jie.lian AT biotronik DOT com)                            */
/*                                                                            */      
/*   This program is free software; you can redistribute it and/or modify     */ 
/*   it under the terms of the GNU General Public License as published by     */ 
/*   the Free Software Foundation; either version 2 of the License, or        */ 
/*   (at your option) any later version.                                      */ 
/*                                                                            */ 
/*   This program is distributed in the hope that it will be useful,          */ 
/*   but WITHOUT ANY WARRANTY; without even the implied warranty of           */ 
/*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */ 
/*   GNU General Public License for more details.                             */
/*                                                                            */      
/*   You should have received a copy of the GNU General Public License        */
/*   along with this program; if not, write to the Free Software Foundation   */
/*   Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA            */ 
/*                                                                            */ 
/*  afvp.c   and its dependents are freely available from Physionet -         */ 
/*  http://www.physionet.org/ - please report any bugs to the author above.   */

/* Last modified: Jul 12th 2006 by Gari Clifford:                             */
/*   - Adapted random number calls to make it ANSI C compliant                */ 
/*   - Added random number seeding and extra command line arguments           */ 
/*   - Removed BI (standby ventricular pacing basic interval) command line    */
/*     option                                                                 */ 


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>    // For random()

#include "afvp.h"    // header file for the AF-VP model
#include "conio.h"   // sometimes a missing system file


// definition of global model variables
EnvStruct  env;     // simulation environment
AtrStruct  atr;     // atrium model
AvjStruct  avj;     // AVJ model
VtrStruct  vtr;     // ventricle model
VeleStruct vele;    // ventricular electrode model

unsigned int rseed; // seed for random number generator 
                    // -can be passed as an argument

int main(int argc, char *argv[])
{
  time_t utc;

  if (argc<2){
    fprintf(stdout,"Usage: %s <config_file> <random number seed>\n",argv[0]);
    fprintf(stdout,"\n config file example: http://www.physionet.org/physiotools/afvp/config.txt \n");
    //    fprintf(stdout," BI = standby ventricular pacing basic interval in ms (default = 800 ms) \n");
    fprintf(stdout," rseed = seed for random number generator - specify a long int for repeatability \n\n");
    exit(1);
  }

  ReadConfig(argv[1]);        // read model configuration file
  //if (argc==3) vele.BI = atoi(argv[2])/1000.0;   // overwrite BI 
  if (argc>2){
    rseed   = atol(argv[2]);          // random seed
  }
  else {
   /* use the system clock to give it a random number */
   time(&utc);
   rseed = (long int) utc;
  }
  if (argc>3)  fprintf(stderr,"Ignoring input after 2nd argument\n");

  ModelInit();                // initialize model parameters

  // Main loop, keep running until MAX_RR beats or MAX_TIME timeout
  while ((env.beat<env.MAX_RR) && (env.t<env.MAX_TIME)) {

    UpdateAtTs();             // model update at sampling time

    if (vtr.tmA[0]/vtr.AntDly+vtr.tmR[0]/vtr.RetDly>=1.0)
      VtrFusion();            // detection of ventricular fusion
    if (atr.t>=atr.NextAA)
      AnteHitAvj();           // antegrade AF impluse hits AVJ
    if (vtr.nR>0 && vtr.tmR[0]>=vtr.RetDly)
      RetrHitAvj();           // VP-induced retrograde impulse hits AVJ
    if (vtr.nA>0 && vtr.tmA[0]>=vtr.AntDly)
      VtrSense();             // there is an ante conducted Vs
    if (vele.t>=vele.NextVV)
      VtrPace();              // timeout for a scheduled Vp
    if (avj.tmA[0][0]>=avj.tmA[0][1])
      AnteEscAvj();           // ante AV delay timeout and invade vtr
    if (avj.tmR[0][0]>=avj.tmR[0][1])
      RetrEscAvj();           // retr AV delay timeout and invade atr

    if (avj.phase==PHASE4) {  // check if AVJ depol & starts ref
      if (avj.Vm>=avj.Vt) {
        ActivateAvj();
        StartAvjRef();
      }
    }
    else {                    // check if AVJ ref ends & return to phase 4
      if (avj.t>=avj.ref)
        StartAvjPh4();
    }
  }

  ModelExit();
} /*** end main ***/

/*********************************************************************
ErrMsg: display error message and abort simulation
Input: msg - message to be displayed
       errcode - exit error code
*********************************************************************/
void ErrMsg(char *msg, int errcode)
{
  printf("%s\n",msg);
  exit(errcode);
}

/*********************************************************************
LogEvent: log event time and type
Input: fp - log file pointer
       msg - message to be logged      
*********************************************************************/
void LogEvent(FILE *fp, char *msg)
{
  fprintf(fp,"T=%.3lf '%s ",env.t,msg);
  if (avj.phase==PHASE4) fprintf(fp,"@PH4' ");
  else fprintf(fp,"@REF' ");
}

/*********************************************************************
LogStat: log statistics after simulation
Input: fp - log file pointer
*********************************************************************/
void LogStat(FILE *fp)
{
	long sum;

  fprintf(fp,"\n----------------------------------------------------\n");
        fprintf(fp,"Random seed = %d\n",rseed);
	fprintf(fp,"Total time (s): %.3lf\n", env.t);
	fprintf(fp,"Total beat: %ld\n", env.beat);
	fprintf(fp,"Vs Count = %ld [%.2lf%%]\n",
	       vele.cnt[SENSE],100.0*vele.cnt[SENSE]/env.beat);
	fprintf(fp,"Vp Count = %ld [%.2lf%%]\n",
	       vele.cnt[PACE], 100.0*vele.cnt[PACE]/env.beat);
	sum = atr.nAf[REFRACTORY]+atr.nAf[PHASE4];
	fprintf(fp,"Af Count = %ld(ref)+%ld(ph4) = %ld\n",
	       atr.nAf[REFRACTORY], atr.nAf[PHASE4], sum);
	fprintf(fp,"Mean RR interval = %.3lf <%dppm>\n",
	       env.t/env.beat, (int)(env.beat*60/env.t));
	fprintf(fp,"Mean PP interval = %.3lf <%dppm>\n",
	       env.t/sum, (int)(sum*60/env.t));
	sum = vtr.nR+avj.nR;
	fprintf(fp,"Active retr waves = %ld [%.2lf%%]\n",
	      sum, 100.0*sum/env.beat);
	sum = vtr.nF[REFRACTORY]+vtr.nF[PHASE4];
	fprintf(fp,"Vtr fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
	       vtr.nF[REFRACTORY], vtr.nF[PHASE4], sum, 100.0*sum/env.beat);
	sum = avj.nF[REFRACTORY]+avj.nF[PHASE4];
	fprintf(fp,"AVJ fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
	       avj.nF[REFRACTORY], avj.nF[PHASE4], sum, 100.0*sum/env.beat);
	sum = atr.nF[REFRACTORY]+atr.nF[PHASE4];
	fprintf(fp,"Atr fusions = %ld(ref)+%ld(ph4) = %ld [%.2lf%%]\n",
	       atr.nF[REFRACTORY], atr.nF[PHASE4], sum, 100.0*sum/env.beat);
	sum = avj.nB[RETROGRADE]+vtr.nB[RETROGRADE];
	fprintf(fp,"Retrgrade Block = %ld(avj)+%ld(vtr) = %ld [%.2lf%%]\n",
	       avj.nB[RETROGRADE], vtr.nB[RETROGRADE], sum, 100.0*sum/env.beat);
}

/*********************************************************************
ReadConfig: read configuration file
Input: fn - config filename
*********************************************************************/
void ReadConfig(char *fn)
{
  double d;
  FILE *fp;
  char str[BUFLIM];

  if ((fp=fopen(fn,"rt"))==NULL) ErrMsg("Cannot open config file",1);

  // read parameters of the simulation environment
  ReadLine(fp,-1,env.fnRR);
  ReadLine(fp,-1,env.fnAA);
  ReadLine(fp,-1,env.fnAV);
  ReadLine(fp,-1,env.fnLOG);
  env.MAX_RR = (long)ReadLine(fp,500,str);
  env.MAX_TIME = ReadLine(fp,1000,str);
  env.Ts = ReadLine(fp,0.001,str);
  env.RR0 = ReadLine(fp,1.0,str);

  // read parameters of the atrium model
  atr.AA_MODEL = (int)ReadLine(fp,0,str);
  atr.lambda = ReadLine(fp,10,str);
  atr.AAstd = ReadLine(fp,0,str);
  atr.dVmean = ReadLine(fp,33,str);
  atr.dVstd = ReadLine(fp,0,str);
  atr.AtrDly = ReadLine(fp,0.03,str);
  atr.S1S2 = ReadLine(fp,0.4,str);
  atr.S2S3 = ReadLine(fp,0.4,str);

  // read parameters of the AVJ model
  avj.Vt = ReadLine(fp,-40,str);
  avj.Vr = ReadLine(fp,-90,str);
  avj.dVdt = ReadLine(fp,30,str);
  avj.MinAVDa = ReadLine(fp,0.07,str);
  avj.MinAVDr = ReadLine(fp,0.07,str);
  avj.alpha = ReadLine(fp,0.13,str);
  avj.tau_c = ReadLine(fp,0.1,str);
  avj.MinRef = ReadLine(fp,0.05,str);
  avj.beta = ReadLine(fp,0.25,str);
  avj.tau_r = ReadLine(fp,0.5,str);
  avj.Ref_std = ReadLine(fp,0.0,str);
  avj.delta = ReadLine(fp,10,str);
  avj.theta = ReadLine(fp,10,str);

  // read parameters of the ventricle model
  vtr.AntDly = ReadLine(fp,0.05,str);
  vtr.RetDly = ReadLine(fp,0.15,str);
  vtr.ref = ReadLine(fp,0.1,str);

  // read parameters of the RV electrode model
  vele.VP_MODEL = (int)ReadLine(fp,0,str);
  vele.BI = ReadLine(fp,1.0,str);

  fclose(fp);

  LogConfig(stdout);    // display model parameters on terminal
}

/*********************************************************************
ReadLine: read a line from the config file
Input:  fp - pointer to the config file
        v  - default value
        str - char pointer (string read after '=' sign)
Output: 'val' read from the line if found, otherwise return 'v'
Note: Comment line starts with '%' sign
      Value is separated from text by token '='
*********************************************************************/
double ReadLine(FILE *fp, double v, char *str)
{
  char buf[BUFLIM], *p;
  char tok[] = "=";                        // token followed by 'val'
  double val;                              // value read from the line

  *str = '\0';                             // default: null string
  do {
    p = fgets(buf,sizeof(buf),fp);         // read a line
    if (p==NULL) return v;                 // force return on EOF
  } while (strlen(buf)<=1 || buf[0]=='%'); // skip comment or blank lines
  p = strstr(buf,tok);                     // search for "=" sign
  if (p!=NULL) {                           // if found "="
    p += strlen(tok);                      // skip "=" sign
    sscanf(p,"%s",str);                    // read string after "=" sign
    if (sscanf(p,"%lf",&val)>0) return val;// read value after "=" sign
  }                                        // otherwise
  return v;                                // return default value
}

/*********************************************************************
LogConfig: log model configuration
Input: fp - output file pointer
*********************************************************************/
void LogConfig(FILE *fp)
{
  fprintf(fp,"----------------------------------------------------\n");
  fprintf(fp,"*** Simulation Environment ***\n");
  fprintf(fp,"fnRR=%s\nfnAA=%s\nfnLog=%s\nfnAV=%s\n",
         env.fnRR,env.fnAA,env.fnAV,env.fnLOG);
  fprintf(fp,"MAX_RR=%ld\nMAX_TIME=%lf\nTs=%lf\nRR0=%lf\n",
         env.MAX_RR,env.MAX_TIME,env.Ts,env.RR0);
  fprintf(fp,"\n*** Atrium Model ***\n");
  fprintf(fp,"AA_MODEL=%d\nlambda=%lf\nAAstd=%lf\n",
         atr.AA_MODEL,atr.lambda,atr.AAstd);
  fprintf(fp,"dVmean=%lf\ndVstd=%lf\nAtrDly=%lf\nS1S2=%lf\nS2S3=%lf\n",
         atr.dVmean,atr.dVstd,atr.AtrDly,atr.S1S2,atr.S2S3);
  fprintf(fp,"\n*** AVJ Model ***\n");
  fprintf(fp,"Vt=%lf\nVr=%lf\ndVdt=%lf\nMinAVDa=%lf\nMinAVDr=%lf\n",
         avj.Vt,avj.Vr,avj.dVdt,avj.MinAVDa,avj.MinAVDr);
  fprintf(fp,"alpha=%lf\ntau_c=%lf\nMinRef=%lf\nbeta=%lf\n",
         avj.alpha,avj.tau_c,avj.MinRef,avj.beta);
  fprintf(fp,"tau_r=%lf\nRef_std=%lf\ndelta=%lf\ntheta=%lf\n",
         avj.tau_r,avj.Ref_std,avj.delta,avj.theta);
  fprintf(fp,"\n*** Ventricle Model ***\n");
  fprintf(fp,"AntDly=%lf\nRetDly=%lf\nref=%lf\n",
         vtr.AntDly,vtr.RetDly,vtr.ref);
  fprintf(fp,"\n*** Ventricular Electrode Model ***\n");
  fprintf(fp,"VP_MODEL=%d\nBI=%lf\n",vele.VP_MODEL,vele.BI);
  fprintf(fp,"----------------------------------------------------\n\n");
}

/*********************************************************************
ModelInit: initialization of model parameters
*********************************************************************/
void ModelInit()
{
  int i;

  if ((env.fprr=fopen(env.fnRR,"wt"))==NULL)  // output file of RR intervals
    ErrMsg("Cannot open output file of RR intervals",1);
  if ((env.fpaa=fopen(env.fnAA,"wt"))==NULL)  // output file of AA intervals
    ErrMsg("Cannot open output file of AA intervals",1);
  if ((env.fpav=fopen(env.fnAV,"wt"))==NULL)  // output file of AVJ status
    ErrMsg("Cannot open output file of AVJ status",1);
  if ((env.fplog=fopen(env.fnLOG,"wt"))==NULL)// output file of event log
    ErrMsg("Cannot open output file of event log",1);

  env.beat = 0;                            // clear beat counter
  env.t = 0;                               // reset simulation time clock
  avj.phase = PHASE4;                      // set AVJ phase = PHASE4
  avj.Vm = avj.Vr;                         // set Vm to Vr
  avj.V4 = avj.dVdt * env.Ts;              // initialize V4 (mV/ms)
  avj.RT0 = 0;                             // AVJ start recovery

  avj.nA = avj.nR = 0;                     // no activation wave in AVJ
  vtr.nA = vtr.nR = 0;                     // no activation wave in vtr
  for (i=0;i<MAX_AVJ_TM;i++) {             // disable AVJ activation timers
    avj.tmA[i][0] = -1;
    avj.tmR[i][0] = -1;
  }
  for (i=0;i<MAX_VENT_TM;i++) {            // disable vtr activation timers
    vtr.tmA[i] = -1;
    vtr.tmR[i] = -1;
  }

  vele.cnt[PACE] = 0;                      // clear Vp counter
  vele.cnt[SENSE] = 0;                     // clear Vs counter
  atr.nAf[0] = atr.nAf[1] = 0;             // clear AF event counters
  atr.nF[0] = atr.nF[1] = 0;               // clear atrial fusion counters
  avj.nF[0] = avj.nF[1] = 0;               // clear AVJ fusion counters
  vtr.nF[0] = vtr.nF[1] = 0;               // clear ventricle fusion counters

  vtr.tmA0 = 0;                            // clear last ante AV timeout time
  vtr.tmR0 = 0;                            // clear last Vp delivery time

  // Get the first ventricular pacing interval
  vele.NextVV = GetVpIntv(env.RR0,vele.BI,vele.BI,vele.VP_MODEL);

  // randomize();                             // Initialize randomization
  srandom(rseed);  // Initialize randomization with seeded random number

  // Get the first random interval for next AF depol.
  atr.NextAA = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL);
  fprintf(env.fpaa,"%lf\n",atr.NextAA);
  
  LogConfig(env.fplog);                    // log config parameters
}

/*********************************************************************
ModelExit: exit simulation
*********************************************************************/
void ModelExit()
{
  LogStat(stdout);                   // output statistics to stdout
  LogStat(env.fplog);                // output statistics to log file
  fclose(env.fprr);
  fclose(env.fpaa);
  fclose(env.fpav);
  fclose(env.fplog);
}

/*********************************************************************
UpdateAtTs: model update at sampling time
*********************************************************************/
void UpdateAtTs()
{
  int i;

  env.t += env.Ts;          // update simulation timer
  atr.t += env.Ts;          // update atrial timer
  vtr.t += env.Ts;          // update ventricle timer
  vele.t += env.Ts;         // update Vp timer

  // update AVJ and ventricle activation timers
  for (i=0;i<avj.nA;i++) avj.tmA[i][0] += env.Ts;
  for (i=0;i<avj.nR;i++) avj.tmR[i][0] += env.Ts;
  for (i=0;i<vtr.nA;i++) vtr.tmA[i] += env.Ts;
  for (i=0;i<vtr.nR;i++) vtr.tmR[i] += env.Ts;

  if (avj.phase==PHASE4) {        // during phase 4
    avj.Vm += avj.V4;             // Vm linear rise from Vr to Vt
    if (avj.Vm>=avj.Vt) avj.dir = ANTEGRADE;
  }
  else                            // during refractory period
    avj.t += env.Ts;              // update AVJ refractory timer
}

/*********************************************************************
AnteHitAvj: event service when antegrade impulse hits AVJ
*********************************************************************/
void AnteHitAvj()
{
  double ext;
  long r;

  LogEvent(env.fplog,"AnteHitAvj"); 
  
  atr.t = 0;                // reset atrium timer
  atr.nAf[avj.phase]++;     // increment AF impulse counter
  // Get the random interval for next AF depolarization
  atr.NextAA = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL);
  fprintf(env.fpaa,"%lf\n",atr.NextAA);
  if (avj.phase==PHASE4) {  // AVJ in phase IV
    //srandom(-1000);
    r = (long) random();      // simulate variation of dV
    atr.dV = gaussian(atr.dVmean,atr.dVstd,ran2(&r),ran2(&r));
    if (atr.dV>0) avj.Vm += atr.dV;  // delta increment of Vm
    if (avj.Vm>=avj.Vt) avj.dir = ANTEGRADE;
    fprintf(env.fplog,      // log after AnteHitAvj() service @PH4
      "NextAA=%.3lf, Vm=%.2lf+%.2lf=%.2lf\n",
      atr.NextAA, avj.Vm-atr.dV, atr.dV, avj.Vm);
  }
  else {                    // AVJ in refractory period
    r = atr.dV/(avj.Vt-avj.Vr); // relative strength of the AF impulse
    if (r>1) r = 1.0;       // with upper limit 1.0
    r = pow(avj.t/avj.ref,avj.theta)*pow(r,avj.delta);
    ext = avj.MinRef*r;
    avj.ref += ext;         // ref extension due to concealed AF
    fprintf(env.fplog,      // log after AnteHitAvj() service @REF
      "NextAA=%.3lf, ref=%.3lf+%.3lf=%.3lf\n",
      atr.NextAA, avj.ref-ext, ext, avj.ref);
  }
}

/*********************************************************************
RetrHitAvj: event service when retrograde impulse hits AVJ
*********************************************************************/
void RetrHitAvj()
{
  int i;
  double r,ext;

  LogEvent(env.fplog,"RetrHitAvj"); 

  StopVtrTm(RETROGRADE);               // stop a retr vtr activation timer
  
  if (avj.phase==PHASE4) {             // retr invasion in phase 4
    if (avj.Vm<avj.Vt) {               // retr AVJ depolarization
      avj.Vm = avj.Vt;                 // suprathreshold AVJ immediately
      avj.dir = RETROGRADE;            // retrograde activation
      fprintf(env.fplog,               // log after RetHitAvj() service @PH4
        "vtr.nR=%d, Vm=%.2lf Retr AVJ activation\n", vtr.nR, avj.Vm);
    }
    else {        // simultaneous ante/retr depol (rare), no wave in AVJ
      avj.dir = BIDIRECTION;
      fprintf(env.fplog,               // log after RetHitAvj() service @PH4
        "vtr.nR=%d, Vm=%.2lf Bidir AVJ activation\n", vtr.nR, avj.Vm);
    }
  }
  else {                               // retr invasion in refractory period
    if (avj.nA>0) {                    // presence of ante wave in AVJ
      StopAvjTm(ANTEGRADE);            // stop the ante AVJ timer
      avj.nF[avj.phase]++;             // incr AVJ fusion counter
    }
    avj.nB[RETROGRADE]++;              // incr counter of blocked retr waves
    r = pow(avj.t/avj.ref,avj.theta);  // strong impulse
    ext = avj.MinRef*r;
    avj.ref += ext;                    // ref extension by concealed impulse
    fprintf(env.fplog,                 // log after RetHitAvj() service @REF
      "vtr.nR=%d, avj.nA=%d, nB=%d, nF=%d, ref=%.2lf+%.2lf=%.2lf\n", 
      vtr.nR,avj.nA,avj.nB[RETROGRADE],avj.nF[REFRACTORY],
      avj.ref-ext,ext,avj.ref);
  }
}

/*********************************************************************
ActivateAvj: AVJ depolarization service
*********************************************************************/
void ActivateAvj()
{
  LogEvent(env.fplog,"ActivateAvj"); 

  avj.RT = env.t-avj.RT0;             // calculated AVJ recovery time 
  avj.AVDa = avj.MinAVDa+avj.alpha*exp(-avj.RT/avj.tau_c); // ante AV delay
  avj.AVDr = avj.MinAVDr+avj.alpha*exp(-avj.RT/avj.tau_c); // retr AV delay

  if (avj.dir==ANTEGRADE) {           // antegrade AVJ activation
    if (avj.nR==0) {                  // no existing retr activation in AVJ
      StartAvjTm(ANTEGRADE);          // start a new timer for ante AV delay
      fprintf(env.fpav,               // record AVJ status
        "%.3lf %d %.3lf %.3lf ", env.t, avj.dir, avj.RT, avj.AVDa);
      fprintf(env.fplog,              // log after ActivateAvj() service
        "start AVDa = %.3lf, avj.nA=%d\n", avj.AVDa, avj.nA);     
    }
    else {                            // collide w/ existing retr wave
      StopAvjTm(RETROGRADE);          // stop a retr AVJ activation timer
      avj.nF[avj.phase]++;            // incr AVJ collision counter (Phase 4)
      fprintf(env.fpav,               // record AVJ status
        "%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
      fprintf(env.fplog,              // log after ActivateAvj() service
        "collide w/ retr wave, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
    }
  }
  else if (avj.dir==RETROGRADE) {     // retrograde AVJ activation
    if (avj.nA==0) {                  // no existing ante activation in AVJ
      StartAvjTm(RETROGRADE);         // start a new timer for retr AV delay
      fprintf(env.fpav,               // record AVJ status
        "%.3lf %d %.3lf %.3lf ", env.t, avj.dir, avj.RT, avj.AVDr);
      fprintf(env.fplog,              // log after ActivateAvj() service
        "start AVDr = %.3lf, avj.nR=%d\n", avj.AVDr, avj.nR);     
    }
    else {                            // collide w/ existing ante wave in AVJ
      StopAvjTm(ANTEGRADE);           // stop an ante AVJ activation timer
      avj.nF[avj.phase]++;            // incr AVJ collision counter (phase 4)
      fprintf(env.fpav,               // record AVJ status
        "%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
      fprintf(env.fplog,              // log after ActivateAvj() service
        "collide w/ ante wave, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
    }
  }
  else {                              // simultaneous bidirectional depol
    avj.nF[avj.phase]++;              // incr AVJ collision counter (phase 4)
    fprintf(env.fpav,                 // record AVJ status
      "%.3lf %d %.3lf %d ", env.t, avj.dir, avj.RT, -1);
    fprintf(env.fplog,                // log after ActivateAvj() service
      "bidir depolarization, avj.nF(%d)=%ld\n", avj.phase, avj.nF[avj.phase]);
  }  
}

/*********************************************************************
StartAvjRef: AVJ starts refractory period
*********************************************************************/
void StartAvjRef()
{
  long r;

  // calculate AVJ refractory period
  avj.ref = avj.MinRef+avj.beta*(1-exp(-avj.RT/avj.tau_r));
  // Assume AVJ refractory period does not exceed the AVD 
  if (avj.nA>0 && avj.ref<avj.AVDa) avj.ref = avj.AVDa;  // limited by AVD
  if (avj.nR>0 && avj.ref<avj.AVDr) avj.ref = avj.AVDr;  // limited by AVD
  //srandom(-1000);
  r = (long) random();
  avj.ref = gaussian(avj.ref,avj.Ref_std,ran2(&r),ran2(&r));
  avj.phase = REFRACTORY;             // update AVJ phase to REFRACTORY
  avj.t = 0;                          // start AVJ refractory timer
  
  fprintf(env.fpav,"%.3lf\n", avj.ref); // record AVJ status
  fprintf(env.fplog,                  // log after StartAvjRef() service
    "T=%.3lf 'Start AVJ REF' RT=%.3lf, ref=%.3lf\n", 
    env.t, avj.RT, avj.ref);  
}

/*********************************************************************
StartAvjPh4: AVJ starts phase 4
*********************************************************************/
void StartAvjPh4()
{
  avj.phase = PHASE4;                 // update AVJ phase to PHASE4
  avj.Vm = avj.Vr;                    // set Vm to Vr
  avj.RT0 = env.t;                    // AVJ starts recovery

  fprintf(env.fplog,                  // log after StartAvjPh4() service
    "T=%.3lf 'Start AVJ PH4' Vm=%.2lf, avj.nA=%d, avj.nR=%d\n",
    env.t, avj.Vm, avj.nA, avj.nR);  
}

/*********************************************************************
StartAvjTm: start a new activation timer in AVJ
Input:  dir - direction of the activation
*********************************************************************/
void StartAvjTm(int dir)
{
  if (dir==ANTEGRADE) {              // a new timer for ante AV delay
    avj.tmA[avj.nA][0] = 0;          // start the timer
    avj.tmA[avj.nA][1] = avj.AVDa;   // set ante AV delay time
    avj.nA++;                        // allow multiple ante waves in AVJ
  }
  else {                             // a new timer for retr AV delay
    avj.tmR[avj.nR][0] = 0;          // start the timer
    avj.tmR[avj.nR][1] = avj.AVDr;   // set retr AV delay time
    avj.nR++;                        // allow multiple retr waves in AVJ
  }
}

/*********************************************************************
StopAvjTm: stop an activation timer in AVJ
Input:  dir - direction of the activation
*********************************************************************/
void StopAvjTm(int dir)
{
  int i;

  if (dir==ANTEGRADE) {              // stop an ante AV delay timer
    for (i=0;i<avj.nA-1;i++) {       // update AVJ ante activation queue
      avj.tmA[i][0] = avj.tmA[i+1][0]; // activation timer
      avj.tmA[i][1] = avj.tmA[i+1][1]; // conduction delay
    }
    avj.nA--;                        // decr ante AVJ activation counter
    avj.tmA[avj.nA][0] = -1;         // disable an ante AVJ timer
  }
  else {                             // stop a retr AV delay timer
    for (i=0;i<avj.nR-1;i++) {       // update AVJ retr activation queue
      avj.tmR[i][0] = avj.tmR[i+1][0]; // activation timer
      avj.tmR[i][1] = avj.tmR[i+1][1]; // conduction delay
    }
    avj.nR--;                        // decr retro AVJ activation counter
    avj.tmR[avj.nR][0] = -1;         // disable a retr AVJ timer
  }
}

/*********************************************************************
StartVtrTm: start a new activation timer in ventricle
Input:  dir - direction of the activation
*********************************************************************/
void StartVtrTm(int dir)
{
  if (dir==ANTEGRADE) {              // a new timer for ante vtr delay
    vtr.tmA[vtr.nA] = 0;             // starts an ante activation in vtr
    vtr.nA++;                        // incr ante vtr activation counter
    if (vtr.nA>MAX_VENT_TM) ErrMsg("Too many ante waves in ventricle!",2);
  }
  else {                             // a new timer for retr vtr delay
    vtr.tmR[vtr.nR] = 0;             // starts a retr activation in vtr
    vtr.nR++;                        // incr retr vtr activation counter
    if (vtr.nR>MAX_VENT_TM) ErrMsg("Too many retr waves in ventricle!",2);
  }
}

/*********************************************************************
StopVtrTm: stop an activation timer in ventricle
Input:  dir - direction of the activation
*********************************************************************/
void StopVtrTm(int dir)
{
  int i;

  if (dir==ANTEGRADE) {              // stop an ante vtr delay timer
    for (i=0;i<vtr.nA-1;i++)         // update ante vtr activation queue
      vtr.tmA[i] = vtr.tmA[i+1];
    vtr.nA--;                        // decr number of ante waves in vtr
    vtr.tmA[vtr.nA] = -1;            // stop an ante vtr activation timer
  }
  else {                             // stop a retr vtr delay timer
    for (i=0;i<vtr.nR-1;i++)         // update vtr retr activation queue
      vtr.tmR[i] = vtr.tmR[i+1];
    vtr.nR--;                        // decr retr vtr activation counter
    vtr.tmR[vtr.nR] = -1;            // stop a retr vtr activation timer
  }
}

/*********************************************************************
AnteEscAvj: event service when ante wave escapes AVJ and invade vtr
*********************************************************************/
void AnteEscAvj()
{
  LogEvent(env.fplog,"AnteEscAvj"); 
  
  StopAvjTm(ANTEGRADE);              // stop an ante AVJ activation timer
  if (avj.nA==0) avj.RT0 = env.t;    // AVJ starts recovery

  if (env.t-vtr.tmA0<=vtr.ref) {     // non-capture because vtr refr
    vtr.nB[ANTEGRADE]++;             // will not generate ante wave in vtr
    fprintf(env.fplog,               // log after AnteEscAvj() service
      "avj.nA=%d, vtr.nB(a)=%d(V-Blk)\n", avj.nA, vtr.nB[ANTEGRADE]);
    return;                          // force return
  }                                  // otherwise, capture vtr and ...
  StartVtrTm(ANTEGRADE);             // starts an ante activation in vtr
  vtr.tmA0 = env.t;                  // update last ante captured time
  
  fprintf(env.fplog,                 // log after AnteEscAvj() service
    "avj.nA=%d, vtr.nA=%d\n", avj.nA, vtr.nA);
}

/*********************************************************************
RetrEscAvj: event service when retr wave escapes AVJ and invade atrium
*********************************************************************/
void RetrEscAvj()
{
  double t;

  LogEvent(env.fplog,"RetrEscAvj"); 
  
  atr.nF[avj.phase]++;               // incr atrial fusion counter
  if (atr.NextAA-atr.t > atr.AtrDly) { // invade SA node before next impulse
    t = GetAfIntv(atr.lambda,atr.AAstd,atr.AA_MODEL); // reset SA node
    fprintf(env.fpaa,"%lf\n",t);
    atr.NextAA += t;                 // update arrival time of next AF
  }                                  // otherwise, simply atrial fusion
  StopAvjTm(RETROGRADE);             // stop a retr AVJ activation timer
  if (avj.nR==0) avj.RT0 = env.t;    // AVJ starts recovery

  fprintf(env.fplog,                 // log after RetrEscAvj() service
    "avj.nR=%d, atr.nF(%d)=%ld, nextAA(%.3lf)+%.3lf=%.3lf\n", 
    avj.nR, avj.phase, atr.nF[avj.phase], atr.NextAA-t, t, atr.NextAA);
}

/*********************************************************************
VtrFusion: event service when detecting ventricular fusion
*********************************************************************/
void VtrFusion()
{
  int i;

  LogEvent(env.fplog,"VtrFusion"); 

  vtr.nF[avj.phase]++;               // incr ventricular fusion counter
  StopVtrTm(ANTEGRADE);              // stop an ante vtr activation timer
  StopVtrTm(RETROGRADE);             // stop a retr vtr activation timer
  
  fprintf(env.fplog,                 // log after VtrFusion() service
    "vtr.nF(%d)=%ld, vtr.nA=%d, vtr.nR=%d\n", 
    avj.phase, vtr.nF[avj.phase], vtr.nA, vtr.nR);
}

/*********************************************************************
RegEventV:  register ventricular event (Vs/Vp) and update parameters
Input: event - PACE or SENSE
*********************************************************************/
void RegEventV(int event)
{
  env.beat++;               // increment beat counter
  vele.VpFlg = event;       // flag for sensed/paced beat
  vele.cnt[event]++;        // increment Vs/Vp counter
  vele.RR = vtr.t;          // record RR interval ended with Vs/Vp
  // Get the next ventricular pacing interval
  vele.NextVV = GetVpIntv(vele.RR,vele.BI,vele.NextVV,vele.VP_MODEL);
  vele.t = 0;               // reset Vp clock timer (inhibit Vp)
  vtr.t = 0;                // clear vtr time to measure RR interval
  fprintf(env.fprr,         // output RR interval and Vp/Vs flag
    "%.3lf %d\n",vele.RR,event);
}

/*********************************************************************
VtrSense: event service for ventricular sense
*********************************************************************/
void VtrSense()
{
  int i;

  LogEvent(env.fplog,"VtrSense"); 

  RegEventV(SENSE);         // update parameters on Vs event
  for (i=0;i<vtr.nA-1;i++)  // update ante vtr activation timers
    vtr.tmA[i] = vtr.tmA[i+1];
  vtr.nA--;
  vtr.tmA[vtr.nA] = -1;     // stop an ante vtr activation timer

  fprintf(env.fplog,        // log after VtrSense() service
    "vtr.nA=%d, nextVV=%.3lf, Vs#%ld, Beat#%ld, RR=%.3lf\n", 
    vtr.nA, vele.NextVV, vele.cnt[SENSE], env.beat, vele.RR);
}

/*********************************************************************
VtrPace: event service for ventricular pace
*********************************************************************/
void VtrPace()
{
  LogEvent(env.fplog,"VtrPace"); 

  RegEventV(PACE);          // update parameters on Vp event
  if (env.t-vtr.tmR0<=vtr.ref) { // non-capture because vtr refr
    vtr.nB[RETROGRADE]++;   // will not generate a retr wave in vtr
    return;                 // force return
  }                         // otherwise, Vp capture vtr and ...
  StartVtrTm(RETROGRADE);   // starts a retr activation in vtr
  vtr.tmR0 = env.t;         // update last captured Vp time
  
  fprintf(env.fplog,        // log after VtrPace() service
    " vtr.nR=%d, nextVV=%.3lf, Vp#%ld, Beat#%ld, RR=%.3lf\n", 
    vtr.nR, vele.NextVV, vele.cnt[PACE], env.beat, vele.RR);
}

/*********************************************************************
GetAfIntv: get next Afib interval
*********************************************************************/
double GetAfIntv(double lambda, double sd, int method)
{
 long init;
 long k;
 double next_afib_int, mean, next_afib_rate;

 k = atr.nAf[0]+atr.nAf[1];
 mean = 1/lambda;
 //srandom(-1000);
 init = (long) random();
 switch (method) {
   case 0:      // random value from truncated exponential distribution
//   next_afib_int = expdev(mean,&init);
     next_afib_int = 0.05+expdev(mean-0.05,&init); // limit AA >= 50ms
     break;
   case 1:      // inverse of random value from Poisson distribution
     next_afib_rate = poidev(lambda,&init);
     next_afib_rate++;     // Ensure next_afib_rate>0 (avoid divided by 0)
     next_afib_int = 1.0/next_afib_rate;
     break;
   case 2:      // random value from uniform distribution
     next_afib_int = mean+sd*mean*ran2(&init);
     break;
   case 3:      // random value from Gaussian distribution
     next_afib_int = gaussian(mean,sd*mean,ran2(&init),ran2(&init));
     break;
   case 4:      // atrial pacing protocl with programmed S1S2
     if (k < env.MAX_RR/2) 
       next_afib_int = 1.0/lambda;
     else 
       next_afib_int = atr.S1S2;
     break;
   case 5:      // atrial pacing protocl with programmed S1S2/S2S3
     if (k < env.MAX_RR/2) 
       next_afib_int = 1.0/lambda;
     else  {  
       if (k%2==0)                // even beat
         next_afib_int = atr.S1S2;
       else                       // odd beat
         next_afib_int = atr.S2S3;
     }
     break;
   default:     // fixed rate, atrial flutter protocol
     next_afib_int = 1.0/lambda;
     break;
  }

  if (next_afib_int < 0.05) next_afib_int = 0.05; // ensure interval >= 50ms
  return next_afib_int;
}