/* 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 */
/*PrtLinePres print line radiation pressures for current conditions */
#include "cddefines.h"
#include "pressure.h"
#include "taulines.h"
#include "iso.h"
#include "dense.h"
#include "lines_service.h"
#include "h2.h"
#include "rt.h"
#include "prt.h"
/* faintest pressure we will bother with */
#define	THRESH	0.05

void PrtLinePres(void)
{
	long int i, 
	  ip,
	  ipLo,
	  ipHi,
	  nelem;
	int ier;
	double RadPres1;

/*	this is limit to number of lines we can store at one time */
#	define	NLINE	100
	/* labels, wavelengths, and fraction of total pressure, for important lines */
	char chLab[NLINE][5];
	float wl[NLINE] , frac[NLINE] ;
	long int iperm[NLINE];

	/* will be used to check on size of opacity, was capped at this value */
	float smallfloat=SMALLFLOAT*10.f;

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

	/* this routine only called if printout of contributors to line
	 * radiation pressure is desired */

	/* compute radiation pressure due to special H and He lines */
	ip = 0;

	for(i=0; i<NLINE; ++i)
	{
		frac[i] = FLT_MAX;
		wl[i] = FLT_MAX;
	}

	/* this wil be H L alpha
	strcpy( chLab[ip] , "Lya " );
	frac[ip] = pressure.HLalpha/MAX2(smallfloat,(float)pressure.PresRadCurr);
	wl[ip] = 1216.;
	++ip; */

	if( pressure.PresRadCurr > 1e-30 )
	{
		long int ipISO = 0;
		/*TODO	1	make this and eval rad pressure same routine, with flag saying to 
		 * print contributors - copy code from other routine - this code has been
		 * left behind */
		/* line radiation pressure for hydrogen lines */
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( dense.IonHigh[nelem] == nelem + 1 )
			{
				/* >>chng 01 jan 08, do not include highest level since maser can occur,
				 * and pops are set to small number in this case */
				for( ipLo=ipH1s; ipLo < (iso.numLevels[ipH_LIKE][nelem] - 2); ipLo++ ) 
				{
					/* this loop does not include Lya since done just above, 
					 * and set into plints[0] */
					for( ipHi=ipLo+1; ipHi < iso.numLevels[ipH_LIKE][nelem]-1; ipHi++ )
					{
						if( EmisLines[ipISO][nelem][ipHi][ipLo].Aul > 10. &&
							EmisLines[ipISO][nelem][ipHi][ipLo].PopHi > smallfloat &&
							EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc > smallfloat &&
							/* >>chng 01 nov 01, add test that have not overrun optical depth scale */
							( (EmisLines[ipISO][nelem][ipHi][ipLo].TauTot - EmisLines[ipISO][nelem][ipHi][ipLo].TauIn) > smallfloat ) &&
							EmisLines[ipISO][nelem][ipHi][ipLo].ipCont>0)
						{
							RadPres1 = 5.551e-2*
							  (powi(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN/1.e6,4))*
							  (EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi/
							  iso.stat[ipH_LIKE][nelem][ipHi])/
							  (EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc/
							  iso.stat[ipH_LIKE][nelem][ipLo])*
							  /* RT_LineWidth gets line width in terms of RT effects */
							  RT_LineWidth(&EmisLines[ipH_LIKE][nelem][ipHi][ipLo]);
							if( RadPres1 > pressure.PresRadCurr*THRESH )
							{
								wl[ip] = EmisLines[ipH_LIKE][nelem][ipHi][ipLo].WLAng;

								/* put null terminated line label into chLab using line array*/
								chIonLbl(chLab[ip], &EmisLines[ipH_LIKE][nelem][ipHi][ipLo]);
								frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

								ip = MIN2((long)NLINE-1,ip+1);
							}
						}
					}
				}
			}
		}

		/* line radiation pressure for he-like lines */
		ipISO = 1;
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( dense.IonHigh[nelem] >= nelem )
			{
				for( ipHi=1; ipHi <iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
				{
					for( ipLo=0; ipLo < ipHi; ipLo++ ) 
					{
						/*NB - this code must be kep parrallel with code in prt_linepres */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].Aul > 10. &&
							EmisLines[ipISO][nelem][ipHi][ipLo].PopHi > smallfloat &&
							EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc > smallfloat &&
							/* >>chng 01 nov 01, add test that have not overrun optical depth scale */
							( (EmisLines[ipISO][nelem][ipHi][ipLo].TauTot - EmisLines[ipISO][nelem][ipHi][ipLo].TauIn) > smallfloat ) &&
							EmisLines[ipISO][nelem][ipHi][ipLo].ipCont>0)
						{
							RadPres1 = 5.551e-2*
							  (powi(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/1.e6,4))*
							  (EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi/
							  iso.stat[ipHE_LIKE][nelem][ipHi])/
							  (EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc/
							  iso.stat[ipHE_LIKE][nelem][ipLo])*
							  /* RT_LineWidth gets line width in terms of RT effects */
							  RT_LineWidth(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
							if( RadPres1 > pressure.PresRadCurr*THRESH )
							{
								wl[ip] = EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng;

								/* put null terminated line label into chLab using line array*/
								chIonLbl(chLab[ip], &EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
								frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

								ip = MIN2((long)NLINE-1,ip+1);
							}
						}
					}
				}
			}
		}

		/* line radiation pressure from large set of level 1 lines */
		for( i=1; i <= nLevel1; i++ )
		{
			if( TauLines[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(TauLines[i].EnergyWN/
				  1.e6,4))*(TauLines[i].PopHi/TauLines[i].gHi)/
				  (TauLines[i].PopLo/TauLines[i].gLo)*
				  RT_LineWidth(&TauLines[i]);
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.PresRadCurr*THRESH )
			{
				wl[ip] = TauLines[i].WLAng;

				/* put null terminated line label into chLab using line array*/
				chIonLbl(chLab[ip], &TauLines[i]);
				frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

				ip = MIN2((long)NLINE-1,ip+1);
			}
		}

		for( i=0; i < nWindLine; i++ )
		{
			if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
			{
				if( TauLine2[i].PopHi > 1e-30 )
				{
					RadPres1 = 5.551e-2*(powi(TauLine2[i].EnergyWN/
					1.e6,4))*(TauLine2[i].PopHi/TauLine2[i].gHi)/
					(TauLine2[i].PopLo/TauLine2[i].gLo)*
					RT_LineWidth(&TauLine2[i]);
				}
				else
				{
					RadPres1 = 0.;
				}

				if( RadPres1 > pressure.PresRadCurr*THRESH )
				{
					wl[ip] = TauLine2[i].WLAng;

					/* put null terminated line label into chLab using line array*/
					chIonLbl(chLab[ip], &TauLine2[i]);
					frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

					ip = MIN2((long)NLINE-1,ip+1);
				}
			}
		}

		for( i=0; i < nHFLines; i++ )
		{
			if( HFLines[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(HFLines[i].EnergyWN/
				  1.e6,4))*(HFLines[i].PopHi/HFLines[i].gHi)/
				  (HFLines[i].PopLo/HFLines[i].gLo)*
				  RT_LineWidth(&HFLines[i]);
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.PresRadCurr*THRESH )
			{
				wl[ip] = HFLines[i].WLAng;

				/* put null terminated line label into chLab using line array*/
				chIonLbl(chLab[ip], &HFLines[i]);
				frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

				ip = MIN2((long)NLINE-1,ip+1);
			}
		}

		/* radiation pressure due to H2 */
		RadPres1 = H2_RadPress();
		if( RadPres1 > pressure.PresRadCurr*THRESH )
		{
			wl[ip] = 0;

			/* put null terminated 4 char line label into chLab using line array*/
			strcpy(chLab[ip], "H2  ");
			frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

			ip = MIN2((long)NLINE-1,ip+1);
		}

		for( i=0; i < nCORotate; i++ )
		{
			if( C12O16Rotate[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(C12O16Rotate[i].EnergyWN/
				  1.e6,4))*(C12O16Rotate[i].PopHi/C12O16Rotate[i].gHi)/
				  (C12O16Rotate[i].PopLo/C12O16Rotate[i].gLo)*
				  RT_LineWidth(&C12O16Rotate[i]);
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.PresRadCurr*THRESH )
			{
				wl[ip] = C12O16Rotate[i].WLAng;

				/* put null terminated line label into chLab using line array*/
				chIonLbl(chLab[ip], &C12O16Rotate[i]);
				frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

				ip = MIN2((long)NLINE-1,ip+1);
			}
		}

		for( i=0; i < nCORotate; i++ )
		{
			if( C13O16Rotate[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(C13O16Rotate[i].EnergyWN/
				  1.e6,4))*(C13O16Rotate[i].PopHi/C13O16Rotate[i].gHi)/
				  (C13O16Rotate[i].PopLo/C13O16Rotate[i].gLo)*
				  RT_LineWidth(&C13O16Rotate[i]);
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.PresRadCurr*THRESH )
			{
				wl[ip] = C13O16Rotate[i].WLAng;

				/* put null terminated line label into chLab using line array*/
				chIonLbl(chLab[ip], &C13O16Rotate[i]);
				frac[ip] = (float)(RadPres1/pressure.PresRadCurr);

				ip = MIN2((long)NLINE-1,ip+1);
			}
		}

		/* this works since we forced Lya into stack as first element */
		ASSERT( ip> 0);

		/* sort from strong to weak, then only print strongest */
		spsort(
			/* input array to be sorted */
			frac, 
			/* number of values in x */
			ip, 
			/* permutation output array */
			iperm, 
			/* flag saying what to do - 1 sorts into increasing order, not changing
			* the original routine */
			-1, 
			/* error condition, should be 0 */
			&ier);

		/* this works since we forced Lya into stack as first element */
		ASSERT( ip> 0);

		/* now print six strongest contributors to radiation pressure */
		fprintf( ioQQQ, " P(Lines):" );
		for( i=0; i < MIN2(6,ip); i++ )
		{
			int ipline = iperm[i];
			fprintf( ioQQQ, "(%4.4s ", chLab[ipline]);
			prt_wl(ioQQQ, wl[ipline] );
			fprintf(ioQQQ," %.2f) ",frac[ipline]);
		}

		/* finally end the line */
		fprintf( ioQQQ, "\n" );
	}

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

