/*lines main routine to put emission line intensities into line stack */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "linesave.h"
#include "ionfracs.h"
#include "trace.h"
#include "dhla.h"
#include "compton.h"
#include "hmihet.h"
#include "tcool.h"
#include "opac.h"
#include "heat.h"
#include "rfield.h"
#include "phycon.h"
#include "nh.h"
#include "extramax.h"
#include "holod.h"
#include "elementnames.h"
#include "reccno.h"
#include "hphoto.h"
#include "hfbnet.h"
#include "linadd.h"
#include "lindst.h"
#include "linefit.h"
#include "rdsum.h"
#include "lines.h"
#include "radius.h"
#include "pntforline.h"
#include "putline.h"

void lines(void)
{
	char chLabel[5];
	static int lgRecOn;
	long int i, 
	  ipnt,
	  ipZ;
	double BigstExtra, 
	  ExtraCool, 
	  bac, 
	  f1, 
	  f2; 
	double sum;

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

	/* major routines used here:
	 *
	 * PutLine( tarray )
	 * this uses information in tarray to give various
	 * contributions to lines, and their intensities
	 *
	 * PutExtra( extra )
	 * extra is some extra intensity to add to the line
	 * it will go into the totl contribution put out by PutLine,
	 * and this contriibution should be indicated by independent
	 * call to linadd
	 * */

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

	heat.power += heat.htot*radius.dVeff;
	dhla.dstpow += dhla.dstotl*radius.dVeff;

	/* total compton heating - cooling */
	Compton.comtot += Compton.cmheat*radius.dVeff;
	hmihetCom.hmitot += hmihetCom.hmihet*radius.dVeff;
	tcool.totcol += tcool.ctot*radius.dVeff;
	tcool.GrossHeat += holod.feheat*radius.dVeff;

	/* up up induced recom cooling */
	for( ipZ=0; ipZ<LIMELM; ++ipZ )
	{
		HPhoto.cintot += HFBNet.crcind[ipZ]*(float)radius.dVeff;
	}

	/* nsum is line pointer within large stack of line intensities */
	LineSave.nsum = 0;
	LineSave.ndsum = 0;

	/* last arg in call to lindst and linadd is info on line
	 * info is char variable indicating type of line this is
	 * 'c' cooling
	 * 'h' heating
	 * 'i' information only
	 * 'r' recombination line
	 *
	 * all components of lines are entered into the main line stack here
	 * when printing, filters exist to not print Inwd component */
	 
	/* initialize routine that can generate pointers for forbidden lines  */
	PntForLine(0.,"FILL",&i);

	/* evaluate rec coef for rec lines of C, N, O
	 * some will be used in LineSet2 and then zeroed out,
	 * others left alone and used below */
	LineFit(phycon.te,RecCNO.RecCoefCNO);

	/* punch energetics, H and He lines */
	LineSet1();

	/* do heavies, lithium through neon */
	LineSet2();

	/* do heavies, sodium through argon */
	LineSet3();

	/* do heavies, potasium through zinc */
	LineSet4();

	rdsum(&sum," SUM");
	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* optional sum of certain emission lines, set with "print sum" */
	linadd(sum/radius.dVeff,0,"Stoy",'i');

	/***********************************************************************
	 * stuff in Bac ratio - continuum above the BJ 
	 * this is trick, zeroing out saved continuum integrated so far,
	 * and adding the current version, so that the line array gives the 
	 * value in the final continuum 
	 *
	 * reflec continuum is different from others since relative to inner
	 * radius, others for for this radius 
	 *************************************************************************/

	/***************************************************************************
	 * "Bac " , 3646,  this is total continuum at peak of Balmer Jump          *
	 ***************************************************************************/
	f1 = (rfield.outcon[nh.ipHn[0][2]-1] + 
		rfield.ConOutNoInter[nh.ipHn[0][2]-1] + 
		rfield.refcon[nh.ipHn[0][2]-1]/radius.r1r0sq )/
	  rfield.widflx[nh.ipHn[0][2]-1]/opac.tmn[nh.ipHn[0][2]-1]/radius.GeoDil;

	/* extrapolated continuum below head */
	f2 = (rfield.outcon[nh.ipHn[0][2]-2]+ 
		rfield.ConOutNoInter[nh.ipHn[0][2]-2] + 
		rfield.refcon[nh.ipHn[0][2]-2]/radius.r1r0sq )/
	  rfield.widflx[nh.ipHn[0][2]-2]/opac.tmn[nh.ipHn[0][2]-2]/radius.GeoDil;

	/* convert to nuFnu units */
	bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;
	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* residual flux at head of Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"Bac ",'i');

	/*********************************************************************
	 * "cout" , 3646,  this is outward continuum at peak of Balmer Jump  *
	 * equal to total in spherical geometry, half in opt thin open geo   *
	 *********************************************************************/
	f1 = (rfield.outcon[nh.ipHn[0][2]-1]+ rfield.ConOutNoInter[nh.ipHn[0][2]-1])/
		rfield.widflx[nh.ipHn[0][2]-1]/opac.tmn[nh.ipHn[0][2]-1]/radius.GeoDil;

	f2 = (rfield.outcon[nh.ipHn[0][2]-2]+ rfield.ConOutNoInter[nh.ipHn[0][2]-2])/
		rfield.widflx[nh.ipHn[0][2]-2]/opac.tmn[nh.ipHn[0][2]-2]/radius.GeoDil;

	/* net balmer jump */
	bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"cout",'i');

	/*********************************************************************
	 * "cref" , 3646,  this is reflected continuum at peak of Balmer Jump*
	 * equal to zero in spherical geometry, half of total in op thin opn *
	 *********************************************************************/
	f1 = rfield.refcon[nh.ipHn[0][2]-1]/
		rfield.widflx[nh.ipHn[0][2]-1]/opac.tmn[nh.ipHn[0][2]-1]/radius.GeoDil;

	f2 = rfield.refcon[nh.ipHn[0][2]-2]/
		rfield.widflx[nh.ipHn[0][2]-2]/opac.tmn[nh.ipHn[0][2]-2]/radius.GeoDil;

	/* net balmer jump */
	bac = (f1 - f2)*0.250*0.250*EN1RYD;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"cref",'i');

	/*********************************************************************
	 * "thin" , 3646, tot optically thin continuum at peak of Balmer Jump*/
	f1 = rfield.DiffLocal[nh.ipHn[0][2]-1]/rfield.widflx[nh.ipHn[0][2]-1];
	f2 = rfield.DiffLocal[nh.ipHn[0][2]-2]/rfield.widflx[nh.ipHn[0][2]-2];
	bac = (f1 - f2)*0.250*0.250*EN1RYD;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac),3646,"thin",'i');


	/***********************************************************************
	 * large number of C, N, and O recombination lines                     *
	 *************************************************************************/

	if( LineSave.ipass <= 0 )
	{
		/* initialize to TRUE, may be set false if density is very high below */
		lgRecOn = TRUE;
	}

	for( i=0; i < 471; i++ )
	{
		/* generate label for the line */
		strcpy( chLabel, elementnames.chElementSym[(long)(RecCNO.RecCoefCNO[0][i])-1] );
		strcat( chLabel, elementnames.chIonStage[(long)(RecCNO.RecCoefCNO[0][i]-
			RecCNO.RecCoefCNO[1][i]+1.01)-1] );

		/* number of rec per unit vol
		 * >>chng 97 aug 28, do not predict rec intensities at high densities
		 * blr.in and oldblr.in go nuts if this is added in outward beam */
		if( phycon.eden < 1e8 && lgRecOn )
		{
			f2 = RecCNO.RecCoefCNO[3][i]*phycon.eden*xIonFracs[(long)(RecCNO.RecCoefCNO[0][i]-
			  RecCNO.RecCoefCNO[1][i]+2)][(long)(RecCNO.RecCoefCNO[0][i])-1];

			/* convert to intensity */
			f2 = f2*1.99e-8/RecCNO.RecCoefCNO[2][i];
		}
		else
		{
			lgRecOn = FALSE;
			f2 = 0.;
		}
		/* finally stuff it into the stack */
		PntForLine(RecCNO.RecCoefCNO[2][i],chLabel,&ipnt);
		lindst(f2,(long)(RecCNO.RecCoefCNO[2][i]+.5),chLabel,ipnt, 'r',TRUE);
	}

	/* add in all the other level 2 wind lines
	 * Dima'a 6k lines */
	ExtraCool = 0.;
	BigstExtra = 0.;
	for( i=0; i < nWindLine; i++ )
	{
		PutLine(&TauLine2[i]);
		if( TauLine2[i].cool > BigstExtra )
		{
			BigstExtra = TauLine2[i].cool;
			ExtraMax.ipMaxExtra = i+1;
		}
		ExtraCool += TauLine2[i].cool;
	}
	ExtraMax.GBarMax = (float)MAX2(ExtraMax.GBarMax,ExtraCool/tcool.ctot);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " lines returns\n" );
	}

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