cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
hydroreccool.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*HydroRecCool hydrogen recombination cooling, called by iso_cool */
4 #include "cddefines.h"
5 #include "hydrogenic.h"
6 #include "phycon.h"
7 #include "iso.h"
8 #include "freebound.h"
9 
10 double HydroRecCool(
11  /* n is the prin quantum number on the physical scale */
12  long int n ,
13  /* nelem is the charge on the C scale, 0 is hydrogen */
14  long int nelem)
15 {
16  long int nm1; /* save n - 1*/
17 
18  double fac,
19  hclf_v,
20  x;
21  static const double a[15]={-26.6446988,-26.9066506,-27.0619832,-27.1826903,
22  -27.2783527,-27.3595949,-27.569406,-27.611159,-27.654748,-27.70479,
23  -27.745398,-27.776126,-27.807261,-27.833093,-27.866134};
24  static const double b[15]={-0.40511045,-0.41644707,-0.45834359,-0.49137946,
25  -0.51931762,-0.54971231,-0.18555807,-0.29204736,-0.36741085,
26  -0.45843009,-0.49753713,-0.51418739,-0.52287028,-0.52445456,
27  -0.52292803};
28  static const double c[15]={11.29232731,11.71035693,12.89309608,13.85569087,
29  14.67354775,15.56090323,6.147461,9.0304953,11.094731,13.630431,
30  14.721959,15.172335,15.413946,15.458123,15.428761};
31  static const double d[15]={.067257375,.07638384,.089925637,.102252192,
32  .111016071,.119518918,0.0093832482,0.028119606,0.039357697,0.050378417,
33  0.051644049,0.051367182,0.04938724,0.050139066,0.043085968};
34  static const double e[15]={-1.99108378,-2.26898352,-2.65163846,-3.02333001,
35  -3.29462338,-3.56633674,-1.0019228,-1.5128672,-1.8327058,-2.1866371,
36  -2.2286257,-2.1932699,-2.1205891,-2.1317169,-1.9175186};
37  static const double f[15]={-0.0050802618,-0.005849291,-0.0074854405,-0.0085677543,
38  -0.0093067267,-0.0098455637,0.040903604,0.037491802,0.035618861,
39  0.034132954,0.032418252,0.02947883,0.027393564,0.027607009,0.02433868};
40  static const double g[15]={.166267838,.196780541,.242675042,.282237211,
41  .310204623,.335160025,-0.81087376,-0.73435108,-0.69164333,-0.64907209,
42  -0.61216299,-0.55239109,-0.51048669,-0.51963194,-0.4504203};
43  static const double h[15]={.00020528663,.00027588492,.00033980652,.00041445515,
44  .00046423276,.0005121808,-0.0011986559,-0.0011333973,-0.0010992935,
45  -0.0010878727,-0.0010412678,-0.00095539899,-0.00089141547,-0.00089294364,
46  -0.00079179756};
47  static const double i[15]={-0.0071357493,-0.0099630604,-0.01178647,-0.014696455,
48  -0.01670318,-0.01887373,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.};
49 
50  DEBUG_ENTRY( "HydroRecCool()" );
51 
52  /* confirm that n is 1 or greater. */
53  ASSERT( n > 0 );
54 
55  /* this is log of (temperature divided by charge squared) since sqlogz is log10 Z^2 */
56  x = phycon.telogn[0] - phycon.sqlogz[nelem];
57 
58  /* this will be called for very silly low temperatures, due to high charge.
59  * evaluate at lowest fitted temp, and then extrapolate using known form.
60  * this is ok since these very high charge species do not contribute much
61  * recombination cooling */
62 
63  if( n > 15 || x < 0.2 )
64  {
65  /*double oldfac;*/
66  /* use scale factor from actual recombination rate for this level, this ion,
67  * fac accounts for decreasing efficiency high up
68  * fac is 0.38 at n=15, and 0.32 at 25 */
69  /*oldfac = 0.219 + (2.58 - 2.586/POW2((double)n));
70  oldfac = 0.35*0.69* pow( (15./(double)n ) , 0.67 );
71  fprintf(ioQQQ,"%e %e %e\n", oldfac , fac , fac/oldfac );*/
72  /* >>chng 00 dec 07 use LaMothe & Ferland rates */
73  /* >>refer H recom cool LaMothe, J., & Ferland, G.J., 2001, PASP, in press */
74  fac = HCoolRatio( phycon.te*POW2((double)n) /POW2(nelem+1.) );
75  hclf_v = iso_sp[ipH_LIKE][nelem].fb[n].RadRecomb[ipRecRad]*phycon.te*BOLTZMANN* fac;
76  return( hclf_v );
77  }
78 
79  /* bail if te too high (if this ever happens, change logic so that HCoolRatio is
80  * used in this limit - the process must be small in this case and routine is
81  * well bounded at high-energy end)*/
82  if( x > phycon.TEMP_LIMIT_HIGH_LOG )
83  {
84  fprintf( ioQQQ, " HydroRecCool called with invalid temperature=%e nelem=%li\n",
85  phycon.te , nelem );
87  }
88 
89  /* convert onto c array for n*/
90  nm1 = n - 1;
91 
92  if( nelem == 0 )
93  {
94  /* simple case, most often called, for hydrogen itself */
95  fac = (a[nm1] +
96  c[nm1]*phycon.telogn[0] +
97  e[nm1]*phycon.telogn[1] +
98  g[nm1]*phycon.telogn[2] +
99  i[nm1]*phycon.telogn[3])/
100  (1. + b[nm1]*phycon.telogn[0] +
101  d[nm1]*phycon.telogn[1] +
102  f[nm1]*phycon.telogn[2] +
103  h[nm1]*phycon.telogn[3]);
104  }
105  else
106  {
107  /* hydrogenic ions, expand as powers in t-z2 */
108  fac = (a[nm1] +
109  c[nm1]*x +
110  e[nm1]*POW2( x) +
111  g[nm1]*POW3( x) +
112  i[nm1]*powi( x,4) ) /
113  (1. + b[nm1]*x +
114  d[nm1]*POW2( x ) +
115  f[nm1]*POW3( x ) +
116  h[nm1]*powi( x ,4) );
117  }
118 
119  hclf_v = exp10(fac)*POW3(nelem+1.);
120  return( hclf_v );
121 }
122 
123 /* this function returns the ratio of cooling to recombination as
124  * derived in
125  * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
126 double HCoolRatio(
127  /* the scaled temperature, Tn^2/Z^2 */
128  double t )
129 {
130  double gamma;
131 
132  DEBUG_ENTRY( "HCoolRatio()" );
133 
134  if( t< 1e0 )
135  {
136  gamma = 1.;
137  }
138  else if( t < 7.4e5 )
139  {
140  double y;
141  double x1,x2,x3,x4;
142  x1=t;
143  x2=t*sqrt(t);
144  x3=t*t;
145  x4=t*t*log(t);
146  y=1.000285197084355-7.569939287228937E-06*x1
147  +2.791888685624040E-08*x2-1.289820289839189E-10*x3
148  +7.829204293134294E-12*x4;
149  gamma = y;
150  }
151  else if( t < 5e10 )
152  {
153  double y;
154  double x1,x2,x3,x4,xl;
155  xl = log(t);
156  x1=t;
157  x2=xl*xl;
158  x3=1.0/sqrt(t);
159  x4=xl/(t*t);
160  y=0.2731170438382388+6.086879204730784E-14*x1
161  -0.0003748988159766978*x2+270.2454763661910*x3
162  -1982634355.349780*x4;
163  gamma = y;
164  }
165  else if( t < 3e14 )
166  {
167  double y;
168  double x1,x2;
169  x1=sqrt(t);
170  x2=log(t);
171  y=-17.02819709397900+4.516090033327356E-05*x1
172  +1.088324678258230*x2;
173  gamma = 1/y;
174  }
175  else
176  {
177  /*gamma = 3.85e11 * pow(t , -1. );*/
178  gamma = 1.289e11 * pow(t , -0.9705 );
179  }
180 
181  gamma = MIN2( 1.0, gamma );
182  gamma = MAX2( 0.0, gamma );
183 
184  return( gamma );
185 }
double HydroRecCool(long int n, long int ipZ)
static double x2[63]
double exp10(double x)
Definition: cddefines.h:1368
static double x1[83]
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
double sqlogz[LIMELM]
Definition: phycon.h:86
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
#define POW2
Definition: cddefines.h:979
#define EXIT_FAILURE
Definition: cddefines.h:168
const double TEMP_LIMIT_HIGH_LOG
Definition: phycon.h:123
#define cdEXIT(FAIL)
Definition: cddefines.h:482
double powi(double, long int)
Definition: service.cpp:594
const int ipRecRad
Definition: cddefines.h:333
#define ASSERT(exp)
Definition: cddefines.h:613
const int ipH_LIKE
Definition: iso.h:64
double HCoolRatio(double t)
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
double telogn[7]
Definition: phycon.h:86
#define POW3
Definition: cddefines.h:986
double te
Definition: phycon.h:21