/*CoolCarb evaluate total cooling due to carbon */
#include "cddefines.h"
#include "coolheavy.h"
#include "embesq.h"
#include "h.h"
#include "tepowers.h"
#include "phycon.h"
#include "taulines.h"
#include "hmi.h"
#include "r1335.h"
#include "tcool.h"
#include "edsqte.h"
#include "poplevls.h"
#include "ligbar.h"
#include "coladd.h"
#include "level2.h"
#include "putcs.h"
#include "level3.h"
#include "beseq.h"
#include "pop3.h"
#include "coolmetals.h"
#include "he.h"
#include "ionfracs.h"

void CoolCarb()
{
	double SaveAbun, 
	  a21, 
	  a31, 
	  a32, 
	  cs, 
	  cs01, 
	  cs02, 
	  cs12, 
	  cs13, 
	  cs23, 
	  cs2s2p, 
	  cs2s3p ,
	  popup;

	/* added to impliment Peter van Hoof additions for new ground term
	 * atomic collision data */
	double cse01,cse12,cse02,csh01,csh12,csh02,csh201,csh212,csh202 ,
	  csh2p01,csh2p12,csh2p02,csh2o01,csh2o12,csh2o02;
	  float pciexc ;

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

	/* subroutine level3( t10,t21,t20)
	 *
	 * Carbon cooling
	 *
	 * C I 1656, collision strength from transition prob */
	/*PutCS(7.3,t1656);
	level2(t1656);*/
	PutCS(7.3, &TauLines[ipT1656] );
	level2(&TauLines[ipT1656]);

	/* C I fine structure lines data originally from
	 * >>refer Tielens, A.G.G., & Hollenbach, D. 1985, ApJ, 291, 722
	 * >>chng 99 Jun 01, to more recent ground term collision data
	 * by Peter van Hoof */
	/* ======================================================== */

	/* effective collision strength of C I(3P) with e
	 * >>refer Johnson, C.T., Burke, P.G., Kingston, A.E. 1987, JPhB, 20, 2553
	 * these data are valid for 7.5K <= Te <= 10,000K*/
	if( phycon.te<=3.0e3 )
	{
		/* the first fit is valid for 10K <= Te <= 300K, 
		 * the second 300K <= Te <= 3000K*/
		cse01 = MAX2(4.80E-06*phycon.te*TePowers.te20/TePowers.te03,
			8.24E-07*TePowers.te32/TePowers.te01);

		cse12 = MAX2(7.67E-05*phycon.te/TePowers.te10/TePowers.te03,
			1.47E-06*TePowers.te32*TePowers.te10/TePowers.te03);

		cse02 = MAX2(4.72E-05*TePowers.te70*TePowers.te03,
			3.05E-07*TePowers.te32*TePowers.te10);
	}
	else
	{
		/* the first fit is valid for 300K <= Te <= 3000K, 
		 * the second upto 10,000K */
		cse01 = MIN2(8.24E-07*TePowers.te32/TePowers.te01,
			0.0035*TePowers.sqrte*TePowers.te01);

		cse12 = MIN2(1.47E-06*TePowers.te32*TePowers.te10/TePowers.te03,
			0.0088*TePowers.sqrte*TePowers.te01*TePowers.te005);

		cse02 = MIN2(3.05E-07*TePowers.te32*TePowers.te10,
			0.00448*TePowers.sqrte/TePowers.te10*TePowers.te03*TePowers.te005);
	}

	if( phycon.te<=1.0e3 )
	{
		/* rate coefficients for collisional de-excitation of C I(3P) with neutral H(2S1/2)
		 * >>refer Launay & Roueff 1977, AA 56, 289
		 * the first fit is for Te <= 100K, the second for Te >= 250K
		 * these data are valid for 4K <= Te <= 1000K*/
		csh01 = MAX2(1.61e-10,5.66e-11*TePowers.te20);

		/* these data are valid for 7K <= Te <= 1000K*/
		csh12 = MAX2(1.93e-10*TePowers.te05*TePowers.te03,
			5.64e-11*TePowers.te30*TePowers.te02);

		/* these data are valid for 10K <= Te <= 1000K*/
		csh02 = MAX2(1.08e-10/TePowers.te03,
			1.67e-11*TePowers.te30*TePowers.te02*TePowers.te02);
	}
	else
	{
		csh01 = 0.;
		csh12 = 0.;
		csh02 = 0.;
	}

	if( phycon.te<=1.2e3 )
	{
		/* rate coefficients for collisional de-excitation of O I(3P) with H2(J=1,0)
		 * >>refer Schroeder et al. 1991, J.Phys.B 24, 2487
		 * these data are valid for 10K <= Te <= 1200K
		 * the first entry is contribution from ortho H2, the second para H2.*/
		if( phycon.te<=30. )
		{
			csh2p01 = MIN2(8.38E-11*TePowers.te05*TePowers.te01,
				2.12e-10/TePowers.te20/TePowers.te05/TePowers.te01);

			csh2o01 = MIN2(5.17E-11*TePowers.te10*TePowers.te05,
				1.07e-10/TePowers.te10*TePowers.te01);
		}
		else if( phycon.te<=150. )
		{
			csh2p01 = MAX2(6.60e-11,
				2.12e-10/TePowers.te20/TePowers.te05/TePowers.te01);

			csh2o01 = MAX2(7.10e-11,
				1.07e-10/TePowers.te10*TePowers.te01);
		}
		else
		{
			csh2p01 = MAX2(6.60e-11,3.38e-11*TePowers.te10*TePowers.te03);
			csh2p01 = MIN2(8.10e-11,csh2p01);

			csh2o01 = MAX2(7.1e-11,3.37e-11*TePowers.te10*TePowers.te02*TePowers.te02);
			csh2o01 = MIN2(8.57e-11,csh2o01);
		}

		/* assume standard mixture of ortho and para H2*/
		csh201 = 0.75*csh2o01 + 0.25*csh2p01;
		if( phycon.te<=30. )
		{
			csh2p12 = MIN2(1.48E-10*TePowers.te05*TePowers.te02,
				2.25e-10/TePowers.te03/TePowers.te03);
		}
		else if( phycon.te <= 100. )
		{
			csh2p12 = MAX2(1.75e-10,
				2.25e-10/TePowers.te03/TePowers.te03);
		}
		else
		{
			csh2p12 = MAX2(1.75e-10,6.23e-11*TePowers.te20*TePowers.te01);
			csh2p12 = MIN2(2.61e-10,csh2p12);
		}

		csh2o12 = MIN2(2.83e-10,4.46e-11*TePowers.te30/TePowers.te03);
		csh212 = 0.75*csh2o12 + 0.25*csh2p12;

		if( phycon.te<=30 )
		{
			csh2p02 = MIN2(8.67E-11*TePowers.te02*TePowers.te02,
				1.35e-10/TePowers.te10);
		}
		else if( phycon.te<=150. )
		{
			csh2p02 = MAX2(8.40e-11,
				1.35e-10/TePowers.te10);
		}
		else
		{
			csh2p02 = MAX2(8.4e-11,4.04e-11*TePowers.te10*TePowers.te02*TePowers.te02);
			csh2p02 = MIN2(1.04e-10,csh2p02);
		}

		csh2o02 = MIN2(1.11e-10,3.16e-11*TePowers.te20/TePowers.te02);
		csh202 = 0.75*csh2o02 + 0.25*csh2p02;
	}
	else
	{
		csh201 = 0.;
		csh212 = 0.;
		csh202 = 0.;
	}

	/* assume CS for HeI is the same as HI*/
	cs01 = cse01 + 3.*(csh01*(h.hi+he.hei) + csh201*hmi.htwo)/edsqte.cdsqte;
	cs12 = cse12 + 5.*(csh12*(h.hi+he.hei) + csh212*hmi.htwo)/edsqte.cdsqte;
	cs02 = cse02 + 5.*(csh02*(h.hi+he.hei) + csh202*hmi.htwo)/edsqte.cdsqte;

	PutCS( cs01 , &TauLines[ipT610] );
	PutCS( cs12 , &TauLines[ipT370] );
	PutCS( cs02 , &TauDummy );

	/* ======================================================== */
	/* end changes 99 Jun 01, by Peter van Hoot */
	level3( 
		&TauLines[ipT610],
		&TauLines[ipT370],
		&TauDummy);

	/* C I 9850, 8727, A from 
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143 */
	if( xIonFracs[1][5] > 0. && phycon.te < 40000. )
	{
		cs12 = 1.156e-4*phycon.te*(1.09 - 7.5e-6*phycon.te - 2.1e-10*
		  phycon.te*phycon.te);
		cs13 = 2.8e-3*TePowers.sqrte;
		cs23 = 2.764e-3*TePowers.sqrte;

		a21 = 3.26e-4*TauLines[ipT9830].Pesc;
		a31 = 2.73e-3;
		a32 = 0.528*TauLines[ipT8727].Pesc;
		/* 666 error! change to level3 */
		CoolHeavy.c8727 = pop3(9.,5.,1.,cs12,cs13,cs23,a21,a31,a32,
		  1.417e4,1.255e4,&pciexc,xIonFracs[1][5],0.)*a32*
		  2.28e-12;
		TauLines[ipT9830].PopOpc = xIonFracs[1][5];
		TauLines[ipT9830].PopLo = xIonFracs[1][5];
		TauLines[ipT9830].PopHi = 0.;
		TauLines[ipT9830].cs = (float)cs12;
		TauLines[ipT8727].PopOpc = (float)(CoolHeavy.c8727/(a32*2.28e-12));
		TauLines[ipT8727].PopLo = (float)(CoolHeavy.c8727/(a32*2.28e-12));
		TauLines[ipT8727].PopHi = 0.;
		TauLines[ipT8727].cs = (float)cs23;

		CoolHeavy.c9850 = pciexc*a21*2.02e-12;
		tcool.dCooldT += CoolHeavy.c9850*(1.468e4*tcool.tsq1 + tcool.halfte);
		tcool.dCooldT += CoolHeavy.c8727*(1.255e4*tcool.tsq1 + tcool.halfte);

		/* C I 9850 correction for deexcitation, needed for rec line */
		r1335.r9850 = (float)(a21/(a21 + cs12/5.*8.629e-6/TePowers.sqrte*phycon.eden));

	}
	else
	{
		CoolHeavy.c9850 = 0.;
		CoolHeavy.c8727 = 0.;
		r1335.r9850 = 0.;
		TauLines[ipT9830].PopOpc = 0.;
		TauLines[ipT9830].PopLo = 0.;
		TauLines[ipT9830].PopHi = 0.;
		TauLines[ipT8727].PopOpc = 0.;
		TauLines[ipT8727].PopLo = 0.;
		TauLines[ipT8727].PopHi = 0.;
	}
	coladd("C  1",8727,CoolHeavy.c8727);
	coladd("C  1",9850,CoolHeavy.c9850);

	/* C II 158 micron emission, A=
	 * >>refer Froese Fischer, C. 1983, J.Phys. B, 16, 157
	 * CS From 
	 * >>refer Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425
	 * neutral collision data from 
	 * >>refer Tielens, A.G.G., & Hollenbach, D. 1985, ApJ, 291, 722
	 * >>chng 96aug 01, better fit to cs  */
	cs = MIN2(2.20,0.403*TePowers.te20/TePowers.te02*TePowers.te001*
	  TePowers.te001) + 5.8e-10*TePowers.te02/edsqte.cdsqte*4.*(h.hi + 
	  hmi.htwo);
	PutCS(cs,&TauLines[ipT157]);
	level2(&TauLines[ipT157]);

	/* following used for pumping
	 * 666 error! put all these in cooling */
	TauLines[ipT386].PopOpc = xIonFracs[3][5];
	TauLines[ipT386].PopLo = xIonFracs[3][5];
	TauLines[ipT386].PopHi = 0.;
	TauLines[ipT310].PopOpc = xIonFracs[3][5];
	TauLines[ipT310].PopLo = xIonFracs[3][5];
	TauLines[ipT310].PopHi = 0.;
	TauLines[ipT291].PopOpc = xIonFracs[3][5];
	TauLines[ipT291].PopLo = xIonFracs[3][5];
	TauLines[ipT291].PopHi = 0.;
	TauLines[ipT280].PopOpc = xIonFracs[3][5];
	TauLines[ipT280].PopLo = xIonFracs[3][5];
	TauLines[ipT280].PopHi = 0.;
	TauLines[ipT274].PopOpc = xIonFracs[3][5];
	TauLines[ipT274].PopLo = xIonFracs[3][5];
	TauLines[ipT274].PopHi = 0.;
	TauLines[ipT270].PopOpc = xIonFracs[3][5];
	TauLines[ipT270].PopLo = xIonFracs[3][5];
	TauLines[ipT270].PopHi = 0.;

	/* all collision strengths and A'S from 
	 * >>refer Lennon, D.J., Dufton, P.L., Hibbert, A., Kingston, A.E. 1985,
	 * >>refer ApJ, 294, 200
	 * 3134=guess
	 * CS and trans prob from 
	 * >>refer Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 */
	PutCS(2.51,&TauLines[ipT2326]);
	cs = MIN2(6.73,2.316*TePowers.te10);
	PutCS(cs,&TauLines[ipT1335]);
	PutCS(5.39,&TauDummy);
	/* subroutine level3( t10,t21,t20) */
	level3(
		&TauLines[ipT2326],
		&TauDummy,
		&TauLines[ipT1335]);

	/* assumed to be optically thin
	 * 666 error! transfer this line */
	CoolHeavy.c3134 = 0.;

	/* C III  1909
	 * A for 1909 itself from 
	 * >>refer Kwong, V., Fang, Z., Gibbons, T.T., Parkinson, W.H., Smith, P.L.
	 * >>refer 1993, ApJ, 411, 431
	 * experimental value of 121 is larger than old NS 96, cs from
	 * >>refer Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
	 * >>refer At. Data Nucl. Data Tables, 33, 195
	 * BESEQ(CS23,CS24,CS34,tarray,A41) */
	cs = MIN2(1.1,2.67/TePowers.te10);
	PutCS(cs,&TauLines[ipT1909]);
	/* C1909 = BESEQ(.96,.73,2.8 , T1909 ,5.19E-3 )
	 * A's 
	 * >>refer Fleming, J., Bell, K.L, Hibbert, A., Vaeck, N., Godefroid, M.R.
	 * >>refer 1996, MNRAS, 279 , 1289 */
	beseq(.96,.73,2.8,&TauLines[ipT1909],5.149e-3);
	/* call dumpline( t1909 ) */
	embesq.em1908 = (float)(PopLevls.PopLevels[3]*5.19e-3*1.04e-11);

	/* CIII 1175 excited state line 
	 * following were computed by previous call to beseq  */
	popup = PopLevls.PopLevels[1] + PopLevls.PopLevels[2] + PopLevls.PopLevels[3];
	SaveAbun = xIonFracs[3][5];
	xIonFracs[3][5] = (float)popup;
	/* cs 
	 * >>refer Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
	 * >>refer At. Data Nucl. Data Tables, 33, 195 */
	cs = MIN2(30.,4.806*TePowers.te10*TePowers.te05/TePowers.te01/
	  TePowers.te003);
	PutCS(18.45,&TauLines[ipc31175]);
	level2(&TauLines[ipc31175]);
	xIonFracs[3][5] = (float)SaveAbun;

	/* C III 977, cs from 
	 * >>refer Berrington, K.A. 1985, J.Phys. B, 18, L395 */
	cs = MIN2(7.0,1.556*TePowers.te10);
	PutCS(cs,&TauLines[ipT977]);
	level2(&TauLines[ipT977]);

	/* CIV 1548, 1550 doublet
	 * >>refer Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(
		6,
		&TauLines[ipT1548],
		&TauLines[ipT312],
		&cs2s2p,&cs2s3p);
	PutCS(cs2s2p,&TauLines[ipT1548]);
	PutCS(cs2s2p*0.5,&TauLines[ipT1550]);
	PutCS(1.0,&TauDummy);
	level3(
		&TauLines[ipT1550],
		&TauDummy,
		&TauLines[ipT1548]);

	PutCS(cs2s3p,&TauLines[ipT312]);
	level2(&TauLines[ipT312]);


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