/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*lines_hydro put H-iso seq into line intensity stack */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "secondaries.h"
#include "lines_service.h"
#include "hydrogenic.h"
#include "dense.h"
#include "atmdat.h"
#include "rfield.h"
#include "phycon.h"
#include "geometry.h"
#include "radius.h"
#include "lines.h"

void lines_hydro(void)
{
	long int i, nelem, ipHi, ipLo;
	double hbetab, 
		em , 
		EmisFac, 
		pump ,
		caseb;
	char chLabel[5];

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

	/* next come the H-iso sequence lines */
	i = StuffComment( "H-like iso-sequence" );
	linadd( 0., (float)i , "####", 'i');

	/* total collisional cooling due to all hydrogen lines */
	linadd(MAX2(0.,iso.xLineTotCool[ipH_LIKE][ipHYDROGEN]),912,"Clin",'c');

	/* total collisional heating due to all hydrogen lines */
	linadd(MAX2(0.,-iso.xLineTotCool[ipH_LIKE][ipHYDROGEN]),912,"Hlin",'h');
	/*fprintf(ioQQQ," debugg\t%.2e\t%.2e\t%.2e\n", 
		radius.drad,
		iso.xLineTotCool[ipH_LIKE][ipHYDROGEN] , 
		iso.cLya_cool[ipH_LIKE][ipHYDROGEN]);*/

	/*>>chng 95 jun 25 changed from info to cooling to pick this up in primal.in */
	/* collisionally excited La cooling */
	linadd(MAX2(0.,iso.cLya_cool[ipH_LIKE][ipHYDROGEN]),1216,"Cool",'i');

	/* collisionally de-excited La heating */
	linadd(MAX2(0.,-iso.cLya_cool[ipH_LIKE][ipHYDROGEN]),1216,"Heat",'i');

	/* cooling due to n>2 Lyman lines */
	linadd(MAX2(0.,iso.cLyrest_cool[ipH_LIKE][ipHYDROGEN]),960,"Crst",'i');

	/* heating due to n>2 Lyman lines */
	linadd(MAX2(0.,-iso.cLyrest_cool[ipH_LIKE][ipHYDROGEN]),960,"Hrst",'i');

	/* cooling due to n>3 Balmer lines */
	linadd(MAX2(0.,iso.cBal_cool[ipH_LIKE][ipHYDROGEN]),4861,"Crst",'i');

	/* heating due to n>3 Balmer lines */
	linadd(MAX2(0.,-iso.cBal_cool[ipH_LIKE][ipHYDROGEN]),4861,"Hrst",'i');

	/* cooling due to higher Paschen lines */
	linadd(MAX2(0.,iso.cRest_cool[ipH_LIKE][ipHYDROGEN]),0,"Crst",'i');

	/* heating due to higher Paschen lines */
	linadd(MAX2(0.,-iso.cRest_cool[ipH_LIKE][ipHYDROGEN]),0,"Hrst",'i');

	/* remember largest fractional ionization of H due to secondaries */
	secondaries.SecHIonMax = MAX2( secondaries.SecHIonMax , secondaries.sec2total ) ;

	/* remember fraction of H ionizations due to ct */
	atmdat.HIonFracMax = MAX2( atmdat.HIonFracMax, atmdat.HIonFrac);

	/* remember largest fraction of thermal collisional ionization of H ground state */
	hydro.HCollIonMax = 
		(float)MAX2( hydro.HCollIonMax , hydro.H_ion_frac_collis );

	/* la contribution from suprathermal secondaries from ground */
	linadd(secondaries.x12tot*dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*1.634e-11,1216,"LA X" ,'i');

	/* factor of 0.4836 is ratio of A(4-2)/(A(4-3)+A(4-2))
	 * the IPLNPUMP is the actual pumping rate per atom */
	/* H-beta produced by continuum pumping in optically thin ld limit */
	pump = (double)(EmisLines[ipH_LIKE][ipHYDROGEN][4][ipH1s].pump*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]*
	  dense.xIonDense[ipHYDROGEN][1]*4.09e-12*0.4836);
	linadd(pump,4861,"Pump",'r');

	/* CHION is HEX*HII (in COOLR) and is net col ionz-3 body heat */
	/* collision ionization cooling of hydrogen */
	linadd(MAX2(0.,iso.coll_ion[ipH_LIKE][ipHYDROGEN]),0,"CION",'c');

	/* this is the heating due to 3-body recombination */
	linadd(MAX2(-iso.coll_ion[ipH_LIKE][ipHYDROGEN],0.),0,"3bHt",'h');

	/* Stark broadening contribution to line */
	linadd(dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]*0.*hydro.pestrk[ipH2p][ipH1s]*1.634e-11,
	  1216,"Strk",'i');

	/* Stark broadening contribution to line */
	linadd(dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][3]*hydro.pestrk[3][2]*3.025e-12,
	  6563,"Strk",'i');

	/* Stark broadening contribution to line */
	linadd(dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][4]*hydro.pestrk[4][2]*4.084e-12,
	  4861,"Strk",'i');

	/* Stark broadening contribution to line */
	linadd(dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][4]*hydro.pestrk[4][3]*1.059e-12,
	  18751,"Strk",'i');

	/* pestrk[5,4] is A[4,5]*pest[4,5] 
	 * Stark broadening contribution to line */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 5 )
		linadd(dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][5]*hydro.pestrk[5][4]*4.900e-13,
	  40512,"Strk",'i');

	/* this can fail if RT_line_all never updates the ots rates, a logic error,
	 * but only assert this during actual calculation (ipass>0), */
	ASSERT( LineSave.ipass  <1 ||
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ots>= 0.);
	/* portion of line lost due to absorp by background opacity */
	linadd(EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ots*EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg, 1216,"Dest",'i');

	/* portion of line lost due to absorp by background opacity */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 3 )
		linadd(EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH2s].ots*EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH2s].EnergyErg, 6563,"Dest",'i');

	/* portion of line lost due to absorp by background opacity */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 5 )
		linadd(EmisLines[ipH_LIKE][ipHYDROGEN][5][4].ots*EmisLines[ipH_LIKE][ipHYDROGEN][5][4].EnergyErg,40516, "Dest",'i');

	/* portion of line lost due to absorp by background opacity */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 4 )
		linadd(EmisLines[ipH_LIKE][ipHYDROGEN][4][ipH2s].ots*EmisLines[ipH_LIKE][ipHYDROGEN][4][ipH2s].EnergyErg, 4861,"Dest",'i');

	/* portion of line lost due to absorp by background opacity */
	if( iso.numLevels[ipH_LIKE][ipHYDROGEN] > 4 )
		linadd(EmisLines[ipH_LIKE][ipHYDROGEN][4][3].ots*EmisLines[ipH_LIKE][ipHYDROGEN][4][3].EnergyErg ,18751, "Dest",'i');

	/* Ly-alpha destroyed by overlap with FeII */
	linadd(iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]*dense.xIonDense[ipHYDROGEN][1]*EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul*
	  hydro.dstfe2lya*EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg , 1216 , "Fe 2" , 'i' );

	/* case b intensity of Ly-alpha, no two photon */
	/* do we want intensity/luminosity or emissivity */
	if( hydro.lgHydEmiss )
	{
		/* option to print emissivity, not intensity */
		EmisFac = 1.;
	}
	else
	{
		EmisFac = dense.xIonDense[ipHYDROGEN][1]*dense.eden;
	}

	/* simple high-density case b intensity of Ly-alpha, no two photon */
	linadd(iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*EmisFac*1.64e-11,1216,"Ca B",'i');

	/* H-beta computed from Q(H) and specified covering factor */
	if( nzone == 1 )
	{
		/* evaluate the case b emissivity by interpolating on the hummer & storey tables */
		caseb = rfield.qhtot*
			atmdat_HS_caseB( 4 , 2 , 1 , phycon.te , dense.eden, 'b' ) / iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN];
		/* the atmdat_HS_caseB returned -1 if the physical conditions were outside range of validity.  
		 * In this case use simple approximation with no temperature or density dependence */
		if( caseb < 0 )
		{
			caseb = rfield.qhtot*4.75e-13;
		}
		LineSv[LineSave.nsum].sumlin = 0.;
	}
	else
	{
		caseb = 0.;
	}
	/* >>chng 02 nov 05, better approximation for Lya for temperature of first zone */
	/* H-beta computed from Q(H) and specified covering factor */
	linadd( caseb/radius.dVeff*geometry.covgeo , 4861 , "Q(H)" , 'i' );

	if( nzone == 1 )
	{
		caseb = rfield.qhtot*EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].EnergyErg;
		LineSv[LineSave.nsum].sumlin = 0.;
	}
	else
	{
		caseb = 0.;
	}
	/* >>chng 02 nov 05, better approximation for Lya for temperature of first zone */
	/* Ly-alpha from Q(H), high-dens lim, specified covering factor */
	linadd( caseb/radius.dVeff*geometry.covgeo , 1216 , "Q(H)" , 'i' );

	/* this is main loop creating hydrogenic intensity or emissivity */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( hydro.lgHydEmiss )
		{
			/* option set with Hydrogen emissivity command
			 * will produce emissivity not intensity - luminosity */
			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem]  )
				{
					for( ipHi=ipH2s; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
					{
						for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
						{
							/* this is in real units not emissivit*/
							EmisLines[ipH_LIKE][nelem][ipHi][ipLo].phots = 
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc*
								dense.xIonDense[nelem][nelem+1];

							/* now find line intensity  */
							EmisLines[ipH_LIKE][nelem][ipHi][ipLo].xIntensity = 
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc*
							    (double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg/dense.eden;

						}
					}
				}
			}
		}

		else
		{
			/* this is the default case - luminosity - intensity */
			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				if( dense.lgElmtOn[nelem]  )
				{
					for( ipHi=ipH2s; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
					{
						for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
						{
							EmisLines[ipH_LIKE][nelem][ipHi][ipLo].phots = 
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
								EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc*
								dense.xIonDense[nelem][nelem+1];

							/* now find line intensity and line ots rates */
							/* set hydrogen line emission and ots fields */
							/* >>chng 01 jan 15, add double cast to work at low densities */
							EmisLines[ipH_LIKE][nelem][ipHi][ipLo].xIntensity = 
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc*
								(double)dense.xIonDense[nelem][nelem+1]*
								(double)EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg;

						}
					}
				}
			}
		}
	}

	nelem = 0;
	ipLo = 0;
	ipHi = 1;
