/*******************************************************************
  pbpk2cpp.cpp

  Revision 1.0

  initial: 020/04/08


********************************************************************/
/* Copyright (C) 2008  Doug Neubauer
*
* 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.
*
* 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.
*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <conio.h>

const char *precisionName = "double";
const  char *expectedVersion = "1.0";
       
const  char *eulersName  =  "eulers";
const  char *rungeName  =  "rungekutta4";
   
#define  INPUTBUFFERSIZE    (8192)

typedef  struct  _PlotSubstance{
    char *substance ;
    struct _PlotSubstance  *next;
    int  bloodSubstance;
} PlotSubstance ;           

          
typedef  struct  _CompartmentSubstance{
    char *substance ;  
    double  partitionCoefficient;
    double  clearance ;
    char *massBalance ;
    struct _CompartmentSubstance  *next;
} CompartmentSubstance ;               
          

typedef enum {
    LEVEL0,
    LEVELCONSTANTS,
    LEVELSUBSTANCE,
    LEVELCOMPARTMENT,
    LEVELCOMPARTMENTSUBSTANCE,
    LEVELDOSAGE,
    LEVELDOSAGE1,
    LEVELPLOTTING,
}  ;


class  Substance {
    public:  // accessible by all
        Substance(){}
        static  void  AddSubstance(char *name1, char *units1)
        {
            Substance *pointer = new  Substance() ;
            pointer->name = name1 ;
            pointer->quantityUnits = units1;
            pointer->next = substanceHead ;
            substanceHead = pointer;
        }
        static  Substance *ReadSubstanceHead(){return(substanceHead) ;}  
        char *ReadName(){return(name) ;}
        char *ReadQuantityUnits(){return(quantityUnits) ;}
        Substance *ReadNext(){return(next) ;}  
        static  char *FindSubstanceUnits(char *substanceName)  
        {
            Substance *pointer;
            for (pointer = substanceHead; pointer; pointer = pointer->next) {
                if (stricmp(pointer->name, substanceName) == 0) {
                    return(pointer->quantityUnits) ;
                } /* if */     
            } /* forLoop */     
            return("units?? (substance not found)") ;
        }       

    protected:  // accessible by subclass

    private:  // accessible by this class only
        char *name;
        char *quantityUnits;
        Substance *next;
        static  Substance  *substanceHead;
};

Substance *Substance::substanceHead = 0;




                         
class  Dosage {
    public:  // accessible by all
        Dosage(){}
        static  void  AddDosage(char *name1,char *substance1, char *compartment1,double  amount1,double  startTime1,char    *type1)
        {
            Dosage *pointer = new  Dosage() ;
            pointer->name = name1 ;
            pointer->substance = substance1;
            pointer->compartment = compartment1;
            pointer->amount = amount1;
            pointer->startTime = startTime1;
            pointer->type = type1;
            pointer->next = dosageHead ;
            dosageHead = pointer;
        }

        static  void  PrintDosageReport( FILE *filePointer, char *volumeUnits,char *timeUnits)
        {                 
            Dosage *pointer;
            fprintf(filePointer,"<table border=\"1\" cellpadding=\"2\">\n");
            fprintf(filePointer,"<tr><th>Dosage Name</th><th>Substance</th><th>Compartment</th><th>Amount</th><th>startTime</th><th>Type</th></tr>\n");
            for (pointer = dosageHead; pointer; pointer = pointer->next) {
                fprintf(filePointer,"<tr><td>%s</td><td>%s</td><td>%s</td><td>%g %s</td><td>%g</td><td>%s</td></tr>\n",pointer->name,pointer->substance,pointer->compartment,pointer->amount, Substance::FindSubstanceUnits(pointer->substance)  ,pointer->startTime,pointer->type);
            } /* forLoop */
            fprintf(filePointer,"</table>\n");
        }         

        static  Dosage *ReadDosageHead(){return(dosageHead) ;}  
        char *ReadSubstance(){return(substance) ;}
        char *ReadCompartment(){return(compartment) ;}
        double  ReadAmount(){return(amount) ;}
        double  ReadStartTime(){return(startTime) ;}
        Dosage *ReadNext(){return(next) ;}

    protected:  // accessible by subclass

    private:  // accessible by this class only
        char *name;
        char *substance;
        char *compartment;
        double  amount;
        double  startTime;
        char    *type;
        Dosage *next;
        static  Dosage  *dosageHead;
};

Dosage *Dosage::dosageHead = 0;


class  Compartment {
    public:  // accessible by all
        Compartment(){}

        char *ReadName(){return(name) ;}
        char *ReadInputCompartment(){return(inputCompartment) ;}
        char *ReadOutputCompartment(){return(outputCompartment) ;}
        double  ReadVolume(){return(volume) ;}
        double  ReadBloodFlow(){return(bloodFlow) ;}
        double  ReadBloodFlowOutput(){return(bloodFlowOutput) ;}
        Compartment *ReadNext(){return(next) ;}
        CompartmentSubstance *ReadCompartmentSubstancePointer(){return(compartmentSubstancePointer) ;} 
        static  Compartment *ReadCompartmentHead(){return(compartmentHead) ;} 

        static  void  AddCompartment(  char *name1, char *inputCompartment1 , char *outputCompartment1 , double  volume1, double  bloodFlow1, double  bloodFlowOutput1, CompartmentSubstance *compartmentSubstancePointer1 )
        {
            Compartment *pointer = new  Compartment() ;
            pointer->name                 = name1 ;
            pointer->inputCompartment     = inputCompartment1 ;
            pointer->outputCompartment    = outputCompartment1 ;
            pointer->volume               = volume1;
            pointer->bloodFlow               = bloodFlow1;
            pointer->bloodFlowOutput               = bloodFlowOutput1;
            pointer->compartmentSubstancePointer = compartmentSubstancePointer1;
            pointer->next = compartmentHead;
            compartmentHead = pointer;
        }

