00001
00002
00003
00004
00005
00006
00007
00008 #include "cddefines.h"
00009 #include "dynamics.h"
00010 #include "atmdat.h"
00011 #include "conv.h"
00012 #include "cosmology.h"
00013 #include "dense.h"
00014 #include "elementnames.h"
00015 #include "helike.h"
00016 #include "helike_recom.h"
00017 #include "hydrogenic.h"
00018 #include "ionbal.h"
00019 #include "iso.h"
00020 #include "opacity.h"
00021 #include "phycon.h"
00022 #include "physconst.h"
00023 #include "prt.h"
00024 #include "save.h"
00025 #include "thermal.h"
00026 #include "thirdparty.h"
00027 #include "trace.h"
00028 #include "rt.h"
00029 #include "taulines.h"
00030
00031
00032
00033 static double ****RRCoef;
00034 static long **NumLevRecomb;
00035 static double ***TotalRecomb;
00036
00037
00038 static double TeRRCoef[N_ISO_TE_RECOMB];
00039 static double kTRyd,global_EthRyd;
00040 static long int globalZ,globalISO;
00041 static long int globalN, globalL, globalS;
00042
00043 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements );
00044 STATIC double iso_recomb_integrand(double EE);
00045 STATIC void iso_put_recomb_error( long ipISO, long nelem );
00046
00047 double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
00048 {
00049 double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
00050 double change[5] = {0.,0.,0.,0.,0.};
00051
00052 DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" );
00053
00054 if( ipISO==ipH_LIKE && ipLo == 0 )
00055 return t_ADfA::Inst().H_rad_rec(nelem+1,ipLo,temp);
00056
00057 global_EthRyd = iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd;
00058
00059
00060 b = MILNE_CONST * iso_sp[ipISO][nelem].st[ipLo].g() * pow(temp,-1.5);
00061
00062 if( ipISO==ipH_LIKE )
00063 b /= 2.;
00064 else if( ipISO==ipHE_LIKE )
00065 b /= 4.;
00066
00067
00068 kTRyd = temp / TE1RYD;
00069 globalISO = ipISO;
00070 globalZ = nelem;
00071 globalN = N_(ipLo);
00072 globalL = L_(ipLo);
00073 globalS = S_(ipLo);
00074
00075
00076
00077 E1 = global_EthRyd;
00078
00079 if( ipISO==ipH_LIKE )
00080 step = MIN2( 0.125*kTRyd, 0.5*E1 );
00081 else if( ipISO==ipHE_LIKE )
00082 step = MIN2( 0.25*kTRyd, 0.5*E1 );
00083 else
00084 TotalInsanity();
00085
00086 E2 = E1 + step;
00087
00088 RecomIntegral = qg32( E1, E2, iso_recomb_integrand);
00089
00090
00091
00092 do
00093 {
00094 OldRecomIntegral = RecomIntegral;
00095 E1 = E2;
00096 step *= 1.25;
00097 E2 = E1 + step;
00098 RecomIntegral += qg32( E1, E2, iso_recomb_integrand);
00099 change[4] = change[3];
00100 change[3] = change[2];
00101 change[2] = change[1];
00102 change[1] = change[0];
00103 change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
00104 TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
00105
00106
00107
00108
00109 } while ( ((E2-global_EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );
00110
00111
00112 alpha = b * RecomIntegral;
00113
00114 alpha = MAX2( alpha, SMALLDOUBLE );
00115
00116 return alpha;
00117 }
00118
00119
00120 STATIC double iso_recomb_integrand(double ERyd)
00121 {
00122 double x1, temp;
00123
00124
00125 x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - global_EthRyd ) / kTRyd);
00126 temp = iso_cross_section( ERyd , global_EthRyd, globalN, globalL, globalS, globalZ, globalISO );
00127 x1 *= temp;
00128
00129 return x1;
00130 }
00131
00132 double iso_cross_section( double EgammaRyd , double EthRyd, long n, long l, long S, long globalZ, long globalISO )
00133 {
00134 double cross_section;
00135 DEBUG_ENTRY( "iso_cross_section()" );
00136
00137 if( globalISO == ipH_LIKE )
00138 cross_section = H_cross_section( EgammaRyd , EthRyd, n, l, globalZ );
00139 else if( globalISO == ipHE_LIKE )
00140 cross_section = He_cross_section( EgammaRyd , EthRyd, n, l, S, globalZ );
00141 else
00142 TotalInsanity();
00143
00144 return cross_section;
00145 }
00146
00147
00148
00149 void iso_radiative_recomb(
00150 long ipISO,
00151
00152 long nelem )
00153 {
00154 long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel;
00155 double topoff, TotMinusExplicitResolved,
00156 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
00157 double RecExplictLevels, TotalRadRecomb, RecCollapsed;
00158 static double TeUsed[NISO][LIMELM]={
00159 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00160 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00161 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00162 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00163 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00164 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
00165
00166 DEBUG_ENTRY( "iso_radiative_recomb()" );
00167
00168 iso_put_recomb_error( ipISO, nelem );
00169
00170
00171
00172
00173
00174 RecExplictLevels = 0.;
00175
00176
00177 ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_local-iso_sp[ipISO][nelem].nCollapsed_local;
00178
00179 ASSERT( ipFirstCollapsed == iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local );
00180 if( !fp_equal(phycon.te,TeUsed[ipISO][nelem]) || iso_sp[ipISO][nelem].lgMustReeval || !conv.nTotalIoniz || cosmology.lgDo )
00181 {
00182 TeUsed[ipISO][nelem] = phycon.te;
00183
00184 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
00185 {
00186
00187 double RadRec;
00188
00189 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
00190 {
00191 RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel );
00192 }
00193 else
00194 {
00195 RadRec = iso_radrecomb_from_cross_section( ipISO, phycon.te, nelem, ipLevel);
00196 }
00197 ASSERT( RadRec > 0. );
00198
00199 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] = RadRec;
00200
00201 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] > 0. );
00202
00203 RecExplictLevels += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
00204
00205 if( iso_ctrl.lgDielRecom[ipISO] )
00206 {
00207
00208 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
00209 Total_DR_Added += iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
00210 }
00211 }
00212
00213
00214
00215
00216 RecCollapsed = 0.;
00217 for( ipLevel=ipFirstCollapsed; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
00218 {
00219
00220 double RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, N_(ipLevel), phycon.te);
00221
00222
00223 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] = RadRec;
00224
00225 RecCollapsed += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
00226
00227 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] > 0. );
00228
00229 if( iso_ctrl.lgDielRecom[ipISO] )
00230 {
00231
00232 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
00233 Total_DR_Added += iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
00234 }
00235 }
00236 for( ipLevel=iso_sp[ipISO][nelem].numLevels_local; ipLevel<iso_sp[ipISO][nelem].numLevels_max;++ipLevel )
00237 {
00238 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = 0.;
00239 }
00240
00241 for( ipLevel = 0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00242 {
00243 if( N_(ipLevel) == ThisN )
00244 {
00245 TotRRThisN += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
00246 }
00247 else
00248 {
00249 ASSERT( iso_ctrl.lgRandErrGen[ipISO] || (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) || (ThisN==iso_sp[ipISO][nelem].n_HighestResolved_local+1) );
00250 LastN = ThisN;
00251 ThisN = N_(ipLevel);
00252 TotRRLastN = TotRRThisN;
00253 TotRRThisN = iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
00254
00255 {
00256
00257
00258 enum {DEBUG_LOC=false};
00259
00260 static long RUNONCE = false;
00261
00262 if( !RUNONCE && DEBUG_LOC )
00263 {
00264 static long FIRSTTIME = true;
00265
00266 if( FIRSTTIME )
00267 {
00268 fprintf( ioQQQ,"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n",
00269 ipISO, nelem, phycon.te);
00270 FIRSTTIME= false;
00271 }
00272
00273 fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN );
00274 }
00275 RUNONCE = true;
00276 }
00277 }
00278 }
00279
00280
00281 if( iso_ctrl.lgNoRecombInterp[ipISO] )
00282 {
00283
00284 TotalRadRecomb = RecCollapsed + RecExplictLevels;
00285
00286
00287 for( long nLo = N_(ipFirstCollapsed) + iso_sp[ipISO][nelem].nCollapsed_max; nLo < NHYDRO_MAX_LEVEL; nLo++ )
00288 {
00289 TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te);
00290 }
00291
00292 for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ )
00293 {
00294 TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo );
00295 }
00296 }
00297 else if( iso_sp[ipISO][nelem].lgLevelsLowered )
00298 {
00299 TotalRadRecomb = 0.;
00300 for( ipLevel = 0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00301 TotalRadRecomb += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
00302 }
00303 else
00304 {
00305
00306 TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max );
00307 }
00308 if(ipISO==0 && nelem==0 )
00309 {
00310
00311 TeUsed[ipISO][nelem] = phycon.te*0.999;
00312
00313
00314
00315 }
00316
00317
00318 if( iso_ctrl.lgRandErrGen[ipISO] )
00319 {
00320
00321 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,iso_sp[ipISO][nelem].numLevels_max,IPRAD,0.0001f,0.0001f);
00322
00323
00324
00325
00326
00327
00328
00329
00330 }
00331
00332
00333 iso_sp[ipISO][nelem].RadRec_caseB = TotalRadRecomb - iso_sp[ipISO][nelem].fb[0].RadRecomb[ipRecRad];
00334 ASSERT( iso_sp[ipISO][nelem].RadRec_caseB > 0.);
00335
00336
00337 for( unsigned i = 0; i < iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
00338 {
00339 long index = iso_sp[ipISO][nelem].numLevels_max-1;
00340 opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] = iso_sp[ipISO][nelem].HighestLevelOpacStack[i];
00341 }
00342
00343
00344
00345
00346 if( !iso_sp[ipISO][nelem].lgLevelsLowered )
00347 {
00348
00349
00350
00351 TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
00352
00353 topoff = TotMinusExplicitResolved - RecCollapsed;
00354
00355
00356
00357
00358 if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) &&
00359 fabs( topoff/TotalRadRecomb ) > 0.01 )
00360 {
00361 fprintf(ioQQQ," PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",
00362 nelem, topoff/TotalRadRecomb, phycon.te );
00363 }
00364
00365
00366 if( !iso_ctrl.lgTopoff[ipISO] )
00367 topoff *= 1E-20;
00368
00369 topoff = MAX2( 0., topoff );
00370 double scale_factor = 1. + topoff/iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].RadRecomb[ipRecRad];
00371 ASSERT( scale_factor >= 1. );
00372
00373
00374 for( unsigned i = 0; i < iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
00375 {
00376 long index = iso_sp[ipISO][nelem].numLevels_max-1;
00377 opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] *= scale_factor;
00378 }
00379
00380
00381 iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].RadRecomb[ipRecRad] += topoff;
00382
00383
00384 if( Total_DR_Added > TotalRadRecomb/100. )
00385 {
00386 if( Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 )
00387 {
00388 fprintf(ioQQQ," PROBLEM negative DR topoff for %li, tot1, tot2 = %.2e\t%.2e rel error was %.2e, temp was %.2f\n",
00389 nelem,
00390 Total_DR_Added,
00391 ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO],
00392 Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0,
00393 phycon.te );
00394 }
00395 }
00396
00397 ASSERT( iso_sp[ipISO][nelem].numLevels_max == iso_sp[ipISO][nelem].numLevels_local );
00398
00399 if( iso_ctrl.lgDielRecom[ipISO] && iso_ctrl.lgTopoff[ipISO] )
00400 {
00401
00402
00403 iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].DielecRecomb +=
00404 MAX2( 0., ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added );
00405 }
00406 }
00407
00408 }
00409
00410
00411
00412
00413
00414
00415 iso_sp[ipISO][nelem].RadRec_effec = 0.;
00416
00417 for( ipLevel=0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
00418 {
00419
00420 if( opac.lgCaseB && ipLevel==0 )
00421 {
00422 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 1e-10;
00423 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 1e-10;
00424 }
00425 else if( cosmology.lgDo && ipLevel==0 )
00426 {
00427 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 1e-30;
00428 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 1e-30;
00429 }
00430 else
00431 {
00432 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] =
00433 RT_recom_effic(iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon);
00434
00435
00436 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] =
00437 MIN2( 1.,
00438 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] +
00439 (1.-iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc])*
00440 iso_sp[ipISO][nelem].fb[ipLevel].ConOpacRatio );
00441 }
00442
00443 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] >= 0. );
00444 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] >= 0. );
00445 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] <= 1. );
00446
00447
00448 iso_sp[ipISO][nelem].RadRec_effec += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad]*
00449 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc];
00450 }
00451
00452
00453 for( ipLevel=iso_sp[ipISO][nelem].numLevels_local; ipLevel<iso_sp[ipISO][nelem].numLevels_max; ++ipLevel )
00454 {
00455
00456 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 0.;
00457
00458 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 0.;
00459 }
00460
00461
00462 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
00463 {
00464 fprintf( ioQQQ, " iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
00465
00466
00467 fprintf( ioQQQ, " iso_radiative_recomb recomb effic" );
00468 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00469 {
00470 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] ));
00471 }
00472 fprintf( ioQQQ, "\n" );
00473
00474
00475 fprintf( ioQQQ, " iso_radiative_recomb recomb net effic" );
00476 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00477 {
00478 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc]) );
00479 }
00480 fprintf( ioQQQ, "\n" );
00481
00482
00483 fprintf( ioQQQ, " iso_radiative_recomb in optic dep" );
00484 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00485 {
00486 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[0][iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon-1] ));
00487 }
00488 fprintf( ioQQQ, "\n" );
00489
00490
00491 fprintf( ioQQQ, " iso_radiative_recomb out op depth" );
00492 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00493 {
00494 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[1][iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon-1] ));
00495 }
00496 fprintf( ioQQQ, "\n" );
00497
00498
00499 fprintf( ioQQQ, " iso_radiative_recomb rad rec coef " );
00500 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00501 {
00502 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad]) );
00503 }
00504 fprintf( ioQQQ, "\n" );
00505 }
00506
00507 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
00508 {
00509 fprintf( ioQQQ, " iso_radiative_recomb total rec coef" );
00510 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].RadRec_effec ));
00511 fprintf( ioQQQ, " case A=" );
00512 fprintf( ioQQQ,PrintEfmt("%10.3e",
00513 iso_sp[ipISO][nelem].RadRec_caseB + iso_sp[ipISO][nelem].fb[ipH1s].RadRecomb[ipRecRad] ) );
00514 fprintf( ioQQQ, " case B=");
00515 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].RadRec_caseB ));
00516 fprintf( ioQQQ, "\n" );
00517 }
00518
00519
00520
00521
00522 {
00523 bool lgOK = true;
00524 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
00525 {
00526 if( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] <= 0. )
00527 {
00528 fprintf( ioQQQ,
00529 " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
00530 ipISO, nelem, ipLevel, iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] , phycon.te);
00531 lgOK = false;
00532 }
00533 }
00534
00535 if( !lgOK )
00536 {
00537 ShowMe();
00538 cdEXIT(EXIT_FAILURE);
00539 }
00540
00541 }
00542
00543
00544 ASSERT( iso_sp[ipISO][nelem].fb[0].RadRecomb[ipRecRad] > 0. );
00545 ASSERT( iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_local-1].RadRecomb[ipRecRad] > 0. );
00546
00547
00548 if( save.lgioRecom )
00549 {
00550
00551 fprintf( save.ioRecom, "%s %s %2li ",
00552 iso_ctrl.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 );
00553 fprintf( save.ioRecom,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].RadRec_effec ));
00554 fprintf( save.ioRecom, "\n" );
00555 }
00556
00557 return;
00558 }
00559
00560 STATIC void iso_put_recomb_error( long ipISO, long nelem )
00561 {
00562 long level;
00563
00564
00565 realnum He_errorOpti[31] = {
00566 0.0000f,
00567 0.0009f,0.0037f,0.0003f,0.0003f,0.0003f,0.0018f,
00568 0.0009f,0.0050f,0.0007f,0.0003f,0.0001f,0.0007f,
00569 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f,
00570 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f };
00571
00572 realnum He_errorPess[31] = {
00573 0.0100f,
00574 0.0100f,0.0060f,0.0080f,0.0080f,0.0080f,0.0200f,
00575 0.0200f,0.0200f,0.0200f,0.0600f,0.0600f,0.0080f,
00576 0.0200f,0.0200f,0.0070f,0.0100f,0.0100f,0.0020f,0.0030f,0.0070f,
00577 0.0300f,0.0300f,0.0100f,0.0200f,0.0200f,0.0200f,0.0200f,0.0010f,0.0004f,0.0090f };
00578
00579
00580 for( long ipHi=0; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00581 {
00582 if( ipISO==ipHE_LIKE && nelem==ipHELIUM )
00583 {
00584 if( ipHi >= iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max )
00585 level = 21;
00586 else if( ipHi<=30 )
00587 level = ipHi;
00588 else
00589 level = iso_sp[ipISO][nelem].QuantumNumbers2Index[5][ MIN2(L_(ipHi),4) ][S_(ipHi)];
00590
00591 ASSERT( level <=30 );
00592
00593 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,ipHi,IPRAD, He_errorOpti[level], He_errorPess[level]);
00594 }
00595 else
00596 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,ipHi,IPRAD,0.1f,0.1f);
00597 }
00598
00599 return;
00600 }
00601
00602 void iso_radiative_recomb_effective( long ipISO, long nelem )
00603 {
00604 DEBUG_ENTRY( "iso_radiative_recomb_effective()" );
00605
00606
00607 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
00608 {
00609 iso_sp[ipISO][nelem].fb[ipHi].RadEffec = 0.;
00610
00611
00612 for( long ipHigher=ipHi; ipHigher < iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
00613 {
00614 ASSERT( iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] >= 0. );
00615 ASSERT( iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad] >= 0. );
00616
00617 iso_sp[ipISO][nelem].fb[ipHi].RadEffec += iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] *
00618 iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad];
00619 }
00620 }
00621
00622
00623
00624
00625 {
00626 enum {DEBUG_LOC=false};
00627
00628 if( DEBUG_LOC )
00629 {
00630 const int maxPrt=10;
00631
00632 fprintf( ioQQQ,"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem, phycon.te );
00633 fprintf( ioQQQ, "N\tL\tS\tRadEffec\tLifetime\n" );
00634
00635 for( long i=0; i<maxPrt; i++ )
00636 {
00637 fprintf( ioQQQ, "%li\t%li\t%li\t%e\t%e\n", N_(i), L_(i), S_(i),
00638 iso_sp[ipISO][nelem].fb[i].RadEffec,
00639 MAX2( 0., iso_sp[ipISO][nelem].st[i].lifetime() ) );
00640 }
00641 fprintf( ioQQQ, "\n" );
00642 }
00643 }
00644
00645
00646 if( iso_ctrl.lgRandErrGen[ipISO] )
00647 {
00648 dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
00649
00650
00651 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
00652 {
00653 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = 0.;
00654
00655
00656 for( long ipHigher=ipHi; ipHigher < iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
00657 {
00658 ASSERT( iso_sp[ipISO][nelem].ex[iso_sp[ipISO][nelem].numLevels_max][ipHigher].Error[IPRAD] >= 0. );
00659 ASSERT( iso_sp[ipISO][nelem].ex[ipHigher][ipHi].SigmaCascadeProb >= 0. );
00660
00661
00662 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec += pow( iso_sp[ipISO][nelem].ex[iso_sp[ipISO][nelem].numLevels_max][ipHigher].Error[IPRAD] *
00663 iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] * iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad], 2.) +
00664 pow( iso_sp[ipISO][nelem].ex[ipHigher][ipHi].SigmaCascadeProb * iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad], 2.);
00665 }
00666
00667 ASSERT( iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec >= 0. );
00668 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec );
00669
00670 for( long ipLo = 0; ipLo < ipHi; ipLo++ )
00671 {
00672 if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) )
00673 {
00674 double EnergyInRydbergs = iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd - iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd;
00675 realnum wavelength = (realnum)(RYDLAM/MAX2(1E-8,EnergyInRydbergs));
00676 double emissivity = iso_sp[ipISO][nelem].fb[ipHi].RadEffec * iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] * EN1RYD * EnergyInRydbergs;
00677 double sigma_emiss = 0., SigmaBranchRatio = 0.;
00678
00679 if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (N_(ipHi)<=5) )
00680 {
00681 SigmaBranchRatio = iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] * sqrt(
00682 pow( (double)iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD], 2. ) +
00683 pow( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) );
00684
00685 sigma_emiss = EN1RYD * EnergyInRydbergs * sqrt(
00686 pow( (double)iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec * iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo], 2.) +
00687 pow( SigmaBranchRatio * iso_sp[ipISO][nelem].fb[ipHi].RadEffec, 2.) );
00688
00689
00690 dprintf( ioQQQ, "%li\t%li\t", ipHi, ipLo );
00691 prt_wl( ioQQQ, wavelength );
00692 fprintf( ioQQQ, "\t%e\t%e\t%e\t%e\t%e\t%e\n",
00693 emissivity,
00694 sigma_emiss,
00695 iso_sp[ipISO][nelem].fb[ipHi].RadEffec,
00696 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec,
00697 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo],
00698 SigmaBranchRatio);
00699 }
00700 }
00701 }
00702 }
00703 }
00704
00705 return;
00706 }
00707
00708 double iso_RRCoef_Te( long ipISO, long nelem , long n )
00709 {
00710 double rate = 0.;
00711
00712 DEBUG_ENTRY( "iso_RRCoef_Te()" );
00713
00714 ASSERT( !iso_ctrl.lgNoRecombInterp[ipISO] );
00715
00716
00717 if( n == iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max )
00718 {
00719 rate = TempInterp( TeRRCoef, TotalRecomb[ipISO][nelem], N_ISO_TE_RECOMB );
00720 }
00721 else
00722 {
00723 rate = TempInterp( TeRRCoef, RRCoef[ipISO][nelem][n], N_ISO_TE_RECOMB );
00724 }
00725
00726
00727 rate = pow( 10. , rate );
00728
00729 return rate;
00730 }
00731
00732
00733
00734 double iso_recomb_check( long ipISO, long nelem, long level, double temperature )
00735 {
00736 double RecombRelError ,
00737 RecombInterp,
00738 RecombCalc,
00739 SaveTemp;
00740
00741 DEBUG_ENTRY( "iso_recomb_check()" );
00742
00743
00744 SaveTemp = phycon.te;
00745
00746
00747 TempChange(temperature);
00748
00749
00750
00751 RecombCalc = iso_radrecomb_from_cross_section( ipISO, temperature , nelem , level );
00752
00753
00754 RecombInterp = iso_RRCoef_Te( ipISO, nelem , level );
00755
00756
00757 TempChange(SaveTemp);
00758
00759 RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc );
00760
00761 return RecombRelError;
00762 }
00763
00764
00765 void iso_recomb_malloc( void )
00766 {
00767 DEBUG_ENTRY( "iso_recomb_malloc()" );
00768
00769 NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO );
00770 TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO );
00771 RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO );
00772
00773 for( long ipISO=0; ipISO<NISO; ipISO++ )
00774 {
00775 TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
00776 RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
00777
00778 NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM );
00779
00780 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
00781 {
00782 long int MaxLevels, maxN;
00783
00784 TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
00785
00786 if( nelem == ipISO )
00787 maxN = RREC_MAXN;
00788 else
00789 maxN = LIKE_RREC_MAXN( nelem );
00790
00791 NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 );
00792
00793 if( nelem == ipISO || dense.lgElmtOn[nelem] )
00794 {
00795
00796
00797 MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso_sp[ipISO][nelem].numLevels_max );
00798
00799
00800
00801
00802 RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) );
00803
00804 for( long ipLo=0; ipLo < MaxLevels;++ipLo )
00805 {
00806 RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
00807 }
00808 }
00809 }
00810 }
00811
00812 for(long i = 0; i < N_ISO_TE_RECOMB; i++)
00813 {
00814
00815 TeRRCoef[i] = 0.25*(i);
00816 }
00817
00818
00819
00820 TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f;
00821
00822 return;
00823 }
00824
00825 void iso_recomb_auxiliary_free( void )
00826 {
00827 DEBUG_ENTRY( "iso_recomb_auxiliary_free()" );
00828
00829 for( long i = 0; i< NISO; i++ )
00830 {
00831 free( NumLevRecomb[i] );
00832 }
00833 free( NumLevRecomb );
00834
00835 return;
00836 }
00837
00838 void iso_recomb_setup( long ipISO )
00839 {
00840 double RadRecombReturn;
00841 long int i, i1, i2, i3, i4, i5;
00842 long int ipLo, nelem;
00843
00844 const int chLine_LENGTH = 1000;
00845 char chLine[chLine_LENGTH];
00846
00847 const char* chFilename[NISO] =
00848 { "h_iso_recomb.dat", "he_iso_recomb.dat" };
00849
00850 FILE *ioDATA;
00851 bool lgEOL;
00852
00853 DEBUG_ENTRY( "iso_recomb_setup()" );
00854
00855
00856 if( iso_ctrl.lgCompileRecomb[ipISO] )
00857 {
00858 iso_ctrl.lgNoRecombInterp[ipISO] = false;
00859 }
00860
00861 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
00862 {
00863
00865
00866
00867 if( !iso_ctrl.lgCompileRecomb[ipISO] )
00868 {
00869 if( trace.lgTrace )
00870 fprintf( ioQQQ," iso_recomb_setup opening %s:", chFilename[ipISO] );
00871
00872
00873 try
00874 {
00875 ioDATA = open_data( chFilename[ipISO], "r" );
00876 }
00877 catch( cloudy_exit )
00878 {
00879 fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
00880 fprintf( ioQQQ, " but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
00881 ioDATA = NULL;
00882 }
00883 if( ioDATA == NULL )
00884 {
00885
00886 for( nelem = ipISO; nelem < LIMELM; nelem++ )
00887 {
00888 if( dense.lgElmtOn[nelem] )
00889 {
00890
00891 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00892 {
00893 TotalRecomb[ipISO][nelem][i] = 0.;
00894 }
00895
00896
00897
00898
00899 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
00900 {
00901
00902 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00903 {
00904
00905 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
00906 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
00907 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
00908 }
00909 }
00910 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00911 {
00912 for( i1 = iso_sp[ipISO][nelem].n_HighestResolved_max+1; i1< NHYDRO_MAX_LEVEL; i1++ )
00913 {
00914 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
00915 }
00916 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
00917 {
00918 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
00919 }
00920 TotalRecomb[ipISO][nelem][i] = log10( TotalRecomb[ipISO][nelem][i] );
00921 }
00922 }
00923 }
00924 }
00925
00926 else
00927 {
00928
00929 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00930 {
00931 fprintf( ioQQQ, " iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
00932 cdEXIT(EXIT_FAILURE);
00933 }
00934 i = 1;
00935 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00936 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00937 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00938 i4 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00939 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
00940 {
00941 fprintf( ioQQQ,
00942 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
00943 fprintf( ioQQQ,
00944 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
00945 RECOMBMAGIC ,
00946 NumLevRecomb[ipISO][ipISO],
00947 NumLevRecomb[ipISO][ipISO+1],
00948 N_ISO_TE_RECOMB,
00949 i1 , i2 , i3, i4 );
00950 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00951 fprintf( ioQQQ,
00952 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00953 cdEXIT(EXIT_FAILURE);
00954 }
00955
00956 i5 = 1;
00957
00958 for( nelem = ipISO; nelem < LIMELM; nelem++ )
00959 {
00960 for( ipLo=0; ipLo <= NumLevRecomb[ipISO][nelem]; ipLo++ )
00961 {
00962 i5++;
00963
00964 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00965 {
00966 fprintf( ioQQQ, " iso_recomb_setup could not read line %li of %s.\n", i5,
00967 chFilename[ipISO] );
00968 cdEXIT(EXIT_FAILURE);
00969 }
00970
00971 i3 = 1;
00972 i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
00973 i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
00974
00975 if( i1!=nelem || i2!=ipLo )
00976 {
00977 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
00978 fprintf( ioQQQ,
00979 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
00980 cdEXIT(EXIT_FAILURE);
00981 }
00982
00983
00984 for(i = 0; i < N_ISO_TE_RECOMB; i++)
00985 {
00986 double ThisCoef = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
00987
00988 if( nelem == ipISO || dense.lgElmtOn[nelem] )
00989 {
00990
00991 if( ipLo == NumLevRecomb[ipISO][nelem] )
00992 {
00993 TotalRecomb[ipISO][nelem][i] = ThisCoef;
00994 }
00995 else
00996 RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
00997 }
00998
00999 if( lgEOL )
01000 {
01001 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
01002 fprintf( ioQQQ,
01003 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
01004 cdEXIT(EXIT_FAILURE);
01005 }
01006 }
01007 }
01008
01009
01010
01011
01012 if( nelem == ipISO || dense.lgElmtOn[nelem] )
01013 {
01014 for( ipLo=NumLevRecomb[ipISO][nelem]; ipLo < iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
01015 {
01016 for(i = 0; i < N_ISO_TE_RECOMB; i++)
01017 {
01018
01019 RRCoef[ipISO][nelem][ipLo][i] = log10(iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo));
01020 }
01021 }
01022 }
01023 }
01024
01025
01026 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01027 {
01028 fprintf( ioQQQ, " iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
01029 cdEXIT(EXIT_FAILURE);
01030 }
01031 i = 1;
01032 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01033 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01034 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01035 i4 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01036
01037 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
01038 {
01039 fprintf( ioQQQ,
01040 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
01041 fprintf( ioQQQ,
01042 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
01043 RECOMBMAGIC ,
01044 NumLevRecomb[ipISO][ipISO],
01045 NumLevRecomb[ipISO][ipISO+1],
01046 N_ISO_TE_RECOMB,
01047 i1 , i2 , i3, i4 );
01048 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01049 fprintf( ioQQQ,
01050 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
01051 cdEXIT(EXIT_FAILURE);
01052 }
01053
01054
01055 fclose( ioDATA );
01056 }
01057 }
01058
01059 else if( iso_ctrl.lgCompileRecomb[ipISO] )
01060 {
01061
01062
01063 FILE *ioRECOMB;
01064
01065 ASSERT( !iso_ctrl.lgNoRecombInterp[ipISO] );
01066
01067
01068
01069
01070 ioRECOMB = open_data( chFilename[ipISO], "w", AS_LOCAL_ONLY );
01071 fprintf(ioRECOMB,"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient H-LIke [or HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
01072 RECOMBMAGIC ,
01073 NumLevRecomb[ipISO][ipISO],
01074 NumLevRecomb[ipISO][ipISO+1],
01075 N_ISO_TE_RECOMB,
01076 iso_ctrl.chISO[ipISO],
01077 NumLevRecomb[ipISO][ipISO],
01078 elementnames.chElementSym[ipISO],
01079 NumLevRecomb[ipISO][ipISO+1],
01080 N_ISO_TE_RECOMB );
01081
01082 for( nelem = ipISO; nelem < LIMELM; nelem++ )
01083 {
01084
01085 ASSERT( NumLevRecomb[ipISO][nelem] <= iso_sp[ipISO][nelem].numLevels_max );
01086
01087
01088 for(i = 0; i < N_ISO_TE_RECOMB; i++)
01089 {
01090 TotalRecomb[ipISO][nelem][i] = 0.;
01091 }
01092
01093 for( ipLo=ipHe1s1S; ipLo < NumLevRecomb[ipISO][nelem]; ipLo++ )
01094 {
01095 fprintf(ioRECOMB, "%li\t%li", nelem, ipLo );
01096
01097 for(i = 0; i < N_ISO_TE_RECOMB; i++)
01098 {
01099
01100 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
01101 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
01102 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
01103 fprintf(ioRECOMB, "\t%f", RRCoef[ipISO][nelem][ipLo][i] );
01104 }
01105 fprintf(ioRECOMB, "\n" );
01106 }
01107
01108
01109
01110
01111 fprintf(ioRECOMB, "%li\t%li", nelem, NumLevRecomb[ipISO][nelem] );
01112 for(i = 0; i < N_ISO_TE_RECOMB; i++)
01113 {
01114 for( i1 = ( (nelem == ipISO) ? (RREC_MAXN + 1) : (LIKE_RREC_MAXN( nelem ) + 1) ); i1< NHYDRO_MAX_LEVEL; i1++ )
01115 {
01116 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
01117 }
01118 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
01119 {
01120 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
01121 }
01122 fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[ipISO][nelem][i] ) );
01123 }
01124 fprintf(ioRECOMB, "\n" );
01125 }
01126
01127 fprintf(ioRECOMB,"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient [H-LIke/HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
01128 RECOMBMAGIC ,
01129 NumLevRecomb[ipISO][ipISO],
01130 NumLevRecomb[ipISO][ipISO+1],
01131 N_ISO_TE_RECOMB,
01132 iso_ctrl.chISO[ipISO],
01133 NumLevRecomb[ipISO][ipISO],
01134 elementnames.chElementSym[ipISO],
01135 NumLevRecomb[ipISO][ipISO+1],
01136 N_ISO_TE_RECOMB );
01137
01138 fclose( ioRECOMB );
01139
01140 fprintf( ioQQQ, "iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
01141 fprintf( ioQQQ, "The compilation is completed successfully.\n");
01142 cdEXIT(EXIT_SUCCESS);
01143 }
01144 }
01145
01146 return;
01147 }
01148
01149 double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo )
01150 {
01151 double rate;
01152 long ipTe, i;
01153 double TeDRCoef[NUM_DR_TEMPS];
01154 const freeBound *fb = &iso_sp[ipISO][nelem].fb[ipLo];
01155 const double Te_over_Z1_Squared[NUM_DR_TEMPS] = {
01156 1.00000, 1.30103, 1.69897, 2.00000, 2.30103, 2.69897, 3.00000,
01157 3.30103, 3.69897, 4.00000, 4.30103, 4.69897, 5.00000, 5.30103,
01158 5.69897, 6.00000, 6.30103, 6.69897, 7.00000 };
01159
01160 DEBUG_ENTRY( "iso_dielec_recomb_rate()" );
01161
01162
01163 ASSERT( ipISO == ipHE_LIKE );
01164 ASSERT( ipLo >= 0 );
01165
01166
01167 for( i=0; i<NUM_DR_TEMPS; i++ )
01168 {
01169 TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (double) nelem );
01170 }
01171
01172 if( ipLo == ipHe1s1S )
01173 {
01174 rate = 0.;
01175 }
01176 else if( ipLo<iso_sp[ipISO][nelem].numLevels_max )
01177 {
01178 if( phycon.alogte <= TeDRCoef[0] )
01179 {
01180
01181 rate = fb->DielecRecombVsTemp[0];
01182 }
01183 else if( phycon.alogte >= TeDRCoef[NUM_DR_TEMPS-1] )
01184 {
01185
01186 rate = fb->DielecRecombVsTemp[NUM_DR_TEMPS-1] *
01187 pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] - phycon.alogte ) ) ;
01188 }
01189 else
01190 {
01191
01192 ipTe = hunt_bisect( TeDRCoef, NUM_DR_TEMPS, phycon.alogte );
01193
01194 ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1) );
01195
01196 if( fb->DielecRecombVsTemp[ipTe+1] == 0. )
01197 rate = 0.;
01198 else if( fb->DielecRecombVsTemp[ipTe] == 0. )
01199 rate = fb->DielecRecombVsTemp[ipTe+1];
01200 else
01201 {
01202
01203 rate = log10( fb->DielecRecombVsTemp[ipTe]) +
01204 (phycon.alogte-TeDRCoef[ipTe])*
01205 (log10(fb->DielecRecombVsTemp[ipTe+1])-log10(fb->DielecRecombVsTemp[ipTe]))/
01206 (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]);
01207
01208 rate = pow( 10., rate );
01209 }
01210 }
01211 }
01212 else
01213 {
01214 rate = 0.;
01215 }
01216
01217 ASSERT( rate >= 0. && rate < 1.0e-12 );
01218
01219 return rate*iso_ctrl.lgDielRecom[ipISO];
01220 }
01221
01222
01224 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements )
01225 {
01226 static long int ipTe=-1;
01227 double rate = 0.;
01228 long i0;
01229
01230 DEBUG_ENTRY( "TempInterp()" );
01231
01232 ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 );
01233
01234 if( ipTe < 0 )
01235 {
01236
01237 if( ( phycon.alogte < TempArray[0] ) ||
01238 ( phycon.alogte > TempArray[NumElements-1] ) )
01239 {
01240 fprintf(ioQQQ," TempInterp called with te out of bounds \n");
01241 cdEXIT(EXIT_FAILURE);
01242 }
01243 ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );
01244 }
01245 else if( phycon.alogte < TempArray[ipTe] )
01246 {
01247
01248 ASSERT( phycon.alogte > TempArray[0] );
01249
01250 while( ( phycon.alogte < TempArray[ipTe] ) && ipTe > 0)
01251 --ipTe;
01252 }
01253 else if( phycon.alogte > TempArray[ipTe+1] )
01254 {
01255
01256 ASSERT( phycon.alogte <= TempArray[NumElements-1] );
01257
01258 while( ( phycon.alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1)
01259 ++ipTe;
01260 }
01261
01262 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
01263
01264
01265 ASSERT( ( phycon.alogte >= TempArray[ipTe] )
01266 && ( phycon.alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );
01267
01268 if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. )
01269 {
01270 rate = 0.;
01271 }
01272 else
01273 {
01274
01275 const int ORDER = 3;
01276 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
01277 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
01278 }
01279
01280 return rate;
01281 }