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