cloudy  trunk
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] */
33 static double ****RRCoef/*[LIMELM][NumLevRecomb[ipISO][nelem]][N_ISO_TE_RECOMB]*/;
34 static long **NumLevRecomb;
35 static double ***TotalRecomb; /*[ipISO][nelem][i]*/
36 
37 /* the array of logs of temperatures at which RRCoef was defined */
38 static double TeRRCoef[N_ISO_TE_RECOMB];
39 static double kTRyd,global_EthRyd;
40 static long int globalZ,globalISO;
41 static long int globalN, globalL, globalS;
42 
43 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements );
44 STATIC double iso_recomb_integrand(double EE);
45 STATIC void iso_put_recomb_error( long ipISO, long nelem );
46 
47 double 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
84  TotalInsanity();
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. */
120 STATIC 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);
126  temp = iso_cross_section( ERyd , global_EthRyd, globalN, globalL, globalS, globalZ, globalISO );
127  x1 *= temp;
128 
129  return x1;
130 }
131 
132 double 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
142  TotalInsanity();
143 
144  return cross_section;
145 }
146 
147 /*=======================================================*/
148 /* iso_radiative_recomb get rad recomb rate coefficients for iso sequences */
149 void iso_radiative_recomb(
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();
538  cdEXIT(EXIT_FAILURE);
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 
560 STATIC 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 
602 void 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 );
691  prt_wl( ioQQQ, wavelength );
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 */
708 double 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 */
734 double 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 */
765 void iso_recomb_malloc( void )
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 
825 void iso_recomb_auxiliary_free( void )
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 
838 void 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]);
932  cdEXIT(EXIT_FAILURE);
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" ,
945  RECOMBMAGIC ,
946  NumLevRecomb[ipISO][ipISO],
947  NumLevRecomb[ipISO][ipISO+1],
948  N_ISO_TE_RECOMB,
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" );
953  cdEXIT(EXIT_FAILURE);
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] );
968  cdEXIT(EXIT_FAILURE);
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" );
980  cdEXIT(EXIT_FAILURE);
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" );
1004  cdEXIT(EXIT_FAILURE);
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]);
1029  cdEXIT(EXIT_FAILURE);
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],
1046  N_ISO_TE_RECOMB,
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" );
1051  cdEXIT(EXIT_FAILURE);
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],
1075  N_ISO_TE_RECOMB,
1076  iso_ctrl.chISO[ipISO],
1077  NumLevRecomb[ipISO][ipISO],
1078  elementnames.chElementSym[ipISO],
1079  NumLevRecomb[ipISO][ipISO+1],
1080  N_ISO_TE_RECOMB );
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],
1131  N_ISO_TE_RECOMB,
1132  iso_ctrl.chISO[ipISO],
1133  NumLevRecomb[ipISO][ipISO],
1134  elementnames.chElementSym[ipISO],
1135  NumLevRecomb[ipISO][ipISO+1],
1136  N_ISO_TE_RECOMB );
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");
1142  cdEXIT(EXIT_SUCCESS);
1143  }
1144  }
1145 
1146  return;
1147 }
1148 
1149 double 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 */
1224 STATIC 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");
1241  cdEXIT(EXIT_FAILURE);
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 }
thermal.h
chLine_LENGTH
#define chLine_LENGTH
freeBound::DielecRecombVsTemp
double DielecRecombVsTemp[NUM_DR_TEMPS]
Definition: freebound.h:35
t_ADfA::H_rad_rec
double H_rad_rec(long int iz, long int n, double t)
Definition: atmdat_adfa.cpp:712
t_trace::lgIsoTraceFull
bool lgIsoTraceFull[NISO]
Definition: trace.h:88
t_iso_sp::HighestLevelOpacStack
vector< double > HighestLevelOpacStack
Definition: iso.h:588
helike.h
globalS
static long int globalS
Definition: iso_radiative_recomb.cpp:41
t_trace::lgHeBug
bool lgHeBug
Definition: trace.h:82
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
iso_put_error
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
SumUpToThisN
#define SumUpToThisN
Definition: iso.h:104
SMALLDOUBLE
const double SMALLDOUBLE
Definition: cpu.h:195
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
dense
t_dense dense
Definition: dense.cpp:24
t_opac::lgCaseB
bool lgCaseB
Definition: opacity.h:161
NHYDRO_MAX_LEVEL
const int NHYDRO_MAX_LEVEL
Definition: cddefines.h:266
elementnames.h
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:175
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:51
TeRRCoef
static double TeRRCoef[N_ISO_TE_RECOMB]
Definition: iso_radiative_recomb.cpp:38
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:281
t_isoCTRL::chISO
const char * chISO[NISO]
Definition: iso.h:330
FFmtRead
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition: service.cpp:381
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
t_isoCTRL::lgTopoff
bool lgTopoff[NISO]
Definition: iso.h:408
t_iso_sp::RadRec_effec
double RadRec_effec
Definition: iso.h:517
iso_put_recomb_error
STATIC void iso_put_recomb_error(long ipISO, long nelem)
Definition: iso_radiative_recomb.cpp:560
AS_LOCAL_ONLY
@ AS_LOCAL_ONLY
Definition: cpu.h:208
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
L_
#define L_(A_)
Definition: iso.h:21
thirdparty.h
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_isoCTRL::lgRandErrGen
bool lgRandErrGen[NISO]
Definition: iso.h:403
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
iso_recomb_integrand
STATIC double iso_recomb_integrand(double EE)
Definition: iso_radiative_recomb.cpp:120
cloudy_exit
Definition: cddefines.h:396
t_iso_sp::CascadeProb
multi_arr< double, 2 > CascadeProb
Definition: iso.h:450
ex
static double * ex
Definition: species2.cpp:28
iso_dielec_recomb_rate
double iso_dielec_recomb_rate(long ipISO, long nelem, long ipLo)
Definition: iso_radiative_recomb.cpp:1149
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:82
dynamics.h
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
freeBound
Definition: freebound.h:9
Recomb_Seaton59
double Recomb_Seaton59(long nelem, double temp, long n)
Definition: helike_recom.cpp:194
opac
t_opac opac
Definition: opacity.cpp:5
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
atmdat.h
t_ionbal::DR_Badnell_rate_coef
double ** DR_Badnell_rate_coef
Definition: ionbal.h:205
x1
static double x1[83]
Definition: atmdat_3body.cpp:28
MIN2
#define MIN2
Definition: cddefines.h:761
iso_RRCoef_Te
double iso_RRCoef_Te(long ipISO, long nelem, long n)
Definition: iso_radiative_recomb.cpp:708
t_isoCTRL::lgNoRecombInterp
bool lgNoRecombInterp[NISO]
Definition: iso.h:385
hunt_bisect
long hunt_bisect(const T x[], long n, T xval)
Definition: thirdparty.h:270
iso_radiative_recomb_effective
void iso_radiative_recomb_effective(long ipISO, long nelem)
Definition: iso_radiative_recomb.cpp:602
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
RT_recom_effic
double RT_recom_effic(long int ip)
Definition: rt_recom_effic.cpp:11
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
globalL
static long int globalL
Definition: iso_radiative_recomb.cpp:41
iso_recomb_check
double iso_recomb_check(long ipISO, long nelem, long level, double temperature)
Definition: iso_radiative_recomb.cpp:734
t_isoCTRL::lgCompileRecomb
bool lgCompileRecomb[NISO]
Definition: iso.h:380
H_cross_section
double H_cross_section(double EgammaRyd, double EthRyd, long n, long l, long nelem)
Definition: hydro_recom.cpp:22
dense.h
globalN
static long int globalN
Definition: iso_radiative_recomb.cpp:41
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
trace
t_trace trace
Definition: trace.cpp:5
iso_recomb_malloc
void iso_recomb_malloc(void)
Definition: iso_radiative_recomb.cpp:765
cddefines.h
kTRyd
static double kTRyd
Definition: iso_radiative_recomb.cpp:39
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
NUM_DR_TEMPS
#define NUM_DR_TEMPS
Definition: freebound.h:7
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
lagrange
double lagrange(const double x[], const double y[], long n, double xval)
Definition: thirdparty_interpolate.cpp:433
N_ISO_TE_RECOMB
#define N_ISO_TE_RECOMB
Definition: iso.h:100
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
RREC_MAXN
#define RREC_MAXN
Definition: iso.h:95
t_opac::OpacStack
double * OpacStack
Definition: opacity.h:151
MILNE_CONST
const UNUSED double MILNE_CONST
Definition: physconst.h:233
IPRAD
#define IPRAD
Definition: iso.h:86
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
LIKE_RREC_MAXN
#define LIKE_RREC_MAXN(A_)
Definition: iso.h:98
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
LIMELM
const int LIMELM
Definition: cddefines.h:258
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
iso_radrecomb_from_cross_section
double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
Definition: iso_radiative_recomb.cpp:47
t_cosmology::lgDo
bool lgDo
Definition: cosmology.h:44
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
ipRecRad
const int ipRecRad
Definition: cddefines.h:283
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:487
save.h
cosmology
t_cosmology cosmology
Definition: cosmology.cpp:11
iso_radiative_recomb
void iso_radiative_recomb(long ipISO, long nelem)
Definition: iso_radiative_recomb.cpp:149
iso_recomb_setup
void iso_recomb_setup(long ipISO)
Definition: iso_radiative_recomb.cpp:838
cosmology.h
S
#define S(I_, J_)
Definition: optimize_subplx.cpp:1835
prt.h
hydrogenic.h
He_cross_section
double He_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long nelem)
Definition: helike_recom.cpp:44
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:513
iso_recomb_auxiliary_free
void iso_recomb_auxiliary_free(void)
Definition: iso_radiative_recomb.cpp:825
t_save::ioRecom
FILE * ioRecom
Definition: save.h:345
rt.h
RECOMBMAGIC
#define RECOMBMAGIC
Definition: iso.h:106
EE
const UNUSED double EE
Definition: physconst.h:23
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
t_save::lgioRecom
bool lgioRecom
Definition: save.h:346
ionbal.h
t_iso_sp::lgMustReeval
bool lgMustReeval
Definition: iso.h:481
globalISO
static long int globalISO
Definition: iso_radiative_recomb.cpp:40
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
TempInterp
STATIC double TempInterp(double *TempArray, double *ValueArray, long NumElements)
Definition: iso_radiative_recomb.cpp:1224
min
long min(int a, long b)
Definition: cddefines.h:723
physconst.h
t_iso_sp::BranchRatio
multi_arr< double, 2 > BranchRatio
Definition: iso.h:451
conv
t_conv conv
Definition: conv.cpp:5
t_trace::ipIsoTrace
long int ipIsoTrace[NISO]
Definition: trace.h:91
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1472
TotalRecomb
static double *** TotalRecomb
Definition: iso_radiative_recomb.cpp:35
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
t_phycon::alogte
double alogte
Definition: phycon.h:82
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
read_whole_line
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition: service.cpp:70
taulines.h
NumLevRecomb
static long ** NumLevRecomb
Definition: iso_radiative_recomb.cpp:34
dprintf
int dprintf(FILE *fp, const char *format,...)
Definition: service.cpp:1009
cross_section
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
Definition: helike_recom.cpp:69
phycon.h
t_iso_sp::n_HighestResolved_local
long int n_HighestResolved_local
Definition: iso.h:507
t_trace::lgHBug
bool lgHBug
Definition: trace.h:85
ShowMe
void ShowMe(void)
Definition: service.cpp:181
ipH1s
const int ipH1s
Definition: iso.h:27
globalZ
static long int globalZ
Definition: iso_radiative_recomb.cpp:40
TE1RYD
const UNUSED double TE1RYD
Definition: physconst.h:183
global_EthRyd
static double global_EthRyd
Definition: iso_radiative_recomb.cpp:39
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
qg32
double qg32(double, double, double(*)(double))
Definition: service.cpp:1053
S_
#define S_(A_)
Definition: iso.h:22
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:279
RRCoef
static double **** RRCoef
Definition: iso_radiative_recomb.cpp:33
iso_get_total_num_levels
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
Definition: iso_create.cpp:1465
EN1RYD
const UNUSED double EN1RYD
Definition: physconst.h:179
iso_cross_section
double iso_cross_section(double EgammaRyd, double EthRyd, long n, long l, long S, long globalZ, long globalISO)
Definition: iso_radiative_recomb.cpp:132
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
t_iso_sp::QuantumNumbers2Index
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:461
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
N_
#define N_(A_)
Definition: iso.h:20
helike_recom.h
wavelength
static realnum * wavelength
Definition: monitor_results.cpp:70
Generated by   doxygen 1.8.17