/*timer time dependent models */
#include <string.h>
#include "cddefines.h"
#include "linesave.h"
#include "pressure.h"
#include "timefc.h"
#include "timed.h"
#include "tepowers.h"
#include "he.h"
#include "h.h"
#include "heat.h"
#include "phycon.h"
#include "tcool.h"
#include "zonecnt.h"
#include "radius.h"
#include "rfield.h"
#include "stopcalc.h"
#include "called.h"
#include "sexp.h"
#include "converge.h"
#include "prtzone.h"
#include "lines.h"
#include "prtfinal.h"
#include "timer.h"

void timer(void)
{
	long int i;

	double c31, 
	  coll, 
	  dion, 
	  fac, 
	  psave, 
	  rad, 
	  rec, 
	  time;
	static double change = 40.;

#	ifdef DEBUG_FUN
	fputs( "<+>timer()\n", debug_fp );
#	endif

	/* time dependent rec-col behind shock with no radiation.
	 * initial conditions set by physical conditions in last zone */
	rfield.nflux = -1;
	time = 0.;
	/* this is the flag that will be used to say time-dep model */
	timed.itime = -1;
	rad = 0.;
	PrtZone();
	ZoneCnt.nzone = 1;
	timefcCom.timefc = 1.;
	heat.power = 0.;
	tcool.totcol = 0.;
	tcool.GrossHeat = 0.;
	heat.htot = TotalPressure()*1.5;
	psave = heat.htot;

	if( called.lgTalk )
	{
		fprintf( ioQQQ, "1                              TIME DEPENDENT MODEL OF ZONE%3ld THERMAL ENERGY=%10.3e\n", 
		  timed.itime, psave );
	}

	/* clear out line array */
	for( i=0; i < LineSave.nsum; i++ )
	{
		LineSv[i].sumlin = 0.;
	}

	/*-------------------------------------------------------------------
	 * choose dt to 1/change of cooling time
	 * continue until TE drops below TEND */
	while( phycon.te > StopCalc.tend )
	{
		timed.dt = (float)(heat.htot/tcool.ctot/change);

		/* hydrogen
		 * soln of rad trans like eqn for ioniz of h */
		coll = 1.03e-11*sexp(1.579e5/phycon.te)*pow(phycon.te,0.652);

		if( phycon.te < 2e4 )
		{
			rec = 2.9e-10/(TePowers.te70*TePowers.te10/TePowers.te03);
		}
		else
		{
			rec = 1.31e-8/(phycon.te*TePowers.te10*TePowers.te03);
		}

		fac = sexp(phycon.eden*(rec+coll)*timed.dt);
		h.hii = (float)(phycon.hden*coll/(coll + rec)*(1. - fac) + h.hii*fac);
		h.hi = (float)phycon.hden - h.hii;

		/* helium + => 0 */
		rec = 9.96e-11*pow(phycon.te,-0.669);
		c31 = 1.074e-5/TePowers.sqrte*exp(-1.15e4/phycon.te);
		he.hei3 = (float)(rec*phycon.eden*he.heiii/(c31*phycon.eden + 1.27e-4));
		dion = he.heii*(1. - exp(-phycon.eden*timed.dt*rec*1.3333));
		he.hei1 = (float)((he.hei1 + dion)*timefcCom.timefc);
		he.hei = (float)((he.hei + dion)*timefcCom.timefc);
		he.heii = (float)((he.heii - dion)*timefcCom.timefc);

		/* helium ++ => + */
		dion = he.heiii*(1. - exp(-phycon.eden*timed.dt*2.44e-9/(TePowers.te70*
		  TePowers.te10)));
		he.heiii = (float)((he.heiii - dion)*timefcCom.timefc);
		he.heii = (float)((he.heii + dion)*timefcCom.timefc);

		/* heavy elements
		 * call carbon
		 * call nitrog
		 * call oxygen
		 * call neon
		 * call magnes
		 * call alumin
		 * call silicn
		 * call sulphr
		 * call clrine
		 * call argoni
		 * call calcm
		 * call iron
		 * call cobalt
		 * call nickel
		 *
		 * get new electron density, find emission spectrum and rATE */

		radius.drad = timed.dt;
		lines();

		radius.depth = time;

		if( (ZoneCnt.nzone/ZoneCnt.nprint)*ZoneCnt.nprint == ZoneCnt.nzone )
			PrtZone();

		ZoneCnt.nzone += 1;
		heat.htot += -tcool.ctot*timed.dt;

		if( heat.htot < 0. )
		{
			if( called.lgTalk )
				{
				fprintf( ioQQQ, " htot is <0, htot, ctot, dt=%12.3e%12.3e%12.3e\n", 
				  heat.htot, tcool.ctot, timed.dt );
				}
			puts( "[Stop in timer]" );
			exit(1);
		}

		phycon.te = phycon.te*(float)(heat.htot/(heat.htot + tcool.ctot*timed.dt));
		rad += tcool.ctot*timed.dt;
		time += timed.dt;

		/* following for density changing with time */
		timefcCom.timefc = 1.;

		/* if( tmpwr.eq.0. ) go to4 */
		if( timed.tmpwr != 0. )
		{
			timefcCom.timefc = (float)pow((timed.tmt1 + timed.dt)/timed.tmt1,timed.tmpwr);
			strcpy( pressure.chCPres, "TIME" );
			PresChng();
			timed.tmt1 += timed.dt;
			heat.htot *= timefcCom.timefc;
		}
	}

	heat.power = psave;
	tcool.totcol = rad;
	PrtZone();
	PrtFinal();


#	ifdef DEBUG_FUN
	fputs( " <->timer()\n", debug_fp );
#	endif
	return;
}