cloudy trunk
Loading...
Searching...
No Matches
zero.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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 "co.h"
22#include "coolheavy.h"
23#include "lines.h"
24#include "dynamics.h"
25#include "carb.h"
26#include "mean.h"
27#include "atomfeii.h"
28#include "iso.h"
29#include "conv.h"
30#include "geometry.h"
31#include "timesc.h"
32#include "peimbt.h"
33#include "ionbal.h"
34#include "continuum.h"
35#include "atmdat.h"
36#include "mole.h"
37#include "ca.h"
38#include "input.h"
39#include "atoms.h"
40#include "pressure.h"
41#include "numderiv.h"
42#include "colden.h"
43#include "yield.h"
44#include "hmi.h"
45#include "rfield.h"
46#include "abund.h"
47#include "radius.h"
48#include "opacity.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#include "dark_matter.h"
59
60// //////////////////////////////////////////////////////////////////////////
61//
62//
63// NB DO NOT ADD VARIABLES TO THIS FILE! THE GOAL IS TO REMOVE THIS FILE
64// initialization of variables should be done in one of the ini_*.cpp routines
65//
66//
67// //////////////////////////////////////////////////////////////////////////
68
69/* zero out or initialize variables, called by cdInit, but also by optimize_func
70 * during optimization, called before command parser, one time per model,
71 * in a grid one time per grid point (so called nGrid times),
72 * only one time in multi-iteration models */
73void zero(void)
74{
75 //long int i;
76
77 /* this is used to signify the first call to this routine. At that
78 * stage some memory has not been allocated so must not initialize,
79 * set false at very end of this routine */
80 static bool lgFirstCall = true;
81
82 DEBUG_ENTRY( "zero()" );
83
84 /* this routine is called exactly one time at the start of
85 * the calculation of a single model. When the code is used as a subroutine
86 * this routine is called one time for each model. It is called before
87 * the boundary conditions are read in, and is never called again
88 * during that calculation of the one model.
89 * All default variables that must be initialized before a calculation starts
90 * must appear in the routine. In a grid they are reset for each model
91 */
92
93 /* parameters having to do with magnetic field */
95
96 /* set all initial abundances */
98
99 /* zero out parameters needed by large FeII atom */
100 FeIIZero();
101
102 /* zero out warnings, cautions, notes, etc */
103 wcnint();
104
105 /* this is number of iterations that have been malloced - we could
106 * increase this if more iterations are needed */
107 iterations.iter_malloc = 200;
108 /* >>chng 06 jun 27, only malloc on first call - memory leak */
109 if( lgFirstCall)
110 {
111 iterations.IterPrnt = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) );
112 geometry.nend = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) );
113 radius.StopThickness = (double*)MALLOC( (size_t)iterations.iter_malloc*sizeof(double) );
114 radius.StopRadius = (double*)MALLOC( (size_t)iterations.iter_malloc*sizeof(double) );
115 }
116 for( long i=0; i < iterations.iter_malloc; i++ )
117 {
118 iterations.IterPrnt[i] = 10000;
119 }
120 iterations.itermx = 0;
121 /* this implements set coverage command */
122 iterations.lgConverge_set = false;
123 iteration = 0;
124
125 /* limits for highest and lowest stages of ionization in TrimStage */
126 ionbal.trimhi = 1e-6;
127 ionbal.lgTrimhiOn = true;
128 ionbal.trimlo = 1e-10;
129
130 hyperfine.lgLya_pump_21cm = true;
131
132 /* variable to do with geometry */
133 geometry.nprint = 1000;
134 geometry.lgZoneSet = false;
135 geometry.lgZoneTrp = false;
136 geometry.lgEndDflt = true;
137
138 /* some variables for saving the codes' state */
139 state.lgGet_state = false;
140 state.lgPut_state = false;
141 state.lgState_print = false;
142
143 /* this is default number of zones
144 * >>chng 96 jun 5, from 400 to 500 for thickest corners4 grid */
145 /* >>chng 04 jan 30, from 600 to 800, code uses finer zoning today */
146 /* >>chng 04 dec 24, from 800 to 1400, so that HII region - molecular cloud
147 * sims do not need set nend - all sims in test suite will run ok without set nend */
148 geometry.nEndDflt = 1400;
149
150 for( long i=0; i < iterations.iter_malloc; i++ )
151 {
152 geometry.nend[i] = geometry.nEndDflt;
153 /*>>chng 03 nov 13, from 1e30 to 1e31, because default inner radius raised to 1e30 */
154 radius.StopThickness[i] = 1e31;
155 radius.StopRadius[i] = -1.;
156 }
157
158 geometry.fiscal = 1.;
159 geometry.FillFac = 1.;
160 geometry.filpow = 0.;
161
162 /* default is open geometry, not sphere */
163 geometry.lgSphere = false;
164 /* the radiative transport covering factor */
165 geometry.covrt = 0.;
166 /* the geometric covering factor */
167 geometry.covgeo = 1.;
168 /* default is expanding when geometry set */
169 geometry.lgStatic = false;
170 /* option to tell code not to complain when geometry static done without iterating,
171 * set with (OK) option on geometry command */
172 geometry.lgStaticNoIt = false;
173 /* this is exponent for emissivity contributing to observed luminosity, r^2.
174 * set to 1 with aperture slit, to 0 with aperture beam command */
175 geometry.iEmissPower = 2;
176
177 /* this counts number of times ionize is called by PressureChange, in current zone
178 * these are reset here, so that we count from first zone not search */
179 conv.nPres2Ioniz = 0;
180
181 /* clear flag indicating the ionization convergence failures
182 * have ever occurred in current zone
183 conv.lgConvIonizThisZone = false; */
184
185 conv.resetCounters();
186
187 /* general abort flag */
188 lgAbort = false;
189
190 /* cooling tolerance heating tolerance - allowed error in heating - cooling balance */
191 /*conv.HeatCoolRelErrorAllowed = 0.02f;*/
192 /* >>chng 04 sep 25, change te tolerance from 0.02 to 4x smaller, 0.005, drove instabilities
193 * in chemistry */
194 conv.HeatCoolRelErrorAllowed = 0.005f;
195
196 /* this is the default allowed relative error in the electron density */
197 conv.EdenErrorAllowed = 1e-2;
198
199 conv.IonizErrorAllowed = 1e-2;
200
201 conv.dCmHdT = 0.;
202
203 conv.LimFail = 20;
204 conv.lgMap = false;
205
206 /* this counts how many times ionize is called in this model after startr,
207 * and is flag used by ionize to understand it is being called the first time*/
208 conv.nTotalIoniz = 0;
209 /* these are variables to remember the biggest error in the
210 * electron density, and the zone in which it occurred */
211 conv.BigEdenError = 0.;
212 conv.AverEdenError = 0.;
213 conv.BigHeatCoolError = 0.;
214 conv.AverHeatCoolError = 0.;
215 conv.BigPressError = 0.;
216 conv.AverPressError = 0.;
217 strcpy( conv.chSolverEden, "vWDB" );
218 strcpy( conv.chSolverTemp, "vWDB" );
219 strcpy( conv.chNotConverged, "none" );
220 strcpy( conv.chConvEden, "none" );
221 conv.resetConvIoniz();
222 /* iterate to convergence flag */
223 conv.lgAutoIt = false;
224 /* convergence criteria */
225 conv.autocv = 0.20f;
226 conv.lgConvTemp = true;
227 conv.lgConvPres = true;
228 conv.lgConvEden = true;
229 conv.lgUpdateCouplings = false;
230 /* >>chng 04 jan 25, only set lgConvIoniz true where used in ConvXXX path */
231 /*conv.lgConvIoniz = true;*/
232
233 /* this option, use the new atmdat_rad_rec recombination rates */
235
236 /* age of the cloud, to check for time-steady */
237 timesc.CloudAgeSet = -1.f;
238 /* some timescale for CO and H2 */
239 timesc.time_H2_Dest_longest = 0.;
240 timesc.time_H2_Form_longest = 0.;
241 /* remains neg if not evaluated */
242 timesc.time_H2_Dest_here = -1.;
243 timesc.time_H2_Form_here = 0.;
244
245 timesc.BigCOMoleForm = 0.;
246
247 timesc.TimeH21cm = 0.;
248 timesc.sound_speed_isothermal = 0.;
249
250 peimbt.tsqden = 1e7;
251
252 /* CO related variables */
253 co.codfrc = 0.;
254 co.codtot = 0.;
255 co.CODissHeat = 0.;
256
257 NumDeriv.lgNumDeriv = false;
258
259 /* index within the line in the line stack
260 * default is Hbeta total - the third line in the stack
261 * 0th is a zero for sanity, 1st is unit, 2nd is a comment */
262 /* >>chng 02 apr 22 from 2 to 3 since added unit at 1 */
263 /* >>chng 06 mar 11, from 3 to -1 will now set to "H 1" 4861 */
264 LineSave.ipNormWavL = -1;
265 LineSave.WavLNorm = 4861.36f;
266 LineSave.lgNormSet = false;
267 LineSave.sig_figs = 4;
268
269 /* the label for the normalization line */
270 strcpy( LineSave.chNormLab, " " );
271
272 /* the scale factor for the normalization line */
273 LineSave.ScaleNormLine = 1.;
274
275 /* this is scale factor, reset with set resolution command, for setting
276 * the continuum resolution. Setting to 0.1 will increase resolution by 10x.
277 * this multiplies the resolution contained in the continuum_mesh.ini file */
278 continuum.ResolutionScaleFactor = 1.;
279
280 continuum.lgCoStarInterpolationCaution = false;
281 continuum.lgCon0 = false;
282
283 /* upper limit to energies of inner shell opacities in Ryd
284 * this is 1 MeV by default */
285 continuum.EnergyKshell = 7.35e4;
286
287 /* free free heating, cooling, net */
288 CoolHeavy.lgFreeOn = true;
289 CoolHeavy.brems_cool_h = 0.;
290 CoolHeavy.colmet = 0.;
291
292 CoolHeavy.brems_cool_net = 0.;
293 hydro.cintot = 0.;
294
295 /* option to print emissivity instead of intensity/luminosity */
296 hydro.lgHiPop2 = false;
297 hydro.pop2mx = 0.;
298
299 /* flag for Lya masing */
300 hydro.HCollIonMax = 0.;
301
302 /* type of hydrogen atom top off, options are " add" and "scal"
303 * in versions 90 this was " add", but was "scal" in 91
304 * >>chng 99 jan 16, changed back to " add"*/
305 /*strcpy( hydro.chHTopType, "scal" );*/
306 strcpy( hydro.chHTopType, " add" );
307
308 /* Lya excitation temperature, counter for hotter than gas */
309 hydro.TexcLya = 0.;
310 hydro.TLyaMax = 0.;
311 hydro.nLyaHot = 0;
312
313 /* option to kill damping wings of Lya */
314 hydro.DampOnFac = 1.;
315
316 /* is continuum pumping of H Lyman lines included? yes, but turned off
317 * with atom h-like Lyman pumping off command */
318 hydro.lgLymanPumping = true;
319
320 /* multiplicative factor for all continuum pumping of H I Lyman lines,
321 * account for possible emission in the line */
322 hydro.xLymanPumpingScaleFactor = 1.f;
323
324 /* >>refer abundance D/H Pettini, M., & Bowen, D.V., 2001, ApJ, 560, 41 */
325 /* quoted error is +/- 0.35 */
326 hydro.D2H_ratio = 1.65e-5;
327
328 /* zero fractions of He0 destruction due to 23S */
329 he.nzone = 0;
330 he.frac_he0dest_23S = 0.;
331 he.frac_he0dest_23S_photo = 0.;
332
333 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
334 {
335 /* option to disable continuum lowering */
336 iso_ctrl.lgContinuumLoweringEnabled[ipISO] = true;
337
338 /* flag set by compile he-like command, says to regenerate table of recombination coef */
339 iso_ctrl.lgCompileRecomb[ipISO] = false;
340 iso_ctrl.lgNoRecombInterp[ipISO] = false;
341
342 /* how the gbar cs will be treated - set with atom he-like gbar command */
344 iso_ctrl.lgCS_Vriens[ipISO] = true;
345 iso_ctrl.lgCS_Vrinceanu[ipISO] = true;
346
347 fixit(); /* make this the default for ipH_LIKE if not too slow. */
348 iso_ctrl.lgCS_Vrinceanu[ipH_LIKE] = false;
349
350 iso_ctrl.lgCS_therm_ave[ipISO] = false;
351 iso_ctrl.lgCS_None[ipISO] = false;
352 /* when set try actually set to 1 or 2, depending on which fit is to be used,
353 * 1 is the broken power law fit */
354 /* >>chng 02 dec 21, change to broken power law fit */
355 iso_ctrl.nCS_new[ipISO] = 1;
356 /* This flag says whether the density is high enough that helium is sufficiently l-mixed. */
357 iso_ctrl.lgCritDensLMix[ipISO] = true;
358 /* flag saying whether to include fine-structure mixing in spontaneous decays
359 * set with ATOM HE-LIKE FSM command */
360 iso_ctrl.lgFSM[ipISO] = 0;
361 /* This is the flag saying whether to generate errors. false means don't. */
362 iso_ctrl.lgRandErrGen[ipISO] = false;
363 /* this is the flag saying whether we should include excess recombination in the
364 * helike sequence. Should only be off if testing effect of top off approximations. */
365 iso_ctrl.lgTopoff[ipISO] = true;
366 /* Dielectronic recombination for helike ions is on by default. */
367 iso_ctrl.lgDielRecom[ipISO] = true;
368
369 /* number of Lyman lines to include in opacities, this can be vastly larger
370 * than the number of actual levels in the model atom */
371 iso_ctrl.nLyman[ipISO] = 100;
372 iso_ctrl.nLyman_malloc[ipISO] = 100;
373
374 /* controls whether l-mixing and collisional ionization included */
375 iso_ctrl.lgColl_l_mixing[ipISO] = true;
376 iso_ctrl.lgColl_excite[ipISO] = true;
377 iso_ctrl.lgColl_ionize[ipISO] = true;
378 iso_ctrl.lgLTE_levels[ipISO] = false;
379 iso_ctrl.lgPrintNumberOfLevels = false;
380 }
381
382 /* Dielectronic recombination forming hydrogen-like ions does not exist. */
383 iso_ctrl.lgDielRecom[ipH_LIKE] = false;
384
385 /* smallest transition probability allowed */
386 iso_ctrl.SmallA = 1e-30f;
387
388 /* reset with SET IND2 command, turns on/off induced two photon */
389 iso_ctrl.lgInd2nu_On = false;
390
391 /* hydrogen redistribution functions */
392 iso_ctrl.ipLyaRedist[ipH_LIKE] = ipPRD;
393 iso_ctrl.ipResoRedist[ipH_LIKE] = ipCRD;
394 iso_ctrl.ipSubRedist[ipH_LIKE] = ipCRDW;
395
396 /* this is the upper level for each Lya, which uses the special ipLY_A */
397 iso_ctrl.nLyaLevel[ipH_LIKE] = ipH2p;
398 iso_ctrl.nLyaLevel[ipHE_LIKE] = ipHe2p1P;
399
400 /* he-like redistribution functions */
401 iso_ctrl.ipLyaRedist[ipHE_LIKE] = ipPRD;
402 iso_ctrl.ipResoRedist[ipHE_LIKE] = ipCRD;
403 iso_ctrl.ipSubRedist[ipHE_LIKE] = ipCRDW;
404
405 iso_ctrl.lgPessimisticErrors = false;
406
407 /* do not average collision strengths - evaluate at kT
408 * set true with command SET COLLISION STRENGHTS AVERAGE */
409 iso_ctrl.lgCollStrenThermAver = false;
410
411
412 /**********************************************************************
413 * all parameters having to do with secondary ionization
414 * by suprathermal electrons
415 **********************************************************************/
416 secondaries.SetCsupra = 0.;
417 secondaries.lgCSetOn = false;
418 secondaries.lgSecOFF = false;
419 secondaries.SecHIonMax = 0.;
420
421 secondaries.HeatEfficPrimary = 1.;
422 secondaries.SecIon2PrimaryErg = 0.;
423 secondaries.SecExcitLya2PrimaryErg = 0.;
424 secondaries.x12tot = 0.;
425 secondaries.sec2total = 0.;
426
427 if( lgFirstCall )
428 {
429 /* malloc space for supra[nelem][ion] */
430 secondaries.csupra = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) );
431 secondaries.csupra_effic = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) );
432 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
433 {
434 secondaries.csupra[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) );
435 secondaries.csupra_effic[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) );
436 }
437 }
438 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
439 {
440 for( long ion=0; ion<nelem+1; ++ion )
441 {
442 /* secondary ionization rate for each species */
443 secondaries.csupra[nelem][ion] = 0.;
444 /* the rate of each species relative to H0 */
445 secondaries.csupra_effic[nelem][ion] = 1.f;
446 }
447 }
448 /* this scale factor is from table 10 of Tielens & Hollenbach 1985 */
449 secondaries.csupra_effic[ipHELIUM][0] = 1.08f;
450
451 /* on first call, these arrays do not exist, only zero here on
452 * second and later calls, on first call, create them */
453 if( lgFirstCall )
454 {
455 /* these will save bound electron recoil information data */
456 ionbal.ipCompRecoil =
457 (long**)MALLOC(sizeof(long*)*(unsigned)LIMELM );
458 ionbal.CompRecoilIonRate =
459 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
460 ionbal.CompRecoilIonRateSave =
461 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
462 ionbal.CompRecoilHeatRate =
463 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
464 ionbal.CompRecoilHeatRateSave =
465 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
466 ionbal.PhotoRate_Shell =
467 (double****)MALLOC(sizeof(double***)*(unsigned)LIMELM );
468 ionbal.CollIonRate_Ground =
469 (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM );
470 ionbal.UTA_ionize_rate =
471 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
472 ionbal.UTA_heat_rate =
473 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
474
475 /* these are source and sink terms for heavy element ionization balance from the
476 * chemistry */
477 mole.source =
478 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
479 mole.sink =
480 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
481 mole.xMoleChTrRate =
482 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
483
484 /* space for ionization recombination arrays */
485 ionbal.RateIoniz = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
486 ionbal.RateRecomTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
487 ionbal.RateRecomIso = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
488 ionbal.RR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
489 ionbal.RR_Verner_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
490
491 /* rate coefficients [cm3 s-1] for Badnell DR recombination */
492 ionbal.DR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
493 ionbal.RR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
494 ionbal.CX_recomb_rate_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM );
495
496 /* create arrays for ions */
497 for( long nelem=0; nelem<LIMELM; ++nelem )
498 {
499 ionbal.DR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
500 ionbal.RR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
501 ionbal.CX_recomb_rate_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
502
503 ionbal.RateIoniz[nelem] = (double **)MALLOC(sizeof(double *)*(unsigned)(nelem+1) );
504 ionbal.RateRecomTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
505
506 for( long ion=0; ion<nelem+1; ++ion )
507 {
508 ionbal.RateIoniz[nelem][ion] = (double *)MALLOC(sizeof(double )*(unsigned)(nelem+2) );
509 for( long ion2=0; ion2<nelem+2; ++ion2 )
510 ionbal.RateIoniz[nelem][ion][ion2] = 0.;
511 }
512
513 ionbal.RR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
514 ionbal.RR_Verner_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
515 ionbal.UTA_ionize_rate[nelem] =
516 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
517 ionbal.UTA_heat_rate[nelem] =
518 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
519 ionbal.ipCompRecoil[nelem] =
520 (long*)MALLOC(sizeof(long)*(unsigned)(nelem+1) );
521 ionbal.CompRecoilIonRate[nelem] =
522 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
523 ionbal.CompRecoilIonRateSave[nelem] =
524 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
525 ionbal.CompRecoilHeatRate[nelem] =
526 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
527 ionbal.CompRecoilHeatRateSave[nelem] =
528 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) );
529 ionbal.PhotoRate_Shell[nelem] =
530 (double***)MALLOC(sizeof(double**)*(unsigned)(nelem+1) );
531 ionbal.CollIonRate_Ground[nelem] =
532 (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) );
533 /* chemistry source and sink terms for ionization ladders */
534 mole.source[nelem] =
535 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
536 mole.sink[nelem] =
537 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
538 mole.xMoleChTrRate[nelem] =
539 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
540 for( long ion=0; ion<nelem+2; ++ion )
541 {
542 mole.xMoleChTrRate[nelem][ion] =
543 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
544 }
545 ionbal.RateRecomIso[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(NISO) );
546 for( long ipISO=0; ipISO<NISO; ++ipISO )
547 {
548 ionbal.RateRecomIso[nelem][ipISO] = 0.;
549 }
550
551 for( long ion=0; ion<nelem+1; ++ion )
552 {
553 /* >>chng 03 aug 09, set these to impossible values */
554 ionbal.RateRecomTot[nelem][ion] = -1.;
555 ionbal.UTA_ionize_rate[nelem][ion] = -1.;
556 ionbal.UTA_heat_rate[nelem][ion] = -1.;
557 ionbal.ipCompRecoil[nelem][ion] = -99;
558 ionbal.CompRecoilIonRate[nelem][ion] = -1.;
559 ionbal.CompRecoilIonRateSave[nelem][ion] = -1.;
560 ionbal.CompRecoilHeatRate[nelem][ion] = -1.;
561 ionbal.CompRecoilHeatRateSave[nelem][ion] = -1.;
562
563 /* finish mallocing space */
564 ionbal.PhotoRate_Shell[nelem][ion] =
565 (double**)MALLOC(sizeof(double*)*(unsigned)NSHELLS );
566 ionbal.CollIonRate_Ground[nelem][ion] =
567 (double*)MALLOC(sizeof(double)*(unsigned)2 );
568 for( long ns=0; ns<NSHELLS; ++ns )
569 {
570 ionbal.PhotoRate_Shell[nelem][ion][ns] =
571 (double*)MALLOC(sizeof(double)*(unsigned)3 );
572 }
573
574 /* now set to impossible values */
575 ionbal.ipCompRecoil[nelem][ion] = -100000;
576 ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.;
577 ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.;
578 }
579
580 set_NaN( ionbal.RR_rate_coef_used[nelem], nelem+1 );
581 set_NaN( ionbal.RR_Verner_rate_coef[nelem], nelem+1 );
582 set_NaN( ionbal.CX_recomb_rate_used[nelem], nelem+1 );
583 }
584 }
585
586 /* now zero out these arrays */
587 for( long nelem=0; nelem< LIMELM; ++nelem )
588 {
589 for( long ion=0; ion<nelem+1; ++ion )
590 {
591
592 ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
593 ionbal.CompRecoilIonRate[nelem][ion] = 0.;
594 ionbal.UTA_ionize_rate[nelem][ion] = 0.;
595 ionbal.UTA_heat_rate[nelem][ion] = 0.;
596 ionbal.CollIonRate_Ground[nelem][ion][0] = 0.;
597 ionbal.CollIonRate_Ground[nelem][ion][1] = 0.;
598 ionbal.RateRecomTot[nelem][ion] = 0.;
599 for( long ns=0; ns < NSHELLS; ++ns )
600 {
601 /* must be zero since ion routines use these when
602 * not yet defined */
603 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 0.;
604 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
605 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
606 }
607 }
608 /* these have one more ion than above */
609 for( long ion=0; ion<nelem+2; ++ion )
610 {
611 /* zero out the source and sink arrays */
612 mole.source[nelem][ion] = 0.;
613 mole.sink[nelem][ion] = 0.;
614 for( long ion2=0; ion2<nelem+2; ++ion2 )
615 {
616 mole.xMoleChTrRate[nelem][ion][ion2] = 0.;
617 }
618 }
619 }
620
621 ionbal.lgPhotoIoniz_On = true;
622 ionbal.lgCompRecoil = true;
623
624 /* these three adjust the treatment of UTA ionization */
625 ionbal.lgInnerShellLine_on = true;
626 ionbal.lgInnerShell_Kisielius = true;
627 ionbal.lgInnerShell_Gu06 = true;
628
629 /* default condition is burgess suppressed, Nussbaumer and Storey not */
630 ionbal.lgSupDie[0] = true;
631 ionbal.lgSupDie[1] = false;
632
633 ionbal.lgNoCota = false;
634 for( long nelem = 0; nelem < LIMELM; ++nelem )
635 {
636 ionbal.CotaRate[nelem] = 0.;
637 }
638 ionbal.ilt = 0;
639 ionbal.iltln = 0;
640 ionbal.ilthn = 0;
641 ionbal.ihthn = 0;
642 ionbal.ifail = 0;
643 ionbal.lgGrainIonRecom = true;
644
645 /* option to print recombination coefficient then exit */
646 ionbal.lgRecom_Badnell_print = false;
647 ionbal.guess_noise = 0.;
648
649 /**********************************************************************
650 * these are options to print errors to special window,
651 * set with print errors command,
652 * output will go to standard error
653 * defined in cdInit
654 **********************************************************************/
655 lgPrnErr = false;
656 ioPrnErr = stderr;
657
658 /* main arrays to save ionization fractions*/
659 dense.zero();
660 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
661 {
662 dense.SetGasPhaseDensity( nelem, 0. );
663 for( long ion=0; ion < LIMELM+1; ion++ )
664 {
665 dense.xIonDense[nelem][ion] = 0.;
666 }
667 }
668 dense.xMassTotal = 0.;
669
670 /* this is the simple Fred Hamann FeII atom */
672
673 /* zero out volume and column density save arrays */
674 mean.MeanZero();
675
676 /* zero out heating rates */
677 HeatZero();
678
679 /* some parameters dealing with calcium */
680 ca.Ca2RmLya = 0.;
681 ca.popca2ex = 0.;
682 ca.Ca3d = 0.;
683 ca.Ca4p = 0.;
684 ca.dstCala = 0.;
685
686 /* this is the default allowed relative error in the pressure */
687 conv.PressureErrorAllowed = 0.01f;
688
689 conv.MaxFractionalDensityStepPerIteration = 0.03;
690
691 /* default error in total gas-phase density of each element, including molecules */
692 conv.GasPhaseAbundErrorAllowed = 1e-5f;
693
694 /* this is abort option set with SET PRESIONIZ command */
695 conv.limPres2Ioniz = 3000;
696
697 conv.nTeFail = 0;
698 conv.nTotalFailures = 0;
699 conv.nPreFail = 0;
700 conv.failmx = 0.;
701 conv.nIonFail = 0;
702 conv.nPopFail = 0;
703 conv.nNeFail = 0;
704 conv.nGrainFail = 0;
705 conv.dCmHdT = 0.;
706
707 /* some titles and line images */
708 for( long i=0; i<74; ++i)
709 {
710 input.chTitle[i] = ' ';
711 }
712 input.chTitle[75] = '\0';
713
714 /* velocity field information */
715 /* the turbulent velocity at illuminated face, internally in cm/s,
716 * but entered with turbulence command in km/s */
717 DoppVel.TurbVel = 0.;
718 /* is a turbulent gradient imposed? Default is no. */
719 DoppVel.lgTurbLawOn = false;
720 /* the log of the turbulence gradient power law. Default is zero. */
721 DoppVel.TurbVelLaw = 0.;
722 /* the parameter F in eq 34 of
723 *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, 624, 773 */
724 DoppVel.Heiles_Troland_F = 0.;
725 /* is TurbVel included in pressure? - can be done two ways, with the velocity
726 * being set of with equipartition - true when TurbVel set if not equipartition,
727 * false with NO PRESSURE option on turbulence command */
728 DoppVel.lgTurb_pressure = true;
729 /* The scale in cm over which the turbulence is dissipated. Normally 0,
730 * only set if dissipate keyword appears on turbulence command */
731 DoppVel.DispScale = 0.;
732 /* equipartition option on turbulence command, to set turbulence from B */
733 DoppVel.lgTurbEquiMag = false;
734
735 /* pressure related variables */
736
737 pressure.RhoGravity_dark = 0.;
738 pressure.RhoGravity_self = 0.;
739 pressure.RhoGravity_external = 0.;
740 pressure.RhoGravity = 0.;
741 pressure.IntegRhoGravity = 0.;
742 pressure.gravity_symmetry = -1;
743 pressure.self_mass_factor = 1.;
744
745 pressure.PresRamCurr = 0.;
746 pressure.pres_radiation_lines_curr = 0.;
747 pressure.lgPradCap = false;
748 pressure.lgPradDen = false;
749 pressure.lgLineRadPresOn = true;
750 /* normally abort when radiation pressure exceeds gas pressure in const pres mod,
751 * this is option to keep going, set with NO ABORT on constant pressure command */
752 pressure.lgRadPresAbortOK = true;
753 /* Ditto for whether to stop at sonic point, this gets set to false
754 * for some of the dynamics pressure modes (strongd, shock, antishock)*/
755 pressure.lgSonicPointAbortOK = true;
756 /* this flag will say we hit the sonic point */
757 pressure.lgSonicPoint = false;
758 /* True when no physical solution for desired pressure in strong D fronts */
759 pressure.lgStrongDLimbo = false;
760
761 pressure.RadBetaMax = 0.;
762 pressure.ipPradMax_nzone = -1;
763 pressure.PresMax = 0.;
764
765 /* initial and current pressures */
766 pressure.PresTotlInit = 0.;
767 pressure.PresTotlCurr = 0.;
768
769 /* zero out some dynamics variables */
770 DynaZero();
771
772 phycon.lgPhysOK = true;
773 /* largest relative changes in Te, ne, H+, H2, and CO in structure
774 * this is computed as part of prtcomment so does not exist when code not talking,
775 * set to zero in zero and still zero if prtcomment not called */
776 phycon.BigJumpTe = 0.;
777 phycon.BigJumpne = 0.;
778 phycon.BigJumpH2 = 0.;
779 phycon.BigJumpCO = 0.;
780
781 dense.xNucleiTotal = 1.;
782 /* WJH */
783 dense.xMassDensity0 = -1.0f;
784
785 // needed for TempChange to work but arrays needed for EdenChange to
786 // work are not yet defined
787 dense.eden = 1.;
788
789 /* now set physical conditions array
790 * following will force updating all temperature - density variables */
791 TempChange( 1e4 , true);
792
793 /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the
794 * electron density to account for collisions with atomic H. it is an order of
795 * magnitude guess, so this command provides ability to see whether it affects results */
796 dense.HCorrFac = 1.f;
797
798 dark.lgNFW_Set = false;
799 dark.r_200 = 0.;
800 dark.r_s = 0.;
801
802 atoms.nNegOI = 0;
803 for( long i=0; i< N_OI_LEVELS; ++i )
804 {
805 atoms.popoi[i] = 0.;
806 }
807 atoms.popmg2 = 0.;
808
809 /* do we want to save negative opacities */
810 opac.lgNegOpacIO = false;
811
812 opac.otsmin = 0.;
813
814 /* this flag says to use the static opacities,
815 * only evaluate them at start of each new zone.
816 * when set false with
817 * no static opacities
818 * command, always reevaluate them */
819 opac.lgOpacStatic = true;
820
821 /* set true in radinc if negative opacities ever occur */
822 opac.lgOpacNeg = false;
823
824 /* can turn of scattering opacities for some geometries */
825 opac.lgScatON = true;
826
827 /* variables having to do with compiling and/or using the
828 * ancillary file of stored opacities */
829 opac.lgCompileOpac = false;
830 /* "no file opacity" command sets following var false, says not to use file */
832 opac.lgUseFileOpac = false;
833
834 dynamics.Upstream_density = 0.;
835 /* effects of fast neutrons */
836 hextra.frcneu = 0.;
837 hextra.effneu = 1.;
838 hextra.totneu = 0.;
839 hextra.lgNeutrnHeatOn = false;
840 hextra.CrsSecNeutron = 4e-26;
841
842 opac.stimax[0] = 0.;
843 opac.stimax[1] = 0.;
844
845
846 hmi.H2_total = 0.;
847 hmi.H2_total_f = 0.f;
848 hmi.HD_total = 0.;
849 hmi.H2_frac_abund_set = 0.;
850 hmi.hmihet = 0.;
851 hmi.h2plus_heat = 0.;
852 hmi.HeatH2Dish_used = 0.;
853 hmi.HeatH2Dexc_used = 0.;
854 hmi.HeatH2Dish_TH85 = 0.;
855 hmi.HeatH2Dexc_TH85 = 0.;
856 hmi.UV_Cont_rel2_Draine_DB96_face = 0.;
857 hmi.UV_Cont_rel2_Draine_DB96_depth = 0.;
858 hmi.UV_Cont_rel2_Habing_TH85_face = 0.;
859 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.;
860 hmi.HeatH2DexcMax = 0.;
861 hmi.CoolH2DexcMax = 0.;
862 hmi.hmitot = 0.;
863 hmi.H2Opacity = 0.;
864 hmi.hmidep = 1.;
865 hmi.h2dep = 1.;
866 hmi.h2pdep = 1.;
867 hmi.h3pdep = 1.;
868
869 /* option to kill effects of H2 in CO chemistry - set with
870 * set Leiden hack h2* off */
871 hmi.lgLeiden_Keep_ipMH2s = true;
872 hmi.lgLeidenCRHack = true;
873
874 /* flag to turn off molecular network */
875 mole_global.lgNoMole = false;
876 mole_global.lgNoHeavyMole = false;
877 /* capture of molecules onto grain surfaces - formation of ices
878 * flag says to include this process - turned off with the
879 * command NO GRAIN MOLECULES */
880 mole_global.lgGrain_mole_deplete = true;
881 /* flag saying that H2O water destruction rate went to zero */
882 mole_global.lgH2Ozer = false;
883 /* option to turn on the UMIST rates, naturally this will be 1, set to zero
884 with the set UMIST rates command */
885 mole_global.lgLeidenHack = false;
886 /* option to use diffuse cloud chemical rates from Table 8 of
887 * >> refer Federman, S. R. & Zsargo, J. 2003, ApJ, 589, 319
888 * By default, this is false - changed with set chemistry command */
889 mole_global.lgFederman = true;
890 /* option to use effective temperature as defined in
891 * >> refer Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319
892 * By default, this is false - changed with set chemistry command */
893 mole_global.lgNonEquilChem = false;
897 mole_global.lgProtElim = true;
901 mole_global.lgNeutrals = true;
902 /* option to use H2 continuum dissociation cross sections computed by P.C. Stancil
903 * By default, this is true - changed with "set H2 continuum dissociation xxx" command
904 * options are "Stancil" or "AD69" */
905 mole_global.lgStancil = false;
906 // all isotopes are currently disabled by default
907 mole_global.lgTreatIsotopes.resize( LIMELM );
908 fill( mole_global.lgTreatIsotopes.begin(), mole_global.lgTreatIsotopes.end(), false );
909
910 /* this says which estimate of the rate of the Solomon process to use,
911 * default is Tielens & Hollenbach 1985a, changed with
912 * set h2 Solomon command, options are TH85 and BD96,
913 * the second for the Bertoldi & Draine rates - they
914 * differ by 1 dex. when large H2 turned on this is ignored */
915 /* the Tielens & Hollenbach 1985 treatment */
916 hmi.chH2_small_model_type = 'T';
917 /* the improved H2 formalism given by
918 *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M
919 >>refcon 1990, ApJ, 365, 620 */
920 hmi.chH2_small_model_type = 'H';
921 /* the Bertoldi & Draine 1996 treatment */
922 /* >>chng 03 nov 15, change default to BD96 */
923 hmi.chH2_small_model_type = 'B';
924 /* >>chng 05 dec 08, use the Elwert et al. approximations as the default */
925 hmi.chH2_small_model_type = 'E';
926
927 /* set NaN */
928 set_NaN( hmi.HeatH2Dish_used );
929 set_NaN( hmi.HeatH2Dish_TH85 );
930 set_NaN( hmi.HeatH2Dish_BD96 );
931 set_NaN( hmi.HeatH2Dish_BHT90 );
932 set_NaN( hmi.HeatH2Dish_ELWERT );
933
936 set_NaN( hmi.HeatH2Dexc_used );
937 set_NaN( hmi.HeatH2Dexc_TH85 );
938 set_NaN( hmi.HeatH2Dexc_BD96 );
939 set_NaN( hmi.HeatH2Dexc_BHT90 );
940 set_NaN( hmi.HeatH2Dexc_ELWERT );
941
943 set_NaN( hmi.deriv_HeatH2Dexc_used );
944 set_NaN( hmi.deriv_HeatH2Dexc_TH85 );
945 set_NaN( hmi.deriv_HeatH2Dexc_BD96 );
946 set_NaN( hmi.deriv_HeatH2Dexc_BHT90 );
947 set_NaN( hmi.deriv_HeatH2Dexc_ELWERT );
948
949 set_NaN( hmi.H2_Solomon_dissoc_rate_used_H2g );
950 set_NaN( hmi.H2_Solomon_dissoc_rate_TH85_H2g );
951 set_NaN( hmi.H2_Solomon_dissoc_rate_BHT90_H2g );
952 set_NaN( hmi.H2_Solomon_dissoc_rate_BD96_H2g );
953 set_NaN( hmi.H2_Solomon_dissoc_rate_ELWERT_H2g );
954
955 set_NaN( hmi.H2_Solomon_dissoc_rate_used_H2s );
956 set_NaN( hmi.H2_Solomon_dissoc_rate_TH85_H2s );
957 set_NaN( hmi.H2_Solomon_dissoc_rate_BHT90_H2s );
958 set_NaN( hmi.H2_Solomon_dissoc_rate_BD96_H2s );
959 set_NaN( hmi.H2_Solomon_dissoc_rate_ELWERT_H2s );
960
965 set_NaN( hmi.H2_photodissoc_used_H2g );
966 set_NaN( hmi.H2_photodissoc_used_H2s );
967 set_NaN( hmi.H2_photodissoc_ELWERT_H2g );
968 set_NaN( hmi.H2_photodissoc_ELWERT_H2s );
969 set_NaN( hmi.H2_photodissoc_TH85 );
970 set_NaN( hmi.H2_photodissoc_BHT90 );
971
972 /* default grain formation pumping - Takahashi 2001 */
973 hmi.chGrainFormPump = 'T';
974
975 /* set which approximation for Jura rate - Cazaux & Tielens
976 * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */
977 hmi.chJura = 'C';
978
979 /* scale factor to multiply Jura rate, set Jura rate command */
980 hmi.ScaleJura = 1.f;
981
982 /* binding energy for change in H2 population while on grain surface,
983 * set with "set h2 Tad" command */
984 hmi.Tad = 800.;
985
986 hmi.lgH2_Thermal_BigH2 = true;
987 hmi.lgH2_Chemistry_BigH2 = true;
988
989 /* zero out some column densities */
990 for( long i=0; i < NCOLD; i++ )
991 {
992 colden.colden[i] = 0.;
993 }
994 colden.He123S = 0.;
995 colden.coldenH2_ov_vel = 0.;
996
997 /* F=0 and F=1 column densities of H0*/
998 colden.H0_21cm_upper =0;
999 colden.H0_21cm_lower =0;
1000
1001 for( long i=0; i < 5; i++ )
1002 {
1003 colden.C2Pops[i] = 0.;
1004 colden.C2Colden[i] = 0.;
1005 /* pops and column density for SiII atom */
1006 colden.Si2Pops[i] = 0.;
1007 colden.Si2Colden[i] = 0.;
1008 }
1009 for( long i=0; i < 3; i++ )
1010 {
1011 /* pops and column density for CI atom */
1012 colden.C1Pops[i] = 0.;
1013 colden.C1Colden[i] = 0.;
1014 /* pops and column density for OI atom */
1015 colden.O1Pops[i] = 0.;
1016 colden.O1Colden[i] = 0.;
1017 /* pops and column density for CIII atom */
1018 colden.C3Pops[i] = 0.;
1019 }
1020 for( long i=0; i < 4; i++ )
1021 {
1022 /* pops and column density for CIII atom */
1023 colden.C3Colden[i] = 0.;
1024 }
1025
1026 /* variables to do with Jeans mass and radius */
1027 colden.TotMassColl = 0.;
1028 colden.tmas = 0.;
1029 colden.wmas = 0.;
1030 colden.rjnmin = FLT_MAX;
1031 colden.ajmmin = FLT_MAX;
1032
1033 /* variables dealing with the radius */
1034 radius.rinner = 0.;
1035 radius.distance = 0.;
1036 radius.Radius = 0.;
1037 radius.Radius_mid_zone = 0.;
1038 radius.depth = DEPTH_OFFSET;
1039 radius.depth_mid_zone = DEPTH_OFFSET/2.;
1040 radius.depth_x_fillfac = 0.;
1041 radius.lgRadiusKnown = false;
1042 radius.drad = 0.;
1043 radius.drad_mid_zone = 0.;
1044 radius.r1r0sq = 1.;
1045 /* this is changed with the roberto command, to go from out to in */
1046 radius.dRadSign = 1.;
1047
1048 /* RDFALT is log of default starting radius (cm) */
1049 /* >>chng 03 nov 12, from 25 to 30 for Lya clouds */
1050 /*radius.rdfalt = 25.;*/
1051 radius.rdfalt = 30.;
1052
1053 /* set default cylinder thickness */
1054 radius.CylindHigh = 1e35f;
1055 radius.lgCylnOn = false;
1056
1057 radius.drad_x_fillfac = 1.;
1058 radius.darea_x_fillfac = 1.;
1059 radius.dVeffVol = 1.;
1060 radius.dVeffAper = 1.;
1061 radius.drNext = 1.;
1062 radius.dRNeff = 1.;
1063 radius.lgdR2Small = false;
1064
1065 radius.sdrmin = SMALLFLOAT;
1066 radius.lgSdrminRel = false;
1067 radius.sdrmax = 1e30;
1068 radius.lgSdrmaxRel = false;
1069 radius.lgSMinON = false;
1070 radius.lgDrMnOn = true;
1071 radius.lgFixed = false;
1072 radius.sdrmin_rel_depth = 1e-5;
1073
1074 radius.lgDrMinUsed = false;
1075
1076 rfield.lgHabing = false;
1077
1078 /* flag to turn off Lya ots */
1079 rfield.lgLyaOTS = true;
1080 /* HeII rec and Lya ots */
1081 rfield.lgHeIIOTS = true;
1082 rfield.lgKillOTSLine = false;
1083 rfield.lgKillOutLine = false;
1084 rfield.lgKillOutCont = false;
1085
1086 /* rfield.DiffPumpOn is unity unless process disabled by setting to 1
1087 * with no diffuse line pumping command */
1088 rfield.DiffPumpOn = 1.;
1089
1090 /* 02 jun 13, by Ryan...added this line */
1091 rfield.lgCompileGauntFF = false;
1092
1093 /* >>chng 03 nov 28, add option to not do line transfer */
1094 rfield.lgDoLineTrans = true;
1095
1096 /* flag saying whether to constantly reevaluated opacities -
1097 * set false with no opacity reevaluate command */
1098 rfield.lgOpacityReevaluate = true;
1099
1100 /* flag saying whether to constantly reevaluated ionization -
1101 * set false with no ionization reevaluate command */
1102 rfield.lgIonizReevaluate = true;
1103 /* this element is default for choosing line width */
1104 rfield.fine_opac_nelem = ipIRON;
1105 /* there will be this many resolution elements in one FWHM for this element,
1106 * at the lowest temperature to be considered */
1107 rfield.fine_opac_nresolv = 1;
1108 /* continuum scale factor for case of time varying continuum */
1109 rfield.time_continuum_scale = 1.;
1110 /* will fine optical depths be punched? */
1111 rfield.lgSaveOpacityFine = false;
1112
1113 /* first is set true if one of the incident continua needs to have
1114 * H-ionizing radiation blocked. Second is set true is it is blocked
1115 * with extinguish command - want both true if first is true */
1116 rfield.lgMustBlockHIon = false;
1117 rfield.lgBlockHIon = false;
1118
1119 /* reset some variable related to cooling */
1120 CoolZero();
1121
1122 thermal.lgCNegChk = true;
1123 thermal.CoolHeatMax = 0.;
1124 thermal.wlCoolHeatMax = 0;
1125 thermal.totcol = 0.;
1126 thermal.heatl = 0.;
1127 thermal.coolheat = 0.;
1128 thermal.lgCExtraOn = false;
1129 thermal.CoolExtra = 0.;
1130 thermal.ctot = 1.;
1131
1132 thermal.htot = 1.;
1133 thermal.power = 0.;
1134 thermal.FreeFreeTotHeat = 0.;
1135 thermal.char_tran_cool = 0.;
1136 thermal.char_tran_heat = 0.;
1137
1138 fnzone = 0.;
1139 nzone = 0;
1140 /* save initial condition for talk in case PRINT LAST used */
1141 called.lgTalkSave = called.lgTalk;
1142
1143 oxy.poiii2 = 0.;
1144 oxy.poiii3 = 0.;
1145 oxy.poiexc = 0.;
1146
1147 oxy.d5007r = 0.;
1148 oxy.d5007t = 0.;
1149 oxy.d4363 = 0.;
1150 oxy.d6300 = 0.;
1151
1152 atmdat.nsbig = 0;
1153
1154 /***************************************************
1155 * charge transfer ionization and recombination
1156 ***************************************************/
1157 /* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone
1158 * or cooling due to ct */
1159 atmdat.HCharHeatMax = 0.;
1160 atmdat.HCharCoolMax = 0.;
1161
1162 for ( long nelem=0; nelem < t_atmdat::NCX; ++nelem)
1163 {
1164 atmdat.CharExcIonTotal[nelem] = 0.;
1165 atmdat.CharExcRecTotal[nelem] = 0.;
1166 }
1167 atmdat.HIonFrac = 0.;
1168 atmdat.HIonFracMax = 0.;
1169 /* option to turn off all charge transfer, turned off with no charge transfer command */
1170 atmdat.lgCTOn = true;
1171
1172 /* flag saying that charge transfer heating should be included,
1173 * turned off with no CTHeat commmand */
1174 atmdat.HCharHeatOn = 1.;
1175 for( long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
1176 {
1177 for( long nelem=0; nelem< LIMELM; ++nelem )
1178 {
1179 for( long ion=0; ion<LIMELM; ++ion )
1180 {
1181 atmdat.CharExcIonOf[nelem1][nelem][ion] = 0.;
1182 atmdat.CharExcRecTo[nelem1][nelem][ion] = 0.;
1183 }
1184 }
1185 }
1186
1187 /* >>chng 97 jan 6, from 0 to 8.5e-10*q as per Alex Dalgarno e-mail
1188 * >>chng 97 feb 6, from 8.5e-10*q 1.92e-9x as per Alex Dalgarno e-mail */
1189 atmdat.HCTAlex = 1.92e-9;
1190
1191 for( long nelem=0; nelem < LIMELM; nelem++ )
1192 {
1193 /* these are depletion scale factors */
1194 abund.depset[nelem] = 1.;
1195 /*begin sanity check */
1196 if( abund.depset[nelem] == 0. )
1197 {
1198 fprintf( ioQQQ, " ZERO finds insane abundance or depletion.\n" );
1199 fprintf( ioQQQ, " atomic number=%6ld abundance=%10.2e depletion=%10.2e\n",
1200 nelem, abund.solar[nelem], abund.depset[nelem] );
1201 ShowMe();
1203 }
1204 /*end sanity check */
1205 }
1206
1207 /* typical ISM depletion factors, subjective mean of Cowie and Songaila
1208 * and Jenkins
1209 * */
1210 abund.Depletion[0] = 1.;
1211 abund.Depletion[1] = 1.;
1212 abund.Depletion[2] = .16f;
1213 abund.Depletion[3] = .6f;
1214 abund.Depletion[4] = .13f;
1215 abund.Depletion[5] = 0.4f;
1216 abund.Depletion[6] = 1.0f;
1217 abund.Depletion[7] = 0.6f;
1218 abund.Depletion[8] = .3f;
1219 abund.Depletion[9] = 1.f;
1220 abund.Depletion[10] = 0.2f;
1221 abund.Depletion[11] = 0.2f;
1222 abund.Depletion[12] = 0.01f;
1223 abund.Depletion[13] = 0.03f;
1224 abund.Depletion[14] = .25f;
1225 abund.Depletion[15] = 1.0f;
1226 abund.Depletion[16] = 0.4f;
1227 abund.Depletion[17] = 1.0f;
1228 abund.Depletion[18] = .3f;
1229 abund.Depletion[19] = 1e-4f;
1230 abund.Depletion[20] = 5e-3f;
1231 abund.Depletion[21] = 8e-3f;
1232 abund.Depletion[22] = 6e-3f;
1233 abund.Depletion[23] = 6e-3f;
1234 abund.Depletion[24] = 5e-2f;
1235 abund.Depletion[25] = 0.01f;
1236 abund.Depletion[26] = 0.01f;
1237 abund.Depletion[27] = 0.01f;
1238 abund.Depletion[28] = .1f;
1239 abund.Depletion[29] = .25f;
1240
1241 abund.lgDepln = false;
1242 abund.ScaleMetals = 1.;
1243
1244 /* this tells the code to use standard Auger yields */
1246
1247 rt.dTauMase = 0.;
1248 rt.lgMaserCapHit = false;
1249 rt.lgMaserSetDR = false;
1250
1251 rt.DoubleTau = 1.;
1252 rt.lgFstOn = true;
1253 rt.lgElecScatEscape = true;
1254
1255 /* there was a call to TestCode */
1256 lgTestCodeCalled = false;
1257 /* test code enabled with set test command */
1258 lgTestCodeEnabled = false;
1259
1260 /* zero out some grain variables */
1261 GrainZero();
1262
1263 /* this is flag saying whether this is very first call,
1264 * a time when space has not been allocated */
1265 lgFirstCall = false;
1266 return;
1267}
t_abund abund
Definition abund.cpp:5
void AbundancesZero(void)
t_atmdat atmdat
Definition atmdat.cpp:6
@ PHFIT96
Definition atmdat.h:276
void FeIIZero(void)
t_atoms atoms
Definition atoms.cpp:5
const int N_OI_LEVELS
Definition atoms.h:236
t_ca ca
Definition ca.cpp:5
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioPrnErr
Definition cddefines.cpp:9
bool lgTestCodeCalled
Definition cddefines.cpp:11
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgTestCodeEnabled
Definition cddefines.cpp:12
bool lgPrnErr
Definition cddefines.cpp:13
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
const int ipIRON
Definition cddefines.h:330
#define MALLOC(exp)
Definition cddefines.h:501
const double DEPTH_OFFSET
Definition cddefines.h:272
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
const int ipCRDW
Definition cddefines.h:294
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipCRD
Definition cddefines.h:292
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
const int ipPRD
Definition cddefines.h:290
static t_ADfA & Inst()
Definition cddefines.h:175
void set_version(phfit_version val)
Definition atmdat.h:318
void zero_opacity()
void reset_yield()
Definition yield.h:79
t_co co
Definition co.cpp:5
t_colden colden
Definition colden.cpp:5
#define NCOLD
Definition colden.h:9
STATIC void fill(double fenlo, double fenhi, double resolv, long int *n0, long int *ipnt, bool lgCount)
t_continuum continuum
Definition continuum.cpp:5
t_conv conv
Definition conv.cpp:5
void CoolZero(void)
Definition cool_etc.cpp:50
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
void HeatZero(void)
Definition heat_sum.cpp:928
void set_NaN(sys_float &x)
Definition cpu.cpp:682
const realnum SMALLFLOAT
Definition cpu.h:191
t_dark_matter dark
t_dense dense
Definition dense.cpp:24
t_DoppVel DoppVel
Definition doppvel.cpp:5
t_dynamics dynamics
Definition dynamics.cpp:44
void DynaZero(void)
t_geometry geometry
Definition geometry.cpp:5
void GrainZero(void)
Definition grains.cpp:500
t_he he
Definition he.cpp:5
t_hextra hextra
Definition hextra.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_input input
Definition input.cpp:12
t_ionbal ionbal
Definition ionbal.cpp:5
#define NSHELLS
Definition ionbal.h:75
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe2p1P
Definition iso.h:49
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_iterations iterations
Definition iterations.cpp:5
t_LineSave LineSave
Definition lines.cpp:5
void Magnetic_init(void)
Definition magnetic.cpp:133
t_mean mean
Definition mean.cpp:17
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
t_NumDeriv NumDeriv
Definition numderiv.cpp:5
t_opac opac
Definition opacity.cpp:5
t_oxy oxy
Definition oxy.cpp:5
t_peimbt peimbt
Definition peimbt.cpp:5
t_phycon phycon
Definition phycon.cpp:6
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
t_secondaries secondaries
t_state state
Definition state.cpp:19
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
void wcnint(void)
Definition warnings.cpp:13
void zero(void)
Definition zero.cpp:73