/* 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_general put general information and energetics into line intensity stack */
/*GetMaxhLine find the strongest heating line */
#include "cddefines.h"
#include "taulines.h"
#include "coolheavy.h"
#include "hydrogenic.h"
#include "dense.h"
#include "thermal.h"
#include "continuum.h"
#include "geometry.h"
#include "dynamics.h"
#include "iso.h"
#include "rfield.h"
#include "trace.h"
#include "ionbal.h"
#include "lines_service.h"
#include "radius.h"
#include "lines.h"
/*GetMaxhLine find the strongest heating line */
static void GetMaxhLine(void);

void lines_general(void)
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  nelem, 
	  ipnt;

	double 
	  EmisFac ,/* used for hydrogenic emissivity or intensity */
	  hbetac, 
	  HeatMetal ,
	  ee511, 
	  hlalph;

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

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

	/* start of general properties */
	i = StuffComment( "general properties" );
	linadd( 0., (float)i , "####", 'i');

	/* do we want intensity/luminosity or emissivity */
	if( hydro.lgHydEmiss )
	{
		/* option to print emissivity, not intensity */
		EmisFac = 1./dense.eden;
	}
	else
	{
		EmisFac = dense.xIonDense[ipHYDROGEN][1];
	}
	/* total H-beta from multi-level atom */
	nelem = 0;
	ipHi = 4;
	ipLo = ipH2p;
	hbetac = 
		(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul * 
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Pesc +
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul *
		EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Pesc ) *
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg*
		EmisFac ;
	/*fprintf(ioQQQ,"nzonehbeta\t%li\t%.3e\n", nzone, hbetac);*/
	/* total H-beta from multi-level atom */
	/*hbetac = EmisLines[ipH_LIKE][ipHYDROGEN][4][ipH2s].xIntensity + */
		/*EmisLines[ipH_LIKE][ipHYDROGEN][4][ipH2p].xIntensity;*/

	/* these lines added to outlin in metdif - following must be false 
	 * this passes array index for line energy in continuum mesh - in rest
	 * of code this is set by a previous call to PntForLine, this index
	 * is on the f not c scale */
	/*linadd(hbetac,4861,"TOTL",'i' );*/
	/* total H-beta from multi-level atom */
	lindst(hbetac,4861,"TOTL",EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].ipCont,'i',FALSE );

	/* total Ly-a from multi-level atom */
	ipHi = ipH2p;
	ipLo = ipH1s;
	hlalph = 
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul* 
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc*
		EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg*
		EmisFac ;
	/*hlalph = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].xIntensity;*/
	linadd(hlalph,1216,"TOTL",'i');

	continuum.totlsv = continuum.totlsv/radius.dVeff*geometry.covgeo;

	/* H 21 cm, A from 
	 * >>refer	H1	A	Gould, ApJ 423, 522
	 * =2.88426e-15 s-1 */
	
	/* total luminosity in incident continuum */
	linadd(continuum.totlsv,0,"Inci",'i');
	continuum.totlsv = continuum.totlsv*radius.dVeff/geometry.covgeo;

	/* ipass is flag to indicate whether to only set up line array
	 * (ipass=0) or actually evaluate lines intensities (ipass=1) */
	if( LineSave.ipass > 0 )
	{
		continuum.totlsv = 0.;
	}

	/* total heating, all forms, information since individuals added later */
	linadd(thermal.htot,0,"TotH",'i');
	/* total cooling, all forms, information since individuals added later */
	linadd(thermal.ctot,0,"TotC",'i');

	/* hydrogen photoionization heating, ground state only */
	linadd(thermal.heating[0][0],0,"BFH1",'h');

	/* net hydrogen photoionization heating less rec cooling, all excited states
	 * normally zero, positive if excited states are net heating */
	linadd(thermal.heating[0][1],0,"BFHx",'h');

	/* heating due to induced lines absorption of continuum */
	linadd(thermal.heating[0][22],0,"Line",'h');
	if( thermal.htot > 0. )
	{
		if( thermal.heating[0][22]/thermal.htot > thermal.HeatLineMax )
		{
			thermal.HeatLineMax = (float)(thermal.heating[0][22]/thermal.htot);
			/* finds the strongest heating line */
			GetMaxhLine();
		}
	}

	/* total helium photoionization heating, all stages */
	linadd(thermal.heating[1][0]+thermal.heating[1][1]+thermal.heating[1][2],
	  0,"BFHe",'h');

	HeatMetal = 0.;
	/* some sums that will be printed in the stack */
	for( nelem=2; nelem<LIMELM; ++nelem)
	{
		/* we now have final solution for this element */
		for( i=dense.IonLow[nelem]; i < dense.IonHigh[nelem]; i++ )
		{
			/* total metal photo heating for LINES */
			HeatMetal += thermal.heating[nelem][i];
		}
	}
	/* total heavy element photoionization heating, all stages */
	linadd(HeatMetal,0,"TotM",'h');

	/* heating due to pair production */
	linadd(thermal.heating[0][21],0,"pair",'h');

	/* ipass is flag to indicate whether to only set up line array
	 * (ipass=0) or actually evaluate lines intensities (ipass=1) */
	if( LineSave.ipass > 0 )
	{
		/* this will be max local heating due to bound compton */
		ionbal.CompHeating_Max = MAX2( ionbal.CompHeating_Max , ionbal.CompRecoilHeatLocal/thermal.htot);
	}
	else
	{
		ionbal.CompHeating_Max = 0.;
	}

	/* heating due to bound compton scattering */
	linadd(ionbal.CompRecoilHeatLocal,0,"Cbnd",'h');

	/* Compton heating */
	linadd(rfield.cmheat,0,"ComH",'h');

	/* total Compton cooling */
	linadd(CoolHeavy.tccool,0,"ComC",'c');

	/* record max local heating due to advection */
	dynamics.HeatMax = MAX2( dynamics.HeatMax , dynamics.Heat /thermal.htot );
	/* record max local cooling due to advection */
	dynamics.CoolMax = MAX2( dynamics.CoolMax , dynamics.Cool /thermal.htot );

	/* cooling due to advection */
	linadd(dynamics.Cool  , 0 , "advC" , 'i' );
	/* heating due to advection */
	linadd(dynamics.Heat , 0 , "advH" , 'i' );

	/* heating due to charge transfer */
	linadd(MAX2(0.,thermal.heating[0][24]),0,"CT H",'h');
	/* cooling due to charge transfer */
	linadd(MAX2(0.,-thermal.heating[0][24]),0,"CT C",'c');

	/* cosmic ray heating */
	linadd(thermal.heating[1][6],0,"CR H",'h');

	/* "extra" heat added to this zone, from HEXTRA command */
	linadd(thermal.heating[0][20],0,"extH",'h');

	/* "extra" cooling added to this zone, from CEXTRA command */
	linadd(CoolHeavy.cextxx,0,"extC",'c');

	/* 511keV annihilation line */
	ee511 = (dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM])*ionbal.PairProducPhotoRate[0]*2.*8.20e-7;
	PntForLine(2.427e-2,"e-e+",&ipnt);
	lindst(ee511,0,"e-e+",ipnt,'i',TRUE );

	/* expansion cooling, only non-zero for wind */
	linadd(CoolHeavy.expans,0,"Expn",'c');

	/* H radiative recombination cooling */
	linadd(iso.RadRecCool[ipH_LIKE][ipHYDROGEN],0,"H FB",'i');

	/* net free-bound cooling */
	linadd(MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBc",'c');
	/* net free-bound heating */
	linadd(MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBh",'h');

	/* cooling due to induced rec of hydrogen */
	linadd(iso.RecomInducCool_Rate[ipH_LIKE][ipHYDROGEN],0,"Hind",'c');

	/* cooling due to induced rec of fully ionized helium */
	linadd(CoolHeavy.c3ind,0,"3He2",'c');

	/* cyclotron cooling */
	linadd(CoolHeavy.cyntrn,0,"Cycn",'c');

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

/*GetMaxhLine find the strongest heating line */
static void GetMaxhLine(void)
{
	long int i;
	double strong;

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

	/* routine called to find which is the strongest heating line */
	strong = 0.;

	/* possible for levlmax to remain 0 if induced processes turned off */
	thermal.levlmax = 0;

	for( i=1; i <= nLevel1; i++ )
	{
		if( TauLines[i].heat > strong )
		{
			strong = TauLines[i].heat;
			thermal.levlmax = 1;
			thermal.ipHeatlmax = i;
		}
	}
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			if( TauLine2[i].heat > strong )
			{
				strong = TauLine2[i].heat;
				thermal.levlmax = 2;
				thermal.ipHeatlmax = i;
			}
		}
	}
	for( i=0; i < nHFLines; i++ )
	{
		if( HFLines[i].heat > strong )
		{
			strong = HFLines[i].heat;
			thermal.levlmax = 2;
			thermal.ipHeatlmax = i;
		}
	}
	for( i=0; i < nCORotate; i++ )
	{
		if( C12O16Rotate[i].heat > strong )
		{
			strong = C12O16Rotate[i].heat;
			thermal.levlmax = 3;
			thermal.ipHeatlmax = i;
		}
		if( C13O16Rotate[i].heat > strong )
		{
			strong = C13O16Rotate[i].heat;
			thermal.levlmax = 4;
			thermal.ipHeatlmax = i;
		}
	}

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

