cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
zero.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 /*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization,
4  * this is called before any commands are parsed, called one time per model, at very start */
5 /*rfield_optac_zero zero out rfield arrays between certain limits */
6 #include "cddefines.h"
7 #include "physconst.h"
8 #include "iterations.h"
9 #include "hydrogenic.h"
10 #include "oxy.h"
11 #include "doppvel.h"
12 #include "dense.h"
13 #include "hextra.h"
14 #include "grains.h"
15 #include "magnetic.h"
16 #include "state.h"
17 #include "rt.h"
18 #include "he.h"
19 #include "struc.h"
20 #include "h2.h"
21 #include "coolheavy.h"
22 #include "lines.h"
23 #include "dynamics.h"
24 #include "carb.h"
25 #include "mean.h"
26 #include "atomfeii.h"
27 #include "iso.h"
28 #include "conv.h"
29 #include "geometry.h"
30 #include "timesc.h"
31 #include "peimbt.h"
32 #include "ionbal.h"
33 #include "continuum.h"
34 #include "atmdat.h"
35 #include "mole.h"
36 #include "ca.h"
37 #include "input.h"
38 #include "atoms.h"
39 #include "pressure.h"
40 #include "numderiv.h"
41 #include "colden.h"
42 #include "yield.h"
43 #include "hmi.h"
44 #include "rfield.h"
45 #include "abund.h"
46 #include "radius.h"
47 #include "opacity.h"
48 #include "broke.h"
49 #include "secondaries.h"
50 #include "called.h"
51 #include "phycon.h"
52 #include "warnings.h"
53 #include "thermal.h"
54 #include "cooling.h"
55 #include "fe.h"
56 #include "hyperfine.h"
57 #include "init.h"
58 
59 /*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization,
60  * this is called before any commands are parsed, called one time per model, at very start */
61 void zero(void)
62 {
63  long int i,
64  ion,
65  ipISO ,
66  nelem,
67  ns;
68 
69  /* this is used to signify the first call to this routine. At that
70  * stage some memory has not been allocated so must not initialize,
71  * set false at very end of this routine */
72  static bool lgFirstCall = true;
73 
74  DEBUG_ENTRY( "zero()" );
75 
76  /* this routine is called exactly one time at the start of
77  * the calculation of a single model. When the code is used as a subroutine
78  * this routine is called one time for each model. It is called before
79  * the boundary conditions are read in, and is never called again
80  * during that calculation of the one model.
81  * All default variables that must be initialized before a calculation starts
82  * must appear in the routine. In a grid they are reset for each model
83  */
84 
85  /* parameters having to do with magnetic field */
86  Magnetic_init();
87 
88  /* set all initial abundances */
90 
91  /* zero out parameters needed by large FeII atom */
92  FeIIZero();
93 
94  /* zero out warnings, cautions, notes, etc */
95  wcnint();
96 
97  /* these tell the molecular solver what zone and iteration it have
98  * been evaluated on */
99  hmole_init();
100  H2_Init();
101 
102  /* this is number of iterations that have been malloced - we could
103  * increase this if more iterations are needed */
104  iterations.iter_malloc = 200;
105  /* >>chng 06 jun 27, only malloc on first call - memory leak */
106  if( lgFirstCall)
107  {
108  iterations.IterPrnt = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) );
109  geometry.nend = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) );
110  radius.router = (double*)MALLOC( (size_t)iterations.iter_malloc*sizeof(double) );
111  }
112  for( i=0; i < iterations.iter_malloc; i++ )
113  {
114  iterations.IterPrnt[i] = 10000;
115  }
116  iterations.itermx = 0;
117  /* this implements set coverage command */
118  iterations.lgConverge_set = false;
119  iteration = 0;
120 
121  /* limits for highest and lowest stages of ionization in TrimStage */
122  ionbal.trimhi = 1e-6;
123  ionbal.lgTrimhiOn = true;
124  ionbal.trimlo = 1e-10;
125 
126  hyperfine.lgLya_pump_21cm = true;
127 
128  /* variable to do with geometry */
129  geometry.nprint = 1000;
130  geometry.lgZoneSet = false;
131  geometry.lgZoneTrp = false;
132  geometry.lgEndDflt = true;
133 
134  /* some variables for saving the codes' state */
135  state.lgGet_state = false;
136  state.lgPut_state = false;
137  state.lgState_print = false;
138 
139  /* this is default number of zones
140  * >>chng 96 jun 5, from 400 to 500 for thickest corners4 grid */
141  /* >>chng 04 jan 30, from 600 to 800, code uses finer zoning today */
142  /* >>chng 04 dec 24, from 800 to 1400, so that HII region - molecular cloud
143  * sims do not need set nend - all sims in test suite will run ok without set nend */
144  geometry.nEndDflt = 1400;
145 
146  for( i=0; i < iterations.iter_malloc; i++ )
147  {
149  /*>>chng 03 nov 13, from 1e30 to 1e31, because default inner radius raised to 1e30 */
150  radius.router[i] = 1e31;
151  }
152 
153  /* change angle of illumination
154  * this is angle away from the normal, so 0 is a normal ray, the default*/
157  geometry.fiscal = 1.;
158  geometry.FillFac = 1.;
159  geometry.filpow = 0.;
160 
161  /* default is open geometry, not sphere */
162  geometry.lgSphere = false;
163  /* the radiative transport covering factor */
164  geometry.covrt = 0.;
165  /* the geometric covering factor */
166  geometry.covgeo = 1.;
167  /* default is expanding when geometry set */
168  geometry.lgStatic = false;
169  /* option to tell code not to complain when geometry static done without iterating,
170  * set with (OK) option on geometry command */
171  geometry.lgStaticNoIt = false;
172  /* this is exponent for emissivity contributing to observed luminosity, r^2.
173  * set to 1 with aperture slit, to 0 with aperture beam command */
174  geometry.iEmissPower = 2;
175 
176  carb.p1909 = 0.;
177  carb.p2326 = 0.;
178  oxy.p1666 = 0.;
179  oxy.p1401 = 0.;
180 
181  /* this counts number of times ionize is called by PressureChange, in current zone
182  * these are reset here, so that we count from first zone not search */
183  conv.nPres2Ioniz = 0;
184 
185  /* clear flag indicating the ionization convergence failures
186  * have ever occurred in current zone
187  conv.lgConvIonizThisZone = false; */
188 
189  /* general abort flag */
190  lgAbort = false;
191 
192  /* cooling tolerance heating tolerance - allowed error in heating - cooling balance */
193  /*conv.HeatCoolRelErrorAllowed = 0.02f;*/
194  /* >>chng 04 sep 25, change te tolerance from 0.02 to 4x smaller, 0.005, drove instabilities
195  * in chemistry */
196  conv.HeatCoolRelErrorAllowed = 0.005f;
197 
198  /* this is the default allowed relative error in the electron density */
199  conv.EdenErrorAllowed = 0.01;
200 
201  conv.LimFail = 20;
202  conv.lgMap = false;
203 
204  /* true if level 2 lines were contributors to the ots rates, set in dimacool */
205  conv.lgLevel2_OTS_Imp = true;
206  /* true if level 2 lines were contributors to the cooling, set in dimacool */
207  conv.lgLevel2_Cool_Imp = true;
208 
209  /* this counts how many times ionize is called in this model after startr,
210  * and is flag used by ionize to understand it is being called the first time*/
211  conv.nTotalIoniz = 0;
212  /* these are variables to remember the biggest error in the
213  * electron density, and the zone in which it occurred */
214  conv.BigEdenError = 0.;
215  conv.AverEdenError = 0.;
216  conv.BigHeatCoolError = 0.;
217  conv.AverHeatCoolError = 0.;
218  conv.BigPressError = 0.;
219  conv.AverPressError = 0.;
220  strcpy( conv.chSolverEden , "simple" );
221  /* >>chng 03 mar 24, turning this on resulted in host of errors,
222  * nearly all due to more iterations per zone, test suite took longer time
223  * see edensolver.txt in tsuite directory for full list. Of these, only
224  * orion_hii_pdr, hizqso, and globule, actually had eden faults
225  * make initial steps to solution larger, not dependent on tolerance */
226  /* this did not work - conserv.in failed
227  * this should work - used for dynamics
228  strcpy( conv.chSolverEden , "new" );*/
229 
230  strcpy( conv.chSolverTemp , "simple" );
231  /* most models failed with impossible bracket to solver with this on,
232  * check blister.in for one, this does not really work
233  strcpy( conv.chSolverTemp , "new" );*/
234 
235  /* iterate to convergence flag */
236  conv.lgAutoIt = false;
237  /* convergence criteria */
238  conv.autocv = 0.20f;
239  conv.lgConvTemp = true;
240  conv.lgConvPres = true;
241  conv.lgConvEden = true;
242  /* >>chng 04 jan 25, only set lgConvIoniz true where used in ConvXXX path */
243  /*conv.lgConvIoniz = true;*/
244 
245  /* this option, use the new atmdat_rad_rec recombination rates */
246  t_ADfA::Inst().set_version( PHFIT96 );
247 
248  /* age of the cloud, to check for time-steady */
249  timesc.CloudAgeSet = -1.f;
250  /* some timescale for CO and H2 */
253  /* remains neg if not evaluated */
256 
257  timesc.BigCOMoleForm = 0.;
258 
259  timesc.TimeH21cm = 0.;
261 
262  peimbt.tsqden = 1e7;
263 
264  /* CO related variables */
265  co.codfrc = 0.;
266  co.codtot = 0.;
267  co.CODissHeat = 0.;
268  /* the CO molecule */
269  co.COCoolBigFrac = 0.;
270  co.lgCOCoolCaped = false;
271 
272  /* flag to turn off molecular network */
273  co.lgNoCOMole = false;
274 
275  NumDeriv.lgNumDeriv = false;
276  /* nsum is line pointer within large stack of line intensities */
277  LineSave.nsum = 0;
278 
279  /* index within the line in the line stack
280  * default is Hbeta total - the third line in the stack
281  * 0th is a zero for sanity, 1st is unit, 2nd is a comment */
282  /* >>chng 02 apr 22 from 2 to 3 since added unit at 1 */
283  /* >>chng 06 mar 11, from 3 to -1 will now set to "H 1" 4861 */
284  LineSave.ipNormWavL = -1;
285  LineSave.WavLNorm = 4861.36f;
286  LineSave.lgNormSet = false;
287  LineSave.sig_figs = 4;
288 
289  /* the label for the normalization line */
290  strcpy( LineSave.chNormLab, " " );
291 
292  /* the scale factor for the normalization line */
293  LineSave.ScaleNormLine = 1.;
294 
295  /* this says to work with intrinsic spectrum, set true or 0 for emergent spectrum */
296  LineSave.lgLineEmergent = false;
297 
298  /* this is scale factor, reset with set resolution command, for setting
299  * the continuum resolution. Setting to 0.1 will increase resolution by 10x.
300  * this multiplies the resolution contained in the continuum_mesh.ini file */
302 
304  continuum.lgCon0 = false;
305 
306  /* upper limit to energies of inner shell opacities in Ryd
307  * this is 1 MeV by default */
308  continuum.EnergyKshell = 7.35e4;
309 
310  /* free free heating, cooling, net */
311  CoolHeavy.lgFreeOn = true;
312  CoolHeavy.brems_cool_h = 0.;
313  CoolHeavy.colmet = 0.;
314 
316  hydro.cintot = 0.;
317 
318  /* option to print emissivity instead of intensity/luminosity */
319  hydro.lgHiPop2 = false;
320  hydro.pop2mx = 0.;
321  hydro.lgReevalRecom = true;
322 
323  /* flag for Lya masing */
324  hydro.HCollIonMax = 0.;
325 
326  /* this is flag indicating which type of model atom to use,
327  * set with atom hydrogen populations, departure, lowt */
328  /* >>chng 02 mar 13, this had been set to POPU which had the effect
329  * of forcing the populations solution even for very extreme low ionization.
330  * a few models did have neg pops as a result. There was no reason given for the
331  * change nor was a time stamp present. So, the code had been running without
332  * this protection for some time. Changed back to none here, to enable the
333  * test in hydrolevel, and also made limits there much more stringent (since
334  * they had (unintentionally?) been very stringent for some time now). */
335  /*strcpy( iso.chTypeAtomSet , "POPU" );*/
336  for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
337  {
338  strcpy( iso.chTypeAtomSet[ipISO] , "none" );
339  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
340  {
341  /* what was actually used */
342  strcpy( iso.chTypeAtomUsed[ipISO][nelem] , "none" );
343  }
344  }
345 
346  /* type of hydrogen atom top off, options are " add" and "scal"
347  * in versions 90 this was " add", but was "scal" in 91
348  * >>chng 99 jan 16, changed back to " add"*/
349  /*strcpy( hydro.chHTopType, "scal" );*/
350  strcpy( hydro.chHTopType, " add" );
351 
352  /* Lya excitation temperature, counter for hotter than gas */
353  hydro.TexcLya = 0.;
354  hydro.TLyaMax = 0.;
355  hydro.nLyaHot = 0;
356 
357  /* option to kill damping wings of Lya */
358  hydro.DampOnFac = 1.;
359 
360  /* is continuum pumping of H Lyman lines included? yes, but turned off
361  * with atom h-like Lyman pumping off command */
362  hydro.lgLymanPumping = true;
363 
364  /* multiplicative factor for all continuum pumping of H I Lyman lines,
365  * account for possible emission in the line */
367 
368  /* >>refer abundance D/H Pettini, M., & Bowen, D.V., 2001, ApJ, 560, 41 */
369  /* quoted error is +/- 0.35 */
370  hydro.D2H_ratio = 1.65e-5;
371 
372  /* zero fractions of He0 destruction due to 23S */
373  he.nzone = 0;
374  he.frac_he0dest_23S = 0.;
376 
377  for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
378  {
379  /* flag set by compile he-like command, says to regenerate table of recombination coef */
380  iso.lgCompileRecomb[ipISO] = false;
381  iso.lgNoRecombInterp[ipISO] = false;
382 
383  /* how the gbar cs will be treated - set with atom he-like gbar command */
385  iso.lgCS_Vriens[ipISO] = true;
386  iso.lgCS_Vrinceanu[ipISO] = true;
387 
388  fixit(); /* make this the default for ipH_LIKE if not too slow. */
389  iso.lgCS_Vrinceanu[ipH_LIKE] = false;
390 
391  iso.lgCS_therm_ave[ipISO] = false;
392  iso.lgCS_None[ipISO] = false;
393  /* when set try actually set to 1 or 2, depending on which fit is to be used,
394  * 1 is the broken power law fit */
395  /* >>chng 02 dec 21, change to broken power law fit */
396  iso.nCS_new[ipISO] = 1;
397  /* This flag says whether the density is high enough that helium is sufficiently l-mixed. */
398  iso.lgCritDensLMix[ipISO] = true;
399  /* flag saying whether to include fine-structure mixing in spontaneous decays
400  * set with ATOM HE-LIKE FSM command */
401  iso.lgFSM[ipISO] = 0;
402  /* This is the flag saying whether to generate errors. false means don't. */
403  iso.lgRandErrGen[ipISO] = false;
404  /* this is the flag saying whether we should include excess recombination in the
405  * helike sequence. Should only be off if testing effect of top off approximations. */
406  iso.lgTopoff[ipISO] = true;
407  /* Dielectronic recombination for helike ions is on by default. */
408  iso.lgDielRecom[ipISO] = true;
409 
410  /* number of Lyman lines to include in opacities, this can be vastly larger
411  * than the number of actual levels in the model atom */
412  iso.nLyman[ipISO] = 100;
413  iso.nLyman_malloc[ipISO] = 100;
414 
415  /* controls whether l-mixing and collisional ionization included */
416  iso.lgColl_l_mixing[ipISO] = true;
417  iso.lgColl_excite[ipISO] = true;
418  iso.lgColl_ionize[ipISO] = true;
419  iso.lgPrintNumberOfLevels = false;
420 
421  /* this will become a check on how well case b worked */
422  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
423  {
424  /* option to print departure coefficient */
425  iso.lgPrtDepartCoef[ipISO][nelem] = false;
426  /* print hydrogenic level populations */
427  iso.lgPrtLevelPops[ipISO][nelem] = false;
428  iso.CaseBCheck[ipISO][nelem] = -FLT_MAX;
429  iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.;
430  /* a first guess at the recombination coefficients */
431  iso.RadRec_caseB[ipISO][nelem] = 1e-13;
432  iso.lgLevelsLowered[ipISO][nelem] = false;
433  iso.lgLevelsEverLowered[ipISO][nelem] = false;
434  /* error generation done yet? false means not done. */
435  iso.lgErrGenDone[ipISO][nelem] = false;
436  }
437  }
438 
439  /* Dielectronic recombination forming hydrogen-like ions does not exist. */
440  iso.lgDielRecom[ipH_LIKE] = false;
441 
442  /* smallest transition probability allowed */
443  iso.SmallA = 1e-30f;
444 
445  /* reset with SET IND2 command, turns on/off induced two photon */
446  iso.lgInd2nu_On = false;
447 
448  /* hydrogen redistribution functions */
452 
453  /* this is the upper level for each Lya, which uses the special ipLY_A */
456 
457  /* he-like redistribution functions */
462 
463  /* do not average collision strengths - evaluate at kT
464  * set true with command SET COLLISION STRENGHTS AVERAGE */
465  iso.lgCollStrenThermAver = false;
466 
467 
468  /**********************************************************************
469  * all parameters having to do with secondary ionization
470  * by suprathermal electrons
471  **********************************************************************/
472  secondaries.SetCsupra = 0.;
473  secondaries.lgCSetOn = false;
474  secondaries.lgSecOFF = false;
475  secondaries.SecHIonMax = 0.;
476 
480  secondaries.x12tot = 0.;
481  secondaries.sec2total = 0.;
482 
483  if( lgFirstCall )
484  {
485  /* malloc space for supra[nelem][ion] */
486  secondaries.csupra = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) );
487  secondaries.csupra_effic = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) );
488  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
489  {
490  secondaries.csupra[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) );
491  secondaries.csupra_effic[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) );
492  }
493  }
494  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
495  {
496  for( ion=0; ion<nelem+1; ++ion )
497  {
498  /* secondary ionization rate for each species */
499  secondaries.csupra[nelem][ion] = 0.;
500  /* the rate of each species relative to H0 */
501  secondaries.csupra_effic[nelem][ion] = 1.f;
502  }
503  }
504  /* this scale factor is from table 10 of Tielens & Hollenbach 1985 */
505  secondaries.csupra_effic[ipHELIUM][0] = 1.08f;
506 
507  /* on first call, these arrays do not exist, only zero here on
508  * second and later calls, on first call, create them */
509  if( lgFirstCall )
510  {
511  /* these will save bound electron recoil information data */
513  (long**)MALLOC(sizeof(long*)*(unsigned)LIMELM );
515  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
517  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
519  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
521  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
523  (double****)MALLOC(sizeof(double***)*(unsigned)LIMELM );
525  (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
527  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
529  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
530 
531  /* these are source and sink terms for heavy element ionization balance from the
532  * chemistry */
533  mole.source =
534  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
535  mole.sink =
536  (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
538  (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
539 
540  /* space for ionization recombination arrays */
541  ionbal.RateIonizTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
542  ionbal.RateRecomTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
543  ionbal.RR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
544  ionbal.DR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
545  ionbal.RR_Verner_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
546 
547  /* rate coefficients [cm3 s-1] for Badnell DR recombination */
548  ionbal.DR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
549  ionbal.DR_Badnell_rate_coef_mean_ion = (double *)MALLOC(sizeof(double)*(unsigned)LIMELM );
550  /* do these rate coefficients exist? */
551  ionbal.lgDR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM );
552  ionbal.RR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
553  /* do these rate coefficients exist? */
554  ionbal.lgRR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM );
555  /* rate coefficients [cm3 s-1] for older DR recombination */
556  ionbal.DR_old_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
557 
558  /* create arrays for ions */
559  for(nelem=0; nelem<LIMELM; ++nelem )
560  {
561  ionbal.DR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
562  ionbal.lgDR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) );
563  ionbal.RR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
564  ionbal.lgRR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) );
565  ionbal.DR_old_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
566 
567  ionbal.RateIonizTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
568  ionbal.RateRecomTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
569  ionbal.RR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
570  ionbal.DR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
571  ionbal.RR_Verner_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
572  ionbal.UTA_ionize_rate[nelem] =
573  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
574  ionbal.UTA_heat_rate[nelem] =
575  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
576  ionbal.ipCompRecoil[nelem] =
577  (long*)MALLOC(sizeof(long)*(unsigned)(nelem+1) );
578  ionbal.CompRecoilIonRate[nelem] =
579  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
580  ionbal.CompRecoilIonRateSave[nelem] =
581  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
582  ionbal.CompRecoilHeatRate[nelem] =
583  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
585  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
586  ionbal.PhotoRate_Shell[nelem] =
587  (double***)MALLOC(sizeof(double**)*(unsigned)(nelem+1) );
588  ionbal.CollIonRate_Ground[nelem] =
589  (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
590  /* chemistry source and sink terms for ionization ladders */
591  mole.source[nelem] =
592  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
593  mole.sink[nelem] =
594  (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
595  mole.xMoleChTrRate[nelem] =
596  (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
597  for( ion=0; ion<nelem+2; ++ion )
598  {
599  mole.xMoleChTrRate[nelem][ion] =
600  (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
601  }
602 
603  for( ion=0; ion<nelem+1; ++ion )
604  {
605  /* >>chng 03 aug 09, set these to impossible values */
606  ionbal.RateIonizTot[nelem][ion] = -1.;
607  ionbal.RateRecomTot[nelem][ion] = -1.;
608  ionbal.UTA_ionize_rate[nelem][ion] = -1.;
609  ionbal.UTA_heat_rate[nelem][ion] = -1.;
610  ionbal.ipCompRecoil[nelem][ion] = -99;
611  ionbal.CompRecoilIonRate[nelem][ion] = -1.;
612  ionbal.CompRecoilIonRateSave[nelem][ion] = -1.;
613  ionbal.CompRecoilHeatRate[nelem][ion] = -1.;
614  ionbal.CompRecoilHeatRateSave[nelem][ion] = -1.;
615 
616  /* finish mallocing space */
617  ionbal.PhotoRate_Shell[nelem][ion] =
618  (double**)MALLOC(sizeof(double*)*(unsigned)NSHELLS );
619  ionbal.CollIonRate_Ground[nelem][ion] =
620  (double*)MALLOC(sizeof(double)*(unsigned)2 );
621  for( ns=0; ns<NSHELLS; ++ns )
622  {
623  ionbal.PhotoRate_Shell[nelem][ion][ns] =
624  (double*)MALLOC(sizeof(double)*(unsigned)3 );
625  }
626 
627  /* now set to impossible values */
628  ionbal.ipCompRecoil[nelem][ion] = -100000;
629  ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
630  ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
631  ionbal.DR_old_rate_coef[nelem][ion] = 0.;
632  }
633 
634  set_NaN( ionbal.DR_rate_coef_used[nelem], nelem+1 );
635  set_NaN( ionbal.RR_rate_coef_used[nelem], nelem+1 );
636  set_NaN( ionbal.RR_Verner_rate_coef[nelem], nelem+1 );
637  }
638  }
639 
640  for(ion=0; ion<LIMELM; ++ion )
641  {
643  }
644 
645  /* now zero out these arrays */
646  for( nelem=0; nelem< LIMELM; ++nelem )
647  {
648  for( ion=0; ion<nelem+1; ++ion )
649  {
650 
651  ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
652  ionbal.CompRecoilIonRate[nelem][ion] = 0.;
653  ionbal.UTA_ionize_rate[nelem][ion] = 0.;
654  ionbal.UTA_heat_rate[nelem][ion] = 0.;
655  ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
656  ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
657  ionbal.RateRecomTot[nelem][ion] = 0.;
658  for( ns=0; ns < NSHELLS; ++ns )
659  {
660  /* must be zero since ion routines use these when
661  * not yet defined */
662  ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 0.;
663  ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
664  ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
665  }
666  }
667  /* these have one more ion than above */
668  for( ion=0; ion<nelem+2; ++ion )
669  {
670  long int ion2;
671  /* zero out the source and sink arrays */
672  mole.source[nelem][ion] = 0.;
673  mole.sink[nelem][ion] = 0.;
674  for( ion2=0; ion2<nelem+2; ++ion2 )
675  {
676  mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
677  }
678  }
679  }
680 
681  /* capture of molecules onto grain surfaces - formation of ices
682  * flag says to include this process - turned off with the
683  * command NO GRAIN MOLECULES */
684  mole.lgGrain_mole_deplete = true;
685 
686  ionbal.lgPhotoIoniz_On = true;
687  ionbal.lgCompRecoil = true;
688 
689  /* these three adjust the treatment of UTA ionization */
691  /*>>chng 07 jan 25, turn Kisielius off by default - rates show very
692  * ragged trend with Z, use Gu et al. 06 which are more regular */
694  ionbal.lgInnerShell_Gu06 = true;
695 
696  for( i=0; i<4; ++i )
697  {
698  ionbal.GuessDiel[i] = 1.;
699  }
700 
701  /* default condition is burgess suppressed, Nussbaumer and Storey not */
702  ionbal.lgSupDie[0] = true;
703  ionbal.lgSupDie[1] = false;
704 
705  ionbal.lgNoCota = false;
706  for( i=0; i < LIMELM; i++ )
707  {
708  ionbal.CotaRate[i] = 0.;
709  }
710  ionbal.ilt = 0;
711  ionbal.iltln = 0;
712  ionbal.ilthn = 0;
713  ionbal.ihthn = 0;
714  ionbal.ifail = 0;
715  ionbal.lgGrainIonRecom = true;
716 
717  /*>>chng 06 nov 29, use Badnell DR by default */
718  /* do not turn on Badnell DR by default */
720  /* turn on Badnell DR by default */
722 
723  /* >>chng 06 nov 23, use RR by default */
724  /* do not turn on Badnell RR by default */
726  /* turn on Badnell RR by default */
728 
729  /* option to print recombination coefficient then exit */
731 
732  /* this flag says how to generate S DR -
733  * 0 == use mix of real Badnell DR and guess */
734  ionbal.nDR_S_guess = 0;
735 
736  /* set true to use Badnell mean in place of existing hacks */
739 
740  /* setting this non-zero will turn on guess for dr for entire range of ions */
741  /* >>chng 03 nov 22, change from false to true - turn on process */
742  /* >>refer dr guess Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ in press */
743  ionbal.lg_guess_coef = true;
744  ionbal.guess_noise = 0.;
745 
746  /* do O H charge transfer in chemistry by default, changes with
747  * set HO charge transfer */
748  ionbal.lgHO_ct_chem = true;
749 
750  /**********************************************************************
751  * these are options to print errors to special window,
752  * set with print errors command,
753  * output will go to standard error
754  * defined in cdInit
755  **********************************************************************/
756  lgPrnErr = false;
757  ioPrnErr = stderr;
758 
759  /* the code is ok at startup */
760  broke.lgBroke = false;
761  broke.lgFixit = false;
762  broke.lgCheckit = false;
763 
764  /* main arrays to save ionization fractions*/
765  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
766  {
767  dense.gas_phase[nelem] = 0.;
768  dense.xMolecules[nelem] = 0.;
769  for( ion=0; ion < LIMELM+1; ion++ )
770  {
771  dense.xIonDense[nelem][ion] = 0.;
772  }
773  }
774  dense.xMassTotal = 0.;
775 
776  /* >> chng 05 26, 2006 NPA. Define the mass of each molecule for
777  * use in non-equilibrium chemistry calculations */
778 
779  for(i=0;i<N_H_MOLEC;++i)
780  {
781  co.hmole_mass[i] = 0;
782  }
783 
784  /* Now define the molar mass of each molecule in the hydrogen
785  * chemistry. This also includes He+, which also is a reactant
786  * in the heavy element chemistry */
787  co.hmole_mass[0] = 1.0*ATOMIC_MASS_UNIT;
788  co.hmole_mass[1] = 1.0*ATOMIC_MASS_UNIT;
789  co.hmole_mass[2] = 1.0*ATOMIC_MASS_UNIT;
790  co.hmole_mass[3] = 2.0*ATOMIC_MASS_UNIT;
791  co.hmole_mass[4] = 2.0*ATOMIC_MASS_UNIT;
792  co.hmole_mass[5] = 3.0*ATOMIC_MASS_UNIT;
793  co.hmole_mass[6] = 2.0*ATOMIC_MASS_UNIT;
794  co.hmole_mass[7] = 5.0*ATOMIC_MASS_UNIT;
795  co.hmole_mass[8] = 4.0*ATOMIC_MASS_UNIT;
796 
797  /* this is the simple Fred Hamann FeII atom */
799 
800  /* zero out volume and column density save arrays */
801  MeanZero();
802 
803  /* zero out heating rates */
804  HeatZero();
805 
806  /* some parameters dealing with calcium */
807  ca.Ca2RmLya = 0.;
808  ca.popca2ex = 0.;
809  ca.oldcla = 0.;
810  ca.Ca3d = 0.;
811  ca.Ca4p = 0.;
812  ca.dstCala = 0.;
813 
814  /* C12/C13 isotope ratio, sets the ratio of C12O16 to C13O16 and
815  * C13/C12 1909 */
817 
818  /* flag saying that H2O water destruction rate went to zero */
819  co.lgH2Ozer = false;
820 
821  /* extra electron density, set with eden command */
822  dense.EdenExtra = 0.;
823 
824  /* forced electron density, set with set eden command */
825  dense.EdenSet = 0.;
826 
827  /* this is the default allowed relative error in the pressure */
828  conv.PressureErrorAllowed = 0.01f;
829 
830  /* this is abort option set with SET PRESIONIZ command */
831  conv.limPres2Ioniz = 100000000;
832 
833  /* some titles and line images */
834  for( i=0; i<74; ++i)
835  {
836  input.chTitle[i] = ' ';
837  }
838  input.chTitle[75] = '\0';
839 
840  /* velocity field information */
841  /* the turbulent velocity at illuminated face, internally in cm/s,
842  * but entered with turbulence command in km/s */
843  DoppVel.TurbVel = 0.;
844  /* the parameter F in eq 34 of
845  *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, 624, 773 */
847  /* is TurbVel included in pressure? - can be done two ways, with the velocity
848  * being set of with equipartition - true when TurbVel set if not equipartition,
849  * false with NO PRESSURE option on turbulence command */
850  DoppVel.lgTurb_pressure = true;
851  /* The scale in cm over which the turbulence is dissipated. Normally 0,
852  * only set if dissipate keyword appears on turbulence command */
853  DoppVel.DispScale = 0.;
854  /* equipartition option on turbulence command, to set turbulence from B */
855  DoppVel.lgTurbEquiMag = false;
856 
857  /* pressure related variables */
858  pressure.PresInteg = 0.;
859  pressure.pinzon = 0.;
860 
861  pressure.PresRamCurr = 0.;
863  pressure.lgPradCap = false;
864  pressure.lgPradDen = false;
865  pressure.lgLineRadPresOn = true;
866  /* normally abort when radiation pressure exceeds gas pressure in const pres mod,
867  * this is option to keep going, set with NO ABORT on constant pressure command */
868  pressure.lgRadPresAbortOK = true;
869  /* Ditto for whether to stop at sonic point, this gets set to false
870  * for some of the dynamics pressure modes (strongd, shock, antishock)*/
872  /* this flag will say we hit the sonic point */
873  pressure.lgSonicPoint = false;
874  /* True when no physical solution for desired pressure in strong D fronts */
875  pressure.lgStrongDLimbo = false;
876 
877  pressure.RadBetaMax = 0.;
879  pressure.PresMax = 0.;
880 
881  /* initial and current pressures */
882  pressure.PresTotlInit = 0.;
883  pressure.PresTotlCurr = 0.;
884 
885  /* zero out some dynamics variables */
886  DynaZero();
887 
888  phycon.lgPhysOK = true;
889  /* largest relative changes in Te, ne, H+, H2, and CO in structure
890  * this is computed as part of prtcomment so does not exist when code not talking,
891  * set to zero in zero and still zero if prtcomment not called */
892  phycon.BigJumpTe = 0.;
893  phycon.BigJumpne = 0.;
894  phycon.BigJumpH2 = 0.;
895  phycon.BigJumpCO = 0.;
896 
897  dense.xNucleiTotal = 1.;
898  /* WJH */
899  dense.xMassDensity0 = -1.0f;
900 
901  dense.eden = 1.;
902 
903  /* now set physical conditions array
904  * following will force updating all temperature - density variables */
905  TempChange( 1e4 , true);
906 
907 
908  /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
909  * electron density to account for collisions with atomic H. it is an order of
910  * magnitude guess, so this command provides ability to see whether it affects results */
911  dense.HCorrFac = 1.f;
912 
913  dense.H_sum_in_CO = 0.;
914 
915  atoms.nNegOI = 0;
916  for( i=0; i< N_OI_LEVELS; ++i )
917  {
918  atoms.popoi[i] = 0.;
919  }
920  atoms.popmg2 = 0.;
921 
922  /* do we want to punch negative opacities */
923  opac.lgNegOpacIO = false;
924 
925  /* this is flag for turning on case b */
926  opac.lgCaseB = false;
927 
928  /* this is separate flag for turning off collisions from n=2 */
929  opac.lgCaseB_HummerStorey = false;
930 
931  /* this is separate flag for turning off excited state photoionization */
932  opac.lgCaseB_no_photo = false;
933  /* another case b option, turn off background opacities, no Pdest */
934  opac.lgCaseB_no_pdest = false;
935 
936  /* smallest allowed line and Lya optical depths, reset with
937  * Case B command */
938  opac.taumin = 1e-20f;
939  opac.tlamin = 1e-20f;
940 
941  opac.otsmin = 0.;
942 
943  /* this flag says to use the static opacities,
944  * only evaluate them at start of each new zone.
945  * when set false with
946  * no static opacities
947  * command, always reevaluate them */
948  opac.lgOpacStatic = true;
949 
950  /* set true in radinc if negative opacities ever occur */
951  opac.lgOpacNeg = false;
952 
953  /* can turn of scattering opacities for some geometries */
954  opac.lgScatON = true;
955 
956  /* variables having to do with compiling and/or using the
957  * ancillary file of stored opacities */
958  opac.lgCompileOpac = false;
959  /* "no file opacity" command sets following var false, says not to use file */
961  opac.lgUseFileOpac = false;
962 
963  opac.telec = opac.taumin;
964  opac.thmin = opac.taumin;
965 
966  opac.tpcah[0] = opac.taumin;
967  opac.tpcah[1] = 1e20f;
968 
969  dynamics.dDensityDT = 0.;
970  /* effects of fast neutrons */
971  hextra.frcneu = 0.;
972  hextra.effneu = 1.;
973  hextra.totneu = 0.;
974  hextra.lgNeutrnHeatOn = false;
975  hextra.CrsSecNeutron = 4e-26;
976 
977  opac.stimax[0] = 0.;
978  opac.stimax[1] = 0.;
979 
980  /* this is limit on ionization we shall check for zoning
981  * and prtcomment */
982  struc.dr_ionfrac_limit = 1e-3f;
983  struc.nzone = -1;
984  for(i=0;i<N_H_MOLEC;i++)
985  {
986  /* number of protons in each molecule */
987  int hmi_protons[N_H_MOLEC] = {1,1,1,2,2,3,2,1};
988  /* >>chng 03 nov 28, add electrons to hmi struc */
989  /* number of electrons in each molecule, THAT IS NOT COUNTED AS ION,
990  * used to get H mole electron sum */
991  int hmi_electrons[N_H_MOLEC] = {0,0,-1,0,1,1,0,1};
992 
993  hmi.Hmolec[i] = 0.;
994  hmi.nProton[i] = hmi_protons[i];
995  /* number of free electrons contributed, -1 for H- */
996  hmi.nElectron[i] = hmi_electrons[i];
997  }
998 
999  hmi.H2_total = 0.;
1000  hmi.H2_frac_abund_set = 0.;
1001  hmi.hmihet = 0.;
1002  hmi.h2plus_heat = 0.;
1003  hmi.HeatH2Dish_used = 0.;
1004  hmi.HeatH2Dexc_used = 0.;
1005  hmi.HeatH2Dish_TH85 = 0.;
1006  hmi.HeatH2Dexc_TH85 = 0.;
1007  hmi.HeatH2Dish_BigH2 = 0.;
1008  hmi.HeatH2Dexc_BigH2 = 0.;
1013  hmi.HeatH2DexcMax = 0.;
1014  hmi.CoolH2DexcMax = 0.;
1015  hmi.hmitot = 0.;
1016  hmi.H2Opacity = 0.;
1017  hmi.hmidep = 1.;
1018  hmi.h2dep = 1.;
1019  hmi.h2pdep = 1.;
1020  hmi.h3pdep = 1.;
1021  hmi.BiggestH2 = -1.f;
1022  /* option to turn off H molecules */
1023  hmi.lgNoH2Mole = false;
1024 
1025  /* option to kill effects of H2 in CO chemistry - set with
1026  * set Leiden hack h2* off */
1027  hmi.lgLeiden_Keep_ipMH2s = true;
1028  hmi.lgLeidenCRHack = true;
1029 
1030  /* option to turn on the UMIST rates, naturally this will be 1, set to zero
1031  with the set UMIST rates command */
1032  co.lgUMISTrates = true;
1033 
1034  /* option to use diffuse cloud chemical rates from Table 8 of
1035  * >> refer Federman, S. R. & Zsargo, J. 2003, ApJ, 589, 319
1036  * By default, this is false - changed with set chemistry command */
1037  co.lgFederman = true;
1038  /* option to use effective temperature as defined in
1039  * >> refer Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
1040  * By default, this is false - changed with set chemistry command */
1041  co.lgNonEquilChem = false;
1045  co.lgProtElim = true;
1049  co.lgNeutrals = true;
1050 
1051  /* this says which estimate of the rate of the Solomon process to use,
1052  * default is Tielens & Hollenbach 1985a, changed with
1053  * set h2 Solomon command, options are TH85 and BD96,
1054  * the second for the Bertoldi & Draine rates - they
1055  * differ by 1 dex. when large H2 turned on this is ignored */
1056  /* the Tielens & Hollenbach 1985 treatment */
1057  hmi.chH2_small_model_type = 'T';
1058  /* the improved H2 formalism given by
1059  *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
1060  >>refcon 1990, ApJ, 365, 620 */
1061  hmi.chH2_small_model_type = 'H';
1062  /* the Bertoldi & Draine 1996 treatment */
1063  /* >>chng 03 nov 15, change default to BD96 */
1064  hmi.chH2_small_model_type = 'B';
1065  /* >>chng 05 dec 08, use the Elwert et al. approximations as the default */
1066  hmi.chH2_small_model_type = 'E';
1067 
1068  /* set NaN */
1075 
1084 
1092 
1099 
1106 
1119 
1120  /* default grain formation pumping - Takahashi 2001 */
1121  hmi.chGrainFormPump = 'T';
1122 
1123  /* set which approximation for Jura rate - Cazaux & Tielens
1124  * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
1125  hmi.chJura = 'C';
1126 
1127  /* scale factor to multiply Jura rate, set Jura rate command */
1128  hmi.ScaleJura = 1.f;
1129 
1130  /* binding energy for change in H2 population while on grain surface,
1131  * set with "set h2 Tad" command */
1132  hmi.Tad = 800.;
1133 
1134  /* some vars in the h2 molecule */
1135  H2_Zero();
1136 
1137  /* zero out some column densities */
1138  for( i=0; i < NCOLD; i++ )
1139  {
1140  colden.colden[i] = 0.;
1141  }
1142  colden.He123S = 0.;
1143  colden.coldenH2_ov_vel = 0.;
1144  /* ortho and para column densities */
1145  h2.ortho_colden = 0.;
1146  h2.para_colden = 0.;
1147 
1148  /* F=0 and F=1 column densities of H0*/
1149  colden.H0_21cm_upper =0;
1150  colden.H0_21cm_lower =0;
1151 
1152  for( i=0; i < 5; i++ )
1153  {
1154  colden.C2Pops[i] = 0.;
1155  colden.C2Colden[i] = 0.;
1156  /* pops and column density for SiII atom */
1157  colden.Si2Pops[i] = 0.;
1158  colden.Si2Colden[i] = 0.;
1159  }
1160  for( i=0; i < 3; i++ )
1161  {
1162  /* pops and column density for CI atom */
1163  colden.C1Pops[i] = 0.;
1164  colden.C1Colden[i] = 0.;
1165  /* pops and column density for OI atom */
1166  colden.O1Pops[i] = 0.;
1167  colden.O1Colden[i] = 0.;
1168  /* pops and column density for CIII atom */
1169  colden.C3Pops[i] = 0.;
1170  }
1171  for( i=0; i < 4; i++ )
1172  {
1173  /* pops and column density for CIII atom */
1174  colden.C3Colden[i] = 0.;
1175  }
1176 
1177  /* variables to do with Jeans mass and radius */
1178  colden.TotMassColl = 0.;
1179  colden.tmas = 0.;
1180  colden.wmas = 0.;
1181  colden.rjnmin = FLT_MAX;
1182  colden.ajmmin = FLT_MAX;
1183 
1184  /* >>chng 01 jul 26, moved next statement from below loop to avoid bug in gcc 2.95.3, PvH */
1185  /* default diffuse fields is outward only */
1186  strcpy( rfield.chDffTrns, "OU2" );
1187 
1188  /* three arrays that are needed to save tabulated continua, LIMSPC is
1189  * limit to number of continua that can be entered. This is very large
1190  * for simplicity, the second dimension for the array will only be
1191  * created in parse table when needed, so not waste memory */
1192  if( lgFirstCall )
1193  {
1194  rfield.tNuRyd = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) );
1195  rfield.tslop = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) );
1196  rfield.tFluxLog = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) );
1197  rfield.lgContMalloc = (int*)MALLOC((size_t)(LIMSPC*sizeof(int)) );
1198  /* this flag will say that the above vector has been malloced into the continuum array
1199  * this can be returned once the continuum is generated */
1200  /* >>chng 06 jul 16, move inside if-statement -> this fixes memory leak when optimizing, PvH */
1201  for( i=0; i < LIMSPC; ++i )
1202  {
1203  rfield.lgContMalloc[i] = false;
1204  }
1205  }
1206 
1208 
1209  /* strcpy( rfield.chDffTrns, "OU2" ); */
1210  rfield.lgOutOnly = true;
1211  rfield.lgUSphON = false;
1212 
1213  /* flags for whether continuum is defined over all energies */
1214  rfield.lgMMok = true;
1215  rfield.lgHPhtOK = true;
1216  rfield.lgXRayOK = true;
1217  rfield.lgGamrOK = true;
1218 
1219  /* these are the low and high energy bounds of the continuum */
1220  rfield.emm = 1.001e-8f;
1221  rfield.egamry = 7.354e6f;
1222 
1224 
1225  /* where cloud is optically thick to brems */
1227  rfield.EnergyBremsThin = 0.;
1228 
1229  /* this is the faintest the high-energy tail of the continuum be */
1230  rfield.FluxFaint = 1e-10f;
1231 
1232  for( i=0; i < LIMSPC; i++ )
1233  {
1234  /* this is set true if particular continuum source can vary with time
1235  * set true if TIME appears on intensity / luminosity command line */
1236  rfield.lgTimeVary[i] = false;
1237  /* most continua enter as a beam rather than isotropic */
1238  rfield.lgBeamed[i] = true;
1239  /* default energy range is H-ionizing radiation */
1240  rfield.range[i][0] = HIONPOT;
1241  rfield.range[i][1] = rfield.egamry;
1242  rfield.ioTableRead[i].clear();
1243  }
1244  rfield.comtot = 0.;
1245  rfield.comoff = 1.;
1246  rfield.cmcool = 0.;
1247  rfield.cinrat = 0.;
1252  rfield.EnerGammaRay = 7676.;
1253 
1254  /* variables dealing with the radius */
1255  radius.rinner = 0.;
1256  radius.distance = 0.;
1257  radius.Radius = 0.;
1258  radius.Radius_mid_zone = 0.;
1261  radius.depth_x_fillfac = 0.;
1262  radius.lgRadiusKnown = false;
1263  radius.drad = 0.;
1264  radius.r1r0sq = 1.;
1265  /* this is changed with the roberto command, to go from out to in */
1266  radius.dRadSign = 1.;
1267  radius.lgDrNeg = false;
1268 
1269  /* RDFALT is log of default starting radius (cm) */
1270  /* >>chng 03 nov 12, from 25 to 30 for Lya clouds */
1271  /*radius.rdfalt = 25.;*/
1272  radius.rdfalt = 30.;
1273 
1274  /* set default cylinder thickness */
1275  radius.CylindHigh = 1e35f;
1276  radius.lgCylnOn = false;
1277 
1278  radius.drad_x_fillfac = 1.;
1280  radius.dVeff = 1.;
1281  radius.drNext = 1.;
1282  radius.dRNeff = 1.;
1283  radius.lgdR2Small = false;
1284 
1285  radius.glbdst = 0.;
1286  radius.glbrad = 0.;
1287 
1289  radius.sdrmax = 1e30;
1290  radius.lgSMinON = false;
1291  radius.lgDrMnOn = true;
1292 
1293  radius.lgDrMinUsed = false;
1294 
1295  /* drChange was reset to get orion flux in h-beta correct
1296  * drChange is really tau of current zone */
1297  radius.drChange = 0.15f;
1298 
1299  rfield.lgHabing = false;
1300  /* this flag to make sure more than one table read does not occur in same run */
1301  rfield.lgTableRead = false;
1302 
1303  /* flag to turn off Lya ots */
1304  rfield.lgLyaOTS = true;
1305  /* HeII rec and Lya ots */
1306  rfield.lgHeIIOTS = true;
1307  rfield.lgKillOTSLine = false;
1308  rfield.lgKillOutLine = false;
1309  rfield.lgKillOutCont = false;
1310 
1311  /* rfield.DiffPumpOn is unity unless process disabled by setting to 1
1312  * with no diffuse line pumping command */
1313  rfield.DiffPumpOn = 1.;
1314 
1315  rfield.lgInducProcess = true;
1316  /* 02 jun 13, by Ryan...added this line */
1317  rfield.lgCompileGauntFF = false;
1318 
1319  /* >>chng 03 nov 28, add option to not do line transfer */
1320  rfield.lgDoLineTrans = true;
1321 
1322  /* flag saying whether to constantly reevaluated opacities -
1323  * set false with no opacity reevaluate command */
1324  rfield.lgOpacityReevaluate = true;
1325 
1326  /* flag saying whether to constantly reevaluated ionization -
1327  * set false with no ionization reevaluate command */
1328  rfield.lgIonizReevaluate = true;
1329 
1330  /* this flag says that CMB has been set */
1331  rfield.lgCMB_set = false;
1332 
1333  /* line overlap opacity, turn off with no fine opacity command */
1334  rfield.lgOpacityFine = true;
1335  /* this element is default for choosing line width */
1337  /* there will be this many resolution elements in one FWHM for this element,
1338  * at the lowest temperature to be considered */
1340  /* continuum scale factor for case of time varying continuum */
1342  /* will fine optical depths be punched? */
1343  rfield.lgPunchOpacityFine = false;
1344 
1345  /* first is set true if one of the incident continua needs to have
1346  * H-ionizing radiation blocked. Second is set true is it is blocked
1347  * with extinguish command - want both true if first is true */
1348  rfield.lgMustBlockHIon = false;
1349  rfield.lgBlockHIon = false;
1350 
1351  /* reset some variable related to cooling */
1352  CoolZero();
1353 
1354  thermal.lgCNegChk = true;
1355  thermal.CoolHeatMax = 0.;
1356  thermal.wlCoolHeatMax = 0;
1357  thermal.totcol = 0.;
1358  thermal.heatl = 0.;
1359  thermal.coolheat = 0.;
1360  thermal.lgCExtraOn = false;
1361  thermal.CoolExtra = 0.;
1362  thermal.ctot = 1.;
1363 
1364  thermal.HeatNet = 0.;
1365  thermal.htot = 1.;
1366  thermal.power = 0.;
1367  thermal.FreeFreeTotHeat = 0.;
1368  thermal.char_tran_cool = 0.;
1369  thermal.char_tran_heat = 0.;
1370 
1371  fnzone = 0.;
1372  nzone = 0;
1373  /* save initial condition for talk in case PRINT LAST used */
1375 
1376  oxy.poiii2 = 0.;
1377  oxy.poiii3 = 0.;
1378  oxy.poiexc = 0.;
1379 
1380  oxy.d5007r = 0.;
1381  oxy.d5007t = 0.;
1382  oxy.d4363 = 0.;
1383  oxy.d6300 = 0.;
1384 
1385  atmdat.nsbig = 0;
1386 
1387  /* option to turn off collisional ionization with "no collisional ionization" cmmnd */
1388  atmdat.lgCollIonOn = true;
1389 
1390  /***************************************************
1391  * charge transfer ionization and recombination
1392  ***************************************************/
1393  /* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone
1394  * or cooling due to ct */
1395  atmdat.HCharHeatMax = 0.;
1396  atmdat.HCharCoolMax = 0.;
1397 
1398  atmdat.HIonFrac = 0.;
1399  atmdat.HCharExcIonTotal = 0.;
1400  atmdat.HIonFracMax = 0.;
1401  atmdat.HCharExcRecTotal = 0.;
1402  /* option to turn off all charge transfer, turned off with no charge transfer command */
1403  atmdat.lgCTOn = true;
1404 
1405  /* flag saying that charge transfer heating should be included,
1406  * turned off with no CTHeat commmand */
1407  atmdat.HCharHeatOn = 1.;
1408  for( nelem=0; nelem< LIMELM; ++nelem )
1409  {
1410  for( ion=0; ion<LIMELM; ++ion )
1411  {
1412  atmdat.HCharExcIonOf[nelem][ion] = 0.;
1413  atmdat.HCharExcRecTo[nelem][ion] = 0.;
1414  }
1415  }
1416 
1417  /* >>chng 97 jan 6, from 0 to 8.5e-10*q as per Alex Dalgarno e-mail
1418  * >>chng 97 feb 6, from 8.5e-10*q 1.92e-9x as per Alex Dalgarno e-mail */
1419  atmdat.HCTAlex = 1.92e-9;
1420 
1421  for( nelem=0; nelem < LIMELM; nelem++ )
1422  {
1423  /* these are depletion scale factors */
1424  abund.depset[nelem] = 1.;
1425  /*begin sanity check */
1426  if( abund.depset[nelem] == 0. )
1427  {
1428  fprintf( ioQQQ, " ZERO finds insane abundance or depletion.\n" );
1429  fprintf( ioQQQ, " atomic number=%6ld abundance=%10.2e depletion=%10.2e\n",
1430  nelem, abund.solar[nelem], abund.depset[nelem] );
1431  ShowMe();
1432  cdEXIT(EXIT_FAILURE);
1433  }
1434  /*end sanity check */
1435  }
1436 
1437  /* typical ISM depletion factors, subjective mean of Cowie and Songaila
1438  * and Jenkins
1439  * */
1440  abund.Depletion[0] = 1.;
1441  abund.Depletion[1] = 1.;
1442  abund.Depletion[2] = .16f;
1443  abund.Depletion[3] = .6f;
1444  abund.Depletion[4] = .13f;
1445  abund.Depletion[5] = 0.4f;
1446  abund.Depletion[6] = 1.0f;
1447  abund.Depletion[7] = 0.6f;
1448  abund.Depletion[8] = .3f;
1449  abund.Depletion[9] = 1.f;
1450  abund.Depletion[10] = 0.2f;
1451  abund.Depletion[11] = 0.2f;
1452  abund.Depletion[12] = 0.01f;
1453  abund.Depletion[13] = 0.03f;
1454  abund.Depletion[14] = .25f;
1455  abund.Depletion[15] = 1.0f;
1456  abund.Depletion[16] = 0.4f;
1457  abund.Depletion[17] = 1.0f;
1458  abund.Depletion[18] = .3f;
1459  abund.Depletion[19] = 1e-4f;
1460  abund.Depletion[20] = 5e-3f;
1461  abund.Depletion[21] = 8e-3f;
1462  abund.Depletion[22] = 6e-3f;
1463  abund.Depletion[23] = 6e-3f;
1464  abund.Depletion[24] = 5e-2f;
1465  abund.Depletion[25] = 0.01f;
1466  abund.Depletion[26] = 0.01f;
1467  abund.Depletion[27] = 0.01f;
1468  abund.Depletion[28] = .1f;
1469  abund.Depletion[29] = .25f;
1470 
1471  abund.lgDepln = false;
1472  abund.ScaleMetals = 1.;
1473 
1474  /* this tells the code to use standard Auger yields */
1476 
1477  rt.dTauMase = 0.;
1478  rt.lgMaserCapHit = false;
1479  rt.lgMaserSetDR = false;
1480 
1481  rt.DoubleTau = 1.;
1482  rt.lgFstOn = true;
1483  rt.lgElecScatEscape = true;
1484 
1485  /* there was a call to TestCode */
1486  lgTestCodeCalled = false;
1487  /* test code enabled with set test command */
1488  lgTestCodeEnabled = false;
1489 
1490  /* zero out some grain variables */
1491  GrainZero();
1492 
1493  /* following should never be set non-zero */
1494  FeII.fe21406 = 0.;
1495  FeII.fe21507 = 0.;
1496  FeII.fe21508 = 0.;
1497  FeII.fe21609 = 0.;
1498  FeII.fe21308 = 0.;
1499  FeII.fe21207 = 0.;
1500  FeII.fe21106 = 0.;
1501  FeII.fe21006 = 0.;
1502  FeII.fe21204 = 0.;
1503  FeII.fe21103 = 0.;
1504  FeII.fe21104 = 0.;
1505  FeII.fe21001 = 0.;
1506  FeII.fe21002 = 0.;
1507  FeII.fe20201 = 0.;
1508  FeII.fe20302 = 0.;
1509  FeII.fe20706 = 0.;
1510  FeII.fe20807 = 0.;
1511  FeII.fe20908 = 0.;
1512  FeII.fe21007 = 0.;
1513  FeII.fe21107 = 0.;
1514  FeII.fe21108 = 0.;
1515  FeII.fe21110 = 0.;
1516  FeII.fe21208 = 0.;
1517  FeII.fe21209 = 0.;
1518  FeII.fe21211 = 0.;
1519 
1520  fe.Fe4CoolTot = 0.;
1521  fe.fe40401 = 0.;
1522  fe.fe42836 = 0.;
1523  fe.fe42829 = 0.;
1524  fe.fe42567 = 0.;
1525  fe.fe41207 = 0.;
1526  fe.fe41206 = 0.;
1527  fe.fe41106 = 0.;
1528  fe.fe41007 = 0.;
1529  fe.fe41008 = 0.;
1530  fe.fe40906 = 0.;
1531 
1532  fe.Fe7CoolTot = 0.;
1533 # if 0
1534  fe.Fe7_6087 = 0.;
1535  fe.Fe7_5721 = 0.;
1536  fe.Fe7_6601 = 0.;
1537  fe.Fe7_3760 = 0.;
1538  fe.Fe7_3588 = 0.;
1539 # endif
1540 
1541  /* this is flag saying whether this is very first call,
1542  * a time when space has not been allocated */
1543  lgFirstCall = false;
1544  return;
1545 }
1546 
1547 /*rfield_opac_zero zero out rfield arrays between certain limits */
1549  /* index for first element in arrays to be set to zero */
1550  long lo ,
1551  /* array index for highest element to be set */
1552  long ihi )
1553 {
1554  long int i;
1555 
1556  /* >>chng 01 aug 19, space not allocated yet,
1557  * following code must also be present in contcreatemesh where
1558  * space allocated for the first time */
1559  if( lgRfieldMalloced )
1560  {
1561  unsigned long n=(unsigned long)(ihi-lo+1);
1562  memset(&rfield.OccNumbDiffCont[lo] , 0 , n*sizeof(realnum) );
1563  memset(&rfield.OccNumbContEmitOut[lo] , 0 , n*sizeof(realnum) );
1564  memset(&rfield.ContBoltz[lo] , 0 , n*sizeof(double) );
1565  /*>>chng 06 aug 15, this is now 2D array, saving diffuse continuum
1566  * over all zones for use in exact RT */
1567  /*memset(&rfield.ConEmitLocal[lo] , 0 , n*sizeof(realnum) );*/
1568  memset(&rfield.ConEmitReflec[lo] , 0 , n*sizeof(realnum) );
1569  memset(&rfield.ConEmitOut[lo] , 0 , n*sizeof(realnum) );
1570  memset(&rfield.reflin[lo] , 0 , n*sizeof(realnum) );
1571  memset(&rfield.ConRefIncid[lo] , 0 , n*sizeof(realnum) );
1572  memset(&rfield.SummedCon[lo] , 0 , n*sizeof(realnum) );
1573  memset(&rfield.OccNumbBremsCont[lo] , 0 , n*sizeof(realnum) );
1574  memset(&rfield.convoc[lo] , 0 , n*sizeof(realnum) );
1575  memset(&rfield.flux[lo] , 0 , n*sizeof(realnum) );
1576  memset(&rfield.flux_beam_const_save[lo] , 0 , n*sizeof(realnum) );
1577  memset(&rfield.flux_time_beam_save[lo] , 0 , n*sizeof(realnum) );
1578  memset(&rfield.flux_isotropic_save[lo] , 0 , n*sizeof(realnum) );
1579  memset(&rfield.SummedOcc[lo] , 0 , n*sizeof(realnum) );
1580  memset(&rfield.SummedDif[lo] , 0 , n*sizeof(realnum) );
1581  memset(&rfield.flux_accum[lo] , 0 , n*sizeof(realnum) );
1582  memset(&rfield.otslin[lo] , 0 , n*sizeof(realnum) );
1583  memset(&rfield.otscon[lo] , 0 , n*sizeof(realnum) );
1584  memset(&rfield.ConInterOut[lo] , 0 , n*sizeof(realnum) );
1585  memset(&rfield.outlin[lo] , 0 , n*sizeof(realnum) );
1586  memset(&rfield.outlin_noplot[lo] , 0 , n*sizeof(realnum) );
1587  memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(realnum) );
1588  memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(realnum) );
1589  memset(&rfield.otsconNew[lo] , 0 , n*sizeof(realnum) );
1590  memset(&rfield.otslinNew[lo] , 0 , n*sizeof(realnum) );
1591  memset(&opac.OldOpacSave[lo] , 0 , n*sizeof(double) );
1592  memset(&opac.opacity_abs[lo] , 0 , n*sizeof(double) );
1593  memset(&opac.opacity_sct[lo] , 0 , n*sizeof(double) );
1594  memset(&opac.albedo[lo] , 0 , n*sizeof(double) );
1595  memset(&opac.FreeFreeOpacity[lo] , 0 , n*sizeof(double) );
1596 
1597  /* these are not defined on first iteration */
1598  memset( &opac.E2TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
1599  memset( &opac.E2TauAbsOut[lo] , 0 , n*sizeof(realnum) );
1600  memset( &opac.TauAbsTotal[lo] , 0 , n*sizeof(realnum) );
1601 
1602  for( i=lo; i <= ihi; i++ )
1603  {
1604  opac.TauTotalGeo[0][i] = opac.taumin;
1605  opac.TauAbsGeo[0][i] = opac.taumin;
1606  opac.TauScatGeo[0][i] = opac.taumin;
1607  opac.tmn[i] = 1.;
1608  opac.ExpZone[i] = 1.;
1609  opac.E2TauAbsFace[i] = 1.;
1610  opac.ExpmTau[i] = 1.;
1611  opac.OpacStatic[i] = 1.;
1612  }
1613  /* also zero out fine opacity fine grid fine mesh array */
1614  memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
1615  /* also zero out fine opacity array */
1616  memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) );
1617  }
1618  return;
1619 }
1620 

Generated for cloudy by doxygen 1.8.1.1