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