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