cloudy: /home66/gary/public_html/cloudy/c08_branch/source/init

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 }
/home66/gary/public_html/cloudy/c08_branch/source/init_coreload.cpp
