cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_ots.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 /*RT_OTS compute diffuse fields due to H, He atoms, ion, triplets, metal recombination,
4  * called by ConvBase */
5 /*RT_OTS_AddLine add local destruction of lines to ots field */
6 /*RT_OTS_AddCont add local destruction of continuum to ots field */
7 /*RT_OTS_Update sum flux, otscon, otslin, ConInterOut, outlin, to form SummeDif, SummedCon SummedOcc */
8 /*RT_OTS_Zero - zero out some vectors -
9  * this is only called when code initialized by ContSetIntensity */
10 /*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
11 #include "cddefines.h"
12 #include "taulines.h"
13 #include "opacity.h"
14 #include "dense.h"
15 #include "iso.h"
16 #include "hmi.h"
17 #include "h2.h"
18 #include "rfield.h"
19 #include "conv.h"
20 #include "rt.h"
21 #include "heavy.h"
22 #include "he.h"
23 #include "trace.h"
24 #include "mole.h"
25 #include "freebound.h"
26 #include "two_photon.h"
27 #include "cosmology.h"
28 
29 /* this flag may be used for debugging ots rate changes */
30 static int nOTS_Line_type = -1;
31 static int nOTS1=-1 , nOTS2=-1;
32 /*add local destruction of continuum to ots field */
34  /* the ots rate itself */
35  realnum ots,
36  /* pointer to continuum cell for ots, on f scale */
37  long int ip);
38 
39 /* =================================================================== */
40 
41 void RT_OTS(void)
42 {
43  long int
44  ipla,
45  ipISO ,
46  nelem,
47  n;
48  realnum
49  difflya,
50  esc,
51  ots;
52 
53  /* the Bowen HeII yield
54  * >>chng 06 aug 08, from 0.3 to 0.4, update with netzer */
55  const realnum BOWEN = 0.5f;
56  long int ipHi,
57  ipLo;
58 
59  double bwnfac;
60  double ots660;
61  realnum cont_phot_destroyed;
62  double save_lya_dest,
63  save_he2lya_dest;
64 
65  double save_he2rec_dest;
66 
67  /*static long int nCall=0;
68  ++nCall;
69  fprintf(ioQQQ,"debugggtos enter %li\n", nCall );*/
70 
71  DEBUG_ENTRY( "RT_OTS()" );
72 
73  for( long i=0; i < rfield.nflux; i++ )
74  {
75  rfield.otslin[i] = 0.;
76  rfield.otscon[i] = 0.;
77  }
78 
79  /**************************************************************************
80  *
81  * the bowen HeII - OIII fluorescence problem
82  *
83  **************************************************************************/
84  nOTS_Line_type = 0;
85  nelem = ipHELIUM;
86  if( dense.lgElmtOn[nelem] )
87  {
88  /* conversion per unit atom to OIII, at end of sub we divide by it,
89  * to fix lines back to proper esc/dest probs */
90  bwnfac = BOWEN * MAX2(0.f,1.f- iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Pesc_total());
91 
92  /* the last factor accounts for fact that two photons are produced,
93  * and the branching ratio */
94  ots660 = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Aul()*
95  iso_sp[ipH_LIKE][nelem].st[ipH2p].Pop()*
96  /*>>chng 06 aug 08, rm 0.8 factor from below, renorm aft discuss with Netzer */
97  iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Pdest() *BOWEN*2.0;
98 
99  /* now add this to the ots field */
100  if( ots660 > SMALLFLOAT )
101  RT_OTS_AddLine(ots660 , he.ip660 );
102 
103  /* decrease the destruction prob by the amount we will add elsewhere,
104  * ok since dest probs updated on every iteration*/
105  iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Pdest() *= (realnum)bwnfac;
106  ASSERT( iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().Pdest() >= 0. );
107  {
108  /* debugging code for line oscillation problems
109  * most often Lya OTS oscillations*/
110  enum {DEBUG_LOC=false};
111  if( DEBUG_LOC )
112  {
113  fprintf(ioQQQ,"DEBUG HeII Bowen nzone %li bwnfac:%.2e bwnfac esc:%.2e ots660 %.2e\n",
114  nzone,
115  bwnfac ,
116  bwnfac/BOWEN ,
117  ots660 );
118  }
119  }
120  }
121 
122  else
123  {
124  bwnfac = 1.;
125  }
126 
127  /* save Lya loss rates so we can reset at end */
128  save_lya_dest = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest();
129 
130  /* this is option to kill Lya ots rates,
131  * rfield.lgLyaOTS is usually true (1), and set false (0) with
132  * no lya ots command */
134 
135  /* option to kill heii lya and rec continua ots */
136  save_he2lya_dest = iso_sp[ipH_LIKE][ipHELIUM].trans(ipH2p,ipH1s).Emis().Pdest();
138 
139  /* option to kill heii lya and rec continua ots */
140  save_he2rec_dest = iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].RadRecomb[ipRecRad];
142 
143  nOTS_Line_type = 1;
144 
145  /* make ots fields due to lines and continua of species treated with unified
146  * isoelectronic sequence */
147  /* loop over all elements */
148  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
149  {
150  for( nelem=ipISO; nelem < LIMELM; nelem++ )
151  {
152  nOTS1 = ipISO;
153  nOTS2 = nelem;
154  /* do if this stage exists */
158  if( (dense.IonHigh[nelem] >= nelem+1-ipISO ) )
159  {
160  /* generate line ots rates */
161  /* now loop over all possible levels, but skip non-radiative
162  * since there is no pointer to this continuum */
163  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
164  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
165  {
166  for( ipLo=0; ipLo < ipHi; ipLo++ )
167  {
168  /* this signifies a fake line */
169  /* >>chng 03 may 19, DEST0 is the smallest possible
170  * dest prob - not a real number, don't add to ots field */
171  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() < 1 ||
172  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest()<= DEST0 )
173  continue;
174 
175  /* ots rates, the destp prob was set in hydropesc */
176  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots() =
177  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul()*
178  iso_sp[ipISO][nelem].st[ipHi].Pop()*
179  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest();
180 
181  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots() >= 0. );
182  /* way to kill lyalpha ots rates
183  if( nelem==ipHYDROGEN && ipHi==ipH2p && ipLo==ipH1s )
184  {
185  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots() = 0.;
186  } */
187 
188  /* finally dump the ots rate into the stack */
189  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots() > SMALLFLOAT )
190  RT_OTS_AddLine(iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots(),
191  iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() );
192  }
193  }
194  {
195  /* debugging code for line oscillation problems
196  * most often Lya OTS oscillations*/
197  /*@-redef@*/
198  enum {DEBUG_LOC=false};
199  /*@+redef@*/
200  if( DEBUG_LOC )
201  {
202  long ip;
203  if( ipISO==0 && nelem==0 && nzone>500 )
204  {
205  ipHi = 2;
206  ipLo = 0;
207  ip = iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont();
208  fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\tEdenTrue\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
209  fnzone,
210  dense.EdenTrue ,
211  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().ots(),
212  opac.opacity_abs[ip-1],
213  iso_sp[ipISO][nelem].st[ipHi].Pop(),
214  dense.xIonDense[nelem][nelem+1-ipISO],
215  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest(),
216  rfield.otslin[ip-1]);
217  }
218  }
219  }
220 
221  /**************************************************************************
222  *
223  * ots recombination bound-free b-f continua continuum
224  *
225  **************************************************************************/
226 
227  /* put in OTS continuum */
228  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
229  for( n=0; n < iso_sp[ipISO][nelem].numLevels_local; n++ )
230  {
231  cont_phot_destroyed = (realnum)(iso_sp[ipISO][nelem].fb[n].RadRecomb[ipRecRad]*
232  (1. - iso_sp[ipISO][nelem].fb[n].RadRecomb[ipRecEsc])*dense.eden*
233  dense.xIonDense[nelem][nelem+1-ipISO]);
234  ASSERT( cont_phot_destroyed >= 0. );
235 
236  /* continuum energy index used in this routine is decremented by one there */
237  RT_OTS_AddCont(cont_phot_destroyed,iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon);
238  /* debugging code for rec continua */
239  {
240  /*@-redef@*/
241  enum {DEBUG_LOC=false};
242  /*@+redef@*/
243  if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 )
244  {
245  long ip = iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1;
246  fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
247  fnzone,
248  n ,
249  iso_sp[ipISO][nelem].st[n].Pop(),
250  cont_phot_destroyed,
251  cont_phot_destroyed/opac.opacity_abs[ip],
252  rfield.otscon[ip] ,
253  opac.opacity_abs[ip] ,
254  findspecieslocal("H-")->den ,
256  }
257  }
258  }
259  }
260  }
261  }
262  /* more debugging code for rec continua */
263  {
264  enum {DEBUG_LOC=false};
265  if( DEBUG_LOC )
266  {
267  nelem = 0;
268  fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
269  nzone,
270  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[0].ipIsoLevNIonCon-1],
271  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[1].ipIsoLevNIonCon-1],
272  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[3].ipIsoLevNIonCon-1],
273  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[4].ipIsoLevNIonCon-1],
274  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[5].ipIsoLevNIonCon-1],
275  rfield.otscon[iso_sp[ipH_LIKE][nelem].fb[6].ipIsoLevNIonCon-1],
276  opac.opacity_abs[iso_sp[ipH_LIKE][nelem].fb[6].ipIsoLevNIonCon-1]);
277  }
278  }
279 
280  /* now reset Lya dest prob in case is was clobbered by rfield.lgHeIIOTS */
281  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest() = (realnum)save_lya_dest;
282  iso_sp[ipH_LIKE][ipHELIUM].trans(ipH2p,ipH1s).Emis().Pdest() = (realnum)save_he2lya_dest;
283  iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].RadRecomb[ipRecRad] = save_he2rec_dest;
284 
285  nelem = ipHELIUM;
286  if( dense.lgElmtOn[nelem] && bwnfac > 0. )
287  {
288  /* increase the destruction prob by the amount we decreased it above */
290  }
291 
292  if( trace.lgTrace )
293  {
294  fprintf(ioQQQ," RT_OTS Pdest %.2e ots rate %.2e in otslin %.2e con opac %.2e\n",
295  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
296  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().ots() ,
299  );
300  }
301 
302  nOTS_Line_type = 2;
303  /* recombination Lya for all elements not yet converted into std isoelectronc form */
304  for( nelem=NISO; nelem < LIMELM; nelem++ )
305  {
306  long int ion;
307  /* do not include species treated in iso-electronic fashion in the following,
308  * these were treated above */
309  for( ion=0; ion < nelem+1-NISO; ion++ )
310  {
311  if( dense.xIonDense[nelem][ion+1] > 0. )
312  {
313  nOTS1 = nelem;
314  nOTS2 = ion;
315  /* now do the recombination Lya */
316  ipla = Heavy.ipLyHeavy[nelem][ion];
317  ASSERT( ipla>0 );
318  esc = opac.ExpmTau[ipla-1];
319  /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
320  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
321  /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
322  * setting the ots rates to zero */
323  ots = difflya*MAX2(0.f,1.f-esc);
324  /*if( nelem==6 && ion==2 )
325  fprintf(ioQQQ," debugggnly\t %.2f\t%.2e\n",fnzone, ots );*/
326  ASSERT( ots >= 0.);
327  /*if( iteration == 2 && nzone>290 && ipla== 2339 )
328  fprintf(ioQQQ,"recdebugg1 %.2e %li %li %.2e %.2e \n",
329  ots, nelem, ion,
330  esc , dense.xIonDense[nelem][ion+1]);*/
331  if( ots > SMALLFLOAT )
332  RT_OTS_AddLine(ots,ipla);
333 
334  /* now do the recombination balmer lines */
335  ipla = Heavy.ipBalHeavy[nelem][ion];
336  esc = opac.ExpmTau[ipla-1];
337  /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
338  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
339  /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
340  * setting the ots rates to zero */
341  ots = difflya*MAX2(0.f,1.f-esc);
342  ASSERT( ots >= 0.);
343  /*if( iteration == 2 &&nzone==294 && ipla== 2339 )
344  fprintf(ioQQQ,"recdebugg2 %.2e %li %li\n", ots, nelem, ion );*/
345  if( ots > SMALLFLOAT )
346  RT_OTS_AddLine(ots,ipla);
347  }
348  }
349  }
350 
351  nOTS_Line_type = 5;
352  /* do the level2 level 2 lines */
353  for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ )
354  {
355  if( (*TauLine2[nOTS1].Hi()).IonStg() < (*TauLine2[nOTS1].Hi()).nelem()+1-NISO )
356  {
357  TauLine2[nOTS1].Emis().ots() = (*TauLine2[nOTS1].Hi()).Pop() * TauLine2[nOTS1].Emis().Aul() * TauLine2[nOTS1].Emis().Pdest();
358  if( TauLine2[nOTS1].Emis().ots() > SMALLFLOAT )
359  RT_OTS_AddLine( TauLine2[nOTS1].Emis().ots() , TauLine2[nOTS1].ipCont());
360  }
361  }
362 
363  nOTS_Line_type = 6;
364  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
365  {
366  if( dBaseSpecies[ipSpecies].lgActive )
367  {
368  for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
369  tr != dBaseTrans[ipSpecies].end(); ++tr)
370  {
371  int ipHi = (*tr).ipHi();
372  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
373  continue;
374  (*tr).Emis().ots() = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() * (*tr).Emis().Pdest();
375  RT_OTS_AddLine( (*tr).Emis().ots() , (*tr).ipCont());
376  }
377  }
378  }
379 
380  nOTS_Line_type = 7;
381  /* the large H2 molecule */
382  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
383  (*diatom)->H2_RT_OTS();
384 
385  return;
386 }
387 
388 /* =================================================================== */
389 
390 void RT_OTS_AddLine(double ots,
391  /* pointer on the f scale */
392  long int ip )
393 {
394 
395  DEBUG_ENTRY( "RT_OTS_AddLine()" );
396 
397  /* add ots due to line destruction to radiation field */
398 
399  /* return if outside bounds of this continuum source, ip > rfield.nflux
400  * first case ip==0 happens when called with dummy line */
401  if( ip==0 || ip > rfield.nflux )
402  {
403  return;
404  }
405 
406  /*the local ots rate must be non-negative */
407  ASSERT( ots >= 0. );
408  /* continuum pointer must be positive */
409  ASSERT( ip > 0 );
410 
411  /* add locally destroyed flux of photons to line OTS array */
412  /* check whether local gas opacity (units cm-1) is positive, if so
413  * convert line destruction rate into ots rate by dividing by it */
414  if( opac.opacity_abs[ip-1] > 0. )
415  {
416  rfield.otslin[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
417  }
418  /* first iteration is 1, second is two */
419  {
420  enum {DEBUG_LOC=false};
421  if( DEBUG_LOC && /*iteration == 2 && nzone>294 &&*/ ip== 2363 /*&& ots > 1e16*/)
422  {
423  fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t",
424  ots ,
425  opac.opacity_abs[ip-1],
426  ots/opac.opacity_abs[ip-1] );
427  fprintf(ioQQQ,"iteration %li type %i %i %i \n",
428  iteration,
430  nOTS1,nOTS2 );
431  }
432  }
433  return;
434 }
435 
436 /* =================================================================== */
437 
438 /*add local destruction of continuum to ots field */
440  /* the ots rate itself */
441  realnum ots,
442  /* pointer to continuum cell for ots, on f scale */
443  long int ip)
444 {
445 
446  DEBUG_ENTRY( "RT_OTS_AddCont()" );
447 
448  /*
449  * routine called to add ots due to continuum destruction to
450  * radiation field
451  */
452 
453  /* check if outside bounds of this continuum source */
454  if( ip > rfield.nflux )
455  {
456  return;
457  }
458 
459  ASSERT( ip > 0 );
460  ASSERT( ots >= 0. );
461  ASSERT( ip <= rfield.nflux_with_check );
462 
463  /* add locally destroyed flux of photons to continuum OTS array */
464  /* check whether local gas opacity (units cm-1) is positive, if so
465  * convert continuum destruction rate into ots rate by dividing by it */
466  if( opac.opacity_abs[ip-1] > 0. )
467  {
468  rfield.otscon[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
469  }
470  return;
471 }
472 
473 /* =================================================================== */
474 
475 /*RT_OTS_Update update ots rates, called in ConvBase */
476 void RT_OTS_Update(double *SumOTS) /* summed ots rates */
477 {
478  long int i;
479 
480  DEBUG_ENTRY( "RT_OTS_Update()" );
481 
482  /* option to kill ots rates with no ots lines command */
483  if( rfield.lgKillOTSLine )
484  {
485  for( i=0; i < rfield.nflux; i++ )
486  {
487  rfield.otslin[i] = 0.;
488  }
489  }
490 
491  memset(rfield.ConOTS_local_photons , 0 , (unsigned)rfield.nflux_with_check*sizeof(realnum) );
492  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
493  {
494  /* >>chng 01 sep 23, rewrote for iso sequences */
495  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
496  {
497  if( dense.IonHigh[nelem] >= nelem+1-ipISO )
498  {
499  t_iso_sp* sp = &iso_sp[ipISO][nelem];
500 
501  for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
502  {
504  for( long nu=0; nu < tnu->ipTwoPhoE; nu++ )
505  {
506  rfield.ConOTS_local_photons[nu] += tnu->local_emis[nu] * (1.f - opac.ExpmTau[nu]);
507  }
508  }
509  }
510  }
511  }
512 
513  /* remember largest change in ots rates */
514  *SumOTS = 0.;
515  /* now update new ots rates */
516  for( i=0; i < rfield.nflux; ++i )
517  {
518  double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );
519 
520  /* this is local ots continuum created by destroyed diffuse continua,
521  * currently only two-photon */
522  rfield.ConOTS_local_OTS_rate[i] = (realnum)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);
523 
524  /* remember sum of ots rates for convergence criteria */
525  *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
526 
527  if( !cosmology.lgDo )
528  {
532  }
533 
534  rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
536  }
537 
539 
540  /* sum of accumulated flux from particular frequency to infinity */
541  rfield.flux_accum[rfield.nflux-1] = 0.;
542  for( i=1; i < rfield.nflux; i++ )
543  {
546  }
547 
548  /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
549  * between plasma frequency and energy where brems is thin */
550  ASSERT( rfield.ipPlasma > 0 );
551 
552  /* all radiation fields are zero below plasma frequency */
553  for( i=0; i < rfield.ipPlasma-1; i++ )
554  {
555  rfield.otscon[i] = 0.;
558  rfield.otslin[i] = 0.;
559  rfield.SummedDif[i] = 0.;
560  rfield.SummedCon[i] = 0.;
561  rfield.SummedOcc[i] = 0.;
562  rfield.ConInterOut[i] = 0.;
563  }
564  return;
565 }
566 
567 /* =================================================================== */
568 
569 /*RT_OTS_Zero zero out some vectors - this is only called when code
570  * initialized by ContSetIntensity */
571 void RT_OTS_Zero( void )
572 {
573  long int i;
574 
575  DEBUG_ENTRY( "RT_OTS_Zero()" );
576 
577  /* this loop goes up to nflux itself (<=) since the highest cell
578  * will be used to pass unity through the code to verify integrations */
579  for( i=0; i <= rfield.nflux; i++ )
580  {
581  rfield.SummedDif[i] = 0.;
582  /* the main ots vectors */
583  rfield.otscon[i] = 0.;
584  rfield.otslin[i] = 0.;
585 
586  rfield.ConInterOut[i] = 0.;
587  rfield.outlin[0][i] = 0.;
588  rfield.outlin_noplot[i] = 0.;
589  rfield.SummedDif[i] = 0.;
590  /* "zero" for the summed con will be just the incident radiation field */
591  rfield.SummedCon[i] = rfield.flux[0][i];
595  }
596 
598 
600  return;
601 }
602 
603 /* =================================================================== */
604 
605 /*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
606 void RT_OTS_ChkSum(long int ipPnt)
607 {
608  static long int nInsane=0;
609  long int i;
610  double phisig;
611  const int LIM_INSAME_PRT = 30;
612 
613  DEBUG_ENTRY( "RT_OTS_ChkSum()" );
614 
615  /* this entire sub is a sanity check */
616  /* >>chng 02 jul 23, lower bound from 0 to rfield.ipEnergyBremsThin - since now
617  * set radiation field to black body below this energy */
618  for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ )
619  {
620  phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly +
621  rfield.outlin[0][i]+
624  /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero whild
625  * phisig is just above small float */
626  if( phisig > 0. && rfield.SummedDif[i]> 0.)
627  {
628  if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 )
629  {
630  ++nInsane;
631  /* limit amount of printout - in many failures would print entire
632  * continuum array */
633  if( nInsane < LIM_INSAME_PRT )
634  {
635  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n",
636  rfield.anu(i), rfield.SummedDif[i]/phisig - 1., i );
637  fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n",
638  rfield.SummedDif[i], phisig );
639  fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n",
641  rfield.outlin[0][i]+rfield.outlin_noplot[i] );
642  fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n",
643  rfield.chLineLabel[i].c_str(), rfield.chContLabel[i].c_str() );
644  }
645  }
646  }
647 
648  phisig += rfield.flux[0][i];
649  /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero when
650  * phisig is just above small float */
651  if( phisig > 0. && rfield.SummedDif[i]> 0. )
652  {
653  if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 )
654  {
655  ++nInsane;
656  /* limit amount of printout - in many failures would print entire
657  * continuum array */
658  if( nInsane < LIM_INSAME_PRT )
659  {
660  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum %3ld, insane SummedCon at energy %.5e error=%.2e i=%ld\n",
661  ipPnt, rfield.anu(i), rfield.SummedCon[i]/phisig - 1., i );
662  fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n",
663  rfield.SummedCon[i], phisig );
664  fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n",
666  rfield.outlin[0][i]+rfield.outlin_noplot[i], rfield.flux[0][i] );
667  fprintf( ioQQQ, " line continuum here are %s %s\n",
668  rfield.chLineLabel[i].c_str(), rfield.chContLabel[i].c_str()
669  );
670  }
671  }
672  }
673  }
674 
675  if( nInsane > 0 )
676  {
677  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n");
678  /* TotalInsanity exits after announcing a problem */
679  TotalInsanity();
680  }
681  return;
682 }
683 
684 /* =================================================================== */
685 
686 /*RT_OTS_PrtRate print continuum and line ots rates when trace ots is on */
688  /* weakest rate to print */
689  double weak ,
690  /* flag, 'c' continuum, 'l' line, 'b' both */
691  int chFlag )
692 {
693  long int i;
694 
695  DEBUG_ENTRY( "RT_OTS_PrtRate()" );
696 
697  /* arg must be one of these three */
698  ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' );
699 
700  /*
701  * both printouts have cell number (on C array scale)
702  * energy in ryd
703  * the actual value of the ots rate
704  * the ots rate relative to the continuum at that energy
705  * rate times opacity
706  * all are only printed if greater than weak
707  */
708 
709  /*===================================================================*/
710  /* first print ots continua */
711  /*===================================================================*/
712  if( chFlag == 'c' || chFlag == 'b' )
713  {
714  fprintf( ioQQQ, " DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n",
715  weak,fnzone ,TorF(conv.lgConvIoniz()) );
716 
717  for( i=0; i < rfield.nflux_with_check; i++ )
718  {
719  if( rfield.otscon[i]*opac.opacity_abs[i] > weak )
720  {
721  fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s \n",
722  i,
723  rfield.anu(i),
724  rfield.otscon[i],
725  opac.opacity_abs[i],
726  rfield.otscon[i]*opac.opacity_abs[i],
727  rfield.chContLabel[i].c_str());
728 
729  }
730  }
731  }
732 
733  /*===================================================================*/
734  /* second print ots line rates */
735  /*===================================================================*/
736  if( chFlag == 'l' || chFlag == 'b' )
737  {
738  fprintf( ioQQQ, "DEBUG density He %.2e He+2 %.2e O+2 %.2e\n",
740  dense.xIonDense[ipOXYGEN][2] );
741  fprintf( ioQQQ, " DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n",
742  weak,fnzone,TorF(conv.lgConvIoniz()) );
743 
744  for( i=0; i < rfield.nflux_with_check; i++ )
745  {
746  if( rfield.otslin[i]*opac.opacity_abs[i] > weak )
747  {
748  fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n",
749  i,
750  rfield.anu(i),
751  rfield.otslin[i],
752  opac.opacity_abs[i],
754  rfield.chLineLabel[i].c_str() ,
755  rfield.line_count[i] );
756  }
757  }
758  }
759  return;
760 }
void RT_OTS_ChkSum(long int ipPnt)
Definition: rt_ots.cpp:606
#define DEST0
Definition: rt.h:154
void RT_OTS_Update(double *SumOTS)
Definition: rt_ots.cpp:476
long int * line_count
Definition: rfield.h:61
double & ots() const
Definition: emission.h:678
realnum * ConOTS_local_OTS_rate
Definition: rfield.h:170
double HMinus_photo_rate
Definition: hmi.h:53
double * opacity_abs
Definition: opacity.h:103
qList st
Definition: iso.h:482
long int ipEnergyBremsThin
Definition: rfield.h:228
void setTrimming()
Definition: rfield.cpp:92
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
STATIC void RT_OTS_AddCont(realnum ots, long int ip)
Definition: rt_ots.cpp:439
t_opac opac
Definition: opacity.cpp:5
realnum ** flux
Definition: rfield.h:70
t_Heavy Heavy
Definition: heavy.cpp:5
vector< string > chContLabel
Definition: rfield.h:215
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
realnum * outlin_noplot
Definition: rfield.h:191
long int IonHigh[LIMELM+1]
Definition: dense.h:130
char TorF(bool l)
Definition: cddefines.h:753
void RT_OTS_PrtRate(double weak, int chFlag)
Definition: rt_ots.cpp:687
const int ipOXYGEN
Definition: cddefines.h:355
const double SMALLDOUBLE
Definition: cpu.h:250
realnum xLyaHeavy[LIMELM][LIMELM]
Definition: heavy.h:21
t_conv conv
Definition: conv.cpp:5
static int nOTS2
Definition: rt_ots.cpp:31
double * SummedCon
Definition: rfield.h:163
realnum * SummedOcc
Definition: rfield.h:165
static int nOTS1
Definition: rt_ots.cpp:31
bool lgKillOTSLine
Definition: rfield.h:421
realnum ** outlin
Definition: rfield.h:191
FILE * ioQQQ
Definition: cddefines.cpp:7
static int nOTS_Line_type
Definition: rt_ots.cpp:30
molezone * findspecieslocal(const char buf[])
long int nzone
Definition: cddefines.cpp:14
bool lgConvIoniz() const
Definition: conv.h:108
realnum * ConOTS_local_photons
Definition: rfield.h:170
vector< freeBound > fb
Definition: iso.h:481
TransitionList TauLine2("TauLine2",&AnonStates)
bool lgDo
Definition: cosmology.h:44
double anu(size_t i) const
Definition: mesh.h:111
long int nSpecies
Definition: taulines.cpp:22
realnum * flux_accum
Definition: rfield.h:79
bool lgHeIIOTS
Definition: rfield.h:412
t_dense dense
Definition: global.cpp:15
void resetCoarseTransCoef()
Definition: rfield.h:485
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nflux_with_check
Definition: rfield.h:51
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
void RT_OTS_Zero(void)
Definition: rt_ots.cpp:571
bool lgOutOnly
Definition: rfield.h:224
vector< two_photon > TwoNu
Definition: iso.h:598
const int ipH1s
Definition: iso.h:29
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
long int ipPlasma
Definition: rfield.h:434
EmissionList::reference Emis() const
Definition: transition.h:447
bool lgLyaOTS
Definition: rfield.h:408
double EdenTrue
Definition: dense.h:232
t_rfield rfield
Definition: rfield.cpp:9
long & ipCont() const
Definition: transition.h:489
realnum * convoc
Definition: rfield.h:115
realnum * ConInterOut
Definition: rfield.h:156
float realnum
Definition: cddefines.h:124
vector< diatomics * > diatoms
Definition: h2.cpp:8
realnum * otscon
Definition: rfield.h:185
const int ipRecRad
Definition: cddefines.h:332
void RT_OTS_AddLine(double ots, long int ip)
Definition: rt_ots.cpp:390
const int ipRecEsc
Definition: cddefines.h:328
long nWindLine
Definition: cdinit.cpp:19
vector< string > chLineLabel
Definition: rfield.h:212
bool lgElmtOn[LIMELM]
Definition: dense.h:160
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
Definition: iso.h:470
bool lgInd2nu_On
Definition: iso.h:375
realnum gas_phase[LIMELM]
Definition: dense.h:76
const int ipH2p
Definition: iso.h:31
long int ipLyHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:11
realnum & Pdest() const
Definition: emission.h:588
#define ASSERT(exp)
Definition: cddefines.h:617
long int ipBalHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:11
long int ip660
Definition: he.h:21
t_he he
Definition: he.cpp:5
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
t_cosmology cosmology
Definition: cosmology.cpp:8
const int ipHELIUM
Definition: cddefines.h:349
vector< species > dBaseSpecies
Definition: taulines.cpp:15
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:828
realnum * ExpmTau
Definition: opacity.h:144
bool lgInducProcess
Definition: rfield.h:235
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
realnum * SummedDif
Definition: rfield.h:164
t_hmi hmi
Definition: hmi.cpp:5
double fnzone
Definition: cddefines.cpp:15
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
void CalcTwoPhotonEmission(two_photon &tnu, bool lgDoInduced)
Definition: two_photon.cpp:125
const int ipHYDROGEN
Definition: cddefines.h:348
long int nflux
Definition: rfield.h:48
realnum & Aul() const
Definition: emission.h:668
void RT_OTS(void)
Definition: rt_ots.cpp:41
long int numLevels_local
Definition: iso.h:529
EmissionList & Emis()
Definition: transition.h:363
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13