        static  void  PrintPhysiologicalReport( FILE *filePointer, char *volumeUnits,char *timeUnits, double  blood2PlasmaRatio)
        {
            char  buffer[INPUTBUFFERSIZE];
            int  substanceCount,i;
            Compartment *pointer;
            Substance *substancePointer;
            CompartmentSubstance  *compartmentSubstancePointer;
            double  totalVolume, totalFlow, temporary, adjustedTemporary;
            totalFlow = totalVolume = 0.0;  
            fprintf(filePointer,"Physiological Data\n");
            fprintf(filePointer,"<table border=\"1\" cellpadding=\"2\">\n");
            fprintf(filePointer,"<tr><th>Compartment</th><th>Volume (%% of Total)</th><th>Blood flow rate (%% of CO)</th><th>Input Compartment</th><th>Output Compartment</th></tr>\n" );
            for (pointer = compartmentHead; pointer; pointer = pointer->next) {
                buffer[0] = 0;
                if (pointer->bloodFlowOutput > 0.0) {
                    sprintf(buffer,"<br />(Output: %g)",pointer->bloodFlowOutput);
                    pointer->CheckBloodOutput();
                } /* if */     
                fprintf(filePointer,"<tr><td>%s</td><td>%g</td><td>%g %s</td><td>%s</td><td>%s</td></tr>\n",pointer->name,pointer->volume,pointer->bloodFlow,buffer,pointer->inputCompartment,pointer->outputCompartment );
                totalVolume += pointer->volume ;
                totalFlow += pointer->bloodFlow;
            } /* forLoop */                  
            fprintf(filePointer,"<tr><td>Total</td><td>%g</td><td>%g</td><td> - </td><td> - </td></tr>\n",totalVolume,totalFlow);
            fprintf(filePointer,"</table>\n");

            fprintf(filePointer,"<br />\n");
            fprintf(filePointer,"<br />\n");
                                                   
            fprintf(filePointer,"Substance Dependent Physiological Data\n");
            fprintf(filePointer,"<table border=\"1\" cellpadding=\"2\">\n");
            fprintf(filePointer,"<tr><th rowspan=\"2\">Compartment</th>");
            for (substancePointer = Substance::ReadSubstanceHead(), substanceCount = 0; substancePointer; substancePointer = substancePointer->ReadNext(), substanceCount++) {
                fprintf(filePointer,"<th colspan=\"4\">%s (%s)</th>",substancePointer->ReadName(),substancePointer->ReadQuantityUnits());
            } /* forLoop */     
            fprintf(filePointer,"</tr>\n");
            fprintf(filePointer,"<tr>");
            for (i = 0; i < substanceCount; i++) {
                fprintf(filePointer,"<th>partitionCoefficient<br />(Ctissue/Cplasma)</th><th>Adjusted partitionCoefficient<br />(Ctissue/Cbtotal)</th><th>clearance (%s/%s)</th><th>massBalance model</th>",volumeUnits,timeUnits);
            } /* iForLoop */     
            fprintf(filePointer,"</tr>\n");

            for (pointer = compartmentHead; pointer; pointer = pointer->next) {
                fprintf(filePointer,"<tr><td>%s</td>",pointer->name);
                for (substancePointer = Substance::ReadSubstanceHead(); substancePointer; substancePointer = substancePointer->ReadNext()) {
                    compartmentSubstancePointer = pointer->FindSubstance(substancePointer->ReadName());
                    if (compartmentSubstancePointer) {
                        temporary = compartmentSubstancePointer->partitionCoefficient;
                        if (temporary < 0) {
                            temporary = 1.0 ;
                            adjustedTemporary = 1.0;
                        } /* if */     
                        else {
                            adjustedTemporary = temporary * blood2PlasmaRatio;
                        } /* else */            
                        fprintf(filePointer,"<td>%g</td><td>%g</td><td>%g</td><td>%s</td>", temporary,adjustedTemporary, compartmentSubstancePointer->clearance, compartmentSubstancePointer->massBalance);
                    } /* if */     
                    else {
                        fprintf(filePointer,"<td>-</td><td>-</td><td>-</td><td>-</td>");
                    } /* else */                                              
                } /* forLoop */     
                fprintf(filePointer,"</tr>\n");
            } /* forLoop */                  
            fprintf(filePointer,"</table>\n");              
            fprintf(filePointer,"note: Adjusted partitionCoefficient = partitionCoefficient * blood2PlasmaRatio\n<br />\n");              
        }                     
        CompartmentSubstance *FindSubstance(char *substanceName)
        {   
            CompartmentSubstance  *pointer;
            for (pointer = compartmentSubstancePointer; pointer; pointer = pointer->next) {
                if (strcmp(pointer->substance, substanceName) == 0) {
                    return(pointer) ;
                } /* if */     
            } /* forLoop */     
            return(0) ;
        }                     

    protected:  // accessible by subclass

    private:  // accessible by this class only
        void  Compartment::CheckBloodOutput ()
        {
            double  total;
            Compartment *pointer;
        
            total = bloodFlow;
            for (pointer = compartmentHead; pointer; pointer = pointer->next) {
                if (stricmp(name, pointer->outputCompartment) == 0) {
                    total += pointer->bloodFlow;
                } /* if */                                    
            } /* forLoop */                  
            if (bloodFlowOutput != total) {
                fprintf(stderr,"CheckBloodOutput: %s  expected:%g  found:%g\n",name,bloodFlowOutput,total);
                exit(1);
            } /* if */     
        }    

        char *name;
        char *inputCompartment ;
        char *outputCompartment ;
        double  volume;
        double  bloodFlow;
        double  bloodFlowOutput;
        Compartment  *next;
        CompartmentSubstance  *compartmentSubstancePointer; 
        static  Compartment  *compartmentHead;
};


                
Compartment  *Compartment::compartmentHead;



       

char *AddText( char *s )
{
    char *rt = new char[strlen(s) + 1];
    if (rt == 0)
        exit(1);
    strcpy(rt,s);
    return(rt);
}


class  World {
    public:  // accessible by all
        World(){InitializeWorld();}
        void  World::ProcessInput ();
        void  PrintWorld(); 
        void  World::PrintReport (char *filename);
        void  World::PrintPlotHeader (char *filename, char *outputFilename);
        void  World::GenerateSimulationCore ();
    protected:  // accessible by subclass

