/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*CoolOxyg evaluate total cooling due to oxygen */
#include "cddefines.h"
#include "coolheavy.h"
#include "dense.h"
#include "taulines.h"
#include "h2.h"
#include "phycon.h"
#include "embesq.h"
#include "hmi.h"
#include "oxy.h"
#include "colden.h"
#include "mole.h"
#include "ligbar.h"
#include "thermal.h"
#include "lines_service.h"
#include "atoms.h"
#include "cooling.h"
/* oi3pcs make collision strengths with electrons for 3P O I ground term.
 *>>refer o1	cs	Berrington, K.A. 1988, J.Phys.B, 21, 1083  for Te > 3000K
 *>>refer	o1	cs	Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K
 * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom)
 * written by Kirk Korista, 29 may 96
 * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data)
 * all data above 3000K have been adjusted down by a constant factor to make
 * them line up with Bell et al. data. */
static void oi3Pcs(double *cs21, 
	  double *cs20, 
	  double *cs10);

void CoolOxyg()
{
	double a21, 
	  a31, 
	  a32, 
	  a41, 
	  a42, 
	  a43, 
	  a51, 
	  a52, 
	  a53, 
	  a54, 
	  a6300, 
	  a6363, 
	  aeff, 
	  cs, 
	  cs2s2p, 
	  cs2s3p ,
	  cs01, 
	  cs02, 
	  cs12, 
	  cs13, 
	  cs21, 
	  cs23, 
	  cs31, 
	  cs32, 
	  cs41, 
	  cs42, 
	  cs43, 
	  cs51, 
	  cs52, 
	  cs53, 
	  cs54, 
	  csoi, 
	  Te_bounded, 
	  Te_bounded_log ,
	  o3cs23, 
	  p[5],
	  r12 , 
	  r13,
	  pump_rate;
	double cs_ab ,cs_ac ;
	/* these were added by Peter van Hoof to update the collision
	 * data within the OI ground term */
	double cse01,cse12,cse02 ,
	  csh01,cshe01,csh201,csh12,cshe12,csh212,csh02,cshe02,csh202 ,
	  csh2o01,csh2o02,csh2o12,csh2p01,csh2p02,csh2p12,csp01,csp02,
	  csp12;
	static double go2[5]={4.,6.,4.,4.,2.};
	static double exo2[4]={26808.4,21.0,13637.5,1.5};
	/* these will be used to update change in cs wrt temperature,
	 * and its affect on cooling deriv */
	static double oi_cs_tsave=-1. , oi_te_tsave=-1. , oi_dcdt_tsave=-1.;
	long int i;
	static int lgFirst=TRUE;
	static long int *ipO4Pump=NULL,
		nO4Pump=0;

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

	/* following does the OI Bowen Ly-bet pumping problem */
	atom_oi_calc(&CoolHeavy.coolOi);
	CoolAdd("o  1",8446,CoolHeavy.coolOi);

	/* O I 6300, 6363, A from
	 * >>refer	all	all	Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
	 * only electrons assumed; ok if ionized fraction GT 1E-4
	 * >>refer	oi	cs	Federman, S.R., & Shipsey, E.J. 1983, ApJ, 269, 791 */
	csoi = 2.68e-5*phycon.te*(1. + 1.67e-6*phycon.te - 2.95e-10*phycon.te*
	  phycon.te);
	csoi = MAX2(0.1,csoi);
	cs23 = 1.05e-3*phycon.sqrte;
	cs13 = 3.24e-6*phycon.te;
	a6300 = TauLines[ipT6300].Aul*TauLines[ipT6300].Pesc;
	TauLines[ipT6300].PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./5.);
	TauLines[ipT6300].PopLo = (dense.xIonDense[ipOXYGEN][0]*5./5.);
	TauLines[ipT6300].PopHi = 0.;
	TauLines[ipT6300].cs = (float)(csoi*5./9.);
	TauLines[ipT6363].PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./3.);
	TauLines[ipT6363].PopLo = (dense.xIonDense[ipOXYGEN][0]*5./3.);
	TauLines[ipT6363].PopHi = 0.;
	TauLines[ipT6363].cs = (float)(csoi*3./9.);
	a6363 = TauLines[ipT6363].Aul*TauLines[ipT6363].Pesc;

	a21 = a6300 + a6363 + oxy.d6300;
	a32 = TauLines[ipT5577].Aul*TauLines[ipT5577].Pesc;

	PutCS(cs23,&TauLines[ipT5577]);

	/* rate of new populations of O^0 due to dissociation of OH,
	 * co.rate_OH_dissoc is rate OH -> O + H [cm-3 s-1],
	 * must make it per unit O atom, so this rate is s-1 excitations per O atom */
	r12 = co.rate_OH_dissoc*0.55/SDIV( dense.xIonDense[ipOXYGEN][0] );
	r13 = co.rate_OH_dissoc*0.05/SDIV( dense.xIonDense[ipOXYGEN][0] );

	/* below is correction for fraction of excitations the produce emission */
	CoolHeavy.c6300_frac_emit = (a6300+a6363)/(a6300+a6363+csoi*dense.cdsqte/5.) ;
	CoolHeavy.c5577_frac_emit = (a32)/(a32+cs23*dense.cdsqte/3.);

	/* d6300 is the photoionization rate from the excited level
	 * was computed when ionization balance done */
	CoolHeavy.c5577 = atom_pop3(9.,5.,1.,csoi,cs13,cs23,a21,7.56e-2,a32,
	  22590.,25775.7,&oxy.poiexc,dense.xIonDense[ipOXYGEN][0],0.,r12,r13)*a32*
	  3.57e-12;

	TauLines[ipT5577].PopOpc = oxy.poiexc;
	TauLines[ipT5577].PopLo = oxy.poiexc;
	TauLines[ipT5577].PopHi = 0.;

	/* convert poiexc to relative abundances */
	/* >>chng 04 apr 20, include correction for fraction of 6300 due to OH pump */
	CoolHeavy.c5577 *= (1.-r13/(SDIV(atoms.c13)));

	CoolHeavy.c6300 = oxy.poiexc*a6300*TauLines[ipT6300].EnergyErg *
		(1.-r12/(SDIV(atoms.c12)));
	CoolHeavy.c6363 = oxy.poiexc*a6363*TauLines[ipT6363].EnergyErg *
		(1.-r12/(SDIV(atoms.c12)));
	/* must introduce correction for fraction of 6300 that is photo produced */
	thermal.dCooldT += CoolHeavy.c6300*(2.28e4*thermal.tsq1 + thermal.halfte) *
		/* note that atoms.c12 has all ra tes 1->2 */
		(1.-r12/(SDIV(atoms.c12)));

	oxy.poiexc /= (float)MAX2(1e-20,dense.xIonDense[ipOXYGEN][0]);
	CoolAdd("o  1",5577,CoolHeavy.c5577);
	CoolAdd("o  1",6300,CoolHeavy.c6300);
	CoolAdd("o  1",6363,CoolHeavy.c6363);

	/* O I fine structure lines rad data from 
	 * >>refer	all	cs	Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
	 * >>refer	o1	as	Berrington, K.A. 1988, J.Phys. B, 21, 1083
	 * hydrogen collisions from 
	 * >>refer	oi	cs	Tielens, A.G.G., & Hollenbach, D. 1985, ApJ, 291, 722
	 * factor in () is their rate coef
	 * assume H2 and H are same
	 * CDSQTE = 8.629E-6*EDEN/SQRTE
	 * cs01 = 9.8e-6*te + (4.2e-12*te70/te03) / cdsqte * 3. * hdcor
	 * cs12 = 2.6e-6*te  + (1.5e-10*sqrte/te03/te03)/cdsqte*hdcor
	 * cs02 = 2.9e-6*te  + (1.1e-12*te70*te10)/cdsqte*hdcor
	 * evaluate fits to OI electron rates, indices on var do not agree
	 * with Kirk's in sub, but are OK */
	/*==============================================================*/
	/* >>>chng 99 jun 01,
	 * following changed to be parallel to Peter van Hoof's changes 
	 * in the fortran C90.05 version 
	 * following is collisions with electrons */
	oi3Pcs(&cse01,&cse02,&cse12);
	/* remember the electron part of the cs */
	cs = cse01;

	/* rate coefficients for collisional de-excitation of O^o(3P) with H^o(2S1/2)
	 * NOTE: due to lack of data these relations are extrapolated to higher Te !
	 * >>refer	o1	cs	Launay & Roueff 1977, AA 56, 289
	 * the first fit is for Te <= 300K, the second for higher temps
	 * these data are valid for 50K <= Te <= 1000K*/
	csh12 = MAX2(2.00e-11*phycon.te30*phycon.te05*phycon.te02,
		7.67e-12*phycon.sqrte*phycon.te03);

	/* these data are valid for 100K <= Te <= 1000K */
	csh01 = MIN2(3.12e-12*phycon.te70*phycon.te02*phycon.te02,
		7.51e-12*phycon.sqrte*phycon.te05*phycon.te03) ;

	/* these data are valid for 150K <= Te <= 1000K*/
	csh02 = MIN2(6.96e-13*phycon.te/phycon.te10/phycon.te02,
		1.49e-12*phycon.te70*phycon.te05);

	/* rate coefficients for collisional de-excitation of O^o(3P) with He^o(1S)
	 * NOTE: due to lack of data these relations are extrapolated to higher Te !
	 * >>refer	oi	cs	Monteiro & Flower 1987, MNRAS 228, 101
	 * the first fit is for Te <= 300K, the second for higher temps
	 * these data are valid for 100K <= Te <= 1000K */
	cshe12 = MIN2(8.09e-16*phycon.te32*phycon.te10*phycon.te03,
		9.72e-15*phycon.te*phycon.te20);

	cshe01 = 1.57e-12*phycon.te70/phycon.te01;

	cshe02 = MIN2(1.80e-12*phycon.te70*phycon.te05,
		4.45e-12*phycon.te70/phycon.te10);

	if( phycon.te<=1.5e3 )
	{
		/* >>chng 04 mar 15, use explicit ortho-para densities */
		double ortho_frac = h2.ortho_density/SDIV(hmi.H2_total);
		/* rate coefficients for collisional de-excitation of O^o(3P) with H2(J=1,0)
		 * >>refer	oi	cs	Jaquet et al. 1992, J.Phys.B 25, 285
		 * these data are valid for 40K <= Te <= 1500K
		 * the first entry is contribution from ortho H2, the second para H2.*/
		csh2o12 = 2.21e-14*phycon.te*phycon.te10/phycon.te01;
		csh2p12 = 9.45e-15*phycon.te*phycon.te20/phycon.te01;
		csh212 = ortho_frac*csh2o12 + (1.-ortho_frac)*csh2p12;

		csh2o01 = 2.37e-11*phycon.te30*phycon.te10/phycon.te02;
		csh2p01 = 3.40e-11*phycon.te30*phycon.te02;
		csh201 = ortho_frac*csh2o01 + (1.-ortho_frac)*csh2p01;

		csh2o02 = 4.86e-11*phycon.te30*phycon.te02*phycon.te02;
		csh2p02 = 6.82e-11*phycon.te30/phycon.te02;
		csh202 = ortho_frac*csh2o02 + (1.-ortho_frac)*csh2p02 ;
	}
	else
	{
		csh212 = 0.;
		csh201 = 0.;
		csh202 = 0.;
	}

	/* effecitive collision strength of O^o(3P) with p
	 * >>refer	oi	cs	Pequignot, D. 1990, A&A 231, 499
	 * original data:
	 * >>refer	oi	cs	Chambaud et al., 1980, J.Phys.B, 13, 4205 (upto 5000K)
	 * >>refer	oi	cs	Roueff, private communication (10,000K and 20,000K)*/
	if( phycon.te<=1000. )
	{
		csp01 = MAX2(2.22e-5*phycon.te/phycon.te10,
			/* >>chng 05 jul 05, eden to dense.EdenHCorr */
			/*2.69e-6*phycon.te*phycon.te30)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
			2.69e-6*phycon.te*phycon.te30)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;

		csp02 = MIN2(7.07e-8*phycon.te32*phycon.te10,2.46e-7*
			/* >>chng 05 jul 05, eden to dense.EdenHCorr */
			/*phycon.te32/phycon.te10)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
			phycon.te32/phycon.te10)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;
	}
	else
	{
		csp01 = MIN2(2.69e-6*phycon.te*phycon.te30,9.21e-5*phycon.te/phycon.te10/
			/* >>chng 05 jul 05, eden to dense.EdenHCorr */
			/*phycon.te03)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
			phycon.te03)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;

		csp02 = MIN2(2.46e-7*phycon.te32/phycon.te10,5.21e-5*phycon.te/
			/* >>chng 05 jul 05, eden to dense.EdenHCorr */
			/*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
			phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;
	}

	csp12 = MIN2(2.35e-6*phycon.te*phycon.te05*phycon.te01,3.98e-5*
		/* >>chng 05 jul 05, eden to dense.EdenHCorr */
		/*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
		/*phycon.te70/phycon.te01)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
		phycon.te70/phycon.te01)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr;

	/*cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.Hmolec[ipMH2g])/
		dense.cdsqte;
	cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.Hmolec[ipMH2g])/
		dense.cdsqte;
	cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.Hmolec[ipMH2g])/
		dense.cdsqte;*/

	cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.H2_total)/
		dense.cdsqte;
	cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.H2_total)/
		dense.cdsqte;
	cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.H2_total)/
		dense.cdsqte;

	/* end changes */
	/*==============================================================*/
	PutCS(cs01,&TauLines[ipT63]);
	PutCS(cs12,&TauLines[ipT146]);
	PutCS(cs02,&TauDummy);
	atom_level3(&TauLines[ipT63],&TauLines[ipT146],&TauDummy);

	/* now save pops to add col den in radinc */
	for( i=0; i<3; ++i)
	{
		/* >>chng 02 oct 23, bug - had been O1Colden rather than O1Pops */
		colden.O1Pops[i] = (float)atoms.PopLevels[i];
	}
	/*fprintf(ioQQQ,"colllll\t%li\t%.3e\t%.3e\n",
		nzone, 
		colden.O1Pops[0], 
		colden.O1Pops[1] );*/

	/* >>>chng 99 apr 29, for orion pdr.in calc, temp oscillated due to
	 * bad dC/dT estimate, due to rapid change in cs with T.  added following*/
	/* when neutral collisions onto oi dominate the cooling the derivative
	 * produced by atom_level3 is vastly too small because of the strong temperature
	 * dependence of the neutral collision strenght.  Increase cooling deriv by
	 * ratio of cs due to electron (nearly constant) and total */
	/* following was correction for this between 00 apr 29 and 02 jul 25 */
