/home66/gary/public_html/cloudy/c08_branch/source/cool_oxyg.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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);

#if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
#pragma optimization_level 1
#endif
void CoolOxyg(void)
{
        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, 
          Te_bounded, 
          Te_bounded_log ,
          o3cs23, 
          p[5],
          r12 , 
          r13,
          pump_rate;
        /*double cs_ab ,cs_ac;(Commented out- Humeshkar Nemala)*/
        /* these were added by Peter van Hoof to update the collision
         * data within the OI ground term */
        double cse01=-1.,cse12=-1.,cse02 =-1.,
          csh01=-1.,cshe01=-1.,csh201=-1.,csh12=-1.,cshe12=-1.,csh212=-1.,csh02=-1.,cshe02=-1.,csh202 =-1.,
          csh2o01=-1.,csh2o02=-1.,csh2o12=-1.,csh2p01=-1.,csh2p02=-1.,csh2p12=-1.,csp01=-1.,csp02=-1.,
          csp12=-1.;
        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 derivative */
        static double oi_cs_tsave=-1. , oi_te_tsave=-1. , oi_dcdt_tsave=-1.;
        long int i;
        static bool lgFirst=true;
        static long int *ipO4Pump=NULL,
                nO4Pump=0;
        double rate_OH_dissoc;

        DEBUG_ENTRY( "CoolOxyg()" );

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

        /* [OI] 6300, 6363, 5575, etc 
         * >>chng 06 oct 02, Humeshkar Nemala incorporate Barklem cs data for OI
         * largest difference is 3P - 1D (6300+6363) which is now roughly 3x smaller */
        /* This is the transition from 3P(J=2,1,0;the levels are reversed) to 1D2
         * The rate coefficients were converted to CS
         *>>refer       oi      cs      Barklem,P.S.,2006,A&A (astroph 0609684) 
         * Data pts are avilable over 1000,3000,5000,8000,12000,20000 and 50000
         * Fits between 1000K and 20000K are reliable;this is not the case 
         * between 20000K & 50000K*/
        /* this was old calculation of cs12 that was replaced by Humeshkar */
#if 0
        cs12 = 2.68e-5*phycon.te*(1. + 1.67e-6*phycon.te - 2.95e-10*phycon.te*
          phycon.te);
        cs12 = MAX2(0.1,cs12);
#endif

        /* [OI] electron 6300, 6363 collision strength
         * >>chng 06 oct 02, Humeshkar Nemala code, based on Barklem paper */
        if(phycon.te < 8E3)
        {
                cs12 = (4E-08)*(phycon.te*phycon.te70*phycon.te05);
        }
        else if(phycon.te < 2E4)
        {
                cs12 = (4.630155E-05)*((phycon.te/phycon.te04)*phycon.te005*phycon.te0001);
        }
        else
        {
                /* this branch T > 20,000K */
                cs12 = (1.5394E-03)*(phycon.sqrte*phycon.te10*phycon.te01*phycon.te001*phycon.te0003);
        }

        /* this block adds on collisional excitation by H0 */
        {
                /* >>chng 06 aug 18, add atomic hydrogen collisional processes using rates from
                 * >>refer      O1      coll    Krems, R.V., Jameson, M.J. & Dalgarno, A. 2006, ApJ, 647, 1531 
                 * their equation 10 - the deecxiation rate coefficient, cm3 s-1
                 * >>referold   oi      cs      Federman, S.R., & Shipsey, E.J. 1983, ApJ, 269, 791 */
                double te_scale = phycon.te / 6000.;
                double rate_H0 = (1.74*te_scale + 0.6)*1e-12*sexp(0.47*te_scale) / sqrt(te_scale ) * 
                        dense.xIonDense[ipHYDROGEN][0];
                /*fprintf(ioQQQ,"DEBUG OI H0 rate %.2e " , cs12 );*/
                cs12 += ConvRate2CS( 5.f , (realnum)rate_H0 ); 
                /*fprintf(ioQQQ,"%.2e\n" , cs12 );*/
        }
        /* following was old fit to 2-3 5755 transition */
        /*cs23 = 1.05e-3*phycon.sqrte;*/
        /*>>chng 06 oct 02, Humeshkar Nemala's new fit to Barklem paper */
        if(phycon.te < 5E3)
        {
                cs23 = (7E-08)*(phycon.te*phycon.sqrte*phycon.te10*phycon.te007*phycon.te0001);
        }
        else if(phycon.te<2E4)
        {
                cs23 = (1.98479e-04)*((phycon.te70/phycon.te03)*phycon.te003*phycon.te0007);
        }
        else
        {
                cs23 = (9.31E-04)*(phycon.sqrte*phycon.te01*phycon.te007*phycon.te0005*phycon.te0001);
        }

        /* following was old fit to 1,2,3 -> 5 transition */
        /*cs13 = 3.24e-6*phycon.te;*/
        /* following is Humeshkar's new fit to Barklem paper */
        if(phycon.te < 2E4)
        {
                cs13 = (2E-05)*(phycon.sqrte*phycon.te30*phycon.te05*phycon.te01*(phycon.te004/phycon.te0002));
        }
        else
        {
                cs13 = (1.054E-03)*(phycon.sqrte/phycon.te04)*phycon.te003*phycon.te0005;
        }

        /* 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 */
        a6300 = TauLines[ipT6300].Emis->Aul*TauLines[ipT6300].Emis->Pesc;
        TauLines[ipT6300].Emis->PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./5.);
        TauLines[ipT6300].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./5.);
        TauLines[ipT6300].Hi->Pop = 0.;
        TauLines[ipT6300].Coll.cs = (realnum)(cs12*5./9.);
        TauLines[ipT6363].Emis->PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./3.);
        TauLines[ipT6363].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./3.);
        TauLines[ipT6363].Hi->Pop = 0.;
        TauLines[ipT6363].Coll.cs = (realnum)(cs12*3./9.);
        a6363 = TauLines[ipT6363].Emis->Aul*TauLines[ipT6363].Emis->Pesc;

        a21 = a6300 + a6363 + oxy.d6300;
        a32 = TauLines[ipT5577].Emis->Aul*TauLines[ipT5577].Emis->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 */
        rate_OH_dissoc = CO_findrate("PHOTON,OH=>O,H");
        r12 = rate_OH_dissoc*0.55/SDIV( dense.xIonDense[ipOXYGEN][0] );
        r13 = 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+cs12*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.,cs12,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].Emis->PopOpc = oxy.poiexc;
        TauLines[ipT5577].Lo->Pop = oxy.poiexc;
        TauLines[ipT5577].Hi->Pop = 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 /= (realnum)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.;
        }

        /* effective 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] = (realnum)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 strength.  Increase cooling derivative 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 derivative due to change in collision strength, which is capped to be
         * less than 4x the thermal derivative just above */
        thermal.dCooldT += TauLines[ipT63].Coll.cool*oi_dcdt_tsave;
        /* this is a bebug print statement for the numerical cs derivative */
        {
                /*@-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].Coll.cool , 
                                TauLines[ipT63].Coll.cs , 
                                TauLines[ipT63].Coll.cool*TauLines[ipT63].EnergyK*thermal.tsq1, 
                                TauLines[ipT63].Coll.cool*oi_dcdt_tsave );
                }
        }
        /* [O II]
         * two sets of A's, Zeippen 1982 being the default
         * >>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 
         * these are selected with the set atomic data ion oxygen 2 command */
        if( dense.lgAsChoose[ipOXYGEN][1] )
        {
                /* not used unless requested with set atomic data command */
                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;
        }
        else
        {
                /* default - the Zeippen (1982) values 
                *>>refer        O2      As      Zeippen, C. J. 1982, MNRAS, 198, 111
                 * these are recommended by 
                 *>>refer       O2      As      Wang et al. 2004 A&A, 427, 873
                 >>refer        O2      As      Chen, Qing & Li Phys Rev A, 2007, 76, 042507
                 */
                 a21 = 3.82e-5;
                 a31 = 1.65e-4;
                 a41 = 5.64e-2;
                 a51 = 2.3202e-2;
                 a32 = 1.20e-7;
                 a42 = 0.11678;
                 a52 = 0.0615;
                 a43 = 0.0614;
                 a53 = 0.10175;
                 a54 = 2.08e-11;
        }
        /* >>chng 01 may 18, update cs to
         * >>referold   o2      cs      McLaughlin, B.M., & Bell, K.L. 1998, J Phys B 31, 4317 
         * >>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 
         * >>chng 04 nov 01, statistical weights were reversed, caught by Kevin Blagrave 
         * >>refer      o2      cs      A.K.Pradhan, M.Montenegro,S.N.Nahar & W.Eissner,2006,MNRAS,366,L6-L9
         */
        cs21 = (realnum)(0.8229*(phycon.te007*phycon.te0005*phycon.te0001));
        cs31 = (realnum)(0.5976/phycon.te002);
        cs32 = (realnum)(2.6084/(phycon.te07/(phycon.te003*phycon.te0001)));
        cs41 = (realnum)(0.2471*phycon.te02*(phycon.te007/phycon.te0004));
        cs42 = (realnum)(0.7075*phycon.te03*phycon.te004*phycon.te0002);
        cs43 = (realnum)(0.415*phycon.te04*(phycon.te004/phycon.te0004));
        cs51 = (realnum)(0.1294*(phycon.te02/(phycon.te001*phycon.te0004)));
        cs52 = (realnum)(0.2841*phycon.te04*phycon.te0004);
        cs53 = (realnum)(0.2737*phycon.te04*phycon.te003*phycon.te0001);
        cs54 = (realnum)(0.2037*phycon.te04*(phycon.te002/phycon.te0004));
        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 = (realnum)(p[1]*a21*5.34e-12);
        CoolHeavy.O3726 = (realnum)(p[2]*a31*5.34e-12);
        CoolHeavy.O2471 = (realnum)((p[3]*a41 + p[4]*a51)*8.05e-12);
        CoolHeavy.O7323 = (realnum)((p[3]*a42 + p[4]*a52)*2.72e-12);
        CoolHeavy.O7332 = (realnum)((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. 
         * >>referold   o3      cs      Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
        /* >>chng 06 jun 30- Humeshkar Nemala */
        /*Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /*Data available in the temperature range 2500 K to 2E6 K*/
        /*fits at temperatures below 30000K and at temperatures above 30000K*/
        if(phycon.te < 30000)
        {
                cs = (realnum)(0.2519*(phycon.te07*phycon.te02*phycon.te007*phycon.te001*phycon.te0002));
        }
        else
        {
                cs = (realnum)(6.166388*(1/(phycon.te20*phycon.te01*phycon.te002)));
        }
        /*PutCS(0.756,&TauLines[ipT1666]);*/
        PutCS(cs ,&TauLines[ipT1666]);
        /*>>chng 06 jun 30- Humeshkar Nemala*/
        /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /*Data available in the temperature range 2500 K to 2E6 K*/
        /*fits at temperatures below 30000K and at temperatures above 30000K*/
        if(phycon.te < 30000)
        {
                cs=(realnum)((0.1511)*(phycon.te07*phycon.te02*phycon.te007*phycon.te001*phycon.te0002));
        }
        else
        {
                cs=(realnum)((3.668474)*(1/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002)));
        }
        /*PutCS(0.454,&TauLines[ipT1661]);*/
        PutCS(cs,&TauLines[ipT1661]);
        /*this line correspond to the transition 3P0-5S^o2*/
        /*>>chng 06 jun 30- Humeshkar Nemala*/
        /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /*Data available in the temperature range 2500 K to 2E6 K*/
        /*fits at temperatures below 30000K and at temperatures above 30000K*/
        if(phycon.te < 30000)
        {
                cs = (realnum)(0.0504*((phycon.te10/phycon.te01)*phycon.te005*phycon.te003*phycon.te0002));
        }
        else
        {
                cs = (realnum)(1.223633/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002));
        }
        /*PutCS(1.3,&TauDummy);*/
        PutCS(cs,&TauDummy);
        atom_level3(&TauDummy,&TauLines[ipT1666],&TauLines[ipT1661]);
        /* if( nLev3Fail.gt.0 )write(qq,'('' nlev3fail='',i3)') nLev3Fail
         *
         * OIII 304 not in cooling function */
        TauLines[ipT304].Emis->PopOpc = dense.xIonDense[ipOXYGEN][2];
        TauLines[ipT304].Lo->Pop = dense.xIonDense[ipOXYGEN][2];
        TauLines[ipT304].Hi->Pop = 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 */

        /* >>refer      o3      cs      Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /* >>chng 06 jul 24- Humeshkar Nemala */
        /*cs are available over two temperature ranges: below 30000K and above 30000K*/
        /*The old values seemed to saturate at around 100,000K.The new values are around 
        10-15% different from the old values*/
        /*The cs of the transitions 3P0,1,2 to 1D2 are added together to give oxy.o3cs12 */
        /*the cs of the transition 1D2-1S0 is mentioned as o3cs23*/
        /*The cs of the transitions 3P0,1,2 to 1S0 are added together to give oxy.o3cs13*/
        if(phycon.te < 3E4)
        {
                oxy.o3cs12 = (realnum)(0.9062*(phycon.te10/phycon.te002));
                o3cs23 = (realnum)(0.0995*phycon.te10*phycon.te07*(phycon.te007/phycon.te0002));
                oxy.o3cs13 = (realnum)(0.1237*phycon.te07*phycon.te02*phycon.te005*phycon.te0005);

        }
        else if(phycon.te < 6E4)
        {       
                oxy.o3cs12 = (realnum)(1.710073*(phycon.te04/phycon.te004)*phycon.te0004);
                oxy.o3cs13 = (realnum)(0.1963109*phycon.te05*phycon.te0007);
                o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001));
        }
        else
        {
                oxy.o3cs12 = (realnum)(6.452638/((phycon.te10/phycon.te02)*phycon.te004*phycon.te0003));
                oxy.o3cs13 = (realnum)(0.841578/(phycon.te07*phycon.te01*(phycon.te002/phycon.te0004)));
                o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001));
        }
