/* 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 */
/*ParseAtomISO parse information from the atom XX-like command line */
#include "cddefines.h"
#include "elementnames.h"
#include "optimize.h"
#include "hydrogenic.h"
#include "input.h"
#include "iso.h"
#include "parser.h"
#include "phycon.h"
#include "rfield.h"
#include "taulines.h"
#include "thirdparty.h"

/*ParseAtomISO parse parameters off the XX-like command */
void ParseAtomISO(long ipISO, Parser &p )
{
	long int numLevels;

	DEBUG_ENTRY( "ParseAtomISO()" );

	/* look for the name of an element - if we don't find one do the entire
	 * iso sequence - returns negative number if element not found */
	long int nelem = p.GetElem( );

	/* H-like Helium is not possible */
	if( ipISO==ipHE_LIKE && nelem==ipHYDROGEN )
	{
		fprintf(ioQQQ," Sorry, H-like He does not exist.\n");
		cdEXIT(EXIT_FAILURE);
	}

	/* collisions - don't pick up the levels collapsed command */
	if( p.nMatch("COLL") && !p.nMatch("LEVE"  ) )
	{
		/* option to turn collisions off, all are
		 * set to 1 in zero.  command can accept only one option at a time */
		if( p.nMatch("EXCI") )
		{
			/* turn off collisional excitation */
			iso_ctrl.lgColl_excite[ipISO] = false;
			phycon.lgPhysOK = false;
		}
		else if( p.nMatch("IONI") )
		{
			/* turn off collisional ionization */
			iso_ctrl.lgColl_ionize[ipISO] = false;
			phycon.lgPhysOK = false;
		}

		else if( p.nMatch("2S2P") || ( p.nMatch("2P2S") && ipISO == ipH_LIKE ) )
		{
			/* >>chng 02 feb 07, change from 2s2p to l-mixing */
			/* this is the atom h-like collision l-mixing command */
			fprintf(ioQQQ,"This command changed to ATOM H-LIKE COLLISIONS L-MIXING\n");
			fprintf(ioQQQ,"I will parse it for now, but may not in the future.\n");
			/* turn off 2s - 2p collisions */
			iso_ctrl.lgColl_l_mixing[ipISO] = false;
			phycon.lgPhysOK = false;
		}

		else if( p.nMatch("L-MI") )
		{
			if( ipISO == ipH_LIKE )
			{
				/* >>chng 02 feb 07, change from 2s2p to l-mixing */
				/* this is the atom h-like collision l-mixing command */
				iso_ctrl.lgColl_l_mixing[ipISO] = false;
				phycon.lgPhysOK = false;
			}
			else if( p.nMatch("THER") )
			{
				/* use l-mix from 
				 *>>refer	l-mix	all Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */
				if( p.nMatch("NO T") )			
				{
					/* This is the "NO Thermal average" command.  It
					 * causes collisions strengths to be evaluated at kT rather than
					 * integrated over a Maxwellian peaked at kT. */
					iso_ctrl.lgCS_therm_ave[ipISO] = false;
				}
				else
				{
					iso_ctrl.lgCS_therm_ave[ipISO] = true;
				}
			}
			else if( p.nMatch("PENG") )
			{
				iso_ctrl.lgCS_Vrinceanu[ipISO] = false;
			}
			else if( p.nMatch(" OFF"  ) )
			{
				/* this is the atom xx-like collision l-mixing command */
				/* turn off same-n collisions */
				iso_ctrl.lgColl_l_mixing[ipISO] = false;
				phycon.lgPhysOK = false;
				iso_ctrl.lgCS_Vrinceanu[ipISO] = false;
				iso_ctrl.lgCS_therm_ave[ipISO] = false;
			}
			else
			{
				fprintf( ioQQQ, " needs parameter\n" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		else if( p.nMatch(" OFF"  ) )
		{
			/* turn everything off, since no keyword given */
			iso_ctrl.lgColl_excite[ipISO] = false;
			iso_ctrl.lgColl_ionize[ipISO] = false;
			iso_ctrl.lgColl_l_mixing[ipISO] = false;
			phycon.lgPhysOK = false;
		}
		else
		{
			fprintf( ioQQQ, " needs parameter\n" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( p.nMatch("CONT") && p.nMatch("LOWE") )
	{
		/* disable continuum lowering for this isoelectronic sequence */
		if( p.nMatch("OFF") )
			iso_ctrl.lgContinuumLoweringEnabled[ipISO] = false;
		else
			iso_ctrl.lgContinuumLoweringEnabled[ipISO] = true;
	}

	else if( p.nMatch("DAMP") )
	{
		if( ipISO == ipHE_LIKE )
		{
			fprintf(ioQQQ," Sorry, the DAMPING option is not implemented for the he-like sequence.\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* turn off absorption due to Ly alpha damping wings */
		hydro.DampOnFac = 0.;
	}

	else if( p.nMatch("DIEL") )
	{
		if( ipISO == ipH_LIKE )
		{
			fprintf(ioQQQ," Sorry, but dielectronic recombination onto the h-like sequence is not possible.\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* This sets which set of data to use for dielectronic recombination.	*/
		if( p.nMatch(" OFF") )
		{
			iso_ctrl.lgDielRecom[ipISO] = false;
		}
		else 
			iso_ctrl.lgDielRecom[ipISO] = true;
	}

	else if( p.nMatch("LEVE") )
	{
		/* the number of levels read in is n, the principal quantum number
		 * only lines with upper levels less than n will be printed */

		/* number of levels for iso-sequence */
		/* there are two options here,
		 * when keyword ELEMENT appears, scan off element name and change levels only for
		 * that one.
		 * when there is no ELEMENT then set all in iso to same number */

		/* lgHydroMalloc is false at start of calculation, set true when space 
		 * allocated for the hydrogen and helium lines.  Once done we must ignore all 
		 * future changes in the number of levels */
		if( p.nMatch("LTE") )
		{
			/* force level ratios to LTE */
			iso_ctrl.lgLTE_levels[ipISO] = true;
		}
		else if( p.nMatch("PRIN") )
		{
			/* only print - do not change levels */
			iso_ctrl.lgPrintNumberOfLevels = true;
		}
		else if( !lgHydroMalloc )
		{
			numLevels = (long int)p.FFmtRead();

			if( !p.lgEOL() )
			{
				if( ipISO == ipH_LIKE && numLevels > NHYDRO_MAX_LEVEL-2 )
				{
					fprintf( ioQQQ, " Not possible to set nhlvl to >NHYDRO_MAX_LEVEL-2= %i\n",
					  NHYDRO_MAX_LEVEL-2 );
					fprintf( ioQQQ, " change NHYDRO_MAX_LEVEL\n");
					cdEXIT(EXIT_FAILURE);
				}

				/* check that alpha transition of highest level is within energy bounds of continuum */
				if( !p.nMatch("COLL") && ipISO == ipH_LIKE &&
					( 2. / POW3((double)numLevels) < rfield.emm ) )
				{
					fprintf( ioQQQ, " Not possible to set iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max to such a high value, since "
						"alpha transition not within energy bounds of code\n");

					fprintf( ioQQQ, " lowest energy is %e and corresponding highest level is %li\n" ,
						rfield.emm, (long)pow(2./rfield.emm, 0.3333) );
					cdEXIT(EXIT_FAILURE);
				}
			}

			if( p.lgEOL() )
			{
				int LevelsResolved=-1 , LevelsCollapsed=10;
				/* no number, so check for either large or small */
				if( p.nMatch("LARG") )
				{
					/* includes all levels with tabulated rec coefficient */
					LevelsResolved = RREC_MAXN;
				}

				/* this is small or compact keyword */
				else if( p.nMatch("SMAL") || p.nMatch("COMP") )
				{
					if( ipISO == ipH_LIKE )
						LevelsResolved = 5;
					else if( ipISO == ipHE_LIKE )
						LevelsResolved = 3;
				}
				else
					/* punch out if no number */
					p.NoNumb("levels");

				if( nelem<0 )
				{
					// element not specified, do entire sequence
					for( nelem=ipISO; nelem<LIMELM; ++nelem )
					{
						iso_sp[ipISO][nelem].nCollapsed_max =
								MIN2( iso_sp[ipISO][nelem].nCollapsed_max , LevelsCollapsed );
						iso_sp[ipISO][nelem].n_HighestResolved_max =
								MIN2( iso_sp[ipISO][nelem].n_HighestResolved_max , LevelsResolved );
						iso_update_num_levels( ipISO, nelem );
					}
				}
				else
				{
					iso_sp[ipISO][nelem].nCollapsed_max = LevelsCollapsed;
					iso_sp[ipISO][nelem].n_HighestResolved_max = LevelsResolved;
					iso_update_num_levels( ipISO, nelem );
				}
			}

			else if( p.nMatch("COLLAP") )
			{
				// set number of collapsed levels
				if( numLevels < 1 )
				{
					fprintf( ioQQQ, "There must be at least one collapsed level.\n");
					cdEXIT(EXIT_FAILURE);
				}

				if( nelem<0 )
				{
					// element not specified, do entire sequence
					for( nelem=ipISO; nelem<LIMELM; ++nelem )
					{
						iso_sp[ipISO][nelem].nCollapsed_max = numLevels;
						iso_update_num_levels( ipISO, nelem );
					}
				}
				else
				{
					iso_sp[ipISO][nelem].nCollapsed_max = numLevels;
					iso_update_num_levels( ipISO, nelem );
				}
			}
			else if( p.nMatch("RESOLV") )
			{
				// number of resolved levels
				if( ( numLevels < 3 ) && !p.nMatch("COLL") )
				{
					fprintf( ioQQQ, " cannot have fewer than 3 resolved levels, the requested number was %li\n" ,
						numLevels  );
					fprintf( ioQQQ, " Sorry.\n" );
					cdEXIT(EXIT_FAILURE);
				}

				if( nelem<0 )
				{
					// element not specified, do entire sequence
					for( nelem=ipISO; nelem<LIMELM; ++nelem )
					{
						iso_sp[ipISO][nelem].n_HighestResolved_max = numLevels;
						iso_update_num_levels( ipISO, nelem );
					}
				}
				else
				{
					iso_sp[ipISO][nelem].n_HighestResolved_max = numLevels;
					iso_update_num_levels( ipISO, nelem );
				}
			}
			else
			{
				fprintf(ioQQQ, "I did not recognize a keyword on this atom xx-like levels command."
					"  Should be COLLAPSED or RESOLVED.\n Sorry.\n\n");
					cdEXIT(EXIT_FAILURE);
			}
		}
	}

	else if( p.nMatch("ERRO") && p.nMatch("GENE"  ) )
	{
		/* Rates will be modified by a randomly generated error that falls within
		 * the range specifically set for each rate (or set of rates).	*/
		iso_ctrl.lgRandErrGen[ipISO] = true;
		iso_ctrl.modelRank[ipISO] = (int)p.FFmtRead();

		iso_ctrl.modelRank[ipISO] = MAX2( 0, iso_ctrl.modelRank[ipISO] );
		if( p.lgEOL() )
			/* Changed to avoid lint complaint */
			/* iso_ctrl.modelRank[ipISO] = (unsigned)time(NULL); */
			iso_ctrl.modelRank[ipISO] = abs((int)time(NULL));

		/* this allows a seed that is dependent upon the processor rank
		 * in a parallel run. */
		/* We add 2 so that the seed is always greater than 1, which would reset the generator. */
		init_genrand( (unsigned)iso_ctrl.modelRank[ipISO] + 2);

		if( p.nMatch("PESS") )
			iso_ctrl.lgPessimisticErrors = true;
		else
			iso_ctrl.lgPessimisticErrors = false;
	}

	else if( p.nMatch(" FSM") )
	{
		if( ipISO == ipH_LIKE )
		{
			fprintf(ioQQQ," Sorry, but fine-structure mixing can only be implemented for the He-like sequence.\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* turn on fine structure mixing of spontaneous decays.  
		 * >>refer	Helike	FSM	Bauman, R., MacAdams, K., and Ferland, G. (2003).	*/
		if( p.nMatch(" OFF") )
			iso_ctrl.lgFSM[ipISO] = false;
		else
			iso_ctrl.lgFSM[ipISO] = true;
	}

	else if( p.nMatch("GBAR") )
	{
		if( ipISO == ipH_LIKE )
		{
			fprintf(ioQQQ," Sorry, the GBAR option is only implemented for the He-like sequence.\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* the HEGBAR command - to change cs of higher levels */
		/* first turn all off, one will be turned back on */
		iso_ctrl.lgCS_Vriens[ipISO] = false;
		iso_ctrl.lgCS_None[ipISO] = false;
		iso_ctrl.nCS_new[ipISO] = false;

		/* now turn one on */
		if( p.nMatch("VRIE") )
		{
			/* option to change how collisions for unknown levels affect things */
			iso_ctrl.lgCS_Vriens[ipISO] = true;
		}
		else if( p.nMatch(" NEW") )
		{
			/* option to change how collisions for unknown levels affect things */
			iso_ctrl.nCS_new[ipISO] = (int)p.FFmtRead();

			/* there are two options, 1 and 2, for which fit - default (no number)
			 * will be 1, the broken power law fit */
			if( p.lgEOL() )
				iso_ctrl.nCS_new[ipISO] = 1;

			ASSERT( iso_ctrl.nCS_new[ipISO] );
		}
		else if( p.nMatch(" OFF") )
		{
			/* option to change how collisions for unknown levels affect things */
			iso_ctrl.lgCS_None[ipISO] = true;
		}
		else
		{
			fprintf( ioQQQ, " needs parameter\n" );
			cdEXIT(EXIT_FAILURE);
		}
	}


	else if( p.nMatch("LYMA") )
	{
		if( ipISO == ipH_LIKE && p.nMatch("PUMP") )
		{
			/* >>chng 05 jul 08, separate out Lyman pump commands */
			if( p.nMatch(" OFF") )
			{
				/* option to turn off all continuum pumping of Lyman lines */
				hydro.lgLymanPumping = false;
			}
			else if( p.nMatch("SCALE") )
			{
				/* multiplicative factor for all continuum pumping of H I Lyman lines,
				 * account for possible emission in the line - only affects H I
				 * not entire H-like iso sequence */
				hydro.xLymanPumpingScaleFactor = 
					(realnum)p.FFmtRead();
				/* scale factor is log if <=0, 
				 * represents line in emission if >1
				 * LOG keyword forces interpretation as a log */
				if( hydro.xLymanPumpingScaleFactor <= 0. ||
					p.nMatch(" LOG") )
				{
					hydro.xLymanPumpingScaleFactor = 
						(realnum)pow((realnum)10.f , hydro.xLymanPumpingScaleFactor );
				}

				/* vary option */
				if( optimize.lgVarOn )
				{
					optimize.nvarxt[optimize.nparm] = 1;
					strcpy( optimize.chVarFmt[optimize.nparm], "ATOM H-LIKE LYMAN PUMPING SCALE %f LOG" );

					/*  pointer to where to write */
					optimize.nvfpnt[optimize.nparm] = input.nRead;

					/*  current parameters - always log so steps are log  */
					optimize.vincr[optimize.nparm] = 0.1f;
					optimize.vparm[0][optimize.nparm] = (realnum)log10(hydro.xLymanPumpingScaleFactor);
					++optimize.nparm;
				}
			}
			else
			{
				fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM H-LIKE LYMAN PUMP command.\n");
				fprintf(ioQQQ," The options are \" OFF\", and \"SCALE\".\n");  
				cdEXIT(EXIT_FAILURE);
			}
		}
		else if( p.nMatch("EXTRA") )
		{
			/* option to set number of "extra" Lyman lines, used for optical depths only */
			iso_ctrl.nLyman[ipISO] = (long int)p.FFmtRead();
			if( p.lgEOL() )
				p.NoNumb("'extra' Lyman lines");
			if( iso_ctrl.nLyman[ipISO] < 2 )
			{
				// Code does not elsewhere protect against values less than 2.
				fprintf(ioQQQ," Sorry, the value on this ATOM xx-LIKE LYMAN command must be at least 2.\n"); 
				cdEXIT(EXIT_FAILURE);
			}
		}
		else
		{
			fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM xx-LIKE LYMAN command.\n");
			fprintf(ioQQQ," The options are \"PUMP\", and \"EXTRA\".\n");
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* don't interpolate on look-up tables but compute recombination on the fly instead */
	else if( p.nMatch("RECO") &&
		p.nMatch(" NO ") && p.nMatch("INTE") )
	{
		/* flag set by atom xx-like no recombination interp command, 
		 * says to generate recombination coefficients
		 * on the fly */
		iso_ctrl.lgNoRecombInterp[ipISO] = true;
	}

	else if( p.nMatch("REDI") )
	{
		int ipRedis=0;
		/* there are three functions, PRD_, CRD_, and CRDW,
		 * representing partial redistribution, 
		 * complete redistribution (doppler core only, no wings)
		 * and complete with wings */
		/* partial redistribution */
		if( p.nMatch(" PRD") )
		{
			ipRedis = ipPRD;
		}
		/* complete redistribution no wings */
		else if( p.nMatch(" CRD") )
		{
			ipRedis = ipCRD;
		}
		/* complete redistribution with wings */
		else if( p.nMatch("CRDW") )
		{
			ipRedis = ipCRDW;
		}

		/* if not SHOW option (handled below) then we have a problem */
		else if( !p.nMatch("SHOW") )
		{
			fprintf(ioQQQ," There should have been a second keyword on this command.\n");
			fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space).  Sorry.\n");
			cdEXIT(EXIT_FAILURE);
		}

		/* resonance lines - not Lya*/
		if( p.nMatch("ALPH") )
		{
			iso_ctrl.ipLyaRedist[ipISO] = ipRedis;
		}
		/* Lya itself */
		else if( p.nMatch("RESO") )
		{
			iso_ctrl.ipResoRedist[ipISO] = ipRedis;
		}
		/* subordinate lines */
		else if( p.nMatch("SUBO") )
		{
			iso_ctrl.ipSubRedist[ipISO] = ipRedis;
		}
		/* the show option, say what we are assuming */
		else if( p.nMatch("SHOW") )
		{
			fprintf(ioQQQ," Ly a is ");
			if( iso_ctrl.ipLyaRedist[ipISO] ==ipCRDW )
			{
				fprintf(ioQQQ,"complete redistribution with wings\n");
			}
			else if( iso_ctrl.ipLyaRedist[ipISO] ==ipCRD )
			{
				fprintf(ioQQQ,"complete redistribution with core only.\n");
			}
			else if( iso_ctrl.ipLyaRedist[ipISO] ==ipPRD )
			{
				fprintf(ioQQQ,"partial redistribution.\n");
			}
			else if( iso_ctrl.ipLyaRedist[ipISO] ==ipLY_A )
			{
				fprintf(ioQQQ,"special Lya.\n");
			}
			else
			{
				fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipLyaRedist.\n");
				TotalInsanity();
			}

			fprintf(ioQQQ," Other %s resonance lines are ",
				elementnames.chElementSym[ipISO] );

			if( iso_ctrl.ipResoRedist[ipISO] ==ipCRDW )
			{
				fprintf(ioQQQ,"complete redistribution with wings\n");
			}
			else if( iso_ctrl.ipResoRedist[ipISO] ==ipCRD )
			{
				fprintf(ioQQQ,"complete redistribution with core only.\n");
			}
			else if( iso_ctrl.ipResoRedist[ipISO] ==ipPRD )
			{
				fprintf(ioQQQ,"partial redistribution.\n");
			}
			else
			{
				fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipResoRedist.\n");
				TotalInsanity();
			}

			fprintf(ioQQQ," %s subordinate lines are ",
				elementnames.chElementSym[ipISO] );

			if( iso_ctrl.ipSubRedist[ipISO] ==ipCRDW )
			{
				fprintf(ioQQQ,"complete redistribution with wings\n");
			}
			else if( iso_ctrl.ipSubRedist[ipISO] ==ipCRD )
			{
				fprintf(ioQQQ,"complete redistribution with core only.\n");
			}
			else if( iso_ctrl.ipSubRedist[ipISO] ==ipPRD )
			{
				fprintf(ioQQQ,"partial redistribution.\n");
			}
			else
			{
				fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipSubRedist.\n");
				TotalInsanity();
			}
		}
		else
		{
			fprintf(ioQQQ," here should have been another keyword on this command.\n");
			fprintf(ioQQQ," Options are ALPHA, RESONANCE, SUBORDINATE.  Sorry.\n");
			cdEXIT(EXIT_FAILURE);
		}
	}

	else if( p.nMatch("TOPO") )
	{
		if( p.nMatch(" OFF") )
		{
			iso_ctrl.lgTopoff[ipISO] = false;
			fprintf( ioQQQ, "ISO %li TOPOFF is OFF\n", ipISO );
		}
		else
			iso_ctrl.lgTopoff[ipISO] = true;
	}


	else
	{
		fprintf( ioQQQ, " There should have been a keyword on this ATOM H-LIKE or HE-LIKE command.\n Sorry.\n" );
		cdEXIT(EXIT_FAILURE);
	}
	return;
}