text: /* c the CCl4 inhalation model
   -------- ex_ccl4model.c -> ex_ccl4model.dll ------
   compile in R with: system("gcc -shared -o ex_ccl4model.dll ex_ccl4model.c")
   or with system("R CMD SHLIB ex_ccl4model.c")
  */
  
  #include <R.h>
  static double parms[21];
  #define BW     parms[0]
  #define QP     parms[1]
  #define QC     parms[2]
  #define VFC    parms[3]
  #define VLC    parms[4]
  #define VMC    parms[5]
  #define QFC    parms[6]
  #define QLC    parms[7]
  #define QMC    parms[8]
  #define PLA    parms[9]
  #define PFA    parms[10]
  #define PMA    parms[11]
  #define PTA    parms[12]
  #define PB     parms[13]
  #define MW     parms[14]
  #define VMAX   parms[15]
  #define KM     parms[16]
  #define CONC   parms[17]
  #define KL     parms[18]
  #define RATS   parms[19]
  #define VCHC   parms[20]
  
  double V[5], P[4], AI0, VTC, Q[4];
  
  #define DOSE   out[0]
  #define MASS   out[1]
  #define CP     out[2]
  
  /* 
  c=======================================================================
  c=======================================================================
  c                      Model initialisation
  c=======================================================================
  c=======================================================================
  
  c=======================================================================
  2c Initialise primary parameter common block
  c=======================================================================
   */
  
  void initccl4(void (* odeparms)(int *, double *))
  {
      void derived();
      
      int N=21;
      odeparms(&N, parms);
      derived();           
  
  }
  
  /*=======================================================================
   In this "event", state variable 1 is increased with 1. DOES NOT WORK...
  =======================================================================
  */
  
  void eventfun(int *n, double *t, double *y)
  {
        y[0] = y[0] + 1;
  }
        
  /*=======================================================================
  c Calculate derived parameters from primary parameters
  c=======================================================================
  */
  
  void derived ()
  {
  //     Fraction viscera (kg/(kg BW))
        VTC = 0.91 - (VLC+VFC+VMC);
  
  //     Net chamber volume
  
        V[0] = VCHC - RATS*BW; 
        V[1] = VMC*BW;
        V[2] = VTC*BW;
        V[3] = VFC*BW;
        V[4] = VLC*BW;
  
  //     Initial amt. in chamber (mg)         
        AI0 = CONC*V[0]*MW/24450.; 
        P[0]  = PMA/PB;
        P[1]  = PTA/PB;
        P[2]  = PFA/PB;
        P[3]  = PLA/PB;
        
        Q[2] = QFC*QC;
        Q[3] = QLC*QC;
        Q[0] = QMC*QC;
      
        Q[1] = QC - (Q[0]+Q[3]+Q[2]);
  
  }
  
  /*=======================================================================
  c The dynamic model
  c=======================================================================
  */
  void derivsccl4 (int *neq, double *t, double *y, double *ydot,
               double *out, int *ip)
  {             
        double vconc[5], tconc[5], CA, CX, RAM;
        int i;
        
        if (ip[0] < 3) error("nout should be at least 3"); 
        
  /*c y = AI, AAM, AT, AF, AL CLT, AM
   where clt = the area under the concentration-time curve in the liver 
         AM  = total amount metabolised
  
       concentrations */
        for (i =0; i<5; i++) {
          tconc[i] = y[i]/V[i];
        }
     
  /*     vconc(1) is conc in mixed venous blood */
        vconc[0] = 0.0;
        for (i = 1; i<5; i++){
         vconc[i] = tconc[i]/P[i-1];
         vconc[0] = vconc[0] + vconc[i]*Q[i-1]/QC ;
        }
        
  /*     CA is conc in arterial blood  */
        CA = (QC * vconc[0] + QP * tconc[0])/ (QC + QP/PB);
  
  /*     Exhaled chemical   */
        CX = CA/PB;
  
  /*     metabolisation rate   */
        RAM = VMAX*vconc[4]/(KM + vconc[4]);
    
  /*      the rate of change    */
        ydot[0] = RATS*QP*(CX - tconc[0]) - KL*y[0];
        for ( i = 1; i<5; i++)
          ydot[i] = Q[i-1]*(CA-vconc[i]);      
  
        ydot[4] = ydot[4] - RAM;
        ydot[5] = tconc[4];
        ydot[6] = RAM;
        
  /*  the mass balance (MASS=AAM+AT+AF+AL+AM), should be constant */
        DOSE =  AI0 - y[0];
        MASS = (y[1]+y[2]+y[3]+y[4]+y[6])*RATS;
        CP   = tconc[0]*24450.0/MW;
  }
  