cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*iso_radiative_recomb find state specific creation and destruction rates for iso sequences */ 00004 /*iso_RRCoef_Te - evaluate radiative recombination coef at some temperature */ 00005 /*iso_recomb_check - called by SanityCheck to confirm that recombination coef are ok, 00006 * return value is relative error between new calculation of recom, and interp value */ 00007 00008 #include "cddefines.h" 00009 #include "atmdat.h" 00010 #include "conv.h" 00011 #include "dense.h" 00012 #include "elementnames.h" 00013 #include "helike.h" 00014 #include "helike_recom.h" 00015 #include "hydrogenic.h" 00016 #include "ionbal.h" 00017 #include "iso.h" 00018 #include "opacity.h" 00019 #include "phycon.h" 00020 #include "physconst.h" 00021 #include "prt.h" 00022 #include "punch.h" 00023 #include "thermal.h" 00024 #include "thirdparty.h" 00025 #include "trace.h" 00026 #include "rt.h" 00027 00028 /* this will save log of radiative recombination rate coefficients at N_ISO_TE_RECOMB temperatures. 00029 * there will be NumLevRecomb[ipISO][nelem] levels saved in RRCoef[nelem][level][temp] */ 00030 static double ****RRCoef/*[LIMELM][NumLevRecomb[ipISO][nelem]][N_ISO_TE_RECOMB]*/; 00031 static long **NumLevRecomb; 00032 static double ***TotalRecomb; /*[ipISO][nelem][i]*/ 00033 00034 /* the array of logs of temperatures at which RRCoef was defined */ 00035 static double TeRRCoef[N_ISO_TE_RECOMB]; 00036 00037 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements ); 00038 00039 double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo) 00040 { 00041 double alpha; 00042 00043 DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" ); 00044 00045 if( ipISO == ipH_LIKE ) 00046 alpha = hlike_radrecomb_from_cross_section( temp, nelem, ipLo); 00047 else if( ipISO == ipHE_LIKE ) 00048 alpha = helike_radrecomb_from_cross_section( temp, nelem, ipLo); 00049 else 00050 TotalInsanity(); 00051 00052 return alpha; 00053 } 00054 00055 /*=======================================================*/ 00056 /* iso_radiative_recomb get rad recomb rate coefficients for iso sequences */ 00057 void iso_radiative_recomb( 00058 long ipISO, 00059 /* nelem on the c scale, He is 1 */ 00060 long nelem ) 00061 { 00062 long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel; 00063 double topoff, TotMinusExplicitResolved, 00064 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.; 00065 double RecExplictLevels, TotalRadRecomb, RecCollapsed; 00066 static double TeUsed[NISO][LIMELM]={ 00067 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00068 0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00069 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, 00070 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00071 0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 00072 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}}; 00073 00074 DEBUG_ENTRY( "iso_radiative_recomb()" ); 00075 00076 /* evaluate recombination escape probability for all levels */ 00077 00078 /* define radiative recombination rates for all levels */ 00079 /* this will be the sum of all levels explicitly included in the model atom */ 00080 RecExplictLevels = 0.; 00081 00082 /* number of resolved levels, so first collapsed level is [ipFirstCollapsed] */ 00083 ipFirstCollapsed = iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem]; 00084 00085 if( (fabs(1.-TeUsed[ipISO][nelem]/phycon.te)> 0.001) || hydro.lgReevalRecom || !conv.nTotalIoniz ) 00086 { 00087 TeUsed[ipISO][nelem] = phycon.te; 00088 00089 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel ) 00090 { 00091 /* this is radiative recombination rate coefficient */ 00092 double RadRec; 00093 00094 if( !iso.lgNoRecombInterp[ipISO] ) 00095 { 00096 RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel ); 00097 } 00098 else 00099 { 00100 RadRec = iso_radrecomb_from_cross_section( ipISO, phycon.te, nelem, ipLevel); 00101 } 00102 00103 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec; 00104 00105 if( iso.lgRandErrGen[ipISO] ) 00106 { 00107 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 00108 /* this has to be from iso.numLevels_max instead of iso.numLevels_local because 00109 * the error factors for rrc are always stored at iso.numLevels_max, regardless of 00110 * continuum lowering effects. */ 00111 iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD]; 00112 } 00113 00114 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. ); 00115 00116 RecExplictLevels += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 00117 00118 if( iso.lgDielRecom[ipISO] ) 00119 { 00120 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */ 00121 iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel ); 00122 Total_DR_Added += iso.DielecRecomb[ipISO][nelem][ipLevel]; 00123 } 00124 } 00125 00126 /**************************************************/ 00127 /*** Add on recombination to collapsed levels ***/ 00128 /**************************************************/ 00129 RecCollapsed = 0.; 00130 for( ipLevel=ipFirstCollapsed; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel ) 00131 { 00132 /* use hydrogenic for collapsed levels */ 00133 double RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, N_(ipLevel), phycon.te); 00134 00135 /* this is radiative recombination rate coefficient */ 00136 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec; 00137 00138 /* This kills recombination into the collapsed level so that the forced 00139 * statistical weighting can be bypassed */ 00140 if( !iso.lgTopoff[ipISO] ) 00141 { 00142 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 1E-10; 00143 } 00144 00145 if( iso.lgRandErrGen[ipISO] ) 00146 { 00147 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 00148 /* this has to be from iso.numLevels_max instead of iso.numLevels_local because 00149 * the error factors for rrc are always stored at iso.numLevels_max, regardless of 00150 * continuum lowering effects. */ 00151 iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD]; 00152 } 00153 00154 RecCollapsed += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 00155 00156 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. ); 00157 00158 if( iso.lgDielRecom[ipISO] ) 00159 { 00160 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */ 00161 iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel ); 00162 Total_DR_Added += iso.DielecRecomb[ipISO][nelem][ipLevel]; 00163 } 00164 } 00165 00166 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 00167 for( ipLevel = 0; ipLevel<iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00168 { 00169 if( N_(ipLevel) == ThisN ) 00170 { 00171 TotRRThisN += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 00172 } 00173 else 00174 { 00175 ASSERT( (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) ); 00176 LastN = ThisN; 00177 ThisN = N_(ipLevel); 00178 TotRRLastN = TotRRThisN; 00179 TotRRThisN = iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]; 00180 00181 { 00182 /* Print the sum of recombination coefficients per n at current temp. */ 00183 /*@-redef@*/ 00184 enum {DEBUG_LOC=false}; 00185 /*@+redef@*/ 00186 static long RUNONCE = false; 00187 00188 if( !RUNONCE && DEBUG_LOC ) 00189 { 00190 static long FIRSTTIME = true; 00191 00192 if( FIRSTTIME ) 00193 { 00194 fprintf( ioQQQ,"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n", 00195 ipISO, nelem, phycon.te); 00196 FIRSTTIME= false; 00197 } 00198 00199 fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN ); 00200 } 00201 RUNONCE = true; 00202 } 00203 } 00204 } 00205 00206 /* Get total recombination into all levels, including those not explicitly considered. */ 00207 if( iso.lgNoRecombInterp[ipISO] ) 00208 { 00209 /* We are not interpolating, must calculate total now, as sum of resolved and collapsed levels... */ 00210 TotalRadRecomb = RecCollapsed + RecExplictLevels; 00211 00212 /* Plus additional levels out to a predefined limit... */ 00213 for( long nLo = N_(ipFirstCollapsed) + iso.nCollapsed_max[ipISO][nelem]; nLo < NHYDRO_MAX_LEVEL; nLo++ ) 00214 { 00215 TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te); 00216 } 00217 /* Plus a bunch more for good measure */ 00218 for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ ) 00219 { 00220 TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo ); 00221 } 00222 } 00223 else 00224 { 00225 /* We are interpolating, and total was calculated here in iso_recomb_setup */ 00226 TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ); 00227 } 00228 00229 /* If generating random error, apply one to total recombination */ 00230 if( iso.lgRandErrGen[ipISO] ) 00231 { 00232 /* this has to be from iso.numLevels_max instead of iso.numLevels_local because 00233 * the error factors for rrc are always stored at iso.numLevels_max, regardless of 00234 * continuum lowering effects. */ 00235 TotalRadRecomb *= 00236 iso.ErrorFactor[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][iso.numLevels_max[ipISO][nelem]][IPRAD]; 00237 } 00238 00239 /* this is case B recombination, sum without the ground included */ 00240 iso.RadRec_caseB[ipISO][nelem] = TotalRadRecomb - iso.RadRecomb[ipISO][nelem][0][ipRecRad]; 00241 ASSERT( iso.RadRec_caseB[ipISO][nelem] > 0.); 00242 00243 /**********************************************************************/ 00244 /*** Add topoff (excess) recombination to top level. This is only ***/ 00245 /*** done if atom is not full size. ***/ 00246 /**********************************************************************/ 00247 if( !iso.lgLevelsLowered[ipISO][nelem] ) 00248 { 00249 /* at this point we have RecExplictLevels, the sum of radiative recombination 00250 * to all levels explicitly included in the model atom the total 00251 * recombination rate. The difference is the amount of "topoff" we will need to do */ 00252 TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels; 00253 00254 topoff = TotMinusExplicitResolved - RecCollapsed; 00255 00256 /* the t_ADfA::Inst().rad_rec fits are too small at high temperatures, so this atom is 00257 * better than the topoff. Only a problem for helium itself, at high temperatures. 00258 * complain if the negative topoff is not for this case */ 00259 if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) && 00260 fabs( topoff/TotalRadRecomb ) > 0.01 ) 00261 { 00262 fprintf(ioQQQ," PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n", 00263 nelem, topoff/TotalRadRecomb, phycon.te ); 00264 } 00265 00266 if( !iso.lgTopoff[ipISO] ) 00267 topoff *= 1E-20; 00268 00269 /* We always have at least one collapsed level if continuum is not lowered. Put topoff there. */ 00270 iso.RadRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1][ipRecRad] += MAX2(0., topoff ); 00271 00272 /* check for negative DR topoff, but only if Total_DR_Added is not negligible compared with TotalRadRecomb */ 00273 if( Total_DR_Added > TotalRadRecomb/100. ) 00274 { 00275 if( Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 ) 00276 { 00277 fprintf(ioQQQ," PROBLEM negative DR topoff for %li, rel error was %.2e, temp was %.2f\n", 00278 nelem, 00279 Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0, 00280 phycon.te ); 00281 } 00282 } 00283 00284 if( iso.lgDielRecom[ipISO] ) 00285 { 00286 /* \todo 2 suppress this total rate for continuum lowering using factors from Jordan (1969). */ 00287 /* put extra DR in top level */ 00288 iso.DielecRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] = 00289 ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added; 00290 } 00291 } 00292 00293 } 00294 00295 /**************************************************************/ 00296 /*** Stuff escape probabilities, and effective rad recomb ***/ 00297 /**************************************************************/ 00298 00299 /* total effective radiative recombination, initialize to zero */ 00300 iso.RadRec_effec[ipISO][nelem] = 0.; 00301 00302 for( ipLevel=0; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel ) 00303 { 00304 /* option for case b conditions, kill ground state recombination */ 00305 if( opac.lgCaseB && ipLevel==0 ) 00306 { 00307 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 1e-10; 00308 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 1e-10; 00309 } 00310 else 00311 { 00312 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 00313 RT_recom_effic(iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]); 00314 00315 /* net escape prob includes dest by background opacity */ 00316 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 00317 MIN2( 1., 00318 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] + 00319 (1.-iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc])* 00320 iso.ConOpacRatio[ipISO][nelem][ipLevel] ); 00321 } 00322 00323 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] >= 0. ); 00324 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] >= 0. ); 00325 ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] <= 1. ); 00326 00327 /* sum of all effective rad rec */ 00328 iso.RadRec_effec[ipISO][nelem] += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]* 00329 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]; 00330 } 00331 00332 /* zero out escape probabilities of levels above numLevels_local */ 00333 for( ipLevel=iso.numLevels_local[ipISO][nelem]; ipLevel<iso.numLevels_max[ipISO][nelem]; ++ipLevel ) 00334 { 00335 /* this is escape probability */ 00336 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 0.; 00337 /* net escape prob includes dest by background opacity */ 00338 iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 0.; 00339 } 00340 00341 /* trace escape probabilities */ 00342 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) ) 00343 { 00344 fprintf( ioQQQ, " iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem ); 00345 00346 /* print continuum escape probability */ 00347 fprintf( ioQQQ, " iso_radiative_recomb recomb effic" ); 00348 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00349 { 00350 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] )); 00351 } 00352 fprintf( ioQQQ, "\n" ); 00353 00354 /* net recombination efficiency factor, including background opacity*/ 00355 fprintf( ioQQQ, " iso_radiative_recomb recomb net effic" ); 00356 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00357 { 00358 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]) ); 00359 } 00360 fprintf( ioQQQ, "\n" ); 00361 00362 /* inward continuous optical depths */ 00363 fprintf( ioQQQ, " iso_radiative_recomb in optic dep" ); 00364 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00365 { 00366 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 00367 } 00368 fprintf( ioQQQ, "\n" ); 00369 00370 /* outward continuous optical depths*/ 00371 fprintf( ioQQQ, " iso_radiative_recomb out op depth" ); 00372 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00373 { 00374 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] )); 00375 } 00376 fprintf( ioQQQ, "\n" ); 00377 00378 /* print radiative recombination coefficients */ 00379 fprintf( ioQQQ, " iso_radiative_recomb rad rec coef " ); 00380 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00381 { 00382 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]) ); 00383 } 00384 fprintf( ioQQQ, "\n" ); 00385 } 00386 00387 if( trace.lgTrace ) 00388 { 00389 fprintf( ioQQQ, " iso_radiative_recomb total rec coef" ); 00390 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipISO][nelem] )); 00391 fprintf( ioQQQ, " case A=" ); 00392 fprintf( ioQQQ,PrintEfmt("%10.3e", 00393 iso.RadRec_caseB[ipISO][nelem] + iso.RadRecomb[ipISO][nelem][ipH1s][ipRecRad] ) ); 00394 fprintf( ioQQQ, " case B="); 00395 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipISO][nelem] )); 00396 fprintf( ioQQQ, "\n" ); 00397 } 00398 00399 /****************************/ 00400 /*** begin sanity check ***/ 00401 /****************************/ 00402 { 00403 bool lgOK = true; 00404 for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ ) 00405 { 00406 if( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] <= 0. ) 00407 { 00408 fprintf( ioQQQ, 00409 " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n", 00410 ipISO, nelem, ipLevel, iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] , phycon.te); 00411 lgOK = false; 00412 } 00413 } 00414 /* bail if we found problems */ 00415 if( !lgOK ) 00416 { 00417 ShowMe(); 00418 cdEXIT(EXIT_FAILURE); 00419 } 00420 /*end sanity check */ 00421 } 00422 00423 /* confirm that we have good rec coef at bottom and top of atom/ion */ 00424 ASSERT( iso.RadRecomb[ipISO][nelem][0][ipRecRad] > 0. ); 00425 ASSERT( iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecRad] > 0. ); 00426 00427 /* set true when punch recombination coefficients command entered */ 00428 if( punch.lgioRecom ) 00429 { 00430 /* this prints Z on physical, not C, scale */ 00431 fprintf( punch.ioRecom, "%s %s %2li ", 00432 iso.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 ); 00433 fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipISO][nelem] )); 00434 fprintf( punch.ioRecom, "\n" ); 00435 } 00436 00437 return; 00438 } 00439 00440 void iso_radiative_recomb_effective( long ipISO, long nelem ) 00441 { 00442 DEBUG_ENTRY( "iso_radiative_recomb_effective()" ); 00443 00444 /* Find effective recombination coefficients */ 00445 for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00446 { 00447 iso.RadEffec[ipISO][nelem][ipHi] = 0.; 00448 00449 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 00450 for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ ) 00451 { 00452 ASSERT( iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. ); 00453 ASSERT( iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad] >= 0. ); 00454 00455 iso.RadEffec[ipISO][nelem][ipHi] += iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] * 00456 iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad]; 00457 } 00458 } 00459 00460 /**************************************************************/ 00461 /*** option to print effective recombination coefficients ***/ 00462 /**************************************************************/ 00463 { 00464 enum {DEBUG_LOC=false}; 00465 00466 if( DEBUG_LOC ) 00467 { 00468 const int maxPrt=10; 00469 00470 fprintf( ioQQQ,"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem, phycon.te ); 00471 fprintf( ioQQQ, "N\tL\tS\tRadEffec\tLifetime\n" ); 00472 00473 for( long i=0; i<maxPrt; i++ ) 00474 { 00475 fprintf( ioQQQ, "%li\t%li\t%li\t%e\t%e\n", N_(i), L_(i), S_(i), 00476 iso.RadEffec[ipISO][nelem][i], 00477 MAX2( 0., StatesElem[ipISO][nelem][i].lifetime ) ); 00478 } 00479 fprintf( ioQQQ, "\n" ); 00480 } 00481 } 00482 00483 /* If we have the variable set, find errors in rad rates */ 00484 if( iso.lgRandErrGen[ipISO] ) 00485 { 00486 dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" ); 00487 00488 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 00489 for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00490 { 00491 iso.SigmaRadEffec[ipISO][nelem][ipHi] = 0.; 00492 00493 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */ 00494 for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ ) 00495 { 00496 ASSERT( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] >= 0. ); 00497 ASSERT( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. ); 00498 00499 /* iso.RadRecomb has to appear here because iso.Error is only relative error */ 00500 iso.SigmaRadEffec[ipISO][nelem][ipHi] += pow( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] * 00501 iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.) + 00502 pow( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.); 00503 } 00504 00505 ASSERT( iso.SigmaRadEffec[ipISO][nelem][ipHi] >= 0. ); 00506 iso.SigmaRadEffec[ipISO][nelem][ipHi] = sqrt( iso.SigmaRadEffec[ipISO][nelem][ipHi] ); 00507 00508 for( long ipLo = 0; ipLo < ipHi; ipLo++ ) 00509 { 00510 if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) ) 00511 { 00512 double EnergyInRydbergs = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]; 00513 realnum wavelength = (realnum)(RYDLAM/MAX2(1E-8,EnergyInRydbergs)); 00514 double emissivity = iso.RadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * EN1RYD * EnergyInRydbergs; 00515 double sigma_emiss = 0., SigmaBranchRatio = 0.; 00516 00517 if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (N_(ipHi)<=5) ) 00518 { 00519 SigmaBranchRatio = iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * sqrt( 00520 pow( (double)iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD], 2. ) + 00521 pow( iso.SigmaAtot[ipISO][nelem][ipHi]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) ); 00522 00523 sigma_emiss = EN1RYD * EnergyInRydbergs * sqrt( 00524 pow( (double)iso.SigmaRadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo], 2.) + 00525 pow( SigmaBranchRatio * iso.RadEffec[ipISO][nelem][ipHi], 2.) ); 00526 00527 /* \todo 2 make this a trace option */ 00528 dprintf( ioQQQ, "%li\t%li\t", ipHi, ipLo ); 00529 prt_wl( ioQQQ, wavelength ); 00530 fprintf( ioQQQ, "\t%e\t%e\t%e\t%e\t%e\t%e\n", 00531 emissivity, 00532 sigma_emiss, 00533 iso.RadEffec[ipISO][nelem][ipHi], 00534 iso.SigmaRadEffec[ipISO][nelem][ipHi], 00535 iso.BranchRatio[ipISO][nelem][ipHi][ipLo], 00536 SigmaBranchRatio); 00537 } 00538 } 00539 } 00540 } 00541 } 00542 00543 return; 00544 } 00545 /*iso_RRCoef_Te evaluated radiative recombination coef at some temperature */ 00546 double iso_RRCoef_Te( long ipISO, long nelem , long n ) 00547 { 00548 double rate = 0.; 00549 00550 DEBUG_ENTRY( "iso_RRCoef_Te()" ); 00551 00552 ASSERT( !iso.lgNoRecombInterp[ipISO] ); 00553 00554 /* if n is equal to the number of levels, return the total recomb, else recomb for given level. */ 00555 if( n == iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ) 00556 { 00557 rate = TempInterp( TeRRCoef, TotalRecomb[ipISO][nelem], N_ISO_TE_RECOMB ); 00558 } 00559 else 00560 { 00561 rate = TempInterp( TeRRCoef, RRCoef[ipISO][nelem][n], N_ISO_TE_RECOMB ); 00562 } 00563 00564 /* that was the log, now make linear */ 00565 rate = pow( 10. , rate ); 00566 00567 return rate; 00568 } 00569 00570 /*iso_recomb_check called by SanityCheck to confirm that recombination coef are ok 00571 * return value is relative error between new calculation of recom, and interp value */ 00572 double iso_recomb_check( long ipISO, long nelem, long level, double temperature ) 00573 { 00574 double RecombRelError , 00575 RecombInterp, 00576 RecombCalc, 00577 SaveTemp; 00578 00579 DEBUG_ENTRY( "iso_recomb_check()" ); 00580 00581 /* first set temp to desired value */ 00582 SaveTemp = phycon.te; 00583 TempChange(temperature); 00584 00585 /* actually calculate the recombination coefficient from the Milne relation, 00586 * normally only done due to compile he-like command */ 00587 RecombCalc = iso_radrecomb_from_cross_section( ipISO, temperature , nelem , level ); 00588 00589 /* interpolate the recombination coefficient, this is the usual method */ 00590 RecombInterp = iso_RRCoef_Te( ipISO, nelem , level ); 00591 00592 /* reset temp */ 00593 TempChange(SaveTemp); 00594 00595 RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc ); 00596 00597 return RecombRelError; 00598 } 00599 00600 /* malloc space needed for iso recombination tables */ 00601 void iso_recomb_malloc( void ) 00602 { 00603 DEBUG_ENTRY( "iso_recomb_malloc()" ); 00604 00605 NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO ); 00606 TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO ); 00607 RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO ); 00608 00609 for( long ipISO=0; ipISO<NISO; ipISO++ ) 00610 { 00611 TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00612 RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ); 00613 /* The number of recombination coefficients to be read from file for each element. */ 00614 NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ); 00615 00616 for( long nelem=ipISO; nelem < LIMELM; ++nelem ) 00617 { 00618 long int MaxLevels, maxN; 00619 00620 TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB ); 00621 00622 if( nelem == ipISO ) 00623 maxN = RREC_MAXN; 00624 else 00625 maxN = LIKE_RREC_MAXN( nelem ); 00626 00627 NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 ); 00628 00629 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 00630 { 00631 /* must always have at least NumLevRecomb[ipISO][nelem] levels since that is number 00632 * that will be read in from he rec data file, but possible to require more */ 00633 MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso.numLevels_max[ipISO][nelem] ); 00634 00635 /* always define this */ 00636 /* >>chng 02 jan 24, RRCoef will be iso.numLevels_max[ipISO][nelem] levels, not iso.numLevels_max, 00637 * code will stop if more than this is requested */ 00638 RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) ); 00639 00640 for( long ipLo=0; ipLo < MaxLevels;++ipLo ) 00641 { 00642 RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB ); 00643 } 00644 } 00645 } 00646 } 00647 00648 for(long i = 0; i < N_ISO_TE_RECOMB; i++) 00649 { 00650 /* this is the vector of temperatures */ 00651 TeRRCoef[i] = 0.25*(i); 00652 } 00653 00654 /* >>chng 06 jun 06, NP, assert thrown at T == 1e10 K, just bump the 00655 * high temperature end slightly. */ 00656 TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f; 00657 00658 return; 00659 } 00660 00661 void iso_recomb_auxiliary_free( void ) 00662 { 00663 DEBUG_ENTRY( "iso_recomb_auxiliary_free()" ); 00664 00665 for( long i = 0; i< NISO; i++ ) 00666 { 00667 free( NumLevRecomb[i] ); 00668 } 00669 free( NumLevRecomb ); 00670 00671 return; 00672 } 00673 00674 void iso_recomb_setup( long ipISO ) 00675 { 00676 double RadRecombReturn; 00677 long int i, i1, i2, i3, i4, i5; 00678 long int ipLo, nelem; 00679 00680 const int chLine_LENGTH = 1000; 00681 char chLine[chLine_LENGTH]; 00682 /* this must be longer than data path string, set in path.h/cpu.cpp */ 00683 const char* chFilename[NISO] = 00684 { "h_iso_recomb.dat", "he_iso_recomb.dat" }; 00685 00686 FILE *ioDATA; 00687 bool lgEOL; 00688 00689 DEBUG_ENTRY( "iso_recomb_setup()" ); 00690 00691 /* if we are compiling the recombination data file, we must interpolate in temperature */ 00692 if( iso.lgCompileRecomb[ipISO] ) 00693 { 00694 iso.lgNoRecombInterp[ipISO] = false; 00695 } 00696 00697 if( !iso.lgNoRecombInterp[ipISO] ) 00698 { 00699 /******************************************************************/ 00701 /******************************************************************/ 00702 /* This flag says we are not compiling the data file */ 00703 if( !iso.lgCompileRecomb[ipISO] ) 00704 { 00705 if( trace.lgTrace ) 00706 fprintf( ioQQQ," iso_recomb_setup opening %s:", chFilename[ipISO] ); 00707 00708 /* Now try to read from file...*/ 00709 try 00710 { 00711 ioDATA = open_data( chFilename[ipISO], "r" ); 00712 } 00713 catch( cloudy_exit ) 00714 { 00715 fprintf( ioQQQ, " Defaulting to on-the-fly computation, "); 00716 fprintf( ioQQQ, " but the code runs much faster if you compile %s!\n", chFilename[ipISO]); 00717 ioDATA = NULL; 00718 } 00719 if( ioDATA == NULL ) 00720 { 00721 /* Do on the fly computation of R.R. Coef's instead. */ 00722 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 00723 { 00724 if( dense.lgElmtOn[nelem] ) 00725 { 00726 /* Zero out the recombination sum array. */ 00727 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00728 { 00729 TotalRecomb[ipISO][nelem][i] = 0.; 00730 } 00731 00732 /* NumLevRecomb[ipISO][nelem] corresponds to n = 40 for H and He and 20 for ions, at present */ 00733 /* There is no need to fill in values for collapsed levels, because we do not need to 00734 * interpolate for a given temperature, just calculate it directly with a hydrogenic routine. */ 00735 for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ ) 00736 { 00737 /* loop over temperatures to produce array of recombination coefficients */ 00738 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00739 { 00740 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */ 00741 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo); 00742 TotalRecomb[ipISO][nelem][i] += RadRecombReturn; 00743 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn); 00744 } 00745 } 00746 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00747 { 00748 for( i1 = iso.n_HighestResolved_max[ipISO][nelem]+1; i1< NHYDRO_MAX_LEVEL; i1++ ) 00749 { 00750 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i])); 00751 } 00752 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ ) 00753 { 00754 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 ); 00755 } 00756 TotalRecomb[ipISO][nelem][i] = log10( TotalRecomb[ipISO][nelem][i] ); 00757 } 00758 } 00759 } 00760 } 00761 /* Data file is present and readable...begin read. */ 00762 else 00763 { 00764 /* check that magic number is ok */ 00765 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00766 { 00767 fprintf( ioQQQ, " iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]); 00768 cdEXIT(EXIT_FAILURE); 00769 } 00770 i = 1; 00771 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00772 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00773 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00774 i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00775 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB ) 00776 { 00777 fprintf( ioQQQ, 00778 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] ); 00779 fprintf( ioQQQ, 00780 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" , 00781 RECOMBMAGIC , 00782 NumLevRecomb[ipISO][ipISO], 00783 NumLevRecomb[ipISO][ipISO+1], 00784 N_ISO_TE_RECOMB, 00785 i1 , i2 , i3, i4 ); 00786 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00787 fprintf( ioQQQ, 00788 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" ); 00789 cdEXIT(EXIT_FAILURE); 00790 } 00791 00792 i5 = 1; 00793 /* now read in the data */ 00794 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 00795 { 00796 for( ipLo=0; ipLo <= NumLevRecomb[ipISO][nelem]; ipLo++ ) 00797 { 00798 i5++; 00799 /* get next line image */ 00800 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00801 { 00802 fprintf( ioQQQ, " iso_recomb_setup could not read line %li of %s.\n", i5, 00803 chFilename[ipISO] ); 00804 cdEXIT(EXIT_FAILURE); 00805 } 00806 /* each line starts with element and level number */ 00807 i3 = 1; 00808 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 00809 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL); 00810 /* check that these numbers are correct */ 00811 if( i1!=nelem || i2!=ipLo ) 00812 { 00813 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]); 00814 fprintf( ioQQQ, 00815 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" ); 00816 cdEXIT(EXIT_FAILURE); 00817 } 00818 00819 /* loop over temperatures to produce array of recombination coefficients */ 00820 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00821 { 00822 double ThisCoef = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL); 00823 00824 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 00825 { 00826 /* The last line for each element is the total recombination for each temp. */ 00827 if( ipLo == NumLevRecomb[ipISO][nelem] ) 00828 { 00829 TotalRecomb[ipISO][nelem][i] = ThisCoef; 00830 } 00831 else 00832 RRCoef[ipISO][nelem][ipLo][i] = ThisCoef; 00833 } 00834 00835 if( lgEOL ) 00836 { 00837 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]); 00838 fprintf( ioQQQ, 00839 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" ); 00840 cdEXIT(EXIT_FAILURE); 00841 } 00842 } 00843 } 00844 00845 /* following loop only executed if we need more levels than are 00846 * stored in the recom coef data set 00847 * do not do collapsed levels since will use H-like recom there */ 00848 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 00849 { 00850 for( ipLo=NumLevRecomb[ipISO][nelem]; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ ) 00851 { 00852 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00853 { 00854 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */ 00855 RRCoef[ipISO][nelem][ipLo][i] = log10(iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo)); 00856 } 00857 } 00858 } 00859 } 00860 00861 /* check that ending magic number is ok */ 00862 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00863 { 00864 fprintf( ioQQQ, " iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]); 00865 cdEXIT(EXIT_FAILURE); 00866 } 00867 i = 1; 00868 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00869 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00870 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00871 i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00872 00873 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB ) 00874 { 00875 fprintf( ioQQQ, 00876 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] ); 00877 fprintf( ioQQQ, 00878 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" , 00879 RECOMBMAGIC , 00880 NumLevRecomb[ipISO][ipISO], 00881 NumLevRecomb[ipISO][ipISO+1], 00882 N_ISO_TE_RECOMB, 00883 i1 , i2 , i3, i4 ); 00884 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00885 fprintf( ioQQQ, 00886 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" ); 00887 cdEXIT(EXIT_FAILURE); 00888 } 00889 00890 /* close the data file */ 00891 fclose( ioDATA ); 00892 } 00893 } 00894 /* We are compiling the he_iso_recomb.dat data file. */ 00895 else if( iso.lgCompileRecomb[ipISO] ) 00896 { 00897 /* option to create table of recombination coefficients, 00898 * executed with the compile he-like command */ 00899 FILE *ioRECOMB; 00900 00901 ASSERT( !iso.lgNoRecombInterp[ipISO] ); 00902 00903 /*RECOMBMAGIC the magic number for the table of recombination coefficients */ 00904 /*NumLevRecomb[ipISO][nelem] the number of levels that will be done */ 00905 /* create recombination coefficients */ 00906 ioRECOMB = open_data( chFilename[ipISO], "w", AS_LOCAL_ONLY ); 00907 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", 00908 RECOMBMAGIC , 00909 NumLevRecomb[ipISO][ipISO], 00910 NumLevRecomb[ipISO][ipISO+1], 00911 N_ISO_TE_RECOMB, 00912 iso.chISO[ipISO], 00913 NumLevRecomb[ipISO][ipISO], 00914 elementnames.chElementSym[ipISO], 00915 NumLevRecomb[ipISO][ipISO+1], 00916 N_ISO_TE_RECOMB ); 00917 00918 for( nelem = ipISO; nelem < LIMELM; nelem++ ) 00919 { 00920 /* this must pass since compile xx-like command reset numLevels to the macro */ 00921 ASSERT( NumLevRecomb[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 00922 00923 /* Zero out the recombination sum array. */ 00924 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00925 { 00926 TotalRecomb[ipISO][nelem][i] = 0.; 00927 } 00928 00929 for( ipLo=ipHe1s1S; ipLo < NumLevRecomb[ipISO][nelem]; ipLo++ ) 00930 { 00931 fprintf(ioRECOMB, "%li\t%li", nelem, ipLo ); 00932 /* loop over temperatures to produce array of recombination coefficients */ 00933 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00934 { 00935 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */ 00936 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo); 00937 TotalRecomb[ipISO][nelem][i] += RadRecombReturn; 00938 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn); 00939 fprintf(ioRECOMB, "\t%f", RRCoef[ipISO][nelem][ipLo][i] ); 00940 } 00941 fprintf(ioRECOMB, "\n" ); 00942 } 00943 00944 /* Store one additional line in XX_iso_recomb.dat that gives the total recombination, 00945 * as computed by the sum so far, plus levels up to NHYDRO_MAX_LEVEL using Verner's fits, 00946 * plus levels up to SumUpToThisN using Seaton 59, for each element and each temperature. */ 00947 fprintf(ioRECOMB, "%li\t%li", nelem, NumLevRecomb[ipISO][nelem] ); 00948 for(i = 0; i < N_ISO_TE_RECOMB; i++) 00949 { 00950 for( i1 = ( (nelem == ipISO) ? (RREC_MAXN + 1) : (LIKE_RREC_MAXN( nelem ) + 1) ); i1< NHYDRO_MAX_LEVEL; i1++ ) 00951 { 00952 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i])); 00953 } 00954 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ ) 00955 { 00956 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 ); 00957 } 00958 fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[ipISO][nelem][i] ) ); 00959 } 00960 fprintf(ioRECOMB, "\n" ); 00961 } 00962 /* end the file with the same information */ 00963 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", 00964 RECOMBMAGIC , 00965 NumLevRecomb[ipISO][ipISO], 00966 NumLevRecomb[ipISO][ipISO+1], 00967 N_ISO_TE_RECOMB, 00968 iso.chISO[ipISO], 00969 NumLevRecomb[ipISO][ipISO], 00970 elementnames.chElementSym[ipISO], 00971 NumLevRecomb[ipISO][ipISO+1], 00972 N_ISO_TE_RECOMB ); 00973 00974 fclose( ioRECOMB ); 00975 00976 fprintf( ioQQQ, "iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] ); 00977 fprintf( ioQQQ, "The compilation is completed successfully.\n"); 00978 cdEXIT(EXIT_SUCCESS); 00979 } 00980 } 00981 00982 return; 00983 } 00984 00985 double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo ) 00986 { 00987 double rate; 00988 long ipTe, i; 00989 double TeDRCoef[NUM_DR_TEMPS]; 00990 const double Te_over_Z1_Squared[NUM_DR_TEMPS] = { 00991 1.00000, 1.30103, 1.69897, 2.00000, 2.30103, 2.69897, 3.00000, 00992 3.30103, 3.69897, 4.00000, 4.30103, 4.69897, 5.00000, 5.30103, 00993 5.69897, 6.00000, 6.30103, 6.69897, 7.00000 }; 00994 00995 DEBUG_ENTRY( "iso_dielec_recomb_rate()" ); 00996 00997 /* currently only two iso sequences and only he-like is applicable. */ 00998 ASSERT( ipISO == ipHE_LIKE ); 00999 ASSERT( ipLo >= 0 ); 01000 01001 /* temperature grid is nelem^2 * constant temperature grid above. */ 01002 for( i=0; i<NUM_DR_TEMPS; i++ ) 01003 { 01004 TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (double) nelem ); 01005 } 01006 01007 if( ipLo == ipHe1s1S ) 01008 { 01009 rate = 0.; 01010 } 01011 else if( ipLo<iso.numLevels_max[ipISO][nelem] ) 01012 { 01013 if( phycon.alogte <= TeDRCoef[0] ) 01014 { 01015 /* Take lowest tabulated value for low temperature end. */ 01016 rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][0]; 01017 } 01018 else if( phycon.alogte >= TeDRCoef[NUM_DR_TEMPS-1] ) 01019 { 01020 /* use T^-1.5 extrapolation at high temperatures. */ 01021 rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][NUM_DR_TEMPS-1] * 01022 pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] - phycon.alogte ) ) ; 01023 } 01024 else 01025 { 01026 /* find temperature in tabulated values. */ 01027 ipTe = hunt_bisect( TeDRCoef, NUM_DR_TEMPS, phycon.alogte ); 01028 01029 ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1) ); 01030 01031 if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1] == 0. ) 01032 rate = 0.; 01033 else if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe] == 0. ) 01034 rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1]; 01035 else 01036 { 01037 /* interpolate between tabulated points */ 01038 rate = log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]) + 01039 (phycon.alogte-TeDRCoef[ipTe])* 01040 (log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1])-log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]))/ 01041 (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]); 01042 01043 rate = pow( 10., rate ); 01044 } 01045 } 01046 } 01047 else 01048 { 01049 rate = 0.; 01050 } 01051 01052 ASSERT( rate >= 0. && rate < 1.0e-12 ); 01053 01054 return rate*iso.lgDielRecom[ipISO]; 01055 } 01056 01057 /* TempInterp - interpolate on an array */ 01059 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements ) 01060 { 01061 static long int ipTe=-1; 01062 double rate = 0.; 01063 long i0; 01064 01065 DEBUG_ENTRY( "TempInterp()" ); 01066 01067 ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 ); 01068 01069 if( ipTe < 0 ) 01070 { 01071 /* te totally unknown */ 01072 if( ( phycon.alogte < TempArray[0] ) || 01073 ( phycon.alogte > TempArray[NumElements-1] ) ) 01074 { 01075 fprintf(ioQQQ," TempInterp called with te out of bounds \n"); 01076 cdEXIT(EXIT_FAILURE); 01077 } 01078 ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte ); 01079 } 01080 else if( phycon.alogte < TempArray[ipTe] ) 01081 { 01082 /* temp is too low, must also lower ipTe */ 01083 ASSERT( phycon.alogte > TempArray[0] ); 01084 /* decrement ipTe until it is correct */ 01085 while( ( phycon.alogte < TempArray[ipTe] ) && ipTe > 0) 01086 --ipTe; 01087 } 01088 else if( phycon.alogte > TempArray[ipTe+1] ) 01089 { 01090 /* temp is too high */ 01091 ASSERT( phycon.alogte <= TempArray[NumElements-1] ); 01092 /* increment ipTe until it is correct */ 01093 while( ( phycon.alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1) 01094 ++ipTe; 01095 } 01096 01097 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) ); 01098 01099 /* ipTe should now be valid */ 01100 ASSERT( ( phycon.alogte >= TempArray[ipTe] ) 01101 && ( phycon.alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) ); 01102 01103 if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. ) 01104 { 01105 rate = 0.; 01106 } 01107 else 01108 { 01109 /* Do a four-point interpolation */ 01110 const int ORDER = 3; /* order of the fitting polynomial */ 01111 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0); 01112 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte ); 01113 } 01114 01115 return rate; 01116 }