00001
00002
00003 #include "cddefines.h"
00004 #include "atmdat.h"
00005 #include "conv.h"
00006 #include "dense.h"
00007 #include "heavy.h"
00008 #include "helike_cs.h"
00009 #include "hydrogenic.h"
00010 #include "hydro_vs_rates.h"
00011 #include "ionbal.h"
00012 #include "iso.h"
00013 #include "opacity.h"
00014 #include "phycon.h"
00015 #include "physconst.h"
00016 #include "rfield.h"
00017 #include "secondaries.h"
00018 #include "trace.h"
00019 #include "taulines.h"
00020
00021
00022 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
00023
00024 void iso_collisional_ionization( long ipISO, long nelem )
00025 {
00026 ASSERT( ipISO < NISO );
00027
00028 DEBUG_ENTRY( "iso_collisional_ionization()" );
00029
00030
00031 iso.ColIoniz[ipISO][nelem][0] = iso.lgColl_ionize[ipISO] *
00032 t_ADfA::Inst().coll_ion( nelem+1, 1+ipISO, phycon.te );
00033
00034 iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],0,IPCOLLIS,0.20f,0.20f);
00035
00036 for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00037 {
00038 if( nelem == ipISO )
00039 {
00040
00041
00042 iso.ColIoniz[ipISO][nelem][ipHi] = hydro_vs_ioniz( iso.xIsoLevNIonRyd[ipISO][nelem][ipHi], phycon.te );
00043 }
00044 else
00045 {
00046
00047
00048
00049
00050
00051 iso.ColIoniz[ipISO][nelem][ipHi] =
00052 Hion_coll_ioniz_ratecoef( ipISO, nelem, N_(ipHi), iso.xIsoLevNIonRyd[ipISO][nelem][ipHi], phycon.te );
00053 }
00054
00055
00056
00057 if( ipHi < iso.numLevels_max[ipISO][nelem] - 1 )
00058 iso.ColIoniz[ipISO][nelem][ipHi] *= iso.lgColl_ionize[ipISO];
00059
00060 iso_put_error(ipISO,nelem,iso.numLevels_max[ipISO][nelem],ipHi,IPCOLLIS,0.20f,0.20f);
00061 }
00062
00063
00064
00065
00066
00067
00068
00069
00070 if( !iso.lgLevelsLowered[ipISO][nelem] )
00071 {
00072 iso.ColIoniz[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] *= 100.;
00073 }
00074
00075 return;
00076 }
00077
00078 void iso_suprathermal( long ipISO, long nelem )
00079 {
00080 DEBUG_ENTRY( "iso_suprathermal()" );
00081
00082
00083 ASSERT( ipISO < NISO );
00084 ASSERT( nelem >= ipISO );
00085 ASSERT( nelem < LIMELM );
00086
00087
00088
00089
00090
00091
00092
00093 for( long i=1; i < iso.numLevels_max[ipISO][nelem]; i++ )
00094 {
00095 if( Transitions[ipISO][nelem][i][0].ipCont > 0 )
00096 {
00097
00098
00099
00100
00101 secondaries.Hx12[ipISO][nelem][i] = secondaries.x12tot *
00102 (Transitions[ipISO][nelem][i][0].Emis->gf/
00103 Transitions[ipISO][nelem][i][0].EnergyWN) /
00104 (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].Emis->gf/
00105 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].EnergyWN);
00106 }
00107 else
00108 secondaries.Hx12[ipISO][nelem][i] = 0.;
00109 }
00110
00111 return;
00112 }
00113
00114
00115
00116 void iso_collide( long ipISO, long nelem )
00117 {
00118 long ipHi, ipLo;
00119 double factor, ConvLTEPOP;
00120
00121
00122
00123 static double TeUsed[NISO][LIMELM]={
00124 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
00125 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
00126
00127
00128
00129
00130 static double TeUsedForCS[NISO][LIMELM]={
00131 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
00132 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
00133
00134 DEBUG_ENTRY( "iso_collide()" );
00135
00136
00137 ASSERT( ipISO < NISO );
00138 ASSERT( nelem >= ipISO );
00139 ASSERT( nelem < LIMELM );
00140
00141
00142
00143
00144 if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz && !iso.lgMustReeval[ipISO][nelem] )
00145 {
00146 ASSERT( Transitions[ipISO][nelem][ iso.nLyaLevel[ipISO] ][0].Coll.ColUL >= 0. );
00147
00148 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
00149 {
00150 fprintf( ioQQQ,
00151 " iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n",
00152 iso.chISO[ipISO], nelem );
00153 }
00154 }
00155 else
00156 {
00157 TeUsed[ipISO][nelem] = phycon.te;
00158
00159 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
00160 {
00161 fprintf( ioQQQ,
00162 " iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n",
00163 iso.chISO[ipISO], nelem );
00164 }
00165
00166
00167
00168
00169
00170
00171
00172
00173 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00174 {
00175 for( ipLo=0; ipLo<ipHi; ipLo++ )
00176 {
00177 iso.Boltzmann[ipISO][nelem][ipHi][ipLo] =
00178 sexp( Transitions[ipISO][nelem][ipHi][ipLo].EnergyK / phycon.te );
00179 }
00180 }
00181
00182
00183 factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
00184 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
00185
00186
00187 ConvLTEPOP = pow(factor,1.5)/(2.*iso.stat_ion[ipISO])/phycon.te32;
00188
00189 iso.lgPopLTE_OK[ipISO][nelem] = true;
00190
00191
00192
00193
00194 #define MAX_POP_LTE (MAX_DENSITY/dense.density_low_limit/dense.density_low_limit)
00195
00196
00197 for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]; ipLo++ )
00198 {
00199
00200 StatesElemNEW[nelem][nelem-ipISO][ipLo].ConBoltz =
00201 dsexp(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]/phycon.te_ryd);
00202
00203
00204
00205
00206
00207
00208
00209 if( StatesElemNEW[nelem][nelem-ipISO][ipLo].ConBoltz > SMALLDOUBLE )
00210 {
00211
00212 iso.PopLTE[ipISO][nelem][ipLo] =
00213 StatesElemNEW[nelem][nelem-ipISO][ipLo].g / StatesElemNEW[nelem][nelem-ipISO][ipLo].ConBoltz * ConvLTEPOP;
00214 ASSERT( iso.PopLTE[ipISO][nelem][ipLo] < BIGDOUBLE );
00215 }
00216 else
00217 {
00218 iso.PopLTE[ipISO][nelem][ipLo] = 0.;
00219 }
00220
00221 iso.PopLTE[ipISO][nelem][ipLo] = MIN2( iso.PopLTE[ipISO][nelem][ipLo], MAX_POP_LTE );
00222
00223
00224 if( iso.PopLTE[ipISO][nelem][ipLo] <= 0. )
00225 {
00226 iso.lgPopLTE_OK[ipISO][nelem] = false;
00227 }
00228 }
00229
00230 iso_collisional_ionization( ipISO, nelem );
00231
00232
00233
00234
00235
00236
00237
00238 if( iso.lgColl_excite[ipISO] &&
00239 ( iso.lgRandErrGen[ipISO] ||
00240 ! fp_equal( phycon.te, TeUsedForCS[ipISO][nelem] ) ) )
00241 {
00242
00243 TeUsedForCS[ipISO][nelem] = phycon.te;
00244
00245 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00246 {
00247 for( ipLo=0; ipLo < ipHi; ipLo++ )
00248 {
00249 for( long ipCollider = ipELECTRON; ipCollider <= ipALPHA; ipCollider++ )
00250 {
00251 double cs_temp = 0.;
00252
00253 if( N_(ipHi) == N_(ipLo) && !iso.lgColl_l_mixing[ipISO] )
00254 cs_temp = 0.;
00255 else if( N_(ipHi)-N_(ipLo) > 2 && ipCollider > ipELECTRON )
00256 cs_temp = 0.;
00257 else if( ipISO == ipH_LIKE )
00258 cs_temp = HydroCSInterp( nelem , ipHi , ipLo, ipCollider );
00259 else if( ipISO == ipHE_LIKE )
00260 cs_temp = HeCSInterp( nelem , ipHi , ipLo, ipCollider );
00261 else
00262 TotalInsanity();
00263
00264 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipCollider] = (realnum) cs_temp;
00265
00266
00267 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipELECTRON] >= 0. );
00268 }
00269
00270 if( N_(ipHi) <= 5 && N_(ipLo) <= 2 )
00271 iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.10f, 0.30f );
00272 else
00273 iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.20f, 0.30f );
00274
00275
00276 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_str =
00277 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipELECTRON];
00278
00279 }
00280 }
00281 }
00282
00283
00284
00285
00286
00287
00288 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00289 {
00290 for( ipLo=0; ipLo<ipHi; ipLo++ )
00291 {
00292 double reduced_mass_proton = dense.AtomicWeight[nelem]*ColliderMass[ipPROTON]/
00293 (dense.AtomicWeight[nelem]+ColliderMass[ipPROTON])*ATOMIC_MASS_UNIT;
00294
00295 double reduced_mass_heplus = dense.AtomicWeight[nelem]*ColliderMass[ipHE_PLUS]/
00296 (dense.AtomicWeight[nelem]+ColliderMass[ipHE_PLUS])*ATOMIC_MASS_UNIT;
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL = (realnum)(
00312 (
00313
00314 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_str +
00315
00316 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipPROTON]*
00317 (realnum)(dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr)*
00318 pow(ELECTRON_MASS/reduced_mass_proton, 1.5) +
00319
00320 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipHE_PLUS]*
00321 (realnum)(dense.xIonDense[ipHELIUM][1]/dense.EdenHCorr)*
00322 pow(ELECTRON_MASS/reduced_mass_heplus, 1.5) +
00323
00324 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipALPHA]*
00325 (realnum)(dense.xIonDense[ipHELIUM][2]/dense.EdenHCorr)*
00326 pow(ELECTRON_MASS/reduced_mass_heplus, 1.5)
00327 ) / phycon.sqrte*COLL_CONST/(double)StatesElemNEW[nelem][nelem-ipISO][ipHi].g );
00328
00329 if( ipISO == ipH_LIKE )
00330 {
00331 if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] &&
00332 N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] )
00333 {
00334 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *=
00335 (8.f/3.f)*(log((realnum)N_(ipHi))+2.f);
00336 }
00337 }
00338 else if( ipISO == ipHE_LIKE )
00339 {
00340 fixit();
00341
00342
00343
00344 if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] &&
00345 N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] )
00346 {
00347 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *=
00348 (2.f/3.f)*(log((realnum)N_(ipHi))+2.f);
00349 }
00350 }
00351
00352 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *= iso.lgColl_excite[ipISO];
00353 }
00354 }
00355
00356 if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
00357 {
00358 fprintf( ioQQQ, " iso_collide: %s Z=%li de-excitation rates coefficients\n", iso.chISO[ipISO], nelem + 1 );
00359 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00360 {
00361 fprintf( ioQQQ, " %li\t", ipHi );
00362 for( ipLo=0; ipLo < ipHi; ipLo++ )
00363 {
00364 fprintf( ioQQQ,PrintEfmt("%10.3e", Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL ));
00365 }
00366 fprintf( ioQQQ, "\n" );
00367 }
00368
00369 fprintf( ioQQQ, " iso_collide: %s Z=%li collisional ionization coefficients\n", iso.chISO[ipISO], nelem + 1 );
00370 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00371 {
00372 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ColIoniz[ipISO][nelem][ipHi] ));
00373 }
00374 fprintf( ioQQQ, "\n" );
00375
00376 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 );
00377 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00378 {
00379 fprintf( ioQQQ,PrintEfmt("%10.3e", StatesElemNEW[nelem][nelem-ipISO][ipHi].ConBoltz ));
00380 }
00381 fprintf( ioQQQ, "\n" );
00382
00383 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 );
00384 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00385 {
00386 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.PopLTE[ipISO][nelem][ipHi] ));
00387 }
00388 fprintf( ioQQQ, "\n" );
00389 }
00390
00391
00392
00393 if( opac.lgCaseB_HummerStorey )
00394 {
00395 for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]-1; ipLo++ )
00396 {
00397 if( N_(ipLo)>=3 )
00398 break;
00399
00400 iso.ColIoniz[ipISO][nelem][ipLo] = 0.;
00401
00402 for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00403 {
00404
00405 if( N_(ipLo)==2 && N_(ipHi)==2 )
00406 continue;
00407
00408 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL = 0.;
00409 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_str = 0.;
00410 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipPROTON] = 0.;
00411 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipHE_PLUS] = 0.;
00412 Transitions[ipISO][nelem][ipHi][ipLo].Coll.col_stri[ipALPHA] = 0.;
00413 }
00414 }
00415 }
00416 }
00417
00418 iso_suprathermal( ipISO, nelem );
00419
00420
00421
00422 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] =
00423 iso.ColIoniz[ipISO][nelem][0]*dense.EdenHCorr;
00424
00425
00426 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] =
00427 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]*
00428 rfield.anu[Heavy.ipHeavy[nelem][nelem-ipISO]-1]*EN1RYD;
00429
00430 return;
00431 }