/* 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 */
/*lines main routine to put emission line intensities into line stack,
 * calls lineset1, 2, 3, 4 */
/*Drive_cdLine do the drive cdLine command */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "thermal.h"
#include "yield.h"
#include "ionbal.h"
#include "cddrive.h"
#include "trace.h"
#include "dense.h"
#include "prt.h"
#include "coolheavy.h"
#include "rfield.h"
#include "phycon.h"
#include "elementnames.h"
#include "iso.h"
#include "hyperfine.h"
#include "hydrogenic.h"
#include "lines_service.h"
#include "atmdat.h"
#include "lines.h"
#include "radius.h"
static void Drive_cdLine( void );

void lines(void)
{
	char chLabel[5];
	static int lgRecOn;
	long int i, 
	  ipnt,
	  nelem;
	double BigstExtra, 
	  ExtraCool,  
	  f2, sum; 

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

	/* LineSave.ipass
	 * -1 - space not yet allocated - just count number of lines entered into stack
	 *  0 - first call with space allocated - must create labels and add in wavelengths
	 * +1 - later calls - just add intensity 
	 */

	/* major routines used here:
	 *
	 * PutLine( tarray )
	 * this uses information in tarray to give various
	 * contributions to lines, and their intensities
	 *
	 * PutExtra( extra )
	 * extra is some extra intensity to add to the line
	 * it will go into the totl contribution put out by PutLine,
	 * and this contriibution should be indicated by independent
	 * call to linadd
	 * */

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " lines called\n" );
	}

	/* the drive cdline command, which checks that all lines can be pulled out by cdLine */
	if( trace.lgDrv_cdLine  && LineSave.ipass > 0 )
		Drive_cdLine();

	/* total luminosity radiated by this model - will be compared with energy in incident
	 * continuum when calculation is complete */
	thermal.power += thermal.htot*radius.dVeff;

	/* remember the total free-free heating */
	thermal.FreeFreeTotHeat += CoolHeavy.brems_heat_total*radius.dVeff;

	/* total compton heating - cooling */
	rfield.comtot += rfield.cmheat*radius.dVeff;
	thermal.totcol += thermal.ctot*radius.dVeff;

	/* up up induced recom cooling */
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		hydro.cintot += iso.RecomInducCool_Rate[ipH_LIKE][nelem]*radius.dVeff;
	}

	/* nsum is line pointer within large stack of line intensities */
	LineSave.nsum = 0;
	LineSave.ndsum = 0;
	LineSave.nComment = 0;

	/* last arg in call to lindst and linadd is info on line
	 * info is char variable indicating type of line this is
	 * 'c' cooling
	 * 'h' heating
	 * 'i' information only
	 * 'r' recombination line
	 *
	 * all components of lines are entered into the main line stack here
	 * when printing, filters exist to not print Inwd component */
	 
	/* initialize routine that can generate pointers for forbidden lines,
	 * these are lines that are not transferred otherwise,
	 * in following routines there will be pairs of calls, first to
	 * PntForLine to get pointer, then lindst to add to stack */
	PntForLine(0.,"FILL",&i);

	/* evaluate rec coef for rec lines of C, N, O
	 * some will be used in LineSet2 and then zeroed out,
	 * others left alone and used below */
	atmdat_rec_lines(phycon.te,LineSave.RecCoefCNO);

	/* initialize ExtraInten with a zero */
	PutExtra(0.);

	/* put in something impossible in element 0 of line stack */
	linadd(0.f,0,"zero",'i');

	/* this is compared with true volume in final.  The number can't
	 * actualy be unity since this would overflow on a 32 bit machine */
	/* integrate the volume as a sanity check */
	linadd( 1.e-10 , 1 , "Unit" , 'i' );

	/* initial set of general properties */
	lines_general();

	/* do all continua */
	lines_continuum();

	/* now do all molecules */
	lines_molecules();

	/* information about grains */
	lines_grains();

	/* do all hydrogenic ions */
	lines_hydro();

	/* enter He-iso sequence lines */
	lines_helium();

