/*PrtLinePres print line radiation pressures for current conditions */
#include <string.h>
#include <float.h>
#include <limits.h>
#include "cddefines.h"
#include "powi.h"
#include "pressure.h"
#include "taulines.h"
#include "chionlbl.h"
#include "iwavlen.h"
#include "widthline.h"
#include "prtlinepres.h"
/* this is limit to number of lines we can store at one time */
#define	NLINE	20
/* faintest pressure we will bother with */
#define	THRESH	0.05

void PrtLinePres(void)
{
	char chLab[NLINE][5];
	long int i, 
	  ip, 
	  iwl[NLINE];
	double RadPres1, 
	  frac[NLINE];

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

	/* this routine only called if printout of contributors to line
	 * radiation pressure is desired
	 *
	 * plines(1)=l-alpha,          in hlevel(ionize)
	 * plines(2) = hign n lyman ;   plines(3)=excit h lines iN HLEVEL
	 * plines(4) = he i 10830
	 * plines(5)= HeII 304 plines(11)=HeI La (in helium) 
	 * plines(6)= HeI 584 */
	/* compute radiation pressure due to special H and He lines */
	ip = 0;

	for(i=0; i<NLINE; ++i)
	{
		frac[i] = DBL_MAX;
		iwl[i] = LONG_MAX;
	}

	if( pressure.RadPres > 1e-30 )
	{
		for( i=0; i < NLINE; i++ )
		{
			if( pressure.plines[i] > pressure.RadPres*THRESH )
			{

				if( i == 0 )
				{
					iwl[ip] = 1216;
					strcpy( chLab[ip], "H  1" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}

				else if( i == 1 )
				{
					iwl[ip] = 1034;
					strcpy( chLab[ip], "H  1" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}
				else if( i == 2 )
				{
					iwl[ip] = 6563;
					strcpy( chLab[ip], "H  1" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}
				else if( i == 3 )
				{
					iwl[ip] = 10830;
					strcpy( chLab[ip], "He 1" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}
				else if( i == 4 )
				{
					iwl[ip] = 304;
					strcpy( chLab[ip], "He 2" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}
				else if( i == 5 )
				{
					iwl[ip] = 584;
					strcpy( chLab[ip], "He 1" );
					frac[ip] = pressure.plines[i]/pressure.RadPres;
				}
				else
				{
					fprintf( ioQQQ, " PrtLinePres insanity\n" );
					puts( "[Stop in prtlinepres]" );
					exit(1);
				}
				++ip;
			}
		}

		/* 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)*
				  WidthLine(&TauLines[i])*2.;
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.RadPres*THRESH )
			{
				iwl[ip] = iWavLen(&TauLines[i]);

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

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

		for( i=0; i < nWindLine; i++ )
		{
			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)*
				  WidthLine(&TauLine2[i])*2.;
			}
			else
			{
				RadPres1 = 0.;
			}

			if( RadPres1 > pressure.RadPres*THRESH )
			{
				iwl[ip] = iWavLen(&TauLine2[i]);

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

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

		/* now print them */
		fprintf( ioQQQ, " P(Lines):" );
		for( i=0; i < ip; i++ )
		{
			fprintf( ioQQQ, "%8.2f %4.4s%5ld,", 
				frac[i], chLab[i], iwl[i] );
		}

		/* possible that dominant contributor not in above, in which case ip is 0
		 * and we will prin contributor that came from TotalPressure */
		if( ip <1 )
		{
			fprintf( ioQQQ, " %s ", pressure.chLineRadPres);
		}
		fprintf( ioQQQ, "\n" );
	}

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