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