#	if 0
	fprintf(ioQQQ,"%.2e %.2e %.2e %.2e %.2e %.2e \n",
		iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2s] , EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi,
		8.23, EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc,
		1.634e-11, EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg );
	fprintf(ioQQQ,"%.2e %.2e \n",
		dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2s]*8.23*1.634e-11, 
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].xIntensity );
#	endif

	/* create emissivity or intensity for hydrogenic species,
	 * first combine/bring balmer series together */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( dense.IonHigh[nelem] == nelem + 1 )
		{
			/* bring n=> 2s + 2p together into one line with 2s lower level
			 * must be 2s since it comes first in line stack */
			for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][nelem]; ipHi++ )
			{
				EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].xIntensity += 
					EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].xIntensity;

				/* kill 2p */
				EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].xIntensity = 0.;
			}
		}
	}

	/* H beta recombination, assuming old case B */
	hbetab = (double)((pow(10.,-20.89 - 0.10612*POW2(phycon.alogte - 4.4)))/
	  phycon.te);
	/* convert to emission per unit then dump into line stack */
	if( !hydro.lgHydEmiss )
	{
		hbetab *= dense.xIonDense[ipHYDROGEN][1]*dense.eden;
	}
	/* this is old case b, had been in code for very long time */
	linadd(hbetab,4861,"CaBo",'i');

	/* 1640 1640 1640 */
	em = 2.03e-20/(phycon.te70*phycon.te10*phycon.te03);
	if( !hydro.lgHydEmiss )
	{
		em *= dense.xIonDense[ipHELIUM][2]*dense.eden;
	}

	/* old prediction of He II 1640, case B at low densities */
	linadd(em,1640,"CaBo",'i');

	/* hydrogenic helium */
	/* old prediction of He II 4686, case B */
	em = 2.52e-20/(pow(phycon.te,1.05881f));
	if( !hydro.lgHydEmiss )
	{
		em *= dense.xIonDense[ipHELIUM][2]*dense.eden;
	}
	linadd(em,4686,"CaBo",'i');

	/* predict case b intensities of hydrogen lines */

	if( LineSave.ipass <= 0 )
	{
		for(nelem=0; nelem<HS_NZ; ++nelem )
		{
			atmdat.lgHCaseBOK[0][nelem] = TRUE;
			atmdat.lgHCaseBOK[1][nelem] = TRUE;
		}
	}
	/* this is the main printout, where line intensities are entered into the stack */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			/* HS_NZ is limit to HS predictions, now 8 */
			/* but don't do the minor elements - these were not read in and so should not be
			 * printed - remove equivalent if statement in createdata to read them in */
			if( nelem < HS_NZ && (nelem<2 || nelem>4) )
			{
				int iCase;
				for( iCase=0; iCase<2; ++iCase )
				{
					char chAB[2]={'A','B'};
					char chLab[5]="Ca  ";
					/*for( ipLo=ipH2p; ipLo<MIN2(6,iso.numLevels[ipH_LIKE][nelem]); ++ipLo )*/
					/* adding iCase means start from n=1 for case A, n=2 for Case B,
					 * note that prin quant number is on physics scale, not C */
					for( ipLo=1+iCase; ipLo<MIN2(6,iso.numLevels[ipH_LIKE][nelem]); ++ipLo )
					{
						for( ipHi=ipLo+1; ipHi< MIN2(ipLo+5,iso.numLevels[ipH_LIKE][nelem]); ++ipHi )
						{
							float wl;
							/* Put case b predictions into line stack
							 * NB NB NB each Hummer & Storey case b line must be 
							 * explicitly clobbered by hand in routine final if 
							 * atmdat.lgHCaseBOK[iCase][nelem] flag is set false
							 * since this indicates that we exceeded bounds of table,
							 * DO NOT want to print lines in that case */

							/* first do case b emssivity of balmer lines */

							/* get HS predictions */
							hbetab = atmdat_HS_caseB( ipHi,ipLo , nelem+1, phycon.te , dense.eden, chAB[iCase] );
							if( hbetab<=0. )
							{
								atmdat.lgHCaseBOK[iCase][nelem] = FALSE;
								hbetab = 0.;
							}

							/* convert to emission per unit vol then dump into line stack */
							if( !hydro.lgHydEmiss )
							{
								hbetab *= dense.xIonDense[nelem][nelem+1]*dense.eden;
							}

							/* this is wavelength of interpolated case b from HS tables */
							if( iCase==0 && ipLo==1 )
							{
								wl = EmisLines[ipH_LIKE][nelem][ipHi][0].WLAng;
							}
							else
							{
								wl = EmisLines[ipH_LIKE][nelem][ipHi][ipLo].WLAng;
							}

							/* make label either Ca A or Ca B */
							chLab[3] = chAB[iCase];
							/* case a or case b from Hummer & Storey tables */
							linadd(hbetab,wl,chLab,'i' );
						}
					}
				}
			}
			/* do not print more than first 50 levels */
			/* NB NB - low and high must be in this order so that all balmer, paschen,
			 * etc series line up correctly in final printout */

			/* >>chng 01 sep 10, do two-photon here */
			if( LineSave.ipass == 0 )
			{
				/* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" 
				 * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
				/* two photon emission */
				chIonLbl(chLabel, &EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s]);
			}
			linadd( EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].xIntensity , 0,chLabel,'r');

			/* induced two photon emission */
			linadd(
				iso.Pop2Ion[ipH_LIKE][nelem][ipH2s]*
				dense.xIonDense[nelem][nelem+1]*
				iso.TwoNu_induc_dn[ipH_LIKE][nelem]*
				EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].EnergyErg,
				22, chLabel ,'i');

			for( ipLo=ipH1s; ipLo < MIN2(50,iso.numLevels[ipH_LIKE][nelem]-1); ipLo++ )
			{
				/* don't bother with decays to 2p ince we set them to zero above */
				if( ipLo==ipH2p )
					continue;

				/* >>chng 01 feb 23, do not print very highest level since edge
				 * effects distort the intensity */
				for( ipHi=ipLo+1; ipHi < MIN2(50,iso.numLevels[ipH_LIKE][nelem]-1); ipHi++ )
				{
					/* skip the 2s-2p and two photon transitions */
					/* >>chng 02 feb 10, had been explicit levels, change to gneral ipCont */
					if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
						continue;
					/* this will enter .xIntensity into the line stack */
					PutLine(&EmisLines[ipH_LIKE][nelem][ipHi][ipLo]);
				}
			}
		}
	}

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

