cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_hyperfine.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 /* HyperfineCreat establish space for hf arrays, reads atomic data from hyperfine.dat */
4 /* HyperfineCS - returns collision strengths for hyperfine struc transitions */
5 /*H21cm computes rate for H 21 cm from upper to lower excitation by atomic hydrogen */
6 /*h21_t_ge_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */
7 /*h21_t_lt_20 compute rate for H 21 cm from upper to lower excitation by atomic hydrogen */
8 /*H21cm_electron compute H 21 cm rate from upper to lower excitation by electrons - call by CoolEvaluate */
9 /*H21cm_H_atom - evaluate H atom spin changing collision rate, called by CoolEvaluate */
10 /*H21cm_proton - evaluate proton spin changing H atom collision rate, */
11 #include "cddefines.h"
12 #include "lines_service.h"
13 #include "phycon.h"
14 #include "dense.h"
15 #include "rfield.h"
16 #include "taulines.h"
17 #include "iso.h"
18 #include "trace.h"
19 #include "hyperfine.h"
20 #include "physconst.h"
21 
22 /* H21_cm_pops - fine level populations for 21 cm with Lya pumping included
23  * called in CoolEvaluate */
24 void H21_cm_pops( void )
25 {
26  /*atom_level2( HFLines[0] );*/
27  /*return;*/
28  /*
29  things we know on entry to this routine:
30  total population of 2p: StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop
31  total population of 1s: StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop
32  continuum pumping rate (lo-up) inside 21 cm line: HFLines[0].pump
33  upper to lower collision rate inside 21 cm line: HFLines[0].cs*dense.cdsqte
34  occupation number inside Lya: OccupationNumberLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] )
35 
36  level populations (cm-3) must be computed:
37  population of upper level of 21cm: HFLines[0].Hi->Pop
38  population of lower level of 21cm: HFLines[0].Lo->Pop
39  stimulated emission corrected population of lower level: HFLines[0].Emis->PopOpc
40  */
41 
42  double x,
43  PopTot;
44  double a32,a31,a41,a42,a21, occnu_lya ,
45  rate12 , rate21 , pump12 , pump21 , coll12 , coll21,
46  texc , occnu_lya_23 , occnu_lya_13,occnu_lya_24,occnu_lya_14, texc1, texc2;
47  a31 = 2.08e8; /* Einstein co-efficient for transition 1p1/2 to 0s1/2 */
48  a32 = 4.16e8; /* Einstein co-efficient for transition 1p1/2 to 1s1/2 */
49  a41 = 4.16e8; /* Einstein co-efficient for transition 1p3/2 to 0s1/2 */
50  a42 = 2.08e8; /* Einstein co-efficient for transition 1p3/2 to 1s1/2 */
51  /* These A values are determined from eqn. 17.64 of "The theory of Atomic structure
52  * and Spectra" by R. D. Cowan
53  * A hyperfine level has degeneracy Gf=(2F + 1)
54  * a2p1s = 6.24e8; Einstein co-efficient for transition 2p to 1s */
55  a21 = 2.85e-15; /* Einstein co-efficient for transition 1s1/2 to 0s1/2 */
56 
57  /* above is spontaneous rate - the net rate is this times escape and destruction
58  * probabilities */
59  a21 *= (HFLines[0].Emis->Pdest + HFLines[0].Emis->Pesc + HFLines[0].Emis->Pelec_esc);
60 
61  /* >>chng 04 dec 01, add hyperfine.lgLya_pump_21cm, option to turn off Lya pump
62  * of 21 cm, with no 21cm lya pump */
65 
66  /* >>chng 05 apr 20, GS, the lya occupation number for the hyperfine levels 0S1/2 and 1S1/2 are different*/
67  texc = TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
68  /* >>chng 05 apr 21, GS, Energy difference between 2p1/2 and 2p3/2 is taken from NSRDS */
69  if( texc > 0. )
70  {
71  /* convert to boltz factor, which will applied to occupation number of hiher energy transition */
72  texc1 = sexp(0.068/texc);
73  texc2 = sexp(((82259.272-82258.907)*T1CM)/texc);
74  }
75  else
76  {
77  texc1 = 0.;
78  texc2 = 0.;
79  }
80 
81  /* the continuum within Lya seen by the two levels is not exactly the same brightness. They
82  * differ by the exp when Lya is on Wein tail of black body, which must be true if 21 cm is important */
83 
84  occnu_lya_23 = occnu_lya;
85  occnu_lya_13 = occnu_lya*texc1;
86  occnu_lya_14 = occnu_lya_13*texc2;
87  occnu_lya_24 = occnu_lya*texc2;
88 
89  /* this is the 21 cm upward continuum pumping rate [s-1] for the attenuated incident and
90  * local continuum and including line optical depths */
91  pump12 = HFLines[0].Emis->pump;
92  pump21 = pump12 * HFLines[0].Lo->g / HFLines[0].Hi->g;
93 
94  /* collision rates s-1 within 1s,
95  * were multiplied by collider density when evaluated in CoolEvaluate */
96  /* ContBoltz is Boltzmann factor for wavelength of line */
98  coll21 = HFLines[0].Coll.cs*dense.cdsqte/HFLines[0].Hi->g;
99 
100  /* set up rate (s-1) equations
101  * all process out of 1 that eventually go to 2 */
102  rate12 =
103  /* collision rate (s-1) from 1 to 2 */
104  coll12 +
105  /* direct external continuum pumping (s-1) in 21 cm line - usually dominated by CMB */
106  pump12 +
107  /* pump rate (s-1) up to 3, times fraction that decay to 2, hence net 1-2 */
108  3.*a31*occnu_lya_13 *a32/(a31+a32)+
109  /* pump rate (s-1) up to 4, times fraction that decay to 2, hence net 1-2 */
110  /* >>chng 05 apr 04, GS, degeneracy corrected from 6 to 3 */
111  3.*a41*occnu_lya_14 *a42/(a41+a42);
112 
113  /* set up rate (s-1) equations
114  * all process out of 2 that eventually go to 1 */
115  /* spontaneous + induced 2 -> 1 by external continuum inside 21 cm line */
116  /* >>chng 04 dec 03, do not include spontaneous decay, for numerical stability */
117  rate21 =
118  /* collisional deexcitation */
119  coll21 +
120  /* net spontaneous decay plus external continuum pumping in 21 cm line */
121  pump21 +
122  /* rate from 2 to 3 time fraction that go back to 1, hence net 2 - 1 */
123  /* >>chng 05 apr 04,GS, degeneracy corrected from 2 to unity */
124  occnu_lya_23*a32 * a31/(a31+a32)+
125  occnu_lya_24*a42*a41/(a41+a42);
126 
127  /* x = HFLines[0].Hi->Pop/HFLines[0].Lo->Pop */
128  x = rate12 / SDIV(a21 + rate21);
129 
131  /* the Transitions term is the total population of 1s */
132  HFLines[0].Hi->Pop = (x/(1.+x))* PopTot;
133  HFLines[0].Lo->Pop = (1./(1.+x))* PopTot;
134 
135  /* the population with correction for stimulated emission */
136  HFLines[0].Emis->PopOpc = HFLines[0].Lo->Pop*((3*rate21- rate12) + 3*a21)/SDIV(3*(a21+ rate21));
137 
138  /* number of escaping line photons, used elsewhere for outward beam */
139  HFLines[0].Emis->phots = MAX2(0.,a21*HFLines[0].Hi->Pop);
140 
141  /* intensity of line */
142  HFLines[0].Emis->xIntensity = HFLines[0].Emis->phots*HFLines[0].EnergyErg;
143 
144  /* ratio of collisional to total (collisional + pumped) excitation */
145  HFLines[0].Emis->ColOvTot = coll12 / rate12;
146 
147  /* finally save the spin temperature */
148  if( HFLines[0].Hi->Pop > SMALLFLOAT )
149  {
150  hyperfine.Tspin21cm = TexcLine( &HFLines[0] );
151  /* this line must be non-zero - it does strongly mase in limit_compton_hi_t sim -
152  * in that sim pop ratio goes to unity for a float and TexcLine ret zero */
153  if( hyperfine.Tspin21cm == 0. )
155  }
156  else
157  {
159  }
160 
161  return;
162 }
163 
164 /*H21cm_electron computes rate for H 21 cm from upper to lower excitation by electrons - call by CoolEvaluate
165  * >>refer H1 cs Smith, F.J., 1966, Planet. Space Sci 14, 929 */
166 double H21cm_electron( double temp )
167 {
168  double hold;
169  temp = MIN2(1e4 , temp );
170  /* following fit is from */
171  /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */
172 
173  hold = -9.607 + log10( sqrt(temp)) * sexp( pow(log10(temp) , 4.5 ) / 1800. );
174  hold = pow(10.,hold );
175  return( hold );
176 }
177 
178 /* computes rate for H 21 cm from upper to lower excitation by atomic hydrogen
179  * from
180  * >>refer H1 cs Allison, A.C., & Dalgarno A., 1969, ApJ 158, 423 */
181 /* the following is the best current survey of 21 cm excitation */
182 /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */
183 #if 0
184 STATIC double h21_t_ge_20( double temp )
185 {
186  double y;
187  double x1,
188  teorginal = temp;
189  /* data go up to 1,000K must not go above this */
190  temp = MIN2( 1000.,temp );
191  x1 =1.0/sqrt(temp);
192  y =-21.70880995483007-13.76259674006133*x1;
193  y = exp(y);
194 
195  /* >>chng 02 feb 14, extrapolate above 1e3 K as per Liszt 2001 recommendation
196  * page 699 of */
197  /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */
198  if( teorginal > 1e3 )
199  {
200  y *= pow(teorginal/1e3 , 0.33 );
201  }
202 
203  return( y );
204 }
205 
206 /* this branch for T < 20K, data go down to 1 K */
207 STATIC double h21_t_lt_20( double temp )
208 {
209  double y;
210  double x1;
211 
212  /* must not go below 1K */
213  temp = MAX2( 1., temp );
214  x1 =temp*log(temp);
215  y =9.720710314268267E-08+6.325515312006680E-08*x1;
216  return(y*y);
217 }
218 #endif
219 
220 /* >> chng 04 dec 15, GS. The fitted rate co-efficients (cm3s-1) in the temperature range 1K to 300K is from
221  * >>refer H1 cs Zygelman, B. 2005, ApJ preprint doi:10.1086/'427682".
222  * The rate is 4/3 times the Dalgarno (1969) rate for the
223  temperature range 300K to 1000K. Above 1000K, the rate is extrapolated according to Liszt 2001.*/
224 STATIC double h21_t_ge_10( double temp )
225 {
226  double y;
227  double x1,x2,x3,
228  teorginal = temp;
229  /* data go up to 300K */
230  temp = MIN2( 300., temp );
231  x1 =temp;
232  y =1.4341127e-9+9.4161077e-15*x1-9.2998995e-9/(log(x1))+6.9539411e-9/sqrt(x1)+1.7742293e-8*(log(x1))/pow(x1,2);
233  if( teorginal > 300. )
234  {
235  /* data go up to 1000*/
236  x3 = MIN2( 1000., teorginal );
237  x2 =1.0/sqrt(x3);
238  y =-21.70880995483007-13.76259674006133*x2;
239  y = 1.236686*exp(y);
240 
241  }
242  if( teorginal > 1e3 )
243  {
244  /*data go above 1000*/
245  y *= pow(teorginal/1e3 , 0.33 );
246  }
247  return( y );
248 }
249 /* this branch for T < 10K, data go down to 1 K */
250 STATIC double h21_t_lt_10( double temp )
251 {
252  double y;
253  double x1;
254 
255  /* must not go below 1K */
256  temp = MAX2(1., temp );
257  x1 =temp;
258  y =8.5622857e-10+2.331358e-11*x1+9.5640586e-11*pow((log(x1)),2)-4.6220869e-10*sqrt(x1)-4.1719545e-10/sqrt(x1);
259  return(y);
260 }
261 
262 /*H21cm_H_atom - evaluate H atom spin changing H atom collision rate,
263  * called by CoolEvaluate
264  * >>refer H1 cs Allison, A.C. & Dalgarno, A., 1969, ApJ 158, 423
265  */
266 double H21cm_H_atom( double temp )
267 {
268  double hold;
269  if( temp >= 10. )
270  {
271  hold = h21_t_ge_10( temp );
272  }
273  else
274  {
275  hold = h21_t_lt_10( temp );
276  }
277 
278  return hold;
279 }
280 
281 /*H21cm_proton - evaluate proton spin changing H atom collision rate,
282 * called by CoolEvaluate */
283 double H21cm_proton( double temp )
284 {
285  /*>>refer 21cm p coll Furlanetto, S.R. & Furlanetto, M.R. 2007,
286  *>>refcon MNRAS, doi:10.1111/j.1365-2966.2007.11921.x
287  * previously had used proton rate, which is 3.2 times H0 rate according to
288  *>>refer 21cm cs Liszt, H. A&A, 371, 698 */
289  /* fit to table 1 of first paper */
290  /*--------------------------------------------------------------*
291  TableCurve Function: c:\storage\litera~1\21cm\atomic~1\p21cm.c Jun 20, 2007 3:37:50 PM
292  proton coll deex
293  X= temperature (K)
294  Y= rate coefficient (1e-9 cm3 s-1)
295  Eqn# 4419 y=a+bx+cx^2+dx^(0.5)+elnx/x
296  r2=0.9999445384690351
297  r2adj=0.9999168077035526
298  StdErr=5.559328579039901E-12
299  Fstat=49581.16793656295
300  a= 9.588389834316704E-11
301  b= -5.158891920816405E-14
302  c= 5.895348443553458E-19
303  d= 2.05304960232429E-11
304  e= 9.122617940315725E-10
305  *--------------------------------------------------------------*/
306 
307  double hold;
308  /* only fit this range, did not include T = 1K point which
309  * causes an inflection */
310  temp = MAX2( 2. , temp );
311  temp = MIN2( 2e4 , temp );
312 
313  /* within range of fitted rate coefficients */
314  double x1,x2,x3,x4;
315  x1 = temp;
316  x2 = temp*temp;
317  x3 = sqrt(temp);
318  x4 = log(temp)/temp;
319  hold =9.588389834316704E-11 - 5.158891920816405E-14*x1
320  +5.895348443553458E-19*x2 + 2.053049602324290E-11*x3
321  +9.122617940315725E-10*x4;
322 
323  return hold;
324 }
325 
326 /*
327  * HyperfineCreate, HyperfineCS written July 2001
328  * William Goddard for Gary Ferland
329  * This code calculates line intensities for known
330  * hyperfine transitions.
331  */
332 
333 /* two products, the transition structure HFLines, which contains all information for the lines,
334  * and nHFLines, the number of these lines.
335  *
336  * these are in taulines.h
337  *
338  * info to create them contained in hyperfine.dat
339  *
340  * abundances of nuclei are also in hyperfine.dat, stored in
341  */
342 
343 /* Ion contains twelve varying temperatures, specified above, used for */
344 /* calculating collision strengths. */
345 typedef struct
346 {
347  double strengths[12];
348 } Ion;
349 
350 static Ion *Strengths;
351 
352 /* HyperfineCreat establish space for hf arrays, reads atomic data from hyperfine.dat */
353 void HyperfineCreate(void)
354 {
355  FILE *ioDATA;
356  char chLine[INPUT_LINE_LENGTH];
357  bool lgEOL;
358  /*double c, h, k, N, Ne, q12, q21, upsilon, x;*/
359  realnum spin, wavelength;
360  long int i, j, mass, nelec, ion, nelem;
361 
362  DEBUG_ENTRY( "HyperfineCreate()" );
363 
364  /* list of ion collision strengths for the temperatures listed in table */
365  /* HFLines containing all the data in Hyperfine.dat, and transition is */
366  /* defined in cddefines.h */
367 
368  /*transition *HFLines;*/
369 
370  /* get the line data for the hyperfine lines */
371  if( trace.lgTrace )
372  fprintf( ioQQQ," Hyperfine opening hyperfine.dat:");
373 
374  ioDATA = open_data( "hyperfine.dat", "r" );
375 
376  /* first line is a version number and does not count */
377  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
378  {
379  fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
380  cdEXIT(EXIT_FAILURE);
381  }
382  /* count how many lines are in the file, ignoring all lines
383  * starting with '#' */
384  nHFLines = 0;
385  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
386  {
387  /* we want to count the lines that do not start with #
388  * since these contain data */
389  if( chLine[0] != '#')
390  ++nHFLines;
391  }
392 
393  /* MALLOC the transition HFLines array */
394  HFLines = (transition *)MALLOC( (size_t)(nHFLines)*sizeof(transition) );
395  /* initialize array to impossible values to make sure eventually done right */
396  for( i=0; i< nHFLines; ++i )
397  {
398  TransitionJunk( &HFLines[i] );
399  HFLines[i].Hi = AddState2Stack();
400  HFLines[i].Lo = AddState2Stack();
401  HFLines[i].Emis = AddLine2Stack( true );
402  }
403 
404  Strengths = (Ion *)MALLOC( (size_t)(nHFLines)*sizeof(Ion) );
405  hyperfine.HFLabundance = (realnum *)MALLOC( (size_t)(nHFLines)*sizeof(realnum) );
406 
407  /* now rewind the file so we can read it a second time*/
408  if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
409  {
410  fprintf( ioQQQ, " Hyperfine could not rewind hyperfine.dat.\n");
411  cdEXIT(EXIT_FAILURE);
412  }
413 
414  /* check that magic number is ok, read the line */
415  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
416  {
417  fprintf( ioQQQ, " Hyperfine could not read first line of hyperfine.dat.\n");
418  cdEXIT(EXIT_FAILURE);
419  }
420 
421  /* check that magic number is ok, scan it in */
422  i = 1;
423  nelem = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
424  nelec = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
425  ion = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
426 
427  {
428  /* the following is the set of numbers that appear at the start of hyperfine.dat 01 07 12 */
429  static int iYR=2, iMN=10, iDY=28;
430  if( ( nelem != iYR ) || ( nelec != iMN ) || ( ion != iDY ) )
431  {
432  fprintf( ioQQQ,
433  " Hyperfine: the version of hyperfine.dat in the data directory is not the current version.\n" );
434  fprintf( ioQQQ,
435  " I expected to find the number %i %i %i and got %li %li %li instead.\n" ,
436  iYR, iMN , iDY ,
437  nelem , nelec , ion );
438  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
439  cdEXIT(EXIT_FAILURE);
440  }
441  }
442 
443  /*
444  * scan the string taken from Hyperfine.dat, parsing into
445  * needed variables.
446  * nelem is the atomic number.
447  * IonStg is the ionization stage. Atom = 1, Z+ = 2, Z++ = 3, etc.
448  * Aul is used to find the einstein A in the function GetGF.
449  * most of the variables are floats.
450  */
451 
452  /* this will count the number of lines */
453  j = 0;
454 
455  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
456  {
457  /* only look at lines without '#' in first col */
458  /* make sure line in data file is < 140 char */
459  if( chLine[0] != '#')
460  {
461  double Aul;
462  ASSERT( j < nHFLines );
463 
464  double help[3];
465  sscanf(chLine, "%li%i%le%i%le%le%le%le%le%le%le%le%le%le%le%le%le%le%le",
466  &mass,
467  &HFLines[j].Hi->nelem,
468  &help[0],
469  &HFLines[j].Hi->IonStg,
470  &help[1],
471  &help[2],
472  &Aul,
473  &Strengths[j].strengths[0],
474  &Strengths[j].strengths[1],
475  &Strengths[j].strengths[2],
476  &Strengths[j].strengths[3],
477  &Strengths[j].strengths[4],
478  &Strengths[j].strengths[5],
479  &Strengths[j].strengths[6],
480  &Strengths[j].strengths[7],
481  &Strengths[j].strengths[8],
482  &Strengths[j].strengths[9],
483  &Strengths[j].strengths[10],
484  &Strengths[j].strengths[11]);
485  spin = (realnum)help[0];
486  hyperfine.HFLabundance[j] = (realnum)help[1];
487  wavelength = (realnum)help[2];
488  HFLines[j].Emis->Aul = (realnum)Aul;
489  HFLines[j].Hi->g = (realnum) (2*(spin + .5) + 1);
490  HFLines[j].Lo->g = (realnum) (2*(spin - .5) + 1);
491  HFLines[j].WLAng = wavelength * 1e8f;
492  HFLines[j].EnergyWN = (realnum) (1. / wavelength);
493  HFLines[j].Emis->gf = (realnum)(GetGF(HFLines[j].Emis->Aul, HFLines[j].EnergyWN, HFLines[j].Hi->g));
494  /*fprintf(ioQQQ,"HFLinesss1\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
495  j,HFLines[j].opacity , HFLines[j].Emis->gf , HFLines[j].Emis->Aul , HFLines[j].EnergyWN,HFLines[j].Lo->g);*/
496 
497  ASSERT(HFLines[j].Emis->gf > 0.);
498  /* increment the counter */
499  j++;
500  }
501 
502  }
503 
504  /* closing the file */
505  fclose(ioDATA);
506 
507 # if 0
508  /* for debugging and developing only */
509  /* calculating the luminosity for each isotope */
510  for(i = 0; i < nHFLines; i++)
511  {
512  N = dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1];
513  Ne = dense.eden;
514 
515  h = 6.626076e-27; /* erg * sec */
516  c = 3e10; /* cm / sec */
517  k = 1.380658e-16; /* erg / K */
518 
519  upsilon = HyperfineCS(i);
520  /*statistical weights must still be identified */
521  q21 = COLL_CONST * upsilon / (phycon.sqrte * HFLines[i].Hi->g);
522 
523  q12 = HFLines[i].Hi->g/ HFLines[i].Lo->g * q21 * exp(-1 * h * c * HFLines[i].EnergyWN / (k * phycon.te));
524 
525  x = Ne * q12 / (HFLines[i].Emis->Aul * (1 + Ne * q21 / HFLines[i].Aul));
526  HFLines[i].xIntensity = N * HFLines[i].Emis->Aul * x / (1.0 + x) * h * c / (HFLines[i].WLAng / 1e8);
527 
528  }
529 # endif
530 
531  return;
532 }
533 
534 
535 /*HyperfineCS returns interpolated collision strength for element nelem and ion ion */
536 double HyperfineCS( long i )
537 {
538 
539  /*Search HFLines to find i of ion */
540  int /*i = 0,*/ j = 0;
541 # define N_TE_TABLE 12
542  double slope, upsilon, TemperatureTable[N_TE_TABLE] = {.1e6, .15e6, .25e6, .4e6, .6e6,
543  1.0e6, 1.5e6, 2.5e6, 4e6, 6e6, 10e6, 15e6};
544 
545  DEBUG_ENTRY( "HyperfineCS()" );
546 
547  /*while(i < nHFLines+1 && HFLines[i].Hi->nelem != nelem && HFLines[i].Hi->IonStg != ion)
548  ++i;*/
549  ASSERT( i >= 0. && i <= nHFLines );
550 
551  /*j = temperature */
552  /* calculate actual cooling rate for isotope. */
553  /* The temperature-dependent collision strength must first be calculated. */
554  /* phycon.te is compared to the first, last and intermediate values of strengths[]*/
555  /* and interpolated by taking the log of the collision strength */
556  /* and temperature. */
557 
558  if( phycon.te <= TemperatureTable[0])
559  {
560  /* temperature below bounds of table */
561  upsilon = Strengths[i].strengths[0];
562  }
563  else if( phycon.te >= TemperatureTable[N_TE_TABLE-1])
564  {
565  /* temperature above bounds of table */
566  upsilon = Strengths[i].strengths[N_TE_TABLE-1];
567  }
568  else
569  {
570  j = 1;
571  /* want Table[j-1] < te < Table[j] */
572  /*while ( phycon.te >= TemperatureTable[j] && j < N_TE_TABLE )*/
573  /*while ( TemperatureTable[j] < phycon.te && j < N_TE_TABLE )*/
574  while ( j < N_TE_TABLE && TemperatureTable[j] < phycon.te )
575  j++;
576 
577  ASSERT( j >= 0 && j < N_TE_TABLE);
578  ASSERT( (TemperatureTable[j-1] <= phycon.te ) && (TemperatureTable[j] >= phycon.te) );
579 
580  if( fp_equal( phycon.te, TemperatureTable[j] ) )
581  {
582  upsilon = Strengths[i].strengths[j];
583  }
584  /*phycon.te must be less than TemperatureTable[j], greater than TemperatureTable[j-1][0]*/
585  else if( phycon.te < TemperatureTable[j])
586  {
587  slope = (log10(Strengths[i].strengths[j-1]) - log10(Strengths[i].strengths[j])) /
588  (log10(TemperatureTable[j-1]) - log10(TemperatureTable[j]));
589 
590  upsilon = (log10((phycon.te)) - (log10(TemperatureTable[j-1])-6.))*slope + log10(Strengths[i].strengths[j-1]);
591  upsilon = pow(10., upsilon);
592  }
593  else
594  {
595  upsilon = Strengths[i].strengths[j-1];
596  }
597  }
598 
599  return upsilon;
600 }

Generated for cloudy by doxygen 1.8.3.1