/* 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 */
/*CoolDima compute cooling due to level 2 lines */
/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "dense.h"
#include "phycon.h"
#include "converge.h"
#include "thermal.h"
#include "opacity.h"
#include "lines_service.h"
#include "rfield.h"
#include "mewecoef.h"
#include "atoms.h"
#include "cooling.h"

/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
static double ColStrGBar(EmLine * t , float cs1 );

void CoolDima(void)
{
	long int i, 
	  ion,
	  nelem;
	double cs;
	double OTSLevel2,
		sumots,
		CoolLevel2;
	static float TeEvalCS = -1.f;
	static long int nzoneEval=0;

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

	/* >>chng 03 nov 29, add this option to skip rest of routine */
	/* no level2 command sets nWindLine to -1 */
	if( nWindLine<0 )
		return;

	/* only force evaluation of cooling when significant, or first call
	 * in this zone, or not into first zone */
	if( nzone != nzoneEval || conv.lgLevel2_Cool_Imp || !nzone )
	{
		nzoneEval = nzone ;

		/* check whether we need to reevaluate the collision strenghts.
		 * Do so if large change in temperature  or any stage of ionization has been lowered. */
		if( (conv.lgSearch || conv.lgIonStageTrimed || 
			fabs(phycon.te-TeEvalCS)/phycon.te > 0.05) )
		{
			for( i=0; i < nWindLine; i++ )
			{
				if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
				{
					/* only evaluate cs if positive abundance */
					ion = TauLine2[i].IonStg ;
					nelem = TauLine2[i].nelem ;

					if( dense.xIonDense[nelem-1][ion-1] > 0. )
					{
						/* now generate the collision strength */
						cs = ColStrGBar(&TauLine2[i] , cs1_flag_lev2[i] );
					}
					else
					{
						cs = 1.;
					}
					/* now put the cs into the line array */
					PutCS(cs,&TauLine2[i] );
				}
			}
			TeEvalCS = phycon.te;
		}

		/* now always update the cooling since this adds ots flux */
		for( i=0; i < nWindLine; i++ )
		{
			/* we need to call this even if nothing present since then sets zero
			 * ignore all lines with negative atomic number
			 * >>chng 97 aug 30, get rid of test for non-pos ipLnNelem
			 * pointer tells atom_level2 to use existing WineEtc labels */
			/* only call non-hydrogenic, non-he-like lines */
			 /* >>chng 02 sug 11, do not include he-like in sum */
			if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
			{
				atom_level2(&TauLine2[i] );
			}
		}
	}
	else
	{
		/* now use old cooling and ots flux */
		for( i=0; i < nWindLine; i++ )
		{
			if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
			{
				/* recall that these should be trivial cooling, ots rates */
				CoolAdd( "    ", TauLine2[i].WLAng , TauLine2[i].cool);
				/*>>chng 03 oct 04, move to RT_OTS */
				/*RT_OTS_AddLine( TauLine2[i].ots , TauLine2[i].ipCont);*/
			}
		}
	}

	/* find ratio of level 1 to level 2 ots to see how important level 2 is.
	 * if level 2 ots not important then will update dest probs only first time
	 * in each zone.  If significant then continuously update */
	OTSLevel2 = 0.;
	CoolLevel2 = 0.;
	if( nzone > 1 )
	{
		for( i=0; i < nWindLine; i++ )
		{
			long int ip = TauLine2[i].ipCont-1;
			CoolLevel2 += TauLine2[i].cool;
			if( opac.opacity_abs[ip] > SMALLFLOAT )
			{
				OTSLevel2 += TauLine2[i].ots / ( (float)opac.opacity_abs[ip]) ;
				ASSERT( fabs(OTSLevel2) < 1e37 );
			}
		}

		/* now sum all ots rates */
		sumots = 0.;
		for( i=0; i < rfield.nflux; i++ )
		{
			sumots += rfield.otslin[i];
		}
		/* is the ots contribution significant ? */
		if( OTSLevel2/SDIV(sumots) > 1e-4 )
		{
			/* true if level 2 lines were contributors to the ots rates, set in dimacool */
			conv.lgLevel2_OTS_Imp = TRUE;
		}
		else
		{
			/* true if level 2 lines were contributors to the ots rates, set in dimacool */
			conv.lgLevel2_OTS_Imp = FALSE;
		}

		/* set flags saying if important */
		if( CoolLevel2/SDIV(thermal.ctot) > 1e-4 )
		{
			/* true if level 2 lines were contributors to the cooling, set in dimacool */
			conv.lgLevel2_Cool_Imp = TRUE;
		}
		else
		{
			/* true if level 2 lines were contributors to the cooling, set in dimacool */
			conv.lgLevel2_Cool_Imp = FALSE;
		}

		{
			/* debugging code for Lya problems */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				fprintf(ioQQQ,"%li\t%.2e\t%.2e\n", 
					nzone, OTSLevel2/SDIV(sumots)  , CoolLevel2/SDIV(thermal.ctot) );
			}
		}
	}
	else
	{
		/* true if level 2 lines were contributors to the ots rates, set in dimacool */
		conv.lgLevel2_OTS_Imp = TRUE;
		/* true if level 2 lines were contributors to the cooling, set in dimacool */
		conv.lgLevel2_Cool_Imp = TRUE;
	}

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