    private:  // accessible by this class only
        char *ExtractPlotCompartment(char *text)
        {            
            static  char buffer[INPUTBUFFERSIZE];
            char *s;
            strcpy(buffer, text) ;
            s = strchr(buffer,'.') ;
            if (s == 0) {
                fprintf(stderr,"ExtractPlotCompartment: missing '.' [%s]\n",text);
                exit(1);
            } /* if */     
            *s = 0 ;
            return(buffer) ;
        }                   
        char *ExtractPlotSubstance(char *text)
        {            
            static  char buffer[INPUTBUFFERSIZE];
            char *s;
            strcpy(buffer, text) ;
            s = strchr(buffer,'.') ;
            if (s == 0) {
                fprintf(stderr,"ExtractPlotSubstance: missing '.' [%s]\n",text);
                exit(1);
            } /* if */     
            return(s + 1) ;      
        }                   
        void  World::GeneratePhysiologyCore ();
        void  World::GenerateRungeKutta ();
        void  World::GenerateEulers ();
        void  World::GenerateInitialization ();
        void  World::GenerateSubstanceFunctions ();
        void  World::InitializeWorld ();
        char *World::ExtractParameter (char *inputBuffer);
        char *World::FindParameter (char *inputBuffer, char *parameter);
        int  lineNumber;

        double  blood2PlasmaRatio;
        double  totalVolume;
        char *volumeUnits;
        char *timeUnits;
        double  totalFlow;
        double   startTime;
        double  stopTime;
        double  dTime;
        char *numericalIntegration;

        double  plot_plotTime  ;
        char *plot_plotTitle ;
        char *plot_yAxis     ; 
        char *plot_yMinimum  ; 
        char *plot_yMaximum  ;
        char *plot_yLog      ;
        char *plot_xLog      ;
        char *plot_xAxis     ;
        char *plot_experimentalFile ;
        char *plot_experimentalTitle ;
        double  plot_xScale ;
        PlotSubstance  *plotSubstancePointer;  
};


void  World::InitializeWorld ()
{
    blood2PlasmaRatio = -1.0;
    totalVolume = 0.0;
    volumeUnits = "?";
    timeUnits = "?";
    totalFlow = 0.0;
    startTime = 0.0;
    stopTime = 0.0;
    dTime = 1.0;
    numericalIntegration = "?";
    lineNumber = 0;       
    plot_plotTime          = 0;
    plot_plotTitle         = "?";
    plot_yAxis             = "?";
    plot_yMinimum          = "?";
    plot_yMaximum          = "?";
    plot_yLog              = "?";
    plot_xLog              = "?";
    plot_xAxis             = "?";
    plotSubstancePointer   = 0;        
    plot_experimentalFile = "?" ;
    plot_experimentalTitle = "?" ;
    plot_xScale = 1.0 ;
}




void  World::PrintReport (char *filename)
{
    FILE *filePointer ;
    struct tm *newtime ;
    time_t aclock;
    
    filePointer = fopen(filename,"w") ;
    if (filePointer == 0) {
        fprintf(stderr,"can't open: %s\n",filename);
        exit(1);
    } /* if */     
          
    fprintf(filePointer,"\n");
    fprintf(filePointer,"\n<html>\n");
    fprintf(filePointer,"<head>\n");
    fprintf(filePointer,"</head>\n");
    fprintf(filePointer,"<body>\n");

    time(&aclock);
    newtime = localtime(&aclock);
    fprintf(filePointer,"<br />Report generated by pbpk2cpp.cpp on %s\n", asctime(newtime) );
    fprintf(filePointer,"<br/>\n");
    fprintf(filePointer,"<br/>\n");

    fprintf(filePointer,"<table border=\"1\" cellpadding=\"2\">\n");
    fprintf(filePointer,"<tr><th>Parameter</th><th>Value</th></tr>\n" );
    fprintf(filePointer,"<tr><td>blood2PlasmaRatio    </td><td>%g </td></tr>\n", blood2PlasmaRatio     );
    fprintf(filePointer,"<tr><td>totalVolume          </td><td>%g  %s</td></tr>\n", totalVolume  ,volumeUnits         );
    fprintf(filePointer,"<tr><td>volumeUnits          </td><td>%s </td></tr>\n", volumeUnits           );
    fprintf(filePointer,"<tr><td>timeUnits            </td><td>%s </td></tr>\n", timeUnits           );
    fprintf(filePointer,"<tr><td>totalFlow (cardiac output, CO)           </td><td>%g %s/%s</td></tr>\n", totalFlow, volumeUnits,timeUnits            );
    fprintf(filePointer,"<tr><td>startTime            </td><td>%g %s</td></tr>\n", startTime     ,timeUnits        );
    fprintf(filePointer,"<tr><td>stopTime             </td><td>%g %s</td></tr>\n", stopTime      ,timeUnits        );
    fprintf(filePointer,"<tr><td>dTime                </td><td>%g %s</td></tr>\n", dTime         ,timeUnits        );
    fprintf(filePointer,"<tr><td>numericalIntegration </td><td>%s </td></tr>\n", numericalIntegration  );
    fprintf(filePointer,"</table>\n");
    fprintf(filePointer,"<br/>\n");
    fprintf(filePointer,"<br/>\n");

    Compartment::PrintPhysiologicalReport(filePointer,volumeUnits,timeUnits,blood2PlasmaRatio);

    fprintf(filePointer,"<br />Dosing Schedule\n" );
    Dosage::PrintDosageReport(filePointer,volumeUnits,timeUnits);

    fprintf(filePointer,"</body>\n");
    fprintf(filePointer,"</html>\n\n");


    fclose(filePointer);

}    



char *World::ExtractParameter (char *inputBuffer)
{
    static  char buffer[INPUTBUFFERSIZE];
    char *s;
    char *s1;

    strcpy(buffer, inputBuffer);

    s1 = strchr(buffer,';') ;
    if (s1 == 0) {
        fprintf(stdout,"missing ';' at line %d\n",lineNumber);
        return(0) ;
    } /* if */

    s = buffer ;
    while ((*s == ' ') || (*s == '\t')) {  // scan white
        s++;
    } /* while */

    while ((*s1 == ' ') || (*s1 == '\t') || (*s1 == ';')) {  // scan white
        s1--;
        if (s1 < s) {
            fprintf(stdout,"missing parameter at line %d\n",lineNumber);
            return(0) ;
        } /* if */
    } /* while */
    s1++ ;
    *s1 = 0 ;
    return(s) ;
}



char *World::FindParameter (char *inputBuffer, char *parameter)
{
    int  length = strlen(parameter) ;
    char *s;

    if (strnicmp(inputBuffer, parameter, length) == 0) {
        s = ExtractParameter(inputBuffer + length) ;
        return(s) ;
    } /* if */
    return(0) ;
}



