/*CoolOxyg evaluate total cooling due to oxygen */
#include "cddefines.h"
#include "oi3pcs.h"
#include "o3exc.h"
#include "coolheavy.h"
#include "oii.h"
#include "taulines.h"
#include "tepowers.h"
#include "phycon.h"
#include "embesq.h"
#include "logte.h"
#include "tcool.h"
#include "h.h"
#include "edsqte.h"
#include "hmi.h"
#include "o3data.h"
#include "poplevls.h"
#include "ionfracs.h"
#include "fivel.h"
#include "ligbar.h"
#include "coladd.h"
#include "level2.h"
#include "putcs.h"
#include "beseq.h"
#include "level3.h"
#include "pop3.h"
#include "coolmetals.h"
#include "he.h"
#include "p8446.h"

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, 
	  hold, 
	  o3cs23, 
	  p[5];
	/* 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};

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

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

	/* O I 6300, 6363, A from
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * only electrons assumed; ok if ionized fraction GT 1E-4
	 * >>refer 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*TePowers.sqrte;
	cs13 = 3.24e-6*phycon.te;
	a6300 = TauLines[ipT6300].Aul*TauLines[ipT6300].Pesc;
	TauLines[ipT6300].PopOpc = (float)(xIonFracs[1][7]*5./5.);
	TauLines[ipT6300].PopLo = (float)(xIonFracs[1][7]*5./5.);
	TauLines[ipT6300].PopHi = 0.;
	TauLines[ipT6300].cs = (float)(csoi*5./9.);
	TauLines[ipT6363].PopOpc = (float)(xIonFracs[1][7]*5./3.);
	TauLines[ipT6363].PopLo = (float)(xIonFracs[1][7]*5./3.);
	TauLines[ipT6363].PopHi = 0.;
	TauLines[ipT6363].cs = (float)(csoi*3./9.);
	a6363 = TauLines[ipT6363].Aul*TauLines[ipT6363].Pesc;

	/* d6300 is the photoionization rate from the excited level
	 * was computed when ionization balance done */
	a21 = a6300 + a6363 + o3exc.d6300;
	a32 = TauLines[ipT5577].Aul*TauLines[ipT5577].Pesc;
	PutCS(cs23,&TauLines[ipT5577]);
	CoolHeavy.c5577 = pop3(9.,5.,1.,csoi,cs13,cs23,a21,7.56e-2,a32,
	  22590.,25775.7,&o3exc.poiexc,xIonFracs[1][7],0.)*a32*
	  3.57e-12;
	TauLines[ipT5577].PopOpc = o3exc.poiexc;
	TauLines[ipT5577].PopLo = o3exc.poiexc;
	TauLines[ipT5577].PopHi = 0.;
	/* convert poiexc to relative abundances */
	CoolHeavy.c6300 = o3exc.poiexc*a6300*TauLines[ipT6300].EnergyErg;
	CoolHeavy.c6363 = o3exc.poiexc*a6363*TauLines[ipT6363].EnergyErg;
	tcool.dCooldT += CoolHeavy.c6300*(2.28e4*tcool.tsq1 + tcool.halfte);
	o3exc.poiexc /= (float)MAX2(1e-20,xIonFracs[1][7]);
	coladd("o  1",5577,CoolHeavy.c5577);
	coladd("o  1",6300,CoolHeavy.c6300);
	coladd("o  1",6363,CoolHeavy.c6363);

	/* O I fine structure lines rad data from 
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * >>refer Berrington, K.A. 1988, J.Phys. B, 21, 1083
	 * hydrogen collisions from 
	 * >>refer 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 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*TePowers.te30*TePowers.te05*TePowers.te02,
		7.67e-12*TePowers.sqrte*TePowers.te03);

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

	/* these data are valid for 150K <= Te <= 1000K*/
	csh02 = MIN2(6.96e-13*phycon.te/TePowers.te10/TePowers.te02,
		1.49e-12*TePowers.te70*TePowers.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 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*TePowers.te32*TePowers.te10*TePowers.te03,
		9.72e-15*phycon.te*TePowers.te20);

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

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

	if( phycon.te<=1.5e3 )
	{
		/* rate coefficients for collisional de-excitation of O^o(3P) with H2(J=1,0)
		 * >>refer 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*TePowers.te10/TePowers.te01;
		csh2p12 = 9.45e-15*phycon.te*TePowers.te20/TePowers.te01;

		/* assume standard mixture of ortho and para H2*/
		csh212 = 0.75*csh2o12 + 0.25*csh2p12;
		csh2o01 = 2.37e-11*TePowers.te30*TePowers.te10/TePowers.te02;
		csh2p01 = 3.40e-11*TePowers.te30*TePowers.te02;
		csh201 = 0.75*csh2o01 + 0.25*csh2p01;
		csh2o02 = 4.86e-11*TePowers.te30*TePowers.te02*TePowers.te02;
		csh2p02 = 6.82e-11*TePowers.te30/TePowers.te02;
		csh202 = 0.75*csh2o02 + 0.25*csh2p02 ;
	}
	else
	{
		csh212 = 0.;
		csh201 = 0.;
		csh202 = 0.;
	}

	/* effecitive collision strength of O^o(3P) with p
	 * >>refer Pequignot, 1990, A&A 231, 499
	 * original data: Chambaud et al., 1980, J.Phys.B, 13, 4205 (upto 5000K)
	 * Roueff, private communication (10,000K and 20,000K)*/
	if( phycon.te<=1000. )
	{
		csp01 = MAX2(2.22e-5*phycon.te/TePowers.te10,
			2.69e-6*phycon.te*TePowers.te30)*h.hii/phycon.eden;

		csp02 = MIN2(7.07e-8*TePowers.te32*TePowers.te10,2.46e-7*
			TePowers.te32/TePowers.te10)*h.hii/phycon.eden;
	}
	else
	{
		csp01 = MIN2(2.69e-6*phycon.te*TePowers.te30,9.21e-5*phycon.te/TePowers.te10/
			TePowers.te03)*h.hii/phycon.eden;

		csp02 = MIN2(2.46e-7*TePowers.te32/TePowers.te10,5.21e-5*phycon.te/
			TePowers.te20)*h.hii/phycon.eden;
	}

	csp12 = MIN2(2.35e-6*phycon.te*TePowers.te05*TePowers.te01,3.98e-5*
		TePowers.te70/TePowers.te01)*h.hii/phycon.eden;

	cs01 = cse01+csp01+3.*(csh01*h.hi + cshe01*he.hei + csh201*hmi.htwo)/
		edsqte.cdsqte;
	cs12 = cse12+csp12+(csh12*h.hi + cshe12*he.hei + csh212*hmi.htwo)/
		edsqte.cdsqte;
	cs02 = cse02+csp02+(csh02*h.hi + cshe02*he.hei + csh202*hmi.htwo)/
		edsqte.cdsqte;

