/*
The permission is granted to use the software as is and redistribute it in
its original and complete form. If you find this program useful and publish
results obtained by it, please cite the program as:
Michael Whitbeck, Desert Research Institute, 1991, "REACT - Program to solve
kinetic equations for chemical systems". Version 1.00, this code is
available from the authors or from Dr. Whitbeck.
*/

/*
 rxn    chemical reaction simulation
*/

#include 
#include "rxn.h"

extern int regress();
extern void lsoda();
extern void rxnparse();
extern void reprint();
extern void initial();
extern void *derivs();
extern void setfiles();
double rwork1, rwork5, rwork6, rwork7;
double *atol, *rtol, t, tout;
int iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9;
int itol, itask, istate, iopt, jt;
main (argc,argv) int argc; char **argv;
{
double step, stop, savet;
int n,s;
char *help=" you must specify a file name\n\n\n\
 REACT is a chemical kinetics simulation program.\n\
 This code may be freely used and distributed but must not be sold;\n\
 you may make a reasonable charge for handling/copying for others.\n\n\
 There are two input files required. If you invoke the code by typing\n\
	react foo\n\
 then react opens file foo.m to read a reaction mechanism of the form\n\
	k-value	reactants -> products\n\
 It also opens a file foo.p containing the run parameters of the form\n\
	start-time stop-time output-time-interval tolerance\n\
	species-name initial-value abs-tolerance relative-tolerance\n\
\n\
 You may provide an optional keyword fit or FIT to specify the existance\n\
of a file foo.f containing observed data for fitting the model.\
\n\n\n\
Send corresponance and bug reports to :\n\
___________________________________________________________\n\
|Mike Whitbeck             | whitbeck@unssun.unr.edu      |\n\
|Desert Research Inst.     | whitbeck@wheeler.wrc.unr.edu |\n\
|POB 60220                 | whitbeck@sanjuan.UUCP        |\n\
|RENO, NV 89506 USA        | 702-673-7348                 |\n\
|__________________________|______________________________|\n\
the University of Nevada\n\n";

int FIT;

extern void wr_results();
extern void wr_species();

if (argc <2)  bail(help);
iwork1= iwork2= iwork5= iwork6= iwork7= iwork8= iwork9= 0;
rwork1= rwork5= rwork6= rwork7= 0.0;

FIT = 1; /* just model, no fit required */
if(argc >2 ) {
  if (!strcmp(argv[1],"fit")) {FIT = 2;}
  if (!strcmp(argv[1],"FIT")) {FIT = 2;}
  if (!strcmp(argv[2],"fit")) FIT = 0;
  if (!strcmp(argv[2],"FIT")) FIT = 0;
}
if (FIT == 2) {FIT = 0; argv++;}
if (FIT && (argc >2)) bail(help);


setfiles(argv[1]); /* build file names */
rxnparse(&n,&s); /* input mechanism and parse */
rtol = (double *) malloc((unsigned)((s+1)*sizeof(double)));
atol = (double *) malloc((unsigned)((s+1)*sizeof(double)));
atol[0]  = rtol[0] = 0.0;
initial(s,&t,&tout,&stop,&step, atol, rtol); /* set initial values */
for (iopt = 1; iopt <= NS; iopt++) DY[iopt] = Y[iopt];
wr_species();
savet = t;
#ifdef DEBUG
reprint(n,s);
for (iout=1; iout <= s; iout++) printf("%15.4e %15.4e\n",atol[iout],rtol[iout]);
#endif

itol = 4; itask = 1; istate = 1; iopt = 0; jt = 2;
switch(FIT){
	case 0: /* do not fit model to data */
                regress(fitdata);
	case 1: /* fit model to data */
		istate =1;
		initial(s,&t,&tout,&stop,&step, atol, rtol); /* set initial values */
		if ((Ywrite = fopen(results,"w+")) == NULL) 
		  bail ("Could not open results file");
		wr_results(s,t,Y); /* initial values */
		do {
		  lsoda(derivs,s,Y,&t,tout,itol,rtol,atol,itask,&istate,iopt,jt,
			iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
			rwork1,rwork5,rwork6,rwork7);
		  wr_results(s,t,Y);
		  /*  if (istate <= 0) {
		      printf("error istate = %d\n",istate);
		      exit(0);
		      } */
		  switch (istate+7){
		  case 6: bail("too many steps being taken\nGiving up\n"); 
		  case 5: bail("increase tolerance values"); break;
		  case 3: bail("possible singularity");break;
		  case 4:
		  case 2: bail("not converging"); break;
		  case 1: bail("component vanished");break;
		  case 0: bail("code problem"); break;
		  case 8:
		  case 9:
		    default: break;
		  }
		  tout += step;
		} while(t < stop);
		fclose(Ywrite);
		default: break;
	      }
}