/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  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 "predcont.h"
#include "pressure.h"
#include "radius.h"
#include "called.h"
#include "wind.h"
#include "hcmap.h"

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

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

	DEBUG_ENTRY( "iter_end_check()" );

	/* >>chng 05 nov 22 - NPA.  Stop calculation when fraction of oxygen frozen
	 * out on grains gets too high - 
	 * NB this test is not used since StopCalc.StopDepleteFrac is set to > 1 */
	oxy_in_grains = 0.0f;
	count_ptr<chem_atom> elOxygen = unresolved_atom_list[ipOXYGEN];
	for(i=0;i<mole_global.num_calc;++i)
	{
		/* define the abundance of oxygen frozen out on grain surfaces */
		if( mole_global.list[i]->nAtom.find(elOxygen) != mole_global.list[i]->nAtom.end() && mole_global.list[i]->parentLabel.empty() )
			oxy_in_grains += (1 - mole_global.list[i]->lgGas_Phase)*mole.species[i].den*mole_global.list[i]->nAtom[elOxygen];
	}
	/*fprintf(ioQQQ, "DEBUG oxy in grains %.2e %e %e\n", 
		oxy_in_grains ,
		oxy_in_grains/MAX2(SMALLFLOAT,dense.gas_phase[ipOXYGEN]) , StopCalc.StopDepleteFrac );*/

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " iter_end_check called, zone %li.\n" , nzone);
	}
	ASSERT( hcmap.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" );
		}
		return( lgEndFun_v );
	}

	/* check whether trace is needed for this zone and iteration */
	if( (nzone >= trace.nznbug && iteration >= trace.npsbug) && trace.lgTrOvrd )
	{
		if( trace.nTrConvg==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.nTrConvg = abs( trace.nTrConvg );
	}

	/* option to turn printout on after certain zone; only the master rank talks */
	if( prt.lgPrtStart && prt.nstart == nzone )
	{
		called.lgTalk = cpu.i().lgMPI_talk();
	}

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

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

	/* set temperature floor option  -
	 * go to constant temperature calculation if temperature
	 * falls below floor */
	if( phycon.te < StopCalc.TeFloor )
	{
		thermal.lgTemperatureConstant = true;
		thermal.ConstTemp = (realnum)StopCalc.TeFloor;
		phycon.te = thermal.ConstTemp;
		TempChange(thermal.ConstTemp , false);
		TotalInsanity();
	}

	/* check on radiation pressure - constant pressure unstable if dominated
	 * by radiation pressure */
	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(ioQQQ);
		}
		lgDone = true;
		conv.lgBadStop = true;
		strncpy( StopCalc.chReasonStop, "of radiation pressure.", sizeof(StopCalc.chReasonStop) );
	}

	/* radius and resulting volume too large for this cpu */
	if( radius.drad_x_fillfac*radius.r1r0sq > BIGFLOAT/10.)
	{
		/* too big */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "volume too large for this cpu.", sizeof(StopCalc.chReasonStop) );
	}
	/* supersonic outflowing wind, initial velocity, windv0, was > 0,
	 * but it has coasted to a stop - lgVelPos is false */
	else if( !wind.lgVelPos && wind.lgBallistic() )
	{
		/* wind solution with negative velocity */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "wind veloc too small.", sizeof(StopCalc.chReasonStop) );
	}

	else if( !wind.lgStatic() && fabs(wind.windv) < StopCalc.StopVelocity )
	{
		/* stop if absolute value of velocity falls below value */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "wind V too small.", sizeof(StopCalc.chReasonStop) );
	}

	else if( (StopCalc.nTotalIonizStop>0) && conv.nTotalIoniz> StopCalc.nTotalIonizStop )
	{
		/* stop if exceed number of calls to conv base set with 
		 * stop nTotalIonizStop command */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "nTotalIonizStop reached.", sizeof(StopCalc.chReasonStop) );
	}

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

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

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

	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;
		strncpy( StopCalc.chReasonStop, "mass reached.", sizeof(StopCalc.chReasonStop) );
	}

	/* >>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;
		strncpy( StopCalc.chReasonStop, "sonic point reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( dense.EdenTrue==0 )
	{
		/* calculation failed */
		conv.lgBadStop = true;
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "zero electron density.", sizeof(StopCalc.chReasonStop) );
	}

	else if( radius.lgdR2Small )
	{
		lgDone = true;
		conv.lgBadStop = true;
		strncpy( StopCalc.chReasonStop, "DR small rel to thick.", sizeof(StopCalc.chReasonStop) );
	}

	else if( dense.eden < StopCalc.StopElecDensity )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "lowest EDEN reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( dense.eden/dense.gas_phase[ipHYDROGEN] < StopCalc.StopElecFrac )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "low electron fraction.", sizeof(StopCalc.chReasonStop) );
	}

	/* >>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;
		strncpy( StopCalc.chReasonStop, "freeze out fraction.", sizeof(StopCalc.chReasonStop) );
	}

	/*else if( 2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN] < StopCalc.StopH2MoleFrac )*/
	else if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > StopCalc.StopH2MoleFrac )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "large H_2/H fraction.", sizeof(StopCalc.chReasonStop) );
	}

	else if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 
		StopCalc.StopHPlusFrac )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "low H_+/H fraction.", sizeof(StopCalc.chReasonStop) );
	}

	else if( radius.lgDrMinUsed )
	{
		/* dr became too small */
		conv.lgBadStop = true;
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "DRAD small- set DRMIN.", sizeof(StopCalc.chReasonStop) );
	}

	else if( radius.depth >= radius.StopThickness[iteration-1]/EPS )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "outer radius reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( StopCalc.iptnu >= 0 && 
		opac.TauAbsGeo[0][StopCalc.iptnu-1] >= StopCalc.tauend/EPS )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "optical depth reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( StopCalc.lgStopSpeciesColumn && 
		findspecieslocal( StopCalc.chSpeciesColumn )->column >= StopCalc.col_species/EPS )
	{
		/* StopCalc.col_species default set to COLUMN_INIT == 1e30 */
		lgDone = true;
		sprintf( StopCalc.chReasonStop, "%s column dens reached.", StopCalc.chSpeciesColumn );
		//strncpy( StopCalc.chReasonStop, "H column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( colden.colden[ipCOL_HTOT] >= StopCalc.HColStop/EPS )
	{
		/* StopCalc.HColStop default set to COLUMN_INIT == 1e30 */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "H column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( colden.colden[ipCOL_Hp] >= StopCalc.colpls/EPS )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "H+ column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( (colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]) >= StopCalc.col_h2/EPS )
	{
		/* >>chng 03 apr 15, add molecular hydrogen */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	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;
		strncpy( StopCalc.chReasonStop, "H2+H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( colden.H0_ov_Tspin >= StopCalc.col_H0_ov_Tspin/EPS )
	{
		/* >>chng 05 jan 09, stopping command for N(H0) / Tspin */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "N(H0)/Tspin column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( findspecieslocal("CO")->column >= StopCalc.col_monoxco/EPS )
	{
		/* >>chng 03 oct 27--Nick Abel, add carbon monoxide */
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "CO column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( colden.colden[ipCOL_H0] >= StopCalc.colnut/EPS )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "H0 column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	/* >>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( 
		( (colden.colden[ipCOL_H0]+2.*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])) >= StopCalc.colnut/EPS)  )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "H0-H0+H2 column dens reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( phycon.te > StopCalc.TempHiStopZone )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "highest Te reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( phycon.te < StopCalc.TempLoStopZone )
	{
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "lowest Te reached.", sizeof(StopCalc.chReasonStop) );
	}

	else if( nzone >= geometry.nend[iteration-1] )
	{
		lgDone = true;
		geometry.lgZoneTrp = true;
		strncpy( StopCalc.chReasonStop, "NZONE reached.", sizeof(StopCalc.chReasonStop) );
	}

	/* option to stop calculation when line intensity ratio reaches certain value,
	 * nstpl is number of stop line commands entered */
	else if( StopCalc.nstpl > 0 || StopCalc.ContIndex.size() > 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 H beta */
			if( LineSv[StopCalc.ipStopLin2[i]].SumLine[StopCalc.nEmergent[i]] > 0. )
			{
				char chString[10];
				if( LineSv[StopCalc.ipStopLin1[i]].SumLine[StopCalc.nEmergent[i]]/
					LineSv[StopCalc.ipStopLin2[i]].SumLine[StopCalc.nEmergent[i]] > 
					StopCalc.stpint[i] )
				{
					lgDone = true;
					sprt_wl( chString , StopCalc.StopLineWl1[i] );
					sprintf( StopCalc.chReasonStop, "line %s %s reached", 
						StopCalc.chStopLabel1[i] , chString );
				}
			}
		}
		/* continuum flux exceeded maximum permitted value */
		for( size_t k=0; k < StopCalc.ContIndex.size(); ++k )
		{
			// there are 4 entries for each wavelength: nFnu, nInu, InwT, InwC
			long ind = t_PredCont::Inst().offset() + 4*StopCalc.ContIndex[k];
			double nFnu_model = LineSv[ind].SumLine[0]*pow(10.,radius.Conv2PrtInten);
			if( nFnu_model >= StopCalc.ContNFnu[k].get("erg/s/cm2") )
			{
				char chTemp[10];
				lgDone = true;
				sprt_wl( chTemp, LineSv[ind].wavelength );
				sprintf( StopCalc.chReasonStop, "flux %s %s reached", 
					 LineSv[ind].chALab, chTemp );
			}
		}
	}

	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;
		strncpy( StopCalc.chReasonStop, "internal error - DRAD.", sizeof(StopCalc.chReasonStop) );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}

	else if( lgAbort )
	{
		/* calculation failed */
		conv.lgBadStop = true;
		lgDone = true;
		strncpy( StopCalc.chReasonStop, "calculation aborted.", sizeof(StopCalc.chReasonStop) );
	}

	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.nTrConvg )
		{
			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( !hcmap.lgMapDone && (hcmap.MapZone == 0 || nzone == hcmap.MapZone) )
	{
		/* print last zone if not already done */
		if( !lgPrinted )
		{
			PrtZone();
		}

		/* say that we are doing a map */
		hcmap.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;
		strncpy( StopCalc.chReasonStop, "MAP command used-stop.", sizeof(StopCalc.chReasonStop) );

		/* >>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;

		/* make really sure that the string contained in StopCalc.chReasonStop is properly terminated */
		StopCalc.chReasonStop[sizeof(StopCalc.chReasonStop)-1] = '\0';

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " iter_end_check returns after map.\n" );
		}
		return( lgEndFun_v );
	}

	if( lgEndFun_v && prt.lgOnlyZone )
	{
		cdEXIT(EXIT_FAILURE);
	}

	/* the string contained in StopCalc.chReasonStop must be properly 
	 * terminated -this can't fail - strlen returns the number of characters
	 * in str, excluding the terminal NULL.*/
	if( strlen( StopCalc.chReasonStop ) >= nCHREASONSTOP-1 )
		TotalInsanity();

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " iter_end_check bottom return.\n" );
	}
	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;
	realnum ratio;

	DEBUG_ENTRY( "dmpary()" );

	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 = (realnum)(thermal.cooling[i]/thermal.ctot);
		if( fabs(ratio) > EPS )
		{
			coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
			  ratio,"DOIT");
		}

		ratio = (realnum)(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");
	return;
}