/*CoolDima compute cooling due to level 2 lines */
/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
#include "cddefines.h"
#include "taulines.h"
#include "phycon.h"
#include "ionfracs.h"
#include "converge.h"
#include "level2.h"
#include "putcs.h"
#include "mewecoef.h"
#include "tepowers.h"
#include "sexp.h"
/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
double ColStrGBar(EmLine * t );

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

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

	/* 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++ )
		{
			/* only evaluate cs if positive abundance */
			ion = TauLine2[i].IonStg ;
			nelem = TauLine2[i].nelem ;

			if( xIonFracs[ion][nelem-1] > 0. )
			{
				/* now generate the collision strength */
				cs = ColStrGBar(&TauLine2[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 level2 to use existing WineEtc labels */
		level2(&TauLine2[i] );
	}

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

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

double ColStrGBar(EmLine * t )
{
	long int i, 
	  j;
	double ColStrGBar_v, 
	  a, 
	  b, 
	  c, 
	  d, 
	  e1, 
	  gb, 
	  ttype, 
	  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 Fisher et al. (1996)
	 * Gaetz & Salpeter (1983, ApJS 52, 155) and 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? */
	ttype = t->cs1;

	/* >>chng 99 feb 27, change to assert */
	assert( ttype >= 0.05 );
#if 0
	if( ttype < 0.05 )
	{
		fprintf( ioQQQ, " ColStrGBar called with insane cs1\n" );
		puts( "[Stop in colstrgbar]" );
		exit(1);
	}
#	endif

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

		if( ttype < 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->EnergyRyd* 13.6);
	}
	else
	{
		i = (long int)ttype;

		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->EnergyRyd;
			/* 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((double)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((double)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*TePowers.sqrte*TePowers.te10*
				  TePowers.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 );
}