/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*iter_end_check after each zone by Cloudy, determines whether model is complete */
#include "cddefines.h"
/*  */
#ifdef EPS
#	undef EPS
#endif
#define	EPS	1.00001
#include "lines.h"
#include "mole.h"
#include "conv.h"
#include "rfield.h"
#include "iterations.h"
#include "trace.h"
#include "dense.h"
#include "colden.h"
#include "taulines.h"
#include "hmi.h"
#include "prt.h"
#include "phycon.h"
#include "geometry.h"
#include "stopcalc.h"
#include "opacity.h"
#include "thermal.h"
#include "cooling.h"
#include "pressure.h"
#include "radius.h"
#include "called.h"
#include "wind.h"
#include "reason.h"
#include "map.h"
/* cloudy.h includes definitions of routines only called here */
#include "cloudy.h"

/*dmpary print all coolants for some zone, as from print cooling command */
static void dmpary(void);

int iter_end_check(void)
{
	int lgDone, 
	  lgEndFun_v, 
	  lgPrinted;
	long int i;
	double oxy_in_grains;

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

	/* >>chng 05 nov 22 - NPA.  Stop calculation when fraction of oxygen frozen
	 * out on grains gets too high */
	oxy_in_grains = 0.0f;
	for (i=0;i<NUM_HEAVY_MOLEC;++i)
	{
		/* define the abundance of oxygen frozen out on grain surfaces */
		oxy_in_grains += (1 - co.lgGas_Phase[i])*co.hevmol[i]*co.nOxyg[i]; 
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone);
	}
	ASSERT( map.MapZone >= 00 || !conv.lgSearch );

	/* >>chng 97 jun 09, now called before first zone with nzone 0 */
	if( nzone == 0 )
	{
		lgEndFun_v = FALSE;
		
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " iter_end_check returns, doing nothing since zone 0.\n" );
		}
#		ifdef DEBUG_FUN
		fputs( " <->iter_end_check()\n", debug_fp );
