cloudy  trunk
 All Data Structures 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-2008 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.h"
10 #include "conv.h"
11 #include "dense.h"
12 #include "elementnames.h"
13 #include "helike.h"
14 #include "helike_recom.h"
15 #include "hydrogenic.h"
16 #include "ionbal.h"
17 #include "iso.h"
18 #include "opacity.h"
19 #include "phycon.h"
20 #include "physconst.h"
21 #include "prt.h"
22 #include "punch.h"
23 #include "thermal.h"
24 #include "thirdparty.h"
25 #include "trace.h"
26 #include "rt.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 
37 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements );
38 
39 double iso_radrecomb_from_cross_section(long ipISO, double temp, long nelem, long ipLo)
40 {
41  double alpha;
42 
43  DEBUG_ENTRY( "iso_radrecomb_from_cross_section()" );
44 
45  if( ipISO == ipH_LIKE )
46  alpha = hlike_radrecomb_from_cross_section( temp, nelem, ipLo);
47  else if( ipISO == ipHE_LIKE )
48  alpha = helike_radrecomb_from_cross_section( temp, nelem, ipLo);
49  else
50  TotalInsanity();
51 
52  return alpha;
53 }
54 
55 /*=======================================================*/
56 /* iso_radiative_recomb get rad recomb rate coefficients for iso sequences */
58  long ipISO,
59  /* nelem on the c scale, He is 1 */
60  long nelem )
61 {
62  long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel;
63  double topoff, TotMinusExplicitResolved,
64  TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
65  double RecExplictLevels, TotalRadRecomb, RecCollapsed;
66  static double TeUsed[NISO][LIMELM]={
67  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
68  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
69  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
70  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
71  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
72  0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
73 
74  DEBUG_ENTRY( "iso_radiative_recomb()" );
75 
76  /* evaluate recombination escape probability for all levels */
77 
78  /* define radiative recombination rates for all levels */
79  /* this will be the sum of all levels explicitly included in the model atom */
80  RecExplictLevels = 0.;
81 
82  /* number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
83  ipFirstCollapsed = iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem];
84 
85  if( (fabs(1.-TeUsed[ipISO][nelem]/phycon.te)> 0.001) || hydro.lgReevalRecom || !conv.nTotalIoniz )
86  {
87  TeUsed[ipISO][nelem] = phycon.te;
88 
89  for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
90  {
91  /* this is radiative recombination rate coefficient */
92  double RadRec;
93 
94  if( !iso.lgNoRecombInterp[ipISO] )
95  {
96  RadRec = iso_RRCoef_Te( ipISO, nelem , ipLevel );
97  }
98  else
99  {
100  RadRec = iso_radrecomb_from_cross_section( ipISO, phycon.te, nelem, ipLevel);
101  }
102 
103  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec;
104 
105  if( iso.lgRandErrGen[ipISO] )
106  {
107  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *=
108  /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
109  * the error factors for rrc are always stored at iso.numLevels_max, regardless of
110  * continuum lowering effects. */
111  iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD];
112  }
113 
114  ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. );
115 
116  RecExplictLevels += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
117 
118  if( iso.lgDielRecom[ipISO] )
119  {
120  /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
121  iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
122  Total_DR_Added += iso.DielecRecomb[ipISO][nelem][ipLevel];
123  }
124  }
125 
126  /**************************************************/
127  /*** Add on recombination to collapsed levels ***/
128  /**************************************************/
129  RecCollapsed = 0.;
130  for( ipLevel=ipFirstCollapsed; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel )
131  {
132  /* use hydrogenic for collapsed levels */
133  double RadRec = t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, N_(ipLevel), phycon.te);
134 
135  /* this is radiative recombination rate coefficient */
136  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] = RadRec;
137 
138  /* This kills recombination into the collapsed level so that the forced
139  * statistical weighting can be bypassed */
140  if( !iso.lgTopoff[ipISO] )
141  {
142  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *= 1E-10;
143  }
144 
145  if( iso.lgRandErrGen[ipISO] )
146  {
147  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] *=
148  /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
149  * the error factors for rrc are always stored at iso.numLevels_max, regardless of
150  * continuum lowering effects. */
151  iso.ErrorFactor[ipISO][nelem][ iso.numLevels_max[ipISO][nelem] ][ipLevel][IPRAD];
152  }
153 
154  RecCollapsed += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
155 
156  ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] > 0. );
157 
158  if( iso.lgDielRecom[ipISO] )
159  {
160  /* \todo 2 suppress these rates for continuum lowering using factors from Jordan (1969). */
161  iso.DielecRecomb[ipISO][nelem][ipLevel] = iso_dielec_recomb_rate( ipISO, nelem, ipLevel );
162  Total_DR_Added += iso.DielecRecomb[ipISO][nelem][ipLevel];
163  }
164  }
165 
166  /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
167  for( ipLevel = 0; ipLevel<iso.numLevels_local[ipISO][nelem]; ipLevel++ )
168  {
169  if( N_(ipLevel) == ThisN )
170  {
171  TotRRThisN += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
172  }
173  else
174  {
175  ASSERT( (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (phycon.te>3E7) || (nelem!=ipHELIUM) );
176  LastN = ThisN;
177  ThisN = N_(ipLevel);
178  TotRRLastN = TotRRThisN;
179  TotRRThisN = iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad];
180 
181  {
182  /* Print the sum of recombination coefficients per n at current temp. */
183  /*@-redef@*/
184  enum {DEBUG_LOC=false};
185  /*@+redef@*/
186  static long RUNONCE = false;
187 
188  if( !RUNONCE && DEBUG_LOC )
189  {
190  static long FIRSTTIME = true;
191 
192  if( FIRSTTIME )
193  {
194  fprintf( ioQQQ,"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n",
195  ipISO, nelem, phycon.te);
196  FIRSTTIME= false;
197  }
198 
199  fprintf( ioQQQ,"%li\t%.2e\n",LastN,TotRRLastN );
200  }
201  RUNONCE = true;
202  }
203  }
204  }
205 
206  /* Get total recombination into all levels, including those not explicitly considered. */
207  if( iso.lgNoRecombInterp[ipISO] )
208  {
209  /* We are not interpolating, must calculate total now, as sum of resolved and collapsed levels... */
210  TotalRadRecomb = RecCollapsed + RecExplictLevels;
211 
212  /* Plus additional levels out to a predefined limit... */
213  for( long nLo = N_(ipFirstCollapsed) + iso.nCollapsed_max[ipISO][nelem]; nLo < NHYDRO_MAX_LEVEL; nLo++ )
214  {
215  TotalRadRecomb += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO, nLo, phycon.te);
216  }
217  /* Plus a bunch more for good measure */
218  for( long nLo = NHYDRO_MAX_LEVEL; nLo<=SumUpToThisN; nLo++ )
219  {
220  TotalRadRecomb += Recomb_Seaton59( nelem+1-ipISO, phycon.te, nLo );
221  }
222  }
223  else
224  {
225  /* We are interpolating, and total was calculated here in iso_recomb_setup */
226  TotalRadRecomb = iso_RRCoef_Te( ipISO, nelem, iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] );
227  }
228 
229  /* If generating random error, apply one to total recombination */
230  if( iso.lgRandErrGen[ipISO] )
231  {
232  /* this has to be from iso.numLevels_max instead of iso.numLevels_local because
233  * the error factors for rrc are always stored at iso.numLevels_max, regardless of
234  * continuum lowering effects. */
235  TotalRadRecomb *=
236  iso.ErrorFactor[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][iso.numLevels_max[ipISO][nelem]][IPRAD];
237  }
238 
239  /* this is case B recombination, sum without the ground included */
240  iso.RadRec_caseB[ipISO][nelem] = TotalRadRecomb - iso.RadRecomb[ipISO][nelem][0][ipRecRad];
241  ASSERT( iso.RadRec_caseB[ipISO][nelem] > 0.);
242 
243  /**********************************************************************/
244  /*** Add topoff (excess) recombination to top level. This is only ***/
245  /*** done if atom is not full size. ***/
246  /**********************************************************************/
247  if( !iso.lgLevelsLowered[ipISO][nelem] )
248  {
249  /* at this point we have RecExplictLevels, the sum of radiative recombination
250  * to all levels explicitly included in the model atom the total
251  * recombination rate. The difference is the amount of "topoff" we will need to do */
252  TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
253 
254  topoff = TotMinusExplicitResolved - RecCollapsed;
255 
256  /* the t_ADfA::Inst().rad_rec fits are too small at high temperatures, so this atom is
257  * better than the topoff. Only a problem for helium itself, at high temperatures.
258  * complain if the negative topoff is not for this case */
259  if( topoff < 0. && (nelem!=ipHELIUM || phycon.te < 1e5 ) &&
260  fabs( topoff/TotalRadRecomb ) > 0.01 )
261  {
262  fprintf(ioQQQ," PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",
263  nelem, topoff/TotalRadRecomb, phycon.te );
264  }
265 
266  if( !iso.lgTopoff[ipISO] )
267  topoff *= 1E-20;
268 
269  /* We always have at least one collapsed level if continuum is not lowered. Put topoff there. */
270  iso.RadRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1][ipRecRad] += MAX2(0., topoff );
271 
272  /* check for negative DR topoff, but only if Total_DR_Added is not negligible compared with TotalRadRecomb */
273  if( Total_DR_Added > TotalRadRecomb/100. )
274  {
275  if( Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 )
276  {
277  fprintf(ioQQQ," PROBLEM negative DR topoff for %li, rel error was %.2e, temp was %.2f\n",
278  nelem,
279  Total_DR_Added / ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0,
280  phycon.te );
281  }
282  }
283 
284  if( iso.lgDielRecom[ipISO] )
285  {
286  /* \todo 2 suppress this total rate for continuum lowering using factors from Jordan (1969). */
287  /* put extra DR in top level */
288  iso.DielecRecomb[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] =
289  ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added;
290  }
291  }
292 
293  }
294 
295  /**************************************************************/
296  /*** Stuff escape probabilities, and effective rad recomb ***/
297  /**************************************************************/
298 
299  /* total effective radiative recombination, initialize to zero */
300  iso.RadRec_effec[ipISO][nelem] = 0.;
301 
302  for( ipLevel=0; ipLevel<iso.numLevels_local[ipISO][nelem]; ++ipLevel )
303  {
304  /* option for case b conditions, kill ground state recombination */
305  if( opac.lgCaseB && ipLevel==0 )
306  {
307  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 1e-10;
308  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 1e-10;
309  }
310  else
311  {
312  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] =
313  RT_recom_effic(iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]);
314 
315  /* net escape prob includes dest by background opacity */
316  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] =
317  MIN2( 1.,
318  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] +
319  (1.-iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc])*
320  iso.ConOpacRatio[ipISO][nelem][ipLevel] );
321  }
322 
323  ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] >= 0. );
324  ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] >= 0. );
325  ASSERT( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] <= 1. );
326 
327  /* sum of all effective rad rec */
328  iso.RadRec_effec[ipISO][nelem] += iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]*
329  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc];
330  }
331 
332  /* zero out escape probabilities of levels above numLevels_local */
333  for( ipLevel=iso.numLevels_local[ipISO][nelem]; ipLevel<iso.numLevels_max[ipISO][nelem]; ++ipLevel )
334  {
335  /* this is escape probability */
336  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] = 0.;
337  /* net escape prob includes dest by background opacity */
338  iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc] = 0.;
339  }
340 
341  /* trace escape probabilities */
342  if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
343  {
344  fprintf( ioQQQ, " iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
345 
346  /* print continuum escape probability */
347  fprintf( ioQQQ, " iso_radiative_recomb recomb effic" );
348  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
349  {
350  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecEsc] ));
351  }
352  fprintf( ioQQQ, "\n" );
353 
354  /* net recombination efficiency factor, including background opacity*/
355  fprintf( ioQQQ, " iso_radiative_recomb recomb net effic" );
356  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
357  {
358  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecNetEsc]) );
359  }
360  fprintf( ioQQQ, "\n" );
361 
362  /* inward continuous optical depths */
363  fprintf( ioQQQ, " iso_radiative_recomb in optic dep" );
364  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
365  {
366  fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] ));
367  }
368  fprintf( ioQQQ, "\n" );
369 
370  /* outward continuous optical depths*/
371  fprintf( ioQQQ, " iso_radiative_recomb out op depth" );
372  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
373  {
374  fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipISO][nelem][ipLevel]-1] ));
375  }
376  fprintf( ioQQQ, "\n" );
377 
378  /* print radiative recombination coefficients */
379  fprintf( ioQQQ, " iso_radiative_recomb rad rec coef " );
380  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
381  {
382  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad]) );
383  }
384  fprintf( ioQQQ, "\n" );
385  }
386 
387  if( trace.lgTrace )
388  {
389  fprintf( ioQQQ, " iso_radiative_recomb total rec coef" );
390  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipISO][nelem] ));
391  fprintf( ioQQQ, " case A=" );
392  fprintf( ioQQQ,PrintEfmt("%10.3e",
393  iso.RadRec_caseB[ipISO][nelem] + iso.RadRecomb[ipISO][nelem][ipH1s][ipRecRad] ) );
394  fprintf( ioQQQ, " case B=");
395  fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipISO][nelem] ));
396  fprintf( ioQQQ, "\n" );
397  }
398 
399  /****************************/
400  /*** begin sanity check ***/
401  /****************************/
402  {
403  bool lgOK = true;
404  for( ipLevel=0; ipLevel < iso.numLevels_local[ipISO][nelem]; ipLevel++ )
405  {
406  if( iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] <= 0. )
407  {
408  fprintf( ioQQQ,
409  " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
410  ipISO, nelem, ipLevel, iso.RadRecomb[ipISO][nelem][ipLevel][ipRecRad] , phycon.te);
411  lgOK = false;
412  }
413  }
414  /* bail if we found problems */
415  if( !lgOK )
416  {
417  ShowMe();
418  cdEXIT(EXIT_FAILURE);
419  }
420  /*end sanity check */
421  }
422 
423  /* confirm that we have good rec coef at bottom and top of atom/ion */
424  ASSERT( iso.RadRecomb[ipISO][nelem][0][ipRecRad] > 0. );
425  ASSERT( iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecRad] > 0. );
426 
427  /* set true when punch recombination coefficients command entered */
428  if( punch.lgioRecom )
429  {
430  /* this prints Z on physical, not C, scale */
431  fprintf( punch.ioRecom, "%s %s %2li ",
432  iso.chISO[ipISO], elementnames.chElementSym[nelem], nelem+1 );
433  fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipISO][nelem] ));
434  fprintf( punch.ioRecom, "\n" );
435  }
436 
437  return;
438 }
439 
440 void iso_radiative_recomb_effective( long ipISO, long nelem )
441 {
442  DEBUG_ENTRY( "iso_radiative_recomb_effective()" );
443 
444  /* Find effective recombination coefficients */
445  for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
446  {
447  iso.RadEffec[ipISO][nelem][ipHi] = 0.;
448 
449  /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
450  for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ )
451  {
452  ASSERT( iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. );
453  ASSERT( iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad] >= 0. );
454 
455  iso.RadEffec[ipISO][nelem][ipHi] += iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] *
456  iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad];
457  }
458  }
459 
460  /**************************************************************/
461  /*** option to print effective recombination coefficients ***/
462  /**************************************************************/
463  {
464  enum {DEBUG_LOC=false};
465 
466  if( DEBUG_LOC )
467  {
468  const int maxPrt=10;
469 
470  fprintf( ioQQQ,"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem, phycon.te );
471  fprintf( ioQQQ, "N\tL\tS\tRadEffec\tLifetime\n" );
472 
473  for( long i=0; i<maxPrt; i++ )
474  {
475  fprintf( ioQQQ, "%li\t%li\t%li\t%e\t%e\n", N_(i), L_(i), S_(i),
476  iso.RadEffec[ipISO][nelem][i],
477  MAX2( 0., StatesElem[ipISO][nelem][i].lifetime ) );
478  }
479  fprintf( ioQQQ, "\n" );
480  }
481  }
482 
483  /* If we have the variable set, find errors in rad rates */
484  if( iso.lgRandErrGen[ipISO] )
485  {
486  dprintf( ioQQQ, "ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
487 
488  /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
489  for( long ipHi=0; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
490  {
491  iso.SigmaRadEffec[ipISO][nelem][ipHi] = 0.;
492 
493  /* >>chng 06 aug 17, from numLevels_max to numLevels_local */
494  for( long ipHigher=ipHi; ipHigher < iso.numLevels_local[ipISO][nelem]; ipHigher++ )
495  {
496  ASSERT( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] >= 0. );
497  ASSERT( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] >= 0. );
498 
499  /* iso.RadRecomb has to appear here because iso.Error is only relative error */
500  iso.SigmaRadEffec[ipISO][nelem][ipHi] += pow( iso.Error[ipISO][nelem][iso.numLevels_max[ipISO][nelem]][ipHigher][IPRAD] *
501  iso.CascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.) +
502  pow( iso.SigmaCascadeProb[ipISO][nelem][ipHigher][ipHi] * iso.RadRecomb[ipISO][nelem][ipHigher][ipRecRad], 2.);
503  }
504 
505  ASSERT( iso.SigmaRadEffec[ipISO][nelem][ipHi] >= 0. );
506  iso.SigmaRadEffec[ipISO][nelem][ipHi] = sqrt( iso.SigmaRadEffec[ipISO][nelem][ipHi] );
507 
508  for( long ipLo = 0; ipLo < ipHi; ipLo++ )
509  {
510  if( (( L_(ipLo) == L_(ipHi) + 1 ) || ( L_(ipHi) == L_(ipLo) + 1 )) )
511  {
512  double EnergyInRydbergs = iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - iso.xIsoLevNIonRyd[ipISO][nelem][ipHi];
513  realnum wavelength = (realnum)(RYDLAM/MAX2(1E-8,EnergyInRydbergs));
514  double emissivity = iso.RadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * EN1RYD * EnergyInRydbergs;
515  double sigma_emiss = 0., SigmaBranchRatio = 0.;
516 
517  if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (N_(ipHi)<=5) )
518  {
519  SigmaBranchRatio = iso.BranchRatio[ipISO][nelem][ipHi][ipLo] * sqrt(
520  pow( (double)iso.Error[ipISO][nelem][ipHi][ipLo][IPRAD], 2. ) +
521  pow( iso.SigmaAtot[ipISO][nelem][ipHi]*StatesElem[ipISO][nelem][ipHi].lifetime, 2.) );
522 
523  sigma_emiss = EN1RYD * EnergyInRydbergs * sqrt(
524  pow( (double)iso.SigmaRadEffec[ipISO][nelem][ipHi] * iso.BranchRatio[ipISO][nelem][ipHi][ipLo], 2.) +
525  pow( SigmaBranchRatio * iso.RadEffec[ipISO][nelem][ipHi], 2.) );
526 
527  /* \todo 2 make this a trace option */
528  dprintf( ioQQQ, "%li\t%li\t", ipHi, ipLo );
529  prt_wl( ioQQQ, wavelength );
530  fprintf( ioQQQ, "\t%e\t%e\t%e\t%e\t%e\t%e\n",
531  emissivity,
532  sigma_emiss,
533  iso.RadEffec[ipISO][nelem][ipHi],
534  iso.SigmaRadEffec[ipISO][nelem][ipHi],
535  iso.BranchRatio[ipISO][nelem][ipHi][ipLo],
536  SigmaBranchRatio);
537  }
538  }
539  }
540  }
541  }
542 
543  return;
544 }
545 /*iso_RRCoef_Te evaluated radiative recombination coef at some temperature */
546 double iso_RRCoef_Te( long ipISO, long nelem , long n )
547 {
548  double rate = 0.;
549 
550  DEBUG_ENTRY( "iso_RRCoef_Te()" );
551 
552  ASSERT( !iso.lgNoRecombInterp[ipISO] );
553 
554  /* if n is equal to the number of levels, return the total recomb, else recomb for given level. */
555  if( n == iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] )
556  {
557  rate = TempInterp( TeRRCoef, TotalRecomb[ipISO][nelem], N_ISO_TE_RECOMB );
558  }
559  else
560  {
561  rate = TempInterp( TeRRCoef, RRCoef[ipISO][nelem][n], N_ISO_TE_RECOMB );
562  }
563 
564  /* that was the log, now make linear */
565  rate = pow( 10. , rate );
566 
567  return rate;
568 }
569 
570 /*iso_recomb_check called by SanityCheck to confirm that recombination coef are ok
571  * return value is relative error between new calculation of recom, and interp value */
572 double iso_recomb_check( long ipISO, long nelem, long level, double temperature )
573 {
574  double RecombRelError ,
575  RecombInterp,
576  RecombCalc,
577  SaveTemp;
578 
579  DEBUG_ENTRY( "iso_recomb_check()" );
580 
581  /* first set temp to desired value */
582  SaveTemp = phycon.te;
583  TempChange(temperature);
584 
585  /* actually calculate the recombination coefficient from the Milne relation,
586  * normally only done due to compile he-like command */
587  RecombCalc = iso_radrecomb_from_cross_section( ipISO, temperature , nelem , level );
588 
589  /* interpolate the recombination coefficient, this is the usual method */
590  RecombInterp = iso_RRCoef_Te( ipISO, nelem , level );
591 
592  /* reset temp */
593  TempChange(SaveTemp);
594 
595  RecombRelError = ( RecombInterp - RecombCalc ) / MAX2( RecombInterp , RecombCalc );
596 
597  return RecombRelError;
598 }
599 
600 /* malloc space needed for iso recombination tables */
601 void iso_recomb_malloc( void )
602 {
603  DEBUG_ENTRY( "iso_recomb_malloc()" );
604 
605  NumLevRecomb = (long **)MALLOC(sizeof(long *)*(unsigned)NISO );
606  TotalRecomb = (double ***)MALLOC(sizeof(double **)*(unsigned)NISO );
607  RRCoef = (double ****)MALLOC(sizeof(double ***)*(unsigned)NISO );
608 
609  for( long ipISO=0; ipISO<NISO; ipISO++ )
610  {
611  TotalRecomb[ipISO] = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
612  RRCoef[ipISO] = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
613  /* The number of recombination coefficients to be read from file for each element. */
614  NumLevRecomb[ipISO] = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM );
615 
616  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
617  {
618  long int MaxLevels, maxN;
619 
620  TotalRecomb[ipISO][nelem] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
621 
622  if( nelem == ipISO )
623  maxN = RREC_MAXN;
624  else
625  maxN = LIKE_RREC_MAXN( nelem );
626 
627  NumLevRecomb[ipISO][nelem] = iso_get_total_num_levels( ipISO, maxN, 0 );
628 
629  if( nelem == ipISO || dense.lgElmtOn[nelem] )
630  {
631  /* must always have at least NumLevRecomb[ipISO][nelem] levels since that is number
632  * that will be read in from he rec data file, but possible to require more */
633  MaxLevels = MAX2( NumLevRecomb[ipISO][nelem] , iso.numLevels_max[ipISO][nelem] );
634 
635  /* always define this */
636  /* >>chng 02 jan 24, RRCoef will be iso.numLevels_max[ipISO][nelem] levels, not iso.numLevels_max,
637  * code will stop if more than this is requested */
638  RRCoef[ipISO][nelem] = (double**)MALLOC(sizeof(double*)*(unsigned)(MaxLevels) );
639 
640  for( long ipLo=0; ipLo < MaxLevels;++ipLo )
641  {
642  RRCoef[ipISO][nelem][ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)N_ISO_TE_RECOMB );
643  }
644  }
645  }
646  }
647 
648  for(long i = 0; i < N_ISO_TE_RECOMB; i++)
649  {
650  /* this is the vector of temperatures */
651  TeRRCoef[i] = 0.25*(i);
652  }
653 
654  /* >>chng 06 jun 06, NP, assert thrown at T == 1e10 K, just bump the
655  * high temperature end slightly. */
656  TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f;
657 
658  return;
659 }
660 
662 {
663  DEBUG_ENTRY( "iso_recomb_auxiliary_free()" );
664 
665  for( long i = 0; i< NISO; i++ )
666  {
667  free( NumLevRecomb[i] );
668  }
669  free( NumLevRecomb );
670 
671  return;
672 }
673 
674 void iso_recomb_setup( long ipISO )
675 {
676  double RadRecombReturn;
677  long int i, i1, i2, i3, i4, i5;
678  long int ipLo, nelem;
679 
680  const int chLine_LENGTH = 1000;
681  char chLine[chLine_LENGTH];
682  /* this must be longer than data path string, set in path.h/cpu.cpp */
683  const char* chFilename[NISO] =
684  { "h_iso_recomb.dat", "he_iso_recomb.dat" };
685 
686  FILE *ioDATA;
687  bool lgEOL;
688 
689  DEBUG_ENTRY( "iso_recomb_setup()" );
690 
691  /* if we are compiling the recombination data file, we must interpolate in temperature */
692  if( iso.lgCompileRecomb[ipISO] )
693  {
694  iso.lgNoRecombInterp[ipISO] = false;
695  }
696 
697  if( !iso.lgNoRecombInterp[ipISO] )
698  {
699  /******************************************************************/
701  /******************************************************************/
702  /* This flag says we are not compiling the data file */
703  if( !iso.lgCompileRecomb[ipISO] )
704  {
705  if( trace.lgTrace )
706  fprintf( ioQQQ," iso_recomb_setup opening %s:", chFilename[ipISO] );
707 
708  /* Now try to read from file...*/
709  try
710  {
711  ioDATA = open_data( chFilename[ipISO], "r" );
712  }
713  catch( cloudy_exit )
714  {
715  fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
716  fprintf( ioQQQ, " but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
717  ioDATA = NULL;
718  }
719  if( ioDATA == NULL )
720  {
721  /* Do on the fly computation of R.R. Coef's instead. */
722  for( nelem = ipISO; nelem < LIMELM; nelem++ )
723  {
724  if( dense.lgElmtOn[nelem] )
725  {
726  /* Zero out the recombination sum array. */
727  for(i = 0; i < N_ISO_TE_RECOMB; i++)
728  {
729  TotalRecomb[ipISO][nelem][i] = 0.;
730  }
731 
732  /* NumLevRecomb[ipISO][nelem] corresponds to n = 40 for H and He and 20 for ions, at present */
733  /* There is no need to fill in values for collapsed levels, because we do not need to
734  * interpolate for a given temperature, just calculate it directly with a hydrogenic routine. */
735  for( ipLo=0; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ )
736  {
737  /* loop over temperatures to produce array of recombination coefficients */
738  for(i = 0; i < N_ISO_TE_RECOMB; i++)
739  {
740  /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
741  RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
742  TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
743  RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
744  }
745  }
746  for(i = 0; i < N_ISO_TE_RECOMB; i++)
747  {
748  for( i1 = iso.n_HighestResolved_max[ipISO][nelem]+1; i1< NHYDRO_MAX_LEVEL; i1++ )
749  {
750  TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
751  }
752  for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
753  {
754  TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
755  }
756  TotalRecomb[ipISO][nelem][i] = log10( TotalRecomb[ipISO][nelem][i] );
757  }
758  }
759  }
760  }
761  /* Data file is present and readable...begin read. */
762  else
763  {
764  /* check that magic number is ok */
765  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
766  {
767  fprintf( ioQQQ, " iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
768  cdEXIT(EXIT_FAILURE);
769  }
770  i = 1;
771  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
772  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
773  i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
774  i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
775  if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
776  {
777  fprintf( ioQQQ,
778  " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
779  fprintf( ioQQQ,
780  " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
781  RECOMBMAGIC ,
782  NumLevRecomb[ipISO][ipISO],
783  NumLevRecomb[ipISO][ipISO+1],
785  i1 , i2 , i3, i4 );
786  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
787  fprintf( ioQQQ,
788  " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
789  cdEXIT(EXIT_FAILURE);
790  }
791 
792  i5 = 1;
793  /* now read in the data */
794  for( nelem = ipISO; nelem < LIMELM; nelem++ )
795  {
796  for( ipLo=0; ipLo <= NumLevRecomb[ipISO][nelem]; ipLo++ )
797  {
798  i5++;
799  /* get next line image */
800  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
801  {
802  fprintf( ioQQQ, " iso_recomb_setup could not read line %li of %s.\n", i5,
803  chFilename[ipISO] );
804  cdEXIT(EXIT_FAILURE);
805  }
806  /* each line starts with element and level number */
807  i3 = 1;
808  i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
809  i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
810  /* check that these numbers are correct */
811  if( i1!=nelem || i2!=ipLo )
812  {
813  fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
814  fprintf( ioQQQ,
815  " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
816  cdEXIT(EXIT_FAILURE);
817  }
818 
819  /* loop over temperatures to produce array of recombination coefficients */
820  for(i = 0; i < N_ISO_TE_RECOMB; i++)
821  {
822  double ThisCoef = FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
823 
824  if( nelem == ipISO || dense.lgElmtOn[nelem] )
825  {
826  /* The last line for each element is the total recombination for each temp. */
827  if( ipLo == NumLevRecomb[ipISO][nelem] )
828  {
829  TotalRecomb[ipISO][nelem][i] = ThisCoef;
830  }
831  else
832  RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
833  }
834 
835  if( lgEOL )
836  {
837  fprintf( ioQQQ, " iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
838  fprintf( ioQQQ,
839  " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
840  cdEXIT(EXIT_FAILURE);
841  }
842  }
843  }
844 
845  /* following loop only executed if we need more levels than are
846  * stored in the recom coef data set
847  * do not do collapsed levels since will use H-like recom there */
848  if( nelem == ipISO || dense.lgElmtOn[nelem] )
849  {
850  for( ipLo=NumLevRecomb[ipISO][nelem]; ipLo < iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem]; ipLo++ )
851  {
852  for(i = 0; i < N_ISO_TE_RECOMB; i++)
853  {
854  /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
855  RRCoef[ipISO][nelem][ipLo][i] = log10(iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo));
856  }
857  }
858  }
859  }
860 
861  /* check that ending magic number is ok */
862  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
863  {
864  fprintf( ioQQQ, " iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
865  cdEXIT(EXIT_FAILURE);
866  }
867  i = 1;
868  i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
869  i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
870  i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
871  i4 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
872 
873  if( i1 !=RECOMBMAGIC || i2 !=NumLevRecomb[ipISO][ipISO] || i3 !=NumLevRecomb[ipISO][ipISO+1] || i4 !=N_ISO_TE_RECOMB )
874  {
875  fprintf( ioQQQ,
876  " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
877  fprintf( ioQQQ,
878  " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
879  RECOMBMAGIC ,
880  NumLevRecomb[ipISO][ipISO],
881  NumLevRecomb[ipISO][ipISO+1],
883  i1 , i2 , i3, i4 );
884  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
885  fprintf( ioQQQ,
886  " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
887  cdEXIT(EXIT_FAILURE);
888  }
889 
890  /* close the data file */
891  fclose( ioDATA );
892  }
893  }
894  /* We are compiling the he_iso_recomb.dat data file. */
895  else if( iso.lgCompileRecomb[ipISO] )
896  {
897  /* option to create table of recombination coefficients,
898  * executed with the compile he-like command */
899  FILE *ioRECOMB;
900 
901  ASSERT( !iso.lgNoRecombInterp[ipISO] );
902 
903  /*RECOMBMAGIC the magic number for the table of recombination coefficients */
904  /*NumLevRecomb[ipISO][nelem] the number of levels that will be done */
905  /* create recombination coefficients */
906  ioRECOMB = open_data( chFilename[ipISO], "w", AS_LOCAL_ONLY );
907  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",
908  RECOMBMAGIC ,
909  NumLevRecomb[ipISO][ipISO],
910  NumLevRecomb[ipISO][ipISO+1],
912  iso.chISO[ipISO],
913  NumLevRecomb[ipISO][ipISO],
914  elementnames.chElementSym[ipISO],
915  NumLevRecomb[ipISO][ipISO+1],
916  N_ISO_TE_RECOMB );
917 
918  for( nelem = ipISO; nelem < LIMELM; nelem++ )
919  {
920  /* this must pass since compile xx-like command reset numLevels to the macro */
921  ASSERT( NumLevRecomb[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] );
922 
923  /* Zero out the recombination sum array. */
924  for(i = 0; i < N_ISO_TE_RECOMB; i++)
925  {
926  TotalRecomb[ipISO][nelem][i] = 0.;
927  }
928 
929  for( ipLo=ipHe1s1S; ipLo < NumLevRecomb[ipISO][nelem]; ipLo++ )
930  {
931  fprintf(ioRECOMB, "%li\t%li", nelem, ipLo );
932  /* loop over temperatures to produce array of recombination coefficients */
933  for(i = 0; i < N_ISO_TE_RECOMB; i++)
934  {
935  /* Store log of recombination coefficients, in N_ISO_TE_RECOMB half dec steps */
936  RadRecombReturn = iso_radrecomb_from_cross_section( ipISO, pow( 10.,TeRRCoef[i] ) ,nelem,ipLo);
937  TotalRecomb[ipISO][nelem][i] += RadRecombReturn;
938  RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
939  fprintf(ioRECOMB, "\t%f", RRCoef[ipISO][nelem][ipLo][i] );
940  }
941  fprintf(ioRECOMB, "\n" );
942  }
943 
944  /* Store one additional line in XX_iso_recomb.dat that gives the total recombination,
945  * as computed by the sum so far, plus levels up to NHYDRO_MAX_LEVEL using Verner's fits,
946  * plus levels up to SumUpToThisN using Seaton 59, for each element and each temperature. */
947  fprintf(ioRECOMB, "%li\t%li", nelem, NumLevRecomb[ipISO][nelem] );
948  for(i = 0; i < N_ISO_TE_RECOMB; i++)
949  {
950  for( i1 = ( (nelem == ipISO) ? (RREC_MAXN + 1) : (LIKE_RREC_MAXN( nelem ) + 1) ); i1< NHYDRO_MAX_LEVEL; i1++ )
951  {
952  TotalRecomb[ipISO][nelem][i] += t_ADfA::Inst().H_rad_rec(nelem+1-ipISO,i1, pow(10.,TeRRCoef[i]));
953  }
954  for( i1 = NHYDRO_MAX_LEVEL; i1<=SumUpToThisN; i1++ )
955  {
956  TotalRecomb[ipISO][nelem][i] += Recomb_Seaton59( nelem+1-ipISO, pow(10.,TeRRCoef[i]), i1 );
957  }
958  fprintf(ioRECOMB, "\t%f", log10( TotalRecomb[ipISO][nelem][i] ) );
959  }
960  fprintf(ioRECOMB, "\n" );
961  }
962  /* end the file with the same information */
963  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",
964  RECOMBMAGIC ,
965  NumLevRecomb[ipISO][ipISO],
966  NumLevRecomb[ipISO][ipISO+1],
968  iso.chISO[ipISO],
969  NumLevRecomb[ipISO][ipISO],
970  elementnames.chElementSym[ipISO],
971  NumLevRecomb[ipISO][ipISO+1],
972  N_ISO_TE_RECOMB );
973 
974  fclose( ioRECOMB );
975 
976  fprintf( ioQQQ, "iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
977  fprintf( ioQQQ, "The compilation is completed successfully.\n");
978  cdEXIT(EXIT_SUCCESS);
979  }
980  }
981 
982  return;
983 }
984 
985 double iso_dielec_recomb_rate( long ipISO, long nelem, long ipLo )
986 {
987  double rate;
988  long ipTe, i;
989  double TeDRCoef[NUM_DR_TEMPS];
990  const double Te_over_Z1_Squared[NUM_DR_TEMPS] = {
991  1.00000, 1.30103, 1.69897, 2.00000, 2.30103, 2.69897, 3.00000,
992  3.30103, 3.69897, 4.00000, 4.30103, 4.69897, 5.00000, 5.30103,
993  5.69897, 6.00000, 6.30103, 6.69897, 7.00000 };
994 
995  DEBUG_ENTRY( "iso_dielec_recomb_rate()" );
996 
997  /* currently only two iso sequences and only he-like is applicable. */
998  ASSERT( ipISO == ipHE_LIKE );
999  ASSERT( ipLo >= 0 );
1000 
1001  /* temperature grid is nelem^2 * constant temperature grid above. */
1002  for( i=0; i<NUM_DR_TEMPS; i++ )
1003  {
1004  TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (double) nelem );
1005  }
1006 
1007  if( ipLo == ipHe1s1S )
1008  {
1009  rate = 0.;
1010  }
1011  else if( ipLo<iso.numLevels_max[ipISO][nelem] )
1012  {
1013  if( phycon.alogte <= TeDRCoef[0] )
1014  {
1015  /* Take lowest tabulated value for low temperature end. */
1016  rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][0];
1017  }
1018  else if( phycon.alogte >= TeDRCoef[NUM_DR_TEMPS-1] )
1019  {
1020  /* use T^-1.5 extrapolation at high temperatures. */
1021  rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][NUM_DR_TEMPS-1] *
1022  pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] - phycon.alogte ) ) ;
1023  }
1024  else
1025  {
1026  /* find temperature in tabulated values. */
1027  ipTe = hunt_bisect( TeDRCoef, NUM_DR_TEMPS, phycon.alogte );
1028 
1029  ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1) );
1030 
1031  if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1] == 0. )
1032  rate = 0.;
1033  else if( iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe] == 0. )
1034  rate = iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1];
1035  else
1036  {
1037  /* interpolate between tabulated points */
1038  rate = log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]) +
1039  (phycon.alogte-TeDRCoef[ipTe])*
1040  (log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe+1])-log10(iso.DielecRecombVsTemp[ipISO][nelem][ipLo][ipTe]))/
1041  (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]);
1042 
1043  rate = pow( 10., rate );
1044  }
1045  }
1046  }
1047  else
1048  {
1049  rate = 0.;
1050  }
1051 
1052  ASSERT( rate >= 0. && rate < 1.0e-12 );
1053 
1054  return rate*iso.lgDielRecom[ipISO];
1055 }
1056 
1057 /* TempInterp - interpolate on an array */
1059 STATIC double TempInterp( double* TempArray , double* ValueArray, long NumElements )
1060 {
1061  static long int ipTe=-1;
1062  double rate = 0.;
1063  long i0;
1064 
1065  DEBUG_ENTRY( "TempInterp()" );
1066 
1067  ASSERT( fabs( 1. - (double)phycon.alogte/log10((double)phycon.te) ) < 0.0001 );
1068 
1069  if( ipTe < 0 )
1070  {
1071  /* te totally unknown */
1072  if( ( phycon.alogte < TempArray[0] ) ||
1073  ( phycon.alogte > TempArray[NumElements-1] ) )
1074  {
1075  fprintf(ioQQQ," TempInterp called with te out of bounds \n");
1076  cdEXIT(EXIT_FAILURE);
1077  }
1078  ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );
1079  }
1080  else if( phycon.alogte < TempArray[ipTe] )
1081  {
1082  /* temp is too low, must also lower ipTe */
1083  ASSERT( phycon.alogte > TempArray[0] );
1084  /* decrement ipTe until it is correct */
1085  while( ( phycon.alogte < TempArray[ipTe] ) && ipTe > 0)
1086  --ipTe;
1087  }
1088  else if( phycon.alogte > TempArray[ipTe+1] )
1089  {
1090  /* temp is too high */
1091  ASSERT( phycon.alogte <= TempArray[NumElements-1] );
1092  /* increment ipTe until it is correct */
1093  while( ( phycon.alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1)
1094  ++ipTe;
1095  }
1096 
1097  ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
1098 
1099  /* ipTe should now be valid */
1100  ASSERT( ( phycon.alogte >= TempArray[ipTe] )
1101  && ( phycon.alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );
1102 
1103  if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. )
1104  {
1105  rate = 0.;
1106  }
1107  else
1108  {
1109  /* Do a four-point interpolation */
1110  const int ORDER = 3; /* order of the fitting polynomial */
1111  i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
1112  rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
1113  }
1114 
1115  return rate;
1116 }

Generated for cloudy by doxygen 1.8.3.1