cloudy: /home66/gary/public_html/cloudy/c08_branch/source/ion

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*IonOxyge derive ionization balance for oxygen */
00004 #include "cddefines.h"
00005 #include "opacity.h"
00006 #include "oxy.h"
00007 #include "thermal.h"
00008 #include "dense.h"
00009 #include "iso.h"
00010 #include "trace.h"
00011 #include "rfield.h"
00012 #include "atmdat.h"
00013 #include "atoms.h"
00014 #include "gammas.h"
00015 #include "ionbal.h"
00016 
00017 void IonOxyge(void)
00018 {
00019         const int NDIM = ipOXYGEN+1;
00020 
00021         static const double dicoef[2][NDIM] = {
00022                 {1.11e-3,5.07e-3,1.48e-2,1.84e-2,4.13e-3,1.06e-1,6.23e-2,0.},
00023                 {.0925,.181,.305,.1,.162,.34,.304,0.}
00024         };
00025         static const double dite[2][NDIM] = {
00026                 {1.75e5,1.98e5,2.41e5,2.12e5,1.25e5,6.25e6,7.01e6,0.},
00027                 {1.45e5,3.35e5,2.83e5,2.83e5,2.27e5,1.12e6,1.47e6,0.}
00028         };
00029         static const double ditcrt[NDIM] = {2.7e4,2.2e4,2.4e4,2.5e4,1.6e4,1.0e6,1.5e6,1e20};
00030         static const double aa[NDIM] = {0.,-0.0036,0.,0.0061,-2.8425,0.,0.,0.};
00031         static const double bb[NDIM] = {0.0238,0.7519,21.8790,0.2269,0.2283,0.,0.,0.};
00032         static const double cc[NDIM] = {0.0659,1.5252,16.2730,32.1419,40.4072,0.,0.,0.};
00033         static const double dd[NDIM] = {0.0349,-0.0838,-0.7020,1.9939,-3.4956,0.,0.,0.};
00034         static const double ff[NDIM] = {0.5334,0.2769,1.1899,-0.0646,1.7558,0.,0.,0.};
00035 
00036         bool lgDebug = false;
00037         long int iup;
00038         double aeff, save_rec;
00039 
00040         DEBUG_ENTRY( "IonOxyge()" );
00041 
00042         /* oxygen, atomic number 8 */
00043         if( !dense.lgElmtOn[ipOXYGEN] )
00044         {
00045                 oxy.poiii2Max = 0.;
00046                 oxy.poiii3Max = 0.;
00047                 oxy.r4363Max = 0.;
00048                 oxy.r5007Max = 0.;
00049                 oxy.poiii2 = 0.;
00050                 oxy.p1666 = 0.;
00051                 oxy.AugerO3 = 0.;
00052                 oxy.p1401 = 0.;
00053                 oxy.s3727 = 0.;
00054                 oxy.s7325 = 0.;
00055                 thermal.heating[7][9] = 0.;
00056                 oxy.poimax = 0.;
00057                 return;
00058         }
00059 
00060         ion_zero(ipOXYGEN);
00061 
00062         ion_photo(ipOXYGEN,false);
00063 
00064         /* find collisional ionization rates */
00065         ion_collis(ipOXYGEN);
00066 
00067         /* get recombination coefficients */
00068         /*lint -e64 type mismatch */
00069         ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipOXYGEN);
00070         /*lint +e64 type mismatch */
00071 
00072         /* photoexcitation of O III 1666 and O IV 1401 */
00075         oxy.p1666 = ionbal.PhotoRate_Shell[ipOXYGEN][3][1][0];
00076 
00077         oxy.p1401 = ionbal.PhotoRate_Shell[ipOXYGEN][2][1][0];
00078 
00079         /* photoionization from O++ 1D
00080          *
00081          * estimate gamma function by assuming no frequency dependence
00082          * betwen 1D and O++3P edge */
00083         /* destroy upper level of OIII 5007*/
00084         oxy.d5007r = (realnum)(GammaK(opac.ipo3exc[0],opac.ipo3exc[1],
00085           opac.ipo3exc[2] , 1. ));
00086 
00087         /* destroy upper level of OIII 4363*/
00088         oxy.d4363 = (realnum)(GammaK(opac.ipo3exc3[0],opac.ipo3exc3[1],
00089           opac.ipo3exc3[2] , 1. ));
00090 
00091         /* destroy upper level of OI 6300*/
00092         oxy.d6300 = (realnum)(GammaK(opac.ipo1exc[0],opac.ipo1exc[1],
00093           opac.ipo1exc[2] , 1. ));
00094 
00095         /* A21 = 0.0263 */
00096         aeff = 0.0263 + oxy.d5007r;
00097 
00098         /* 1. as last arg makes this the relative population */
00099         oxy.poiii2 = (realnum)(atom_pop2(2.5,9.,5.,aeff,2.88e4,1.)/aeff);
00100         {
00101                 /*@-redef@*/
00102                 enum {DEBUG_LOC=false};
00103                 /*@+redef@*/
00104                 if( DEBUG_LOC )
00105                 {
00106                         fprintf(ioQQQ,"pop rel  %.1e rate %.1e  grnd rate %.1e\n", 
00107                                 oxy.poiii2 , oxy.d5007r ,ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] );
00108                 }
00109         }
00110 
00111         /* photoionization from excited states */
00112         if( nzone > 0 ) 
00113         {
00114                 /* neutral oxygen destruction */
00115                 ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]*
00116                   (1. - oxy.poiexc) + oxy.d6300*oxy.poiexc;
00117 
00118                 /* doubly ionized oxygen destruction */
00119                 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]*
00120                   (1. - oxy.poiii2 - oxy.poiii3) + oxy.d5007r*oxy.poiii2 + 
00121                   oxy.d4363*oxy.poiii3;
00122 
00123                 if( ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > 1e-30 && dense.IonLow[ipOXYGEN] <= 2 )
00124                 {
00125                         if( (oxy.d5007r*oxy.poiii2 + oxy.d4363*oxy.poiii3)/
00126                           ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > (oxy.r4363Max + 
00127                           oxy.r5007Max) )
00128                         {
00129                                 oxy.poiii2Max = (realnum)(oxy.d5007r*oxy.poiii2/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00130                                 oxy.poiii3Max = (realnum)(oxy.d4363*oxy.poiii3/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00131                         }
00132                         oxy.r4363Max = (realnum)(MAX2(oxy.r4363Max,oxy.d4363));
00133                         oxy.r5007Max = (realnum)(MAX2(oxy.r5007Max,oxy.d5007r));
00134                 }
00135 
00136                 /* ct into excited states */
00137                 if( dense.IonLow[ipOXYGEN] <= 0 && (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] + 
00138                   atmdat.HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]) > 1e-30 )
00139                 {
00140                         oxy.poimax = (realnum)(MAX2(oxy.poimax,oxy.d6300*oxy.poiexc/
00141                           (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]+
00142                           atmdat.HCharExcIonOf[ipOXYGEN][0]* dense.xIonDense[ipHYDROGEN][1])));
00143                 }
00144         }
00145         else
00146         {
00147                 oxy.poiii2Max = 0.;
00148                 oxy.poiii3Max = 0.;
00149                 oxy.r4363Max = 0.;
00150                 oxy.r5007Max = 0.;
00151                 oxy.poimax = 0.;
00152         }
00153 
00154         /* save atomic oxygen photodistruction rate for 3727 creation */
00155         if( dense.IonLow[ipOXYGEN] == 0 && oxy.i2d < rfield.nflux )
00156         {
00157                 oxy.s3727 = (realnum)(GammaK(oxy.i2d,oxy.i2p,opac.iopo2d , 1. ));
00158 
00159                 iup = MIN2(iso.ipIsoLevNIonCon[ipH_LIKE][1][0],rfield.nflux);
00160                 oxy.s7325 = (realnum)(GammaK(oxy.i2d,iup,opac.iopo2d , 1. ));
00161 
00162                 oxy.s7325 -= oxy.s3727;
00163                 oxy.s3727 = oxy.s3727 + oxy.s7325;
00164 
00165                 /* ratio of cross sections */
00166                 oxy.s7325 *= 0.66f;
00167         }
00168         else
00169         {
00170                 oxy.s3727 = 0.;
00171                 oxy.s7325 = 0.;
00172         }
00173 
00174         oxy.AugerO3 = (realnum)ionbal.PhotoRate_Shell[ipOXYGEN][0][0][0];
00175 
00176         /* >>chng 03 sep 29, synch up ion and co solvers.
00177          * the co solver will have a different O/O+ balance that
00178          * is strongly affected by the chemistry if we are in neutral gas
00179          * (hence the test that the highest stage of ionization is <=2 on 
00180          * physics scale)
00181          * in this case use co solver's O/O+ balance */
00182         save_rec = ionbal.RateRecomTot[ipOXYGEN][0];
00183         /*>>chng 04 apr 27, do not test for ionhigh being 1,
00184          * no reason for test on upper stage of ionization,
00185          * if codrive called but not evaluted then hevmol is all zero */
00186         /* >>chng 04 sep 10, rm check on search phase, no reason for it */
00187 #       if 0
00188         if( dense.IonLow[ipOXYGEN]==0  && 
00189                 co.hevmol[ipOP] > SMALLFLOAT && 
00190                 ionbal.RateIonizTot[ipOXYGEN][0]*co.hevmol[ipATO]>0. )
00191         {
00192                 ionbal.RateRecomTot[ipOXYGEN][0] = 
00193                         ionbal.RateIonizTot[ipOXYGEN][0]*
00194                         co.hevmol[ipATO]/co.hevmol[ipOP];
00195         }
00196 
00197 #       endif
00198 
00199         /* solve for ionization balance */
00200         if(0 &&  nzone > 100 )
00201                 lgDebug = true;
00202         else
00203                 lgDebug = false;
00204         ion_solver(ipOXYGEN,lgDebug);
00205         if( lgDebug )
00206                 fprintf(ioQQQ,"DEBUG O\t%.3e\t%.3e\tH\t%.3e\t%.3e\n",
00207                 dense.xIonDense[ipOXYGEN][0],
00208                 dense.xIonDense[ipOXYGEN][1],
00209                 dense.xIonDense[ipHYDROGEN][0],
00210                 dense.xIonDense[ipHYDROGEN][1]);
00211 
00212         /* reset the var we just hosed */
00213         if( save_rec > 0. )
00214                 ionbal.RateRecomTot[ipOXYGEN][0] = save_rec;
00215 
00216         /* 1666 ratio corrected for phot crs at 50ev */
00217         oxy.p1666 *= dense.xIonDense[ipOXYGEN][1]*0.3;
00218         oxy.p1401 *= dense.xIonDense[ipOXYGEN][2]*0.43;
00219         oxy.s3727 *= dense.xIonDense[ipOXYGEN][0];
00220         oxy.s7325 *= dense.xIonDense[ipOXYGEN][0];
00221         oxy.AugerO3 *= dense.xIonDense[ipOXYGEN][0];
00222 
00223         if( trace.lgTrace )
00224         {
00225                 fprintf( ioQQQ, "     IonOxyge returns; frac=" );
00226                 for( int i=1; i <= 9; i++ )
00227                 {
00228                         fprintf( ioQQQ, " %10.3e", dense.xIonDense[ipOXYGEN][i-1]/
00229                           dense.gas_phase[ipOXYGEN] );
00230                 }
00231                 fprintf( ioQQQ, "\n" );
00232         }
00233         return;
00234 }
/home66/gary/public_html/cloudy/c08_branch/source/ion_oxyge.cpp