#		endif
		return( lgEndFun_v );
	}

	/* check whether trace is needed for this zone and iteration */
	if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd )
	{
		if( trace.lgTrConvg==0 )
		{
			geometry.nprint = 1;
			trace.lgTrace = TRUE;
		}
		else
			/* trace convergence entered = but with negative flag = make positive,
			 * abs and not mult by -1 since may trigger more than one time */
			trace.lgTrConvg = abs( trace.lgTrConvg );
	}

	/* option to turn printout on after certain zone */
	if( prt.lgPrtStart && prt.nstart == nzone )
	{
		called.lgTalk = TRUE;
	}

	/* check whether model is done, various criteria used. */
	lgDone = FALSE;

	/* this is flag to chekc whether stopping reason was bad */
	conv.lgBadStop = FALSE;

	/* >>chng 01 jul 27, add set temper floor option  -
	 * go to constant temperature calculation if temperature
	 * falls below floor */
	if( phycon.te <= StopCalc.TeFloor )
	{
		thermal.lgTSetOn = TRUE;
		thermal.ConstTemp = StopCalc.TeFloor;
		phycon.te = thermal.ConstTemp;
		tfidle(FALSE);
	}
	/*fprintf(ioQQQ,"lgendfun zone %li colden %e HColStop/EPS %e radius.drad %e\n",
		nzone , colden.colden[ipCOL_HTOT] , StopCalc.HColStop/EPS ,radius.drad);*/

	/* >>chng 97 jul 23, from last iteration only, to any greater than 1
	 * >>chng 97 jul 23, back again to original */
	if( (pressure.lgPres_radiation_ON && pressure.pbeta > 1.0) && 
		(strcmp(dense.chDenseLaw ,"CPRE") == 0) && 
		/* >>chng 03 aug 20, check on extreme values of pbeta, and abort if true */
		(iterations.lgLastIt||(pressure.pbeta>1000.)) &&
		/* >>chng 03 aug 19, add check on pbeta, and abort even if "no abort"
		 * was specified, since rad pres dominated limit can lead to VERY
		 * small H densities and crash due to underflow */
		(pressure.lgRadPresAbortOK||(pressure.pbeta>1000.)) )
	{
		/* const total pres model; if RadPres>PGAS, then unstable, stop */
		if( called.lgTalk )
		{
			fprintf( ioQQQ, "\n STOP since P(rad)/P(gas)=%7.3f >1.0\n", 
			  pressure.pbeta );

			fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" );
			PrtLinePres();
		}
		lgDone = TRUE;
		conv.lgBadStop = TRUE;
		strcpy( reason.chReason, "of radiation pressure" );

	}

	/* supersonic outflowing wind, initial velocity, windv0, was > 0,
	 * but it has coasted to a stop - lgVelPos is false */
	else if( !wind.lgVelPos && wind.windv0 > 0.)
	{
		/* wind solution with negative velocity */
		lgDone = TRUE;
		strcpy( reason.chReason, "wind veloc too small " );
	}

	/* this flag says that 21cm line optical depth is the stop quantity */
	else if( StopCalc.lgStop21cm && (HFLines[0].TauCon >=  StopCalc.tauend) )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "21 cm optical depth  " );
	}

	else if( rfield.extin_mag_V_extended >= StopCalc.AV_extended )
	{
		/* stop at specified AV for (1-g) in scattering opacity */
		lgDone = TRUE;
		strcpy( reason.chReason, "A_V reached.         " );
	}

	else if( rfield.extin_mag_V_point >= StopCalc.AV_point )
	{
		/* stop at specified AV without (1-g) in scattering opacity */
		lgDone = TRUE;
		strcpy( reason.chReason, "A_V reached.         " );
	}

	else if( StopCalc.xMass!=0. &&
		log10(SDIV(dense.xMassTotal))+1.0992099+ 2.*log10(radius.rinner) >= StopCalc.xMass )
	{
		/* stop at specified AV without (1-g) in scattering opacity */
		lgDone = TRUE;
		strcpy( reason.chReason, "mass reached.        " );
	}

	/* >>chg 02 may 31, added pressure.lgSonicPoint logic */
	/* WJH 19 May 2004: for some models, we want to get through the
	 * sonic point and out the other side */ 
	else if( pressure.lgSonicPoint && pressure.lgSonicPointAbortOK )
	{
		/* D-critical solution reached sonic point */
		lgDone = TRUE;
		strcpy( reason.chReason, "sonic point reached  " );
	}

	else if( dense.EdenTrue==0 )
	{
		/* calculation failed */
		conv.lgBadStop = TRUE;
		lgDone = TRUE;
		strcpy( reason.chReason, "zero electron density" );
	}

	else if( radius.lgdR2Small )
	{
		lgDone = TRUE;
		conv.lgBadStop = TRUE;
		strcpy( reason.chReason, "DR small rel to thick" );
	}

	else if( dense.eden < StopCalc.StopElecDensity )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "lowest EDEN reached. " );
	}

	else if( dense.eden/dense.gas_phase[ipHYDROGEN] < StopCalc.StopElecFrac )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "low electron fraction" );
	}

	/* >>chng 05 nov 22, NA add this stop condition - stop when too many molecules
	 * are ices or solids on grains - the limit StopCalc.StopDepleteFrac is changed
	 * with the stop molecular depletion command */
	else if( dense.lgElmtOn[ipOXYGEN] &&
		(oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN])) > StopCalc.StopDepleteFrac )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "freeze out fraction " );
	}

	/*else if( 2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN] < StopCalc.StopH2MoleFrac )*/
	else if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > StopCalc.StopH2MoleFrac )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "low H_2/H fraction" );
	}

	else if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 
		StopCalc.StopHPlusFrac )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "low H_+/H fraction" );
	}

	else if( radius.lgDrMinUsed )
	{
		/* dr became too small */
		conv.lgBadStop = TRUE;
		lgDone = TRUE;
		strcpy( reason.chReason, "DRAD small- set DRMIN" );
	}

	else if( radius.depth >= radius.router[iteration-1]/EPS || radius.lgDrNeg )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "outer radius reached." );
	}

	else if( StopCalc.iptnu >= 0 && 
		opac.TauAbsGeo[0][StopCalc.iptnu-1] >= StopCalc.tauend/EPS )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "optical depth reached" );
	}

	else if( colden.colden[ipCOL_HTOT] >= StopCalc.HColStop/EPS )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) >= StopCalc.col_h2/EPS )
	{
		/* >>chng 03 apr 15, add molecular hydrogen */
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( (2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) + colden.colden[ipCOL_H0]) >= StopCalc.col_h2_nut/EPS )
	{
		/* >>chng 04 feb 10, stopping command for H2 + H I */
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( colden.H0_ov_Tspin >= StopCalc.col_H0_ov_Tspin/EPS )
	{
		/* >>chng 05 jan 09, stopping command for n(H0) / Tspin */
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( co.hevcol[ipCO] >= StopCalc.col_monoxco/EPS )
	{
		/* >>chng 03 oct 27--Nick Abel, add carbon monoxide */
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	/* >>chng 99 jul 7, when no h2 molecules, include h2 as neutral atomic hydrogen,
	 * hmi.lgNoH2Mole is set false in zero, true with command no hydrogen molecules */
	else if( hmi.lgNoH2Mole &&
		( (colden.colden[ipCOL_H0]+2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])) >= StopCalc.colnut/EPS)  )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "column dens reached. " );
	}

	else if( phycon.te > StopCalc.T2High )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "highest Te reached.  " );
	}

	else if( phycon.te < StopCalc.tend )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "lowest Te reached.   " );
	}

	else if( nzone >= geometry.nend[iteration-1] )
	{
		lgDone = TRUE;
		geometry.lgZoneTrp = TRUE;
		strcpy( reason.chReason, "NZONE reached.       " );
	}

	/* option to stop calculation when line intensity ratio reaches certain value,
	 * nstpl is number of stop line commands entered */
	else if( StopCalc.nstpl > 0 )
	{
		/* line ratio exceeded maximum permitted value
		 * do not consider case where norm line has zero intensity */
		for( i=0; i < StopCalc.nstpl; i++ )
		{
			/* the second line is always set to something, default is hbeta */
			if( LineSv[StopCalc.ipStopLin2[i]].sumlin > 0. )
			{
				char chString[10];
				if( LineSv[StopCalc.ipStopLin1[i]].sumlin/
					LineSv[StopCalc.ipStopLin2[i]].sumlin > StopCalc.stpint[i] )
				{
					lgDone = TRUE;
					sprt_wl( chString , StopCalc.StopLineWl[i] );
					sprintf( reason.chReason, "line %s reached.   ", chString );
				}
			}
		}
	}

	else if( radius.drNext <= 0. )
	{
		/* this cant happen */
		if( called.lgTalk )
		{
			fprintf( ioQQQ, " drNext=%10.2e STOP\n", radius.drNext );
		}
		lgDone = TRUE;
		conv.lgBadStop = TRUE;
		strcpy( reason.chReason, "internal error - DRAD" );
		ShowMe();
	}

	else if( lgAbort )
	{
		/* calculation failed */
		conv.lgBadStop = TRUE;
		lgDone = TRUE;
		strcpy( reason.chReason, "calculation aborted  " );
	}

	lgPrinted = FALSE;
	if( lgDone )
	{
		/* flag to call it quits */
		lgEndFun_v = TRUE;
		PrtZone();
		lgPrinted = TRUE;
	}

	else
	{
		/* passed all the tests, keep going */
		/* check whether this zone should be printed */
		if( ((nzone/geometry.nprint)*geometry.nprint == nzone || 
		  nzone == 1) || trace.lgTrConvg )
		{
			PrtZone();
			lgPrinted = TRUE;
		}
		/* flag to keep going */
		lgEndFun_v = FALSE;
	}

	/* dump cooling arrays for this zone? */
	if( prt.nzdump == nzone || prt.nzdump == 0 )
		dmpary();

	/* do map of cooling function if desired, and not yet done */
	/* >>chng 02 may 29, MapZone < = to <= 0 - map 0 did not work */
	/* >>chng 04 jun 16, change to MapZone = 0 for map of first zone then quit,
	 * -1 is not set, positive, do map of that zone */
	if( !map.lgMapDone && (map.MapZone == 0 || nzone == map.MapZone) )
	{
		/* print last zone if not already done */
		if( !lgPrinted )
		{
			PrtZone();
		}

		/* say that we are doing a map */
		map.lgMapBeingDone = TRUE;

		/* save old output file then redirect to map file */
		/* >>chng 01 mar 28, ioMAP may not be initialized, PvH */
		if( ioMAP != NULL )
			map_do(ioMAP, " map");
		else
			map_do(ioQQQ, " map");

		/* stop after doing map */
		lgEndFun_v = TRUE;
		strcpy( reason.chReason, "MAP command used-stop" );

		/* >>chng 03 jun 06, reset iterations since we want to stop even if
		 * iterate xx is specified, bug caught by Joop Schaye */
		iterations.itermx = 0;

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " iter_end_check returns after map.\n" );
		}