void  World::ProcessInput ()
{
    char inputBuffer[INPUTBUFFERSIZE];      
    char *s;
    char *versionNumber = "?";
    int  stateValue = LEVEL0;
    char *substanceName ;
    char *substanceUnits;
    char *compartment_name;
    char *compartment_inputCompartment ;
    char *compartment_outputCompartment ;
    char *compartment_substanceName;
    double  compartment_partitionCoefficient;
    double  compartment_volume;
    double  compartment_bloodFlow;
    double  compartment_bloodFlowOutput;
    double  compartment_clearance ;
    char *compartment_massBalance ;
    char *dosage_name;
    char *dosage_substance;
    char *dosage_compartment;
    double  dosage_amount;
    double  dosage_startTime;
    char *dosage_type;
    PlotSubstance  *temporaryPointer;
    CompartmentSubstance  *compartmentSubstancePointer; 
    CompartmentSubstance *temporarySubstance;

    while (1) {
        if (gets(inputBuffer) == 0){
            break;
        }
        lineNumber++;
        if (inputBuffer[0] == 0) {
            continue;
        } /* if */
        if (inputBuffer[0] == ' ') {
            continue;
        } /* if */
        if ((inputBuffer[0] == '/') &&  (inputBuffer[1] == '/') ) {
            continue;
        } /* if */

        switch (stateValue) {
        case  LEVEL0:
            if ((s = FindParameter(inputBuffer, "Pbpk2CppVersion:")) !=  0) {
                versionNumber =  AddText( s ) ;
                break;
            } /* if */
            if (strnicmp(inputBuffer, "Constants {", 11) == 0) {
                stateValue = LEVELCONSTANTS;
                break;
            } /* if */
            if (strnicmp(inputBuffer, "Substance {", 11) == 0) {
                stateValue = LEVELSUBSTANCE;
                substanceName = "?" ;
                substanceUnits = "?";
                break;
            } /* if */
            if (strnicmp(inputBuffer, "Compartment {", 13) == 0) {
                stateValue = LEVELCOMPARTMENT;
                compartment_name = "?" ;
                compartment_inputCompartment  = "-" ;
                compartment_outputCompartment  = "-" ;
                compartment_volume = 0 ;
                compartment_bloodFlow = 0 ;
                compartment_bloodFlowOutput = 0 ;
                compartmentSubstancePointer = 0;
                break;
            } /* if */
            if (strnicmp(inputBuffer, "Dosage {", 11) == 0) {
                stateValue = LEVELDOSAGE;
                break;
            } /* if */
            if (strnicmp(inputBuffer, "Plotting {", 10) == 0) {
                stateValue = LEVELPLOTTING;             
                break;
            } /* if */

            fprintf(stdout,"unknown command: [%s] at line %d \n",inputBuffer,lineNumber);
            break;
        case  LEVELCONSTANTS:
            if (inputBuffer[0] == '}') {
                stateValue = LEVEL0;
                break;
            } /* if */

            if ((s = FindParameter(inputBuffer, "blood2PlasmaRatio:")) != 0) {
                blood2PlasmaRatio =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "totalVolume:")) != 0) {
                totalVolume =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "totalFlow:")) != 0) {
                totalFlow =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "startTime:")) != 0) {
                startTime =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "stopTime:")) != 0) {
                stopTime =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "dTime:")) != 0) {
                dTime =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "volumeUnits:")) != 0) {
                volumeUnits =  AddText( s ) ;
                break;
            } /* if */                           
            if ((s = FindParameter(inputBuffer, "timeUnits:")) != 0) {
                timeUnits =  AddText( s ) ;
                break;
            } /* if */                           
            if ((s = FindParameter(inputBuffer, "numericalIntegration:")) != 0) {
                numericalIntegration =  AddText( s ) ;
                break;
            } /* if */
            fprintf(stdout,"unknown constant: [%s] at line %d \n",inputBuffer,lineNumber);
            break;

        case  LEVELSUBSTANCE:
            if (inputBuffer[0] == '}') {
                stateValue = LEVEL0;
                Substance::AddSubstance(substanceName,substanceUnits);
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "name:")) != 0) {
                substanceName =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "units:")) != 0) {
                substanceUnits =  AddText( s ) ;
                break;
            } /* if */
            fprintf(stdout,"unknown substance parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;

        case  LEVELCOMPARTMENT:
            if (inputBuffer[0] == '}') {
                stateValue = LEVEL0;
                Compartment::AddCompartment(compartment_name,compartment_inputCompartment ,compartment_outputCompartment ,  compartment_volume,  compartment_bloodFlow, compartment_bloodFlowOutput, compartmentSubstancePointer   );
                break;
            } /* if */
            if (inputBuffer[0] == '[') {
                stateValue = LEVELCOMPARTMENTSUBSTANCE;
                compartment_substanceName = "?";
                compartment_partitionCoefficient = -1.0 ;
                compartment_clearance  = 0 ;
                compartment_massBalance  = "-" ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "name:")) != 0) {
                compartment_name =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "inputCompartment:")) != 0) {
                compartment_inputCompartment =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "outputCompartment:")) != 0) {
                compartment_outputCompartment =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "volume:")) != 0) {
                compartment_volume =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "bloodFlow:")) != 0) {
                compartment_bloodFlow =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "bloodFlowOutput:")) != 0) {
                compartment_bloodFlowOutput =  atof( s ) ;
                break;
            } /* if */
            fprintf(stdout,"unknown compartment parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;

        case  LEVELCOMPARTMENTSUBSTANCE:
            if (inputBuffer[0] == ']') {
                stateValue = LEVELCOMPARTMENT;
                temporarySubstance = new  CompartmentSubstance() ;
                temporarySubstance->substance            = compartment_substanceName;
                temporarySubstance->partitionCoefficient = compartment_partitionCoefficient ;
                temporarySubstance->clearance            = compartment_clearance ;
                temporarySubstance->massBalance          = compartment_massBalance ;
                temporarySubstance->next = compartmentSubstancePointer ;
                compartmentSubstancePointer = temporarySubstance;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "substance:")) != 0) {
                compartment_substanceName =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "massBalance:")) != 0) {
                compartment_massBalance =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "partitionCoefficient:")) != 0) {
                compartment_partitionCoefficient =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "clearance:")) != 0) {
                compartment_clearance =  atof( s ) ;
                break;
            } /* if */
            fprintf(stdout,"unknown compartment substance parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;
                  
        case  LEVELDOSAGE:
            if (inputBuffer[0] == '}') {
                stateValue = LEVEL0;
                break;
            } /* if */
            if (inputBuffer[0] == '[') {
                stateValue = LEVELDOSAGE1;
                dosage_name = "?" ;
                dosage_substance = "?" ;
                dosage_compartment = "?" ;
                dosage_amount = 0 ;
                dosage_startTime = 0 ;
                dosage_type = "?" ;
                break;
            } /* if */
            fprintf(stdout,"unknown dosage parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;

        case  LEVELDOSAGE1:
            if (inputBuffer[0] == ']') {
                stateValue = LEVELDOSAGE;
                Dosage::AddDosage(dosage_name, dosage_substance, dosage_compartment, dosage_amount, dosage_startTime, dosage_type    );
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "name:")) != 0) {
                dosage_name =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "substance:")) != 0) {
                dosage_substance =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "compartment:")) != 0) {
                dosage_compartment =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "type:")) != 0) {
                dosage_type =  AddText( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "amount:")) != 0) {
                dosage_amount =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "startTime:")) != 0) {
                dosage_startTime =  atof( s ) ;
                break;
            } /* if */
            fprintf(stdout,"unknown dosage parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;
                  
                  
        case  LEVELPLOTTING:
            if (inputBuffer[0] == '}') {
                stateValue = LEVEL0;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "plotTime:")) != 0) {
                plot_plotTime =  atof( s ) ;
                break;
            } /* if */
            if ((s = FindParameter(inputBuffer, "plotTitle:")) != 0) {
                plot_plotTitle =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "xAxis:")) != 0) {
                plot_xAxis =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "experimentalFile:")) != 0) {
                plot_experimentalFile =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "experimentalTitle:")) != 0) {
                plot_experimentalTitle =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "xScale:")) != 0) {
                plot_xScale =  atof( s ) ;
                break;
            } /* if */             

            if ((s = FindParameter(inputBuffer, "yAxis:")) != 0) {
                plot_yAxis =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "yMinimum:")) != 0) {
                plot_yMinimum =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "yMaximum:")) != 0) {
                plot_yMaximum =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "yLog:")) != 0) {
                plot_yLog =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "xLog:")) != 0) {
                plot_xLog =  AddText( s ) ;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "substance:")) != 0) {
                temporaryPointer = new  PlotSubstance() ;
                temporaryPointer->substance = AddText( s ) ;
                temporaryPointer->next = plotSubstancePointer;
                temporaryPointer->bloodSubstance = 0;
                plotSubstancePointer = temporaryPointer;
                break;
            } /* if */             
            if ((s = FindParameter(inputBuffer, "substanceBlood:")) != 0) {
                temporaryPointer = new  PlotSubstance() ;
                temporaryPointer->substance = AddText( s ) ;
                temporaryPointer->next = plotSubstancePointer;
                temporaryPointer->bloodSubstance = 1;
                plotSubstancePointer = temporaryPointer;
                break;
            } /* if */             
            fprintf(stdout,"unknown plotting parameter: [%s] at line %d \n",inputBuffer,lineNumber);
            break;

        default:
            break;
        } /* switch */
    } /* while */     
    
    if (strcmp(versionNumber, expectedVersion) != 0) {
        fprintf(stderr,"FATAL: bad versionNumber: expected:[%s] found:[%s]\n",expectedVersion,versionNumber);
        exit(1);
    } /* if */     

    if ( (stricmp(numericalIntegration, eulersName) != 0) && (stricmp(numericalIntegration, rungeName) != 0) ) {
        fprintf(stderr,"FATAL: bad numericalIntegration: found:[%s]\n",numericalIntegration);
        exit(1);
    } /* if */     
}


