cloudy trunk
Loading...
Searching...
No Matches
iso_radiative_recomb.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*iso_radiative_recomb find state specific creation and destruction rates for iso sequences */
4/*iso_RRCoef_Te - evaluate radiative recombination coef at some temperature */
5/*iso_recomb_check - called by SanityCheck to confirm that recombination coef are ok,
6 * return value is relative error between new calculation of recom, and interp value */
7
8#include "cddefines.h"
9#include "dynamics.h"
10#include "atmdat.h"
11#include "conv.h"
12#include "cosmology.h"
13#include "dense.h"
14#include "elementnames.h"
15#include "helike.h"
16#include "helike_recom.h"
17#include "hydrogenic.h"
18#include "ionbal.h"
19#include "iso.h"
20#include "opacity.h"
21#include "phycon.h"
22#include "physconst.h"
23#include "prt.h"
24#include "save.h"
25#include "thermal.h"
26#include "thirdparty.h"
27#include "trace.h"
28#include "rt.h"
29#include "taulines.h"
30
31/* this will save log of radiative recombination rate coefficients at N_ISO_TE_RECOMB temperatures.
32 * there will be NumLevRecomb[ipISO][nelem] levels saved in RRCoef[nelem][level][temp] */
33static double ****RRCoef/*[LIMELM][NumLevRecomb[ipISO][nelem]][N_ISO_TE_RECOMB]*/;
34static long **NumLevRecomb;
35static double ***TotalRecomb; /*[ipISO][nelem][i]*/
36
37/* the array of logs of temperatures at which RRCoef was defined */
39static double kTRyd,global_EthRyd;
40static long int globalZ,globalISO;
41static long int globalN, globalL, globalS;
42
43STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements );
44STATIC double iso_recomb_integrand(double EE);
45STATIC void iso_put_recomb_error( long ipISO, long nelem );
46
47double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
48{
49 double alpha,RecomIntegral=0.,b,E1,E2,step,OldRecomIntegral,TotChangeLastFive;
50 double change[5] = {0.,0.,0.,0.,0.};
51
52 DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" );
53
54 if( ipISO==ipH_LIKE && ipLo == 0 )
55 return t_ADfA::Inst().H_rad_rec(nelem+1,ipLo,temp);
56
57 global_EthRyd = iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd;
58
59 /* Factors outside integral in Milne relation. */
60 b = MILNE_CONST * iso_sp[ipISO][nelem].st[ipLo].g() * pow(temp,-1.5);
61
62 if( ipISO==ipH_LIKE )
63 b /= 2.;
64 else if( ipISO==ipHE_LIKE )
65 b /= 4.;
66
67 /* kT in Rydbergs. */
68 kTRyd = temp / TE1RYD;
69 globalISO = ipISO;
70 globalZ = nelem;
71 globalN = N_(ipLo);
72 globalL = L_(ipLo);
73 globalS = S_(ipLo);
74
75 /* Begin integration. */
76 /* First define characteristic step */
77 E1 = global_EthRyd;
78
79 if( ipISO==ipH_LIKE )
80 step = MIN2( 0.125*kTRyd, 0.5*E1 );
81 else if( ipISO==ipHE_LIKE )
82 step = MIN2( 0.25*kTRyd, 0.5*E1 );
83 else
85
86 E2 = E1 + step;
87 /* Perform initial integration, from threshold to threshold + step. */
88 RecomIntegral = qg32( E1, E2, iso_recomb_integrand);
89 /* Repeat the integration, adding each new result to the total,
90 * except that the step size is doubled every time, since values away from
91 * threshold tend to fall off more slowly. */
92 do
93 {
94 OldRecomIntegral = RecomIntegral;
95 E1 = E2;
96 step *= 1.25;
97 E2 = E1 + step;
98 RecomIntegral += qg32( E1, E2, iso_recomb_integrand);
99 change[4] = change[3];
100 change[3] = change[2];
101 change[2] = change[1];
102 change[1] = change[0];
103 change[0] = (RecomIntegral - OldRecomIntegral)/RecomIntegral;
104 TotChangeLastFive = change[0] + change[1] + change[2] + change[3] + change[4];
105 /* Continue integration until the upper limit exceeds 100kTRyd, an arbitrary
106 * point at which the integrand component exp(electron energy/kT) is very small,
107 * making neglible cross-sections at photon energies beyond that point,
108 * OR when the last five steps resulted in less than a 1 percent change. */
109 } while ( ((E2-global_EthRyd) < 100.*kTRyd) && ( TotChangeLastFive > 0.0001) );
110
111 /* Calculate recombination coefficient. */
112 alpha = b * RecomIntegral;
113
114 alpha = MAX2( alpha, SMALLDOUBLE );
115
116 return alpha;
117}
118
119/*iso_recomb_integrand, used in Milne relation for iso sequences - the energy is photon Rydbergs. */
120STATIC double iso_recomb_integrand(double ERyd)
121{
122 double x1, temp;
123
124 /* Milne relation integrand */
125 x1 = ERyd * ERyd * exp(-1.0 * ( ERyd - global_EthRyd ) / kTRyd);
127 x1 *= temp;
128
129 return x1;
130}
131
132double iso_cross_section( double EgammaRyd , double EthRyd, long n, long l, long S, long globalZ, long globalISO )
133{
134 double cross_section;
135 DEBUG_ENTRY( "iso_cross_section()" );
136
137 if( globalISO == ipH_LIKE )
138 cross_section = H_cross_section( EgammaRyd , EthRyd, n, l, globalZ );
139 else if( globalISO == ipHE_LIKE )
140 cross_section = He_cross_section( EgammaRyd , EthRyd, n, l, S, globalZ );
141 else
143
144 return cross_section;
145}
146
147/*=======================================================*/
148/* iso_radiative_recomb get rad recomb rate coefficients for iso sequences */
150 long ipISO,
151 /* nelem on the c scale, He is 1 */
152 long nelem )
153{
154 long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel;
155 double topoff, TotMinusExplicitResolved,
156 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
157 double RecExplictLevels, TotalRadRecomb, RecCollapsed;
158 static double TeUsed[NISO][LIMELM]={
159 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
160 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
161 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
162 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
163 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
164 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
165
166 DEBUG_ENTRY( "iso_radiative_recomb()" );
167
168 iso_put_recomb_error( ipISO, nelem );
169
170 /* evaluate recombination escape probability for all levels */
171
172 /* define radiative recombination rates for all levels */
173 /* this will be the sum of all levels explicitly included in the model atom */
174 RecExplictLevels = 0.;
175
176 /* number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
177 ipFirstCollapsed = iso_sp[ipISO][nelem].numLevels_local-iso_sp[ipISO][nelem].nCollapsed_local;
178
179 ASSERT( ipFirstCollapsed == iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local );
180 if( !fp_equal(phycon.te,TeUsed[ipISO][nelem]) || iso_sp[ipISO][nelem].lgMustReeval || !conv.nTotalIoniz || cosmology.lgDo )
181 {
182 TeUsed[ipISO][nelem] = phycon.te;
183
184 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
185 {
186 /* this is radiative recombination rate coefficient */
187 double RadRec;
188
189 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
190 {
191 RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel );
192 }
193 else
194 {
195 RadRec = iso_radrecomb_from_cross_section( ipISO, phycon.te, nelem, ipLevel);
196 }
197 ASSERT( RadRec > 0. );
198
199 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] = RadRec;
200
201 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] > 0. );
202
203 RecExplictLevels += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
204
205 if( iso_ctrl.lgDielRecom[ipISO] )
206 {
207 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
208 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
209 Total_DR_Added += iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
210 }
211 }
212
213 /**************************************************/
214 /*** Add on recombination to collapsed levels ***/
215 /**************************************************/
216 RecCollapsed = 0.;
217 for( ipLevel=ipFirstCollapsed; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
218 {
219 /* use hydrogenic for collapsed levels */
220 double RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, N_(ipLevel), phycon.te);
221
222 /* this is radiative recombination rate coefficient */
223 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] = RadRec;
224
225 RecCollapsed += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
226
227 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] > 0. );
228
229 if( iso_ctrl.lgDielRecom[ipISO] )
230 {
231 /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
232 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
233 Total_DR_Added += iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
234 }
235 }
236 for( ipLevel=iso_sp[ipISO][nelem].numLevels_local; ipLevel<iso_sp[ipISO][nelem].numLevels_max;++ipLevel )
237 {
238 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = 0.;
239 }
240 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
241 for( ipLevel = 0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
242 {
243 if( N_(ipLevel) == ThisN )
244 {
245 TotRRThisN += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
246 }
247 else
248 {
249 ASSERT( iso_ctrl.lgRandErrGen[ipISO] || (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) || (ThisN==iso_sp[ipISO][nelem].n_HighestResolved_local+1) );
250 LastN = ThisN;
251 ThisN = N_(ipLevel);
252 TotRRLastN = TotRRThisN;
253 TotRRThisN = iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
254
255 {
256 /* Print the sum of recombination coefficients per n at current temp. */
257 /*@-redef@*/
258 enum {DEBUG_LOC=false};
259 /*@+redef@*/
260 static long RUNONCE = false;
261
262 if( !RUNONCE && DEBUG_LOC )
263 {
264 static long FIRSTTIME = true;
265
266 if( FIRSTTIME )
267 {
268 fprintf( ioQQQ,"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n",
269 ipISO, nelem, phycon.te);
270 FIRSTTIME= false;
271 }
272
273 fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN );
274 }
275 RUNONCE = true;
276 }
277 }
278 }
279
280 /* Get total recombination into all levels, including those not explicitly considered. */
281 if( iso_ctrl.lgNoRecombInterp[ipISO] )
282 {
283 /* We are not interpolating, must calculate total now, as sum of resolved and collapsed levels... */
284 TotalRadRecomb = RecCollapsed + RecExplictLevels;
285
286 /* Plus additional levels out to a predefined limit... */
287 for( long nLo = N_(ipFirstCollapsed) + iso_sp[ipISO][nelem].nCollapsed_max; nLo < NHYDRO_MAX_LEVEL; nLo++ )
288 {
289 TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te);
290 }
291 /* Plus a bunch more for good measure */
292 for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ )
293 {
294 TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo );
295 }
296 }
297 else if( iso_sp[ipISO][nelem].lgLevelsLowered )
298 {
299 TotalRadRecomb = 0.;
300 for( ipLevel = 0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
301 TotalRadRecomb += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad];
302 }
303 else
304 {
305 /* We are interpolating, and total was calculated here in iso_recomb_setup */
306 TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max );
307 }
308 if(ipISO==0 && nelem==0 )
309 {
310 // insure rec coef will be always evaluated
311 TeUsed[ipISO][nelem] = phycon.te*0.999;
312
313 // if( iteration > dynamics.n_initial_relax+2 && fudge(-1)>0 )
314 // TotalRadRecomb *= fudge(0);
315 }
316
317 /* If generating random error, apply one to total recombination */
318 if( iso_ctrl.lgRandErrGen[ipISO] )
319 {
320 /* Error for total recombination */
321 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,iso_sp[ipISO][nelem].numLevels_max,IPRAD,0.0001f,0.0001f);
322 //iso_sp[ipISO][nelem].ErrorFactor[iso_sp[ipISO][nelem].numLevels_max][iso_sp[ipISO][nelem].numLevels_max][IPRAD] =
323 // (realnum)MyGaussRand( iso_sp[ipISO][nelem].Error[iso_sp[ipISO][nelem].numLevels_max][iso_sp[ipISO][nelem].numLevels_max][IPRAD] );
324
325 /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
326 * the error factors for rrc are always stored at iso.numLevels_max, regardless of
327 * continuum lowering effects. */
328 //TotalRadRecomb *=
329 // iso_sp[ipISO][nelem].ErrorFactor[iso_sp[ipISO][nelem].numLevels_max][iso_sp[ipISO][nelem].numLevels_max][IPRAD];
330 }
331
332 /* this is case B recombination, sum without the ground included */
333 iso_sp[ipISO][nelem].RadRec_caseB = TotalRadRecomb - iso_sp[ipISO][nelem].fb[0].RadRecomb[ipRecRad];
334 ASSERT( iso_sp[ipISO][nelem].RadRec_caseB > 0.);
335
336 // Restore highest levels opacity stack
337 for( unsigned i = 0; i < iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
338 {
339 long index = iso_sp[ipISO][nelem].numLevels_max-1;
340 opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] = iso_sp[ipISO][nelem].HighestLevelOpacStack[i];
341 }
342
343 /* Add topoff (excess) recombination to top level. This is only done if atom is full size,
344 * that is, no pressure lowered of the continuum the current conditions. Radiative recombination
345 * to non-existent states does not occur, those would dominate the topoff */
346 if( !iso_sp[ipISO][nelem].lgLevelsLowered )
347 {
348 /* at this point we have RecExplictLevels, the sum of radiative recombination
349 * to all levels explicitly included in the model atom the total
350 * recombination rate. The difference is the amount of "topoff" we will need to do */
351 TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
352
353 topoff = TotMinusExplicitResolved - RecCollapsed;
354
355 /* the t_ADfA::Inst().rad_rec fits are too small at high temperatures, so this atom is
356 * better than the topoff. Only a problem for helium itself, at high temperatures.
357 * complain if the negative topoff is not for this case */
358 if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) &&
359 fabs( topoff/TotalRadRecomb ) > 0.01 )
360 {
361 fprintf(ioQQQ," PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",
362 nelem, topoff/TotalRadRecomb, phycon.te );
363 }
364
365 // option to turn off topoff with atom xx-like topoff off
366 if( !iso_ctrl.lgTopoff[ipISO] )
367 topoff *= 1E-20;
368
369 topoff = MAX2( 0., topoff );
370 double scale_factor = 1. + topoff/iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].RadRecomb[ipRecRad];
371 ASSERT( scale_factor >= 1. );
372
373 // Scale highest level opacities to be consistent with recombination topoff
374 for( unsigned i = 0; i < iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
375 {
376 long index = iso_sp[ipISO][nelem].numLevels_max-1;
377 opac.OpacStack[iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] *= scale_factor;
378 }
379
380 /* We always have at least one collapsed level if continuum is not lowered. Put topoff there. */
381 iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].RadRecomb[ipRecRad] += topoff;
382
383 /* check for negative DR topoff, but only if Total_DR_Added is not negligible compared with TotalRadRecomb */
384 if( Total_DR_Added > TotalRadRecomb/100. )
385 {
386 if( Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 )
387 {
388 fprintf(ioQQQ," PROBLEM negative DR topoff for %li, tot1, tot2 = %.2e\t%.2e rel error was %.2e, temp was %.2f\n",
389 nelem,
390 Total_DR_Added,
391 ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO],
392 Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0,
393 phycon.te );
394 }
395 }
396
397 ASSERT( iso_sp[ipISO][nelem].numLevels_max == iso_sp[ipISO][nelem].numLevels_local );
398
399 if( iso_ctrl.lgDielRecom[ipISO] && iso_ctrl.lgTopoff[ipISO] )
400 {
401 /* \todo 2 suppress this total rate for continuum lowering using factors from Jordan (1969). */
402 /* put extra DR in top level */
403 iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_max-1].DielecRecomb +=
404 MAX2( 0., ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added );
405 }
406 }
407
408 }
409
410 /**************************************************************/
411 /*** Stuff escape probabilities, and effective rad recomb ***/
412 /**************************************************************/
413
414 /* total effective radiative recombination, initialize to zero */
415 iso_sp[ipISO][nelem].RadRec_effec = 0.;
416
417 for( ipLevel=0; ipLevel<iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
418 {
419 /* option for case b conditions, kill ground state recombination */
420 if( opac.lgCaseB && ipLevel==0 )
421 {
422 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 1e-10;
423 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 1e-10;
424 }
425 else if( cosmology.lgDo && ipLevel==0 )
426 {
427 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 1e-30;
428 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 1e-30;
429 }
430 else
431 {
432 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] =
433 RT_recom_effic(iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon);
434
435 /* net escape prob includes dest by background opacity */
436 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] =
437 MIN2( 1.,
438 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] +
439 (1.-iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc])*
440 iso_sp[ipISO][nelem].fb[ipLevel].ConOpacRatio );
441 }
442
443 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] >= 0. );
444 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] >= 0. );
445 ASSERT( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] <= 1. );
446
447 /* sum of all effective rad rec */
448 iso_sp[ipISO][nelem].RadRec_effec += iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad]*
449 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc];
450 }
451
452 /* zero out escape probabilities of levels above numLevels_local */
453 for( ipLevel=iso_sp[ipISO][nelem].numLevels_local; ipLevel<iso_sp[ipISO][nelem].numLevels_max; ++ipLevel )
454 {
455 /* this is escape probability */
456 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] = 0.;
457 /* net escape prob includes dest by background opacity */
458 iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc] = 0.;
459 }
460
461 /* trace escape probabilities */
462 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
463 {
464 fprintf( ioQQQ, " iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
465
466 /* print continuum escape probability */
467 fprintf( ioQQQ, " iso_radiative_recomb recomb effic" );
468 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
469 {
470 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecEsc] ));
471 }
472 fprintf( ioQQQ, "\n" );
473
474 /* net recombination efficiency factor, including background opacity*/
475 fprintf( ioQQQ, " iso_radiative_recomb recomb net effic" );
476 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
477 {
478 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecNetEsc]) );
479 }
480 fprintf( ioQQQ, "\n" );
481
482 /* inward continuous optical depths */
483 fprintf( ioQQQ, " iso_radiative_recomb in optic dep" );
484 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
485 {
486 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[0][iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon-1] ));
487 }
488 fprintf( ioQQQ, "\n" );
489
490 /* outward continuous optical depths*/
491 fprintf( ioQQQ, " iso_radiative_recomb out op depth" );
492 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
493 {
494 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[1][iso_sp[ipISO][nelem].fb[ipLevel].ipIsoLevNIonCon-1] ));
495 }
496 fprintf( ioQQQ, "\n" );
497
498 /* print radiative recombination coefficients */
499 fprintf( ioQQQ, " iso_radiative_recomb rad rec coef " );
500 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
501 {
502 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad]) );
503 }
504 fprintf( ioQQQ, "\n" );
505 }
506
507 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
508 {
509 fprintf( ioQQQ, " iso_radiative_recomb total rec coef" );
510 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].RadRec_effec ));
511 fprintf( ioQQQ, " case A=" );
512 fprintf( ioQQQ,PrintEfmt("%10.3e",
513 iso_sp[ipISO][nelem].RadRec_caseB + iso_sp[ipISO][nelem].fb[ipH1s].RadRecomb[ipRecRad] ) );
514 fprintf( ioQQQ, " case B=");
515 fprintf( ioQQQ,PrintEfmt("%10.3e", iso_sp[ipISO][nelem].RadRec_caseB ));
516 fprintf( ioQQQ, "\n" );
517 }
518
519 /****************************/
520 /*** begin sanity check ***/
521 /****************************/
522 {
523 bool lgOK = true;
524 for( ipLevel=0; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
525 {
526 if( iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] <= 0. )
527 {
528 fprintf( ioQQQ,
529 " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
530 ipISO, nelem, ipLevel, iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[ipRecRad] , phycon.te);
531 lgOK = false;
532 }
533 }
534 /* bail if we found problems */
535 if( !lgOK )
536 {
537 ShowMe();
539 }
540 /*end sanity check */
541 }
542
543 /* confirm that we have good rec coef at bottom and top of atom/ion */
544 ASSERT( iso_sp[ipISO][nelem].fb[0].RadRecomb[ipRecRad] > 0. );
545 ASSERT( iso_sp[ipISO][nelem].fb[iso_sp[ipISO][nelem].numLevels_local-1].RadRecomb[ipRecRad] > 0. );
546
547 /* set true when save recombination coefficients command entered */
548 if( save.lgioRecom )
549 {
550 /* this prints Z on physical, not C, scale */
551 fprintf( save.ioRecom, "%s %s %2li ",
552 iso_ctrl.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 );
553 fprintf( save.ioRecom,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].RadRec_effec ));
554 fprintf( save.ioRecom, "\n" );
555 }
556
557 return;
558}
559
560STATIC void iso_put_recomb_error( long ipISO, long nelem )
561{
562 long level;
563
564 /* optimistic and pessimistic errors for HeI recombination */
565 realnum He_errorOpti[31] = {
566 0.0000f,
567 0.0009f,0.0037f,0.0003f,0.0003f,0.0003f,0.0018f,
568 0.0009f,0.0050f,0.0007f,0.0003f,0.0001f,0.0007f,
569 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f,
570 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f };
571
572 realnum He_errorPess[31] = {
573 0.0100f,
574 0.0100f,0.0060f,0.0080f,0.0080f,0.0080f,0.0200f,
575 0.0200f,0.0200f,0.0200f,0.0600f,0.0600f,0.0080f,
576 0.0200f,0.0200f,0.0070f,0.0100f,0.0100f,0.0020f,0.0030f,0.0070f,
577 0.0300f,0.0300f,0.0100f,0.0200f,0.0200f,0.0200f,0.0200f,0.0010f,0.0004f,0.0090f };
578
579 /* now put recombination errors into iso.Error[ipISO] array */
580 for( long ipHi=0; ipHi<iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
581 {
582 if( ipISO==ipHE_LIKE && nelem==ipHELIUM )
583 {
584 if( ipHi >= iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max )
585 level = 21;
586 else if( ipHi<=30 )
587 level = ipHi;
588 else
589 level = iso_sp[ipISO][nelem].QuantumNumbers2Index[5][ MIN2(L_(ipHi),4) ][S_(ipHi)];
590
591 ASSERT( level <=30 );
592
593 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,ipHi,IPRAD, He_errorOpti[level], He_errorPess[level]);
594 }
595 else
596 iso_put_error(ipISO,nelem,iso_sp[ipISO][nelem].numLevels_max,ipHi,IPRAD,0.1f,0.1f);
597 }
598
599 return;
600}
601
602void iso_radiative_recomb_effective( long ipISO, long nelem )
603{
604 DEBUG_ENTRY( "iso_radiative_recomb_effective()" );
605
606 /* Find effective recombination coefficients */
607 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
608 {
609 iso_sp[ipISO][nelem].fb[ipHi].RadEffec = 0.;
610
611 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
612 for( long ipHigher=ipHi; ipHigher < iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
613 {
614 ASSERT( iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] >= 0. );
615 ASSERT( iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad] >= 0. );
616
617 iso_sp[ipISO][nelem].fb[ipHi].RadEffec += iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] *
618 iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad];
619 }
620 }
621
622 /**************************************************************/
623 /*** option to print effective recombination coefficients ***/
624 /**************************************************************/
625 {
626 enum {DEBUG_LOC=false};
627
628 if( DEBUG_LOC )
629 {
630 const int maxPrt=10;
631
632 fprintf( ioQQQ,"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem, phycon.te );
633 fprintf( ioQQQ, "N\tL\tS\tRadEffec\tLifetime\n" );
634
635 for( long i=0; i<maxPrt; i++ )
636 {
637 fprintf( ioQQQ, "%li\t%li\t%li\t%e\t%e\n", N_(i), L_(i), S_(i),
638 iso_sp[ipISO][nelem].fb[i].RadEffec,
639 MAX2( 0., iso_sp[ipISO][nelem].st[i].lifetime() ) );
640 }
641 fprintf( ioQQQ, "\n" );
642 }
643 }
644
645 /* If we have the variable set, find errors in rad rates */
646 if( iso_ctrl.lgRandErrGen[ipISO] )
647 {
648 dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
649
650 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
651 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
652 {
653 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = 0.;
654
655 /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
656 for( long ipHigher=ipHi; ipHigher < iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
657 {
658 ASSERT( iso_sp[ipISO][nelem].ex[iso_sp[ipISO][nelem].numLevels_max][ipHigher].Error[IPRAD] >= 0. );
659 ASSERT( iso_sp[ipISO][nelem].ex[ipHigher][ipHi].SigmaCascadeProb >= 0. );
660
661 /* iso.RadRecomb has to appear here because iso.Error is only relative error */
662 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec += pow( iso_sp[ipISO][nelem].ex[iso_sp[ipISO][nelem].numLevels_max][ipHigher].Error[IPRAD] *
663 iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] * iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad], 2.) +
664 pow( iso_sp[ipISO][nelem].ex[ipHigher][ipHi].SigmaCascadeProb * iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[ipRecRad], 2.);
665 }
666
667 ASSERT( iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec >= 0. );
668 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = sqrt( iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec );
669
670 for( long ipLo = 0; ipLo < ipHi; ipLo++ )
671 {
672 if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) )
673 {
674 double EnergyInRydbergs = iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd - iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd;
675 realnum wavelength = (realnum)(RYDLAM/MAX2(1E-8,EnergyInRydbergs));
676 double emissivity = iso_sp[ipISO][nelem].fb[ipHi].RadEffec * iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] * EN1RYD * EnergyInRydbergs;
677 double sigma_emiss = 0., SigmaBranchRatio = 0.;
678
679 if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (N_(ipHi)<=5) )
680 {
681 SigmaBranchRatio = iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] * sqrt(
682 pow( (double)iso_sp[ipISO][nelem].ex[ipHi][ipLo].Error[IPRAD], 2. ) +
683 pow( iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) );
684
685 sigma_emiss = EN1RYD * EnergyInRydbergs * sqrt(
686 pow( (double)iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec * iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo], 2.) +
687 pow( SigmaBranchRatio * iso_sp[ipISO][nelem].fb[ipHi].RadEffec, 2.) );
688
689 /* \todo 2 make this a trace option */
690 dprintf( ioQQQ, "%li\t%li\t", ipHi, ipLo );
692 fprintf( ioQQQ, "\t%e\t%e\t%e\t%e\t%e\t%e\n",
693 emissivity,
694 sigma_emiss,
695 iso_sp[ipISO][nelem].fb[ipHi].RadEffec,
696 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec,
697 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo],
698 SigmaBranchRatio);
699 }
700 }
701 }
702 }
703 }
704
705 return;
706}
707/*iso_RRCoef_Te evaluated radiative recombination coef at some temperature */
708double iso_RRCoef_Te( long ipISO, long nelem , long n )
709{
710 double rate = 0.;
711
712 DEBUG_ENTRY( "iso_RRCoef_Te()" );
713
714 ASSERT( !iso_ctrl.lgNoRecombInterp[ipISO] );
715
716 /* if n is equal to the number of levels, return the total recomb, else recomb for given level. */
717 if( n == iso_sp[ipISO][nelem].numLevels_max - iso_sp[ipISO][nelem].nCollapsed_max )
718 {
719 rate = TempInterp( TeRRCoef, TotalRecomb[ipISO][nelem], N_ISO_TE_RECOMB );
720 }
721 else
722 {
723 rate = TempInterp( TeRRCoef, RRCoef[ipISO][nelem][n], N_ISO_TE_RECOMB );
724 }
725
726 /* that was the log, now make linear */
727 rate = pow( 10. , rate );
728
729 return rate;
730}
731
732/*iso_recomb_check called by SanityCheck to confirm that recombination coef are ok
733 * return value is relative error between new calculation of recom, and interp value */
734double iso_recomb_check( long ipISO, long nelem, long level, double temperature )
735{
736 double RecombRelError ,
737 RecombInterp,
738 RecombCalc,
739 SaveTemp;
740
741 DEBUG_ENTRY( "iso_recomb_check()" );
742
743 /* first set temp to desired value */
744 SaveTemp = phycon.te;
745 // uses overloaded version of TempChange that does not check
746 // on floor since this is just a sanity check
747 TempChange(temperature);
748
749 /* actually calculate the recombination coefficient from the Milne relation,
750 * normally only done due to compile he-like command */
751 RecombCalc = iso_radrecomb_from_cross_section( ipISO, temperature , nelem , level );
752
753 /* interpolate the recombination coefficient, this is the usual method */
754 RecombInterp = iso_RRCoef_Te( ipISO, nelem , level );
755
756 /* reset temp */
757 TempChange(SaveTemp);
758
759 RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc );
760
761 return RecombRelError;
762}
763
764/* malloc space needed for iso recombination tables */
766{
767 DEBUG_ENTRY( "iso_recomb_malloc()" );
768
769 NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO );
770 TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO );
771 RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO );
772
773 for( long ipISO=0; ipISO<NISO; ipISO++ )
774 {
775 TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
776 RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
777 /* The number of recombination coefficients to be read from file for each element. */
778 NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM );
779
780 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
781 {
782 long int MaxLevels, maxN;
783
784 TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
785
786 if( nelem == ipISO )
787 maxN = RREC_MAXN;
788 else
789 maxN = LIKE_RREC_MAXN( nelem );
790
791 NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 );
792
793 if( nelem == ipISO || dense.lgElmtOn[nelem] )
794 {
795 /* must always have at least NumLevRecomb[ipISO][nelem] levels since that is number
796 * that will be read in from he rec data file, but possible to require more */
797 MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso_sp[ipISO][nelem].numLevels_max );
798
799 /* always define this */
800 /* >>chng 02 jan 24, RRCoef will be iso_sp[ipISO][nelem].numLevels_max levels, not iso.numLevels_max,
801 * code will stop if more than this is requested */
802 RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) );
803
804 for( long ipLo=0; ipLo < MaxLevels;++ipLo )
805 {
806 RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
807 }
808 }
809 }
810 }
811
812 for(long i = 0; i < N_ISO_TE_RECOMB; i++)
813 {
814 /* this is the vector of temperatures */
815 TeRRCoef[i] = 0.25*(i);
816 }
817
818 /* >>chng 06 jun 06, NP, assert thrown at T == 1e10 K, just bump the
819 * high temperature end slightly. */
820 TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f;
821
822 return;
823}
824
826{
827 DEBUG_ENTRY( "iso_recomb_auxiliary_free()" );
828
829 for( long i = 0; i< NISO; i++ )
830 {
831 free( NumLevRecomb[i] );
832 }
833 free( NumLevRecomb );
834
835 return;
836}
837
838void iso_recomb_setup( long ipISO )
839{
840 double RadRecombReturn;
841 long int i, i1, i2, i3, i4, i5;
842 long int ipLo, nelem;
843
844 const int chLine_LENGTH = 1000;
845 char chLine[chLine_LENGTH];
846 /* this must be longer than data path string, set in path.h/cpu.cpp */
847 const char* chFilename[NISO] =
848 { "h_iso_recomb.dat", "he_iso_recomb.dat" };
849
850 FILE *ioDATA;
851 bool lgEOL;
852
853 DEBUG_ENTRY( "iso_recomb_setup()" );
854
855 /* if we are compiling the recombination data file, we must interpolate in temperature */
856 if( iso_ctrl.lgCompileRecomb[ipISO] )
857 {
858 iso_ctrl.lgNoRecombInterp[ipISO] = false;
859 }
860
861 if( !iso_ctrl.lgNoRecombInterp[ipISO] )
862 {
863 /******************************************************************/
865 /******************************************************************/
866 /* This flag says we are not compiling the data file */
867 if( !iso_ctrl.lgCompileRecomb[ipISO] )
868 {
869 if( trace.lgTrace )
870 fprintf( ioQQQ," iso_recomb_setup opening %s:", chFilename[ipISO] );
871
872 /* Now try to read from file...*/
873 try
874 {
875 ioDATA = open_data( chFilename[ipISO], "r" );
876 }
877 catch( cloudy_exit )
878 {
879 fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
880 fprintf( ioQQQ, " but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
881 ioDATA = NULL;
882 }
883 if( ioDATA == NULL )
884 {
885 /* Do on the fly computation of R.R. Coef's instead. */
886 for( nelem = ipISO; nelem < LIMELM; nelem++ )
887 {
888 if( dense.lgElmtOn[nelem] )
889 {
890 /* Zero out the recombination sum array. */
891 for(i = 0; i < N_ISO_TE_RECOMB; i++)
892 {
893 TotalRecomb[ipISO][nelem][i] = 0.;
894 }
895
896 /* NumLevRecomb[ipISO][nelem] corresponds to n = 40 for H and He and 20 for ions, at present */
897 /* There is no need to fill in values for collapsed levels, because we do not need to
898 * interpolate for a given temperature, just calculate it directly with a hydrogenic routine. */
899 for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
900 {
901 /* loop over temperatures to produce array of recombination coefficients */
902 for(i = 0; i < N_ISO_TE_RECOMB; i++)
903 {
904 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
905 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
906 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
907 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
908 }
909 }
910 for(i = 0; i < N_ISO_TE_RECOMB; i++)
911 {
912 for( i1 = iso_sp[ipISO][nelem].n_HighestResolved_max+1; i1< NHYDRO_MAX_LEVEL; i1++ )
913 {
914 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
915 }
916 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
917 {
918 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
919 }
920 TotalRecomb[ipISO][nelem][i] = log10( TotalRecomb[ipISO][nelem][i] );
921 }
922 }
923 }
924 }
925 /* Data file is present and readable...begin read. */
926 else
927 {
928 /* check that magic number is ok */
929 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
930 {
931 fprintf( ioQQQ, " iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
933 }
934 i = 1;
935 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
936 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
937 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
938 i4 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
939 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
940 {
941 fprintf( ioQQQ,
942 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
943 fprintf( ioQQQ,
944 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
946 NumLevRecomb[ipISO][ipISO],
947 NumLevRecomb[ipISO][ipISO+1],
949 i1 , i2 , i3, i4 );
950 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
951 fprintf( ioQQQ,
952 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
954 }
955
956 i5 = 1;
957 /* now read in the data */
958 for( nelem = ipISO; nelem < LIMELM; nelem++ )
959 {
960 for( ipLo=0; ipLo <= NumLevRecomb[ipISO][nelem]; ipLo++ )
961 {
962 i5++;
963 /* get next line image */
964 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
965 {
966 fprintf( ioQQQ, " iso_recomb_setup could not read line %li of %s.\n", i5,
967 chFilename[ipISO] );
969 }
970 /* each line starts with element and level number */
971 i3 = 1;
972 i1 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
973 i2 = (long)FFmtRead(chLine,&i3,sizeof(chLine),&lgEOL);
974 /* check that these numbers are correct */
975 if( i1!=nelem || i2!=ipLo )
976 {
977 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
978 fprintf( ioQQQ,
979 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
981 }
982
983 /* loop over temperatures to produce array of recombination coefficients */
984 for(i = 0; i < N_ISO_TE_RECOMB; i++)
985 {
986 double ThisCoef = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
987
988 if( nelem == ipISO || dense.lgElmtOn[nelem] )
989 {
990 /* The last line for each element is the total recombination for each temp. */
991 if( ipLo == NumLevRecomb[ipISO][nelem] )
992 {
993 TotalRecomb[ipISO][nelem][i] = ThisCoef;
994 }
995 else
996 RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
997 }
998
999 if( lgEOL )
1000 {
1001 fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
1002 fprintf( ioQQQ,
1003 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
1005 }
1006 }
1007 }
1008
1009 /* following loop only executed if we need more levels than are
1010 * stored in the recom coef data set
1011 * do not do collapsed levels since will use H-like recom there */
1012 if( nelem == ipISO || dense.lgElmtOn[nelem] )
1013 {
1014 for( ipLo=NumLevRecomb[ipISO][nelem]; ipLo < iso_sp[ipISO][nelem].numLevels_max-iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
1015 {
1016 for(i = 0; i < N_ISO_TE_RECOMB; i++)
1017 {
1018 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
1019 RRCoef[ipISO][nelem][ipLo][i] = log10(iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo));
1020 }
1021 }
1022 }
1023 }
1024
1025 /* check that ending magic number is ok */
1026 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1027 {
1028 fprintf( ioQQQ, " iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
1030 }
1031 i = 1;
1032 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1033 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1034 i3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1035 i4 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1036
1037 if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
1038 {
1039 fprintf( ioQQQ,
1040 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
1041 fprintf( ioQQQ,
1042 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
1043 RECOMBMAGIC ,
1044 NumLevRecomb[ipISO][ipISO],
1045 NumLevRecomb[ipISO][ipISO+1],
1047 i1 , i2 , i3, i4 );
1048 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1049 fprintf( ioQQQ,
1050 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
1052 }
1053
1054 /* close the data file */
1055 fclose( ioDATA );
1056 }
1057 }
1058 /* We are compiling the he_iso_recomb.dat data file. */
1059 else if( iso_ctrl.lgCompileRecomb[ipISO] )
1060 {
1061 /* option to create table of recombination coefficients,
1062 * executed with the compile he-like command */
1063 FILE *ioRECOMB;
1064
1065 ASSERT( !iso_ctrl.lgNoRecombInterp[ipISO] );
1066
1067 /*RECOMBMAGIC the magic number for the table of recombination coefficients */
1068 /*NumLevRecomb[ipISO][nelem] the number of levels that will be done */
1069 /* create recombination coefficients */
1070 ioRECOMB = open_data( chFilename[ipISO], "w", AS_LOCAL_ONLY );
1071 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",
1072 RECOMBMAGIC ,
1073 NumLevRecomb[ipISO][ipISO],
1074 NumLevRecomb[ipISO][ipISO+1],
1076 iso_ctrl.chISO[ipISO],
1077 NumLevRecomb[ipISO][ipISO],
1078 elementnames.chElementSym[ipISO],
1079 NumLevRecomb[ipISO][ipISO+1],
1081
1082 for( nelem = ipISO; nelem < LIMELM; nelem++ )
1083 {
1084 /* this must pass since compile xx-like command reset numLevels to the macro */
1085 ASSERT( NumLevRecomb[ipISO][nelem] <= iso_sp[ipISO][nelem].numLevels_max );
1086
1087 /* Zero out the recombination sum array. */
1088 for(i = 0; i < N_ISO_TE_RECOMB; i++)
1089 {
1090 TotalRecomb[ipISO][nelem][i] = 0.;
1091 }
1092
1093 for( ipLo=ipHe1s1S; ipLo < NumLevRecomb[ipISO][nelem]; ipLo++ )
1094 {
1095 fprintf(ioRECOMB, "%li\t%li", nelem, ipLo );
1096 /* loop over temperatures to produce array of recombination coefficients */
1097 for(i = 0; i < N_ISO_TE_RECOMB; i++)
1098 {
1099 /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
1100 RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
1101 TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
1102 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
1103 fprintf(ioRECOMB, "\t%f", RRCoef[ipISO][nelem][ipLo][i] );
1104 }
1105 fprintf(ioRECOMB, "\n" );
1106 }
1107
1108 /* Store one additional line in XX_iso_recomb.dat that gives the total recombination,
1109 * as computed by the sum so far, plus levels up to NHYDRO_MAX_LEVEL using Verner's fits,
1110 * plus levels up to SumUpToThisN using Seaton 59, for each element and each temperature. */
1111 fprintf(ioRECOMB, "%li\t%li", nelem, NumLevRecomb[ipISO][nelem] );
1112 for(i = 0; i < N_ISO_TE_RECOMB; i++)
1113 {
1114 for( i1 = ( (nelem == ipISO) ? (RREC_MAXN + 1) : (LIKE_RREC_MAXN( nelem ) + 1) ); i1< NHYDRO_MAX_LEVEL; i1++ )
1115 {
1116 TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
1117 }
1118 for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
1119 {
1120 TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
1121 }
1122 fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[ipISO][nelem][i] ) );
1123 }
1124 fprintf(ioRECOMB, "\n" );
1125 }
1126 /* end the file with the same information */
1127 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",
1128 RECOMBMAGIC ,
1129 NumLevRecomb[ipISO][ipISO],
1130 NumLevRecomb[ipISO][ipISO+1],
1132 iso_ctrl.chISO[ipISO],
1133 NumLevRecomb[ipISO][ipISO],
1134 elementnames.chElementSym[ipISO],
1135 NumLevRecomb[ipISO][ipISO+1],
1137
1138 fclose( ioRECOMB );
1139
1140 fprintf( ioQQQ, "iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
1141 fprintf( ioQQQ, "The compilation is completed successfully.\n");
1143 }
1144 }
1145
1146 return;
1147}
1148
1149double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo )
1150{
1151 double rate;
1152 long ipTe, i;
1153 double TeDRCoef[NUM_DR_TEMPS];
1154 const freeBound *fb = &iso_sp[ipISO][nelem].fb[ipLo];
1155 const double Te_over_Z1_Squared[NUM_DR_TEMPS] = {
1156 1.00000, 1.30103, 1.69897, 2.00000, 2.30103, 2.69897, 3.00000,
1157 3.30103, 3.69897, 4.00000, 4.30103, 4.69897, 5.00000, 5.30103,
1158 5.69897, 6.00000, 6.30103, 6.69897, 7.00000 };
1159
1160 DEBUG_ENTRY( "iso_dielec_recomb_rate()" );
1161
1162 /* currently only two iso sequences and only he-like is applicable. */
1163 ASSERT( ipISO == ipHE_LIKE );
1164 ASSERT( ipLo >= 0 );
1165
1166 /* temperature grid is nelem^2 * constant temperature grid above. */
1167 for( i=0; i<NUM_DR_TEMPS; i++ )
1168 {
1169 TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (double) nelem );
1170 }
1171
1172 if( ipLo == ipHe1s1S )
1173 {
1174 rate = 0.;
1175 }
1176 else if( ipLo<iso_sp[ipISO][nelem].numLevels_max )
1177 {
1178 if( phycon.alogte <= TeDRCoef[0] )
1179 {
1180 /* Take lowest tabulated value for low temperature end. */
1181 rate = fb->DielecRecombVsTemp[0];
1182 }
1183 else if( phycon.alogte >= TeDRCoef[NUM_DR_TEMPS-1] )
1184 {
1185 /* use T^-1.5 extrapolation at high temperatures. */
1186 rate = fb->DielecRecombVsTemp[NUM_DR_TEMPS-1] *
1187 pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] - phycon.alogte ) ) ;
1188 }
1189 else
1190 {
1191 /* find temperature in tabulated values. */
1192 ipTe = hunt_bisect( TeDRCoef, NUM_DR_TEMPS, phycon.alogte );
1193
1194 ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1) );
1195
1196 if( fb->DielecRecombVsTemp[ipTe+1] == 0. )
1197 rate = 0.;
1198 else if( fb->DielecRecombVsTemp[ipTe] == 0. )
1199 rate = fb->DielecRecombVsTemp[ipTe+1];
1200 else
1201 {
1202 /* interpolate between tabulated points */
1203 rate = log10( fb->DielecRecombVsTemp[ipTe]) +
1204 (phycon.alogte-TeDRCoef[ipTe])*
1205 (log10(fb->DielecRecombVsTemp[ipTe+1])-log10(fb->DielecRecombVsTemp[ipTe]))/
1206 (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]);
1207
1208 rate = pow( 10., rate );
1209 }
1210 }
1211 }
1212 else
1213 {
1214 rate = 0.;
1215 }
1216
1217 ASSERT( rate >= 0. && rate < 1.0e-12 );
1218
1219 return rate*iso_ctrl.lgDielRecom[ipISO];
1220}
1221
1222/* TempInterp - interpolate on an array */
1224STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements )
1225{
1226 static long int ipTe=-1;
1227 double rate = 0.;
1228 long i0;
1229
1230 DEBUG_ENTRY( "TempInterp()" );
1231
1232 ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 );
1233
1234 if( ipTe < 0 )
1235 {
1236 /* te totally unknown */
1237 if( ( phycon.alogte < TempArray[0] ) ||
1238 ( phycon.alogte > TempArray[NumElements-1] ) )
1239 {
1240 fprintf(ioQQQ," TempInterp called with te out of bounds \n");
1242 }
1243 ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );
1244 }
1245 else if( phycon.alogte < TempArray[ipTe] )
1246 {
1247 /* temp is too low, must also lower ipTe */
1248 ASSERT( phycon.alogte > TempArray[0] );
1249 /* decrement ipTe until it is correct */
1250 while( ( phycon.alogte < TempArray[ipTe] ) && ipTe > 0)
1251 --ipTe;
1252 }
1253 else if( phycon.alogte > TempArray[ipTe+1] )
1254 {
1255 /* temp is too high */
1256 ASSERT( phycon.alogte <= TempArray[NumElements-1] );
1257 /* increment ipTe until it is correct */
1258 while( ( phycon.alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1)
1259 ++ipTe;
1260 }
1261
1262 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
1263
1264 /* ipTe should now be valid */
1265 ASSERT( ( phycon.alogte >= TempArray[ipTe] )
1266 && ( phycon.alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );
1267
1268 if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. )
1269 {
1270 rate = 0.;
1271 }
1272 else
1273 {
1274 /* Do a four-point interpolation */
1275 const int ORDER = 3; /* order of the fitting polynomial */
1276 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
1277 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
1278 }
1279
1280 return rate;
1281}
static double x1[83]
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define PrintEfmt(F, V)
Definition cddefines.h:1472
const int ipRecEsc
Definition cddefines.h:279
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int NHYDRO_MAX_LEVEL
Definition cddefines.h:266
double qg32(double, double, double(*)(double))
Definition service.cpp:1053
const int LIMELM
Definition cddefines.h:258
#define EXIT_SUCCESS
Definition cddefines.h:138
int dprintf(FILE *fp, const char *format,...)
Definition service.cpp:1009
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
const int ipRecNetEsc
Definition cddefines.h:281
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipRecRad
Definition cddefines.h:283
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long min(int a, long b)
Definition cddefines.h:723
void ShowMe(void)
Definition service.cpp:181
static t_ADfA & Inst()
Definition cddefines.h:175
double DielecRecombVsTemp[NUM_DR_TEMPS]
Definition freebound.h:35
double H_rad_rec(long int iz, long int n, double t)
t_conv conv
Definition conv.cpp:5
t_cosmology cosmology
Definition cosmology.cpp:11
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
@ AS_LOCAL_ONLY
Definition cpu.h:208
const double SMALLDOUBLE
Definition cpu.h:195
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
#define NUM_DR_TEMPS
Definition freebound.h:7
#define chLine_LENGTH
double He_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long nelem)
double Recomb_Seaton59(long nelem, double temp, long n)
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
double H_cross_section(double EgammaRyd, double EthRyd, long n, long l, long nelem)
static double kTRyd
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
#define IPRAD
Definition iso.h:86
const int ipH1s
Definition iso.h:27
#define SumUpToThisN
Definition iso.h:104
#define LIKE_RREC_MAXN(A_)
Definition iso.h:98
#define N_(A_)
Definition iso.h:20
#define N_ISO_TE_RECOMB
Definition iso.h:100
const int ipHE_LIKE
Definition iso.h:63
#define S_(A_)
Definition iso.h:22
#define RREC_MAXN
Definition iso.h:95
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
#define RECOMBMAGIC
Definition iso.h:106
const int ipH_LIKE
Definition iso.h:62
#define L_(A_)
Definition iso.h:21
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
static long int globalL
static long int globalISO
void iso_recomb_malloc(void)
static long int globalZ
static double TeRRCoef[N_ISO_TE_RECOMB]
static long int globalS
STATIC double iso_recomb_integrand(double EE)
void iso_recomb_setup(long ipISO)
static long int globalN
static double **** RRCoef
STATIC double TempInterp(double *TempArray, double *ValueArray, long NumElements)
double iso_RRCoef_Te(long ipISO, long nelem, long n)
static double *** TotalRecomb
static long ** NumLevRecomb
double iso_recomb_check(long ipISO, long nelem, long level, double temperature)
double iso_dielec_recomb_rate(long ipISO, long nelem, long ipLo)
STATIC void iso_put_recomb_error(long ipISO, long nelem)
void iso_recomb_auxiliary_free(void)
double iso_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long globalZ, long globalISO)
static double global_EthRyd
void iso_radiative_recomb_effective(long ipISO, long nelem)
void iso_radiative_recomb(long ipISO, long nelem)
static realnum * wavelength
t_opac opac
Definition opacity.cpp:5
#define S(I_, J_)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double MILNE_CONST
Definition physconst.h:233
UNUSED const double EE
Definition physconst.h:23
UNUSED const double TE1RYD
Definition physconst.h:183
UNUSED const double RYDLAM
Definition physconst.h:176
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
double RT_recom_effic(long int ip)
t_save save
Definition save.cpp:5
static double * ex
Definition species2.cpp:28
void TempChange(double TempNew, bool lgForceUpdate)
double lagrange(const double x[], const double y[], long n, double xval)
long hunt_bisect(const T x[], long n, T xval)
Definition thirdparty.h:270
t_trace trace
Definition trace.cpp:5