#if	0
	/* This is Ryan's code for dumping lots of Helium lines according to
	 * quantum number rather than wavelength, principally for comparison with Rob
	 * Bauman's code. */
	{
		long ipHi, ipLo;
		
		fprintf( ioQQQ,"ipHi\tipLo\tnu\tlu\tsu\tnl\tll\tsl\tWL\tintens\n" );
		for( ipHi=5; ipHi<= 419; ipHi++ )
		{
			for( ipLo=0; ipLo<ipHi; ipLo++ )
			{
				if( EmisLines[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng < 1.0E6 &&
					EmisLines[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng > 2.0E3 &&
					ipLo!=3 && ipLo!=4 )
				{
					fprintf( ioQQQ,"%li\t%li\t%li\t%li\t%li\t%li\t%li\t%li\t%e\t",
						ipHi,
						ipLo,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].n,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].l,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipHi].s,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].n,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].l,
						iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].s,
						EmisLines[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng );

					if( ipHi==5 )
					{
						fprintf( ioQQQ,"%e\t", 
							EmisLines[ipHE_LIKE][ipHELIUM][3][ipLo].xIntensity+
							EmisLines[ipHE_LIKE][ipHELIUM][4][ipLo].xIntensity+
							EmisLines[ipHE_LIKE][ipHELIUM][5][ipLo].xIntensity );
					}
					else if( ipLo==5 )
					{
						fprintf( ioQQQ,"%e\t", 
							EmisLines[ipHE_LIKE][ipHELIUM][ipHi][3].xIntensity+
							EmisLines[ipHE_LIKE][ipHELIUM][ipHi][4].xIntensity+
							EmisLines[ipHE_LIKE][ipHELIUM][ipHi][5].xIntensity );
					}
					else
					{
						fprintf( ioQQQ,"%e\t", 
							EmisLines[ipHE_LIKE][ipHELIUM][ipHi][ipLo].xIntensity );
					}

					fprintf( ioQQQ, "\n");
				}
			}
		}
	}
