cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization, 00004 * this is called before any commands are parsed, called one time per model, at very start */ 00005 /*rfield_optac_zero zero out rfield arrays between certain limits */ 00006 #include "cddefines.h" 00007 #include "physconst.h" 00008 #include "iterations.h" 00009 #include "hydrogenic.h" 00010 #include "oxy.h" 00011 #include "doppvel.h" 00012 #include "dense.h" 00013 #include "hextra.h" 00014 #include "grains.h" 00015 #include "magnetic.h" 00016 #include "state.h" 00017 #include "rt.h" 00018 #include "he.h" 00019 #include "struc.h" 00020 #include "h2.h" 00021 #include "coolheavy.h" 00022 #include "lines.h" 00023 #include "dynamics.h" 00024 #include "carb.h" 00025 #include "mean.h" 00026 #include "atomfeii.h" 00027 #include "iso.h" 00028 #include "conv.h" 00029 #include "geometry.h" 00030 #include "timesc.h" 00031 #include "peimbt.h" 00032 #include "ionbal.h" 00033 #include "continuum.h" 00034 #include "atmdat.h" 00035 #include "mole.h" 00036 #include "ca.h" 00037 #include "input.h" 00038 #include "atoms.h" 00039 #include "pressure.h" 00040 #include "numderiv.h" 00041 #include "colden.h" 00042 #include "yield.h" 00043 #include "hmi.h" 00044 #include "rfield.h" 00045 #include "abund.h" 00046 #include "radius.h" 00047 #include "opacity.h" 00048 #include "broke.h" 00049 #include "secondaries.h" 00050 #include "called.h" 00051 #include "phycon.h" 00052 #include "warnings.h" 00053 #include "thermal.h" 00054 #include "cooling.h" 00055 #include "fe.h" 00056 #include "hyperfine.h" 00057 #include "init.h" 00058 00059 /*zero zero out or initialize variables, called by cdInit, but also by optimize_func during optimization, 00060 * this is called before any commands are parsed, called one time per model, at very start */ 00061 void zero(void) 00062 { 00063 long int i, 00064 ion, 00065 ipISO , 00066 nelem, 00067 ns; 00068 00069 /* this is used to signify the first call to this routine. At that 00070 * stage some memory has not been allocated so must not initialize, 00071 * set false at very end of this routine */ 00072 static bool lgFirstCall = true; 00073 00074 DEBUG_ENTRY( "zero()" ); 00075 00076 /* this routine is called exactly one time at the start of 00077 * the calculation of a single model. When the code is used as a subroutine 00078 * this routine is called one time for each model. It is called before 00079 * the boundary conditions are read in, and is never called again 00080 * during that calculation of the one model. 00081 * All default variables that must be initialized before a calculation starts 00082 * must appear in the routine. In a grid they are reset for each model 00083 */ 00084 00085 /* parameters having to do with magnetic field */ 00086 Magnetic_init(); 00087 00088 /* set all initial abundances */ 00089 AbundancesZero(); 00090 00091 /* zero out parameters needed by large FeII atom */ 00092 FeIIZero(); 00093 00094 /* zero out warnings, cautions, notes, etc */ 00095 wcnint(); 00096 00097 /* these tell the molecular solver what zone and iteration it have 00098 * been evaluated on */ 00099 hmole_init(); 00100 H2_Init(); 00101 00102 /* this is number of iterations that have been malloced - we could 00103 * increase this if more iterations are needed */ 00104 iterations.iter_malloc = 200; 00105 /* >>chng 06 jun 27, only malloc on first call - memory leak */ 00106 if( lgFirstCall) 00107 { 00108 iterations.IterPrnt = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) ); 00109 geometry.nend = (long int*)MALLOC( (size_t)iterations.iter_malloc*sizeof(long int) ); 00110 radius.router = (double*)MALLOC( (size_t)iterations.iter_malloc*sizeof(double) ); 00111 } 00112 for( i=0; i < iterations.iter_malloc; i++ ) 00113 { 00114 iterations.IterPrnt[i] = 10000; 00115 } 00116 iterations.itermx = 0; 00117 /* this implements set coverage command */ 00118 iterations.lgConverge_set = false; 00119 iteration = 0; 00120 00121 /* limits for highest and lowest stages of ionization in TrimStage */ 00122 ionbal.trimhi = 1e-6; 00123 ionbal.lgTrimhiOn = true; 00124 ionbal.trimlo = 1e-10; 00125 00126 hyperfine.lgLya_pump_21cm = true; 00127 00128 /* variable to do with geometry */ 00129 geometry.nprint = 1000; 00130 geometry.lgZoneSet = false; 00131 geometry.lgZoneTrp = false; 00132 geometry.lgEndDflt = true; 00133 00134 /* some variables for saving the codes' state */ 00135 state.lgGet_state = false; 00136 state.lgPut_state = false; 00137 state.lgState_print = false; 00138 00139 /* this is default number of zones 00140 * >>chng 96 jun 5, from 400 to 500 for thickest corners4 grid */ 00141 /* >>chng 04 jan 30, from 600 to 800, code uses finer zoning today */ 00142 /* >>chng 04 dec 24, from 800 to 1400, so that HII region - molecular cloud 00143 * sims do not need set nend - all sims in test suite will run ok without set nend */ 00144 geometry.nEndDflt = 1400; 00145 00146 for( i=0; i < iterations.iter_malloc; i++ ) 00147 { 00148 geometry.nend[i] = geometry.nEndDflt; 00149 /*>>chng 03 nov 13, from 1e30 to 1e31, because default inner radius raised to 1e30 */ 00150 radius.router[i] = 1e31; 00151 } 00152 00153 /* change angle of illumination 00154 * this is angle away from the normal, so 0 is a normal ray, the default*/ 00155 geometry.AngleIllumRadian = 0.; 00156 geometry.DirectionalCosin = 1.; 00157 geometry.fiscal = 1.; 00158 geometry.FillFac = 1.; 00159 geometry.filpow = 0.; 00160 00161 /* default is open geometry, not sphere */ 00162 geometry.lgSphere = false; 00163 /* the radiative transport covering factor */ 00164 geometry.covrt = 0.; 00165 /* the geometric covering factor */ 00166 geometry.covgeo = 1.; 00167 /* default is expanding when geometry set */ 00168 geometry.lgStatic = false; 00169 /* option to tell code not to complain when geometry static done without iterating, 00170 * set with (OK) option on geometry command */ 00171 geometry.lgStaticNoIt = false; 00172 /* this is exponent for emissivity contributing to observed luminosity, r^2. 00173 * set to 1 with aperture slit, to 0 with aperture beam command */ 00174 geometry.iEmissPower = 2; 00175 00176 carb.p1909 = 0.; 00177 carb.p2326 = 0.; 00178 oxy.p1666 = 0.; 00179 oxy.p1401 = 0.; 00180 00181 /* this counts number of times ionize is called by PressureChange, in current zone 00182 * these are reset here, so that we count from first zone not search */ 00183 conv.nPres2Ioniz = 0; 00184 00185 /* clear flag indicating the ionization convergence failures 00186 * have ever occurred in current zone 00187 conv.lgConvIonizThisZone = false; */ 00188 00189 /* general abort flag */ 00190 lgAbort = false; 00191 00192 /* cooling tolerance heating tolerance - allowed error in heating - cooling balance */ 00193 /*conv.HeatCoolRelErrorAllowed = 0.02f;*/ 00194 /* >>chng 04 sep 25, change te tolerance from 0.02 to 4x smaller, 0.005, drove instabilities 00195 * in chemistry */ 00196 conv.HeatCoolRelErrorAllowed = 0.005f; 00197 00198 /* this is the default allowed relative error in the electron density */ 00199 conv.EdenErrorAllowed = 0.01; 00200 00201 conv.LimFail = 20; 00202 conv.lgMap = false; 00203 00204 /* true if level 2 lines were contributors to the ots rates, set in dimacool */ 00205 conv.lgLevel2_OTS_Imp = true; 00206 /* true if level 2 lines were contributors to the cooling, set in dimacool */ 00207 conv.lgLevel2_Cool_Imp = true; 00208 00209 /* this counts how many times ionize is called in this model after startr, 00210 * and is flag used by ionize to understand it is being called the first time*/ 00211 conv.nTotalIoniz = 0; 00212 /* these are variables to remember the biggest error in the 00213 * electron density, and the zone in which it occurred */ 00214 conv.BigEdenError = 0.; 00215 conv.AverEdenError = 0.; 00216 conv.BigHeatCoolError = 0.; 00217 conv.AverHeatCoolError = 0.; 00218 conv.BigPressError = 0.; 00219 conv.AverPressError = 0.; 00220 strcpy( conv.chSolverEden , "simple" ); 00221 /* >>chng 03 mar 24, turning this on resulted in host of errors, 00222 * nearly all due to more iterations per zone, test suite took longer time 00223 * see edensolver.txt in tsuite directory for full list. Of these, only 00224 * orion_hii_pdr, hizqso, and globule, actually had eden faults 00225 * make initial steps to solution larger, not dependent on tolerance */ 00226 /* this did not work - conserv.in failed 00227 * this should work - used for dynamics 00228 strcpy( conv.chSolverEden , "new" );*/ 00229 00230 strcpy( conv.chSolverTemp , "simple" ); 00231 /* most models failed with impossible bracket to solver with this on, 00232 * check blister.in for one, this does not really work 00233 strcpy( conv.chSolverTemp , "new" );*/ 00234 00235 /* iterate to convergence flag */ 00236 conv.lgAutoIt = false; 00237 /* convergence criteria */ 00238 conv.autocv = 0.20f; 00239 conv.lgConvTemp = true; 00240 conv.lgConvPres = true; 00241 conv.lgConvEden = true; 00242 /* >>chng 04 jan 25, only set lgConvIoniz true where used in ConvXXX path */ 00243 /*conv.lgConvIoniz = true;*/ 00244 00245 /* this option, use the new atmdat_rad_rec recombination rates */ 00246 t_ADfA::Inst().set_version( PHFIT96 ); 00247 00248 /* age of the cloud, to check for time-steady */ 00249 timesc.CloudAgeSet = -1.f; 00250 /* some timescale for CO and H2 */ 00251 timesc.time_H2_Dest_longest = 0.; 00252 timesc.time_H2_Form_longest = 0.; 00253 /* remains neg if not evaluated */ 00254 timesc.time_H2_Dest_here = -1.; 00255 timesc.time_H2_Form_here = 0.; 00256 00257 timesc.BigCOMoleForm = 0.; 00258 00259 timesc.TimeH21cm = 0.; 00260 timesc.sound_speed_isothermal = 0.; 00261 00262 peimbt.tsqden = 1e7; 00263 00264 /* CO related variables */ 00265 co.codfrc = 0.; 00266 co.codtot = 0.; 00267 co.CODissHeat = 0.; 00268 /* the CO molecule */ 00269 co.COCoolBigFrac = 0.; 00270 co.lgCOCoolCaped = false; 00271 00272 /* flag to turn off molecular network */ 00273 co.lgNoCOMole = false; 00274 00275 NumDeriv.lgNumDeriv = false; 00276 /* nsum is line pointer within large stack of line intensities */ 00277 LineSave.nsum = 0; 00278 00279 /* index within the line in the line stack 00280 * default is Hbeta total - the third line in the stack 00281 * 0th is a zero for sanity, 1st is unit, 2nd is a comment */ 00282 /* >>chng 02 apr 22 from 2 to 3 since added unit at 1 */ 00283 /* >>chng 06 mar 11, from 3 to -1 will now set to "H 1" 4861 */ 00284 LineSave.ipNormWavL = -1; 00285 LineSave.WavLNorm = 4861.36f; 00286 LineSave.lgNormSet = false; 00287 LineSave.sig_figs = 4; 00288 00289 /* the label for the normalization line */ 00290 strcpy( LineSave.chNormLab, " " ); 00291 00292 /* the scale factor for the normalization line */ 00293 LineSave.ScaleNormLine = 1.; 00294 00295 /* this says to work with intrinsic spectrum, set true or 0 for emergent spectrum */ 00296 LineSave.lgLineEmergent = false; 00297 00298 /* this is scale factor, reset with set resolution command, for setting 00299 * the continuum resolution. Setting to 0.1 will increase resolution by 10x. 00300 * this multiplies the resolution contained in the continuum_mesh.ini file */ 00301 continuum.ResolutionScaleFactor = 1.; 00302 00303 continuum.lgCoStarInterpolationCaution = false; 00304 continuum.lgCon0 = false; 00305 00306 /* upper limit to energies of inner shell opacities in Ryd 00307 * this is 1 MeV by default */ 00308 continuum.EnergyKshell = 7.35e4; 00309 00310 /* free free heating, cooling, net */ 00311 CoolHeavy.lgFreeOn = true; 00312 CoolHeavy.brems_cool_h = 0.; 00313 CoolHeavy.colmet = 0.; 00314 00315 CoolHeavy.brems_cool_net = 0.; 00316 hydro.cintot = 0.; 00317 00318 /* option to print emissivity instead of intensity/luminosity */ 00319 hydro.lgHiPop2 = false; 00320 hydro.pop2mx = 0.; 00321 hydro.lgReevalRecom = true; 00322 00323 /* flag for Lya masing */ 00324 hydro.HCollIonMax = 0.; 00325 00326 /* this is flag indicating which type of model atom to use, 00327 * set with atom hydrogen populations, departure, lowt */ 00328 /* >>chng 02 mar 13, this had been set to POPU which had the effect 00329 * of forcing the populations solution even for very extreme low ionization. 00330 * a few models did have neg pops as a result. There was no reason given for the 00331 * change nor was a time stamp present. So, the code had been running without 00332 * this protection for some time. Changed back to none here, to enable the 00333 * test in hydrolevel, and also made limits there much more stringent (since 00334 * they had (unintentionally?) been very stringent for some time now). */ 00335 /*strcpy( iso.chTypeAtomSet , "POPU" );*/ 00336 for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00337 { 00338 strcpy( iso.chTypeAtomSet[ipISO] , "none" ); 00339 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00340 { 00341 /* what was actually used */ 00342 strcpy( iso.chTypeAtomUsed[ipISO][nelem] , "none" ); 00343 } 00344 } 00345 00346 /* type of hydrogen atom top off, options are " add" and "scal" 00347 * in versions 90 this was " add", but was "scal" in 91 00348 * >>chng 99 jan 16, changed back to " add"*/ 00349 /*strcpy( hydro.chHTopType, "scal" );*/ 00350 strcpy( hydro.chHTopType, " add" ); 00351 00352 /* Lya excitation temperature, counter for hotter than gas */ 00353 hydro.TexcLya = 0.; 00354 hydro.TLyaMax = 0.; 00355 hydro.nLyaHot = 0; 00356 00357 /* option to kill damping wings of Lya */ 00358 hydro.DampOnFac = 1.; 00359 00360 /* is continuum pumping of H Lyman lines included? yes, but turned off 00361 * with atom h-like Lyman pumping off command */ 00362 hydro.lgLymanPumping = true; 00363 00364 /* multiplicative factor for all continuum pumping of H I Lyman lines, 00365 * account for possible emission in the line */ 00366 hydro.xLymanPumpingScaleFactor = 1.f; 00367 00368 /* >>refer abundance D/H Pettini, M., & Bowen, D.V., 2001, ApJ, 560, 41 */ 00369 /* quoted error is +/- 0.35 */ 00370 hydro.D2H_ratio = 1.65e-5; 00371 00372 /* zero fractions of He0 destruction due to 23S */ 00373 he.nzone = 0; 00374 he.frac_he0dest_23S = 0.; 00375 he.frac_he0dest_23S_photo = 0.; 00376 00377 for( ipISO=ipH_LIKE; ipISO<NISO; ipISO++ ) 00378 { 00379 /* flag set by compile he-like command, says to regenerate table of recombination coef */ 00380 iso.lgCompileRecomb[ipISO] = false; 00381 iso.lgNoRecombInterp[ipISO] = false; 00382 00383 /* how the gbar cs will be treated - set with atom he-like gbar command */ 00385 iso.lgCS_Vriens[ipISO] = true; 00386 iso.lgCS_Vrinceanu[ipISO] = true; 00387 00388 fixit(); /* make this the default for ipH_LIKE if not too slow. */ 00389 iso.lgCS_Vrinceanu[ipH_LIKE] = false; 00390 00391 iso.lgCS_therm_ave[ipISO] = false; 00392 iso.lgCS_None[ipISO] = false; 00393 /* when set try actually set to 1 or 2, depending on which fit is to be used, 00394 * 1 is the broken power law fit */ 00395 /* >>chng 02 dec 21, change to broken power law fit */ 00396 iso.nCS_new[ipISO] = 1; 00397 /* This flag says whether the density is high enough that helium is sufficiently l-mixed. */ 00398 iso.lgCritDensLMix[ipISO] = true; 00399 /* flag saying whether to include fine-structure mixing in spontaneous decays 00400 * set with ATOM HE-LIKE FSM command */ 00401 iso.lgFSM[ipISO] = 0; 00402 /* This is the flag saying whether to generate errors. false means don't. */ 00403 iso.lgRandErrGen[ipISO] = false; 00404 /* this is the flag saying whether we should include excess recombination in the 00405 * helike sequence. Should only be off if testing effect of top off approximations. */ 00406 iso.lgTopoff[ipISO] = true; 00407 /* Dielectronic recombination for helike ions is on by default. */ 00408 iso.lgDielRecom[ipISO] = true; 00409 00410 /* number of Lyman lines to include in opacities, this can be vastly larger 00411 * than the number of actual levels in the model atom */ 00412 iso.nLyman[ipISO] = 100; 00413 iso.nLyman_malloc[ipISO] = 100; 00414 00415 /* controls whether l-mixing and collisional ionization included */ 00416 iso.lgColl_l_mixing[ipISO] = true; 00417 iso.lgColl_excite[ipISO] = true; 00418 iso.lgColl_ionize[ipISO] = true; 00419 iso.lgPrintNumberOfLevels = false; 00420 00421 /* this will become a check on how well case b worked */ 00422 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00423 { 00424 /* option to print departure coefficient */ 00425 iso.lgPrtDepartCoef[ipISO][nelem] = false; 00426 /* print hydrogenic level populations */ 00427 iso.lgPrtLevelPops[ipISO][nelem] = false; 00428 iso.CaseBCheck[ipISO][nelem] = -FLT_MAX; 00429 iso.TwoNu_induc_dn_max[ipISO][nelem] = 0.; 00430 /* a first guess at the recombination coefficients */ 00431 iso.RadRec_caseB[ipISO][nelem] = 1e-13; 00432 iso.lgLevelsLowered[ipISO][nelem] = false; 00433 iso.lgLevelsEverLowered[ipISO][nelem] = false; 00434 /* error generation done yet? false means not done. */ 00435 iso.lgErrGenDone[ipISO][nelem] = false; 00436 } 00437 } 00438 00439 /* Dielectronic recombination forming hydrogen-like ions does not exist. */ 00440 iso.lgDielRecom[ipH_LIKE] = false; 00441 00442 /* smallest transition probability allowed */ 00443 iso.SmallA = 1e-30f; 00444 00445 /* reset with SET IND2 command, turns on/off induced two photon */ 00446 iso.lgInd2nu_On = false; 00447 00448 /* hydrogen redistribution functions */ 00449 iso.ipLyaRedist[ipH_LIKE] = ipLY_A; 00450 iso.ipResoRedist[ipH_LIKE] = ipCRD; 00451 iso.ipSubRedist[ipH_LIKE] = ipCRDW; 00452 00453 /* this is the upper level for each Lya, which uses the special ipLY_A */ 00454 iso.nLyaLevel[ipH_LIKE] = ipH2p; 00455 iso.nLyaLevel[ipHE_LIKE] = ipHe2p1P; 00456 00457 /* he-like redistribution functions */ 00459 iso.ipLyaRedist[ipHE_LIKE] = ipPRD; 00460 iso.ipResoRedist[ipHE_LIKE] = ipCRD; 00461 iso.ipSubRedist[ipHE_LIKE] = ipCRDW; 00462 00463 /* do not average collision strengths - evaluate at kT 00464 * set true with command SET COLLISION STRENGHTS AVERAGE */ 00465 iso.lgCollStrenThermAver = false; 00466 00467 00468 /********************************************************************** 00469 * all parameters having to do with secondary ionization 00470 * by suprathermal electrons 00471 **********************************************************************/ 00472 secondaries.SetCsupra = 0.; 00473 secondaries.lgCSetOn = false; 00474 secondaries.lgSecOFF = false; 00475 secondaries.SecHIonMax = 0.; 00476 00477 secondaries.HeatEfficPrimary = 1.; 00478 secondaries.SecIon2PrimaryErg = 0.; 00479 secondaries.SecExcitLya2PrimaryErg = 0.; 00480 secondaries.x12tot = 0.; 00481 secondaries.sec2total = 0.; 00482 00483 if( lgFirstCall ) 00484 { 00485 /* malloc space for supra[nelem][ion] */ 00486 secondaries.csupra = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) ); 00487 secondaries.csupra_effic = (realnum **)MALLOC( (unsigned)LIMELM*sizeof(realnum *) ); 00488 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00489 { 00490 secondaries.csupra[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) ); 00491 secondaries.csupra_effic[nelem] = (realnum *)MALLOC( (unsigned)(nelem+1)*sizeof(realnum) ); 00492 } 00493 } 00494 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00495 { 00496 for( ion=0; ion<nelem+1; ++ion ) 00497 { 00498 /* secondary ionization rate for each species */ 00499 secondaries.csupra[nelem][ion] = 0.; 00500 /* the rate of each species relative to H0 */ 00501 secondaries.csupra_effic[nelem][ion] = 1.f; 00502 } 00503 } 00504 /* this scale factor is from table 10 of Tielens & Hollenbach 1985 */ 00505 secondaries.csupra_effic[ipHELIUM][0] = 1.08f; 00506 00507 /* on first call, these arrays do not exist, only zero here on 00508 * second and later calls, on first call, create them */ 00509 if( lgFirstCall ) 00510 { 00511 /* these will save bound electron recoil information data */ 00512 ionbal.ipCompRecoil = 00513 (long**)MALLOC(sizeof(long*)*(unsigned)LIMELM ); 00514 ionbal.CompRecoilIonRate = 00515 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00516 ionbal.CompRecoilIonRateSave = 00517 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00518 ionbal.CompRecoilHeatRate = 00519 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00520 ionbal.CompRecoilHeatRateSave = 00521 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00522 ionbal.PhotoRate_Shell = 00523 (double****)MALLOC(sizeof(double***)*(unsigned)LIMELM ); 00524 ionbal.CollIonRate_Ground = 00525 (double***)MALLOC(sizeof(double**)*(unsigned)LIMELM ); 00526 ionbal.UTA_ionize_rate = 00527 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00528 ionbal.UTA_heat_rate = 00529 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00530 00531 /* these are source and sink terms for heavy element ionization balance from the 00532 * chemistry */ 00533 mole.source = 00534 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00535 mole.sink = 00536 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ); 00537 mole.xMoleChTrRate = 00538 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM ); 00539 00540 /* space for ionization recombination arrays */ 00541 ionbal.RateIonizTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00542 ionbal.RateRecomTot = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00543 ionbal.RR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00544 ionbal.DR_rate_coef_used = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00545 ionbal.RR_Verner_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00546 00547 /* rate coefficients [cm3 s-1] for Badnell DR recombination */ 00548 ionbal.DR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00549 ionbal.DR_Badnell_rate_coef_mean_ion = (double *)MALLOC(sizeof(double)*(unsigned)LIMELM ); 00550 /* do these rate coefficients exist? */ 00551 ionbal.lgDR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM ); 00552 ionbal.RR_Badnell_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00553 /* do these rate coefficients exist? */ 00554 ionbal.lgRR_Badnell_rate_coef_exist = (int **)MALLOC(sizeof(int *)*(unsigned)LIMELM ); 00555 /* rate coefficients [cm3 s-1] for older DR recombination */ 00556 ionbal.DR_old_rate_coef = (double **)MALLOC(sizeof(double *)*(unsigned)LIMELM ); 00557 00558 /* create arrays for ions */ 00559 for(nelem=0; nelem<LIMELM; ++nelem ) 00560 { 00561 ionbal.DR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00562 ionbal.lgDR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) ); 00563 ionbal.RR_Badnell_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00564 ionbal.lgRR_Badnell_rate_coef_exist[nelem] = (int *)MALLOC(sizeof(int)*(unsigned)(nelem+1) ); 00565 ionbal.DR_old_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00566 00567 ionbal.RateIonizTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00568 ionbal.RateRecomTot[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00569 ionbal.RR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00570 ionbal.DR_rate_coef_used[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00571 ionbal.RR_Verner_rate_coef[nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00572 ionbal.UTA_ionize_rate[nelem] = 00573 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00574 ionbal.UTA_heat_rate[nelem] = 00575 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00576 ionbal.ipCompRecoil[nelem] = 00577 (long*)MALLOC(sizeof(long)*(unsigned)(nelem+1) ); 00578 ionbal.CompRecoilIonRate[nelem] = 00579 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00580 ionbal.CompRecoilIonRateSave[nelem] = 00581 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00582 ionbal.CompRecoilHeatRate[nelem] = 00583 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00584 ionbal.CompRecoilHeatRateSave[nelem] = 00585 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+1) ); 00586 ionbal.PhotoRate_Shell[nelem] = 00587 (double***)MALLOC(sizeof(double**)*(unsigned)(nelem+1) ); 00588 ionbal.CollIonRate_Ground[nelem] = 00589 (double**)MALLOC(sizeof(double*)*(unsigned)(nelem+1) ); 00590 /* chemistry source and sink terms for ionization ladders */ 00591 mole.source[nelem] = 00592 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) ); 00593 mole.sink[nelem] = 00594 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) ); 00595 mole.xMoleChTrRate[nelem] = 00596 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) ); 00597 for( ion=0; ion<nelem+2; ++ion ) 00598 { 00599 mole.xMoleChTrRate[nelem][ion] = 00600 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) ); 00601 } 00602 00603 for( ion=0; ion<nelem+1; ++ion ) 00604 { 00605 /* >>chng 03 aug 09, set these to impossible values */ 00606 ionbal.RateIonizTot[nelem][ion] = -1.; 00607 ionbal.RateRecomTot[nelem][ion] = -1.; 00608 ionbal.UTA_ionize_rate[nelem][ion] = -1.; 00609 ionbal.UTA_heat_rate[nelem][ion] = -1.; 00610 ionbal.ipCompRecoil[nelem][ion] = -99; 00611 ionbal.CompRecoilIonRate[nelem][ion] = -1.; 00612 ionbal.CompRecoilIonRateSave[nelem][ion] = -1.; 00613 ionbal.CompRecoilHeatRate[nelem][ion] = -1.; 00614 ionbal.CompRecoilHeatRateSave[nelem][ion] = -1.; 00615 00616 /* finish mallocing space */ 00617 ionbal.PhotoRate_Shell[nelem][ion] = 00618 (double**)MALLOC(sizeof(double*)*(unsigned)NSHELLS ); 00619 ionbal.CollIonRate_Ground[nelem][ion] = 00620 (double*)MALLOC(sizeof(double)*(unsigned)2 ); 00621 for( ns=0; ns<NSHELLS; ++ns ) 00622 { 00623 ionbal.PhotoRate_Shell[nelem][ion][ns] = 00624 (double*)MALLOC(sizeof(double)*(unsigned)3 ); 00625 } 00626 00627 /* now set to impossible values */ 00628 ionbal.ipCompRecoil[nelem][ion] = -100000; 00629 ionbal.DR_Badnell_rate_coef[nelem][ion] = 0.; 00630 ionbal.RR_Badnell_rate_coef[nelem][ion] = 0.; 00631 ionbal.DR_old_rate_coef[nelem][ion] = 0.; 00632 } 00633 00634 set_NaN( ionbal.DR_rate_coef_used[nelem], nelem+1 ); 00635 set_NaN( ionbal.RR_rate_coef_used[nelem], nelem+1 ); 00636 set_NaN( ionbal.RR_Verner_rate_coef[nelem], nelem+1 ); 00637 } 00638 } 00639 00640 for(ion=0; ion<LIMELM; ++ion ) 00641 { 00642 ionbal.DR_Badnell_rate_coef_mean_ion[ion] = 0.; 00643 } 00644 00645 /* now zero out these arrays */ 00646 for( nelem=0; nelem< LIMELM; ++nelem ) 00647 { 00648 for( ion=0; ion<nelem+1; ++ion ) 00649 { 00650 00651 ionbal.CompRecoilHeatRate[nelem][ion] = 0.; 00652 ionbal.CompRecoilIonRate[nelem][ion] = 0.; 00653 ionbal.UTA_ionize_rate[nelem][ion] = 0.; 00654 ionbal.UTA_heat_rate[nelem][ion] = 0.; 00655 ionbal.CollIonRate_Ground[nelem][ion][0] = 0.; 00656 ionbal.CollIonRate_Ground[nelem][ion][1] = 0.; 00657 ionbal.RateRecomTot[nelem][ion] = 0.; 00658 for( ns=0; ns < NSHELLS; ++ns ) 00659 { 00660 /* must be zero since ion routines use these when 00661 * not yet defined */ 00662 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = 0.; 00663 ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.; 00664 ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.; 00665 } 00666 } 00667 /* these have one more ion than above */ 00668 for( ion=0; ion<nelem+2; ++ion ) 00669 { 00670 long int ion2; 00671 /* zero out the source and sink arrays */ 00672 mole.source[nelem][ion] = 0.; 00673 mole.sink[nelem][ion] = 0.; 00674 for( ion2=0; ion2<nelem+2; ++ion2 ) 00675 { 00676 mole.xMoleChTrRate[nelem][ion][ion2] = 0.; 00677 } 00678 } 00679 } 00680 00681 /* capture of molecules onto grain surfaces - formation of ices 00682 * flag says to include this process - turned off with the 00683 * command NO GRAIN MOLECULES */ 00684 mole.lgGrain_mole_deplete = true; 00685 00686 ionbal.lgPhotoIoniz_On = true; 00687 ionbal.lgCompRecoil = true; 00688 00689 /* these three adjust the treatment of UTA ionization */ 00690 ionbal.lgInnerShellLine_on = true; 00691 /*>>chng 07 jan 25, turn Kisielius off by default - rates show very 00692 * ragged trend with Z, use Gu et al. 06 which are more regular */ 00693 ionbal.lgInnerShell_Kisielius = false; 00694 ionbal.lgInnerShell_Gu06 = true; 00695 00696 for( i=0; i<4; ++i ) 00697 { 00698 ionbal.GuessDiel[i] = 1.; 00699 } 00700 00701 /* default condition is burgess suppressed, Nussbaumer and Storey not */ 00702 ionbal.lgSupDie[0] = true; 00703 ionbal.lgSupDie[1] = false; 00704 00705 ionbal.lgNoCota = false; 00706 for( i=0; i < LIMELM; i++ ) 00707 { 00708 ionbal.CotaRate[i] = 0.; 00709 } 00710 ionbal.ilt = 0; 00711 ionbal.iltln = 0; 00712 ionbal.ilthn = 0; 00713 ionbal.ihthn = 0; 00714 ionbal.ifail = 0; 00715 ionbal.lgGrainIonRecom = true; 00716 00717 /*>>chng 06 nov 29, use Badnell DR by default */ 00718 /* do not turn on Badnell DR by default */ 00719 ionbal.lgDR_recom_Badnell_use = false; 00720 /* turn on Badnell DR by default */ 00721 ionbal.lgDR_recom_Badnell_use = true; 00722 00723 /* >>chng 06 nov 23, use RR by default */ 00724 /* do not turn on Badnell RR by default */ 00725 ionbal.lgRR_recom_Badnell_use = false; 00726 /* turn on Badnell RR by default */ 00727 ionbal.lgRR_recom_Badnell_use = true; 00728 00729 /* option to print recombination coefficient then exit */ 00730 ionbal.lgRecom_Badnell_print = false; 00731 00732 /* this flag says how to generate S DR - 00733 * 0 == use mix of real Badnell DR and guess */ 00734 ionbal.nDR_S_guess = 0; 00735 00736 /* set true to use Badnell mean in place of existing hacks */ 00737 ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = false; 00738 ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = true; 00739 00740 /* setting this non-zero will turn on guess for dr for entire range of ions */ 00741 /* >>chng 03 nov 22, change from false to true - turn on process */ 00742 /* >>refer dr guess Kraemer, S.B., Ferland, G.J., & Gabel, J.R. 2004, ApJ in press */ 00743 ionbal.lg_guess_coef = true; 00744 ionbal.guess_noise = 0.; 00745 00746 /* do O H charge transfer in chemistry by default, changes with 00747 * set HO charge transfer */ 00748 ionbal.lgHO_ct_chem = true; 00749 00750 /********************************************************************** 00751 * these are options to print errors to special window, 00752 * set with print errors command, 00753 * output will go to standard error 00754 * defined in cdInit 00755 **********************************************************************/ 00756 lgPrnErr = false; 00757 ioPrnErr = stderr; 00758 00759 /* the code is ok at startup */ 00760 broke.lgBroke = false; 00761 broke.lgFixit = false; 00762 broke.lgCheckit = false; 00763 00764 /* main arrays to save ionization fractions*/ 00765 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00766 { 00767 dense.gas_phase[nelem] = 0.; 00768 dense.xMolecules[nelem] = 0.; 00769 for( ion=0; ion < LIMELM+1; ion++ ) 00770 { 00771 dense.xIonDense[nelem][ion] = 0.; 00772 } 00773 } 00774 dense.xMassTotal = 0.; 00775 00776 /* >> chng 05 26, 2006 NPA. Define the mass of each molecule for 00777 * use in non-equilibrium chemistry calculations */ 00778 00779 for(i=0;i<N_H_MOLEC;++i) 00780 { 00781 co.hmole_mass[i] = 0; 00782 } 00783 00784 /* Now define the molar mass of each molecule in the hydrogen 00785 * chemistry. This also includes He+, which also is a reactant 00786 * in the heavy element chemistry */ 00787 co.hmole_mass[0] = 1.0*ATOMIC_MASS_UNIT; 00788 co.hmole_mass[1] = 1.0*ATOMIC_MASS_UNIT; 00789 co.hmole_mass[2] = 1.0*ATOMIC_MASS_UNIT; 00790 co.hmole_mass[3] = 2.0*ATOMIC_MASS_UNIT; 00791 co.hmole_mass[4] = 2.0*ATOMIC_MASS_UNIT; 00792 co.hmole_mass[5] = 3.0*ATOMIC_MASS_UNIT; 00793 co.hmole_mass[6] = 2.0*ATOMIC_MASS_UNIT; 00794 co.hmole_mass[7] = 5.0*ATOMIC_MASS_UNIT; 00795 co.hmole_mass[8] = 4.0*ATOMIC_MASS_UNIT; 00796 00797 /* this is the simple Fred Hamann FeII atom */ 00798 t_fe2ovr_la::Inst().zero_opacity(); 00799 00800 /* zero out volume and column density save arrays */ 00801 MeanZero(); 00802 00803 /* zero out heating rates */ 00804 HeatZero(); 00805 00806 /* some parameters dealing with calcium */ 00807 ca.Ca2RmLya = 0.; 00808 ca.popca2ex = 0.; 00809 ca.oldcla = 0.; 00810 ca.Ca3d = 0.; 00811 ca.Ca4p = 0.; 00812 ca.dstCala = 0.; 00813 00814 /* C12/C13 isotope ratio, sets the ratio of C12O16 to C13O16 and 00815 * C13/C12 1909 */ 00816 co.C12_C13_isotope_ratio = 30.; 00817 00818 /* flag saying that H2O water destruction rate went to zero */ 00819 co.lgH2Ozer = false; 00820 00821 /* extra electron density, set with eden command */ 00822 dense.EdenExtra = 0.; 00823 00824 /* forced electron density, set with set eden command */ 00825 dense.EdenSet = 0.; 00826 00827 /* this is the default allowed relative error in the pressure */ 00828 conv.PressureErrorAllowed = 0.01f; 00829 00830 /* this is abort option set with SET PRESIONIZ command */ 00831 conv.limPres2Ioniz = 100000000; 00832 00833 /* some titles and line images */ 00834 for( i=0; i<74; ++i) 00835 { 00836 input.chTitle[i] = ' '; 00837 } 00838 input.chTitle[75] = '\0'; 00839 00840 /* velocity field information */ 00841 /* the turbulent velocity at illuminated face, internally in cm/s, 00842 * but entered with turbulence command in km/s */ 00843 DoppVel.TurbVel = 0.; 00844 /* the parameter F in eq 34 of 00845 *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, 624, 773 */ 00846 DoppVel.Heiles_Troland_F = 0.; 00847 /* is TurbVel included in pressure? - can be done two ways, with the velocity 00848 * being set of with equipartition - true when TurbVel set if not equipartition, 00849 * false with NO PRESSURE option on turbulence command */ 00850 DoppVel.lgTurb_pressure = true; 00851 /* The scale in cm over which the turbulence is dissipated. Normally 0, 00852 * only set if dissipate keyword appears on turbulence command */ 00853 DoppVel.DispScale = 0.; 00854 /* equipartition option on turbulence command, to set turbulence from B */ 00855 DoppVel.lgTurbEquiMag = false; 00856 00857 /* pressure related variables */ 00858 pressure.PresInteg = 0.; 00859 pressure.pinzon = 0.; 00860 00861 pressure.PresRamCurr = 0.; 00862 pressure.pres_radiation_lines_curr = 0.; 00863 pressure.lgPradCap = false; 00864 pressure.lgPradDen = false; 00865 pressure.lgLineRadPresOn = true; 00866 /* normally abort when radiation pressure exceeds gas pressure in const pres mod, 00867 * this is option to keep going, set with NO ABORT on constant pressure command */ 00868 pressure.lgRadPresAbortOK = true; 00869 /* Ditto for whether to stop at sonic point, this gets set to false 00870 * for some of the dynamics pressure modes (strongd, shock, antishock)*/ 00871 pressure.lgSonicPointAbortOK = true; 00872 /* this flag will say we hit the sonic point */ 00873 pressure.lgSonicPoint = false; 00874 /* True when no physical solution for desired pressure in strong D fronts */ 00875 pressure.lgStrongDLimbo = false; 00876 00877 pressure.RadBetaMax = 0.; 00878 pressure.ipPradMax_nzone = -1; 00879 pressure.PresMax = 0.; 00880 00881 /* initial and current pressures */ 00882 pressure.PresTotlInit = 0.; 00883 pressure.PresTotlCurr = 0.; 00884 00885 /* zero out some dynamics variables */ 00886 DynaZero(); 00887 00888 phycon.lgPhysOK = true; 00889 /* largest relative changes in Te, ne, H+, H2, and CO in structure 00890 * this is computed as part of prtcomment so does not exist when code not talking, 00891 * set to zero in zero and still zero if prtcomment not called */ 00892 phycon.BigJumpTe = 0.; 00893 phycon.BigJumpne = 0.; 00894 phycon.BigJumpH2 = 0.; 00895 phycon.BigJumpCO = 0.; 00896 00897 dense.xNucleiTotal = 1.; 00898 /* WJH */ 00899 dense.xMassDensity0 = -1.0f; 00900 00901 dense.eden = 1.; 00902 00903 /* now set physical conditions array 00904 * following will force updating all temperature - density variables */ 00905 TempChange( 1e4 , true); 00906 00907 00908 /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the 00909 * electron density to account for collisions with atomic H. it is an order of 00910 * magnitude guess, so this command provides ability to see whether it affects results */ 00911 dense.HCorrFac = 1.f; 00912 00913 dense.H_sum_in_CO = 0.; 00914 00915 atoms.nNegOI = 0; 00916 for( i=0; i< N_OI_LEVELS; ++i ) 00917 { 00918 atoms.popoi[i] = 0.; 00919 } 00920 atoms.popmg2 = 0.; 00921 00922 /* do we want to punch negative opacities */ 00923 opac.lgNegOpacIO = false; 00924 00925 /* this is flag for turning on case b */ 00926 opac.lgCaseB = false; 00927 00928 /* this is separate flag for turning off collisions from n=2 */ 00929 opac.lgCaseB_HummerStorey = false; 00930 00931 /* this is separate flag for turning off excited state photoionization */ 00932 opac.lgCaseB_no_photo = false; 00933 /* another case b option, turn off background opacities, no Pdest */ 00934 opac.lgCaseB_no_pdest = false; 00935 00936 /* smallest allowed line and Lya optical depths, reset with 00937 * Case B command */ 00938 opac.taumin = 1e-20f; 00939 opac.tlamin = 1e-20f; 00940 00941 opac.otsmin = 0.; 00942 00943 /* this flag says to use the static opacities, 00944 * only evaluate them at start of each new zone. 00945 * when set false with 00946 * no static opacities 00947 * command, always reevaluate them */ 00948 opac.lgOpacStatic = true; 00949 00950 /* set true in radinc if negative opacities ever occur */ 00951 opac.lgOpacNeg = false; 00952 00953 /* can turn of scattering opacities for some geometries */ 00954 opac.lgScatON = true; 00955 00956 /* variables having to do with compiling and/or using the 00957 * ancillary file of stored opacities */ 00958 opac.lgCompileOpac = false; 00959 /* "no file opacity" command sets following var false, says not to use file */ 00961 opac.lgUseFileOpac = false; 00962 00963 opac.telec = opac.taumin; 00964 opac.thmin = opac.taumin; 00965 00966 opac.tpcah[0] = opac.taumin; 00967 opac.tpcah[1] = 1e20f; 00968 00969 dynamics.dDensityDT = 0.; 00970 /* effects of fast neutrons */ 00971 hextra.frcneu = 0.; 00972 hextra.effneu = 1.; 00973 hextra.totneu = 0.; 00974 hextra.lgNeutrnHeatOn = false; 00975 hextra.CrsSecNeutron = 4e-26; 00976 00977 opac.stimax[0] = 0.; 00978 opac.stimax[1] = 0.; 00979 00980 /* this is limit on ionization we shall check for zoning 00981 * and prtcomment */ 00982 struc.dr_ionfrac_limit = 1e-3f; 00983 struc.nzone = -1; 00984 for(i=0;i<N_H_MOLEC;i++) 00985 { 00986 /* number of protons in each molecule */ 00987 int hmi_protons[N_H_MOLEC] = {1,1,1,2,2,3,2,1}; 00988 /* >>chng 03 nov 28, add electrons to hmi struc */ 00989 /* number of electrons in each molecule, THAT IS NOT COUNTED AS ION, 00990 * used to get H mole electron sum */ 00991 int hmi_electrons[N_H_MOLEC] = {0,0,-1,0,1,1,0,1}; 00992 00993 hmi.Hmolec[i] = 0.; 00994 hmi.nProton[i] = hmi_protons[i]; 00995 /* number of free electrons contributed, -1 for H- */ 00996 hmi.nElectron[i] = hmi_electrons[i]; 00997 } 00998 00999 hmi.H2_total = 0.; 01000 hmi.H2_frac_abund_set = 0.; 01001 hmi.hmihet = 0.; 01002 hmi.h2plus_heat = 0.; 01003 hmi.HeatH2Dish_used = 0.; 01004 hmi.HeatH2Dexc_used = 0.; 01005 hmi.HeatH2Dish_TH85 = 0.; 01006 hmi.HeatH2Dexc_TH85 = 0.; 01007 hmi.HeatH2Dish_BigH2 = 0.; 01008 hmi.HeatH2Dexc_BigH2 = 0.; 01009 hmi.UV_Cont_rel2_Draine_DB96_face = 0.; 01010 hmi.UV_Cont_rel2_Draine_DB96_depth = 0.; 01011 hmi.UV_Cont_rel2_Habing_TH85_face = 0.; 01012 hmi.UV_Cont_rel2_Habing_TH85_depth = 0.; 01013 hmi.HeatH2DexcMax = 0.; 01014 hmi.CoolH2DexcMax = 0.; 01015 hmi.hmitot = 0.; 01016 hmi.H2Opacity = 0.; 01017 hmi.hmidep = 1.; 01018 hmi.h2dep = 1.; 01019 hmi.h2pdep = 1.; 01020 hmi.h3pdep = 1.; 01021 hmi.BiggestH2 = -1.f; 01022 /* option to turn off H molecules */ 01023 hmi.lgNoH2Mole = false; 01024 01025 /* option to kill effects of H2 in CO chemistry - set with 01026 * set Leiden hack h2* off */ 01027 hmi.lgLeiden_Keep_ipMH2s = true; 01028 hmi.lgLeidenCRHack = true; 01029 01030 /* option to turn on the UMIST rates, naturally this will be 1, set to zero 01031 with the set UMIST rates command */ 01032 co.lgUMISTrates = true; 01033 01034 /* option to use diffuse cloud chemical rates from Table 8 of 01035 * >> refer Federman, S. R. & Zsargo, J. 2003, ApJ, 589, 319 01036 * By default, this is false - changed with set chemistry command */ 01037 co.lgFederman = true; 01038 /* option to use effective temperature as defined in 01039 * >> refer Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319 01040 * By default, this is false - changed with set chemistry command */ 01041 co.lgNonEquilChem = false; 01045 co.lgProtElim = true; 01049 co.lgNeutrals = true; 01050 01051 /* this says which estimate of the rate of the Solomon process to use, 01052 * default is Tielens & Hollenbach 1985a, changed with 01053 * set h2 Solomon command, options are TH85 and BD96, 01054 * the second for the Bertoldi & Draine rates - they 01055 * differ by 1 dex. when large H2 turned on this is ignored */ 01056 /* the Tielens & Hollenbach 1985 treatment */ 01057 hmi.chH2_small_model_type = 'T'; 01058 /* the improved H2 formalism given by 01059 *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M 01060 >>refcon 1990, ApJ, 365, 620 */ 01061 hmi.chH2_small_model_type = 'H'; 01062 /* the Bertoldi & Draine 1996 treatment */ 01063 /* >>chng 03 nov 15, change default to BD96 */ 01064 hmi.chH2_small_model_type = 'B'; 01065 /* >>chng 05 dec 08, use the Elwert et al. approximations as the default */ 01066 hmi.chH2_small_model_type = 'E'; 01067 01068 /* set NaN */ 01069 set_NaN( hmi.HeatH2Dish_used ); 01070 set_NaN( hmi.HeatH2Dish_BigH2 ); 01071 set_NaN( hmi.HeatH2Dish_TH85 ); 01072 set_NaN( hmi.HeatH2Dish_BD96 ); 01073 set_NaN( hmi.HeatH2Dish_BHT90 ); 01074 set_NaN( hmi.HeatH2Dish_ELWERT ); 01075 01078 set_NaN( hmi.HeatH2Dexc_used ); 01079 set_NaN( hmi.HeatH2Dexc_BigH2 ); 01080 set_NaN( hmi.HeatH2Dexc_TH85 ); 01081 set_NaN( hmi.HeatH2Dexc_BD96 ); 01082 set_NaN( hmi.HeatH2Dexc_BHT90 ); 01083 set_NaN( hmi.HeatH2Dexc_ELWERT ); 01084 01086 set_NaN( hmi.deriv_HeatH2Dexc_used ); 01087 set_NaN( hmi.deriv_HeatH2Dexc_BigH2 ); 01088 set_NaN( hmi.deriv_HeatH2Dexc_TH85 ); 01089 set_NaN( hmi.deriv_HeatH2Dexc_BD96 ); 01090 set_NaN( hmi.deriv_HeatH2Dexc_BHT90 ); 01091 set_NaN( hmi.deriv_HeatH2Dexc_ELWERT ); 01092 01093 set_NaN( hmi.H2_Solomon_dissoc_rate_used_H2g ); 01094 set_NaN( hmi.H2_Solomon_dissoc_rate_BigH2_H2g ); 01095 set_NaN( hmi.H2_Solomon_dissoc_rate_TH85_H2g ); 01096 set_NaN( hmi.H2_Solomon_dissoc_rate_BHT90_H2g ); 01097 set_NaN( hmi.H2_Solomon_dissoc_rate_BD96_H2g ); 01098 set_NaN( hmi.H2_Solomon_dissoc_rate_ELWERT_H2g ); 01099 01100 set_NaN( hmi.H2_Solomon_dissoc_rate_used_H2s ); 01101 set_NaN( hmi.H2_Solomon_dissoc_rate_BigH2_H2s ); 01102 set_NaN( hmi.H2_Solomon_dissoc_rate_TH85_H2s ); 01103 set_NaN( hmi.H2_Solomon_dissoc_rate_BHT90_H2s ); 01104 set_NaN( hmi.H2_Solomon_dissoc_rate_BD96_H2s ); 01105 set_NaN( hmi.H2_Solomon_dissoc_rate_ELWERT_H2s ); 01106 01111 set_NaN( hmi.H2_photodissoc_used_H2g ); 01112 set_NaN( hmi.H2_photodissoc_used_H2s ); 01113 set_NaN( hmi.H2_photodissoc_BigH2_H2s ); 01114 set_NaN( hmi.H2_photodissoc_BigH2_H2g ); 01115 set_NaN( hmi.H2_photodissoc_ELWERT_H2g ); 01116 set_NaN( hmi.H2_photodissoc_ELWERT_H2s ); 01117 set_NaN( hmi.H2_photodissoc_TH85 ); 01118 set_NaN( hmi.H2_photodissoc_BHT90 ); 01119 01120 /* default grain formation pumping - Takahashi 2001 */ 01121 hmi.chGrainFormPump = 'T'; 01122 01123 /* set which approximation for Jura rate - Cazaux & Tielens 01124 * >>refer H2 form Cazaux, S., & Tielens, A.G.G.M., 2002, ApJ, 575, L29 */ 01125 hmi.chJura = 'C'; 01126 01127 /* scale factor to multiply Jura rate, set Jura rate command */ 01128 hmi.ScaleJura = 1.f; 01129 01130 /* binding energy for change in H2 population while on grain surface, 01131 * set with "set h2 Tad" command */ 01132 hmi.Tad = 800.; 01133 01134 /* some vars in the h2 molecule */ 01135 H2_Zero(); 01136 01137 /* zero out some column densities */ 01138 for( i=0; i < NCOLD; i++ ) 01139 { 01140 colden.colden[i] = 0.; 01141 } 01142 colden.He123S = 0.; 01143 colden.coldenH2_ov_vel = 0.; 01144 /* ortho and para column densities */ 01145 h2.ortho_colden = 0.; 01146 h2.para_colden = 0.; 01147 01148 /* F=0 and F=1 column densities of H0*/ 01149 colden.H0_21cm_upper =0; 01150 colden.H0_21cm_lower =0; 01151 01152 for( i=0; i < 5; i++ ) 01153 { 01154 colden.C2Pops[i] = 0.; 01155 colden.C2Colden[i] = 0.; 01156 /* pops and column density for SiII atom */ 01157 colden.Si2Pops[i] = 0.; 01158 colden.Si2Colden[i] = 0.; 01159 } 01160 for( i=0; i < 3; i++ ) 01161 { 01162 /* pops and column density for CI atom */ 01163 colden.C1Pops[i] = 0.; 01164 colden.C1Colden[i] = 0.; 01165 /* pops and column density for OI atom */ 01166 colden.O1Pops[i] = 0.; 01167 colden.O1Colden[i] = 0.; 01168 /* pops and column density for CIII atom */ 01169 colden.C3Pops[i] = 0.; 01170 } 01171 for( i=0; i < 4; i++ ) 01172 { 01173 /* pops and column density for CIII atom */ 01174 colden.C3Colden[i] = 0.; 01175 } 01176 01177 /* variables to do with Jeans mass and radius */ 01178 colden.TotMassColl = 0.; 01179 colden.tmas = 0.; 01180 colden.wmas = 0.; 01181 colden.rjnmin = FLT_MAX; 01182 colden.ajmmin = FLT_MAX; 01183 01184 /* >>chng 01 jul 26, moved next statement from below loop to avoid bug in gcc 2.95.3, PvH */ 01185 /* default diffuse fields is outward only */ 01186 strcpy( rfield.chDffTrns, "OU2" ); 01187 01188 /* three arrays that are needed to save tabulated continua, LIMSPC is 01189 * limit to number of continua that can be entered. This is very large 01190 * for simplicity, the second dimension for the array will only be 01191 * created in parse table when needed, so not waste memory */ 01192 if( lgFirstCall ) 01193 { 01194 rfield.tNuRyd = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) ); 01195 rfield.tslop = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) ); 01196 rfield.tFluxLog = (realnum**)MALLOC((size_t)(LIMSPC*sizeof(realnum*)) ); 01197 rfield.lgContMalloc = (int*)MALLOC((size_t)(LIMSPC*sizeof(int)) ); 01198 /* this flag will say that the above vector has been malloced into the continuum array 01199 * this can be returned once the continuum is generated */ 01200 /* >>chng 06 jul 16, move inside if-statement -> this fixes memory leak when optimizing, PvH */ 01201 for( i=0; i < LIMSPC; ++i ) 01202 { 01203 rfield.lgContMalloc[i] = false; 01204 } 01205 } 01206 01207 rfield_opac_zero( 0 , rfield.nupper ); 01208 01209 /* strcpy( rfield.chDffTrns, "OU2" ); */ 01210 rfield.lgOutOnly = true; 01211 rfield.lgUSphON = false; 01212 01213 /* flags for whether continuum is defined over all energies */ 01214 rfield.lgMMok = true; 01215 rfield.lgHPhtOK = true; 01216 rfield.lgXRayOK = true; 01217 rfield.lgGamrOK = true; 01218 01219 /* these are the low and high energy bounds of the continuum */ 01220 rfield.emm = 1.001e-8f; 01221 rfield.egamry = 7.354e6f; 01222 01223 rfield.nflux = rfield.nupper; 01224 01225 /* where cloud is optically thick to brems */ 01226 rfield.ipEnergyBremsThin = 0; 01227 rfield.EnergyBremsThin = 0.; 01228 01229 /* this is the faintest the high-energy tail of the continuum be */ 01230 rfield.FluxFaint = 1e-10f; 01231 01232 for( i=0; i < LIMSPC; i++ ) 01233 { 01234 /* this is set true if particular continuum source can vary with time 01235 * set true if TIME appears on intensity / luminosity command line */ 01236 rfield.lgTimeVary[i] = false; 01237 /* most continua enter as a beam rather than isotropic */ 01238 rfield.lgBeamed[i] = true; 01239 /* default energy range is H-ionizing radiation */ 01240 rfield.range[i][0] = HIONPOT; 01241 rfield.range[i][1] = rfield.egamry; 01242 rfield.ioTableRead[i].clear(); 01243 } 01244 rfield.comtot = 0.; 01245 rfield.comoff = 1.; 01246 rfield.cmcool = 0.; 01247 rfield.cinrat = 0.; 01248 rfield.extin_mag_B_point = 0.; 01249 rfield.extin_mag_V_point = 0.; 01250 rfield.extin_mag_B_extended = 0.; 01251 rfield.extin_mag_V_extended = 0.; 01252 rfield.EnerGammaRay = 7676.; 01253 01254 /* variables dealing with the radius */ 01255 radius.rinner = 0.; 01256 radius.distance = 0.; 01257 radius.Radius = 0.; 01258 radius.Radius_mid_zone = 0.; 01259 radius.depth = DEPTH_OFFSET; 01260 radius.depth_mid_zone = DEPTH_OFFSET/2.; 01261 radius.depth_x_fillfac = 0.; 01262 radius.lgRadiusKnown = false; 01263 radius.drad = 0.; 01264 radius.r1r0sq = 1.; 01265 /* this is changed with the roberto command, to go from out to in */ 01266 radius.dRadSign = 1.; 01267 radius.lgDrNeg = false; 01268 01269 /* RDFALT is log of default starting radius (cm) */ 01270 /* >>chng 03 nov 12, from 25 to 30 for Lya clouds */ 01271 /*radius.rdfalt = 25.;*/ 01272 radius.rdfalt = 30.; 01273 01274 /* set default cylinder thickness */ 01275 radius.CylindHigh = 1e35f; 01276 radius.lgCylnOn = false; 01277 01278 radius.drad_x_fillfac = 1.; 01279 radius.drad_x_fillfac_mean = 1.; 01280 radius.dVeff = 1.; 01281 radius.drNext = 1.; 01282 radius.dRNeff = 1.; 01283 radius.lgdR2Small = false; 01284 01285 radius.glbdst = 0.; 01286 radius.glbrad = 0.; 01287 01288 radius.sdrmin = SMALLFLOAT; 01289 radius.sdrmax = 1e30; 01290 radius.lgSMinON = false; 01291 radius.lgDrMnOn = true; 01292 01293 radius.lgDrMinUsed = false; 01294 01295 /* drChange was reset to get orion flux in h-beta correct 01296 * drChange is really tau of current zone */ 01297 radius.drChange = 0.15f; 01298 01299 rfield.lgHabing = false; 01300 /* this flag to make sure more than one table read does not occur in same run */ 01301 rfield.lgTableRead = false; 01302 01303 /* flag to turn off Lya ots */ 01304 rfield.lgLyaOTS = true; 01305 /* HeII rec and Lya ots */ 01306 rfield.lgHeIIOTS = true; 01307 rfield.lgKillOTSLine = false; 01308 rfield.lgKillOutLine = false; 01309 rfield.lgKillOutCont = false; 01310 01311 /* rfield.DiffPumpOn is unity unless process disabled by setting to 1 01312 * with no diffuse line pumping command */ 01313 rfield.DiffPumpOn = 1.; 01314 01315 rfield.lgInducProcess = true; 01316 /* 02 jun 13, by Ryan...added this line */ 01317 rfield.lgCompileGauntFF = false; 01318 01319 /* >>chng 03 nov 28, add option to not do line transfer */ 01320 rfield.lgDoLineTrans = true; 01321 01322 /* flag saying whether to constantly reevaluated opacities - 01323 * set false with no opacity reevaluate command */ 01324 rfield.lgOpacityReevaluate = true; 01325 01326 /* flag saying whether to constantly reevaluated ionization - 01327 * set false with no ionization reevaluate command */ 01328 rfield.lgIonizReevaluate = true; 01329 01330 /* this flag says that CMB has been set */ 01331 rfield.lgCMB_set = false; 01332 01333 /* line overlap opacity, turn off with no fine opacity command */ 01334 rfield.lgOpacityFine = true; 01335 /* this element is default for choosing line width */ 01336 rfield.fine_opac_nelem = ipIRON; 01337 /* there will be this many resolution elements in one FWHM for this element, 01338 * at the lowest temperature to be considered */ 01339 rfield.fine_opac_nresolv = 1; 01340 /* continuum scale factor for case of time varying continuum */ 01341 rfield.time_continuum_scale = 1.; 01342 /* will fine optical depths be punched? */ 01343 rfield.lgPunchOpacityFine = false; 01344 01345 /* first is set true if one of the incident continua needs to have 01346 * H-ionizing radiation blocked. Second is set true is it is blocked 01347 * with extinguish command - want both true if first is true */ 01348 rfield.lgMustBlockHIon = false; 01349 rfield.lgBlockHIon = false; 01350 01351 /* reset some variable related to cooling */ 01352 CoolZero(); 01353 01354 thermal.lgCNegChk = true; 01355 thermal.CoolHeatMax = 0.; 01356 thermal.wlCoolHeatMax = 0; 01357 thermal.totcol = 0.; 01358 thermal.heatl = 0.; 01359 thermal.coolheat = 0.; 01360 thermal.lgCExtraOn = false; 01361 thermal.CoolExtra = 0.; 01362 thermal.ctot = 1.; 01363 01364 thermal.HeatNet = 0.; 01365 thermal.htot = 1.; 01366 thermal.power = 0.; 01367 thermal.FreeFreeTotHeat = 0.; 01368 thermal.char_tran_cool = 0.; 01369 thermal.char_tran_heat = 0.; 01370 01371 fnzone = 0.; 01372 nzone = 0; 01373 /* save initial condition for talk in case PRINT LAST used */ 01374 called.lgTalkSave = called.lgTalk; 01375 01376 oxy.poiii2 = 0.; 01377 oxy.poiii3 = 0.; 01378 oxy.poiexc = 0.; 01379 01380 oxy.d5007r = 0.; 01381 oxy.d5007t = 0.; 01382 oxy.d4363 = 0.; 01383 oxy.d6300 = 0.; 01384 01385 atmdat.nsbig = 0; 01386 01387 /* option to turn off collisional ionization with "no collisional ionization" cmmnd */ 01388 atmdat.lgCollIonOn = true; 01389 01390 /*************************************************** 01391 * charge transfer ionization and recombination 01392 ***************************************************/ 01393 /* HCharHeatMax, HCharCoolMax are largest fractions of heating in cur zone 01394 * or cooling due to ct */ 01395 atmdat.HCharHeatMax = 0.; 01396 atmdat.HCharCoolMax = 0.; 01397 01398 atmdat.HIonFrac = 0.; 01399 atmdat.HCharExcIonTotal = 0.; 01400 atmdat.HIonFracMax = 0.; 01401 atmdat.HCharExcRecTotal = 0.; 01402 /* option to turn off all charge transfer, turned off with no charge transfer command */ 01403 atmdat.lgCTOn = true; 01404 01405 /* flag saying that charge transfer heating should be included, 01406 * turned off with no CTHeat commmand */ 01407 atmdat.HCharHeatOn = 1.; 01408 for( nelem=0; nelem< LIMELM; ++nelem ) 01409 { 01410 for( ion=0; ion<LIMELM; ++ion ) 01411 { 01412 atmdat.HCharExcIonOf[nelem][ion] = 0.; 01413 atmdat.HCharExcRecTo[nelem][ion] = 0.; 01414 } 01415 } 01416 01417 /* >>chng 97 jan 6, from 0 to 8.5e-10*q as per Alex Dalgarno e-mail 01418 * >>chng 97 feb 6, from 8.5e-10*q 1.92e-9x as per Alex Dalgarno e-mail */ 01419 atmdat.HCTAlex = 1.92e-9; 01420 01421 for( nelem=0; nelem < LIMELM; nelem++ ) 01422 { 01423 /* these are depletion scale factors */ 01424 abund.depset[nelem] = 1.; 01425 /*begin sanity check */ 01426 if( abund.depset[nelem] == 0. ) 01427 { 01428 fprintf( ioQQQ, " ZERO finds insane abundance or depletion.\n" ); 01429 fprintf( ioQQQ, " atomic number=%6ld abundance=%10.2e depletion=%10.2e\n", 01430 nelem, abund.solar[nelem], abund.depset[nelem] ); 01431 ShowMe(); 01432 cdEXIT(EXIT_FAILURE); 01433 } 01434 /*end sanity check */ 01435 } 01436 01437 /* typical ISM depletion factors, subjective mean of Cowie and Songaila 01438 * and Jenkins 01439 * */ 01440 abund.Depletion[0] = 1.; 01441 abund.Depletion[1] = 1.; 01442 abund.Depletion[2] = .16f; 01443 abund.Depletion[3] = .6f; 01444 abund.Depletion[4] = .13f; 01445 abund.Depletion[5] = 0.4f; 01446 abund.Depletion[6] = 1.0f; 01447 abund.Depletion[7] = 0.6f; 01448 abund.Depletion[8] = .3f; 01449 abund.Depletion[9] = 1.f; 01450 abund.Depletion[10] = 0.2f; 01451 abund.Depletion[11] = 0.2f; 01452 abund.Depletion[12] = 0.01f; 01453 abund.Depletion[13] = 0.03f; 01454 abund.Depletion[14] = .25f; 01455 abund.Depletion[15] = 1.0f; 01456 abund.Depletion[16] = 0.4f; 01457 abund.Depletion[17] = 1.0f; 01458 abund.Depletion[18] = .3f; 01459 abund.Depletion[19] = 1e-4f; 01460 abund.Depletion[20] = 5e-3f; 01461 abund.Depletion[21] = 8e-3f; 01462 abund.Depletion[22] = 6e-3f; 01463 abund.Depletion[23] = 6e-3f; 01464 abund.Depletion[24] = 5e-2f; 01465 abund.Depletion[25] = 0.01f; 01466 abund.Depletion[26] = 0.01f; 01467 abund.Depletion[27] = 0.01f; 01468 abund.Depletion[28] = .1f; 01469 abund.Depletion[29] = .25f; 01470 01471 abund.lgDepln = false; 01472 abund.ScaleMetals = 1.; 01473 01474 /* this tells the code to use standard Auger yields */ 01475 t_yield::Inst().reset_yield(); 01476 01477 rt.dTauMase = 0.; 01478 rt.lgMaserCapHit = false; 01479 rt.lgMaserSetDR = false; 01480 01481 rt.DoubleTau = 1.; 01482 rt.lgFstOn = true; 01483 rt.lgElecScatEscape = true; 01484 01485 /* there was a call to TestCode */ 01486 lgTestCodeCalled = false; 01487 /* test code enabled with set test command */ 01488 lgTestCodeEnabled = false; 01489 01490 /* zero out some grain variables */ 01491 GrainZero(); 01492 01493 /* following should never be set non-zero */ 01494 FeII.fe21406 = 0.; 01495 FeII.fe21507 = 0.; 01496 FeII.fe21508 = 0.; 01497 FeII.fe21609 = 0.; 01498 FeII.fe21308 = 0.; 01499 FeII.fe21207 = 0.; 01500 FeII.fe21106 = 0.; 01501 FeII.fe21006 = 0.; 01502 FeII.fe21204 = 0.; 01503 FeII.fe21103 = 0.; 01504 FeII.fe21104 = 0.; 01505 FeII.fe21001 = 0.; 01506 FeII.fe21002 = 0.; 01507 FeII.fe20201 = 0.; 01508 FeII.fe20302 = 0.; 01509 FeII.fe20706 = 0.; 01510 FeII.fe20807 = 0.; 01511 FeII.fe20908 = 0.; 01512 FeII.fe21007 = 0.; 01513 FeII.fe21107 = 0.; 01514 FeII.fe21108 = 0.; 01515 FeII.fe21110 = 0.; 01516 FeII.fe21208 = 0.; 01517 FeII.fe21209 = 0.; 01518 FeII.fe21211 = 0.; 01519 01520 fe.Fe4CoolTot = 0.; 01521 fe.fe40401 = 0.; 01522 fe.fe42836 = 0.; 01523 fe.fe42829 = 0.; 01524 fe.fe42567 = 0.; 01525 fe.fe41207 = 0.; 01526 fe.fe41206 = 0.; 01527 fe.fe41106 = 0.; 01528 fe.fe41007 = 0.; 01529 fe.fe41008 = 0.; 01530 fe.fe40906 = 0.; 01531 01532 fe.Fe7CoolTot = 0.; 01533 # if 0 01534 fe.Fe7_6087 = 0.; 01535 fe.Fe7_5721 = 0.; 01536 fe.Fe7_6601 = 0.; 01537 fe.Fe7_3760 = 0.; 01538 fe.Fe7_3588 = 0.; 01539 # endif 01540 01541 /* this is flag saying whether this is very first call, 01542 * a time when space has not been allocated */ 01543 lgFirstCall = false; 01544 return; 01545 } 01546 01547 /*rfield_opac_zero zero out rfield arrays between certain limits */ 01548 void rfield_opac_zero( 01549 /* index for first element in arrays to be set to zero */ 01550 long lo , 01551 /* array index for highest element to be set */ 01552 long ihi ) 01553 { 01554 long int i; 01555 01556 /* >>chng 01 aug 19, space not allocated yet, 01557 * following code must also be present in contcreatemesh where 01558 * space allocated for the first time */ 01559 if( lgRfieldMalloced ) 01560 { 01561 unsigned long n=(unsigned long)(ihi-lo+1); 01562 memset(&rfield.OccNumbDiffCont[lo] , 0 , n*sizeof(realnum) ); 01563 memset(&rfield.OccNumbContEmitOut[lo] , 0 , n*sizeof(realnum) ); 01564 memset(&rfield.ContBoltz[lo] , 0 , n*sizeof(double) ); 01565 /*>>chng 06 aug 15, this is now 2D array, saving diffuse continuum 01566 * over all zones for use in exact RT */ 01567 /*memset(&rfield.ConEmitLocal[lo] , 0 , n*sizeof(realnum) );*/ 01568 memset(&rfield.ConEmitReflec[lo] , 0 , n*sizeof(realnum) ); 01569 memset(&rfield.ConEmitOut[lo] , 0 , n*sizeof(realnum) ); 01570 memset(&rfield.reflin[lo] , 0 , n*sizeof(realnum) ); 01571 memset(&rfield.ConRefIncid[lo] , 0 , n*sizeof(realnum) ); 01572 memset(&rfield.SummedCon[lo] , 0 , n*sizeof(realnum) ); 01573 memset(&rfield.OccNumbBremsCont[lo] , 0 , n*sizeof(realnum) ); 01574 memset(&rfield.convoc[lo] , 0 , n*sizeof(realnum) ); 01575 memset(&rfield.flux[lo] , 0 , n*sizeof(realnum) ); 01576 memset(&rfield.flux_beam_const_save[lo] , 0 , n*sizeof(realnum) ); 01577 memset(&rfield.flux_time_beam_save[lo] , 0 , n*sizeof(realnum) ); 01578 memset(&rfield.flux_isotropic_save[lo] , 0 , n*sizeof(realnum) ); 01579 memset(&rfield.SummedOcc[lo] , 0 , n*sizeof(realnum) ); 01580 memset(&rfield.SummedDif[lo] , 0 , n*sizeof(realnum) ); 01581 memset(&rfield.flux_accum[lo] , 0 , n*sizeof(realnum) ); 01582 memset(&rfield.otslin[lo] , 0 , n*sizeof(realnum) ); 01583 memset(&rfield.otscon[lo] , 0 , n*sizeof(realnum) ); 01584 memset(&rfield.ConInterOut[lo] , 0 , n*sizeof(realnum) ); 01585 memset(&rfield.outlin[lo] , 0 , n*sizeof(realnum) ); 01586 memset(&rfield.outlin_noplot[lo] , 0 , n*sizeof(realnum) ); 01587 memset(&rfield.ConOTS_local_OTS_rate[lo], 0 , n*sizeof(realnum) ); 01588 memset(&rfield.ConOTS_local_photons[lo] , 0 , n*sizeof(realnum) ); 01589 memset(&rfield.otsconNew[lo] , 0 , n*sizeof(realnum) ); 01590 memset(&rfield.otslinNew[lo] , 0 , n*sizeof(realnum) ); 01591 memset(&opac.OldOpacSave[lo] , 0 , n*sizeof(double) ); 01592 memset(&opac.opacity_abs[lo] , 0 , n*sizeof(double) ); 01593 memset(&opac.opacity_sct[lo] , 0 , n*sizeof(double) ); 01594 memset(&opac.albedo[lo] , 0 , n*sizeof(double) ); 01595 memset(&opac.FreeFreeOpacity[lo] , 0 , n*sizeof(double) ); 01596 01597 /* these are not defined on first iteration */ 01598 memset( &opac.E2TauAbsTotal[lo] , 0 , n*sizeof(realnum) ); 01599 memset( &opac.E2TauAbsOut[lo] , 0 , n*sizeof(realnum) ); 01600 memset( &opac.TauAbsTotal[lo] , 0 , n*sizeof(realnum) ); 01601 01602 for( i=lo; i <= ihi; i++ ) 01603 { 01604 opac.TauTotalGeo[0][i] = opac.taumin; 01605 opac.TauAbsGeo[0][i] = opac.taumin; 01606 opac.TauScatGeo[0][i] = opac.taumin; 01607 opac.tmn[i] = 1.; 01608 opac.ExpZone[i] = 1.; 01609 opac.E2TauAbsFace[i] = 1.; 01610 opac.ExpmTau[i] = 1.; 01611 opac.OpacStatic[i] = 1.; 01612 } 01613 /* also zero out fine opacity fine grid fine mesh array */ 01614 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) ); 01615 /* also zero out fine opacity array */ 01616 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine_malloc*sizeof(realnum) ); 01617 } 01618 return; 01619 } 01620