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