#	if 0
	if( cs01 > SMALLFLOAT )
	{
		thermal.dCooldT += 
			(cs01-cs)/cs01*TauLines[ipT63].cool*
			TauLines[ipT63].EnergyK*thermal.tsq1;
	}
#	endif

	/* >>chng 02 jul 25, keep track of change in cs for 63 micron line */
	if( fabs(phycon.te - oi_te_tsave)/phycon.te > 1e-4 )
	{
		/* very first time we come through, previous values are -1 */
		if(oi_te_tsave>0. )
		{
			oi_dcdt_tsave = (cs01-oi_cs_tsave) / (phycon.te-oi_te_tsave);
		}
		else
		{
			oi_dcdt_tsave = 0.;
		}
		oi_cs_tsave = cs01;
		oi_te_tsave	= phycon.te;

		/* >>chng 03 jan 21, this factor could become very large - it should
		 * always be positive since neutral cs's are, and not much bigger than
		 * the usual derivative */
		/* can't be negative */
		oi_dcdt_tsave = MAX2( 0. , oi_dcdt_tsave);
		/* can't be bigger than several times normal dC/dT */
		oi_dcdt_tsave = MIN2( TauLines[ipT63].EnergyK*thermal.tsq1*4. , oi_dcdt_tsave);
	}
	/* this is only deriv due to change in collision strength, which is capped to be
	 * less than 4x the thermal deriv just above */
	thermal.dCooldT += TauLines[ipT63].cool*oi_dcdt_tsave;

	/* this is a bebug print statement for the numerical cs deriv */
	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf(ioQQQ,"DEBUG OI\t%.2f\tte\t%.5e\tcool\t%.5e\tcs\t%.4e\told\t%.4e\tnew\t%.4e\n",
				fnzone,
				phycon.te, 
				TauLines[ipT63].cool , 
				TauLines[ipT63].cs , 
				TauLines[ipT63].cool*TauLines[ipT63].EnergyK*thermal.tsq1, 
				TauLines[ipT63].cool*oi_dcdt_tsave );
		}
	}
	/* O II lines, CS data from 
	 * A;s from NIST
	 * full five level atom, ex in wavenumbers
	 * more recent cs 
	 * >>refer	o2	cs	McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753
	 * >>chng 97 feb 19, change in last sig fig as per Bob Rubin e-mail feb 11 97
	 * A from NIST compilation
	 * >>refer	o2	as	Wiese, W.L., Fuhr, J.R., Deters, T.M. 1996, J Phys Chem Ref Data,
	 * >>refercon Monograph 7 */
	a21 = 3.058e-5;
	a31 = 1.776e-4;
	a41 = 5.22e-2;
	a51 = 2.12e-2;
	a32 = 1.30e-7;
	a42 = 0.09907;
	a52 = 0.0519;
	a43 = 0.0534;
	a53 = 0.08672;
	a54 = 1.41e-10;
	/* >>chng 97 feb 19, mostly in last sig figs, as per Bob Rubin e-mail and
	 * >>referold	o2	cs	McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
	/* >>chng 01 may 18, update cs to
	 * >>refer	o2	cs	McLaughlin, B.M., & Bell, K.L. 1998, J Phys B 31, 4317 */
	/*cs21 = 0.82;
	cs31 = 0.42;
	cs41 = 0.261;
	cs51 = 0.113;
	cs32 = 1.24;
	cs42 = 0.875;
	cs52 = 0.334;
	cs43 = 0.822;
	cs53 = 0.213;
	cs54 = 0.157;*/
	/* >>chng 02 mar 13, go back to older values as per Seaton/Osterbrock correspondence
	 * this is a mix of the McLauthlin Bell 93 and 98 results, with the expected
	 * branching to levels within the term */
	cs_ab = 1.24;
	cs_ac = 0.374;
	cs21 = cs_ab * 0.6;
	cs31 = cs_ab * 0.4;
	/* >>chng 04 nov 01, statistical weights were reversed, caught by Kevin Blagrave */
	/*cs41 = cs_ac * 0.3333;
	cs51 = cs_ac * 0.6667;*/
	/*TODO	2	add proper temperature dependence when good coll stren calculations become available */
	cs41 = cs_ac * 0.6667;
	cs51 = cs_ac * 0.3333;
	cs32 = 1.24;
	cs42 = 0.875;
	cs52 = 0.334;
	cs43 = 0.822;
	cs53 = 0.213;
	cs54 = 0.157;
	/* subroutine atom_pop5(g,ex,cs12,cs13,cs14,cs15,
	 *  1 cs23,cs24,cs25,cs34,cs35,cs45, a21,a31,a41,a51, a32,a42,a52,
	 *  2 a43,a53,a54, p,abund) */
	atom_pop5(go2,exo2,cs21,cs31,cs41,cs51,cs32,cs42,cs52,cs43,cs53,cs54,
	  a21,a31,a41,a51,a32,a42,a52,a43,a53,a54,p,dense.xIonDense[ipOXYGEN][1]);
	CoolHeavy.O3730 = (float)(p[1]*a21*5.34e-12);
	CoolHeavy.O3726 = (float)(p[2]*a31*5.34e-12);
	CoolHeavy.O2471 = (float)((p[3]*a41 + p[4]*a51)*8.05e-12);
	CoolHeavy.O7323 = (float)((p[3]*a42 + p[4]*a52)*2.72e-12);
	CoolHeavy.O7332 = (float)((p[3]*a43 + p[4]*a53)*2.71e-12);
	CoolHeavy.c3727 = CoolHeavy.O3730 + CoolHeavy.O3726;
	CoolHeavy.c7325 = CoolHeavy.O7323 + CoolHeavy.O7332;
	thermal.dCooldT += CoolHeavy.c3727*(3.86e4*thermal.tsq1 - thermal.halfte);
	CoolAdd("O  2",3727,CoolHeavy.c3727);
	CoolAdd("O  2",7325,CoolHeavy.c7325);
	CoolAdd("O  2",2470,CoolHeavy.c7325*0.78);

	/* remember ratio of radiative to total decays, to use for estimating
	 * recombination contribution in lines_lv1_li_ne */
	if( (p[3] + p[4]) > SMALLFLOAT )
	{
		CoolHeavy.O2_A3_tot = (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) /
			( (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) +
			( p[3]*(cs41+cs42+cs43)/go2[3] + p[4]*(cs51+cs52+cs53)/go2[4]) * 
			dense.cdsqte );
	}
	else
	{
		CoolHeavy.O2_A3_tot = 0.;
	}

	if( (p[1] + p[2]) > SMALLFLOAT )
	{
		CoolHeavy.O2_A2_tot = (p[1]*a21 + p[2]*a31 ) /
			( (p[1]*a21 + p[2]*a31 ) +
			( p[1]*cs21/go2[1] + p[2]*cs31/go2[2]) * 
			dense.cdsqte );
	}
	else
	{
		CoolHeavy.O2_A2_tot = 0.;
	}

	/* O II 833.9A, CS 
	 * >>refer	o2	cs	McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
	/* >>chng 01 aug 10, turn this line back on - no reason given for turning it off. */
	cs = 0.355*phycon.te10*phycon.te10*phycon.te003*phycon.te001*
	  phycon.te001;
	PutCS(cs,&TauLines[ipT834]);
	atom_level2(&TauLines[ipT834]);

	/* O III 1666, crit den=2.6+10, A from
	 * >>refer	all	cs	Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
	 * c.s. 
	 * >>refer	o3	cs	Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
	PutCS(0.756,&TauLines[ipT1666]);
	PutCS(0.454,&TauLines[ipT1661]);
	PutCS(1.3,&TauDummy);
	atom_level3(&TauDummy,&TauLines[ipT1666],&TauLines[ipT1661]);
	/* if( nLev3Fail.gt.0 )write(qq,'('' nlev3fail='',i3)') nLev3Fail
	 *
	 * oiii 304 not in cooling function
	 * TODO	2	put all these in cooling */
	TauLines[ipT304].PopOpc = dense.xIonDense[ipOXYGEN][2];
	TauLines[ipT304].PopLo = dense.xIonDense[ipOXYGEN][2];
	TauLines[ipT304].PopHi = 0.;

	/* o iii 5007+4959, As 96 nist */
	/* >>chng 01 may 04, As updated to
	 * >>refer	o3	as	Storey, P.J., & Zeippen, C.J., 2000, MNRAS, 312, 813-816,
	 * changed by 10 percent! */
	aeff = 0.027205 + oxy.d5007r;
	/* following data used in routine that deduces oiii temp from spectrum
	 *
	 * >>refer	o3	cs	Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273
	 * IP paper Y(ki) differ significantly from those
	 * calculated by
	 * >>refer	o3	cs	Burke, V.M., Lennon, D.J., & Seaton, M.J. 1989, MNRAS, 236, 353
	 * especially for 1D2-1S0.
	 * BLS89 is adopted for 1D2-1S0 and LB94 for 3P2,1-1D2.
	 * NB!!  these cs's must be kept parallel with those in Peimbert analysis */
	if( phycon.te > 3981 && phycon.te <= 1e5 )
	{
		oxy.o3cs12 = (float)(3.0211144 - 101.57536/phycon.sqrte + 355.06905*
			/*really is base e log of tempeature */
		  log(phycon.te)/phycon.te);
		o3cs23 = 0.32412181 + 79.051672/phycon.sqrte - 4374.7816/
		  phycon.te;
	}
	else if( phycon.te < 3981. )
	{
		oxy.o3cs12 = 2.15f;
		o3cs23 = 0.4781;
	}
	else
	{
		oxy.o3cs12 = 2.6594f;
		o3cs23 = 0.5304;
	}
	oxy.o3cs13 = 
		(float)(MIN2(0.36,0.0784*phycon.te10*phycon.te03*phycon.te01*
	  phycon.te003));
	/* >>chng 01 may 04, derive a21 from aeff, had been 0.02431 
	a21 = 0.02431;*/
	a21 = aeff - oxy.d5007r;
	a31 = .215634;
	a32 = 1.71;
	oxy.o3ex23 = 32947.;
	oxy.o3br32 = (float)(a32/(a31 + a32));
	oxy.o3enro = (float)(4.56e-12/3.98e-12);
	/* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
	oxy.poiii3 = 
		(float)(atom_pop3(9.,5.,1.,oxy.o3cs12,oxy.o3cs13,o3cs23,
	  a21,a31,a32,28990.,oxy.o3ex23,&oxy.poiii2,dense.xIonDense[ipOXYGEN][2],
	  oxy.d5007r,0.,0.));
	CoolHeavy.c4363 = oxy.poiii3*a32*4.56e-12;
	CoolHeavy.c5007 = oxy.poiii2*a21*3.98e-12;
	oxy.d5007t = (float)(CoolHeavy.c5007*oxy.d5007r/aeff);
	thermal.dCooldT += CoolHeavy.c5007*(2.88e4*thermal.tsq1 - thermal.halfte);
	thermal.dCooldT += CoolHeavy.c4363*6.20e4*thermal.tsq1;
	CoolAdd("O  3",5007,CoolHeavy.c5007);
	CoolAdd("O  3",4363,CoolHeavy.c4363);
	CoolAdd("O  3",2331,CoolHeavy.c4363*0.236);
	/* >>chng 02 jun 27, prevent crash on Sun with gcc 3.1, PvH */
	/* if( xIonDense[ipOXYGEN][2] > 0. ) */
	if( MAX2(oxy.poiii2,oxy.poiii3) > 0.f && dense.xIonDense[ipOXYGEN][2]>SMALLFLOAT)
	{
		oxy.poiii3 /= dense.xIonDense[ipOXYGEN][2];
		oxy.poiii2 /= dense.xIonDense[ipOXYGEN][2];
	}

	/* oiii ir lines, col str from iron project, 
	 * >>refer	o3	cs	Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
	PutCS(0.5590,&TauLines[ipTO88]);
	PutCS(1.335,&TauLines[ipT52]);
	PutCS(0.2832,&TauDummy);
	atom_level3(&TauLines[ipTO88],&TauLines[ipT52],&TauDummy);

	/* O III 834A, CS 
	 * >>refer	o3	cs	Aggarwal, K.M., 1985 A&A 146, 149. */
	PutCS(5.,&TauLines[ipT835]);
	atom_level2(&TauLines[ipT835]);
	/* call DumpLine( t835 ); */

	/* these cs data only extend over a modest Temp range */
	Te_bounded = MAX2(phycon.te,4500.);
	Te_bounded = MIN2(Te_bounded,450000.);
	Te_bounded_log = log(Te_bounded);

	/* O IV 789A, CS from 
	 * >>refer	o4	cs	Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
	cs = -3.0102462 + 109.22973/Te_bounded_log - 18666.357/Te_bounded;
	PutCS(cs,&TauLines[ipT789]);
	atom_level2(&TauLines[ipT789]);

	/* O IV 26 mic, CS 
	 * >>refer	o4	cs	Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425
	 * A=
	 * >>refer	o4	as	Brage, T., Judge, P.G., & Brekke, P. 1996, ApJ. 464, 1030 */
	/* O IV 26 micron */
	cs = 9.2141614 - 5.1727759e-3*sqrt(Te_bounded) - 58.116447/Te_bounded_log;
	PutCS(cs,&TauLines[ipT26]);


	/* one time initialization if first call, and level 2 lines are on */
	if( lgFirst && nWindLine )
	{
		lgFirst = FALSE;
		nO4Pump = 0;
		for( i=0; i<nWindLine; ++i )
		{
			/* don't test on nelem==ipIRON since lines on physics, not C, scale */
			if( TauLine2[i].nelem ==8 && TauLine2[i].IonStg==4  )
			{
				++nO4Pump;
			}
		}
		if( nO4Pump<0 )
			TotalInsanity();
		else if( nO4Pump > 0 )
			/* create the space - can't malloc 0 bytes */
			if( (ipO4Pump = (long *)MALLOC((unsigned)(nO4Pump)*sizeof(long) ) )==NULL )
				BadMalloc();
		nO4Pump = 0;
		for( i=0; i<nWindLine; ++i )
		{
			/* don't test on nelem==ipIRON since lines on physics, not C, scale */
			if( TauLine2[i].nelem ==8 && TauLine2[i].IonStg==4  )
			{
#				define	PRT_DEBUG	FALSE
#				if	PRT_DEBUG
				DumpLine( &TauLine2[i] );
#				endif
				ipO4Pump[nO4Pump] = i;
				++nO4Pump;
			}
		}
	}
	else
		/* level 2 lines are not enabled */
		nO4Pump = 0;

	/* now sum pump rates */
	pump_rate = 0.;
	for( i=0; i<nO4Pump; ++i )
	{
		pump_rate += TauLine2[ipO4Pump[i]].pump;
#		if	PRT_DEBUG
		fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
			i,
			TauLine2[ipO4Pump[i]].WLAng , TauLine2[ipO4Pump[i]].pump );