#		ifdef DEBUG_FUN
		fputs( " <->iter_end_check()\n", debug_fp );
#		endif
		return( lgEndFun_v );
	}

	if( lgEndFun_v && prt.lgOnlyZone )
	{
		puts( "[Stop in lgendfun since only zone printout desired.]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " iter_end_check bottom return.\n" );
	}
#	ifdef DEBUG_FUN
	fputs( " <->iter_end_check()\n", debug_fp );
#	endif
	return( lgEndFun_v );
}

#ifdef EPS
#	undef EPS
#endif
#define	EPS	0.005
/*dmpary print all coolants for some zone, as from print cooling command */
static void dmpary(void)
{
	long int i;
	float ratio;

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

	fprintf( ioQQQ, 
		" This is a print out of the cooling array for zone number %3ld\n", 
	  nzone );

	fprintf( ioQQQ, 
		" The total heating was HTOT=%10.2e erg/s/cm3, the total cooling was CTOT=%10.2e, and the temperature was%10.3eK.\n", 
	  thermal.htot, thermal.ctot, phycon.te );

	fprintf( ioQQQ, 
		" All coolants greater than%6.2f%% of the total will be printed.\n", 
	  EPS*100. );

	/* flag all significant coolants */
	coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
	for( i=0; i < thermal.ncltot; i++ )
	{
		ratio = (float)(thermal.cooling[i]/thermal.ctot);
		if( fabs(ratio) > EPS )
		{
			coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
			  ratio,"DOIT");
		}

		ratio = (float)(thermal.heatnt[i]/thermal.ctot);
		if( fabs(ratio) > EPS )
		{
			coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
			  ratio,"DOIT");
		}
	}
	coolpr(ioQQQ,"DONE",1,0.,"DONE");

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

