/*level3 compute three level atom, 10, 21, and 20 are line */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "taulines.h"
#include "showme.h"
#include "ionfracs.h"
#include "edsqte.h"
#include "teinv.h"
#include "phycon.h"
#include "tcool.h"
#include "poplevls.h"
#include "lev3fail.h"
#include "chionlbl.h"
#include "insane.h"
#include "addotslin.h"
#include "coladd.h"
#include "chlinelbl.h"
#include "level3.h"

void level3(EmLine * t10, 
  EmLine * t21, 
  EmLine * t20)
{
	char chLab[5], 
	  chLab10[11];
	double AbunxIon, 
	  a, 
	  a10, 
	  a20, 
	  a21, 
	  b, 
	  beta, 
	  bolt01, 
	  bolt02, 
	  bolt12, 
	  c, 
	  ener10, 
	  ener20, 
	  ener21, 
	  g0, 
	  g010, 
	  g020, 
	  g1, 
	  g110, 
	  g121, 
	  g2, 
	  g220, 
	  g221, 
	  o10, 
	  o20, 
	  o21, 
	  p0, 
	  p1, 
	  p2, 
	  pump01, 
	  pump02, 
	  pump12;

	double TotCool, 
	  TotHeat, 
	  TotInten, 
	  alpha, 
	  alpha1, 
	  alpha2, 
	  c01, 
	  c02, 
	  c10, 
	  c12, 
	  c20, 
	  c21, 
	  cnet01, 
	  cnet02, 
	  cnet12, 
	  cool01, 
	  cool02, 
	  cool12, 
	  heat10, 
	  heat20, 
	  heat21, 
	  hnet01, 
	  hnet02, 
	  hnet12, 
	  pump10, 
	  pump20, 
	  pump21, 
	  r01, 
	  r02, 
	  r10, 
	  r12, 
	  r20, 
	  r21, 
	  temp01, 
	  temp02, 
	  temp12;

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

	/* compute three level atom, 10, 21, and 20 are line
	 * arrays for 10, 21, and 20 transitions.
	 * one can be a dummy */
	/* >>chng 96 dec 6, to double precision due to round off problems below */

	/* generalized three level atom for any ion
	 * sum of three levels normalized to total abundance
	 *
	 * stat weights of all three lines
	 * sanity check will confirm ok */
	g010 = t10->gLo;
	g110 = t10->gHi;

	g121 = t21->gLo;
	g221 = t21->gHi;

	g020 = t20->gLo;
	g220 = t20->gHi;

	/* these are statistical weights */
	if( g010 > 0. )
	{
		g0 = g010;
	}

	else if( g020 > 0. )
	{
		g0 = g020;
	}

	else
	{
		strcpy( chLab10, chLineLbl(t10) );
		fprintf( ioQQQ, " level3: insane stat weights g0 :%10.2e%10.2e %10.10s\n", 
		  g010, g020, chLab10 );
		insane();
		ShowMe();
		exit(1);
	}

	if( g110 > 0. )
	{
		g1 = g110;
	}

	else if( g121 > 0. )
	{
		g1 = g121;
	}

	else
	{
		strcpy( chLab10, chLineLbl(t10) );
		fprintf( ioQQQ, " level3: insane stat weights g1 :%10.2e%10.2e %10.10s\n", 
		  g010, g020, chLab );
		insane();
		ShowMe();
		puts( "[Stop in level3]" );
		exit(1);
	}

	if( g220 > 0. )
	{
		g2 = g220;
	}

	else if( g221 > 0. )
	{
		g2 = g221;
	}

	else
	{
		strcpy( chLab10, chLineLbl(t20) );
		fprintf( ioQQQ, " level3: insane stat weights g2 :%10.2e%10.2e %10.10s\n", 
		  g010, g020, chLab10 );
		insane();
		ShowMe();
		puts( "[Stop in level3]" );
		exit(1);
	}

	/* abundances from the elements grid
	 * one of these must be a true line */
	if( g010 > 0. )
	{
	/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, t10);
		AbunxIon = xIonFracs[t10->IonStg][ t10->nelem -1];
	}

	else if( g121 > 0. )
	{
	/* put null terminated line label into chLab using line array*/
		chIonLbl(chLab, t21);
		AbunxIon = xIonFracs[t21->IonStg][t21->nelem -1];
	}

	else
		/* this cannot possibly happen */
	{
		fprintf( ioQQQ, " level3: insanity at g010 g121 branch \n" );
		insane();
		ShowMe();
		puts( "[Stop in level3]" );
		exit(1);
	}

	a = t10->EnergyK*teinvCom.teinv;
	b = t21->EnergyK*teinvCom.teinv;
	c = t20->EnergyK*teinvCom.teinv;

	if( c == 0. )
	{
		c = a + b;
	}

	/* if still neg at end, then success!, in common so possible to
	 * to check why zero returned */
	Lev3Fail.nLev3Fail = -1;
	if( AbunxIon <= 1e-30 || c > 60. )
	{
		Lev3Fail.nLev3Fail = 0;

		/* all populations are zero */
		PopLevls.PopLevels[0] = (float)AbunxIon;
		PopLevls.PopLevels[1] = 0.;
		PopLevls.PopLevels[2] = 0.;

		/* 666 error! these pops ARE NOT defined below */
		PopLevls.DepLTELevels[0] = 1.;
		PopLevls.DepLTELevels[1] = 0.;
		PopLevls.DepLTELevels[2] = 0.;

		/* level populations */
		t21->PopOpc = 0.;
		t10->PopOpc = (float)AbunxIon;
		t20->PopOpc = (float)AbunxIon;
		t21->PopLo = 0.;
		t10->PopLo = (float)AbunxIon;
		t20->PopLo = (float)AbunxIon;
		t21->PopHi = 0.;
		t10->PopHi = 0.;
		t20->PopHi = 0.;

		/* line heating */
		t20->heat = 0.;
		t21->heat = 0.;
		t10->heat = 0.;

		/* intensity of line */
		t21->xIntensity = 0.;
		t10->xIntensity = 0.;
		t20->xIntensity = 0.;

		/* line cooling */
		t20->cool = 0.;
		t21->cool = 0.;
		t10->cool = 0.;

		/* local ots rates */
		t20->ots = 0.;
		t21->ots = 0.;
		t10->ots = 0.;

		/* number of photons in line zero */
		t21->phots = 0.;
		t10->phots = 0.;
		t20->phots = 0.;

		/* fraction that produced lines zero */
		t21->AovTot = 1.;
		t10->AovTot = 1.;
		t20->AovTot = 1.;

		/* ratio coll over total excitation */
		t21->ColOvTot = 0.;
		t10->ColOvTot = 0.;
		t20->ColOvTot = 0.;

		/* add zero to cooling */
		coladd(chLab, (long)t21->WLAng ,0.);
		coladd(chLab, (long)t10->WLAng ,0.);
		coladd(chLab, (long)t20->WLAng ,0.);

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

	/* collision strengths */
	o10 = t10->cs;
	o21 = t21->cs;
	o20 = t20->cs;

	/* begin sanity checks, check statistic weights, 
	 * first check is protection against dummy lines */
	assert( (g010*g020 == 0.) || g010 == g020 );

	assert( (g110*g121 == 0.) || g110 == g121 );

	assert( (g221*g220 == 0.) || g221 == g220 );

	/* both abundances must be same, equal abundance
	 * xIonFracs(nelem,i) is density of ith ionization stage (cm^-3) */
	assert( (t10->IonStg*t21->IonStg == 0) || (t10->IonStg == t21->IonStg) );

	assert( (t20->IonStg*t21->IonStg == 0) || (t20->IonStg == t21->IonStg ) );

	assert( (t10->nelem * t21->nelem == 0) || (t10->nelem == t21->nelem) );

	assert( (t20->nelem * t21->nelem == 0) || (t20->nelem == t21->nelem) );

	assert( o10 > 0. && o21 > 0. && o20 > 0. );

	/*end sanity checks */

	/* net loss of line escape and destruction */
	a21 = t21->Aul * (t21->Pesc + t21->Pdest);
	a10 = t10->Aul * (t10->Pesc + t10->Pdest);
	a20 = t20->Aul * (t20->Pesc + t20->Pdest);

	/* find energies of all transitions - one line could be a dummy
	 * also find boltzmann factors */
	if( a10 == 0. )
	{
		ener20 = t20->EnergyErg;
		ener21 = t21->EnergyErg;
		ener10 = ener20 - ener21;
		bolt12 = exp(-t21->EnergyK/phycon.te);
		bolt02 = exp(-t20->EnergyK/phycon.te);
		bolt01 = bolt02/bolt12;
		temp12 = t21->EnergyK;
		temp02 = t20->EnergyK;
		temp01 = temp02 - temp12;
	}

	else if( a21 == 0. )
	{
		ener10 = t10->EnergyErg;
		ener20 = t20->EnergyErg;
		ener21 = ener20 - ener10;
		bolt01 = exp(-t10->EnergyK/phycon.te);
		bolt02 = exp(-t20->EnergyK/phycon.te);
		bolt12 = bolt02/bolt01;
		temp02 = t20->EnergyK;
		temp01 = t10->EnergyK;
		temp12 = temp02 - temp01;
	}

	else if( a20 == 0. )
	{
		ener10 = t10->EnergyErg;
		ener21 = t21->EnergyErg;
		ener20 = ener21 + ener10;
		bolt01 = exp(-t10->EnergyK/phycon.te);
		bolt12 = exp(-t21->EnergyK/phycon.te);
		bolt02 = bolt01*bolt12;
		temp01 = t10->EnergyK;
		temp12 = t21->EnergyK;
		temp02 = temp01 + temp12;
	}

	else
	{
		/* all lines are ok */
		ener10 = t10->EnergyErg;
		ener20 = t20->EnergyErg;
		ener21 = t21->EnergyErg;
		bolt01 = exp(-t10->EnergyK/phycon.te);
		bolt12 = exp(-t21->EnergyK/phycon.te);
		bolt02 = bolt01*bolt12;
		temp02 = t20->EnergyK;
		temp01 = t10->EnergyK;
		temp12 = t21->EnergyK;
	}

	/* check all energies positive */
	assert( ener10 > 0. && ener20 > 0. && ener21 > 0. );

	/* check if energy order is ok */
	assert( ener10 < ener20 && ener21 < ener20 );

	/* check if energy scale is ok */
	assert( fabs((ener10+ener21)/ener20-1.) < 1e-4 );

	pump01 = t10->pump;
	pump10 = pump01*g0/g1;
	pump12 = t21->pump;
	pump21 = pump12*g1/g2;
	pump02 = t20->pump;
	pump20 = pump02*g0/g2;

	/* cdsqte is 8.629E-6 / sqrte * eden */
	c01 = o10*bolt01*edsqte.cdsqte/g0;
	r01 = c01 + pump01;
	c10 = o10*edsqte.cdsqte/g1;
	r10 = c10 + a10 + pump10;
	c20 = o20*edsqte.cdsqte/g2;
	r20 = c20 + a20 + pump20;
	c02 = o20*bolt02*edsqte.cdsqte/g0;
	r02 = c02 + pump02;
	c12 = o21*bolt12*edsqte.cdsqte/g1;
	r12 = c12 + pump12;
	c21 = o21*edsqte.cdsqte/g2;
	r21 = c21 + a21 + pump21;

	alpha1 = (double)(AbunxIon)*(r01+r02)/(r10+r01+r02);
	alpha2 = (double)(AbunxIon)*(r01)/(r10+r12+r01);
	alpha = alpha1 - alpha2;

	/*  1( DBLE(r01+r02)/DBLE(r10+r01+r02) - DBLE(r01)/DBLE(r10+r12+r01) )
	 * beta is factor with n2 */
	beta = (r21 - r01)/(r10 + r12 + r01) + (r20 + r01 + r02)/(r10 + 
	  r01 + r02);

	if( alpha/MAX2(alpha1,alpha2) < 1e-11 )
	{
		/* this catches both negative and round off */
		p2 = 0.;
		alpha = 0.;
		Lev3Fail.nLev3Fail = 1;
	}

	else
	{
		p2 = alpha/beta;
	}
	PopLevls.PopLevels[2] = (float)p2;

	if( alpha < 0. || beta < 0. )
	{
		fprintf( ioQQQ, " level3: insane n2 pop alf, bet, p2=%10.2e%10.2e%10.2e %10.10s t=%10.2e\n", 
		  alpha, beta, p2, chLab, phycon.te );
		fprintf( ioQQQ, " gs are%5.1f%5.1f%5.1f\n", g0, g1, 
		  g2 );
		fprintf( ioQQQ, " Bolts are%10.2e%10.2e%10.2e\n", 
		  bolt01, bolt12, bolt02 );
		fprintf( ioQQQ, " As are%10.2e%10.2e%10.2e\n", a10, 
		  a21, a20 );
		fprintf( ioQQQ, " Energies are%10.2e%10.2e%10.2e\n", 
		  ener10, ener21, ener20 );
		fprintf( ioQQQ, " 2 terms, dif of alpha are%15.6e%15.6e\n", 
		  (r01 + r02)/(r10 + r01 + r02), r01/(r10 + r12 + r01) );
		ShowMe();
		puts( "[Stop in level3]" );
		exit(1);
	}

	alpha = (double)(AbunxIon)*(r01+r02) - (double)(p2)*(r20+r01+r02);
	/* guard against roundoff - this should really have been zero
	 * >>chng 96 nov 14, protection against round-off to zero
	 * >>chng 96 dec 3, made r01, etc, double, and changed limit to 1e-9 */
	if( fabs(alpha)/(MAX2(AbunxIon*(r01+r02),p2*(r20+r01+r02))) < 1e-9 )
	{
		alpha = 0.;
		Lev3Fail.nLev3Fail = 2;
	}

	beta = r10 + r01 + r02;
	p1 = alpha/beta;
	PopLevls.PopLevels[1] = (float)p1;

	if( p1 < 0. )
	{
		if( p1 > -1e-37 )
		{
			/* slightly negative solution, probably just round-off, zero it */
			p1 = 0.;
			PopLevls.PopLevels[1] = (float)p1;
			Lev3Fail.nLev3Fail = 3;
		}

		else
		{
			/* very negative solution, better punt */
			fprintf( ioQQQ, " level3: insane n1 pop alf, bet, p1=%10.2e%10.2e%10.2e %10.10s%5f\n", 
			  alpha, beta, p1, chLab,  t10->WLAng );
			fprintf( ioQQQ, " local electron density and temperature were%10.2e%10.2e\n", 
			  phycon.eden, phycon.te );
			ShowMe();
			puts( "[Stop in level3]" );
			exit(1);
		}
	}

	p0 = AbunxIon - p2 - p1;

	/* population of lowest level */
	PopLevls.PopLevels[0] = (float)p0;
	if( p0 <= 0. )
	{
		fprintf( ioQQQ, " level3: insane n0 pop p1, 2, abun=%10.2e%10.2e%10.2e \n", 
		  p1, p2, AbunxIon );
		ShowMe();
		puts( "[Stop in level3]" );
		exit(1);
	}

	/* level populations for line opacities */
	t21->PopLo = (float)p1;
	t10->PopLo = (float)p0;
	t20->PopLo = (float)p0;
	t21->PopOpc = (float)(p1 - p2*g1/g2);
	t10->PopOpc = (float)(p0 - p1*g0/g1);
	t20->PopOpc = (float)(p0 - p2*g0/g2);
	t21->PopHi = (float)p2;
	t10->PopHi = (float)p1;
	t20->PopHi = (float)p2;

	/* line emission - net emission in line */
	t21->phots = (float)(t21->Aul * t21->Pesc*p2);
	t21->xIntensity = t21->phots * t21->EnergyErg;
	t21->ots = (float)(p2 * t21->Aul * t21->Pdest) ;

	t20->phots = (float)(t20->Aul * t20->Pesc*p2);
	t20->xIntensity = t20->phots * t20->EnergyErg;
	t20->ots = (float)(p2 * t20->Aul * t20->Pdest) ;

	t10->phots = (float)(t10->Aul * t10->Pesc*p1);
	t10->xIntensity = t10->phots * t10->EnergyErg;
	t10->ots = (float)(p2 * t10->Aul * t10->Pdest) ;

	/*  now add thess lines to ots field */
	AddOTSLin( t21->ots , t21->ipCont);
	AddOTSLin( t20->ots , t20->ipCont);
	AddOTSLin( t10->ots , t10->ipCont);

	/* total intensity used to divide line up - one may be 0 */
	/* >>chng 99 nov 30, rewrite algebra so double prec throughout,
	 * for very low density models */
	/*TotInten = t21->xIntensity + t20->xIntensity + t10->xIntensity;*/
	TotInten = (double)t21->phots * (double)t21->EnergyErg 
		+ (double)t20->phots * (double)t20->EnergyErg + 
		(double)t10->phots * (double)t10->EnergyErg;

	/* fraction that was collisionally excited */
	if( r12 > 0. )
	{
		t21->ColOvTot = (float)(c12/r12);
	}
	else
	{
		t21->ColOvTot = 0.;
	}

	if( r01 > 0. )
	{
		t10->ColOvTot = (float)(c01/r01);
	}
	else
	{
		t10->ColOvTot = 0.;
	}

	if( r02 > 0. )
	{
		t20->ColOvTot = (float)(c02/r02);
	}
	else
	{
		t20->ColOvTot = 0.;
	}

	/* heating or cooling due to each line */
	heat20 = p2*c20*ener20;
	cool02 = p0*c02*ener20;
	heat21 = p2*c21*ener21;
	cool12 = p1*c12*ener21;
	heat10 = p1*c10*ener10;
	cool01 = p0*c01*ener10;

	/* now get net heating or cooling */
	cnet02 = cool02 - heat20*t20->ColOvTot;
	hnet02 = heat20 *(1. - t20->ColOvTot);
	cnet12 = cool12 - heat21*t21->ColOvTot;
	hnet12 = heat21 *(1. - t21->ColOvTot);
	cnet01 = cool01 - heat10*t10->ColOvTot;
	hnet01 = heat10 *(1. - t10->ColOvTot);
	/* >>chng 96 nov 22, very dense cool models, roundoff error
	 *could cause [OI] 63 mic to be dominant heating, cooling, or
	 *just zero
	 * >>chng 96 dec 6, above change caused o iii 1666 cooling to
	 *be zeroed when important in a model in kirk's grid.  was at 1e-6,
	 * >>chng to 1e-7
	 * >>chng 96 dec 17, from 1e-7 to 1e-8, caused temp fail */
	/* >>chng 99 nov 29, had been 1e-30 to prevent div by zero (?),
	 * change to dble max since 1e-30 was very large compared to
	 * cooling for 1e-10 cm-3 den models */
	/*if( fabs(cnet01/MAX2(1e-30,cool01)) < 1e-8 )*/
	if( fabs(cnet01/MAX2(DBL_MIN,cool01)) < 1e-8 )
	{
		Lev3Fail.nLev3Fail = 4;
		cnet02 = 0.;
		hnet02 = 0.;
		cnet12 = 0.;
		hnet12 = 0.;
		cnet01 = 0.;
		hnet01 = 0.;
	}

	TotCool = cnet02 + cnet12 + cnet01;
	TotHeat = hnet02 + hnet12 + hnet01;

	if( TotInten > 0. )
	{
		cool02 = TotCool * t20->phots * (double)t20->EnergyErg/TotInten;
		cool12 = TotCool * t21->phots * (double)t21->EnergyErg/TotInten;
		cool01 = TotCool * t10->phots * (double)t10->EnergyErg/TotInten;
		heat20 = TotHeat * t20->phots * (double)t20->EnergyErg/TotInten;
		heat21 = TotHeat * t21->phots * (double)t21->EnergyErg/TotInten;
		heat10 = TotHeat * t10->phots * (double)t10->EnergyErg/TotInten;
		t20->cool = (float)cool02;
		t21->cool = (float)cool12;
		t10->cool = (float)cool01;
		t20->heat = (float)heat20;
		t21->heat = (float)heat21;
		t10->heat = (float)heat10;
	}
	else
	{
		Lev3Fail.nLev3Fail = 5;
		cool02 = 0.;
		cool12 = 0.;
		cool01 = 0.;
		heat20 = 0.;
		heat21 = 0.;
		heat10 = 0.;
		t20->cool = 0.;
		t21->cool = 0.;
		t10->cool = 0.;
		t20->heat = 0.;
		t21->heat = 0.;
		t10->heat = 0.;
	}

	/* fraction that produced lines */
	t21->AovTot = (float)(t21->Aul*t21->Pesc/r21);
	t10->AovTot = (float)(t10->Aul*t10->Pesc/r10);
	t20->AovTot = (float)(t20->Aul*t20->Pesc/r20);

	/* add cooling due to each line */
	/* >>chng 99 nov 30, rewrite algebra to keep precision on very low density models*/
	/*coladd(chLab, (long)t21->WLAng ,t21->cool);*/
	/*coladd(chLab, (long)t10->WLAng ,t10->cool);*/
	/*coladd(chLab, (long)t20->WLAng ,t20->cool);*/
	coladd(chLab, (long)t21->WLAng ,cool12);
	coladd(chLab, (long)t10->WLAng ,cool01);
	coladd(chLab, (long)t20->WLAng ,cool02);

	/* derivative of cooling function
	 * dC/dT = Cooling * ( T(excitation)/T_e^2 - 1/(2T) )
	 * in following I assume that a 1-2 exciation will have the 0-2
	 * expnential in the dcdt term = NOT A TYPO */

	tcool.dCooldT += t10->cool*(temp01*tcool.tsq1 - tcool.halfte) + 
	  (t20->cool + t21->cool)*(temp02*tcool.tsq1 - tcool.halfte);
	/* two t20->t's above are not a typo!
	 * */

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