void  World::PrintPlotHeader (char *filename, char *outputFilename)
{
    static  char  buffer[INPUTBUFFERSIZE];
    FILE *filePointer ;
    char *s;
    PlotSubstance  *pointer;
    int  j;
    
    filePointer = fopen(filename,"w") ;
    if (filePointer == 0) {
        fprintf(stderr,"can't open: %s\n",filename);
        exit(1);
    } /* if */     
          
    fprintf(filePointer,"set autoscale\n");
    fprintf(filePointer,"set yrange [%s : %s]\n",plot_yMinimum, plot_yMaximum);
    fprintf(filePointer,"set grid\n");
    fprintf(filePointer,"unset log\n");
    fprintf(filePointer,"unset label\n");
    fprintf(filePointer,"set xtics auto\n");
    fprintf(filePointer,"set ytics auto\n");
    fprintf(filePointer,"set mxtics 5\n");
    fprintf(filePointer,"set mytics 5\n");
    if (stricmp(plot_yLog, "yes") == 0) {
        fprintf(filePointer,"set logscale y\n");
    } /* if */     
    if (stricmp(plot_xLog, "yes") == 0) {
        fprintf(filePointer,"set logscale x\n");
    } /* if */     
    
    fprintf(filePointer,"set title \"%s\"\n",plot_plotTitle);
    fprintf(filePointer,"set ylabel \"%s\"\n",plot_yAxis);
    fprintf(filePointer,"set xlabel \"%s\"\n",plot_xAxis);

    fprintf(filePointer,"set terminal png\n");
    sprintf(buffer,"%s.png",outputFilename);        
    fprintf(filePointer,"set output \"%s\"\n",buffer);

    s = "";
    sprintf(buffer,"%s.dat",outputFilename);    
    fprintf(filePointer,"plot ");
    for (j = 0, pointer = plotSubstancePointer; pointer; pointer = pointer->next,j++) {
        fprintf(filePointer,"%s\"%s\" using 1:%d title '%s' with lines",s,buffer,j + 2,pointer->substance);
        s = ", ";
    } /* forLoop */     
    if (plot_experimentalFile[0] != '?') {
        fprintf(filePointer,"%s\"%s\" using 1:2 title '%s' with points 3",s,plot_experimentalFile,plot_experimentalTitle);
    } /* if */     
    fprintf(filePointer,"\n");
    fclose(filePointer);
}    



