/*
 * fd	integration of 2nd order ODE for driven oscillator 
 *	in a viscous medium
 *
 *	models driven motion of a protein side chain
 *	see derivs() at end of this file for eqns for xdot and xdotdot
 */

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

extern void     lsoda();
void     derivs();
void     bail();
void     wr_results();
double	 atof();

double		period, zeta, I, Ks, Kp, psimax, omega, value, kcal_erg;
double		atol_value, rtol_value;

int
main(argc, argv)
	int             argc;
	char          **argv;
{
	double         *y;		/* concentrations */
	int             num_reactions;	/* no. reactions read in */
	int             num_species;	/* no. species used */

	double          rwork1, rwork5, rwork6, rwork7;
	double         *atol, *rtol, t, tout, step, stop;
	int             iwork1, iwork2, iwork5, iwork6;
	int             iwork7, iwork8, iwork9;
	int             itol, itask, istate, iopt, jt;
	int             interval_factor_flag;

	int		i;

	Kp = 50;
	period = 8.192;
	zeta = 0.1092561;
	I = 7.69e-15;
	Ks = 27.87461;
	psimax = 40;
	num_reactions = 2;
	num_species = 2;

	if (argc > 1) {
		period = atof(argv[1]);
	}
	if (argc > 2) {
		Kp = atof(argv[2]);
	}
	if (argc > 3) {
		zeta = atof(argv[3]);
	}
	if (argc > 4) {
		Ks = atof(argv[4]);
	}
	if (argc > 5) {
		psimax = atof(argv[5]);
	}

	omega = 2 * M_PI * 1.e12 / period;
	kcal_erg = 4.184e10;
	Ks = Ks * kcal_erg;
	Kp = Kp * kcal_erg;
	psimax = psimax * M_PI / 180;

	interval_factor_flag = 0;
	t = 0;
	step = period / 16 * 1.e-12;
/* 	step = 1.2589254117; */
/* 	tout = period * 1.e-12 + step; */
	tout = period * 1.e-12;
	while (tout < 1.024e-12) tout += period * 1.e-12;
	stop = tout + step * 256;
	atol_value = 1.e-10;
	rtol_value = 1.e-10;

	y = (double *) malloc(	/* now that num_species set... */
			      (size_t) ((num_species + 1) * sizeof(double)));
	rtol = (double *) malloc(
			     (size_t) ((num_species + 1) * sizeof(double)));
	atol = (double *) malloc(
			     (size_t) ((num_species + 1) * sizeof(double)));
	if (y == NULL || atol == NULL || rtol == NULL)
		bail("insufficient memory");


	y[0] = 0.0;	/* never referenced in lsoda()? */
	y[1] = 0.0;
	y[2] = 0.0;
/* 	y[1] = 10.0; */
/* 	y[2] = -1.e16; */
	atol[0] = rtol[0] = 0.0;	/* never referenced in lsoda()? */
	atol[1] = atol[2] = 0.0;	
	rtol[1] = rtol[2] = 0.0;	
	atol[1] = atol_value;
	atol[2] = atol_value;
	rtol[1] = rtol_value;
	rtol[2] = rtol_value;

	/*
	 * set initial values, including y[]
	 */

	iopt= 1;
	itol = 4;
	istate = 1;
	jt = 2;
	itask = 1;
	iwork1 = iwork2 = iwork5 = 0;
	iwork6 = 10000;
	iwork7 = iwork8 = iwork9 = 0;
	rwork1 = rwork5 = rwork6 = rwork7 = 0.0;

	/*
	 * calculate and write results
	 */
/* 	wr_results(t, y, num_species);	*//* initial values */
	i = 0;
	do {
/*		itask = 4;	*//* set of itask=4, rwork1=tout to prevent */
/*		rwork1 = tout;	*//* overshoot of tcrit=rwork1 */
		lsoda(derivs, num_species, y, &t, tout, itol, rtol, atol,
		      itask, &istate, iopt, jt,
		      iwork1, iwork2, iwork5, iwork6, iwork7, iwork8, iwork9,
		      rwork1, rwork5, rwork6, rwork7);
		switch (istate + 7) {
		case 6:
			/* step limit exceeded -- continue trying */
			fprintf(stderr, "istate=%d, t=%g, tout=%g, flag=%d\n",
				istate, t, tout, interval_factor_flag);
/*
 * must keep istate set at 1 
 * lsoda c version frees malloc memory at each exit 
 * control values used by lsoda are static and persist
 * t and y are stored in calling program
 * etc... 
 * so ok to reset istate = 1 for exit before reach tout
 * and ok to change optional input parms (old istate = 3)
 */
			istate = 1;
			break;
		case 5:
			bail("reduce tolerances");
			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:
			wr_results(t, y, num_species);
			/*
			 * increment tout within switch, else have conflict
			 * with case 6
			 */
			istate = 1;
			i += 1;
			if (interval_factor_flag)
				tout *= step;
			else
				tout += step;
			break;
		}
	} while (tout <= stop && i < 256);

	return 0;
}

/*
 * bail
 */

void
bail(msg)
	char           *msg;
{
	extern void     exit();

	fprintf(stderr, "FATAL: %s\n", msg);
	exit(1);
}

/*
 * write results
 */

void
wr_results(t, y, num_species)
	double          t, *y;
	int             num_species;
{
	int             i;

	fprintf(stdout, "%22.14e", t);
	for (i = 1; i <= num_species; i++)
		fprintf(stdout, " %22.14e", y[i]);
	fprintf(stdout, "\n");
}

/*
 * derivs
 */

void
derivs(neq, t, y, ydot)
	double          t, *y, *ydot;
	int             neq;	/* == num_species */
{
	extern double	zeta, I, Ks, Kp, psimax, omega;

	ydot[1] = y[2];
	ydot[2] = - 1/I * (zeta * y[2] + Ks * y[1] + Kp * sin(y[1] - psimax * sin(omega * t)));
}