cloudy trunk
Loading...
Searching...
No Matches
iter_startend.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/*IterStart set and save values of many variables at start of iteration after initial temperature set*/
4/*IterRestart restart iteration */
5#include "cddefines.h"
6#include "cddrive.h"
7#include "physconst.h"
8#include "iso.h"
9#include "taulines.h"
10#include "hydrogenic.h"
11#include "struc.h"
12#include "dynamics.h"
13#include "prt.h"
14#include "hyperfine.h"
15#include "dense.h"
16#include "magnetic.h"
17#include "continuum.h"
18#include "geometry.h"
19#include "h2.h"
20#include "co.h"
21#include "he.h"
22#include "grains.h"
23#include "atomfeii.h"
24#include "pressure.h"
25#include "stopcalc.h"
26#include "conv.h"
27#include "mean.h"
28#include "ca.h"
29#include "thermal.h"
30#include "atoms.h"
31#include "wind.h"
32#include "opacity.h"
33#include "timesc.h"
34#include "trace.h"
35#include "colden.h"
36#include "secondaries.h"
37#include "hmi.h"
38#include "radius.h"
39#include "phycon.h"
40#include "called.h"
41#include "mole.h"
42#include "rfield.h"
43#include "ionbal.h"
44#include "atmdat.h"
45#include "lines.h"
46#include "molcol.h"
47#include "input.h"
48#include "rt.h"
49#include "iterations.h"
50#include "cosmology.h"
51#include "deuterium.h"
52/* these are the static saved variables, set in IterStart and reset in IterRestart */
53
62
63/* arguments are atomic number, ionization stage*/
66static double HeatSave[LIMELM][LIMELM];
68/* save values of dr and drNext */
69static double drSave , drNextSave;
70
71/* also places to save them between iterations */
72static long int IonLowSave[LIMELM],
74
75/* these are used to reset the level populations of the h and he model atoms */
76/*static realnum hnsav[LIMELM][LMHLVL+1], */
77static bool lgHNSAV = false;
78
79static double ***HOpacRatSav;
80
81static double ortho_save , para_save;
82static double ***hnsav,
84
86static double *den_save;
87
88/*IterStart set and save values of many variables at start of iteration after initial temperature set*/
89void IterStart(void)
90{
91 long int i,
92 ion,
93 ion2,
94 ipISO,
95 n ,
96 nelem;
97 double fhe1nx,
98 ratio;
99
100 DEBUG_ENTRY( "IterStart()" );
101
102 /* allocate two matrices if first time in this core load
103 * these variables are malloced here because they are file static -
104 * used to save values at start of first zone, and reset to this value at
105 * start of new iteration */
106 if( !lgHNSAV )
107 {
108 /* set flag so we never do this again */
109 lgHNSAV = true;
110
111 HOpacRatSav = (double ***)MALLOC(sizeof(double **)*NISO );
112
113 hnsav = (double ***)MALLOC(sizeof(double **)*NISO );
114
115 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
116 {
117
118 HOpacRatSav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );
119
120 hnsav[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM );
121
122 /* now do the second dimension */
123 for( nelem=ipISO; nelem<LIMELM; ++nelem )
124 {
125 HOpacRatSav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );
126
127 hnsav[ipISO][nelem] = (double *)MALLOC(sizeof(double)*(unsigned)(iso_sp[ipISO][nelem].numLevels_max+1) );
128 }
129 }
131 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
132 saveMoleSink =
133 (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM );
135 (realnum***)MALLOC(sizeof(realnum**)*(unsigned)LIMELM );
136 for(nelem=0; nelem<LIMELM; ++nelem )
137 {
138 /* chemistry source and sink terms for ionization ladders */
139 saveMoleSource[nelem] =
140 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
141 saveMoleSink[nelem] =
142 (double*)MALLOC(sizeof(double)*(unsigned)(nelem+2) );
143 SaveMoleChTrRate[nelem] =
144 (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(nelem+2) );
145 for( ion=0; ion<nelem+2; ++ion )
146 {
147 SaveMoleChTrRate[nelem][ion] =
148 (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2) );
149 }
150 }
151 den_save = (double*)MALLOC(sizeof(double)*(mole_global.num_calc) );
152 }
153
154 /* IterStart is called at the start of EVERY iteration */
155 if( trace.lgTrace )
156 {
157 fprintf( ioQQQ, " IterStart called.\n" );
158 }
159
160 /* check if this is the last iteration */
161 if( iteration > iterations.itermx )
162 iterations.lgLastIt = true;
163
164 /* flag set true in RT_line_one_tauinc if maser ever capped */
165 rt.lgMaserCapHit = false;
166 /* flag for remembering if maser ever set dr in any zone */
167 rt.lgMaserSetDR = false;
168
169 /* zero out charge transfer heating and cooling fractions */
170 atmdat.HCharHeatMax = 0.;
171 atmdat.HCharCoolMax = 0.;
172
173 /* these are set false if limits of Case B tables is ever exceeded */
174 for(nelem=0; nelem<HS_NZ; ++nelem )
175 {
176 atmdat.lgHCaseBOK[0][nelem] = true;
177 atmdat.lgHCaseBOK[1][nelem] = true;
178 }
179
180 /* reset the largest number of levels in the database species */
181 for( long ipSpecies=0; ipSpecies<nSpecies; ++ipSpecies )
182 {
183 dBaseSpecies[ipSpecies].numLevels_local = dBaseSpecies[ipSpecies].numLevels_max;
184 dBaseSpecies[ipSpecies].CoolTotal = 0.;
185 }
186
187 /* the reason this calculation stops */
188 strncpy( StopCalc.chReasonStop, "reason not specified.",
189 sizeof(StopCalc.chReasonStop) );
190
191 /* zero fractions of He0 destruction due to 23S */
192 he.nzone = 0;
193 he.frac_he0dest_23S = 0.;
194 he.frac_he0dest_23S_photo = 0.;
195
196 dense.EdenMax = 0.;
197 dense.EdenMin = BIGFLOAT;
198
199 timesc.time_H2_Dest_longest = 0.;
200 timesc.time_H2_Form_longest = 0.;
201 /* remains neg if not evaluated */
202 timesc.time_H2_Dest_here = -1.;
203 timesc.time_H2_Form_here = 0.;
204
205 timesc.BigCOMoleForm = 0.;
206 timesc.TimeH21cm = 0.;
207 thermal.CoolHeatMax = 0.;
208 hydro.HCollIonMax = 0.;
209
210 atmdat.HIonFracMax = 0.;
211
212 /* will save highest n=2p/1s ratio for hydrogen atom*/
213 hydro.pop2mx = 0.;
214 hydro.lgHiPop2 = false;
215
216 hydro.nLyaHot = 0;
217 hydro.TLyaMax = 0.;
218
219 /* evaluate the gas and radiation pressure */
220 /* this sets values of pressure.PresTotlCurr */
221 pressure.PresInteg = 0.;
222 pressure.pinzon = 0.;
223 pressure.PresIntegElecThin = 0.;
225
226 dynamics.HeatMax = 0.;
227 dynamics.CoolMax = 0.;
228 if( dynamics.lgAdvection )
230
231 /* reset these since we now have a valid solution */
232 pressure.pbeta = 0.;
233 pressure.RadBetaMax = 0.;
234 pressure.lgPradCap = false;
235 pressure.lgPradDen = false;
236 /* this flag will say we hit the sonic point */
237 pressure.lgSonicPoint = false;
238 pressure.lgStrongDLimbo = false;
239
240 /* define two timescales for equilibrium, Compton and thermal */
241 timesc.tcmptn = 1.e0/(rfield.qtot*6.65e-25*dense.eden);
242 timesc.time_therm_long = 1.5*dense.pden*BOLTZMANN*phycon.te/thermal.ctot;
243
244 /* this will be total mass in computed structure */
245 dense.xMassTotal = 0.;
246
247 for( nelem=0; nelem < LIMELM; nelem++ )
248 {
249
250 if( dense.lgElmtOn[nelem] )
251 {
252
253 /* now zero out the ionic fractions */
254 for( ion=0; ion < (nelem + 2); ion++ )
255 {
256 xIonFsave[nelem][ion] = dense.xIonDense[nelem][ion];
257 }
258
259 for( ion=0; ion < LIMELM; ion++ )
260 {
261 HeatSave[nelem][ion] = thermal.heating[nelem][ion];
262 }
263 }
264 }
265
266 deutDenseSave0 = deut.xIonDense[0];
267 deutDenseSave1 = deut.xIonDense[1];
268
269 /* >>chng 02 jan 28, add ipISO loop */
270 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
271 {
272 /* >>chng 03 apr 11, from 0 to ipISO */
273 for( nelem=ipISO; nelem < LIMELM; nelem++ )
274 {
275 if( dense.lgElmtOn[nelem] )
276 {
277 for( n=0; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
278 {
279 HOpacRatSav[ipISO][nelem][n] = iso_sp[ipISO][nelem].fb[n].ConOpacRatio;
280 hnsav[ipISO][nelem][n] = iso_sp[ipISO][nelem].st[n].Pop();
281 }
282 }
283 for( vector<two_photon>::iterator tnu = iso_sp[ipISO][nelem].TwoNu.begin(); tnu != iso_sp[ipISO][nelem].TwoNu.end(); ++tnu )
284 tnu->induc_dn_max = 0.;
285 }
286 }
287
288 rfield.ipEnergyBremsThin = 0;
289 rfield.EnergyBremsThin = 0.;
290
291 p2nit = atoms.p2nit;
292 d5200r = atoms.d5200r;
293 atoms.nNegOI = 0;
294 for( i=0; i< N_OI_LEVELS; ++i )
295 atoms.popoi[i] = 0.;
296
297 /* save molecular fractions and ionization range */
298 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
299 {
300 /* save molecular densities */
301 gas_phase_save[nelem] = dense.gas_phase[nelem];
302 IonLowSave[nelem] = dense.IonLow[nelem];
303 IonHighSave[nelem] = dense.IonHigh[nelem];
304 }
305 /*fprintf(ioQQQ,"DEBUG IterStart save set to gas_phase hden %.4e \n",
306 dense.gas_phase[0]);*/
307
308 edsav = dense.eden;
309
310 /* save molecular densities */
311 for( i=0; i<mole_global.num_calc; i++)
312 {
313 den_save[i] = mole.species[i].den;
314 }
315 ortho_save = h2.ortho_density;
316 para_save = h2.para_density;
317
318 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
319 {
320 /* these have one more ion than above */
321 for( ion=0; ion<nelem+2; ++ion )
322 {
323 /* zero out the source and sink arrays */
324 saveMoleSource[nelem][ion] = mole.source[nelem][ion];
325 saveMoleSink[nelem][ion] = mole.sink[nelem][ion];
326 /*>>chng 06 jun 27, move from here, where leak occurs, to
327 * above where xMoleChTrRate first created */
328 /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
329 for( ion2=0; ion2<nelem+2; ++ion2 )
330 {
331 SaveMoleChTrRate[nelem][ion][ion2] = mole.xMoleChTrRate[nelem][ion][ion2];
332 }
333 }
334 }
335
336 hmihet_save = hmi.hmihet;
337 hmitot_save = hmi.hmitot;
338
339 h2plus_heat_save = hmi.h2plus_heat;
340
341 HeatH2Dish_used_save = hmi.HeatH2Dish_used;
342 HeatH2Dexc_used_save = hmi.HeatH2Dexc_used;
343
344 deriv_HeatH2Dexc_used_save = hmi.deriv_HeatH2Dexc_used;
345 H2_Solomon_dissoc_rate_used_H2g_save = hmi.H2_Solomon_dissoc_rate_used_H2g;
346 H2_Solomon_dissoc_rate_used_H2s_save = hmi.H2_Solomon_dissoc_rate_used_H2s;
347
348 H2_H2g_to_H2s_rate_used_save = hmi.H2_H2g_to_H2s_rate_used;
349
350 H2_photodissoc_used_H2s_save = hmi.H2_photodissoc_used_H2s;
351 H2_photodissoc_used_H2g_save = hmi.H2_photodissoc_used_H2g;
352
353 UV_Cont_rel2_Draine_DB96_face = hmi.UV_Cont_rel2_Draine_DB96_face;
354 UV_Cont_rel2_Draine_DB96_depth = hmi.UV_Cont_rel2_Draine_DB96_depth;
355 UV_Cont_rel2_Habing_TH85_face = hmi.UV_Cont_rel2_Habing_TH85_face;
356
357 /* save zone thickness, and next zone thickness */
358 drSave = radius.drad;
359 drNextSave = radius.drNext;
360 /* remember largest and smallest dr over iteration */
361 radius.dr_min_last_iter = BIGFLOAT;
362 radius.dr_max_last_iter = 0.;
363
364 geometry.nprint = iterations.IterPrnt[iteration-1];
365
366 colden.TotMassColl = 0.;
367 colden.tmas = 0.;
368 colden.wmas = 0.;
369 colden.rjnmin = FLT_MAX;
370 colden.ajmmin = FLT_MAX;
371
372 thermal.nUnstable = 0;
373 thermal.lgUnstable = false;
374
375 rfield.extin_mag_B_point = 0.;
376 rfield.extin_mag_V_point = 0.;
377 rfield.extin_mag_B_extended = 0.;
378 rfield.extin_mag_V_extended = 0.;
379
380 /* plus 1 is to zero out point where unit integration occurs */
381 for( i=0; i < rfield.nupper+1; i++ )
382 {
383 /* diffuse fields continuum */
384 /* >>chng 03 nov 08, recover SummedDif */
385 rfield.SummedDifSave[i] = rfield.SummedDif[i];
386 rfield.ConEmitReflec[0][i] = 0.;
387 rfield.ConEmitOut[0][i] = 0.;
388 rfield.ConInterOut[i] = 0.;
389 /* save OTS continuum for next iteration */
390 rfield.otssav[i][0] = rfield.otscon[i];
391 rfield.otssav[i][1] = rfield.otslin[i];
392 rfield.outlin[0][i] = 0.;
393 rfield.outlin_noplot[i] = 0.;
394 rfield.ConOTS_local_OTS_rate[i] = 0.;
395 rfield.ConOTS_local_photons[i] = 0.;
396 rfield.DiffuseEscape[i] = 0.;
397
398 /* save initial absorption, scattering opacities for next iteration */
399 opac.opacity_abs_savzon1[i] = opac.opacity_abs[i];
400 opac.opacity_sct_savzon1[i] = opac.opacity_sct[i];
401
402 /* will accumulate optical depths through cloud */
403 opac.TauAbsFace[i] = opac.taumin;
404 opac.TauScatFace[i] = opac.taumin;
405 /* >>chng 99 dec 04, having exactly 1 zone thickness for first zone caused discontinuity
406 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
407 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
408 * these were not here at all - changed to dr/2 */
409 /* attenuation of flux by optical depths IN THIS ZONE
410 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
411 * option for illumination of slab at an angle */
412 /* >>chng 04 oct 09, from drad to radius.drad_x_fillfac - include fill fac, PvH */
413 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin);
414
415 /* e(-tau) in inward direction, up to illuminated face */
416 opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
417
418 /* e2(tau) in inward direction, up to illuminated face, this is used to get the
419 * recombination escape probability */
420 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i] );
421 }
422
424
425 /* this zeros out arrays to hold mean ionization fractions
426 * later entered by mean
427 * read out by setlim */
428 mean.MeanZero();
429
430 /* zero out the column densities */
431 for( i=0; i < NCOLD; i++ )
432 {
433 colden.colden[i] = 0.;
434 }
435 colden.He123S = 0.;
436 colden.H0_ov_Tspin = 0.;
437 colden.OH_ov_Tspin = 0.;
438 colden.coldenH2_ov_vel = 0.;
439
440 /* upper and lower levels of H0 1s */
441 colden.H0_21cm_upper =0;
442 colden.H0_21cm_lower =0;
443
444 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
445 {
446 (*diatom)->ortho_colden = 0.;
447 (*diatom)->para_colden = 0.;
448 }
449
450 for( i=0; i < mole_global.num_calc; i++ )
451 {
452 mole.species[i].column = 0.;
453 }
454 for( i=0; i < 5; i++ )
455 {
456 colden.C2Pops[i] = 0.;
457 colden.C2Colden[i] = 0.;
458 /* pops and column density for SiII atom */
459 colden.Si2Pops[i] = 0.;
460 colden.Si2Colden[i] = 0.;
461 }
462 for( i=0; i < 3; i++ )
463 {
464 /* pops and column density for CI atom */
465 colden.C1Pops[i] = 0.;
466 colden.C1Colden[i] = 0.;
467 /* pops and column density for OI atom */
468 colden.O1Pops[i] = 0.;
469 colden.O1Colden[i] = 0.;
470 /* pops and column density for CIII atom */
471 colden.C3Pops[i] = 0.;
472 }
473 for( i=0; i < 4; i++ )
474 {
475 /* pops and column density for CIII atom */
476 colden.C3Colden[i] = 0.;
477 }
478
479 /* these are some line of sight emission measures */
480 colden.dlnenp = 0.;
481 colden.dlnenHep = 0.;
482 colden.dlnenHepp = 0.;
483 colden.dlnenCp = 0.;
484 colden.H0_ov_Tspin = 0.;
485 colden.OH_ov_Tspin = 0.;
486
487 // zero column densities of all states
488 for( unsigned i = 0; i < mole.species.size(); ++i )
489 {
490 if( mole.species[i].levels != NULL )
491 {
492 for( qList::iterator st = mole.species[i].levels->begin(); st != mole.species[i].levels->end(); ++st )
493 {
494 (*st).ColDen() = 0.;
495 }
496 }
497 }
498
499 /* now zero heavy element molecules */
500 molcol("ZERO",ioQQQ);
501
502 /* this will be sum of all free free heating over model */
503 thermal.FreeFreeTotHeat = 0.;
504
505 thermal.thist = 0.;
506 thermal.tlowst = 1e20f;
507
508 wind.AccelAver = 0.;
509 wind.acldr = 0.;
510 ionbal.ilt = 0;
511 ionbal.iltln = 0;
512 ionbal.ilthn = 0;
513 ionbal.ihthn = 0;
514 ionbal.ifail = 0;
515
516 secondaries.SecHIonMax = 0.;
517 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
518 {
519 for( ion=0; ion<nelem+1; ++ion )
520 {
521 supsav[nelem][ion] = secondaries.csupra[nelem][ion];
522 }
523 }
524 secondaries.x12sav = secondaries.x12tot;
525 secondaries.hetsav = secondaries.HeatEfficPrimary;
526 secondaries.savefi = secondaries.SecIon2PrimaryErg;
527 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
528 {
529 for( ion=0; ion<nelem+1; ++ion )
530 {
531 ionbal.CompRecoilHeatRateSave[nelem][ion] = ionbal.CompRecoilHeatRate[nelem][ion];
532 ionbal.CompRecoilIonRateSave[nelem][ion] = ionbal.CompRecoilIonRate[nelem][ion];
533 }
534 }
535
536 /* these will keep track of the number of convergence failures that occur */
537 conv.nTeFail = 0;
538 conv.nTotalFailures = 0;
539 conv.nPreFail = 0;
540 conv.failmx = 0.;
541 conv.nIonFail = 0;
542 conv.nPopFail = 0;
543 conv.nNeFail = 0;
544 conv.nGrainFail = 0;
545 conv.nChemFail = 0;
546 conv.dCmHdT = 0.;
547
548 /* abort flag */
549 lgAbort = false;
550
552
553 rfield.comtot = 0.;
554
555 co.codfrc = 0.;
556 co.codtot = 0.;
557
558 hmi.HeatH2DexcMax = 0.;
559 hmi.CoolH2DexcMax = 0.;
560 hmi.h2dfrc = 0.;
561 hmi.h2line_cool_frac = 0.;
562 hmi.h2dtot = 0.;
563 thermal.HeatLineMax = 0.;
564 thermal.GBarMax = 0.;
565 hyperfine.cooling_max = 0.;
566 hydro.cintot = 0.;
567 geometry.lgZoneTrp = false;
568
569 hmi.h2pmax = 0.;
570
571 /************************************************************************
572 *
573 * allocate space for lines arrays
574 *
575 ************************************************************************/
576
577
578
579 /* this was set in call to lines above */
580 ASSERT( LineSave.nsum > 0);
581 ASSERT( LineSave.nsumAllocated == LineSave.nsum );
582
583 /* zero emission line arrays - this has to be done on every iteration */
584 for( i=0; i < LineSave.nsum; i++ )
585 {
586 LineSv[i].SumLine[0] = 0.;
587 LineSv[i].SumLine[1] = 0.;
588 LineSv[i].emslin[0] = 0.;
589 LineSv[i].emslin[1] = 0.;
590 }
591
592 /* now zero out some variables set by last call to LINES */
593 hydro.cintot = 0.;
594 thermal.totcol = 0.;
595 rfield.comtot = 0.;
596 thermal.FreeFreeTotHeat = 0.;
597 thermal.power = 0.;
598 thermal.HeatLineMax = 0.;
599 thermal.GBarMax = 0.;
600 hyperfine.cooling_max = 0.;
601
602 hmi.h2pmax = 0.;
603
604 co.codfrc = 0.;
605 co.codtot = 0.;
606
607 hmi.h2dfrc = 0.;
608 hmi.h2line_cool_frac = 0.;
609 hmi.h2dtot = 0.;
610 timesc.sound = 0.;
611
612 if( LineSave.lgNormSet )
613 {
614 /* set normalization line index to zero from negative initial value so that
615 * we pass the assert in cdLine, in order to get the norm line index */
616 LineSave.ipNormWavL = 0;
617 LineSave.ipNormWavL = cdLine( LineSave.chNormLab , LineSave.WavLNorm , &fhe1nx , &ratio );
618 if( LineSave.ipNormWavL < 0 )
619 {
620 /* did not find the line if return is negative */
621 fprintf( ioQQQ, "PROBLEM could not find the normalisation line.\n");
622 fprintf( ioQQQ,
623 "IterStart could not find a line with a wavelength of %f and a label of %s\n",
624 LineSave.WavLNorm, LineSave.chNormLab );
625 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
626 fprintf( ioQQQ, "Sorry.\n");
628 }
629 }
630 else
631 {
632 /* set normalization line index to zero from negative initial value so that
633 * we pass the assert in cdLine, in order to get the norm line index */
634 LineSave.ipNormWavL = 0;
635 LineSave.ipNormWavL = cdLine( "TOTL" , 4861.36f , &fhe1nx , &ratio );
636 if( LineSave.ipNormWavL < 0 )
637 {
638 /* did not find the line if return is negative */
639 fprintf( ioQQQ, "PROBLEM could not find the default normalisation line.\n");
640 fprintf( ioQQQ,
641 "IterStart could not find a line with a wavelength of 4861 and a label of TOTL\n" );
642 fprintf( ioQQQ, "Please check the emission line output to find the correct line identification.\n");
643 fprintf( ioQQQ, "Sorry.\n");
645 }
646 }
647
648 /* set up stop line command on first iteration
649 * find index for lines and save for future iterations
650 * StopCalc.nstpl is zero (false) if no stop line commands entered */
651 if( iteration == 1 && StopCalc.nstpl )
652 {
653 /* nstpl is number of stop line commands, 0 if none entered */
654 for( long int nStopLine=0; nStopLine < StopCalc.nstpl; nStopLine++ )
655 {
656 double relint, absint ;
657 /* returns array index for line in array stack if we found the line,
658 * return negative of total number of lines as debugging aid if line not found */
659 StopCalc.ipStopLin1[nStopLine] = cdLine( StopCalc.chStopLabel1[nStopLine],
660 /* wavelength of line in angstroms, not format printed by code */
661 StopCalc.StopLineWl1[nStopLine], &relint, &absint );
662
663 if( StopCalc.ipStopLin1[nStopLine]<0 )
664 {
665 fprintf( ioQQQ,
666 " IterStart could not find first line in STOP LINE command, number %ld with label *%s* and wl %.1f\n",
667 StopCalc.ipStopLin1[nStopLine] ,
668 StopCalc.chStopLabel1[nStopLine],
669 StopCalc.StopLineWl1[nStopLine]);
671 }
672
673 StopCalc.ipStopLin2[nStopLine] = cdLine( StopCalc.chStopLabel2[nStopLine],
674 /* wavelength of line in angstroms, not format printed by code */
675 StopCalc.StopLineWl2[nStopLine], &relint, &absint );
676
677 if( StopCalc.ipStopLin2[nStopLine] < 0 )
678 {
679 fprintf( ioQQQ,
680 " IterStart could not find second line in STOP LINE command, line number %ld with label *%s* and wl %.1f\n",
681 StopCalc.ipStopLin1[nStopLine] ,
682 StopCalc.chStopLabel1[nStopLine],
683 StopCalc.StopLineWl1[nStopLine]);
685 }
686
687 if( trace.lgTrace )
688 {
689 fprintf( ioQQQ,
690 " stop line 1 is number %5ld wavelength is %f label is %4.4s\n",
691 StopCalc.ipStopLin1[nStopLine],
692 LineSv[StopCalc.ipStopLin1[nStopLine]].wavelength,
693 LineSv[StopCalc.ipStopLin1[nStopLine]].chALab );
694 fprintf( ioQQQ,
695 " stop line 2 is number %5ld wavelength is %f label is %4.4s\n",
696 StopCalc.ipStopLin2[nStopLine],
697 LineSv[StopCalc.ipStopLin2[nStopLine]].wavelength,
698 LineSv[StopCalc.ipStopLin2[nStopLine]].chALab );
699 }
700 }
701 }
702
703 /* option to only print last iteration */
704 if( prt.lgPrtLastIt )
705 {
706 if( iteration == 1 )
707 {
708 called.lgTalk = false;
709 }
710
711 /* initial condition of TALK may be off if optimization used or not master rank
712 * sec part is for print last command
713 * lgTalkForcedOff is set true when cdTalk is called with false
714 * to turn off printing */
715 if( iterations.lgLastIt && !prt.lgPrtStart && !called.lgTalkForcedOff )
716 {
717 called.lgTalk = called.lgTalkSave;
718 }
719 }
720
721 if( opac.lgCaseB )
722 {
723 if( trace.lgTrace )
724 {
725 fprintf( ioQQQ, " IterStart does not change mid-zone optical depth since CASE B\n" );
726 }
727 }
728
729 /* check if induced recombination can be ignored */
730 hydro.FracInd = 0.;
731 hydro.fbul = 0.;
732
733 /* remember some things about effects of induced rec on H only
734 * don't do ground state since SPHERE turns it off */
735 for( i=ipH2p; i < (iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max - 1); i++ )
736 {
737 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul() <= iso_ctrl.SmallA )
738 continue;
739
740 ratio = iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE/
741 SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RecomInducRate*iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].PopLTE +
742 iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecRad]*
743 iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].RadRecomb[ipRecNetEsc]);
744 if( ratio > hydro.FracInd )
745 {
746 hydro.FracInd = (realnum)ratio;
747 hydro.ndclev = i;
748 }
749
750 ratio = iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump()/
751 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().pump() +
752 iso_sp[ipH_LIKE][ipHYDROGEN].trans(i+1,i).Emis().Aul());
753
754 if( ratio > hydro.fbul )
755 {
756 hydro.fbul = (realnum)ratio;
757 hydro.nbul = i;
758 }
759 }
760
761 if( trace.lgTrace )
762 fprintf( ioQQQ, " IterStart returns.\n" );
763 return;
764}
765
766/*IterRestart reset many variables at the start of a new iteration
767 * called by cloudy after calculation is completed, when more iterations
768 * are needed - the iteration has been incremented when this routine is
769 * called so iteration == 2 after first iteration, we are starting
770 * the second iteration */
771void IterRestart(void)
772{
773 long int i,
774 ion,
775 ion2,
776 ipISO ,
777 n,
778 nelem;
779 double SumOTS;
780
781 DEBUG_ENTRY( "IterRestart()" );
782
783 /* this is case where temperature floor has been set, if it was hit
784 * then we did a constant temperature calculation, and must go back to
785 * a thermal solution
786 * test on thermal.lgTemperatureConstantCommandParsed distinguishes
787 * from temperature floor option, so not reset if constant temperature
788 * was actually set */
789 if( StopCalc.TeFloor > 0. && !thermal.lgTemperatureConstantCommandParsed )
790 {
791 thermal.lgTemperatureConstant = false;
792 thermal.ConstTemp = 0.;
793 }
794
795 lgAbort = false;
796
797 /* reset some parameters needed for magnetic field */
799
800 opac.stimax[0] = 0.;
801 opac.stimax[1] = 0.;
802
803 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
804 (*diatom)->H2_Reset();
805
806 for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
807 {
808 /* these have one more ion than above */
809 for( ion=0; ion<nelem+2; ++ion )
810 {
811 /* zero out the source and sink arrays */
812 mole.source[nelem][ion] = saveMoleSource[nelem][ion];
813 mole.sink[nelem][ion] = saveMoleSink[nelem][ion];
814 /*>>chng 06 jun 27, move from here, where leak occurs, to
815 * above where xMoleChTrRate first created */
816 /*mole.xMoleChTrRate[nelem][ion] = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(nelem+2));*/
817 for( ion2=0; ion2<nelem+2; ++ion2 )
818 {
819 mole.xMoleChTrRate[nelem][ion][ion2] = SaveMoleChTrRate[nelem][ion][ion2];
820 }
821 }
822 }
823
824 /* reset molecular abundances */
825 for( i=0; i<mole_global.num_calc; i++)
826 {
827 mole.species[i].den = den_save[i];
828 }
829 dense.updateXMolecules();
830 hmi.H2_total = findspecieslocal("H2")->den + findspecieslocal("H2*")->den;
831 hmi.HD_total = findspecieslocal("HD")->den;
832 /*fprintf(ioQQQ," IterRestar sets H2 total to %.2e\n",hmi.H2_total );*/
833 h2.ortho_density = ortho_save;
834 h2.para_density = para_save;
835 {
836 hmi.H2_total_f = (realnum)hmi.H2_total;
837 h2.ortho_density_f = (realnum)h2.ortho_density;
838 h2.para_density_f = (realnum)h2.para_density;
839 }
840
841 /* zero out the column densities */
842 for( i=0; i < NCOLD; i++ )
843 {
844 colden.colden[i] = 0.;
845 }
846 colden.He123S = 0.;
847 colden.coldenH2_ov_vel = 0.;
848
849 for( i=0; i < mole_global.num_calc; i++ )
850 {
851 /* column densities */
852 mole.species[i].column = 0.;
853 /* largest fraction of atoms in molecules */
854 mole.species[i].xFracLim = 0.;
855 mole.species[i].nAtomLim = -1;
856 }
857
858
859 /* close out this iteration if dynamics or time dependence is enabled */
860 if( dynamics.lgAdvection )
861 DynaIterEnd();
862
863 rfield.extin_mag_B_point = 0.;
864 rfield.extin_mag_V_point = 0.;
865 rfield.extin_mag_B_extended = 0.;
866 rfield.extin_mag_V_extended = 0.;
867
868 hmi.hmihet = hmihet_save;
869 hmi.hmitot = hmitot_save;
870
871 hmi.h2plus_heat = h2plus_heat_save;
872 hmi.HeatH2Dish_used = HeatH2Dish_used_save;
873 hmi.HeatH2Dexc_used = HeatH2Dexc_used_save;
874
875 hmi.deriv_HeatH2Dexc_used = (realnum)deriv_HeatH2Dexc_used_save;
876 hmi.H2_Solomon_dissoc_rate_used_H2g = H2_Solomon_dissoc_rate_used_H2g_save;
877 hmi.H2_Solomon_dissoc_rate_used_H2s = H2_Solomon_dissoc_rate_used_H2s_save;
878 hmi.H2_H2g_to_H2s_rate_used = H2_H2g_to_H2s_rate_used_save;
879 hmi.H2_photodissoc_used_H2s = H2_photodissoc_used_H2s_save;
880 hmi.H2_photodissoc_used_H2g = H2_photodissoc_used_H2g_save;
881
882 hmi.UV_Cont_rel2_Draine_DB96_face = (realnum)UV_Cont_rel2_Draine_DB96_face;
883 hmi.UV_Cont_rel2_Draine_DB96_depth = (realnum)UV_Cont_rel2_Draine_DB96_depth;
884 hmi.UV_Cont_rel2_Habing_TH85_face = (realnum)UV_Cont_rel2_Habing_TH85_face;
885
886 rfield.ipEnergyBremsThin = 0;
887 rfield.EnergyBremsThin = 0.;
888 rfield.lgUSphON = false;
889
890 radius.glbdst = 0.;
891 /* zone thickness, and next zone thickness */
892 radius.drad = drSave;
893 radius.drad_mid_zone = drSave/2.;
894 radius.drNext = drNextSave;
895
896 /* was min dr ever used? */
897 radius.lgDrMinUsed = false;
898
899 continuum.cn4861 = continuum.sv4861;
900 continuum.cn1216 = continuum.sv1216;
901
902 /* total luminosity */
903 continuum.totlsv = continuum.TotalLumin;
904
905 /* set debugger on now if NZONE desired is 0 */
906 if( (trace.nznbug == 0 && iteration >= trace.npsbug) && trace.lgTrOvrd )
907 {
908 if( trace.nTrConvg==0 )
909 {
910 trace.lgTrace = true;
911 }
912 else
913 /* trace convergence entered = but with negative flag = make positive,
914 * abs and not mult by -1 since may trigger more than one time */
915 trace.nTrConvg = abs( trace.nTrConvg );
916
917 fprintf( ioQQQ, " IterRestart called.\n" );
918 }
919 else
920 {
921 trace.lgTrace = false;
922 }
923
924 /* zero secondary suprathermals variables for first ionization attem */
925 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
926 {
927 for( ion=0; ion<nelem+1; ++ion )
928 {
929 secondaries.csupra[nelem][ion] = supsav[nelem][ion];
930 }
931 }
932 secondaries.x12tot = secondaries.x12sav;
933 secondaries.HeatEfficPrimary = secondaries.hetsav;
934 secondaries.SecIon2PrimaryErg = secondaries.savefi;
935 for( nelem=0; nelem<LIMELM; ++nelem)
936 {
937 for( ion=0; ion<nelem+1; ++ion )
938 {
939 ionbal.CompRecoilHeatRate[nelem][ion] = ionbal.CompRecoilHeatRateSave[nelem][ion];
940 ionbal.CompRecoilIonRate[nelem][ion] = ionbal.CompRecoilIonRateSave[nelem][ion];
941 }
942 }
943
944 wind.lgVelPos = true;
945 wind.AccelMax = 0.;
946 wind.AccelAver = 0.;
947 wind.acldr = 0.;
948 wind.windv = wind.windv0;
949
950 thermal.nUnstable = 0;
951 thermal.lgUnstable = false;
952
953 pressure.pbeta = 0.;
954 pressure.RadBetaMax = 0.;
955 pressure.lgPradCap = false;
956 pressure.lgPradDen = false;
957 /* this flag will say we hit the sonic point */
958 pressure.lgSonicPoint = false;
959 pressure.lgStrongDLimbo = false;
960
961 pressure.PresInteg = 0.;
962 pressure.pinzon = 0.;
963 pressure.PresIntegElecThin = 0.;
964
965 pressure.RhoGravity_dark = 0.;
966 pressure.RhoGravity_self = 0.;
967 pressure.RhoGravity_external = 0.;
968 pressure.RhoGravity = 0.;
969 pressure.IntegRhoGravity = 0.;
970
971 EdenChange( edsav );
972 dense.EdenHCorr = edsav;
973 dense.EdenHCorr_f = (realnum)dense.EdenHCorr;
974 dense.EdenTrue = edsav;
975
976 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
977 {
978 /* reset molecular densities */
979 dense.SetGasPhaseDensity( nelem, gas_phase_save[nelem] );
980 dense.IonLow[nelem] = IonLowSave[nelem];
981 dense.IonHigh[nelem] = IonHighSave[nelem];
982 }
983 /*fprintf(ioQQQ,"DEBUG IterRestart gas_phase set to save hden %.4e\n",
984 dense.gas_phase[0]);*/
985
986 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
987 {
988 for( long nelem=ipISO; nelem < LIMELM; ++nelem )
989 {
990 iso_sp[ipISO][nelem].Reset();
991 }
992 }
993
994 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
995 {
996 for( nelem=ipISO; nelem < LIMELM; nelem++ )
997 {
998 if( dense.lgElmtOn[nelem] )
999 {
1000 for( n=ipH1s; n < iso_sp[ipISO][nelem].numLevels_max; n++ )
1001 {
1002 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = HOpacRatSav[ipISO][nelem][n];
1003 iso_sp[ipISO][nelem].st[n].Pop() = hnsav[ipISO][nelem][n];
1004 }
1005 }
1006 }
1007 }
1008
1009 if( trace.lgTrace && trace.lgNeBug )
1010 {
1011 fprintf( ioQQQ, " EDEN set to%12.4e by IterRestart.\n",
1012 dense.eden );
1013 }
1014
1015 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
1016 {
1017 for( ion=0; ion < (nelem + 2); ion++ )
1018 {
1019 dense.xIonDense[nelem][ion] = xIonFsave[nelem][ion];
1020 }
1021 for( i=0; i < LIMELM; i++ )
1022 {
1023 thermal.heating[nelem][i] = HeatSave[nelem][i];
1024 }
1025 }
1026
1027 deut.xIonDense[0] = deutDenseSave0;
1028 deut.xIonDense[1] = deutDenseSave1;
1029
1031
1032 rfield.resetCoarseTransCoef();
1033
1034 /* continuum was saved in flux_total_incident */
1035 for( i=0; i < rfield.nupper; i++ )
1036 {
1037 /* time-constant part of beamed continuum */
1038 rfield.flux_beam_const[i] = rfield.flux_beam_const_save[i];
1039 /* continuum flux_time_beam_save has the initial value of the
1040 * time-dependent beamed continuum */
1041 rfield.flux_beam_time[i] = rfield.flux_time_beam_save[i]*
1042 rfield.time_continuum_scale;
1043
1044 if( cosmology.lgDo )
1045 {
1046 double slope_ratio;
1047 double fac_old = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_start ));
1048 double fac_new = TE1RYD*rfield.anu[i]/(CMB_TEMP*(1. + cosmology.redshift_current));
1049
1050 if( fac_old > log10(DBL_MAX) )
1051 {
1052 slope_ratio = 0.;
1053 }
1054 else if( fac_new > log10(DBL_MAX) )
1055 {
1056 slope_ratio = 0.;
1057 }
1058 else if( fac_old > 1.e-5 )
1059 {
1060 slope_ratio = (exp(fac_old) - 1.)/(exp(fac_new) - 1.);
1061 }
1062 else
1063 {
1064 slope_ratio = (fac_old*(1. + fac_old/2.))/(fac_new*(1. + fac_new/2.));
1065 }
1066
1067 rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i] * (realnum)slope_ratio;
1068 }
1069 else
1070 {
1071 /* the isotropic continuum source is time steady */
1072 rfield.flux_isotropic[i] = rfield.flux_isotropic_save[i];
1073 }
1074
1075 rfield.flux[0][i] = rfield.flux_isotropic[i] + rfield.flux_beam_time[i] +
1076 rfield.flux_beam_const[i];
1077 rfield.flux_total_incident[0][i] = rfield.flux[0][i];
1078
1079 rfield.SummedDif[i] = rfield.SummedDifSave[i];
1080 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
1081
1082 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
1083 rfield.otscon[i] = rfield.otssav[i][0];
1084 rfield.otslin[i] = rfield.otssav[i][1];
1085 rfield.ConInterOut[i] = 0.;
1086 rfield.OccNumbBremsCont[i] = 0.;
1087 rfield.OccNumbDiffCont[i] = 0.;
1088 rfield.OccNumbContEmitOut[i] = 0.;
1089 rfield.outlin[0][i] = 0.;
1090 rfield.outlin_noplot[i] = 0.;
1091 rfield.ConOTS_local_OTS_rate[i] = 0.;
1092 rfield.ConOTS_local_photons[i] = 0.;
1093 rfield.DiffuseEscape[i] = 0.;
1094
1095 opac.opacity_abs[i] = opac.opacity_abs_savzon1[i];
1096 opac.OldOpacSave[i] = opac.opacity_abs_savzon1[i];
1097 opac.opacity_sct[i] = opac.opacity_sct_savzon1[i];
1098 opac.albedo[i] =
1099 opac.opacity_sct[i]/MAX2(1e-30,opac.opacity_sct[i] + opac.opacity_abs[i]);
1100 opac.tmn[i] = 1.;
1101 /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
1102 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
1103 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
1104 opac.ExpmTau[i] = 1.;
1105 opac.E2TauAbsFace[i] = 1.;*/
1106 /* >>chng 99 dec 04, having exactly 1 for first zone caused discontinuity
1107 * for heating in very high T models in func_map.in. zone 1 and 2 were 20% different,
1108 * since tau in is 1e-20, e2 is 0.9999, and so some H ots
1109 * these were not here at all*/
1110 /* attenuation of flux by optical depths IN THIS ZONE
1111 * DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
1112 * option for illumination of slab at an angle */
1113 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad/2.*geometry.DirectionalCosin);
1114
1115 /* e(-tau) in inward direction, up to illuminated face */
1116 opac.ExpmTau[i] = (realnum)opac.ExpZone[i];
1117
1118 /* e2(tau) in inward direction, up to illuminated face */
1119 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
1120 rfield.reflin[0][i] = 0.;
1121 rfield.ConEmitReflec[0][i] = 0.;
1122 rfield.ConEmitOut[0][i] = 0.;
1123 rfield.ConRefIncid[0][i] = 0.;
1124
1125 /* escape in the outward direction
1126 * on second and later iterations define outward E2 */
1127 if( iteration > 1 )
1128 {
1129 /* e2 from current position to outer edge of shell */
1130 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
1131 opac.E2TauAbsOut[i] = (realnum)e2( tau );
1132 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
1133 }
1134 else
1135 opac.E2TauAbsOut[i] = 1.;
1136
1137 }
1138
1139 /* update continuum */
1140 RT_OTS_Update(&SumOTS);
1141
1142 thermal.FreeFreeTotHeat = 0.;
1143 atoms.p2nit = p2nit;
1144 atoms.d5200r = d5200r;
1145
1146 if( called.lgTalk )
1147 {
1148 fprintf( ioQQQ, "\f\n Start Iteration Number %ld %75.75s\n\n\n",
1149 iteration, input.chTitle );
1150 }
1151
1152 /* reset variable to do with FeII atom, in FeIILevelPops */
1153 FeIIReset();
1155 return;
1156}
1157
1158/* do some work with ending the iteration */
1159void IterEnd(void)
1160{
1161
1162 DEBUG_ENTRY( "IterEnd()" );
1163
1164 if( lgAbort )
1165 {
1166 return;
1167 }
1168
1169 /* give indication of geometry */
1170 double fac = radius.depth/radius.rinner;
1171 if( fac < 0.1 )
1172 {
1173 geometry.lgGeoPP = true;
1174 }
1175 else
1176 {
1177 geometry.lgGeoPP = false;
1178 }
1179
1180 if( iteration > dynamics.n_initial_relax && dynamics.lgTimeDependentStatic )
1181 {
1182 // report cumulative lines per unit mass rather than flux (per unit
1183 // area), so total is meaningful when density isn't constant
1184
1185 double cumulative_factor = (dynamics.timestep/
1186 colden.TotMassColl);
1187 // save cumulative lines
1188 for( long n=0; n<LineSave.nsum; ++n )
1189 {
1190 for( long nEmType=0; nEmType<2; ++nEmType )
1191 {
1192 LineSv[n].SumLine[nEmType+2] += (realnum) LineSv[n].SumLine[nEmType]*cumulative_factor;
1193 }
1194 }
1195 // save cumulative continua
1196 for( long n=0; n<rfield.nflux; ++n)
1197 {
1198
1199 rfield.flux[1][n] += (realnum) rfield.flux[0][n]*cumulative_factor;
1200 rfield.ConEmitReflec[1][n] += (realnum) rfield.ConEmitReflec[0][n]*cumulative_factor;
1201 rfield.ConEmitOut[1][n] += (realnum) rfield.ConEmitOut[0][n]*cumulative_factor;
1202 rfield.ConRefIncid[1][n] += (realnum) rfield.ConRefIncid[0][n]*cumulative_factor;
1203 rfield.flux_total_incident[1][n] += (realnum) rfield.flux_total_incident[0][n]*cumulative_factor;
1204 rfield.reflin[1][n] += (realnum) rfield.reflin[0][n]*cumulative_factor;
1205 rfield.outlin[1][n] += (realnum) rfield.outlin[0][n]*cumulative_factor;
1206 }
1207 }
1208
1209
1210 /* save previous iteration's results */
1211 struc.nzonePreviousIteration = nzone;
1212 for( long i=0; i<struc.nzonePreviousIteration; ++i )
1213 {
1214 struc.depth_last[i] = struc.depth[i];
1215 struc.drad_last[i] = struc.drad[i];
1216 }
1217
1218 /* all continue were attenuated in last zone in radius_increment to represent the intensity
1219 * in the middle of the next zone - this was too much since we now want
1220 * intensity emergent from outer edge of last zone */
1221 for( long i=0; i < rfield.nflux; i++ )
1222 {
1223 {
1224 enum{DEBUG_LOC=false};
1225 if( DEBUG_LOC)
1226 {
1227 fprintf(ioQQQ,"i=%li opac %.2e \n", i,
1228 (double)opac.opacity_abs[i]*radius.drad_x_fillfac/2.*(double)geometry.DirectionalCosin );
1229 }
1230 }
1231 double tau = opac.opacity_abs[i]*radius.drad_x_fillfac/2.*geometry.DirectionalCosin;
1232 ASSERT( tau > 0. );
1233 fac = sexp( tau );
1234
1235 /* >>chng 02 dec 14, add first test to see whether product of ratio is within
1236 * range of a float - ConInterOut can be finite and fac just above a small float,
1237 * so ratio exceeds largest size of a float */
1238 /*if( fac > SMALLFLOAT )*/
1239 if( (realnum)(fac/SDIV(rfield.ConInterOut[i]))>SMALLFLOAT && fac > SMALLFLOAT )
1240 {
1241 rfield.ConInterOut[i] /= (realnum)fac;
1242 rfield.outlin[0][i] /= (realnum)fac;
1243 rfield.outlin_noplot[i] /= (realnum)fac;
1244 }
1245 }
1246
1247 /* remember thickness of previous iteration */
1248 radius.StopThickness[iteration-1] = radius.depth;
1249 return;
1250}
t_atmdat atmdat
Definition atmdat.cpp:6
#define HS_NZ
Definition atmdat.h:125
void FeIIReset(void)
t_atoms atoms
Definition atoms.cpp:5
const int N_OI_LEVELS
Definition atoms.h:236
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define MALLOC(exp)
Definition cddefines.h:501
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 ipRecNetEsc
Definition cddefines.h:281
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
double e2(double t)
Definition service.cpp:299
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
const int ipRecRad
Definition cddefines.h:283
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
LinSv * LineSv
Definition cdinit.cpp:70
static t_fe2ovr_la & Inst()
Definition cddefines.h:175
double den
Definition mole.h:358
ProxyIterator< qStateProxy, qStateConstProxy > iterator
void zero_opacity()
t_co co
Definition co.cpp:5
t_colden colden
Definition colden.cpp:5
#define NCOLD
Definition colden.h:9
t_continuum continuum
Definition continuum.cpp:5
t_conv conv
Definition conv.cpp:5
void EdenChange(double EdenNew)
t_cosmology cosmology
Definition cosmology.cpp:11
#define CMB_TEMP
Definition cosmology.h:10
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dense dense
Definition dense.cpp:24
t_deuterium deut
Definition deuterium.cpp:8
void DynaIterEnd(void)
Definition dynamics.cpp:874
t_dynamics dynamics
Definition dynamics.cpp:44
void DynaIterStart(void)
t_geometry geometry
Definition geometry.cpp:5
void GrainRestartIter(void)
Definition grains.cpp:551
void GrainStartIter(void)
Definition grains.cpp:513
vector< diatomics * > diatoms
Definition h2.cpp:8
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_he he
Definition he.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
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
static double UV_Cont_rel2_Habing_TH85_face
static double H2_Solomon_dissoc_rate_used_H2g_save
static double UV_Cont_rel2_Draine_DB96_face
static realnum xIonFsave[LIMELM][LIMELM+1]
static double ortho_save
static double *** hnsav
static double hmitot_save
static double h2plus_heat_save
void IterRestart(void)
static realnum d5200r
static double HeatH2Dish_used_save
static realnum deutDenseSave1
static double H2_photodissoc_used_H2g_save
static realnum *** SaveMoleChTrRate
void IterEnd(void)
static double * den_save
static double UV_Cont_rel2_Draine_DB96_depth
static double para_save
static long int IonHighSave[LIMELM]
static double hmihet_save
static double H2_photodissoc_used_H2s_save
static double *** HOpacRatSav
static realnum supsav[LIMELM][LIMELM]
static double edsav
static double drSave
static double ** saveMoleSink
static double H2_H2g_to_H2s_rate_used_save
static double drNextSave
static double H2_Solomon_dissoc_rate_used_H2s_save
static long int IonLowSave[LIMELM]
static realnum p2nit
static double ** saveMoleSource
static realnum gas_phase_save[LIMELM]
static realnum deutDenseSave0
static double deriv_HeatH2Dexc_used_save
static bool lgHNSAV
void IterStart(void)
static double HeatH2Dexc_used_save
static double HeatSave[LIMELM][LIMELM]
t_iterations iterations
Definition iterations.cpp:5
t_LineSave LineSave
Definition lines.cpp:5
void Magnetic_reinit(void)
Definition magnetic.cpp:123
t_mean mean
Definition mean.cpp:17
void molcol(const char *chLabel, FILE *ioMEAN)
Definition molcol.cpp:12
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double TE1RYD
Definition physconst.h:183
t_pressure pressure
Definition pressure.cpp:5
void PresTotCurrent(void)
t_prt prt
Definition prt.cpp:10
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
void RT_OTS_Update(double *SumOTS)
Definition rt_ots.cpp:488
t_secondaries secondaries
t_StopCalc StopCalc
Definition stopcalc.cpp:5
t_struc struc
Definition struc.cpp:6
long int nSpecies
Definition taulines.cpp:21
species * dBaseSpecies
Definition taulines.cpp:14
t_thermal thermal
Definition thermal.cpp:5
t_timesc timesc
Definition timesc.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5