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 /*InitCoreload one time initialization of core load, called from cdDrive, this sets 00004 * minimum set of values needed for the code to start - called after 00005 * input lines have been read in and checked for VARY or GRID - so 00006 * known whether single or multiple sims will be run */ 00007 #include "cddefines.h" 00008 #include "optimize.h" 00009 #include "parse.h" 00010 #include "dense.h" 00011 #include "hcmap.h" 00012 #include "h2.h" 00013 #include "version.h" 00014 #include "hextra.h" 00015 #include "mole.h" 00016 #include "extinc.h" 00017 #include "heavy.h" 00018 #include "grid.h" 00019 #include "ionbal.h" 00020 #include "iso.h" 00021 #include "taulines.h" 00022 /* */ 00023 /* following used to enter total predicted continuum into emission line array 00024 * actually entered into line array in lineset1 00025 */ 00026 #include "predcont.h" 00027 /* energies where diffuse continuum is to be entered into line array 00028 * variables are declared in predcont.h, where array dimensioned 00029 * NB - if these numbers change, the wavelength in the printout will change 00030 * too, since the wavelength is derived form this */ 00031 /* >>>chng 99 mar 23, adjusted energies so that wavelength line list is 00032 * the same as it was in C90 - small changes were caused by going over 00033 * to proper Rydberg constant */ 00034 realnum EnrPredCont[NPREDCONT] = { 00035 /* this is radio continuum at 3.4 cm */ 00036 2.680e-06f 00037 , 7.445e-04f 00038 , 1.498e-03f 00039 , 2.211e-03f 00040 , 2.952e-03f 00041 , 3.677e-03f 00042 , 3.7501e-03f /* Ney-Allen */ 00043 , 3.9915e-03f /* Ney-Allen */ 00044 , 4.2543e-03f /* Ney-Allen */ 00045 , 4.314e-03f 00046 , 4.6446e-03f /* Ney-Allen */ 00047 , 5.162e-03f 00048 , 5.2462e-03f /* Ney-Allen */ 00049 , 5.8079e-03f /* Ney-Allen */ 00050 , 6.240e-03f 00051 /* , 7.320e-03 */ 00052 , 7.3312e-03f /* Ney-Allen */ 00053 , 7.9936e-03f /* Ney-Allen */ 00054 , 8.7119e-03f /* Ney-Allen */ 00055 , 9.6125e-03f /* Ney-Allen */ 00056 , 9.77243e-03f 00057 , 1.1099e-02f /* Ney-Allen */ 00058 , 1.2022e-02f /* Ney-Allen */ 00059 , 1.29253e-02f 00060 , 2.2152e-02f 00061 , 3.92044e-02f 00062 , 5.54593e-02f 00063 /* next two either side of n=4 edge of hydrogen, set to 1.5% off either direction*/ 00064 /* >>chng 00 sep 18, had been too close in energy */ 00065 , 6.1563e-02f 00066 , 6.3437e-02f 00067 , 8.1460e-02f 00068 /* >>chng 00 sep 14, changed energies of paschen jump to be farther away as 00069 * per note on BJ */ 00070 , 0.1094f 00071 , 0.1128f 00072 , 0.14675f 00073 , 0.18653f 00074 /* >>chng 00 sep 14, next two energies changed since they were too close to BJ 00075 * and so both ended up shortward of limit*/ 00076 , 0.246f /* these two are the Balmer jump, below and above. */ 00077 , 0.254f /* continuum binning not much better than 1% so offset energies by more */ 00078 , 0.375f /* peak on two photon continuum */ 00079 , 0.38096f 00080 , 0.43994f 00081 , 0.44394f 00082 , 0.50811f 00083 , 0.57489f 00084 , 0.62487f 00085 , 0.67155f 00086 , 0.70244f 00087 , 0.72163f 00088 , 0.74812f 00089 , 0.76172f 00090 , 0.77551f 00091 , 0.79681f 00092 , 0.81859f 00093 , 0.8260f 00094 , 0.84859f 00095 , 0.85618f 00096 , 0.87967f 00097 , 0.911267f /*exactly 1000A */ 00098 /* points on either side of Lyman jump, 00099 * energies changed to be robust when energy grid changes, 00100 * grid resolution is about 1%, so change from 0.99783 and 1.000 00101 * to 1 +/- 1.5% 00102 * >>chng 00 sep 23 change wavelength points for next two */ 00103 , 0.985f 00104 , 1.015f 00105 , 1.199f 00106 , 1.299f 00107 , 1.4984f 00108 , 1.58441f 00109 /* points on either side of Lyman jump, 00110 * energies changed to be robust when energy grid changes, 00111 * grid resolution is about 1%, so change from 1.80433 and 1.809 00112 * to 1.807 +/- 1.5% 00113 * >>chng 00 sep 23 change wavelength points for next two */ 00114 , 1.780f 00115 , 1.834f 00116 , 2.283f 00117 }; 00118 00119 long int ipPredCont[NPREDCONT]; 00120 00121 /*InitCoreload one time initialization of core load, called from cdDrive, this sets 00122 * minimum set of values needed for the code to start - called after 00123 * input lines have been read in and checked for VARY or GRID - so 00124 * known whether single or multiple sims will be run */ 00125 void InitCoreload( void ) 00126 { 00127 static int nCalled=0; 00128 long int nelem; 00129 00130 DEBUG_ENTRY( "InitCoreload()" ); 00131 00132 /* sanity check */ 00133 if( nCalled ) 00134 { 00135 /* return since already called */ 00136 return; 00137 } 00138 ++nCalled; 00139 00140 /* number of simulation sin this coreload, 00141 * 0 while in this routine, 1 for first sim, increments just after 00142 * this routine is called */ 00143 nSimThisCoreload = 0; 00144 00145 strncpy( chOptimFileName , "optimal.in" , sizeof( chOptimFileName ) ); 00146 00147 /* the constant that multiplies the column density to get optical depth at 1 Ryd */ 00148 extinc.cnst_col2optdepth = 6.22e-18f; 00149 /* the power on the energy */ 00150 extinc.cnst_power = -2.43f; 00151 00152 hcmap.lgMapOK = true; 00153 hcmap.lgMapDone = false; 00154 00155 /* will be reset to positive number when map actually done */ 00156 hcmap.nMapAlloc = 0; 00157 hcmap.nmap = 0; 00158 hcmap.lgMapBeingDone = false; 00159 00160 /* name of output file from optimization run */ 00161 strncpy( chOptimFileName , "optimal.in" , sizeof( chOptimFileName ) ); 00162 00163 /* number of electrons in valence shell that can Compton recoil ionize */ 00164 long int nCom[LIMELM] = 00165 { 00166 1 , 2 , /* K 1s shell */ 00167 1 , 2 , /* L 2s shell */ 00168 1 , 2 , 3 , 4 , 5 , 6 , /* L 2p shell */ 00169 1 , 2 , /* M 3s shell */ 00170 1 , 2 , 3 , 4 , 5 , 6 , /* M 3p shell */ 00171 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , /* M 3d shell */ 00172 1 , 2 , /* N 4s shell */ 00173 1 , 2 /* N 4p shell */ 00174 }; 00175 00176 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00177 { 00178 ionbal.nCompRecoilElec[nelem] = nCom[nelem]; 00179 } 00180 00181 /* list of shells, 1 to 7 */ 00182 strcpy( Heavy.chShell[0], "1s" ); 00183 strcpy( Heavy.chShell[1], "2s" ); 00184 strcpy( Heavy.chShell[2], "2p" ); 00185 strcpy( Heavy.chShell[3], "3s" ); 00186 strcpy( Heavy.chShell[4], "3p" ); 00187 strcpy( Heavy.chShell[5], "3d" ); 00188 strcpy( Heavy.chShell[6], "4s" ); 00189 00190 /* variables for H-like sequence */ 00191 /* default number of levels for hydrogen iso sequence */ 00192 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem ) 00193 { 00194 iso.n_HighestResolved_max[ipH_LIKE][nelem] = 5; 00195 iso.nCollapsed_max[ipH_LIKE][nelem] = 2; 00196 } 00197 00198 iso.n_HighestResolved_max[ipH_LIKE][ipHYDROGEN] = 10; 00199 iso.n_HighestResolved_max[ipH_LIKE][ipHELIUM] = 10; 00200 00201 iso.nCollapsed_max[ipH_LIKE][ipHYDROGEN] = 15; 00202 iso.nCollapsed_max[ipH_LIKE][ipHELIUM] = 15; 00203 00204 /* variables for He-like sequence */ 00205 /* "he-like" hydrogen (H-) is treated elsewhere */ 00206 iso.n_HighestResolved_max[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX; 00207 iso.numLevels_max[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX; 00208 iso.numPrintLevels[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX; 00209 iso.nCollapsed_max[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX; 00210 00211 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem ) 00212 { 00213 /* put at least three resolved and 1 collapsed in every element for he-like */ 00214 iso.n_HighestResolved_max[ipHE_LIKE][nelem] = 3; 00215 iso.nCollapsed_max[ipHE_LIKE][nelem] = 1; 00216 } 00217 00218 iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] = 6; 00219 iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] = 20; 00220 00221 /* And n=5 for these because they are most abundant */ 00222 iso.n_HighestResolved_max[ipHE_LIKE][ipCARBON] = 5; 00223 iso.n_HighestResolved_max[ipHE_LIKE][ipNITROGEN] = 5; 00224 iso.n_HighestResolved_max[ipHE_LIKE][ipOXYGEN] = 5; 00225 iso.n_HighestResolved_max[ipHE_LIKE][ipNEON] = 5; 00226 iso.n_HighestResolved_max[ipHE_LIKE][ipSILICON] = 5; 00227 iso.n_HighestResolved_max[ipHE_LIKE][ipMAGNESIUM] = 5; 00228 iso.n_HighestResolved_max[ipHE_LIKE][ipSULPHUR] = 5; 00229 iso.n_HighestResolved_max[ipHE_LIKE][ipIRON] = 5; 00230 /* also set this, for exercising any possible issues with biggest charge models */ 00231 iso.n_HighestResolved_max[ipHE_LIKE][LIMELM-1] = 5; 00232 00233 iso.chISO[ipH_LIKE] = "H-like "; 00234 iso.chISO[ipHE_LIKE] = "He-like"; 00235 00236 max_num_levels = 0; 00237 for( long ipISO = ipH_LIKE; ipISO < NISO; ipISO++ ) 00238 { 00239 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00240 { 00241 /* set this to LONG_MAX, reduce to actual number later, 00242 * then verify number of levels is not increased after initial coreload */ 00243 iso.numLevels_malloc[ipISO][nelem] = LONG_MAX; 00244 iso_update_num_levels( ipISO, nelem ); 00245 } 00246 } 00247 00248 statesAdded = 0; 00249 lgStatesAdded = false; 00250 linesAdded = 0; 00251 lgLinesAdded = false; 00252 00254 /* these are used to set trace levels of output by options on atom h2 trace command 00255 * first is minimum level of trace, keyword is FINAL */ 00256 mole.nH2_trace_final = 1; 00257 /* intermediate level of trace, info per iteration, key ITERATION */ 00258 mole.nH2_trace_iterations = 2; 00259 /* full trace, keyword is FULL */ 00260 mole.nH2_trace_full = 3; 00261 /* print matrices used in solving X */ 00262 mole.nH2_trace_matrix = 4; 00263 00264 /* turn element on if this is first call to this routine, 00265 * but do not reset otherwise since would clobber how elements were set up */ 00266 for( nelem=0; nelem < LIMELM; nelem++ ) 00267 { 00268 /* never turn on element if turned off on first iteration */ 00269 dense.lgElmtOn[nelem] = true; 00270 00271 /* option to set ionization with element ioniz cmnd, 00272 * default (false) is to solve for ionization */ 00273 dense.lgSetIoniz[nelem] = false; 00274 } 00275 00276 /* smallest density to permit in any ion - if smaller then set to zero */ 00277 dense.density_low_limit = SMALLFLOAT * 1e3; 00278 00279 /* default cosmic ray background */ 00280 /* >>chng 99 jun 24, slight change to value 00281 * quoted by 00282 * >>refer cosmic ray ionization rate McKee, C.M., 1999, astro-ph 9901370 00283 * this will produce a total 00284 * secondary ionization rate of 2.5e-17 s^-1, as tested in 00285 * test suite cosmicray.in. If each ionization produces 2.4 eV of heat, 00286 * the background heating rate should be 9.6e-29 * n*/ 00287 /* >>chng 04 jan 26, update cosmic ray ionization rate for H0 to 00288 * >>refer cosmic ray ionization Williams, J.P., Bergin, E.A., Caselli, P., 00289 * >>refercon Myers, P.C., & Plume, R. 1998, ApJ, 503, 689, 00290 * H0 ionization rate of 2.5e-17 s-1 and a H2 ionization rate twice this 00291 * >>chng 04 mar 15, comment said 2.5e-17 which is correct, but code produce 8e-17, 00292 * fix back to correct value 00293 */ 00294 /* NB - the rate is derived from the density. these two are related by the secondary 00295 * ionization efficiency problem. background_rate is only here to provide the relationship 00296 * for predominantly neutral gas. the background_density is the real rate. 00297 hextra.background_density = 1.99e-9f;*/ 00298 /* >>chng 05 apr 16, to get proper ionization rate in ism_set_cr_rate, where 00299 * H is forced to be fully atomic, no molecules, density from 1.99 to 2.15 */ 00300 hextra.background_density = 2.15e-9f; 00301 hextra.background_rate = 2.5e-17f; 00302 00303 /* initialization for punch files - must call after input commands have 00304 * been scanned for grid and vary options. So known if grid or single run 00305 * default punch is different for these */ 00306 grid.lgGridDone = false; 00307 grid.lgStrictRepeat = false; 00308 00309 PunchFilesInit(); 00310 00311 H2_init_coreload(); 00312 00313 return; 00314 }