#endif

	/* do heavies, lithium through neon */
	lines_lv1_li_ne();

	/* do heavies, sodium through argon */
	lines_lv1_na_ar();

	/* do heavies, potasium through zinc */
	lines_lv1_k_zn();

	/* add up line intensities for certain set of lines */
	sum = PrtLineSum( " SUM" );
	/* zero out the location that will receive this information, 
	 * remember that memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;
	/* optional sum of certain emission lines, set with "print sum" */
	linadd(sum/radius.dVeff,0,"Stoy",'i' );

	/* next come some recombination lines */
	i = StuffComment( "recombination" );
	linadd( 0., (float)i , "####", 'i' );

	/***********************************************************************
	 * large number of C, N, and O recombination lines                     *
	 *************************************************************************/

	if( LineSave.ipass <= 0 )
	{
		/* initialize to TRUE, may be set false if density is very high below */
		lgRecOn = TRUE;
	}

	for( i=0; i < 471; i++ )
	{
		/* generate label for the line */
		strcpy( chLabel, elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][i])-1] );
		strcat( chLabel, elementnames.chIonStage[(long)(LineSave.RecCoefCNO[0][i]-
			LineSave.RecCoefCNO[1][i]+1.01)-1] );

		/* number of rec per unit vol
		 * >>chng 97 aug 28, do not predict rec intensities at high densities
		 * blr.in and oldblr.in go nuts if this is added in outward beam */
		if( dense.eden < 1e8 && lgRecOn )
		{
			f2 = LineSave.RecCoefCNO[3][i]*dense.eden*
				dense.xIonDense[(long)(LineSave.RecCoefCNO[0][i])-1][(long)(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+2)-1];

			/* convert to intensity */
			f2 = f2*1.99e-8/LineSave.RecCoefCNO[2][i];
		}
		else
		{
			lgRecOn = FALSE;
			f2 = 0.;
		}
		/* finally stuff it into the stack */
		PntForLine(LineSave.RecCoefCNO[2][i],chLabel,&ipnt);
		lindst(f2,LineSave.RecCoefCNO[2][i],chLabel,ipnt, 'r',TRUE );
	}

	/* next come the atom_level2 lines */
	i = StuffComment( "level2 lines" );
	linadd( 0., (float)i , "####", 'i' );

	/* add in all the other level 2 wind lines
	 * Dima'a 6k lines */
	ExtraCool = 0.;
	BigstExtra = 0.;
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			PutLine(&TauLine2[i]);
			if( TauLine2[i].cool > BigstExtra )
			{
				BigstExtra = TauLine2[i].cool;
				thermal.ipMaxExtra = i+1;
			}
			ExtraCool += TauLine2[i].cool;
		}
	}
	/* keep track of how important this is */
	thermal.GBarMax = MAX2(thermal.GBarMax,(float)(ExtraCool/thermal.ctot));

	/* next come the hyperfine structure lines */
	i = StuffComment( "hyperfine structure" );
	linadd( 0., (float)i , "####", 'i' );

	/* this is total cooling due to all HF lines */
	linadd( hyperfine.cooling_total, 0., "hfin", 'i' );

	/* remember largest local cooling for possible printout in comments */
	hyperfine.cooling_max = (float)MAX2(hyperfine.cooling_max,hyperfine.cooling_total/thermal.ctot);

	/* the hyperfine lines */
	for( i=0; i < nHFLines; i++ )
	{
		PutLine(&HFLines[i]);
	}

	/* next come the inner shell fluorescent lines */
	i = StuffComment( "inner shell" );
	linadd( 0., (float)i , "####", 'i' );

	/* the group of inner shell fluorescent lines */
	for( i=0; i<yield.nfl_lines; ++i )
	{
		double xInten = 
			/* density of parent ion, cm-3 */
			dense.xIonDense[yield.nfl_nelem[i]][yield.nfl_ion[i]] *
			/* photo rate per atom per second, s-1 */
			ionbal.PhotoRate_Shell[yield.nfl_nelem[i]][yield.nfl_ion[i]][yield.nfl_nshell[i]][0]*
			/* fluor yield - dimensionless */
			yield.fl_yield[i] *
			/* photon energy - ryd, converted into ergs */
			yield.fl_energy[i] * EN1RYD;

		/* create label if initializing line stack */
		if( LineSave.ipass == 0 )
		{
			/* only generate the line label if it is going to be used */
			strcpy( chLabel , elementnames.chElementSym[yield.nfl_nelem[i]] );
			strcat( chLabel , elementnames.chIonStage[yield.nfl_ion_emit[i]] );
#			if 0
			/* only print yields for atoms */
			if( yield.nfl_ion[i] == 0 && yield.nfl_nelem[i]==ipIRON )
			fprintf(ioQQQ,"DEBUGyeild\t%s\t%.3f\t%.3e\n",
				/* line designation, energy in eV, yield */
				chLabel , yield.fl_energy[i]*EVRYD, yield.fl_yield[i] );
#			endif
		}

		/* the group of inner shell fluorescent lines */
		lindst(
			/* intensity of line */
			xInten,
			/* wavelength of line in Angstroms */
			(float)RYDLAM / yield.fl_energy[i] ,
			/* label */
			chLabel ,
			/* continuum array offset for line as set in ipoint */
			yield.nfl_ipoint[i], 
			/* type of line - count as a recombination line */
			'r',
			/* include line in continuum? */
			TRUE );
	}

	/* >>chng 06 jan 03, confirm that number of lines never changes once we have
	 * created the labels */
	{
		static long nLineSave=-1 , ndLineSave=-1;
		if( LineSave.ipass == 0 )
		{
			nLineSave = LineSave.nsum;
			ndLineSave = LineSave.ndsum;
		}
		else if( LineSave.ipass > 0 )
		{
			/* this can't happen */
			if( nLineSave<= 0 || ndLineSave < 0 )
				TotalInsanity();

			/* now make sure that we have the same number of lines as we had previously
			 * created labels.  This would not pass if we did not add in exactly the same
			 * number of lines on each pass */
			if( nLineSave != LineSave.nsum )
			{
				fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
				fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and nLineSave is %li\n",
					LineSave.nsum , 
					nLineSave);
				ShowMe();
				puts( "[Stop in lines" );
				cdEXIT(EXIT_FAILURE);
			}
			if( ndLineSave != LineSave.ndsum )
			{
				fprintf( ioQQQ, "DISASTER number of lines in LineSave.ndsum changed between pass 0 and 1 - this is impossible\n" );
				fprintf( ioQQQ, "DISASTER LineSave.ndsum is %li and ndLineSave is %li\n",
					LineSave.ndsum , 
					ndLineSave);
				ShowMe();
				puts( "[Stop in lines" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " lines returns\n" );
	}

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

/*Drive_cdLine do the drive cdLine command */
static void Drive_cdLine( void )
{
	long int j;
	int lgMiss = FALSE;
	double absval , rel ;

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

	for( j=1; j < LineSave.nsum; j++ )
	{
		if( cdLine( LineSv[j].chALab , LineSv[j].wavelength , &absval , &rel ) <= 0 )
		{
			/* print header if first miss */
			if( !lgMiss )
				fprintf(ioQQQ,"n\tlab\twl\n");
			++lgMiss;

			fprintf(ioQQQ,"%li\t%s\t%f\n", j, LineSv[j].chALab , LineSv[j].wavelength );
		}
	}
	fprintf( ioQQQ, " Thanks for checking on the cdLine routine!\n" );
	puts( "[Stop in Drive_cdLine]" );
	cdEXIT(EXIT_FAILURE);
}
