/* 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 */
/*PunchLineData punches selected line data for all lines transferred in code */
/*Punch1LineData punch data for one line */
#include "cddefines.h"
#include "taulines.h"
#include "hmi.h"
#include "iso.h"
#include "phycon.h"
#include "physconst.h"
#include "elementnames.h"
#include "hydrogenic.h"
#include "lines_service.h"
#include "dense.h"
#include "atomfeii.h"
#include "lines.h"
#include "atmdat.h"
#include "prt.h"
#include "mole.h"
#include "h2.h"
#include "thermal.h"
#include "cooling.h"
#include "punch.h"

void PunchLineData(FILE * io)
{

	long int i, 
	  j, 
	  limit ,
	  nelem ,
	  ipHi , 
	  ipLo;

	const long nskip=2; /* number of emission lines per line of output */
	double tot;
	int lgElemOff=FALSE;
	float a , b; /* dummy vars to pass to rotate cooling routine */

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

	/* routine punches out (on unit io) line data
	 * for all recombination lines, and all transitions that are transferred */

	/* say what is happening so we know why we stopped */
	fprintf( ioQQQ, " punching line data, then stopping\n" );

	/* first check that all lines are turned on */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		if( !dense.lgElmtOn[nelem] )
		{
			fprintf(ioQQQ," WARNING - I am punching line data but element %s is turned off.\n",
			elementnames.chElementName[nelem]);
			lgElemOff = TRUE;
		}
	}
	if( lgElemOff )
	{
		fprintf(ioQQQ,"Some elements are turned off and punch line data requested.\n");
		fprintf(ioQQQ,"Code is now designed to do punch line data only with all elements on.\n");
		fprintf(ioQQQ,"Please try again with all elements on.\n");
		puts( "[Stop in PunchLineData]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* evaluate rec coef at constant temperature if this is set, else
	 * use 10,000K */
	if( thermal.lgTSetOn )
	{
		phycon.te = thermal.ConstTemp;
	}
	else
	{
		phycon.te = 1e4;
	}
	tfidle(FALSE);

	/* this is set of Dima's recombination lines */
	atmdat_rec_lines(phycon.te,LineSave.RecCoefCNO);
	fprintf( io, "\n       Recombination lines of C, N, O\n" );
	fprintf( io, "    Ion  WL(A)   Coef     Ion   WL(A)  Coef\n" );
	for( i=0; i<471; i+=nskip)
	{
		/* nskip is set to 2 above */
		limit = MIN2(471,i+nskip);
		fprintf( io, "    " );
		for( j=i; j < limit; j++ )
		{
			fprintf( io, "%2.2s%2ld%6ld%8.3f    ", 
				elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][j])-1], 
			  (long)(LineSave.RecCoefCNO[0][j]-LineSave.RecCoefCNO[1][j]+1.01), 
			  (long)(LineSave.RecCoefCNO[2][j]+0.5), 
			  log10(SDIV(LineSave.RecCoefCNO[3][j]) ) );
		}
		fprintf( io, "    \n" );
	}
	fprintf( io, "\n\n" );

	dense.eden = 1.;
	dense.gas_phase[ipHYDROGEN] = 1.;
	dense.EdenHCorr = 1.;

	/* want very small neutral fractions so get mostly e- cs */
	dense.xIonDense[ipHYDROGEN][1] = 1.e-5f;
	hmi.Hmolec[ipMH2g] = 0.;
	dense.xIonDense[ipHYDROGEN][1] = 1.;
	for( i=1; i <= nLevel1; i++ )
	{
		TauLines[i].PopLo = 1.;
	}

	for( i=0; i < nWindLine; i++ )
	{
		TauLine2[i].PopLo = 1.;
	}

	for( i=0; i < nUTA; i++ )
	{
		UTALines[i].PopLo = 1.;
	}

	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < LIMELM+1; j++ )
		{
			dense.xIonDense[i][j] = 1.;
		}
	}

	/* evaluate cooling, this forces evaluation of collision strengths */
	CoolEvaluate(&tot);

	fprintf( io, "       Level 1 transferred lines\n" );

	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );

	for( i=1; i <= nLevel1; i++ )
	{
		/* chLineLbl generates a 1 char string from the line transfer array info -
		 * "Ne 2  128" the string is null terminated,
		 * in following printout the first 4 char is used first, followed by
		 * an integer, followed by the rest of the array*/
		Punch1LineData( &TauLines[i] , io , TRUE);
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end level 1, start level 2\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			Punch1LineData( &TauLine2[i] , io , TRUE);
		}
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end level 2, start inner shell\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );

	for( i=0; i < nUTA; i++ )
	{
		if( UTALines[i].Aul > 0. )
			Punch1LineData( &UTALines[i] , io , TRUE);
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end inner shell, start h-like iso seq\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );

	/* h-like iso sequence */
	/* the hydrogen like iso-sequence */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		HydroCollid( nelem );
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			/* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */
			for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipH_LIKE][nelem]-1; ipLo++ )
			{
				/* set ENTIRE array to impossible values, in case of bad pointer */
				for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
				{
					Punch1LineData( &EmisLines[ipH_LIKE][nelem][ipHi][ipLo] , io , FALSE );
				}
			}
		}
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end h-like iso seq, start he-like iso seq\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );
	for( nelem=1; nelem < LIMELM; nelem++ )
	{
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			/* arrays are dim'd ipHe2p1P */
			for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][nelem]-1; ipLo++ )
			{
				/* set ENTIRE array to impossible values, in case of bad pointer */
				for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
				{
					Punch1LineData( &EmisLines[ipHE_LIKE][nelem][ipHi][ipLo] , io , FALSE );
				}
			}
		}
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end he-like iso seq, start hyperfine structure lines\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );
	/* fine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		Punch1LineData( &HFLines[i] , io , TRUE);
	}

	/* want very small ionized fractions so get mostly H2 cs */
	dense.eden = 1e-6;
	dense.gas_phase[ipHYDROGEN] = 1e-6f;
	dense.EdenHCorr = 1e-6f;
	dense.xIonDense[ipHYDROGEN][1] = 1.;
	hmi.Hmolec[ipMH2g] = 1.;
	hmi.Hmolec[ipMH2s] = 1.;
	dense.xIonDense[ipHYDROGEN][1] = 1e-6f;

	/* H2 molecule */
	fprintf( io, "\n\n\n" );
	fprintf( io, "       end hyperfine, start H2 lines\n" );
	fprintf( io, "Eu Vu Ju El Vl Jl        WL    gl gu    gf       A       CS   n(crt)\n" );

	/* io unit, and option to print all possible lines - FALSE indicates
	 * punch only significant lines */
	H2_LevelPops();
	H2_Punch_line_data( io , FALSE );

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end H2, start 12CO rotation lines\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );
	CO_PopsEmisCool(&C12O16Rotate, nCORotate ,1. , "12CO",&a,&b );
	for( j=0; j< nCORotate; ++j )
	{
		Punch1LineData( &C12O16Rotate[j] , io , TRUE);
	}

	fprintf( io, "\n\n\n" );
	fprintf( io, "       end 12CO start 13CO rotation lines\n" );
	fprintf( io, "Ion   WL  gl gu    gf       A       CS   n(crt)\n" );
	CO_PopsEmisCool(&C13O16Rotate, nCORotate , 1. ,"13CO",&a,&b );
	for( j=0; j< nCORotate; ++j )
	{
		Punch1LineData( &C13O16Rotate[j] , io , TRUE);
	}

	/* punch FeII data if atom is currently enabled */
	fprintf( io, "\n\n\n" );
	fprintf( io, "       end 13CO rotation lines, start FeII lines\n" );
	fprintf( io, " Lo  Hi  Ion  label   WL  gl gu    gf       A       CS   n(crt)\n" );

	/* io unit, and option to print all possible lines - FALSE indicates
	 * punch only significant lines */
	FeIIPunData( io , FALSE );

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

	/* stop when done, we have done some serious damage to the code */
	puts( "[Normal stop in PunchLineData - no ChkAssert to report]" );
	fprintf( ioQQQ, " Cloudy ends, exited OK\n" );
	cdEXIT(EXIT_SUCCESS);
}