#       if 0
        /*Old code commented -Humeshkar Nemala */
        if( phycon.te > 3981 && phycon.te <= 1e5 )
        {
                oxy.o3cs12 = (realnum)(3.0211144 - 101.57536/phycon.sqrte + 355.06905*
                        /*really is base e log of temperature */
                  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 = 
                (realnum)(MIN2(0.36,0.0784*phycon.te10*phycon.te03*phycon.te01*
          phycon.te003));
        /*End :Old Code Commented Out-Humeshkar Nemala*/
        /* >>chng 01 may 04, derive a21 from aeff, had been 0.02431 
        a21 = 0.02431;*/
#       endif

        a21 = aeff - oxy.d5007r;
        a31 = .215634;
        a32 = 1.71;
        oxy.o3ex23 = 32947.;
        oxy.o3br32 = (realnum)(a32/(a31 + a32));
        oxy.o3enro = (realnum)(4.56e-12/3.98e-12);
        /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
        oxy.poiii3 = 
                (realnum)(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 = (realnum)(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, 
         * >>referold   o3      cs      Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
        /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /*88 microns refers to the 3P0-3P1 transition*/
        if(phycon.te < 3E4)
        {
                cs = (realnum)(0.3993*(phycon.te03/phycon.te001)*phycon.te0001);
        }
        else if(phycon.te < 1E5)
        {
                cs = (realnum)(0.245712*phycon.te07*phycon.te005*phycon.te001*phycon.te0002);
        }
        else
        {
                cs = (realnum)(1.310467/((phycon.te07/phycon.te001)*phycon.te0002));
        }
        /*PutCS(0.5590,&TauLines[ipTO88]);*/
        PutCS(cs,&TauLines[ipTO88]);
        /*52 microns refers to the 3P1-3P2 transition*/
        if(phycon.te < 3E4)
        {
                cs = (realnum)(0.7812*(phycon.te05/phycon.te0001));
        }
        else if(phycon.te < 1.2E5)
        {
                cs = (realnum)(0.516157*phycon.te07*phycon.te02*phycon.te0001);
        }
        else
        {
                cs = (realnum)(4.038402/(phycon.te05*phycon.te03*phycon.te005*phycon.te0007*phycon.te0001));
        }
        /*PutCS(1.335,&TauLines[ipT52]);*/
        PutCS(cs,&TauLines[ipT52]);
        /*TauDummy refers to the 3P0-3P2 transition*/
        if(phycon.te < 3E4)
        {
                cs = (realnum)(0.1337*phycon.te07*phycon.te002*phycon.te0002);
        }
        else if(phycon.te < 1.2E5)
        {
                cs = (realnum)(0.086446*phycon.te10*phycon.te01*phycon.te004*phycon.te0005);
        }
        else
        {
                cs = (realnum)(0.82031/(phycon.te07*phycon.te007*phycon.te0007*phycon.te0002));
        }
        /*PutCS(0.2832,&TauDummy);*/
        PutCS(cs,&TauDummy);
        atom_level3(&TauLines[ipTO88],&TauLines[ipT52],&TauDummy);

        /* O III 834A, CS */
         /* >>referold  o3      cs      Aggarwal, K.M., 1985 A&A 146, 149. */
        /* >>chng 06 jul 22- Humeshkar Nemala */
        /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
        /*the cs of OIII 834A was obtained by summing the cs of the six transitions:
        3P0-3D^o1,3P1-3D^o1,3P2-3D^o1,3P1-3D^o2,3P2-3D^o2,3P2-3D^o3*/
        if(phycon.te < 3E4)
        {
                cs = (realnum)(0.9595*phycon.te04*phycon.te0002);
        }
        else
        {
                cs = (realnum)(0.10798*phycon.te20*phycon.te05*phycon.te002*phycon.te0001);
        }
        /*PutCS(5.,&TauLines[ipT835]);*/
        PutCS(cs,&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 micron, CS 
         * >>referold   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;*/

        /* >>chng 06 jun 30- Humeshkar Nemala */
        /*>>refer o4 cs Zhang,H.L.,Graziani,M. & Pradhan,A.K. 1994,A&A,283,319*/
        /*cs given over the temp range :900K to 450000K

        if(phycon.te < 1.8E4)
        {
                cs = (realnum)(0.5363*phycon.te10*phycon.te05*(phycon.te01/phycon.te0004));
        }
        else if(phycon.te < 5.4E4)
        {
                cs = (realnum)(1.554285*phycon.te05*phycon.te001);
        }
        else
        {
                cs = (realnum)(95.373669/(phycon.te30*phycon.te02*(phycon.te007/phycon.te0002)));
        }

         * >>chng 06 nov 08 - NPA.  Update collision strength to new data from:
         * >>refer      o4      cs      Tayal, S. 2006, ApJS 166, 634 
         * Equation derived by using TableCurve, and goes to zero as 
         * T => 0 and T => infinity */

        cs = (realnum)(exp(3.265723 - 0.00014686984*phycon.alnte*phycon.sqrte
                 -22.035066/phycon.alnte));
        /* cs goes to zero at very high T, which will cause an assert in the
         * n-level solver - don't let this happen */
        cs = MAX2( cs , 1e-4 );
        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].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4  )
                        {
                                ++nO4Pump;
                        }
                }
                if( nO4Pump<0 )
                        TotalInsanity();
                else if( nO4Pump > 0 )
                        /* create the space - can't malloc 0 bytes */
                        ipO4Pump = (long *)MALLOC((unsigned)(nO4Pump)*sizeof(long) );
                nO4Pump = 0;
                for( i=0; i<nWindLine; ++i )
                {
                        /* don't test on nelem==ipIRON since lines on physics, not C, scale */
                        if( TauLine2[i].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4  )
                        {
#                               if      0
                                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]].Emis->pump;
#               if      0
                fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
                        i,
                        TauLine2[ipO4Pump[i]].WLAng , TauLine2[ipO4Pump[i]].pump );
#               endif
        }

        /*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 = (realnum)(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]);
        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;

        DEBUG_ENTRY( "oi3Pcs()" );
        /* 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);
        }

        return;
}

---