void  World::GenerateSubstanceFunctions ()
{   
    Compartment *compartmentPointer;
    Compartment *sourceCompartment;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    Substance *substancePointer;
    char *sourceCompartmentName;

    // make substance functions  (one for each substance)
    for (substancePointer = Substance::ReadSubstanceHead(); substancePointer; substancePointer = substancePointer->ReadNext()) {
        substanceName = substancePointer->ReadName();
        printf("void  ProcessPhysiology_%s (%s  currentTime ", substanceName, precisionName);
        for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
            compartmentSubstancePointer = compartmentPointer->FindSubstance(substanceName);
            if (compartmentSubstancePointer) {                                     
                compartmentName =  compartmentPointer->ReadName();                     
                printf(", %s &%s_out, %s %s_in ", precisionName, compartmentName, precisionName, compartmentName);
                if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                    printf(", %s &%s_outEx ", precisionName, compartmentName);
                } /* if */            
            } /* if */     
        } /* forLoop */     
        printf(")\n{\n");

        for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
            compartmentSubstancePointer = compartmentPointer->FindSubstance(substanceName);
            if (compartmentSubstancePointer) {                                     
                compartmentName =  compartmentPointer->ReadName();
                printf("    %s_out = -(%s_%s_tissue2OutputRatio * %s_in) ",  compartmentName, compartmentName, substanceName, compartmentName);
                if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                    printf(" - (%s_%s_tissue2ExcretionRatio * %s_in) ", compartmentName, substanceName , compartmentName); 
                } /* if */            

                for (sourceCompartment = Compartment::ReadCompartmentHead(); sourceCompartment; sourceCompartment = sourceCompartment->ReadNext()) {
                    if (sourceCompartment == compartmentPointer) {
                        continue;
                    } /* if */     
                    sourceCompartmentName = sourceCompartment->ReadName();
                    if (stricmp(compartmentName, sourceCompartment->ReadOutputCompartment()) == 0) {
                        printf(" + (%s_%s_tissue2OutputRatio * %s_in) ",  sourceCompartmentName, substanceName,sourceCompartmentName);
                    } /* if */                                    
                    else if (stricmp(compartmentPointer->ReadInputCompartment(), sourceCompartmentName) == 0) {
                        printf(" + (%g * %s_%s_tissue2OutputRatio * %s_in) ",compartmentPointer->ReadBloodFlow(), sourceCompartmentName,substanceName,sourceCompartmentName);
                    } /* elseIf */     
                } /* forLoop */     
                printf(" ;\n");                                                                      
                if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                    printf("    %s_outEx = (%s_%s_tissue2ExcretionRatio * %s_in) ;\n",  compartmentName, compartmentName, substanceName, compartmentName);
                } /* if */            
            } /* if */     
        } /* forLoop */     
        printf("}\n\n");                       
    } /* forLoop */                  
}    



void  World::GenerateRungeKutta ()
{   
    char  buffer[INPUTBUFFERSIZE];
    Compartment *compartmentPointer;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    Substance *substancePointer;
    int  i;

    printf("void  ProcessPhysiologyTop (%s  currentTime)     // runge kutta 4 (rk4)\n",precisionName);
    printf("{\n");
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();                     
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s  %s_%s_k1, %s_%s_k2, %s_%s_k3,  %s_%s_k4;\n", precisionName, compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s  %s_%s_Ex_k1, %s_%s_Ex_k2, %s_%s_Ex_k3,  %s_%s_Ex_k4;\n", precisionName, compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  
    printf("\n");
    
    for (i = 1; i <= 4; i++) {
        for (substancePointer = Substance::ReadSubstanceHead(); substancePointer; substancePointer = substancePointer->ReadNext()) {
            substanceName = substancePointer->ReadName();
            printf("    ProcessPhysiology_%s ( currentTime ", substanceName );
            for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
                compartmentSubstancePointer = compartmentPointer->FindSubstance(substanceName);
                if (compartmentSubstancePointer) {                                     
                    compartmentName =  compartmentPointer->ReadName();                     
                    buffer[0] = 0 ;
                    if ((i == 2) || (i == 3)) {              
                        sprintf(buffer," + (%s_%s_k%d * dTime * 0.5)",compartmentName,substanceName,i-1);
                    } /* if */     
                    if (i == 4) {
                        sprintf(buffer," + (%s_%s_k%d * dTime)",compartmentName,substanceName,i-1);
                    } /* if */     
                    printf(", %s_%s_k%d, (%s_%s_currentQuantity %s) ", compartmentName,substanceName,i, compartmentName, substanceName,buffer);
                    if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                        printf(", %s_%s_Ex_k%d ", compartmentName,substanceName,i);
                    } /* if */            
                } /* if */     
            } /* forLoop */     
            printf(");\n");
        } /* forLoop */                  
    } /* iForLoop */     
    printf("\n");
    
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s_%s_currentQuantity += (%s_%s_k1  + 2.0 * %s_%s_k2 + 2.0 * %s_%s_k3 + %s_%s_k4) * dTime / 6.0 ;\n", compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s_%s_excretedQuantity += (%s_%s_Ex_k1  + 2.0 * %s_%s_Ex_k2 + 2.0 * %s_%s_Ex_k3 + %s_%s_Ex_k4) * dTime / 6.0 ;\n", compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  

    printf("}\n");                       
}    


void  World::GenerateEulers ()
{   
    Compartment *compartmentPointer;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    
    printf("void  ProcessPhysiologyTop (%s  currentTime)  // Eulers\n",precisionName);
    printf("{\n");

    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s  %s_%s_dQuantity = 0.0;\n",precisionName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s  %s_%s_dExcretion = 0.0 ; \n",precisionName, compartmentName, substanceName);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  
    printf("\n");                                                                      

    printf("    ProcessPhysiology_%s ( currentTime ", substanceName );
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        compartmentSubstancePointer = compartmentPointer->FindSubstance(substanceName);
        if (compartmentSubstancePointer) {                                     
            compartmentName =  compartmentPointer->ReadName();                     
            printf(", %s_%s_dQuantity, (%s_%s_currentQuantity) ", compartmentName,substanceName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf(", %s_%s_dExcretion ", compartmentName,substanceName);
            } /* if */            
        } /* if */     
    } /* forLoop */     
    printf(");\n\n");
    
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s_%s_currentQuantity += %s_%s_dQuantity * dTime ;\n", compartmentName, substanceName , compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s_%s_excretedQuantity += %s_%s_dExcretion * dTime ;\n", compartmentName, substanceName , compartmentName, substanceName);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  

    printf("}\n");                       
}     


