cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*HydroRecCool hydrogen recombination cooling, called by iso_cool */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "hydrogenic.h" 00007 #include "phycon.h" 00008 #include "iso.h" 00009 00010 double HydroRecCool( 00011 /* n is the prin quantum number on the physical scale */ 00012 long int n , 00013 /* nelem is the charge on the C scale, 0 is hydrogen */ 00014 long int nelem) 00015 { 00016 long int nm1; /* save n - 1*/ 00017 00018 double fac, 00019 hclf_v, 00020 x; 00021 static double a[15]={-26.6446988,-26.9066506,-27.0619832,-27.1826903, 00022 -27.2783527,-27.3595949,-27.569406,-27.611159,-27.654748,-27.70479, 00023 -27.745398,-27.776126,-27.807261,-27.833093,-27.866134}; 00024 static double b[15]={-0.40511045,-0.41644707,-0.45834359,-0.49137946, 00025 -0.51931762,-0.54971231,-0.18555807,-0.29204736,-0.36741085, 00026 -0.45843009,-0.49753713,-0.51418739,-0.52287028,-0.52445456, 00027 -0.52292803}; 00028 static double c[15]={11.29232731,11.71035693,12.89309608,13.85569087, 00029 14.67354775,15.56090323,6.147461,9.0304953,11.094731,13.630431, 00030 14.721959,15.172335,15.413946,15.458123,15.428761}; 00031 static double d[15]={.067257375,.07638384,.089925637,.102252192, 00032 .111016071,.119518918,0.0093832482,0.028119606,0.039357697,0.050378417, 00033 0.051644049,0.051367182,0.04938724,0.050139066,0.043085968}; 00034 static double e[15]={-1.99108378,-2.26898352,-2.65163846,-3.02333001, 00035 -3.29462338,-3.56633674,-1.0019228,-1.5128672,-1.8327058,-2.1866371, 00036 -2.2286257,-2.1932699,-2.1205891,-2.1317169,-1.9175186}; 00037 static double f[15]={-0.0050802618,-0.005849291,-0.0074854405,-0.0085677543, 00038 -0.0093067267,-0.0098455637,0.040903604,0.037491802,0.035618861, 00039 0.034132954,0.032418252,0.02947883,0.027393564,0.027607009,0.02433868}; 00040 static double g[15]={.166267838,.196780541,.242675042,.282237211, 00041 .310204623,.335160025,-0.81087376,-0.73435108,-0.69164333,-0.64907209, 00042 -0.61216299,-0.55239109,-0.51048669,-0.51963194,-0.4504203}; 00043 static double h[15]={.00020528663,.00027588492,.00033980652,.00041445515, 00044 .00046423276,.0005121808,-0.0011986559,-0.0011333973,-0.0010992935, 00045 -0.0010878727,-0.0010412678,-0.00095539899,-0.00089141547,-0.00089294364, 00046 -0.00079179756}; 00047 static double i[15]={-0.0071357493,-0.0099630604,-0.01178647,-0.014696455, 00048 -0.01670318,-0.01887373,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.}; 00049 00050 DEBUG_ENTRY( "HydroRecCool()" ); 00051 00052 /* confirm that n is 1 or greater. */ 00053 ASSERT( n > 0 ); 00054 00055 /* this is log of (temperature divided by charge squared) since sqlogz is log10 Z^2 */ 00056 x = phycon.telogn[0] - phycon.sqlogz[nelem]; 00057 00058 /* this will be called for very silly low temperatures, due to high charge. 00059 * evaluate at lowest fitted temp, and then extrapolate using known form. 00060 * this is ok since these very high charge species do not contribute much 00061 * recombination cooling */ 00062 00063 if( n > 15 || x < 0.2 ) 00064 { 00065 /*double oldfac;*/ 00066 /* use scale factor from actual recombination rate for this level, this ion, 00067 * fac accounts for decreasing efficiency high up 00068 * fac is 0.38 at n=15, and 0.32 at 25 */ 00069 /*oldfac = 0.219 + (2.58 - 2.586/POW2((double)n)); 00070 oldfac = 0.35*0.69* pow( (15./(double)n ) , 0.67 ); 00071 fprintf(ioQQQ,"%e %e %e\n", oldfac , fac , fac/oldfac );*/ 00072 /* >>chng 00 dec 07 use LaMothe & Ferland rates */ 00073 /* >>refer H recom cool LaMothe, J., & Ferland, G.J., 2001, PASP, in press */ 00074 fac = HCoolRatio( phycon.te*POW2((double)n) /POW2(nelem+1.) ); 00075 hclf_v = iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad]*phycon.te*BOLTZMANN* fac; 00076 return( hclf_v ); 00077 } 00078 00079 /* bail if te too high (if this ever happens, change logic so that HCoolRatio is 00080 * used in this limit - the process must be small in this case and routine is 00081 * well bounded at high-energy end)*/ 00082 if( x > phycon.TEMP_LIMIT_HIGH_LOG ) 00083 { 00084 fprintf( ioQQQ, " HydroRecCool called with invalid temperature=%e nelem=%li\n", 00085 phycon.te , nelem ); 00086 cdEXIT(EXIT_FAILURE); 00087 } 00088 00089 /* convert onto c array for n*/ 00090 nm1 = n - 1; 00091 00092 if( nelem == 0 ) 00093 { 00094 /* simple case, most often called, for hydrogen itself */ 00095 fac = (a[nm1] + 00096 c[nm1]*phycon.telogn[0] + 00097 e[nm1]*phycon.telogn[1] + 00098 g[nm1]*phycon.telogn[2] + 00099 i[nm1]*phycon.telogn[3])/ 00100 (1. + b[nm1]*phycon.telogn[0] + 00101 d[nm1]*phycon.telogn[1] + 00102 f[nm1]*phycon.telogn[2] + 00103 h[nm1]*phycon.telogn[3]); 00104 } 00105 else 00106 { 00107 /* hydrogenic ions, expand as powers in t-z2 */ 00108 fac = (a[nm1] + 00109 c[nm1]*x + 00110 e[nm1]*POW2( x) + 00111 g[nm1]*POW3( x) + 00112 i[nm1]*powi( x,4) ) / 00113 (1. + b[nm1]*x + 00114 d[nm1]*POW2( x ) + 00115 f[nm1]*POW3( x ) + 00116 h[nm1]*powi( x ,4) ); 00117 } 00118 00119 hclf_v = pow(10.,fac)*POW3(nelem+1.); 00120 return( hclf_v ); 00121 } 00122 00123 /* this function returns the ratio of cooling to recombination as 00124 * derived in 00125 * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */ 00126 double HCoolRatio( 00127 /* the scaled temperature, Tn^2/Z^2 */ 00128 double t ) 00129 { 00130 double gamma; 00131 00132 DEBUG_ENTRY( "HCoolRatio()" ); 00133 00134 if( t< 1e0 ) 00135 { 00136 gamma = 1.; 00137 } 00138 else if( t < 7.4e5 ) 00139 { 00140 double y; 00141 double x1,x2,x3,x4; 00142 x1=t; 00143 x2=t*sqrt(t); 00144 x3=t*t; 00145 x4=t*t*log(t); 00146 y=1.000285197084355-7.569939287228937E-06*x1 00147 +2.791888685624040E-08*x2-1.289820289839189E-10*x3 00148 +7.829204293134294E-12*x4; 00149 gamma = y; 00150 } 00151 else if( t < 5e10 ) 00152 { 00153 double y; 00154 double x1,x2,x3,x4,xl; 00155 xl = log(t); 00156 x1=t; 00157 x2=xl*xl; 00158 x3=1.0/sqrt(t); 00159 x4=xl/(t*t); 00160 y=0.2731170438382388+6.086879204730784E-14*x1 00161 -0.0003748988159766978*x2+270.2454763661910*x3 00162 -1982634355.349780*x4; 00163 gamma = y; 00164 } 00165 else if( t < 3e14 ) 00166 { 00167 double y; 00168 double x1,x2; 00169 x1=sqrt(t); 00170 x2=log(t); 00171 y=-17.02819709397900+4.516090033327356E-05*x1 00172 +1.088324678258230*x2; 00173 gamma = 1/y; 00174 } 00175 else 00176 { 00177 /*gamma = 3.85e11 * pow(t , -1. );*/ 00178 gamma = 1.289e11 * pow(t , -0.9705 ); 00179 } 00180 00181 return( gamma ); 00182 }