cloudy
trunk
|
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 }