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 /*OpacityCreateAll compute initial set of opacities for all species */ 00004 /*OpacityCreate1Element generate ionic subshell opacities by calling t_ADfA::Inst().phfit */ 00005 /*opacity_more_memory allocate more memory for opacity stack */ 00006 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */ 00007 /*hmiopc derive total H- H minus opacity */ 00008 /*rayleh compute Rayleigh scattering cross section for Lya */ 00009 /*OpacityValenceRescale routine to rescale non-OP valence shell cross sections */ 00010 /*Yan_H2_CS - cross sections for the photoionization of H2 */ 00011 /****************************************************************************** 00012 *NB NB NB NB NB NB NB NB NB NB NB NB NB NB 00013 * everything set here must be written to the opacity store files 00014 * 00015 ****************************************************************************** */ 00016 #include "cddefines.h" 00017 #include "physconst.h" 00018 #include "dense.h" 00019 #include "continuum.h" 00020 #include "iso.h" 00021 #include "hydrogenic.h" 00022 #include "oxy.h" 00023 #include "trace.h" 00024 #include "heavy.h" 00025 #include "rfield.h" 00026 #include "hmi.h" 00027 #include "atmdat.h" 00028 #include "punch.h" 00029 #include "grains.h" 00030 #include "thirdparty.h" 00031 #include "hydro_bauman.h" 00032 #include "opacity.h" 00033 #include "helike_recom.h" 00034 00035 static const int NCSH2P = 10; 00036 00037 /* limit to number of opacity cells available in the opacity stack*/ 00038 static long int ndimOpacityStack = 2100000L; 00039 00040 /*Yan_H2_CS - cross sections for the photoionization of H2 00041 * may be represented analytically in the paper 00042 *>>refer H2 photo cs Yan, M.; Sadeghpour, H. R.; Dalgarno, A. 1998, ApJ, 496, 1044 00043 * This is revised version given in ERRATUM 2001, ApJ, 559, 1194 00044 * return value is photo cs in cm-2 */ 00045 STATIC double Yan_H2_CS( double energy_ryd /* photon energy in ryd */) 00046 { 00047 00048 double energy_keV; /* keV */ 00049 double cross_section; /* barns */ 00050 double x; /* x = E/15.4 */ 00051 double xsqrt , x15 , x2; 00052 double energy = energy_ryd * EVRYD; 00053 00054 DEBUG_ENTRY( "Yan_H2_CS()" ); 00055 00056 /* energy relative to threshold */ 00057 x = energy / 15.4; 00058 energy_keV = energy/1000.0; 00059 xsqrt = sqrt(x); 00060 x15 = x * xsqrt; 00061 x2 = x*x; 00062 00063 if( energy < 15.4 ) 00064 { 00065 cross_section = 0.; 00066 } 00067 00068 else if(energy >= 15.4 && energy < 18.) 00069 { 00070 cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2); 00071 /* this equation is obviously negative for x = 1 */ 00072 cross_section = MAX2( 0. , cross_section ); 00073 } 00074 00075 else if(energy >= 18. && energy <= 30.) 00076 { 00077 cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5); 00078 } 00079 00080 else if(energy > 30. && energy <= 85.) 00081 { 00082 cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5); 00083 } 00084 00085 /* the high-energy tail */ 00086 else 00087 { 00088 cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*x2) - 81.885/(x*x2))/pow(energy_keV,3.5); 00089 } 00090 00091 return( cross_section * 1e-24 ); 00092 } 00093 00094 /*OpacityCreate1Element generate opacities for entire element by calling t_ADfA::Inst().phfit */ 00095 STATIC void OpacityCreate1Element(long int nelem); 00096 00097 /*opacity_more_memory allocate more memory for opacity stack */ 00098 STATIC void opacity_more_memory(void); 00099 00100 /*hmiopc derive total H- H minus opacity */ 00101 STATIC double hmiopc(double freq); 00102 00103 /*rayleh compute Rayleigh scattering cross section for Lya */ 00104 STATIC double rayleh(double ener); 00105 00106 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */ 00107 STATIC double Opacity_iso_photo_cs( realnum energy , long ipISO , long nelem , long n ); 00108 00109 /*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */ 00110 STATIC void OpacityCreateReilMan(long int low, 00111 long int ihi, 00112 const realnum cross[], 00113 long int ncross, 00114 long int *ipop, 00115 const char *chLabl); 00116 00117 static bool lgRealloc = false; 00118 00119 /*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */ 00120 STATIC void OpacityCreatePowerLaw( 00121 /* lower energy limit on continuum mesh */ 00122 long int ilo, 00123 /* upper energy limit on continuum mesh */ 00124 long int ihi, 00125 /* threshold cross section */ 00126 double cross, 00127 /* power law index */ 00128 double s, 00129 /* pointer to opacity offset where this starts */ 00130 long int *ip); 00131 00132 /*ofit compute cross sections for all shells of atomic oxygen */ 00133 STATIC double ofit(double e, 00134 realnum opart[]); 00135 00136 /*OpacityValenceRescale routine to rescale non-OP valence shell cross sections for atom */ 00137 STATIC void OpacityValenceRescale( 00138 /* element number on C scale */ 00139 long int nelem , 00140 /* scale factor, must be >= 0. */ 00141 double scale ) 00142 { 00143 00144 long int ion , nshell , low , ihi , ipop , ip; 00145 00146 DEBUG_ENTRY( "OpacityValenceRescale()" ); 00147 00148 /* return if element is not turned on 00149 * >>chng 05 oct 19, this had not been done, so low in the opacity offset below was 00150 * not set, and opacity index was negative - only problem when K turned off */ 00151 if( !dense.lgElmtOn[nelem] ) 00152 { 00153 return; 00154 } 00155 00156 /* scale better be >= 0. */ 00157 ASSERT( scale >= 0. ); 00158 00159 ion = 0; 00160 /* this is valence shell on C scale */ 00161 nshell = Heavy.nsShells[nelem][ion] - 1; 00162 00163 /* set lower and upper limits to this range */ 00164 low = opac.ipElement[nelem][ion][nshell][0]; 00165 ihi = opac.ipElement[nelem][ion][nshell][1]; 00166 ipop = opac.ipElement[nelem][ion][nshell][2]; 00167 00168 /* loop over energy range of this shell */ 00169 for( ip=low-1; ip < ihi; ip++ ) 00170 { 00171 opac.OpacStack[ip-low+ipop] *= scale; 00172 } 00173 return; 00174 } 00175 00176 void OpacityCreateAll(void) 00177 { 00178 long int i, 00179 ipISO , 00180 n , 00181 need , 00182 nelem; 00183 00184 realnum opart[7]; 00185 00186 double crs, 00187 dx, 00188 eps, 00189 thom, 00190 thres, 00191 x; 00192 00193 /* >>chng 02 may 29, change to lgOpacMalloced */ 00194 /* remember whether opacities have ever been evaluated in this coreload 00195 static bool lgOpEval = false; */ 00196 00197 /* fits to cross section for photo dist of H_2^+ */ 00198 static const realnum csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f, 00199 6.0f,14.28f,2.7f}; 00200 00201 DEBUG_ENTRY( "OpacityCreateAll()" ); 00202 00203 /* H2+ h2plus h2+ */ 00204 00205 /* make and print dust opacities 00206 * fill in dstab and dstsc, totals, zero if no dust 00207 * may be different if different grains are turned on */ 00208 GrainsInit(); 00209 00210 /* flag lgOpacMalloced says whether opacity stack has been generated 00211 * only do this one time per core load */ 00212 /* >>chng 02 may 29, from lgOpEval to lgOpacMalloced */ 00213 if( lgOpacMalloced ) 00214 { 00215 /* this is not the first time code called */ 00216 if( trace.lgTrace ) 00217 { 00218 fprintf( ioQQQ, " OpacityCreateAll called but NOT evaluated since already done.\n" ); 00219 } 00220 return; 00221 } 00222 00223 lgOpacMalloced = true; 00224 00225 /* create the space for the opacity stack */ 00226 opac.OpacStack = (double*)MALLOC((size_t)ndimOpacityStack*sizeof(double)); 00227 00228 if( trace.lgTrace ) 00229 { 00230 fprintf( ioQQQ, " OpacityCreateAll called, evaluating.\n" ); 00231 } 00232 00233 /* zero out opac since this array sometimes addressed before OpacityAddTotal called */ 00234 for( i=0; i < rfield.nupper; i++ ) 00235 { 00236 opac.opacity_abs[i] = 0.; 00237 } 00238 00239 /* nOpacTot is number of opacity cells in OpacStack filled so far by opacity generating routines */ 00240 opac.nOpacTot = 0; 00241 00242 /* photoionization of h, he-like iso electronic sequences */ 00243 for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO ) 00244 { 00245 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00246 { 00247 if( dense.lgElmtOn[nelem] ) 00248 { 00249 long int nupper; 00250 00251 /* this is the opacity offset in the general purpose pointer array */ 00252 /* indices are type, shell. ion, element, so this is the inner shell, 00253 * NB - this works for H and He, but the second index should be 1 for Li */ 00254 opac.ipElement[nelem][nelem-ipISO][0][2] = opac.nOpacTot + 1; 00255 00256 /* gound state goes to high-energy limit of code, 00257 * but excited states only go up to edge for ground */ 00258 nupper = rfield.nupper; 00259 for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ ) 00260 { 00261 /* this is array index to the opacity offset */ 00262 iso.ipOpac[ipISO][nelem][n] = opac.nOpacTot + 1; 00263 00264 /* first make sure that first energy point is at least near the limit */ 00265 /* >>chng 01 sep 23, increased factor form 0.98 to 0.94, needed since cells now go 00266 * so far into radio, where resolution is poor */ 00267 ASSERT( rfield.AnuOrg[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0.94f * 00268 iso.xIsoLevNIonRyd[ipISO][nelem][n] ); 00269 00270 /* number of cells we will need to do this level */ 00271 need = nupper - iso.ipIsoLevNIonCon[ipISO][nelem][n] + 1; 00272 ASSERT( need > 0 ); 00273 00274 if( opac.nOpacTot + need > ndimOpacityStack ) 00275 opacity_more_memory(); 00276 00277 for( i=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; i < nupper; i++ ) 00278 { 00279 opac.OpacStack[i-iso.ipIsoLevNIonCon[ipISO][nelem][n]+iso.ipOpac[ipISO][nelem][n]] = 00280 Opacity_iso_photo_cs( rfield.AnuOrg[i] , ipISO , nelem , n ); 00281 } 00282 00283 opac.nOpacTot += need; 00284 /* for all excited levels high-energy limit is edge for ground state */ 00285 nupper = iso.ipIsoLevNIonCon[ipISO][nelem][0]; 00286 } 00287 } 00288 } 00289 } 00290 00291 /* This check will get us through Klein-Nishina below. */ 00292 /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */ 00293 if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + rfield.nupper > ndimOpacityStack ) 00294 opacity_more_memory(); 00295 00296 /* Lyman alpha damping wings - Rayleigh scattering */ 00297 opac.ipRayScat = opac.nOpacTot + 1; 00298 for( i=0; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ ) 00299 { 00300 opac.OpacStack[i-1+opac.ipRayScat] = rayleh(rfield.AnuOrg[i]); 00301 } 00302 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; 00303 00304 /* ============================================================== 00305 * this block of code defines the electron scattering cross section 00306 * for all energies */ 00307 00308 /* assume Thomson scattering up to ipCKshell, 20.6 Ryd=0.3 keV */ 00309 thom = 6.65e-25; 00310 opac.iopcom = opac.nOpacTot + 1; 00311 for( i=0; i < opac.ipCKshell; i++ ) 00312 { 00313 opac.OpacStack[i-1+opac.iopcom] = thom; 00314 /*fprintf(ioQQQ,"%.3e\t%.3e\n", 00315 rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/ 00316 } 00317 00318 /* Klein-Nishina from eqn 7.5, 00319 * >>refer Klein-Nishina cs Rybicki and Lightman */ 00320 for( i=opac.ipCKshell; i < rfield.nupper; i++ ) 00321 { 00322 dx = rfield.AnuOrg[i]/3.7573e4; 00323 00324 opac.OpacStack[i-1+opac.iopcom] = thom*3.e0/4.e0*((1.e0 + 00325 dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+ 00326 2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0* 00327 dx)/POW3(1.e0 + 2.e0*dx)); 00328 /*fprintf(ioQQQ,"%.3e\t%.3e\n", 00329 rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/ 00330 } 00331 opac.nOpacTot += rfield.nupper - 1 + 1; 00332 00333 /* ============================================================== */ 00334 00335 /* This check will get us through "H- hminus H minus bound-free opacity" below. */ 00336 /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */ 00337 if( opac.nOpacTot + 3*rfield.nupper - opac.ippr + iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 2 > ndimOpacityStack ) 00338 opacity_more_memory(); 00339 00340 /* pair production */ 00341 opac.ioppr = opac.nOpacTot + 1; 00342 for( i=opac.ippr-1; i < rfield.nupper; i++ ) 00343 { 00344 /* pair production heating rate for unscreened H + He 00345 * fit to figure 41 of Experimental Nuclear Physics, 00346 * Vol1, E.Segre, ed */ 00347 00348 x = rfield.AnuOrg[i]/7.512e4*2.; 00349 00350 opac.OpacStack[i-opac.ippr+opac.ioppr] = 5.793e-28* 00351 POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. + 00352 x*(0.130471301 + x*0.000524906))); 00353 /*fprintf(ioQQQ,"%.3e\t%.3e\n", 00354 rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-opac.ippr+opac.ioppr] );*/ 00355 } 00356 opac.nOpacTot += rfield.nupper - opac.ippr + 1; 00357 00358 /* brems (free-free) opacity */ 00359 opac.ipBrems = opac.nOpacTot + 1; 00360 00361 for( i=0; i < rfield.nupper; i++ ) 00362 { 00363 /* missing factor of 1E-20 to avoid underflow 00364 * free free opacity needs g(ff)*(1-exp(hn/kT))/SQRT(T)*1E-20 */ 00365 opac.OpacStack[i-1+opac.ipBrems] = 00366 /*(realnum)(1.03680e-18/POW3(rfield.AnuOrg[i]));*/ 00367 /* >>chng 00 jun 05, descale by 1e10 so that underflow at high-energy 00368 * end does not occur */ 00369 1.03680e-8/POW3(rfield.AnuOrg[i]); 00370 } 00371 opac.nOpacTot += rfield.nupper - 1 + 1; 00372 00373 opac.iphmra = opac.nOpacTot + 1; 00374 for( i=0; i < rfield.nupper; i++ ) 00375 { 00376 /* following is ratio of h minus to neut h bremss opacity */ 00377 opac.OpacStack[i-1+opac.iphmra] = 0.1175*rfield.anusqr[i]; 00378 } 00379 opac.nOpacTot += rfield.nupper - 1 + 1; 00380 00381 opac.iphmop = opac.nOpacTot + 1; 00382 for( i=hmi.iphmin-1; i < iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]; i++ ) 00383 { 00384 /* H- hminus H minus bound-free opacity */ 00385 opac.OpacStack[i-hmi.iphmin+opac.iphmop] = 00386 hmiopc(rfield.AnuOrg[i]); 00387 } 00388 opac.nOpacTot += iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 1; 00389 00390 /* ============================================================== */ 00391 00392 /* This check will get us through "H2 photoionization cross section" below. */ 00393 /* >>chng 07 oct 10, by Ryan. Added this check for allotted memory. */ 00394 if( opac.nOpacTot + rfield.nupper - opac.ipH2_photo_thresh + 1 > ndimOpacityStack ) 00395 opacity_more_memory(); 00396 00397 opac.ipH2_photo_opac_offset = opac.nOpacTot + 1; 00398 for( i=opac.ipH2_photo_thresh-1; i < rfield.nupper; i++ ) 00399 { 00400 /* >>chng 05 nov 24, add H2 photoionization cross section */ 00401 opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] = 00402 Yan_H2_CS(rfield.AnuOrg[i]); 00403 /*fprintf(ioQQQ,"DEBUG H2 photo\t%.3e\t%.3e\t%.3e\n", 00404 rfield.AnuOrg[i], 00405 opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] , 00406 Opacity_iso_photo_cs( rfield.AnuOrg[i] , 0 , 0 , 0 )*2. );*/ 00407 } 00408 opac.nOpacTot += rfield.nupper - opac.ipH2_photo_thresh + 1; 00409 00410 /* H2+ H2P h2plus photoabsorption 00411 * cs from 00412 * >>refer H2+ photodissoc Buckingham, R.A., Reid, S., & Spence, R. 1952, MNRAS 112, 382, 0 K temp */ 00413 OpacityCreateReilMan(opac.ih2pnt[0],opac.ih2pnt[1],csh2p,NCSH2P,&opac.ih2pof, 00414 "H2+ "); 00415 00416 /* This check will get us through "HeI singlets neutral helium ground" below. */ 00417 /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */ 00418 if( opac.nOpacTot + rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1 > ndimOpacityStack ) 00419 opacity_more_memory(); 00420 00421 /* HeI singlets neutral helium ground */ 00422 opac.iophe1[0] = opac.nOpacTot + 1; 00423 opac.ipElement[ipHELIUM][0][0][2] = opac.iophe1[0]; 00424 for( i=iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1; i < rfield.nupper; i++ ) 00425 { 00426 crs = t_ADfA::Inst().phfit(2,2,1,rfield.AnuOrg[i]*EVRYD); 00427 opac.OpacStack[i-iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]+opac.iophe1[0]] = 00428 crs*1e-18; 00429 } 00430 opac.nOpacTot += rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1; 00431 00432 /* these are opacity offset points that would be defined in OpacityCreate1Element, 00433 * but this routine will not be called for H and He 00434 * generate all heavy element opacities, everything heavier than He, 00435 * nelem is on the C scale, so Li is 2 */ 00436 /*>>chng 99 jan 27, do not reevaluate hydrogenic opacity below */ 00437 for( nelem=2; nelem < LIMELM; nelem++ ) 00438 { 00439 if( dense.lgElmtOn[nelem] ) 00440 { 00441 OpacityCreate1Element(nelem); 00442 } 00443 } 00444 00445 /* option to rescale atoms of some elements that were not done by opacity project 00446 * the valence shell - two arguments - element number on C scale, and scale factor */ 00447 /*>>chng 05 sep 26, fudge factor to get atomic K fraction along well defined line of sight 00448 * to be observed value - this is ratio of cross sections, actual value is very uncertain since 00449 * differences betweeen Verner & opacity project are huge */ 00450 OpacityValenceRescale( ipPOTASSIUM , 5. ); 00451 00452 /* now add on some special cases, where exicted states, etc, come in */ 00453 /* Nitrogen 00454 * >>refer n1 photo Henry, R., ApJ 161, 1153. 00455 * photoionization of excited level of N+ */ 00456 OpacityCreatePowerLaw(opac.in1[0],opac.in1[1],9e-18,1.75,&opac.in1[2]); 00457 00458 /* atomic Oxygen 00459 * only do this if 1996 Verner results are used */ 00460 if( dense.lgElmtOn[ipOXYGEN] && t_ADfA::Inst().get_version() == PHFIT96 ) 00461 { 00462 /* This check will get us through this loop. */ 00463 /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */ 00464 if( opac.nOpacTot + opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1 > ndimOpacityStack ) 00465 opacity_more_memory(); 00466 00467 /* integrate over energy range of the valence shell of atomic oxygen*/ 00468 for( i=opac.ipElement[ipOXYGEN][0][2][0]-1; i < opac.ipElement[ipOXYGEN][0][2][1]; i++ ) 00469 { 00470 /* call special routine to evaluate partial cross section for OI shells */ 00471 eps = rfield.AnuOrg[i]*EVRYD; 00472 crs = ofit(eps,opart); 00473 00474 /* this will be total cs of all processes leaving shell 3 */ 00475 crs = opart[0]; 00476 for( n=1; n < 6; n++ ) 00477 { 00478 /* add up table of cross sections */ 00479 crs += opart[n]; 00480 } 00481 /* convert to cgs and overwrite cross sections set by OpacityCreate1Element */ 00482 crs *= 1e-18; 00483 opac.OpacStack[i-opac.ipElement[ipOXYGEN][0][2][0]+opac.ipElement[ipOXYGEN][0][2][2]] = crs; 00484 } 00485 /* >>chng 02 may 09 - this was a significant error */ 00486 /* >>chng 02 may 08, by Ryan. This loop did not update total slots filled. */ 00487 opac.nOpacTot += opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1; 00488 } 00489 00490 /* Henry nubmers for 1S excit state of OI, OP data very sparse */ 00491 OpacityCreatePowerLaw(opac.ipo1exc[0],opac.ipo1exc[1],4.64e-18,0.,&opac.ipo1exc[2]); 00492 00493 /* photoionization of excited level of O2+ 1D making 5007 00494 * fit to TopBase Opacity Project cs */ 00495 OpacityCreatePowerLaw(opac.ipo3exc[0],opac.ipo3exc[1],3.8e-18,0.,&opac.ipo3exc[2]); 00496 00497 /* photoionization of excited level of O2+ 1S making 4363 */ 00498 OpacityCreatePowerLaw(opac.ipo3exc3[0],opac.ipo3exc3[1],5.5e-18,0.01, 00499 &opac.ipo3exc3[2]); 00500 00501 /* This check will get us through the next two steps. */ 00502 /* >>chng 02 may 08, by Ryan. Added this and other checks for allotted memory. */ 00503 if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1 00504 + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1 > ndimOpacityStack ) 00505 opacity_more_memory(); 00506 00507 /* photoionization to excited states of O+ */ 00508 opac.iopo2d = opac.nOpacTot + 1; 00509 thres = rfield.AnuOrg[oxy.i2d-1]; 00510 for( i=oxy.i2d-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]; i++ ) 00511 { 00512 crs = 3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38* 00513 pow(rfield.AnuOrg[i]/thres,-2.5)); 00514 00515 opac.OpacStack[i-oxy.i2d+opac.iopo2d] = crs; 00516 } 00517 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1; 00518 00519 /* magnesium 00520 * photoionization of excited level of Mg+ 00521 * fit to opacity project data Dima got */ 00522 opac.ipOpMgEx = opac.nOpacTot + 1; 00523 for( i=opac.ipmgex-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ ) 00524 { 00525 opac.OpacStack[i-opac.ipmgex+opac.ipOpMgEx] = 00526 (0.2602325880970085 + 00527 445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))* 00528 1e-18; 00529 } 00530 opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1; 00531 00532 ASSERT( opac.nOpacTot < ndimOpacityStack ); 00533 00534 /* Calcium 00535 * excited states of Ca+ */ 00536 OpacityCreatePowerLaw(opac.ica2ex[0],opac.ica2ex[1],4e-18,1.,&opac.ica2op); 00537 00538 if( trace.lgTrace ) 00539 { 00540 fprintf( ioQQQ, 00541 " OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n", 00542 opac.nOpacTot, ndimOpacityStack ); 00543 } 00544 00545 /* option to compile opacities into file for later use 00546 * this is executed if the 'compile opacities' command is entered */ 00547 if( opac.lgCompileOpac ) 00548 { 00549 fprintf( ioQQQ, "The COMPILE OPACITIES command is currently not supported\n" ); 00550 cdEXIT(EXIT_FAILURE); 00551 } 00552 00553 if( lgRealloc ) 00554 fprintf(ioQQQ, 00555 " Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" , opac.nOpacTot); 00556 return; 00557 } 00558 /*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */ 00559 STATIC void OpacityCreatePowerLaw( 00560 /* lower energy limit on continuum mesh */ 00561 long int ilo, 00562 /* upper energy limit on continuum mesh */ 00563 long int ihi, 00564 /* threshold cross section */ 00565 double cross, 00566 /* power law index */ 00567 double s, 00568 /* pointer to opacity offset where this starts */ 00569 long int *ip) 00570 { 00571 long int i; 00572 double thres; 00573 00574 DEBUG_ENTRY( "OpacityCreatePowerLaw()" ); 00575 00576 /* non-positive cross section is unphysical */ 00577 ASSERT( cross > 0. ); 00578 00579 /* place in the opacity stack where we will stuff cross sections */ 00580 *ip = opac.nOpacTot + 1; 00581 ASSERT( *ip > 0 ); 00582 ASSERT( ilo > 0 ); 00583 thres = rfield.anu[ilo-1]; 00584 00585 if( opac.nOpacTot + ihi - ilo + 1 > ndimOpacityStack ) 00586 opacity_more_memory(); 00587 00588 for( i=ilo-1; i < ihi; i++ ) 00589 { 00590 opac.OpacStack[i-ilo+*ip] = cross*pow(rfield.anu[i]/thres,-s); 00591 } 00592 00593 opac.nOpacTot += ihi - ilo + 1; 00594 return; 00595 } 00596 00597 /*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */ 00598 STATIC void OpacityCreateReilMan(long int low, 00599 long int ihi, 00600 const realnum cross[], 00601 long int ncross, 00602 long int *ipop, 00603 const char *chLabl) 00604 { 00605 long int i, 00606 ics, 00607 j, 00608 ncr; 00609 00610 const int NOP = 100; 00611 realnum cs[NOP], 00612 en[NOP], 00613 slope; 00614 00615 DEBUG_ENTRY( "OpacityCreateReilMan()" ); 00616 00617 /* this is the opacity entering routine designed for 00618 * the Reilman and Manson tables. It works with incident 00619 * photon energy (entered in eV) and cross sections in megabarns 00620 * */ 00621 *ipop = opac.nOpacTot + 1; 00622 ASSERT( *ipop > 0 ); 00623 00624 ncr = ncross/2; 00625 if( ncr > NOP ) 00626 { 00627 fprintf( ioQQQ, " Too many opacities were entered into OpacityCreateReilMan. Increase the value of NOP.\n" ); 00628 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl ); 00629 cdEXIT(EXIT_FAILURE); 00630 } 00631 00632 /* the array CROSS has ordered pairs of elements. 00633 * the first is the energy in eV (not Ryd) 00634 * and the second is the cross section in megabarns */ 00635 for( i=0; i < ncr; i++ ) 00636 { 00637 en[i] = cross[i*2]/13.6f; 00638 cs[i] = cross[(i+1)*2-1]*1e-18f; 00639 } 00640 00641 ASSERT( low>0 ); 00642 if( en[0] > rfield.anu[low-1] ) 00643 { 00644 fprintf( ioQQQ, 00645 " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" ); 00646 fprintf( ioQQQ, 00647 " The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n", 00648 rfield.anu[low-1]*EVRYD, en[0]*EVRYD ); 00649 00650 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl ); 00651 fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" ); 00652 fprintf( ioQQQ, " " ); 00653 00654 for( i=0; i < ncross; i++ ) 00655 { 00656 fprintf( ioQQQ, "%11.4e", cross[i] ); 00657 } 00658 00659 fprintf( ioQQQ, "\n" ); 00660 cdEXIT(EXIT_FAILURE); 00661 } 00662 00663 slope = (cs[1] - cs[0])/(en[1] - en[0]); 00664 ics = 1; 00665 00666 if( opac.nOpacTot + ihi - low + 1 > ndimOpacityStack ) 00667 opacity_more_memory(); 00668 00669 /* now fill in the opacities using linear interpolation */ 00670 for( i=low-1; i < ihi; i++ ) 00671 { 00672 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] ) 00673 { 00674 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 00675 en[ics-1]); 00676 } 00677 00678 else 00679 { 00680 ics += 1; 00681 if( ics + 1 > ncr ) 00682 { 00683 fprintf( ioQQQ, " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" ); 00684 fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n", 00685 rfield.anu[i]*13.6, en[ncr-1]*13.6 ); 00686 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl 00687 ); 00688 fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n", 00689 en[0]*13.65 ); 00690 fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n", 00691 rfield.anu[ihi-1]*13.6 ); 00692 fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n", 00693 rfield.anu[low-1]*13.6 ); 00694 cdEXIT(EXIT_FAILURE); 00695 } 00696 00697 slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]); 00698 if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] ) 00699 { 00700 opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 00701 en[ics-1]); 00702 } 00703 else 00704 { 00705 ASSERT( i > 0); 00706 fprintf( ioQQQ, " Internal logical error in OpacityCreateReilMan.\n" ); 00707 fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n", 00708 rfield.anu[i]*13.6, i, en[ics-1], en[ics] ); 00709 00710 fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n", 00711 rfield.anu[i-1]*13.6 ); 00712 00713 fprintf( ioQQQ, " Here comes the energy array. ICS=%4ld\n", 00714 ics ); 00715 00716 for( j=0; j < ncr; j++ ) 00717 { 00718 fprintf( ioQQQ, "%10.2e", en[j] ); 00719 } 00720 fprintf( ioQQQ, "\n" ); 00721 00722 fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl ); 00723 cdEXIT(EXIT_FAILURE); 00724 } 00725 } 00726 } 00727 /* >>chng 02 may 09, this was a significant logcal error */ 00728 /* >>chng 02 may 08, by Ryan. This routine did not update the total slots filled. */ 00729 opac.nOpacTot += ihi - low + 1; 00730 return; 00731 } 00732 00733 00734 /*ofit compute cross sections for all shells of atomic oxygen */ 00735 STATIC double ofit(double e, 00736 realnum opart[]) 00737 { 00738 double otot, 00739 q, 00740 x; 00741 00742 static const double y[7][5] = { 00743 {8.915,3995.,3.242,10.44,0.0}, 00744 {11.31,1498.,5.27,7.319,0.0}, 00745 {10.5,1.059e05,1.263,13.04,0.0}, 00746 {19.49,48.47,8.806,5.983,0.0}, 00747 {50.,4.244e04,0.1913,7.012,4.454e-02}, 00748 {110.5,0.1588,148.3,-3.38,3.589e-02}, 00749 {177.4,32.37,381.2,1.083,0.0} 00750 }; 00751 static const double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.}; 00752 static const long l[7]={1,1,1,0,1,1,0}; 00753 00754 DEBUG_ENTRY( "ofit()" ); 00755 /*compute cross sections for all shells of atomic oxygen 00756 * Photoionization of OI 00757 * Input parameter: e - photon energy, eV 00758 * Output parameters: otot - total photoionization cross section, Mb 00759 * opart(1) - 2p-shell photoionization, goes to 4So 00760 * opart(2) - 2p-shell photoionization, goes to 2Do 00761 * opart(3) - 2p-shell photoionization, goes to 2Po 00762 * opart(4) - 2s-shell photoionization 00763 * opart(5) - double photoionization, goes to O++ 00764 * opart(6) - triple photoionization, goes to O+++ 00765 * opart(7) - 1s-shell photoionization */ 00766 00767 otot = 0.0; 00768 00769 for( int i=0; i < 7; i++ ) 00770 { 00771 opart[i] = 0.0; 00772 } 00773 00774 for( int i=0; i < 7; i++ ) 00775 { 00776 if( e >= eth[i] ) 00777 { 00778 q = 5.5 - 0.5*y[i][3] + l[i]; 00779 00780 x = e/y[i][0]; 00781 00782 opart[i] = (realnum)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/ 00783 pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3])); 00784 00785 otot += opart[i]; 00786 } 00787 } 00788 return(otot); 00789 } 00790 00791 /******************************************************************************/ 00792 00793 /*OpacityCreate1Element generate ionic subshell opacities by calling t_ADfA::Inst().phfit */ 00794 STATIC void OpacityCreate1Element( 00795 /* atomic number on the C scale, lowest ever called will be Li=2 */ 00796 long int nelem) 00797 { 00798 long int ihi, 00799 ip, 00800 ipop, 00801 limit, 00802 low, 00803 need, 00804 nelec, 00805 ion, 00806 nshell; 00807 double cs; 00808 double energy; 00809 00810 DEBUG_ENTRY( "OpacityCreate1Element()" ); 00811 00812 /* confirm range of validity of atomic number, Li=2 should be the lightest */ 00813 ASSERT( nelem >= 2 ); 00814 ASSERT( nelem < LIMELM ); 00815 00816 /*>>chng 99 jan 27, no longer redo hydrogenic opacity here */ 00817 /* for( ion=0; ion <= nelem; ion++ )*/ 00818 for( ion=0; ion < nelem; ion++ ) 00819 { 00820 00821 /* will be used for a sanity check on number of hits in a cell*/ 00822 for( ip=0; ip < rfield.nupper; ip++ ) 00823 { 00824 opac.opacity_abs[ip] = 0.; 00825 } 00826 00827 /* number of bound electrons */ 00828 nelec = nelem+1 - ion; 00829 00830 /* loop over all shells, from innermost K shell to valence */ 00831 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ ) 00832 { 00833 /* this is array index for start of this shell within large opacity stack */ 00834 opac.ipElement[nelem][ion][nshell][2] = opac.nOpacTot + 1; 00835 00836 /* this is continuum upper limit to array index for energy range of this shell */ 00837 limit = opac.ipElement[nelem][ion][nshell][1]; 00838 00839 /* this is number of cells in continuum needed to store opacity */ 00840 need = limit - opac.ipElement[nelem][ion][nshell][0] + 1; 00841 00842 /* check that opac will have enough frequeny cells */ 00843 if( opac.nOpacTot + need > ndimOpacityStack ) 00844 opacity_more_memory(); 00845 00846 /* set lower and upper limits to this range */ 00847 low = opac.ipElement[nelem][ion][nshell][0]; 00848 ihi = opac.ipElement[nelem][ion][nshell][1]; 00849 ipop = opac.ipElement[nelem][ion][nshell][2]; 00850 00851 /* make sure indices are within correct bounds, 00852 * mainly check on logic for detecting missing shells */ 00853 ASSERT( low <= ihi || low<5 ); 00854 00855 /* loop over energy range of this shell */ 00856 for( ip=low-1; ip < ihi; ip++ ) 00857 { 00858 /* photo energy MAX so that we never eval below threshold */ 00859 energy = MAX2(rfield.AnuOrg[ip]*EVRYD , 00860 t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)); 00861 00862 /* the cross section in mega barns */ 00863 cs = t_ADfA::Inst().phfit(nelem+1,nelec,nshell+1,energy); 00864 /* cannot assert that cs is positive since, at edge of shell, 00865 * energy might be slightly below threshold and hence zero, 00866 * due to finite size of continuum bins */ 00867 opac.OpacStack[ip-low+ipop] = cs*1e-18; 00868 00869 /* add this to total opacity, which we will confirm to be greater than zero below */ 00870 opac.opacity_abs[ip] += cs; 00871 } 00872 00873 opac.nOpacTot += ihi - low + 1; 00874 00875 /* punch pointers option */ 00876 if( punch.lgPunPoint ) 00877 { 00878 fprintf( punch.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n", 00879 nelem, ion, nshell, rfield.anu[low-1], rfield.anu[ihi-1], 00880 opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] ); 00881 } 00882 } 00883 00884 ASSERT( Heavy.nsShells[nelem][ion] >= 1 ); 00885 /*confirm that total opacity is greater than zero */ 00886 for( ip=opac.ipElement[nelem][ion][Heavy.nsShells[nelem][ion]-1][0]-1; 00887 ip < continuum.KshellLimit; ip++ ) 00888 { 00889 ASSERT( opac.opacity_abs[ip] > 0. ); 00890 } 00891 00892 } 00893 return; 00894 } 00895 00896 /*opacity_more_memory allocate more memory for opacity stack */ 00897 STATIC void opacity_more_memory(void) 00898 { 00899 00900 DEBUG_ENTRY( "opacity_more_memory()" ); 00901 00902 /* double size */ 00903 ndimOpacityStack *= 2; 00904 opac.OpacStack = (double *)REALLOC( opac.OpacStack , (size_t)ndimOpacityStack*sizeof(double) ); 00905 fprintf( ioQQQ, " NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack, please consider increasing it.\n" ); 00906 fprintf( ioQQQ, " NOTE OpacityCreate1Element doubled memory allocation to %li.\n",ndimOpacityStack ); 00907 lgRealloc = true; 00908 return; 00909 } 00910 00911 /*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */ 00912 STATIC double Opacity_iso_photo_cs( 00913 /* photon energy ryd */ 00914 realnum energy , 00915 /* iso sequence */ 00916 long ipISO , 00917 /* charge, 0 for H */ 00918 long nelem , 00919 /* n - meaning depends on iso */ 00920 long n ) 00921 { 00922 /* >>chng 01 dec 23, from float to double */ 00923 double thres/*,ejected_electron_energy*/; 00924 double crs=-DBL_MAX; 00925 00926 DEBUG_ENTRY( "Opacity_iso_photo_cs()" ); 00927 00928 if( ipISO==ipH_LIKE ) 00929 { 00930 double photon; 00931 00932 if( n==0 ) 00933 { 00934 /* this is the ground state, use Dima's routine, which works in eV 00935 * and returns megabarns */ 00936 thres = MAX2(energy*(realnum)EVRYD , t_ADfA::Inst().ph1(0,0,nelem,0)); 00937 crs = t_ADfA::Inst().phfit(nelem+1,1,1,thres)* 1e-18; 00938 /* make sure cross section is reasonable */ 00939 ASSERT( crs > 0. && crs < 1e-10 ); 00940 } 00941 else if( n < iso.numLevels_max[ipISO][nelem] - iso.nCollapsed_max[ipISO][nelem] ) 00942 { 00943 /* photon energy relative to threshold */ 00944 photon = MAX2( energy/iso.xIsoLevNIonRyd[ipISO][nelem][n], 1. + FLT_EPSILON*2. ); 00945 00946 crs = H_photo_cs( photon , N_(n), L_(n), nelem+1 ); 00947 /* make sure cross section is reasonable */ 00948 ASSERT( crs > 0. && crs < 1e-10 ); 00949 } 00950 else if( N_(n) <= NHYDRO_MAX_LEVEL ) 00951 { 00952 /* for first cell, depending on the current resolution of the energy mesh, 00953 * the center of the first cell can be below the ionization limit of the 00954 * level. do not let the energy fall below this limit */ 00955 /* This will make sure that we don't call epsilon below threshold, 00956 * the factor 1.001 was chosen so that t_ADfA::Inst().hpfit, which works 00957 * in terms of Dima's Rydberg constant, is not tripped below threshold */ 00958 thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.001f ); 00959 00960 crs = t_ADfA::Inst().hpfit(nelem+1,N_(n),thres*EVRYD); 00961 /* make sure cross section is reasonable */ 00962 ASSERT( crs > 0. && crs < 1e-10 ); 00963 } 00964 else 00965 { 00966 /* photon energy relative to threshold */ 00967 photon = MAX2( energy/iso.xIsoLevNIonRyd[ipISO][nelem][n], 1. + FLT_EPSILON*2. ); 00968 00969 /* cross section for collapsed level should be 00970 * roughly equal to cross-section for yrast level, 00971 * so third parameter is n - 1. */ 00972 crs = H_photo_cs( photon , N_(n), N_(n)-1, nelem+1 ); 00973 00974 /* make sure cross section is reasonable */ 00975 ASSERT( crs > 0. && crs < 1e-10 ); 00976 } 00977 } 00978 else if( ipISO==ipHE_LIKE ) 00979 { 00980 thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]); 00981 /* this would be a collapsed level */ 00982 if( n >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) 00983 { 00984 long int nup = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + n + 1 - 00985 (iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem]); 00986 00987 /* this is a collapsed level - this is hydrogenic routine and 00988 * first he-like energy may not agree exactly with threshold for H */ 00989 crs = t_ADfA::Inst().hpfit(nelem,nup ,thres*EVRYD); 00990 /* make sure cross section is reasonable if away from threshold */ 00991 ASSERT( 00992 (energy < iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.02) || 00993 (crs > 0. && crs < 1e-10) ); 00994 } 00995 else 00996 { 00997 00998 /* He_cross_section returns cross section (cm^-2), 00999 * given EgammaRyd, the photon energy in Ryd, 01000 * ipLevel, the index of the level, 0 is ground, 3 within 2 3P, 01001 * nelem is charge, equal to 1 for Helium, 01002 * this is a wrapper for cross_section */ 01003 crs = He_cross_section( thres , n , nelem ); 01004 01005 /* make sure cross section is reasonable */ 01006 ASSERT( crs > 0. && crs < 1e-10 ); 01007 } 01008 } 01009 else 01010 TotalInsanity(); 01011 return(crs); 01012 } 01013 01014 /*hmiopc derive total H- H minus opacity */ 01015 static const int NCRS = 33; 01016 01017 STATIC double hmiopc(double freq) 01018 { 01019 double energy, 01020 hmiopc_v, 01021 x, 01022 y; 01023 static double y2[NCRS]; 01024 static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805, 01025 2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925, 01026 3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898, 01027 1.567,1.233,.912,.629,.39,.19}; 01028 static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351, 01029 0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129, 01030 0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470, 01031 0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001, 01032 0.5520,0.8557,1.7669}; 01033 static bool lgFirst = true; 01034 01035 DEBUG_ENTRY( "hmiopc()" ); 01036 01037 /* bound free cross section (x10**-17 cm^2) from Doughty et al 01038 * 1966, MNRAS 132, 255; good agreement with Wishart MNRAS 187, 59p. */ 01039 01040 /* photoelectron energy, add 0.05552 to get incoming energy (Ryd) */ 01041 01042 01043 if( lgFirst ) 01044 { 01045 /* set up coefficients for spline */ 01046 spline(ener,crs,NCRS,2e31,2e31,y2); 01047 lgFirst = false; 01048 } 01049 01050 energy = freq - 0.05552; 01051 if( energy < ener[0] || energy > ener[NCRS-1] ) 01052 { 01053 hmiopc_v = 0.; 01054 } 01055 else 01056 { 01057 x = energy; 01058 splint(ener,crs,y2,NCRS,x,&y); 01059 hmiopc_v = y*1e-17; 01060 } 01061 return( hmiopc_v ); 01062 } 01063 01064 /*rayleh compute Rayleigh scattering cross section for Lya */ 01065 STATIC double rayleh(double ener) 01066 { 01067 double rayleh_v; 01068 01069 DEBUG_ENTRY( "rayleh()" ); 01070 01072 /* do hydrogen Rayleigh scattering cross sections; 01073 * fits to 01074 *>>refer Ly scattering Gavrila, M., 1967, Physical Review 163, 147 01075 * and Mihalas radiative damping 01076 * 01077 * >>chng 96 aug 15, changed logic to do more terms for each part of 01078 * rayleigh scattering 01079 * if( ener.lt.0.05 ) then 01080 * rayleh = 8.41e-25 * ener**4 * DampOnFac 01081 * */ 01082 if( ener < 0.05 ) 01083 { 01084 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))* 01085 hydro.DampOnFac; 01086 } 01087 01088 else if( ener < 0.646 ) 01089 { 01090 rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) + 01091 4.71e-22*powi(ener,14))*hydro.DampOnFac; 01092 } 01093 01094 else if( ener >= 0.646 && ener < 1.0 ) 01095 { 01096 rayleh_v = fabs(0.74959-ener); 01097 rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v)); 01098 /* typical energy between Ly-a and Ly-beta */ 01099 rayleh_v = MAX2(rayleh_v,1e-24)*hydro.DampOnFac; 01100 } 01101 01102 else 01103 { 01104 rayleh_v = 0.; 01105 } 01106 return( rayleh_v ); 01107 }