#		endif
#		undef	PRT_DEBUG
	}

	/*atom_level2(&TauLines[ipT26]);*/
	/*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
	/* >>refer	o4	cs	Blum, R.D., & Pradhan, A.K., 1992, ApJS 80, 425
	 * >>refer	o4	cs	Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
	AtomSeqBoron(&TauLines[ipT26], 
	  &TauLines[ipO4_1400], 
	  &TauLines[ipO4_1397], 
	  &TauLines[ipO4_1407], 
	  &TauLines[ipO4_1405], 
	  &TauLines[ipO4_1401], 
	  0.1367 , 0.1560 , 1.1564 , 0.0457 , pump_rate,"O  4");

	/* O V 630, CS from 
	 * >>refer	o5	cs	Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
	 * >>refercon	At. Data Nucl. Data Tables, 33, 195 */
	cs = MIN2(4.0,1.317*phycon.te10/phycon.te03);
	PutCS(cs,&TauLines[ipT630]);
	atom_level2(&TauLines[ipT630]);

	/* O V 1218; coll data from 
	 * >>refer	o4	cs	Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
	 * >>refercon	At. Data Nucl. Data Tables, 33, 195 */
	if( phycon.te <= 3.16e4 )
	{
		cs = 3.224/(phycon.te10*phycon.te03*phycon.te03*phycon.te003);
	}
	else
	{
		cs = 10.549/(phycon.te10*phycon.te10*phycon.te10/phycon.te03);
	}
	/* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
	PutCS(cs,&TauLines[ipT1214]);
	/* c1214 = AtomSeqBeryllium( .87,1.05,3.32, t1214, .0216) * 1.64E-11
	 * AtomSeqBeryllium(CS23,CS24,CS34,tarray,A41)
	 * A's 
	 * >>refer	o5	cs	Fleming, J., Bell, K.L, Hibbert, A., Vaeck, N., Godefroid, M.R.
	 * >>refercon	1996, MNRAS, 279, 1289 */
	AtomSeqBeryllium(.87,0.602,2.86,&TauLines[ipT1214],.02198);
	embesq.em1218 = (float)(atoms.PopLevels[3]*0.0216*1.64e-11);

	/* O VI 1035 li seq
	 * generate collision strengths, then stuff them in
	 * >>refer	o6	vs	Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(8,&TauLines[ipT1032],&TauLines[ipT150],&cs2s2p,&cs2s3p);
	PutCS(cs2s2p,&TauLines[ipT1032]);
	PutCS(cs2s2p*0.5,&TauLines[ipT1037]);
	/* no data for the 2-3 transition */
	PutCS(1.0,&TauDummy);
	/* solve the 3 level atom */
	atom_level3(&TauLines[ipT1037],&TauDummy,&TauLines[ipT1032]);

	PutCS(cs2s3p,&TauLines[ipT150]);
	atom_level2(&TauLines[ipT150]);

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

