/*lgEndFun after each zone by Cloudy, determines whether model is complete */
#include <string.h>
#include "cddefines.h"
/*  */
#ifdef EPS
#	undef EPS
#endif
#define	EPS	1.00001
#include "linesave.h"
#include "lotsco.h"
#include "showme.h"
#include "comneg.h"
#include "drneg.h"
#include "badstp.h"
#include "dump.h"
#include "itercnt.h"
#include "trace.h"
#include "colden.h"
#include "nomole.h"
#include "timed.h"
#include "mappar.h"
#include "phycon.h"
#include "zonecnt.h"
#include "stopcalc.h"
#include "opac.h"
#include "pressure.h"
#include "radius.h"
#include "called.h"
#include "only.h"
#include "pstart.h"
#include "wind.h"
#include "drminu.h"
#include "trovrd.h"
#include "reason.h"
#include "prtlinepres.h"
#include "prtzone.h"
#include "punt.h"
#include "dmpary.h"
#include "timer.h"
#include "lgendfun.h"

int lgEndFun(void)
{
	int lgDone, 
	  lgEndFun_v, 
	  lgPrinted;
	long int i;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " lgEndFun called, zone %li.\n" , ZoneCnt.nzone);
	}

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

	/* check whether trace is needed for this zone and iteration */
	if( (ZoneCnt.nzone >= trace.nznbug && IterCnt.iter >= trace.npsbug) && 
	  trovrd.lgTrOvrd )
	{
		ZoneCnt.nprint = 1;
		trace.lgTrace = TRUE;
	}

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

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

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

	/* >>chng 97 july 23, from last iteration only, to any greater than 1
	 * >>chng 97 july 23, back again to original */
	if( (pressure.lgRadPresON && pressure.pbeta > 1.0) && 
		(strcmp(pressure.chCPres ,"CPRE") == 0) && 
		IterCnt.lgLastIt &&
		pressure.lgRadPresAbortOK )
	{
		/* 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 );

			if( pressure.plines[0]/pressure.RadPres > 0.9 )
			{
				fprintf( ioQQQ, " Caused by Ly-alpha.\n" );
			}
			else
			{
				fprintf( ioQQQ, " Line contributors to radiation pressure follows:\n" );
				PrtLinePres();
			}
		}
		lgDone = TRUE;
		badstp.lgBadStop = TRUE;
		strcpy( reason.chReason, "of radiation pressure" );

	}

	else if( !wind.lgVelPos && wind.windv0 > 0.)
	{
		/* wind solution with negative velocity */
		lgDone = TRUE;
		strcpy( reason.chReason, "wind veloc too small " );
	}

	else if( !wind.lgVelPos && wind.windv0 < 0.)
	{
		/* D-critical solution reached sonic point */
		lgDone = TRUE;
		strcpy( reason.chReason, "sonic point reached  " );
	}

	else if( called.lgBusted )
	{
		/* calculation failed */
		badstp.lgBadStop = TRUE;
		lgDone = TRUE;
		strcpy( reason.chReason, "code returned BUSTED " );
	}

	else if( lotsco.lgLotsCO )
	{
		/* all carbon in CO */
		lgDone = TRUE;
		strcpy( reason.chReason, "carbon fully moleculr" );
	}

	else if( comneg.lgComNeg )
	{
		/* heavy element molecule solution with negative population */
		lgDone = TRUE;
		badstp.lgBadStop = TRUE;
		strcpy( reason.chReason, "negative mole abundan" );
	}

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

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

	else if( phycon.eden/phycon.hden < StopCalc.StopElecFrac )
	{
		lgDone = TRUE;
		strcpy( reason.chReason, "low electron fraction" );
	}

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

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

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

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

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

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

	/* >>chng 9 July 7, when no h2 molecules, include h2 as neutral atomic hydrogen */
	else if( nomole.lgNoH2Mole &&
		( (coldenCom.colden[IPCHI-1]+2.*coldenCom.colden[IPCOLH2-1]) >= 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( ZoneCnt.nzone >= ZoneCnt.nend[IterCnt.iter-1] )
	{
		lgDone = TRUE;
		ZoneCnt.lgZoneTrp = TRUE;
		strcpy( reason.chReason, "NZONE reached.       " );
	}

	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++ )
		{
			if( LineSv[StopCalc.ipStopLin2[i]-1].sumlin > 0. )
			{
				if( LineSv[StopCalc.ipStopLin1[i]-1].sumlin/
					LineSv[StopCalc.ipStopLin2[i]-1].sumlin > StopCalc.stpint[i] )
				{
					lgDone = TRUE;
					sprintf( reason.chReason, "line%5ld reached.   ", 
					  StopCalc.LineStopWl[i] );
				}
			}
		}
	}

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

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

	else
	{
		/* passed all the tests, keep going
		 * time dependent model of this zone? */
		if( timed.itime == ZoneCnt.nzone )
		{
			timer();
			puts( "[Stop in lgendfun]" );
			exit(1);
		}

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

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

	/* do map of cooling function if desired, and not yet done */
	if( !MapPar.lgMapDone && (MapPar.MapZone < 0 || ZoneCnt.nzone == MapPar.MapZone) )
	{
		FILE *ioSAVE;
		/* print last zone if not already done */
		if( !lgPrinted )
		{
			PrtZone();
		}

		/* save old output file then redirect to map file */
		ioSAVE = ioQQQ;
		ioQQQ = ioMAP;
		punt(" map");
		ioQQQ = ioSAVE;

		/* stop after doing map */
		lgEndFun_v = TRUE;
		strcpy( reason.chReason, "MAP command used-stop" );
		
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " lgEndFun returns after map.\n" );
		}
#		ifdef DEBUG_FUN
		fputs( " <->lgEndFun()\n", debug_fp );
#		endif
		return( lgEndFun_v );
	}

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

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