/*Punch1LineData punch data for one line */
void Punch1LineData( EmLine * t , FILE * io  , 
	/* flag saying whether to give collision strength too - in multi level atoms
	 * it will be not valid without a great deal more work */
	 int lgCS_2 )
{

	double CritDen;
	/* these are used for parts of the line label */
	char chLbl[11];

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

	/*TODO	1	define lifetime and collision rate for multi-level species so that the
	 * critical density is derived correctly in this routine.  For now the flag lgCS_2
	 * being true means to punch crit den and is only true for two-level systems 
	 * all places where this routine is called with lgCS_2 false need to be fixed */

	/*iWL = iWavLen( t , &chUnits , &chShift );*/
	/* ion label, like C  1 */
	chIonLbl(chLbl , t );
	fprintf(io,"%s\t", chLbl );

	/* this is the second piece of the line label, pick up string after start */

	/* the wavelength */
	prt_wl(io, t->WLAng );

	fprintf( io, " %3ld%3ld",
	  /* lower and upper stat weights */
	  (long)(t->gLo), 
	  (long)(t->gHi) );

	/* oscillator strength */
	fprintf( io,PrintEfmt("%9.2e",  t->gf));

	/* einstein A for transition */
	fprintf( io,PrintEfmt("%9.2e",  t->Aul));

	/* use different formats depending on size of collision strength */
	if( t->cs > 1. )
	{
		fprintf( io, "%7.3f", t->cs );
	}
	else if( t->cs > .01 )
	{
		fprintf( io, "%7.4f", t->cs );
	}
	else if( t->cs > 0.0 )
	{
		fprintf( io, " %.3e", t->cs );
	}
	else
	{
		fprintf( io, "%7.4f", 0. );
	}

	/* now print critical density but only if cs is positive 
	 * >>chng 06 mar 24, add flag lgCS_2 - in multi-level systems do not want
	 * to punch cs since not computed properly */
	if( lgCS_2 && t->cs> 0. )
	{
		CritDen = t->Aul * t->gHi*phycon.sqrte / (t->cs*COLL_CONST);
		CritDen = log10(CritDen);
	}
	else
	{
		CritDen = 0.;
	}
	fprintf( io, "%7.3f\n",CritDen );

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

	return;
}