#if 0
void  World::GenerateEulers ()
{   
    Compartment *compartmentPointer;
    Compartment *sourceCompartment;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    
    printf("void  ProcessPhysiologyTop (%s  currentTime)\n",precisionName);
    printf("{\n");

    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s  %s_%s_dQuantity = 0.0;\n",precisionName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s  %s_%s_dExcretion = 0.0 ; \n",precisionName, compartmentName, substanceName);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  
    printf("\n");                                                                      

    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s_%s_dQuantity = dTime * %s_%s_tissue2OutputRatio * %s_%s_currentQuantity ;\n", compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s_%s_dExcretion = dTime * %s_%s_tissue2ExcretionRatio * %s_%s_currentQuantity ; \n", compartmentName, substanceName , compartmentName, substanceName, compartmentName, substanceName); 
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  

    printf("\n");                                                                      
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    %s_%s_currentQuantity = %s_%s_currentQuantity - %s_%s_dQuantity ",  compartmentName, substanceName, compartmentName, substanceName, compartmentName, substanceName);
            
            for (sourceCompartment = Compartment::ReadCompartmentHead(); sourceCompartment; sourceCompartment = sourceCompartment->ReadNext()) {
                if (sourceCompartment == compartmentPointer) {
                    continue;
                } /* if */     
                if (stricmp(compartmentName, sourceCompartment->ReadOutputCompartment()) == 0) {
                    printf(" + %s_%s_dQuantity ",  sourceCompartment->ReadName(), substanceName);
                } /* if */                                    
                else if (stricmp(compartmentPointer->ReadInputCompartment(), sourceCompartment->ReadName()) == 0) {
                    printf(" + (%g * %s_%s_dQuantity) ",compartmentPointer->ReadBloodFlow(),  sourceCompartment->ReadName(), substanceName);
                } /* elseIf */     
            } /* forLoop */     
                                                                                        
            printf(" ;\n");                                                                      

            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s_%s_currentQuantity -= %s_%s_dExcretion ; \n",  compartmentName, substanceName, compartmentName, substanceName); 
                printf("    %s_%s_excretedQuantity += %s_%s_dExcretion ; \n",  compartmentName, substanceName, compartmentName, substanceName); 
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  
    printf("}\n");                       
}     
#endif



void  World::GenerateInitialization ()
{   
    Compartment *compartmentPointer;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    double  temporary;
    
    printf("void  InitializeSimulation ()\n");
    printf("{\n");
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            temporary =  compartmentPointer->ReadBloodFlowOutput();
            if (temporary == 0) {
                temporary =  compartmentPointer->ReadBloodFlow();
            } /* if */     

            if (compartmentSubstancePointer->partitionCoefficient < 0) {
                printf("    %s_%s_tissue2OutputRatio = totalFlow * %g / %s_volume ;\n", compartmentName, substanceName, temporary, compartmentName);
            } /* if */     
            else {
                printf("    %s_%s_tissue2OutputRatio = totalFlow * %g / (%s_volume * %g / blood2PlasmaRatio) ;\n", compartmentName, substanceName, temporary, compartmentName, compartmentSubstancePointer->partitionCoefficient);
            } /* else */     

            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    %s_%s_tissue2ExcretionRatio = %g / (%s_volume * %g ); \n", compartmentName, substanceName, compartmentSubstancePointer->clearance , compartmentName, compartmentSubstancePointer->partitionCoefficient ); 
                printf("    printf(\"# %s_%s_tissue2ExcretionRatio:  %%g \\n\", %s_%s_tissue2OutputRatio); \n", compartmentName, substanceName,compartmentName, substanceName);
            } /* if */            
            printf("    printf(\"# %s_%s_tissue2OutputRatio:  %%g \\n\", %s_%s_tissue2OutputRatio); \n", compartmentName, substanceName,compartmentName, substanceName);
        } /* forLoop */     
    } /* forLoop */                  
    printf("}\n");
}    



void  World::GeneratePhysiologyCore ()
{
    printf("\n");                        
    if (stricmp(numericalIntegration, rungeName) == 0) {
        GenerateSubstanceFunctions();
        GenerateRungeKutta();
    } /* if */     
    else {
        GenerateSubstanceFunctions();
        GenerateEulers(); 
    } /* else */     
    printf("\n");   
    GenerateInitialization();
}    


