00001
00002
00003
00004 #include "cddefines.h"
00005 #include "yield.h"
00006 #include "heavy.h"
00007 #include "opacity.h"
00008 #include "dense.h"
00009 #include "thermal.h"
00010 #include "conv.h"
00011 #include "grainvar.h"
00012 #include "elementnames.h"
00013 #include "gammas.h"
00014 #include "ionbal.h"
00015 #include "ca.h"
00016 #include "mole.h"
00017 #include "phycon.h"
00018 #include "hmi.h"
00019 #include "rfield.h"
00020 #include "atoms.h"
00021 #include "iso.h"
00022 #include "oxy.h"
00023 #include "atmdat.h"
00024 #include "fe.h"
00025
00026 void ion_photo(
00027
00028 long int nelem ,
00029
00030 bool lgPrintIt )
00031 {
00032 long int ion,
00033 iphi,
00034 iplow,
00035 ipop,
00036 limit_hi,
00037 limit_lo,
00038 ns;
00039
00040 DEBUG_ENTRY( "ion_photo()" );
00041
00042
00043
00044
00045 ASSERT( nelem < LIMELM );
00046 ASSERT( dense.IonLow[nelem] >= 0 );
00047 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
00048 ASSERT( dense.IonHigh[nelem] <= nelem + 1);
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 limit_hi = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );
00063
00064
00065 limit_hi = MAX2( 1 , limit_hi );
00066
00067
00068
00069 if( !conv.nPres2Ioniz && gv.lgDustOn() )
00070 {
00071 limit_lo = 0;
00072 }
00073 else
00074 {
00075 limit_lo = dense.IonLow[nelem];
00076 }
00077
00078
00079
00080
00081 for( ion=limit_lo; ion < limit_hi; ion++ )
00082 {
00083
00084 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00085 {
00086
00087 if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
00088 {
00089
00090 iplow = opac.ipElement[nelem][ion][ns][0];
00091 iphi = opac.ipElement[nelem][ion][ns][1];
00092 ipop = opac.ipElement[nelem][ion][ns][2];
00093
00094 t_phoHeat photoHeat;
00095
00096
00097
00098 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
00099 GammaK(iplow,iphi,
00100 ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
00101 &photoHeat )*ionbal.lgPhotoIoniz_On;
00102
00103
00104
00105
00106
00107 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = photoHeat.HeatLowEnr*ionbal.lgPhotoIoniz_On;
00108 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = photoHeat.HeatHiEnr*ionbal.lgPhotoIoniz_On;
00109 }
00110 }
00111
00112
00113
00114
00115 ns = (Heavy.nsShells[nelem][ion]-1);
00116
00117 ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ionbal.CompRecoilIonRate[nelem][ion];
00118
00119 ionbal.PhotoRate_Shell[nelem][ion][ns][2] += ionbal.CompRecoilHeatRate[nelem][ion];
00120 }
00121
00122
00123 if( lgPrintIt )
00124 {
00125
00126 ns = 5;
00127 ion = 1;
00128 GammaPrt(
00129 opac.ipElement[nelem][ion][ns][0],
00130 opac.ipElement[nelem][ion][ns][1],
00131 opac.ipElement[nelem][ion][ns][2],
00132 ioQQQ,
00133 ionbal.PhotoRate_Shell[nelem][ion][ns][0],
00134 0.05);
00135
00136
00137 for( ns=0; ns < Heavy.nsShells[nelem][0]; ns++ )
00138 {
00139 fprintf( ioQQQ, "\n %s", elementnames.chElementNameShort[nelem] );
00140 fprintf( ioQQQ, " %s" , Heavy.chShell[ns]);
00141
00142 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
00143 {
00144 if( Heavy.nsShells[nelem][ion] > ns )
00145 {
00146 fprintf( ioQQQ, " %8.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
00147 }
00148 else
00149 {
00150 break;
00151 }
00152 }
00153 }
00154 fprintf(ioQQQ,"\n");
00155 }
00156
00157
00158
00159 if( nelem == ipCALCIUM )
00160 {
00161 long ns = 6, ion = 1;
00162 ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ca.dstCala;
00163 }
00164 if( nelem == ipCARBON )
00165 {
00166
00167
00168 if(mole_global.lgLeidenHack)
00169 {
00170 int nelem=ipCARBON , ion=0 , ns=2;
00171 ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
00172 (HMRATE((1e-10)*3.0,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*
00173 exp(-(3.0*rfield.extin_mag_V_point))/1.66));
00174
00175 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
00176 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
00177 }
00178 }
00179 if( nelem == ipNITROGEN )
00180 {
00181
00182 if( dense.xIonDense[ipNITROGEN][0] > 0. )
00183 {
00184 t_phoHeat photoHeat;
00185
00186 atoms.d5200r = (realnum)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.,&photoHeat);
00187
00188 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]*
00189 (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r;
00190 ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]*
00191 (1. - atoms.p2nit) + photoHeat.HeatNet*atoms.p2nit;
00192 }
00193 else
00194 {
00195 atoms.p2nit = 0.;
00196 atoms.d5200r = 0.;
00197 }
00198 }
00199 if ( nelem == ipMAGNESIUM )
00200 {
00201 if( dense.IonLow[ipMAGNESIUM] <= 1 )
00202 {
00203 t_phoHeat dummy;
00204
00205 realnum rmg2l = (realnum)GammaK(opac.ipmgex,
00206 iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,opac.ipOpMgEx,1., &dummy );
00207 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] += rmg2l*atoms.popmg2;
00208
00209 if( nzone <= 1 )
00210 {
00211 atoms.xMg2Max = 0.;
00212 }
00213 else if( ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] > 1e-30 )
00214 {
00215
00216 atoms.xMg2Max = (realnum)(MAX2(atoms.xMg2Max,rmg2l*atoms.popmg2/
00217 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0]));
00218 }
00219 }
00220 else
00221 {
00222 atoms.xMg2Max = 0.;
00223 }
00224 }
00225
00226 if ( nelem == ipOXYGEN )
00227 {
00228
00229 t_phoHeat dummy;
00230
00231 if( !dense.lgElmtOn[ipOXYGEN] )
00232 {
00233 oxy.poiii2Max = 0.;
00234 oxy.poiii3Max = 0.;
00235 oxy.r4363Max = 0.;
00236 oxy.r5007Max = 0.;
00237 oxy.poiii2 = 0.;
00238 oxy.AugerO3 = 0.;
00239 oxy.s3727 = 0.;
00240 oxy.s7325 = 0.;
00241 thermal.heating[7][9] = 0.;
00242 oxy.poimax = 0.;
00243 return;
00244 }
00245 else
00246 {
00247 double aeff;
00248
00249
00250
00251
00252
00253
00254 oxy.d5007r = (realnum)(GammaK(opac.ipo3exc[0],opac.ipo3exc[1],
00255 opac.ipo3exc[2] , 1., &dummy ));
00256
00257
00258 oxy.d4363 = (realnum)(GammaK(opac.ipo3exc3[0],opac.ipo3exc3[1],
00259 opac.ipo3exc3[2] , 1., &dummy ));
00260
00261
00262 oxy.d6300 = (realnum)(GammaK(opac.ipo1exc[0],opac.ipo1exc[1],
00263 opac.ipo1exc[2] , 1., &dummy ));
00264
00265
00266 aeff = 0.0263 + oxy.d5007r;
00267
00268
00269 oxy.poiii2 = (realnum)(atom_pop2(2.5,9.,5.,aeff,2.88e4,1.)/aeff);
00270 {
00271 enum {DEBUG_LOC=false};
00272 if( DEBUG_LOC )
00273 {
00274 fprintf(ioQQQ,"pop rel %.1e rate %.1e grnd rate %.1e\n",
00275 oxy.poiii2 , oxy.d5007r ,ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] );
00276 }
00277 }
00278
00279
00280 if( nzone > 0 )
00281 {
00282
00283 ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]*
00284 (1. - oxy.poiexc) + oxy.d6300*oxy.poiexc;
00285
00286
00287 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]*
00288 (1. - oxy.poiii2 - oxy.poiii3) + oxy.d5007r*oxy.poiii2 +
00289 oxy.d4363*oxy.poiii3;
00290
00291 if( ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > 1e-30 && dense.IonLow[ipOXYGEN] <= 2 )
00292 {
00293 if( (oxy.d5007r*oxy.poiii2 + oxy.d4363*oxy.poiii3)/
00294 ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > (oxy.r4363Max +
00295 oxy.r5007Max) )
00296 {
00297 oxy.poiii2Max = (realnum)(oxy.d5007r*oxy.poiii2/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00298 oxy.poiii3Max = (realnum)(oxy.d4363*oxy.poiii3/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
00299 }
00300 oxy.r4363Max = (realnum)(MAX2(oxy.r4363Max,oxy.d4363));
00301 oxy.r5007Max = (realnum)(MAX2(oxy.r5007Max,oxy.d5007r));
00302 }
00303
00304
00305 if( dense.IonLow[ipOXYGEN] <= 0 && (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] +
00306 atmdat.HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]) > 1e-30 )
00307 {
00308 oxy.poimax = (realnum)(MAX2(oxy.poimax,oxy.d6300*oxy.poiexc/
00309 (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]+
00310 atmdat.HCharExcIonOf[ipOXYGEN][0]* dense.xIonDense[ipHYDROGEN][1])));
00311 }
00312 }
00313 else
00314 {
00315 oxy.poiii2Max = 0.;
00316 oxy.poiii3Max = 0.;
00317 oxy.r4363Max = 0.;
00318 oxy.r5007Max = 0.;
00319 oxy.poimax = 0.;
00320 }
00321 }
00322 long int iup;
00323
00324 if( dense.IonLow[ipOXYGEN] == 0 && oxy.i2d < rfield.nflux )
00325 {
00326 oxy.s3727 = (realnum)(GammaK(oxy.i2d,oxy.i2p,opac.iopo2d , 1., &dummy ));
00327
00328 iup = MIN2(iso_sp[ipH_LIKE][1].fb[0].ipIsoLevNIonCon,rfield.nflux);
00329 oxy.s7325 = (realnum)(GammaK(oxy.i2d,iup,opac.iopo2d , 1., &dummy ));
00330
00331 oxy.s7325 -= oxy.s3727;
00332 oxy.s3727 = oxy.s3727 + oxy.s7325;
00333
00334
00335 oxy.s7325 *= 0.66f;
00336 }
00337 else
00338 {
00339 oxy.s3727 = 0.;
00340 oxy.s7325 = 0.;
00341 }
00342
00343 oxy.AugerO3 = (realnum)ionbal.PhotoRate_Shell[ipOXYGEN][0][0][0];
00344
00345 oxy.s3727 *= dense.xIonDense[ipOXYGEN][0];
00346 oxy.s7325 *= dense.xIonDense[ipOXYGEN][0];
00347 oxy.AugerO3 *= dense.xIonDense[ipOXYGEN][0];
00348 }
00349 if( nelem == ipIRON )
00350 {
00351 if( !dense.lgElmtOn[ipIRON] )
00352 {
00353 fe.fekcld = 0.;
00354 fe.fekhot = 0.;
00355 fe.fegrain = 0.;
00356 }
00357 else
00358 {
00359 const int NDIM = ipIRON+1;
00360
00361 static const double fyield[NDIM+1] = {.34,.34,.35,.35,.36,.37,.37,.38,.39,.40,
00362 .41,.42,.43,.44,.45,.46,.47,.47,.48,.48,.49,.49,.11,.75,0.,0.,0.};
00363
00364 long int i, limit, limit2;
00365
00366
00367 fe.fekcld = 0.;
00368 limit = MIN2(18,dense.IonHigh[ipIRON]);
00369
00370 for( i=dense.IonLow[ipIRON]; i < limit; i++ )
00371 {
00372 ASSERT( i < NDIM + 1 );
00373 fe.fekcld +=
00374 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
00375 fyield[i]);
00376 }
00377
00378
00379 fe.fekhot = 0.;
00380 limit = MAX2(18,dense.IonLow[ipIRON]);
00381
00382 limit2 = MIN2(ipIRON+1,dense.IonHigh[ipIRON]);
00383 ASSERT( limit2 <= LIMELM + 1 );
00384
00385 for( i=limit; i < limit2; i++ )
00386 {
00387 ASSERT( i < NDIM + 1 );
00388 fe.fekhot +=
00389 (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
00390 fyield[i]);
00391 }
00392
00393
00394
00395 i = 0;
00396
00397 fe.fegrain = ( gv.lgWD01 ) ? 0.f : (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*fyield[i]*
00398 gv.elmSumAbund[ipIRON]);
00399 }
00400 }
00401
00402 return;
00403 }