/*"oi3pcs.h" make collision strengths with electrons for 3P O I ground term.
 *>>refer	o1	cs	Berrington, K.A. 1988, J.Phys.B, 21, 1083  for Te > 3000K
 *>>refer	o1	cs	Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K
 * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom)
 * written by Kirk Korista, 29 may 96
 * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data)
 * all data above 3000K have been adjusted down by a constant factor to make
 * them line up with Bell et al. data. */
static void oi3Pcs(double *cs21, 
	  double *cs20, 
	  double *cs10)
{
	double a, 
	  b, 
	  c, 
	  d;

#	ifdef DEBUG_FUN
	fputs( "<+>oi3Pcs()\n", debug_fp );
#	endif
	/* local variables */

	/* 3P2 - 3P1 */
	if( phycon.te <= 3.0e3 )
	{
		*cs21 = 1.49e-4*phycon.sqrte/phycon.te02/phycon.te02;
	}
	else if( phycon.te <= 1.0e4 )
	{
		a = -5.5634127e-04;
		b = 8.3458102e-08;
		c = 2.3068232e-04;
		*cs21 = 0.228f*(a + b*phycon.te32 + c*phycon.sqrte);
	}
	else
	{
		*cs21 = 0.228*MIN2(0.222,5.563e-06*phycon.te*phycon.te05*
		  phycon.te02);
	}

	/* 3P2 - 3P0 */
	if( phycon.te <= 3.0e3 )
	{
		*cs20 = 4.98e-5*phycon.sqrte;
	}
	else if( phycon.te <= 1.0e4 )
	{
		a = -3.7178028e-04;
		b = 2.0569267e-08;
		c = 1.1898539e-04;
		*cs20 = 0.288*(a + b*phycon.te32 + c*phycon.sqrte);
	}
	else
	{
		*cs20 = 0.288*MIN2(0.0589,1.015e-05*phycon.te/phycon.te10/
		  phycon.te02/phycon.te005);
	}

	/* 3P1 - 3P0 */
	if( phycon.te <= 3.0e3 )
	{
		*cs10 = 1.83e-9*phycon.te*phycon.te30*phycon.te05;
	}
	else if( phycon.te <= 1.0e4 )
	{
		a = -5.9364373e-04;
		b = 0.02946867;
		c = 10768.675;
		d = 3871.9826;
		*cs10 = 0.0269*(a + b*exp(-0.5*POW2((phycon.te-c)/d)));
	}
	else
	{
		*cs10 = 0.0269*MIN2(0.074,7.794e-08*phycon.te32/phycon.te10/
		  phycon.te01);
	}


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