void  World::GenerateSimulationCore ()
{
    struct tm *newtime ;
    time_t aclock;
    Compartment *compartmentPointer;
    CompartmentSubstance  *compartmentSubstancePointer;
    char *compartmentName ;
    char *substanceName;
    Dosage *dosagePointer;
    int  first;
    PlotSubstance *plotPointer;

    time(&aclock);
    newtime = localtime(&aclock);
    printf("\n");
    printf("// (WARNING: DO NOT EDIT) Machine generated by pbpk2cpp.cpp on %s\n", asctime(newtime) );
    printf("// Title: %s\n",plot_plotTitle);
    printf("\n");
    printf("#include <stdlib.h>\n");  
    printf("#include <stdio.h> \n");
    printf("#include <math.h>  \n");
    printf("#include <string.h>\n");
    printf("#include <time.h>  \n");
    printf("\n");
    
    printf("// constants\n");
    printf("const  %s  blood2PlasmaRatio = %g ;\n", precisionName, blood2PlasmaRatio     );
    printf("const  %s  totalVolume       = %g ;\n", precisionName,  totalVolume    );
    printf("const  %s  totalFlow         = %g ;\n", precisionName,  totalFlow    );
    printf("const  %s  dTime = %g ;\n", precisionName,    dTime  );
    printf("\n");

    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        printf("const %s  %s_volume = (%g * totalVolume);\n",precisionName, compartmentPointer->ReadName(), compartmentPointer->ReadVolume());
    } /* forLoop */                  
    printf("\n");
                
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            printf("%s  %s_%s_currentQuantity = 0.0;\n",precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->substance);
            printf("%s  %s_%s_tissue2OutputRatio = 0.0;\n",precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->substance);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("%s  %s_%s_excretedQuantity = 0.0;\n",precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->substance);
                printf("%s  %s_%s_tissue2ExcretionRatio = 0.0;\n",precisionName, compartmentPointer->ReadName(), compartmentSubstancePointer->substance);
            } /* if */            
        } /* forLoop */     
    } /* forLoop */                  


    GeneratePhysiologyCore();

    printf("\n");                        
    printf("void  InjectSubstance (%s  currentTime)\n",precisionName);
    printf("{\n");
    printf("    static  int  finished = 0;\n");       
    printf("    if (finished) {\n");
    printf("        return;    \n");
    printf("    } /* if */     \n");
    
    first = 1;
    for (dosagePointer = Dosage::ReadDosageHead(); dosagePointer; dosagePointer = dosagePointer->ReadNext()) {
        printf("    if (currentTime >= %g) {\n", dosagePointer->ReadStartTime());
        printf("        %s_%s_currentQuantity += %g;\n",dosagePointer->ReadCompartment(),dosagePointer->ReadSubstance(), dosagePointer->ReadAmount());
        if (first) {
            first = 0;
            printf("        finished = 1;    \n");
        } /* if */     
        printf("    } /* if */     \n");
        printf("\n");
        
    } /* forLoop */     
    printf("}\n");                       

               
    printf("\n");   
    printf("void  SanityCheck ()\n");
    printf("{\n");
    printf("    %s total = 0.0;\n",precisionName);
    for (compartmentPointer = Compartment::ReadCompartmentHead(); compartmentPointer; compartmentPointer = compartmentPointer->ReadNext()) {
        for (compartmentSubstancePointer = compartmentPointer->ReadCompartmentSubstancePointer(); compartmentSubstancePointer; compartmentSubstancePointer = compartmentSubstancePointer->next) {
            compartmentName =  compartmentPointer->ReadName();
            substanceName =  compartmentSubstancePointer->substance;
            printf("    total += %s_%s_currentQuantity ;\n", compartmentName, substanceName);
            if (stricmp(compartmentSubstancePointer->massBalance, "liverMetabolism") == 0) {
                printf("    total += %s_%s_excretedQuantity ; \n", compartmentName, substanceName);
            } /* if */            
        } /* forLoop */        
    } /* forLoop */                  
    printf("    printf(\"# SanityCheck: %%g\\n\",total);\n");
    printf("}\n");

    printf("\n");                        
    printf("void  PlotPhysiology (%s  currentTime)\n",precisionName);
    printf("{\n");
    printf("    static  double  nextPlot = 0.0;\n");       
                            
    printf("    if (currentTime >= nextPlot) {\n" );
    printf("        nextPlot += %g;    \n",plot_plotTime);
    
    if (plot_xScale == 1.0) {
        printf("        printf(\"%%g \",currentTime);\n");
    } /* if */     
    else {
        printf("        printf(\"%%g \",currentTime / %g);\n",plot_xScale);
    } /* else */     
    for (plotPointer = plotSubstancePointer; plotPointer; plotPointer = plotPointer->next) {
        if (plotPointer->bloodSubstance) {
            printf("        printf(\"%%g \",%s_%s_currentQuantity / (%s_volume * blood2PlasmaRatio ));   // convert Cbloodtotal to Cplasma\n",ExtractPlotCompartment(plotPointer->substance), ExtractPlotSubstance(plotPointer->substance),ExtractPlotCompartment(plotPointer->substance));
        } /* if */     
        else {
            printf("        printf(\"%%g \",%s_%s_currentQuantity / %s_volume);\n",ExtractPlotCompartment(plotPointer->substance), ExtractPlotSubstance(plotPointer->substance),ExtractPlotCompartment(plotPointer->substance));
        } /* else */     
    } /* forLoop */     
    printf("        printf(\"\\n\");\n");
    printf("        SanityCheck();\n");
    printf("    } /* if */     \n");

    printf("}\n");                       


    printf("\n");
    printf("void  RunSimulation ()\n");
    printf("{\n");
    printf("    %s  currentTime;\n",precisionName);
    printf("\n");
    printf("    InitializeSimulation ();\n");
    printf("    for (currentTime = %g; currentTime < %g; currentTime += dTime) {\n",startTime,stopTime);
    printf("        InjectSubstance(currentTime) ;\n");
    printf("        ProcessPhysiologyTop(currentTime) ;\n");
    printf("        PlotPhysiology(currentTime) ;\n");
    printf("    } /* currentTimeForLoop */ \n");
    printf("}\n");
    
    printf("\n");
    printf("void main ()\n");
    printf("{\n");
    printf("    RunSimulation();\n");
    printf("}\n");
}    


void  test ()
{                           
    double  valueInitial = 10.0;
    double  k = -0.15 ;
    double  dt = 0.5 ;
    double  t, exact, exponential, eulers, exponentialError, eulersError ;
    double  rk,k1,k2,k3,k4, rkError;
    double  exponentialConstant = exp(k * dt) - 1.0;
    double  eulersConstant = k * dt;
    exponential = valueInitial ;
    eulers = valueInitial;
    rk = valueInitial;
    for (t = dt; t < (100.0); t += dt) {
        exact = valueInitial * exp(k * t) ;
        exponential += exponential * exponentialConstant;
        exponentialError = 100.0 * fabs((exponential - exact) / exact);
        eulers += eulers * eulersConstant;                     
        eulersError = 100.0 * fabs((eulers - exact) / exact);

        k1 = k * (rk) ; 
        k2 = k * (rk + 0.5 * dt * k1) ;
        k3 = k * (rk + 0.5 * dt * k2)  ;
        k4 = k * (rk +       dt * k3)   ;
        rk += (k1 + 2.0 * k2 + 2.0 * k3 + k4) * dt / 6.0;
        rkError = 100.0 * fabs((rk - exact) / exact);

        printf("t:%g  exact:%g  exponential:%g (%g%%)  eulers:%g (%g%%)  runge:%g (%g%%)\n",t,exact,exponential,exponentialError,eulers,eulersError,rk,rkError);
        
    } /* forLoop */     
    
}    




int main()
{
// test() ;
// exit(0);

    World  process ;

    process.ProcessInput() ;
    process.PrintReport("report.html");
    process.PrintPlotHeader("data.plt", "calc2");

    process.GenerateSimulationCore();
    return(0) ;
}