#	if 0
	/* assume H atoms and molecule have same rate */
	hdcor = h.hi + hmi.htwo;
	/* add on hydrogen atom collision */
	/* remember the electron part of the cs */
	cs = cs01;
	cs01 = cs + (4.2e-12*TePowers.te70/TePowers.te03)/edsqte.cdsqte*3.*
	  hdcor;
	cs12 += (1.5e-10*TePowers.sqrte/TePowers.te03/TePowers.te03)/edsqte.cdsqte*
	  hdcor;
	cs02 += (1.1e-12*TePowers.te70*TePowers.te10)/edsqte.cdsqte*hdcor;
#	endif
	/* end changes */
	/*==============================================================*/
	PutCS(cs01,&TauLines[ipT63]);
	PutCS(cs12,&TauLines[ipT146]);
	PutCS(cs02,&TauDummy);
	level3(&TauLines[ipT63],&TauLines[ipT146],&TauDummy);

	/* >>>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 level3 is vastly too small.  Increase cooling deriv by
	 * ratio of cs due to electron (nearly constant) and total */
	if( cs01 > SMALLFLOAT )
	{
		tcool.dCooldT += 
			(cs01-cs)/cs01*TauLines[ipT63].cool*
			TauLines[ipT63].EnergyK*tcool.tsq1;
	}

	/* O II lines, CS data from 
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * A;s from NIST
	 * full five level atom, ex in wavenumbers
	 * more recent cs 
	 * >>refer 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 Wiese, W.L., Fuhr, J.R., Deters, T.M. 1996, J Phys Chem Ref Data,
	 * >>refer 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
	 * >>refer McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
	cs21 = 0.825;
	cs31 = 0.55;
	cs41 = 0.2765;
	cs51 = 0.1382;
	cs32 = 1.17;
	/* scale factor from
	 * >>refer McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
	cs42 = 0.730*1.207;
	cs52 = 0.295*1.207;
	cs43 = 0.408*1.207;
	cs53 = 0.275*1.207;
	cs54 = 0.287;
	/* subroutine fivel(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) */
	fivel(go2,exo2,cs21,cs31,cs41,cs51,cs32,cs42,cs52,cs43,cs53,cs54,
	  a21,a31,a41,a51,a32,a42,a52,a43,a53,a54,p,xIonFracs[2][7]);
	oii.O3730 = (float)(p[1]*a21*5.34e-12);
	oii.O3727 = (float)(p[2]*a31*5.34e-12);
	oii.O2471 = (float)((p[3]*a41 + p[4]*a51)*8.05e-12);
	oii.O7323 = (float)((p[3]*a42 + p[4]*a52)*2.72e-12);
	oii.O7332 = (float)((p[3]*a43 + p[4]*a53)*2.71e-12);
	CoolHeavy.c3727 = oii.O3730 + oii.O3727;
	CoolHeavy.c7325 = oii.O7323 + oii.O7332;
	tcool.dCooldT += CoolHeavy.c3727*(3.86e4*tcool.tsq1 - tcool.halfte);
	coladd("O  2",3727,CoolHeavy.c3727);
	coladd("O  2",7325,CoolHeavy.c7325);
	coladd("O  2",2470,CoolHeavy.c7325*0.78);

	/* O II 833.9A, CS 
	 * >>refer McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
	cs = 0.355*TePowers.te10*TePowers.te10*TePowers.te003*TePowers.te001*
	  TePowers.te001;
	PutCS(cs,&TauLines[ipT834]);
	level2(&TauLines[ipT834]);

	/* O III 1666, crit den=2.6+10, A from
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * c.s. 
	 * >>refer 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);
	level3(&TauDummy,&TauLines[ipT1666],&TauLines[ipT1661]);
	/* if( nLev3Fail.gt.0 )write(qq,'('' nlev3fail='',i3)') nLev3Fail
	 *
	 * oiii 304 not in cooling function
	 * 666 error! put all these in cooling */
	TauLines[ipT304].PopOpc = xIonFracs[3][7];
	TauLines[ipT304].PopLo = xIonFracs[3][7];
	TauLines[ipT304].PopHi = 0.;

	/* o iii 5007+4959, As 96 nist */
	aeff = 0.0263 + o3exc.d5007r;
	/* following data used in routine that deduces oiii temp from spectrum
	 *
	 * >>refer Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273
	 * IP paper Y(ki) differ significantly from those
	 * calculated by
	 * >>refer 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 )
	{
		o3data.o3cs12 = (float)(3.0211144 - 101.57536/TePowers.sqrte + 355.06905*
		  logte.alogete/phycon.te);
		o3cs23 = 0.32412181 + 79.051672/TePowers.sqrte - 4374.7816/
		  phycon.te;
	}
	else if( phycon.te < 3981. )
	{
		o3data.o3cs12 = 2.15f;
		o3cs23 = 0.4781;
	}
	else
	{
		o3data.o3cs12 = 2.6594f;
		o3cs23 = 0.5304;
	}
	o3data.o3cs13 = 
		(float)(MIN2(0.36,0.0784*TePowers.te10*TePowers.te03*TePowers.te01*
	  TePowers.te003));
	a21 = 0.02431;
	a31 = .215634;
	a32 = 1.71;
	o3data.o3ex23 = 32947.;
	o3data.o3br32 = (float)(a32/(a31 + a32));
	o3data.o3enro = (float)(4.56e-12/3.98e-12);
	/* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
	o3exc.poiii3 = 
		(float)(pop3(9.,5.,1.,o3data.o3cs12,o3data.o3cs13,o3cs23,
	  a21,a31,a32,28990.,o3data.o3ex23,&o3exc.poiii2,xIonFracs[3][7],
	  o3exc.d5007r));
	CoolHeavy.c4363 = o3exc.poiii3*a32*4.56e-12;
	CoolHeavy.c5007 = o3exc.poiii2*a21*3.98e-12;
	o3exc.d5007t = (float)(CoolHeavy.c5007*o3exc.d5007r/aeff);
	tcool.dCooldT += CoolHeavy.c5007*(2.88e4*tcool.tsq1 - tcool.halfte);
	tcool.dCooldT += CoolHeavy.c4363*6.20e4*tcool.tsq1;
	coladd("O  3",5007,CoolHeavy.c5007);
	coladd("O  3",4363,CoolHeavy.c4363);
	coladd("O  3",2331,CoolHeavy.c4363*0.236);
	if( xIonFracs[3][7] > 0. )
	{
		o3exc.poiii3 /= xIonFracs[3][7];
		o3exc.poiii2 /= xIonFracs[3][7];
	}

	/* oiii ir lines, col str from iron project, 
	 * >>refer 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);
	level3(&TauLines[ipTO88],&TauLines[ipT52],&TauDummy);

	/* O III 834A, CS 
	 * >>refer Aggarwal, K.M., 1985 A&A 146, 149. */
	PutCS(5.,&TauLines[ipT835]);
	level2(&TauLines[ipT835]);
	/* call DumpLine( t835 )
	 *
	 * O IV 26 mic, CS 
	 * >>refer Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425
	 * A=
	 * >>refer Brage, T., Judge, P.G., & Brekke, P. 1996, ApJ. 464, 1030 */
	hold = MAX2(phycon.te,4500.);
	hold = MIN2(hold,450000.);
	aeff = log(hold);
	cs = 9.2141614 - 5.1727759e-3*sqrt(hold) - 58.116447/aeff;
	PutCS(cs,&TauLines[ipT26]);
	level2(&TauLines[ipT26]);

	/* O IV 1402, wl from 
	 * >>refer Wills, D., & Netzer, H. 1979, ApJ, 233, 1
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * CS 
	 * >>refer Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
	cs = -13.860726 - 0.14091452*POW2(aeff) + 2.9626895*aeff;
	PutCS(cs,&TauLines[ipT1402]);
	level2(&TauLines[ipT1402]);

	/* O IV 789A, CS from 
	 * >>refer Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
	cs = -3.0102462 + 109.22973/aeff - 18666.357/hold;
	PutCS(cs,&TauLines[ipT789]);
	level2(&TauLines[ipT789]);
	/* call DumpLine( t789 )
	 *
	 * O V 630, CS from 
	 * >>refer Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
	 * >>refer At. Data Nucl. Data Tables, 33, 195 */
	cs = MIN2(4.0,1.317*TePowers.te10/TePowers.te03);
	PutCS(cs,&TauLines[ipT630]);
	level2(&TauLines[ipT630]);

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

	/* O VI 1035 li seq
	 * generate collision strengths, then stuff them in
	 * >>refer 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 */
	level3(&TauLines[ipT1037],&TauDummy,&TauLines[ipT1032]);

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


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