00001
00002
00003
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
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
00058 DEBUG_EXIT( "IonOxyge()" );
00059 return;
00060 }
00061
00062 ion_zero(ipOXYGEN);
00063
00064 ion_photo(ipOXYGEN,false);
00065
00066
00067 ion_collis(ipOXYGEN);
00068
00069
00070
00071
00072 ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipOXYGEN);
00073
00074
00075
00078 oxy.p1666 = ionbal.PhotoRate_Shell[ipOXYGEN][3][1][0];
00079
00080 oxy.p1401 = ionbal.PhotoRate_Shell[ipOXYGEN][2][1][0];
00081
00082
00083
00084
00085
00086
00087 oxy.d5007r = (float)(GammaK(opac.ipo3exc[0],opac.ipo3exc[1],
00088 opac.ipo3exc[2] , 1. ));
00089
00090
00091 oxy.d4363 = (float)(GammaK(opac.ipo3exc3[0],opac.ipo3exc3[1],
00092 opac.ipo3exc3[2] , 1. ));
00093
00094
00095 oxy.d6300 = (float)(GammaK(opac.ipo1exc[0],opac.ipo1exc[1],
00096 opac.ipo1exc[2] , 1. ));
00097
00098
00099 aeff = 0.0263 + oxy.d5007r;
00100
00101
00102 oxy.poiii2 = (float)(atom_pop2(2.5,9.,5.,aeff,2.88e4,1.)/aeff);
00103 {
00104
00105 enum {DEBUG_LOC=false};
00106
00107 if( DEBUG_LOC )
00108 {
00109 fprintf(ioQQQ,"pop rel %.1e rate %.1e grnd rate %.1e\n",
00110 oxy.poiii2 , oxy.d5007r ,ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] );
00111 }
00112 }
00113
00114
00115 if( nzone > 0 )
00116 {
00117
00118 ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]*
00119 (1. - oxy.poiexc) + oxy.d6300*oxy.poiexc;
00120
00121
00122 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]*
00123 (1. - oxy.poiii2 - oxy.poiii3) + oxy.d5007r*oxy.poiii2 +
00124 oxy.d4363*oxy.poiii3;
00125
00126 if( ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > 1e-30 && dense.IonLow[ipOXYGEN] <= 2 )
00127 {
00128 if( (oxy.d5007r*oxy.poiii2 + oxy.d4363*oxy.poiii3)/
00129 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > (oxy.r4363Max +
00130 oxy.r5007Max) )
00131 {
00132 oxy.poiii2Max = (float)(oxy.d5007r*oxy.poiii2/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00133 oxy.poiii3Max = (float)(oxy.d4363*oxy.poiii3/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00134 }
00135 oxy.r4363Max = (float)(MAX2(oxy.r4363Max,oxy.d4363));
00136 oxy.r5007Max = (float)(MAX2(oxy.r5007Max,oxy.d5007r));
00137 }
00138
00139
00140 if( dense.IonLow[ipOXYGEN] <= 0 && (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] +
00141 atmdat.HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]) > 1e-30 )
00142 {
00143 oxy.poimax = (float)(MAX2(oxy.poimax,oxy.d6300*oxy.poiexc/
00144 (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]+
00145 atmdat.HCharExcIonOf[ipOXYGEN][0]* dense.xIonDense[ipHYDROGEN][1])));
00146 }
00147 }
00148 else
00149 {
00150 oxy.poiii2Max = 0.;
00151 oxy.poiii3Max = 0.;
00152 oxy.r4363Max = 0.;
00153 oxy.r5007Max = 0.;
00154 oxy.poimax = 0.;
00155 }
00156
00157
00158 if( dense.IonLow[ipOXYGEN] == 0 && oxy.i2d < rfield.nflux )
00159 {
00160 oxy.s3727 = (float)(GammaK(oxy.i2d,oxy.i2p,opac.iopo2d , 1. ));
00161
00162 iup = MIN2(iso.ipIsoLevNIonCon[ipH_LIKE][1][0],rfield.nflux);
00163 oxy.s7325 = (float)(GammaK(oxy.i2d,iup,opac.iopo2d , 1. ));
00164
00165 oxy.s7325 -= oxy.s3727;
00166 oxy.s3727 = oxy.s3727 + oxy.s7325;
00167
00168
00169 oxy.s7325 *= 0.66f;
00170 }
00171 else
00172 {
00173 oxy.s3727 = 0.;
00174 oxy.s7325 = 0.;
00175 }
00176
00177 oxy.AugerO3 = (float)ionbal.PhotoRate_Shell[ipOXYGEN][0][0][0];
00178
00179
00180
00181
00182
00183
00184
00185 save_rec = ionbal.RateRecomTot[ipOXYGEN][0];
00186
00187
00188
00189
00190 # if 0
00191 if( dense.IonLow[ipOXYGEN]==0 &&
00192 co.hevmol[ipOP] > SMALLFLOAT &&
00193 ionbal.RateIonizTot[ipOXYGEN][0]*co.hevmol[ipATO]>0. )
00194 {
00195 ionbal.RateRecomTot[ipOXYGEN][0] =
00196 ionbal.RateIonizTot[ipOXYGEN][0]*
00197 co.hevmol[ipATO]/co.hevmol[ipOP];
00198 }
00199
00200 # endif
00201
00202
00203 if(0 && nzone > 100 )
00204 lgDebug = true;
00205 else
00206 lgDebug = false;
00207 ion_solver(ipOXYGEN,lgDebug);
00208 if( lgDebug )
00209 fprintf(ioQQQ,"DEBUG O\t%.3e\t%.3e\tH\t%.3e\t%.3e\n",
00210 dense.xIonDense[ipOXYGEN][0],
00211 dense.xIonDense[ipOXYGEN][1],
00212 dense.xIonDense[ipHYDROGEN][0],
00213 dense.xIonDense[ipHYDROGEN][1]);
00214
00215
00216 if( save_rec > 0. )
00217 ionbal.RateRecomTot[ipOXYGEN][0] = save_rec;
00218
00219
00220 oxy.p1666 *= dense.xIonDense[ipOXYGEN][1]*0.3;
00221 oxy.p1401 *= dense.xIonDense[ipOXYGEN][2]*0.43;
00222 oxy.s3727 *= dense.xIonDense[ipOXYGEN][0];
00223 oxy.s7325 *= dense.xIonDense[ipOXYGEN][0];
00224 oxy.AugerO3 *= dense.xIonDense[ipOXYGEN][0];
00225
00226 if( trace.lgTrace )
00227 {
00228 fprintf( ioQQQ, " IonOxyge returns; frac=" );
00229 for( int i=1; i <= 9; i++ )
00230 {
00231 fprintf( ioQQQ, " %10.3e", dense.xIonDense[ipOXYGEN][i-1]/
00232 dense.gas_phase[ipOXYGEN] );
00233 }
00234 fprintf( ioQQQ, "\n" );
00235 }
00236
00237 DEBUG_EXIT( "IonOxyge()" );
00238 return;
00239 }
00240