/*ColStrGBar generate g-bar collision strengths for level 2 line2 */

static double ColStrGBar(EmLine * t , float cs1 )
{
	long int i, 
	  j;
	double ColStrGBar_v, 
	  a, 
	  b, 
	  c, 
	  d, 
	  e1, 
	  gb, 
	  x, 
	  y;
	double xx, 
	  yy;

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

	/* Calculation of the collision strengths of multiplets.
	 * Neutrals are recalculated from 
	 * >>refer cs	gbar	Fisher et al. (1996)
	 * >>refer cs	gbar	Gaetz & Salpeter (1983, ApJS 52, 155) and 
	 * >>refer cs	gbar	Mewe (1972, A&A 20, 215) 
	 * fits for ions. */

	/* zero hydrogenic lines since they are done by iso-sequence */
	if( t->nelem == t->IonStg )
	{
		ColStrGBar_v = 0.0;
		
#		ifdef DEBUG_FUN
		fputs( " <->ColStrGBar()\n", debug_fp );
#		endif
		return( ColStrGBar_v );
	}

	/*was the block data linked in? */
	ASSERT( MeweCoef.g[1][0] != 0.);

	/* which type of transition is this? cs1 is the flag */

	/* >>chng 01 may 30 - cs1 < 0 means a forced collision strength */
	if( cs1 < 0. )
	{
		ColStrGBar_v = -cs1;
		
#		ifdef DEBUG_FUN
		fputs( " <->ColStrGBar()\n", debug_fp );
#		endif
		return( ColStrGBar_v );
	}

	/* >>chng 99 feb 27, change to assert */
	ASSERT( cs1 >= 0.05 );

	/* excitation energy over kT */
	y = t->EnergyK/phycon.te;
	if( cs1 < 1.5 )
	{
		xx = -log10(y);

		if( cs1 < 0.5 )
		{
			yy = (1.398813573838321 + xx*(0.02943050869177121 + xx*
			  (-0.4439783893114510 + xx*(0.2316073358577902 + xx*(0.001870493481643103 + 
			  xx*(-0.008227246351067403))))))/(1.0 + xx*(-0.6064792600526370 + 
			  xx*(0.1958559534507252 + xx*(-0.02110452007196644 + 
			  xx*(0.01348743933722316 + xx*(-0.0001944731334371711))))));
		}

		else
		{
			yy = (1.359675968512206 + xx*(0.04636500015069853 + xx*
			  (-0.4491620298246676 + xx*(0.2498199231048967 + xx*(0.005053803073345794 + 
			  xx*(-0.01015647880244268))))))/(1.0 + xx*(-0.5904799485819767 + 
			  xx*(0.1877833737815317 + xx*(-0.01536634911179847 + 
			  xx*(0.01530712091180953 + xx*(-0.0001909176790831023))))));
		}

		ColStrGBar_v = pow((double)10,yy)*t->gf/(t->EnergyWN * WAVNRYD* 13.6);
	}
	else
	{
		i = (long int)cs1;

		if( i < 26 )
		{
			e1 = log(1.0+1.0/y) - 0.4/POW2(y + 1.0);
			a = MeweCoef.g[i-1][0];
			b = MeweCoef.g[i-1][1];
			c = MeweCoef.g[i-1][2];
			d = MeweCoef.g[i-1][3];
			x = (double)t->nelem - 3.0;

			if( i == 14 )
			{
				a *= 1.0 - 0.5/x;
				b = 1.0 - 0.8/x;
				c *= 1.0 - 1.0/x;
			}

			else if( i == 16 )
			{
				a *= 1.0 - 0.9/x;
				b *= 1.0 - 1.7/x;
				c *= 1.0 - 2.1/x;
			}

			else if( i == 18 )
			{
				a *= 1.0 + 2.0/x;
				b *= 1.0 - 0.7/x;
			}

			gb = a + (b*y - c*y*y + d)*e1 + c*y;

			/*  ipLnRyd is exciation energy in Rydbergs */
			ColStrGBar_v = 14.510395*t->gf*gb/(t->EnergyWN * WAVNRYD);
			/* following i>=26 */
		}

		else
		{
			/* 210 is the dimem of g, so [209] is largest val */
			if( i < 210 )
			{
				j = (long)(MeweCoef.g[i-1][3]);
				if( j == 1 )
				{
					ColStrGBar_v = t->gLo*MeweCoef.g[i-1][0]*
					  pow(phycon.te/pow((float)10,MeweCoef.g[i-1][2]),MeweCoef.g[i-1][1]);
				}
				else
				{
					ColStrGBar_v = t->gLo*MeweCoef.g[i-1][0]*
					  sexp(MeweCoef.g[i-1][1]*(pow((float)10,MeweCoef.g[i-1][2])/
					  phycon.te));
				}
			}

			else
			{
				/* This is for AlII 1670 line only!
				 *    ColStrGBar=0.0125*te**0.603 */
				/* 98 dec 27, this is still in use */
				ColStrGBar_v = 0.0125*phycon.sqrte*phycon.te10*
				  phycon.te003;
			}
		}
	}

	/* following to make sure that negative values not returned */
	ColStrGBar_v = MAX2(ColStrGBar_v,1e-10);

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