/* 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 */
/*cloudy the main routine, this IS Cloudy, ret 0 normal exit, 1 error exit,
 * called by maincl when used as standalone program */
/*BadStart announce that things are so bad the calculation cannot even start */
#include "cddefines.h"
#include "punch.h"
#include "noexec.h"
#include "lines.h"
#include "abund.h"
#include "continuum.h"
#include "warnings.h"
#include "atmdat.h"
#include "prt.h"
#include "conv.h"
#include "parse.h"
#include "opacity.h"
#include "rt.h"
#include "assertresults.h"
#include "zones.h"
#include "iterations.h"
#include "plot.h"
#include "radius.h"
/* cloudy.h includes definitions of routines only called here */
#include "grid.h"
#include "cloudy.h"

/*BadStart announce that things are so bad the calculation cannot even start */
static void BadStart(void);

/* returns 1 if disaster strikes, 0 if everything appears ok */
int cloudy(void)
{
	int lgOK,
		/* will be used to remember why we stopped, 
		 * return error exit in some cases */
		lgBadEnd;

	/* 
	 * this is the main routine to drive the modules that compute a model
	 * when cloudy is used as a stand-alone program 
	 * the main program (maincl) calls cdInit then cdDrive
	 * this sub is called by cdDrive which returns upon exiting
	 *
	 * this routine uses the following variables:
	 *
	 * nzone 
	 * this is the zone number, and is incremented here
	 * logical function iter_end_check returns a true condition if NZONE reaches
	 * NEND(ITER), the limit to the number of zones in this iteration
	 * nzone is totally controlled in this subroutine
	 *
	 * iteration 
	 * this is the iteration number, it is 1 for the first iteration
	 * and is incremented in this subroutine at the end of each iteration
	 *
	 * iterations.itermx
	 * this is the limiit to the number of iterations, and is entered
	 * by the user
	 * This routine returns when iteration > iterations.itermx
	 */

	/* nzone is zero while initial search for conditions takes place */
	nzone = 0;
	fnzone = 0.;

	/* iteration is iteration number, cal complete when iteration > iterations.itermx */
	iteration = 1;

	/* flag for error exit */
	lgBadEnd = FALSE;

	/* scan in and parse input data */
	ParseCommands();

	/* ContCreateMesh calls fill to set up continuum energy mesh if first call, 
	 * otherwise reset to original mesh.
	 * This is AFTER ParseCommands so that
	 * path and mesh size can be set with commands before mesh is set */
	ContCreateMesh();

	/* create several data sets by allocating memory and reading in data files, 
	 * but only if this is first call */
	atmdat_readin();

	/* fix pointers to ionization edges and frequency array
	 * calls iso_create
	 * this routine only returns if this is a later call of code */
	ContCreatePointers();

	/* fix abundances of elements, in abundances.c */
	AbundancesSet();

	/* set continuum normalization, continuum itself, and ionization stage limits */
	if( ContSetIntensity() )
	{
		/* this happens when disaster strikes in the inital setup of the continuum intensity array,
		 * BadStart is in this file, below */
		BadStart();
		return(1);
	}

	/* print header */
	PrtHeader();

	/* this is an option to stop after printing header only */
	if( noexec.lgNoExec || prt.lgOnlyHead )
		return(0);

	/* guess some outward optical depths, set inward optical depths, 
	 * also calls RT_line_all so escap probs ok before printout of trace */
	RT_tau_init();

	/* generate inital set of opacities, but only if this is the first call 
	 * in this coreload
	 * grains only exist AFTER this routine exits */
	OpacityCreateAll();

	/* print all abundances, gas phase and grains, in abundances.c */
	AbundancesPrt();

	/* this checks that various parts of the code still work properly */
	SanityCheck("begin");

	/* find the initial temperature, punt if initial conditions outside range of code,
	 * abort condition set by flag lgAbort */
	if( ConvInitSolution() )
	{
		BadStart();
		return(1);
	}

	/* set thickness of first zone */
	radius_first();

	/* find thickness of next zone */
	if( radius_next() )
	{
		BadStart();
		return(1);
	}

	/* set up some zone variables, correct continuum for sphercity, 
	 * after return, radius is equal to the inner radius, not outer radius
	 * of this zone */
	ZoneStart("init");

	plot("FIRST");

	/* outer loop is over number of iterations
	 * >>chng 05 mar 12, chng from test on true to not aborting */
	while( !lgAbort )
	{
		IterStart();
		nzone = 0;
		fnzone = 0.;

		/* loop over zones across cloud for this iteration, 
		 * iter_end_check checks whether model is complete and this iteration is done
		 * returns TRUE is current iteration is complete */
		while( !iter_end_check() )
		{
			/* the zone number, 0 during search phase, first zone is 1 */
			++nzone;
			/* this is the zone number plus the number of calls to bottom solvers
			 * from top pressure solver, divided by 100 */
			fnzone = (double)nzone;

			/* use changes in opacity, temperature, ionization, to fix new dr for next zone */
			/* >>chng 03 oct 29, move radius_next up to here from below, so that
			 * precise correct zone thickness will be used in current zone, so that
			 * column density and thickness will be exact 
			 * abort condition is possible */
			if( radius_next() )
				continue;

			/* following sets zone thickness, drad, to drnext */
			/* set variable dealing with new radius, in zones.c */
			ZoneStart("incr");

			/* converge the pressure-temperature-ionization solution for this zone 
			 * NB ignoring return value - should be ok (ret 1 for disaster) */
			ConvPresTempEdenIoniz();

			/* generate diffuse emission from this zone, add to outward & reflected beams */
			RT_diffuse();

			/* do work associated with incrementing this radius, 
			 * total attenuation and dilution of radiation fields takes place here,
			 * reflected continuum incremented here
			 * various mean and extremes of solution are also remembered here */
			radius_increment();

			/* increment optical depths */
			RT_tau_inc();

			/* >>chng 99 may 04, move lines here from after RT_diffuse */
			/* fill in emission line array, adds outward lines */
			/* >>chng 99 dec 29, moved to here from below RT_tau_inc, 
			 * lines adds lines to outward beam,
			 * and these are attenuated in radius_increment */
			lines();

			/* possibly punch out some results from this zone */
			PunchDo("MIDL");

			/* do some things to finish out this zone */
			ZoneEnd();
		}
		/* end loop over zones */

		/* close out this iteration, in startenditer.c */
		IterEnd();

		/* print out some comments, generate warning and cautions*/
		PrtComment();

		GridGather();

		/* punch stuff only needed on completion of this iteration */
		PunchDo("LAST" );

		/* second call to plot routine, to complete plots for this iteration */
		plot("SECND");

		/* print out the results */
		PrtFinal();

		/* this is the normal exit, occurs if we reached limit to number of iterations,
		 * or if code has set busted */
		if( iteration > iterations.itermx || lgAbort )
			break;

		/*ConvIterCheck check whether model has converged or whether more iterations
		 * are needed - implements the iter to converg comnd */
		ConvIterCheck();

		/* reset limiting and initial optical depth variables */
		RT_tau_reset();

		/* increment iteration counter */
		++iteration;

		/* reinitialize some variables to initial conditions at previous first zone
		 * routine in startenditer.c */
		IterRestart();

		/* reset zone number to 0 - done here since this is only routine that sets nzone */
		nzone = 0;
		fnzone = 0.;

		/* now redo escape probabilities */
		RT_line_all(TRUE , FALSE);

		ZoneStart("init");

		/* find new initial temperature, punt if initial conditions outside range of code,
		 * abort condition set by flag lgAbort */
		if( ConvInitSolution() )
		{ 
			lgBadEnd = TRUE;
			break;
		}
	}

	ClosePunchFiles();

	/* this checks that various parts of the code worked properly */
	SanityCheck("final");

	/* check whether any asserts were present and failed.  
	 * return is true if ok, false if not.  routine also checks
	 * number of warnings and returns false if not ok */
	lgOK = lgCheckAsserts(ioQQQ);

	if( lgOK && !warnings.lgWarngs && !lgBadEnd)
	{
		/* no failed asserts or warnings */
		return(0);
	}
	else
	{
		/* there were failed asserts or warnings */
		return(1);
	}

}

/*BadStart announce that things are so bad the calculation cannot even start */
static void BadStart(void)
{
	char chLine[INPUT_LINE_LENGTH];

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

	/* initialize the line saver */
	wcnint();
	sprintf( warnings.chRgcln[0], "   Calculation stopped because initial conditions out of bounds." );
	sprintf( chLine, " W-Calculation could not begin." );
	warnin(chLine);

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

