cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
init_coreload.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*InitCoreload one time initialization of core load, called from cdDrive, this sets
4 * minimum set of values needed for the code to start - called after
5 * input lines have been read in and checked for VARY or GRID - so
6 * known whether single or multiple sims will be run */
7 #include "cddefines.h"
8 #include "optimize.h"
9 #include "parse.h"
10 #include "dense.h"
11 #include "hcmap.h"
12 #include "h2.h"
13 #include "version.h"
14 #include "hextra.h"
15 #include "mole.h"
16 #include "extinc.h"
17 #include "heavy.h"
18 #include "grid.h"
19 #include "ionbal.h"
20 #include "iso.h"
21 #include "taulines.h"
22 /* */
23 /* following used to enter total predicted continuum into emission line array
24 * actually entered into line array in lineset1
25 */
26 #include "predcont.h"
27 /* energies where diffuse continuum is to be entered into line array
28 * variables are declared in predcont.h, where array dimensioned
29 * NB - if these numbers change, the wavelength in the printout will change
30 * too, since the wavelength is derived form this */
31 /* >>>chng 99 mar 23, adjusted energies so that wavelength line list is
32 * the same as it was in C90 - small changes were caused by going over
33 * to proper Rydberg constant */
35  /* this is radio continuum at 3.4 cm */
36  2.680e-06f
37  , 7.445e-04f
38  , 1.498e-03f
39  , 2.211e-03f
40  , 2.952e-03f
41  , 3.677e-03f
42  , 3.7501e-03f /* Ney-Allen */
43  , 3.9915e-03f /* Ney-Allen */
44  , 4.2543e-03f /* Ney-Allen */
45  , 4.314e-03f
46  , 4.6446e-03f /* Ney-Allen */
47  , 5.162e-03f
48  , 5.2462e-03f /* Ney-Allen */
49  , 5.8079e-03f /* Ney-Allen */
50  , 6.240e-03f
51  /* , 7.320e-03 */
52  , 7.3312e-03f /* Ney-Allen */
53  , 7.9936e-03f /* Ney-Allen */
54  , 8.7119e-03f /* Ney-Allen */
55  , 9.6125e-03f /* Ney-Allen */
56  , 9.77243e-03f
57  , 1.1099e-02f /* Ney-Allen */
58  , 1.2022e-02f /* Ney-Allen */
59  , 1.29253e-02f
60  , 2.2152e-02f
61  , 3.92044e-02f
62  , 5.54593e-02f
63  /* next two either side of n=4 edge of hydrogen, set to 1.5% off either direction*/
64  /* >>chng 00 sep 18, had been too close in energy */
65  , 6.1563e-02f
66  , 6.3437e-02f
67  , 8.1460e-02f
68  /* >>chng 00 sep 14, changed energies of paschen jump to be farther away as
69  * per note on BJ */
70  , 0.1094f
71  , 0.1128f
72  , 0.14675f
73  , 0.18653f
74  /* >>chng 00 sep 14, next two energies changed since they were too close to BJ
75  * and so both ended up shortward of limit*/
76  , 0.246f /* these two are the Balmer jump, below and above. */
77  , 0.254f /* continuum binning not much better than 1% so offset energies by more */
78  , 0.375f /* peak on two photon continuum */
79  , 0.38096f
80  , 0.43994f
81  , 0.44394f
82  , 0.50811f
83  , 0.57489f
84  , 0.62487f
85  , 0.67155f
86  , 0.70244f
87  , 0.72163f
88  , 0.74812f
89  , 0.76172f
90  , 0.77551f
91  , 0.79681f
92  , 0.81859f
93  , 0.8260f
94  , 0.84859f
95  , 0.85618f
96  , 0.87967f
97  , 0.911267f /*exactly 1000A */
98  /* points on either side of Lyman jump,
99  * energies changed to be robust when energy grid changes,
100  * grid resolution is about 1%, so change from 0.99783 and 1.000
101  * to 1 +/- 1.5%
102  * >>chng 00 sep 23 change wavelength points for next two */
103  , 0.985f
104  , 1.015f
105  , 1.199f
106  , 1.299f
107  , 1.4984f
108  , 1.58441f
109  /* points on either side of Lyman jump,
110  * energies changed to be robust when energy grid changes,
111  * grid resolution is about 1%, so change from 1.80433 and 1.809
112  * to 1.807 +/- 1.5%
113  * >>chng 00 sep 23 change wavelength points for next two */
114  , 1.780f
115  , 1.834f
116  , 2.283f
117 };
118 
120 
121 /*InitCoreload one time initialization of core load, called from cdDrive, this sets
122 * minimum set of values needed for the code to start - called after
123 * input lines have been read in and checked for VARY or GRID - so
124 * known whether single or multiple sims will be run */
125 void InitCoreload( void )
126 {
127  static int nCalled=0;
128  long int nelem;
129 
130  DEBUG_ENTRY( "InitCoreload()" );
131 
132  /* sanity check */
133  if( nCalled )
134  {
135  /* return since already called */
136  return;
137  }
138  ++nCalled;
139 
140  /* number of simulation sin this coreload,
141  * 0 while in this routine, 1 for first sim, increments just after
142  * this routine is called */
143  nSimThisCoreload = 0;
144 
145  strncpy( chOptimFileName , "optimal.in" , sizeof( chOptimFileName ) );
146 
147  /* the constant that multiplies the column density to get optical depth at 1 Ryd */
148  extinc.cnst_col2optdepth = 6.22e-18f;
149  /* the power on the energy */
150  extinc.cnst_power = -2.43f;
151 
152  hcmap.lgMapOK = true;
153  hcmap.lgMapDone = false;
154 
155  /* will be reset to positive number when map actually done */
156  hcmap.nMapAlloc = 0;
157  hcmap.nmap = 0;
158  hcmap.lgMapBeingDone = false;
159 
160  /* name of output file from optimization run */
161  strncpy( chOptimFileName , "optimal.in" , sizeof( chOptimFileName ) );
162 
163  /* number of electrons in valence shell that can Compton recoil ionize */
164  long int nCom[LIMELM] =
165  {
166  1 , 2 , /* K 1s shell */
167  1 , 2 , /* L 2s shell */
168  1 , 2 , 3 , 4 , 5 , 6 , /* L 2p shell */
169  1 , 2 , /* M 3s shell */
170  1 , 2 , 3 , 4 , 5 , 6 , /* M 3p shell */
171  1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , /* M 3d shell */
172  1 , 2 , /* N 4s shell */
173  1 , 2 /* N 4p shell */
174  };
175 
176  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
177  {
178  ionbal.nCompRecoilElec[nelem] = nCom[nelem];
179  }
180 
181  /* list of shells, 1 to 7 */
182  strcpy( Heavy.chShell[0], "1s" );
183  strcpy( Heavy.chShell[1], "2s" );
184  strcpy( Heavy.chShell[2], "2p" );
185  strcpy( Heavy.chShell[3], "3s" );
186  strcpy( Heavy.chShell[4], "3p" );
187  strcpy( Heavy.chShell[5], "3d" );
188  strcpy( Heavy.chShell[6], "4s" );
189 
190  /* variables for H-like sequence */
191  /* default number of levels for hydrogen iso sequence */
192  for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
193  {
194  iso.n_HighestResolved_max[ipH_LIKE][nelem] = 5;
195  iso.nCollapsed_max[ipH_LIKE][nelem] = 2;
196  }
197 
200 
203 
204  /* variables for He-like sequence */
205  /* "he-like" hydrogen (H-) is treated elsewhere */
207  iso.numLevels_max[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX;
208  iso.numPrintLevels[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX;
209  iso.nCollapsed_max[ipHE_LIKE][ipHYDROGEN] = -LONG_MAX;
210 
211  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
212  {
213  /* put at least three resolved and 1 collapsed in every element for he-like */
214  iso.n_HighestResolved_max[ipHE_LIKE][nelem] = 3;
215  iso.nCollapsed_max[ipHE_LIKE][nelem] = 1;
216  }
217 
220 
221  /* And n=5 for these because they are most abundant */
230  /* also set this, for exercising any possible issues with biggest charge models */
231  iso.n_HighestResolved_max[ipHE_LIKE][LIMELM-1] = 5;
232 
233  iso.chISO[ipH_LIKE] = "H-like ";
234  iso.chISO[ipHE_LIKE] = "He-like";
235 
236  max_num_levels = 0;
237  for( long ipISO = ipH_LIKE; ipISO < NISO; ipISO++ )
238  {
239  for( nelem=ipISO; nelem < LIMELM; nelem++ )
240  {
241  /* set this to LONG_MAX, reduce to actual number later,
242  * then verify number of levels is not increased after initial coreload */
243  iso.numLevels_malloc[ipISO][nelem] = LONG_MAX;
244  iso_update_num_levels( ipISO, nelem );
245  }
246  }
247 
248  statesAdded = 0;
249  lgStatesAdded = false;
250  linesAdded = 0;
251  lgLinesAdded = false;
252 
254  /* these are used to set trace levels of output by options on atom h2 trace command
255  * first is minimum level of trace, keyword is FINAL */
256  mole.nH2_trace_final = 1;
257  /* intermediate level of trace, info per iteration, key ITERATION */
259  /* full trace, keyword is FULL */
260  mole.nH2_trace_full = 3;
261  /* print matrices used in solving X */
263 
264  /* turn element on if this is first call to this routine,
265  * but do not reset otherwise since would clobber how elements were set up */
266  for( nelem=0; nelem < LIMELM; nelem++ )
267  {
268  /* never turn on element if turned off on first iteration */
269  dense.lgElmtOn[nelem] = true;
270 
271  /* option to set ionization with element ioniz cmnd,
272  * default (false) is to solve for ionization */
273  dense.lgSetIoniz[nelem] = false;
274  }
275 
276  /* smallest density to permit in any ion - if smaller then set to zero */
278 
279  /* default cosmic ray background */
280  /* >>chng 99 jun 24, slight change to value
281  * quoted by
282  * >>refer cosmic ray ionization rate McKee, C.M., 1999, astro-ph 9901370
283  * this will produce a total
284  * secondary ionization rate of 2.5e-17 s^-1, as tested in
285  * test suite cosmicray.in. If each ionization produces 2.4 eV of heat,
286  * the background heating rate should be 9.6e-29 * n*/
287  /* >>chng 04 jan 26, update cosmic ray ionization rate for H0 to
288  * >>refer cosmic ray ionization Williams, J.P., Bergin, E.A., Caselli, P.,
289  * >>refercon Myers, P.C., & Plume, R. 1998, ApJ, 503, 689,
290  * H0 ionization rate of 2.5e-17 s-1 and a H2 ionization rate twice this
291  * >>chng 04 mar 15, comment said 2.5e-17 which is correct, but code produce 8e-17,
292  * fix back to correct value
293  */
294  /* NB - the rate is derived from the density. these two are related by the secondary
295  * ionization efficiency problem. background_rate is only here to provide the relationship
296  * for predominantly neutral gas. the background_density is the real rate.
297  hextra.background_density = 1.99e-9f;*/
298  /* >>chng 05 apr 16, to get proper ionization rate in ism_set_cr_rate, where
299  * H is forced to be fully atomic, no molecules, density from 1.99 to 2.15 */
300  hextra.background_density = 2.15e-9f;
301  hextra.background_rate = 2.5e-17f;
302 
303  /* initialization for punch files - must call after input commands have
304  * been scanned for grid and vary options. So known if grid or single run
305  * default punch is different for these */
306  grid.lgGridDone = false;
307  grid.lgStrictRepeat = false;
308 
309  PunchFilesInit();
310 
312 
313  return;
314 }

Generated for cloudy by doxygen 1.8.4