cloudy  trunk
 All Data Structures 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-2008 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 "mole_co_atom.h"
22 #include "atomfeii.h"
23 #include "heavy.h"
24 #include "he.h"
25 #include "trace.h"
26 
27 /* this flag may be used for debugging ots rate changes */
28 static int nOTS_Line_type = -1;
29 static int nOTS1=-1 , nOTS2=-1;
30 /*add local destruction of continuum to ots field */
32  /* the ots rate itself */
33  realnum ots,
34  /* pointer to continuum cell for ots, on f scale */
35  long int ip);
36 
37 /* =================================================================== */
38 
39 void RT_OTS(void)
40 {
41  long int
42  ipla,
43  ipISO ,
44  nelem,
45  n;
46  realnum
47  difflya,
48  esc,
49  ots;
50 
51  /* the Bowen HeII yield
52  * >>chng 06 aug 08, from 0.3 to 0.4, update with netzer */
53 # define BOWEN 0.5f
54  long int ipHi,
55  ipLo;
56 
57  double bwnfac;
58  double ots660;
59  realnum cont_phot_destroyed;
60  double save_lya_dest,
61  save_he2lya_dest;
62 
63  double save_he2rec_dest;
64 
65  /*static long int nCall=0;
66  ++nCall;
67  fprintf(ioQQQ,"debugggtos enter %li\n", nCall );*/
68 
69  DEBUG_ENTRY( "RT_OTS()" );
70 
71  /**************************************************************************
72  *
73  * the bowen HeII - OIII fluorescense problem
74  *
75  **************************************************************************/
76  nOTS_Line_type = 0;
77  nelem = ipHELIUM;
78  if( dense.lgElmtOn[nelem] )
79  {
80  /* conversion per unit atom to OIII, at end of sub we divide by it,
81  * to fix lines back to proper esc/dest probs */
82  bwnfac = BOWEN * MAX2(0.f,1.f- Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pesc -
83  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pelec_esc );
84 
85  /* the last factor accounts for fact that two photons are produced,
86  * and the branching ratio */
87  ots660 = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul*
88  StatesElem[ipH_LIKE][nelem][ipH2p].Pop*dense.xIonDense[nelem][nelem+1]*
89  /*>>chng 06 aug 08, rm 0.8 factor from below, renorm aft discuss with Netzer */
90  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *BOWEN*2.0;
91 
92  /* now add this to the ots field */
93  if( ots660 > SMALLFLOAT )
94  RT_OTS_AddLine(ots660 , he.ip660 );
95 
96  /* decrease the destruction prob by the amount we will add elsewhere,
97  * ok since dest probs updated on every iteration*/
98  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *= (realnum)bwnfac;
99  ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest >= 0. );
100  {
101  /* debugging code for line oscillation problems
102  * most often Lya OTS oscillations*/
103  /*@-redef@*/
104  enum {DEBUG_LOC=false};
105  /*@+redef@*/
106  if( DEBUG_LOC )
107  {
108  fprintf(ioQQQ,"DEBUG HeII Bowen nzone %li bwnfac:%.2e bwnfac esc:%.2e ots660 %.2e\n",
109  nzone,
110  bwnfac ,
111  bwnfac/BOWEN ,
112  ots660 );
113  }
114  }
115  }
116 
117  else
118  {
119  bwnfac = 1.;
120  }
121 
122  /* save Lya loss rates so we can reset at end */
123  save_lya_dest = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest;
124 
125  /* >>chng 03 may 28, hydro.dstfe2lya should not go into the ots field since it was
126  * actually lost in exciting FeII
127  * this is not ready to be used - need to recover ots with small feii atom
128  fprintf(ioQQQ,"feii debug\t%.3e\t%.3e\n",
129  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , hydro.dstfe2lya);
130  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = MAX2(0.f,
131  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest-hydro.dstfe2lya);*/
132 
133  /* this is option to kill Lya ots rates,
134  * rfield.lgLyaOTS is usually true (1), and set false (0) with
135  * no lya ots command */
137 
138  /* option to kill heii lya and rec continua ots */
139  save_he2lya_dest = Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest;
141 
142  /* option to kill heii lya and rec continua ots */
143  save_he2rec_dest = iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad];
145 
146  nOTS_Line_type = 1;
147 
148  /* make ots fields due to lines and continua of species treated with unified
149  * isoelectronic sequence */
150  /* loop over all elements */
151  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
152  {
153  for( nelem=ipISO; nelem < LIMELM; nelem++ )
154  {
155  nOTS1 = ipISO;
156  nOTS2 = nelem;
157  /* do if this stage exists */
161  if( (dense.IonHigh[nelem] >= nelem+1-ipISO ) )
162  {
163  /* generate line ots rates */
164  /* now loop over all possible levels, but cannot include two photon
165  * since there is no pointer to this continuum */
166  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
167  for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
168  {
169  for( ipLo=0; ipLo < ipHi; ipLo++ )
170  {
171  /* this signifies a fake line */
172  /* >>chng 03 may 19, DEST0 is the smallest possible
173  * dest prob - not a real number, don't add to ots field */
174  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 1 ||
175  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest<= DEST0 )
176  continue;
177 
178  /* ots rates, the destp prob was set in hydropesc */
179  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots =
180  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul*
181  StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*
182  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest;
183 
184  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots >= 0. );
185  /* way to kill lyalpha ots rates
186  if( nelem==ipHYDROGEN && ipHi==ipH2p && ipLo==ipH1s )
187  {
188  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots = 0.;
189  } */
190 
191  /* finally dump the ots rate into the stack */
192  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots > SMALLFLOAT )
193  RT_OTS_AddLine(Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots,
194  Transitions[ipISO][nelem][ipHi][ipLo].ipCont );
195  }
196  }
197  {
198  /* debugging code for line oscillation problems
199  * most often Lya OTS oscillations*/
200  /*@-redef@*/
201  enum {DEBUG_LOC=false};
202  /*@+redef@*/
203  if( DEBUG_LOC )
204  {
205  long ip;
206  if( ipISO==0 && nelem==0 && nzone>500 )
207  {
208  ipHi = 2;
209  ipLo = 0;
210  ip = Transitions[ipISO][nelem][ipHi][ipLo].ipCont;
211  fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\tEdenTrue\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
212  fnzone,
213  dense.EdenTrue ,
214  Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots,
215  opac.opacity_abs[ip-1],
216  StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO],
217  StatesElem[ipISO][nelem][ipHi].Pop,
218  dense.xIonDense[nelem][nelem+1-ipISO],
219  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest,
220  rfield.otslinNew[ip-1],
221  rfield.otslin[ip-1]);
222  }
223  }
224  }
225 
226  /**************************************************************************
227  *
228  * ots recombination bound-free b-f continua continuum
229  *
230  **************************************************************************/
231 
232  /* put in OTS continuum */
233  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
234  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
235  {
236  cont_phot_destroyed = (realnum)(iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
237  (1. - iso.RadRecomb[ipISO][nelem][n][ipRecEsc])*dense.eden*
238  dense.xIonDense[nelem][nelem+1-ipISO]);
239  ASSERT( cont_phot_destroyed >= 0. );
240 
241  /* continuum energy index used in this routine is decremented by one there */
242  RT_OTS_AddCont(cont_phot_destroyed,iso.ipIsoLevNIonCon[ipISO][nelem][n]);
243  /* debugging code for rec continua */
244  {
245  /*@-redef@*/
246  enum {DEBUG_LOC=false};
247  /*@+redef@*/
248  if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 )
249  {
250  long ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;
251  fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
252  fnzone,
253  n ,
254  StatesElem[ipISO][ipHYDROGEN][2].Pop*dense.xIonDense[ipHYDROGEN][1],
255  cont_phot_destroyed,
256  cont_phot_destroyed/opac.opacity_abs[ip],
257  rfield.otsconNew[ip] ,
258  rfield.otscon[ip] ,
259  opac.opacity_abs[ip] ,
260  hmi.Hmolec[ipMHm] ,
262  }
263  }
264  }
265  }
266  }
267  }
268  /* more debugging code for rec continua */
269  {
270  /*@-redef@*/
271  enum {DEBUG_LOC=false};
272  /*@+redef@*/
273  if( DEBUG_LOC )
274  {
275  nelem = 0;
276  fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
277  nzone,
285  }
286  }
287 
288  /* now reset Lya dest prob in case is was clobbered */
289  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest = (realnum)save_lya_dest;
290  Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest = (realnum)save_he2lya_dest;
291  iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] = save_he2rec_dest;
292 
293  nelem = ipHELIUM;
294  if( dense.lgElmtOn[nelem] && bwnfac > 0. )
295  {
296  /* increase the destruction prob by the amount we decreased it above */
297  Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest /= (realnum)bwnfac;
298  }
299 
300  if( trace.lgTrace )
301  {
302  fprintf(ioQQQ," RT_OTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n",
303  Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest , Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->ots ,
304  rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ,
306  );
307  }
308 
309  nOTS_Line_type = 2;
310  /* recombination Lya for all elements not yet converted into std isoelectronc form */
311  for( nelem=NISO; nelem < LIMELM; nelem++ )
312  {
313  long int ion;
314  /* do not include species treated in iso-electronic fashon in the following,
315  * these were treated above */
316  for( ion=0; ion < nelem+1-NISO; ion++ )
317  {
318  if( dense.xIonDense[nelem][ion+1] > 0. )
319  {
320  nOTS1 = nelem;
321  nOTS2 = ion;
322  /* now do the recombination Lya */
323  ipla = Heavy.ipLyHeavy[nelem][ion];
324  ASSERT( ipla>0 );
325  esc = opac.ExpmTau[ipla-1];
326  /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
327  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
328  /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
329  * setting the ots rates to zero */
330  ots = difflya*MAX2(0.f,1.f-esc);
331  /*if( nelem==6 && ion==2 )
332  fprintf(ioQQQ," debugggnly\t %.2f\t%.2e\n",fnzone, ots );*/
333  ASSERT( ots >= 0.);
334  /*if( iteration == 2 && nzone>290 && ipla== 2339 )
335  fprintf(ioQQQ,"recdebugg1 %.2e %li %li %.2e %.2e \n",
336  ots, nelem, ion,
337  esc , dense.xIonDense[nelem][ion+1]);*/
338  if( ots > SMALLFLOAT )
339  RT_OTS_AddLine(ots,ipla);
340 
341  /* now do the recombination balmer lines */
342  ipla = Heavy.ipBalHeavy[nelem][ion];
343  esc = opac.ExpmTau[ipla-1];
344  /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
345  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
346  /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
347  * setting the ots rates to zero */
348  ots = difflya*MAX2(0.f,1.f-esc);
349  ASSERT( ots >= 0.);
350  /*if( iteration == 2 &&nzone==294 && ipla== 2339 )
351  fprintf(ioQQQ,"recdebugg2 %.2e %li %li\n", ots, nelem, ion );*/
352  if( ots > SMALLFLOAT )
353  RT_OTS_AddLine(ots,ipla);
354  }
355  }
356  }
357 
358  nOTS_Line_type = 3;
359  /* do OTS and outward parts of FeII lines, if large atom is turned on */
360  FeII_OTS();
361 
362  nOTS_Line_type = 4;
363  /* do the main set of level1 lines */
364 # if 0
365  {
366 # include "lines_service.h"
367  if( nzone>290 ) DumpLine( &TauLines[74] );
368  }
369 # endif
370 /*# define FRACNEW 1.0f
371 # undef FRACNEW*/
372  for( nOTS1=1; nOTS1 < nLevel1; nOTS1++ )
373  {
374  /*realnum otsold = TauLines[nOTS1].ots;*/
376  /* * TauLines[nOTS1].ColOvTot;*/
377  /*if( nzone )
378  TauLines[nOTS1].ots = TauLines[nOTS1].ots*FRACNEW + otsold*(1.f-FRACNEW);*/
379  if( TauLines[nOTS1].Emis->ots > SMALLFLOAT )
380  RT_OTS_AddLine( TauLines[nOTS1].Emis->ots , TauLines[nOTS1].ipCont);
381  }
382 
383  nOTS_Line_type = 5;
384  /* do the level2 level 2 lines */
385  for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ )
386  {
387  if( TauLine2[nOTS1].Hi->IonStg < TauLine2[nOTS1].Hi->nelem+1-NISO )
388  {
390  if( TauLine2[nOTS1].Emis->ots > SMALLFLOAT )
391  RT_OTS_AddLine( TauLine2[nOTS1].Emis->ots , TauLine2[nOTS1].ipCont);
392  }
393  }
394  /*fprintf(ioQQQ,"debuggne1 \t%.2e\t%.2e\t%.2e\t%.2e\n",
395  TauLine2[743].ots/SDIV(opac.opacity_abs[743]),
396  TauLine2[744].ots/SDIV(opac.opacity_abs[744]) ,
397  TauLine2[743].ots,opac.opacity_abs[743]);*/
398  /*CO_OTS - add CO lines to ots fields, called by RT_OTS */
399 
400  nOTS_Line_type = 6;
401  CO_OTS();
402 
403  /* the large H2 molecule */
404  nOTS_Line_type = 7;
405  H2_RT_OTS();
406  /*fprintf(ioQQQ,"DEBUG ggtos exit\n");*/
407  return;
408 }
409 
410 /* =================================================================== */
411 
412 void RT_OTS_AddLine(double ots,
413  /* pointer on the f scale */
414  long int ip )
415 {
416 
417  DEBUG_ENTRY( "RT_OTS_AddLine()" );
418 
419  /* add ots due to line destruction to radiation field */
420 
421  /* return if outside bounds of this continuum source, ip > rfield.nflux
422  * first case ip==0 happens when called with dummy line */
423  if( ip==0 || ip > rfield.nflux )
424  {
425  return;
426  }
427 
428  /*the local ots rate must be non-negative */
429  ASSERT( ots >= 0. );
430  /* continuum pointer must be positive */
431  ASSERT( ip > 0 );
432 
433  /* add locally destroyed flux of photons to line OTS array */
434  /* check whether local gas opacity (units cm-1) is positive, if so
435  * convert line destruction rate into ots rate by dividing by it */
436  if( opac.opacity_abs[ip-1] > 0. )
437  {
438  rfield.otslinNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
439  }
440  /* first iteration is 1, second is two */
441  {
442  /*@-redef@*/
443  enum {DEBUG_LOC=false};
444  /*@+redef@*/
445  if( DEBUG_LOC && /*iteration == 2 && nzone>294 &&*/ ip== 2363 /*&& ots > 1e16*/)
446  {
447  fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t",
448  ots ,
449  opac.opacity_abs[ip-1],
450  ots/opac.opacity_abs[ip-1] );
451  fprintf(ioQQQ,"iteration %li type %i %i %i \n",
452  iteration,
454  nOTS1,nOTS2 );
455  }
456  }
457  return;
458 }
459 
460 /* =================================================================== */
461 
462 /*add local destruction of continuum to ots field */
464  /* the ots rate itself */
465  realnum ots,
466  /* pointer to continuum cell for ots, on f scale */
467  long int ip)
468 {
469 
470  DEBUG_ENTRY( "RT_OTS_AddCont()" );
471 
472  /*
473  * routine called to add ots due to continuum destruction to
474  * radiation field
475  */
476 
477  /* check if outside bounds of this continuum source */
478  if( ip > rfield.nflux )
479  {
480  return;
481  }
482 
483  ASSERT( ip> 0 );
484  ASSERT( ots >= 0. );
485  ASSERT( ip <= rfield.nupper );
486 
487  /* add locally destroyed flux of photons to continuum OTS array */
488  /* check whether local gas opacity (units cm-1) is positive, if so
489  * convert continuum destruction rate into ots rate by dividing by it */
490  if( opac.opacity_abs[ip-1] > 0. )
491  {
492  rfield.otsconNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
493  /*if( ots>0. ) fprintf(ioQQQ,
494  "buggg ip %li ots %.2e opac %.2e cont_phot_destroyed %.2e\n",
495  ip, ots , opac.opacity_abs[ip-1] ,rfield.otsconNew[ip-1]);*/
496  }
497  return;
498 }
499 
500 #if 0
501 /* check whether Lya OTS rate is oscillating */
502 STATIC bool lgLyaOTS_oscil( realnum ots_new )
503 {
504  static realnum old , older;
505  bool lgReturn = false;
506 
507  if( !conv.nTotalIoniz || conv.lgSearch )
508  {
509  /* this is very first call on this iteration, zero everything out */
510  old = 0.;
511  older = 0.;
512  }
513 
514  /* is sign of change oscillating? */
515  if( conv.nPres2Ioniz>1 && (ots_new-old) * (old-older) < 0.)
516  lgReturn = true;
517 
518  older = old;
519  old = ots_new;
520  return( lgReturn );
521 }
522 
523  /* >>chng 04 jan 26, check whether Lya is oscillating */
524  if( lgLyaOTS_oscil( rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ) )
525  {
526  fprintf(ioQQQ,"DEBUG lya ots %.2f %.3e\n",
527  fnzone,
529  }
530 #endif
531 /* =================================================================== */
532 
533 /*RT_OTS_Update update ots rates, called in ConvBase */
535  /* summed ots rates */
536  double *SumOTS ,
537  /* array index for constituent that changed the most
538  long int *ipOTSchange, */
539  /* limit on how large a change to allow, 0 for no limit */
540  double BigFrac)
541 {
542  long int i;
543 
544  static bool lgNeedSpace=true;
545  static double *old_OTS_line_x_opac , *old_OTS_cont_x_opac;
546  realnum FacBig , FacSml;
547 
548  DEBUG_ENTRY( "RT_OTS_Update()" );
549 
550  if( BigFrac <= 0. )
551  BigFrac = 10.;
552  FacBig = (realnum)(1. + BigFrac);
553  FacSml = 1.f/FacBig;
554 
555  /* create space to save arrays if needed, one time per coreload */
556  if( lgNeedSpace )
557  {
558  old_OTS_line_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
559  old_OTS_cont_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
560  }
561  lgNeedSpace = false;
562 
563  *SumOTS = 0.;
564 
565  /* option to kill ots rates with no ots lines command */
566  if( rfield.lgKillOTSLine )
567  {
568  for( i=0; i < rfield.nflux; i++ )
569  {
570  rfield.otslinNew[i] = 0.;
571  }
572  }
573 
574  /* during search phase set ots to current values */
575  if( nzone==0 )
576  {
577  for( i=0; i < rfield.nflux; i++ )
578  {
579  old_OTS_line_x_opac[i] = rfield.otslinNew[i] * opac.opacity_abs[i];
580  old_OTS_cont_x_opac[i] = rfield.otsconNew[i] * opac.opacity_abs[i];
581  }
582  }
583 
584  /* remember largest change in ots rates */
585  *SumOTS = 0.;
586  /* now update new ots rates */
587  for( i=0; i < rfield.nflux; ++i )
588  {
589  double OTS_line_x_opac = rfield.otslinNew[i] * opac.opacity_abs[i];
590  double OTS_cont_x_opac = rfield.otsconNew[i] * opac.opacity_abs[i];
591  double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );
592 
593  /* get new ots line rates for each cell, but don't let change be too big */
594  if( OTS_line_x_opac > old_OTS_line_x_opac[i]*FacBig )
595  {
596  rfield.otslin[i] = rfield.otslin[i]*FacBig;
597  }
598  else if( OTS_line_x_opac < old_OTS_line_x_opac[i]*FacSml )
599  {
600  rfield.otslin[i] = rfield.otslin[i]*FacSml;
601  }
602  else
603  {
604  rfield.otslin[i] = rfield.otslinNew[i];
605  }
606 
607  /* get new ots continuum rates for each cell */
608  if( OTS_cont_x_opac > old_OTS_cont_x_opac[i]*FacBig )
609  {
610  rfield.otscon[i] = rfield.otscon[i]*FacBig;
611  }
612  else if( OTS_cont_x_opac < old_OTS_cont_x_opac[i]*FacSml )
613  {
614  rfield.otscon[i] = rfield.otscon[i]*FacSml;
615  }
616  else
617  {
618  rfield.otscon[i] = rfield.otsconNew[i];
619  }
620 
621  /* zero out the new otscon vector*/
622  rfield.otsconNew[i] = 0.;
623 
624  /* zero out the new otslin vector*/
625  rfield.otslinNew[i] = 0.;
626 
627  /* now update to new values */
628  old_OTS_line_x_opac[i] = rfield.otslin[i] * opac.opacity_abs[i];
629  old_OTS_cont_x_opac[i] = rfield.otscon[i] * opac.opacity_abs[i];
630 
631  /* this is local ots continuum created by destroyed diffuse continua,
632  * currently only two-photon */
633  rfield.ConOTS_local_OTS_rate[i] = (realnum)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);
634 
635  /*if( i==14 )
636  {
637  fprintf(ioQQQ,"DEBUG OTS rate %.2e %.2e %.2e %.2e \n", *SumOTS ,
638  rfield.otscon[14] , rfield.otslin[14],opac.opacity_abs[14]);
639  fflush(ioQQQ);
640  }*/
641 
642  /* remember sum of ots rates for convergence criteria */
643  *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
644 
648 
651 
652  }
653 
654 
655  /* >>chng 02 oct 18, add this */
656  /* sum of accumulated flux from particular frequency to infinity */
657  rfield.flux_accum[rfield.nflux-1] = 0;
658  for( i=1; i < rfield.nflux; i++ )
659  {
662  }
663  /* this has to be positive since is sum of all photons in SummedCon */
664  ASSERT( rfield.flux_accum[0]> 0. );
665 
666  /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
667  * between plasma frequency and energy where brems is thin */
668  ASSERT(rfield.ipPlasma>0 );
669 
670  /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */
671  for( i=0; i < rfield.ipPlasma-1; i++ )
672  {
673  rfield.otscon[i] = 0.;
676  rfield.otslin[i] = 0.;
677  rfield.SummedDif[i] = 0.;
678  rfield.OccNumbBremsCont[i] = 0.;
679  rfield.SummedCon[i] = 0.;
680  rfield.SummedOcc[i] = 0.;
681  rfield.ConInterOut[i] = 0.;
682  }
683  /* this loop evaluates occupation number for brems continuum,
684  * only used for induced two photon emission */
685  if( rfield.ipEnergyBremsThin > 0 )
686  {
687  for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
688  {
689  /* this corrects for opacity / optical depth in brems - brems opacity goes as
690  * energy squared. */
691  /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
692  * so need this max */
693  realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);
694 
695  /* occupation number for black body is 1/ (exp hn/kT) -1) */
696  rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
697  }
698  }
699  return;
700 
701 # if 0
702  OTSOld = 0.;
703  OTSNew = 0.;
704  BigOTSNew = 0.;
705  /* find summed ots rates, old and new */
706  for( i=0; i < rfield.nflux; i++ )
707  {
708  OTSOld += rfield.otscon[i] + rfield.otslin[i];
709  OTSNew += rfield.otsconNew[i] + rfield.otslinNew[i];
710  if( BigOTSNew < rfield.otsconNew[i] + rfield.otslinNew[i] )
711  {
712  BigOTSNew = rfield.otsconNew[i] + rfield.otslinNew[i];
713  }
714  }
715 
716  /* we now have old and new rates, what is the ratio, and by how much will
717  * we allow this to change? */
718  if( BigFrac == 0. || conv.lgSearch || OTSOld < SMALLFLOAT )
719  {
720  /* this branch, want a clear update of ots rates, either because
721  * requested with BigFrac == 0, or we are in search phase */
722  FracOld = 0.;
723  }
724  else
725  {
726  if( OTSNew > OTSOld )
727  {
728  chng = fabs(1. - OTSOld/SDIV(OTSNew) );
729  }
730  else
731  {
732  chng = fabs(1. - OTSNew/SDIV(OTSOld) );
733  }
734 
735  if( chng < BigFrac )
736  {
737  /* this branch, old and new do not differ by much, so use new */
738  FracOld = 0.;
739  }
740  else
741  {
742  /* this branch, too large a difference between old and new, cap it to BigFrac */
743  FracOld = (1. - BigFrac / chng);
744  ASSERT( FracOld >= 0. );
745  FracOld = MIN2( 0.25 , FracOld );
746  }
747  }
748 
749  /* fraction old and new ots rates */
750  FracNew = 1. - FracOld;
751  fprintf(ioQQQ," DEBUG FracNew\t%.2e\n", FracNew);
752 
753  /* remember largest change in ots rates */
754  changeOTS = 0.;
755  /*fprintf(ioQQQ," sumcontinuum zone%li 1168=%e\n", nzone,rfield.otslin[1167]);*/
756 
757  for( i=0; i < rfield.nflux; i++ )
758  {
759  /* >>chng 01 feb 01, define inverse opacity in safe manner */
760  double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );
761 
762  /* remember which energy had largest change in ots rates */
763  if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS)
764  {
765  changeOTS =
767  }
768 
769  /* >>chng 01 apr 10, only change ots with means if FracOld not zero,
770  * this is to avoid taking means when this routine called by StartEndIter
771  * to reset sums of rates */
772  if( BigFrac > 0. && conv.nTotalIoniz > 0 )
773  {
774  /* the New vectors are the ones updated by the AddOTS routines,
775  * and have the current OTS rates. The otscon and otslin vectors
776  * have the rates from the previous refresh of the New vectors*/
777  /* here is the refresh. all were initially set to zero in call to
778  * RT_OTS_Zero (below) from zero */
779  rfield.otscon[i] = (realnum)((rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew));
780 
781  /* this is local ots continuum created by destroyed diffuse continua,
782  * currently only two-photon */
783  /* >>chng 00 dec 27, define this continuum and do not count 2-photon in old var */
784  rfield.ConOTS_local_OTS_rate[i] = (realnum)(rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);
785 
786  /* current ots will be average of old and new */
787  rfield.otslin[i] = (realnum)((rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew));
788  }
789 
790  /* zero out the new otscon vector*/
791  rfield.otsconNew[i] = 0.;
792 
793  /* zero out the new otslin vector*/
794  rfield.otslinNew[i] = 0.;
795 
796  /* remember sum of ots rates for convergence criteria */
797  *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
798  {
799  /* following should be set true to print strongest ots contributors */
800  /*@-redef@*/
801  enum {DEBUG_LOC=false};
802  /*@+redef@*/
803  if( DEBUG_LOC && (nzone==378) )
804  {
805  if( conv.nPres2Ioniz > 3 )
806  cdEXIT( EXIT_FAILURE );
807  fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
809  rfield.anu[i],
810  opac.opacity_abs[i],
811  rfield.otscon[i],
812  rfield.otslin[i]);
813  }
814  }
815 
816  /* >>chng 00 dec 27, include ConOTS_local_OTS_rate */
817  /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
821 
824  }
825 
826  ASSERT(*SumOTS >= 0. );
827 
828  /* >>chng 02 oct 18, add this */
829  /* sum of accumulated flux from particular frequency to infinity */
830  rfield.flux_accum[rfield.nflux-1] = 0;
831  for( i=1; i < rfield.nflux; i++ )
832  {
835  }
836  /* this has to be positive since is sum of all photons in SummedCon */
837  ASSERT( rfield.flux_accum[0]> 0. );
838 
839  /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
840  * between plasma frequency and energy where brems is thin */
841  ASSERT(rfield.ipPlasma>0 );
842 
843  /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */
844  for( i=0; i < rfield.ipPlasma-1; i++ )
845  {
846  rfield.otscon[i] = 0.;
849  rfield.otslin[i] = 0.;
850  rfield.SummedDif[i] = 0.;
851  rfield.OccNumbBremsCont[i] = 0.;
852  rfield.SummedCon[i] = 0.;
853  rfield.SummedOcc[i] = 0.;
854  rfield.ConInterOut[i] = 0.;
855  }
856  /* this loop evaluates occupation number for brems continuum,
857  * only used for induced two photon emission */
858  if( rfield.ipEnergyBremsThin > 0 )
859  {
860  for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
861  {
862  /* this corrects for opacity / optical depth in brems - brems opacity goes as
863  * energy squared. */
864  /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
865  * so need this max */
866  realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);
867 
868  /* occupation number for black body is 1/ (exp hn/kT) -1) */
869  rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
870  }
871  }
872 
873  {
874  /* following should be set true to print strongest ots contributors */
875  /*@-redef@*/
876  enum {DEBUG_LOC=false};
877  /*@+redef@*/
878  if( DEBUG_LOC && (nzone>0) )
879  {
880  double BigOTS;
881  int ipOTS=0;
882  BigOTS = -1.;
883  /* find the biggest ots contributor */
884  /*for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]; i < rfield.nflux; i++ )*/
885  for( i=0; i < rfield.nflux; i++ )
886  {
887  if( (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i] > BigOTS )
888  {
889  BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
890  ipOTS = (int)i;
891  }
892  }
893  fprintf(ioQQQ,
894  " sumcontinuunots zone %li SumOTS %.2e %.2eRyd opc:%.2e OTS%.2e lin%.2e con:%.2e %s %s cell%i FracNew %.2f \n",
895  nzone,
896  *SumOTS,
897  rfield.anu[ipOTS] ,
898  opac.opacity_abs[ipOTS],
899  BigOTS ,
900  rfield.otslin[ipOTS],
901  rfield.otscon[ipOTS] ,
902  /* line label */
903  rfield.chLineLabel[ipOTS] ,
904  /* cont label*/
905  rfield.chContLabel[ipOTS] ,
906  ipOTS,
907  FracNew);
908  }
909  }
910  return;
911 # endif
912 }
913 
914 /* =================================================================== */
915 
916 /*RT_OTS_Zero zero out some vectors - this is only called when code
917  * initialized by ContSetIntensity */
918 void RT_OTS_Zero( void )
919 {
920  long int i;
921 
922  DEBUG_ENTRY( "RT_OTS_Zero()" );
923 
924  /* this loop goes up to nflux itself (<=) since the highest cell
925  * will be used to pass unity through the code to verify integrations */
926  for( i=0; i <= rfield.nflux; i++ )
927  {
928  rfield.SummedDif[i] = 0.;
929  /* the four main ots vectors */
930  rfield.otscon[i] = 0.;
931  rfield.otslin[i] = 0.;
932  rfield.otslinNew[i] = 0.;
933  rfield.otsconNew[i] = 0.;
934  /* >>chng 05 mar 03, add this */
935  rfield.trans_coef_zone[i] = 1.f;
936  rfield.trans_coef_total[i] = 1.f;
937 
938  rfield.ConInterOut[i] = 0.;
939  rfield.outlin[i] = 0.;
940  rfield.outlin_noplot[i] = 0.;
941  rfield.SummedDif[i] = 0.;
942  /* "zero" for the summed con will be just the incident radiation field */
943  rfield.SummedCon[i] = rfield.flux[i];
947  }
948  return;
949 }
950 
951 /* =================================================================== */
952 
953 /*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
954 void RT_OTS_ChkSum(long int ipPnt)
955 {
956  static long int nInsane=0;
957  long int i;
958  double phisig;
959 # define LIM_INSAME_PRT 30
960 
961  DEBUG_ENTRY( "RT_OTS_ChkSum()" );
962 
963  /* this entire sub is a sanity check */
964  /* >>chng 02 jul 23, lower bound from 0 to rfield.ipEnergyBremsThin - since now
965  * set radiation field to black body below this energy */
966  for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ )
967  {
968  phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly +
969  rfield.outlin[i]+
972  /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero whild
973  * phisig is just above small float */
974  if( phisig > 0. && rfield.SummedDif[i]> 0.)
975  {
976  if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 )
977  {
978  ++nInsane;
979  /* limit amount of printout - in many failures would print entire
980  * continuum array */
981  if( nInsane < LIM_INSAME_PRT )
982  {
983  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n",
984  rfield.anu[i], rfield.SummedDif[i]/phisig - 1., i );
985  fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n",
986  rfield.SummedDif[i], phisig );
987  fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n",
990  fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n",
992  }
993  }
994  }
995 
996  phisig += rfield.flux[i];
997  /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero when
998  * phisig is just above small float */
999  if( phisig > 0. && rfield.SummedDif[i]> 0. )
1000  {
1001  if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 )
1002  {
1003  ++nInsane;
1004  /* limit amount of printout - in many failures would print entire
1005  * continuum array */
1006  if( nInsane < LIM_INSAME_PRT )
1007  {
1008  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum %3ld, insane SummedCon at energy %.5e error=%.2e i=%ld\n",
1009  ipPnt, rfield.anu[i], rfield.SummedCon[i]/phisig - 1., i );
1010  fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n",
1011  rfield.SummedCon[i], phisig );
1012  fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n",
1015  fprintf( ioQQQ, " line continuum here are %s %s\n",
1017  );
1018  }
1019  }
1020  }
1021  }
1022 
1023  if( nInsane > 0 )
1024  {
1025  fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n");
1026  /* TotalInsanity exits after announcing a problem */
1027  TotalInsanity();
1028  }
1029  return;
1030 }
1031 
1032 /* =================================================================== */
1033 
1034 /*RT_OTS_PrtRate print continuum and line ots rates when trace ots is on */
1036  /* weakest rate to print */
1037  double weak ,
1038  /* flag, 'c' continuum, 'l' line, 'b' both */
1039  int chFlag )
1040 {
1041  long int i;
1042 
1043  DEBUG_ENTRY( "RT_OTS_PrtRate()" );
1044 
1045  /* arg must be one of these three */
1046  ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' );
1047 
1048  /*
1049  * both printouts have cell number (on C array scale)
1050  * energy in ryd
1051  * the actual value of the ots rate
1052  * the ots rate relative to the continuum at that energy
1053  * rate times opacity
1054  * all are only printed if greater than weak
1055  */
1056 
1057  /*===================================================================*/
1058  /* first print ots continua */
1059  /*===================================================================*/
1060  if( chFlag == 'c' || chFlag == 'b' )
1061  {
1062  fprintf( ioQQQ, " DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n",
1063  weak,fnzone ,TorF(conv.lgConvIoniz) );
1064 
1065  for( i=0; i < rfield.nupper; i++ )
1066  {
1067  if( rfield.otscon[i]*opac.opacity_abs[i] > weak )
1068  {
1069  fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s \n",
1070  i,
1071  rfield.anu[i],
1072  rfield.otscon[i],
1073  opac.opacity_abs[i],
1074  rfield.otscon[i]*opac.opacity_abs[i],
1075  rfield.chContLabel[i]);
1076 
1077  }
1078  }
1079  }
1080 
1081  /*===================================================================*/
1082  /* second print ots line rates */
1083  /*===================================================================*/
1084  if( chFlag == 'l' || chFlag == 'b' )
1085  {
1086  fprintf( ioQQQ, " DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n",
1087  weak,fnzone,TorF(conv.lgConvIoniz) );
1088 
1089  for( i=0; i < rfield.nupper; i++ )
1090  {
1091  if( rfield.otslin[i]*opac.opacity_abs[i] > weak )
1092  {
1093  fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n",
1094  i,
1095  rfield.anu[i],
1096  rfield.otslin[i],
1097  opac.opacity_abs[i],
1098  rfield.otslin[i]*opac.opacity_abs[i],
1099  rfield.chLineLabel[i] ,
1100  rfield.line_count[i] );
1101  }
1102  }
1103  }
1104  return;
1105 }

Generated for cloudy by doxygen 1.8.3.1