cloudy trunk
Loading...
Searching...
No Matches
opacity_addtotal.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/*OpacityAddTotal derive total opacity for this position,
4 * called by ConvBase */
5#include "cddefines.h"
6#include "physconst.h"
7#include "iso.h"
8#include "ipoint.h"
9#include "grainvar.h"
10#include "ca.h"
11#include "rfield.h"
12#include "oxy.h"
13#include "mole.h"
14#include "h2.h"
15#include "h2_priv.h"
16#include "hmi.h"
17#include "dense.h"
18#include "atoms.h"
19#include "conv.h"
20#include "ionbal.h"
21#include "trace.h"
22#include "hmi.h"
23#include "phycon.h"
24#include "opacity.h"
25#include "taulines.h"
26
28{
29 long int i,
30 ion,
31 ipHi,
32 ipISO,
33 ipop,
34 limit,
35 low,
36 nelem,
37 n;
38 double DepartCoefInv ,
39 fac,
40 sum;
41 realnum SaveOxygen1 ,
42 SaveCarbon1;
43
44 DEBUG_ENTRY( "OpacityAddTotal()" );
45
46 /* OpacityZero will zero out scattering and absorption opacities,
47 * and set OldOpacSave to opac to save it */
49
50 /* free electron scattering opacity, Compton recoil energy loss */
51 for( i=0; i < rfield.nflux; i++ )
52 {
53 /* scattering part of total opacity */
54 opac.opacity_sct[i] += opac.OpacStack[i-1+opac.iopcom]*
55 dense.eden;
56 }
57
58 /* opacity due to Compton bound recoil ionization */
59 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
60 {
61 if( dense.lgElmtOn[nelem] )
62 {
63 for( ion=0; ion<nelem+1; ++ion )
64 {
65 realnum factor = dense.xIonDense[nelem][ion];
66 /*>>chng 05 nov 26, add molecular hydrogen assuming same
67 * as two free hydrogen atoms - hence mult density by two
68 *>>KEYWORD H2 bound Compton ionization */
69 if( nelem==ipHYDROGEN )
70 factor += hmi.H2_total*2.f;
71 if( factor > 0. )
72 {
73 // loop_min and loop_max are needed to work around a bug in icc 10.0
74 long loop_min = ionbal.ipCompRecoil[nelem][ion]-1;
75 long loop_max = rfield.nflux;
76 /* ionbal.nCompRecoilElec number of electrons in valence shell
77 * that can compton recoil ionize */
78 factor *= ionbal.nCompRecoilElec[nelem-ion];
79 for( i=loop_min; i < loop_max; i++ )
80 {
81 /* add in bound hydrogen electron scattering, treated as absorption */
82 opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor;
83 }
84 }
85 }
86 }
87 }
88
89 /* opacity due to pair production - does not matter what form these
90 * elements are in */
92 sum = dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM];
93 OpacityAdd1Subshell(opac.ioppr,opac.ippr,rfield.nflux,(realnum)sum,'s');
94
95 /* free-free free free brems emission for all ions */
96
97 /* >>chng 02 jul 21, use full set of ions and gaunt factor */
98 /* ipEnergyBremsThin is index to energy where gas goes optically thin to brems,
99 * so this loop is over optically thick frequencies */
100 static double *TotBremsAllIons;
101 static bool lgFirstTime=true;
102 double BremsThisEner,bfac, sumion[LIMELM+1];
103 long int ion_lo , ion_hi;
104
105 if( lgFirstTime )
106 {
107 /* rfield.nupper will not change in one coreload, so just malloc this once */
108 TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
109 lgFirstTime = false;
110 }
111 bfac = (dense.eden/1e20)/phycon.sqrte/1e10;
112 /* gaunt factors depend only on photon energy and ion charge, so do
113 * sum of ions here before entering into loop over photon energy */
114 sumion[0] = 0.;
115 for(ion=1; ion<=LIMELM; ++ion )
116 {
117 sumion[ion] = 0.;
118 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
119 {
120 if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
121 {
122 sumion[ion] += dense.xIonDense[nelem][ion];
123 }
124 }
125 /* now include the charge, density, and temperature */
126 sumion[ion] *= POW2((double)ion)*bfac;
127 }
128
129 /* add molecular ions */
130 for( long ipMol = 0; ipMol<mole_global.num_calc; ipMol++ )
131 {
132 if( !mole_global.list[ipMol]->isMonatomic() && mole_global.list[ipMol]->charge > 0 &&
133 mole_global.list[ipMol]->parentLabel.empty() &&
134 mole_global.list[ipMol]->label != "H2+" &&
135 mole_global.list[ipMol]->label != "H3+" )
136 {
137 ASSERT( mole_global.list[ipMol]->charge < LIMELM+1 );
138 sumion[mole_global.list[ipMol]->charge] += mole.species[ipMol].den * POW2((realnum)mole_global.list[ipMol]->charge)*bfac;
139 }
140 }
141
142 /* now find lowest and highest ion we need to consider - following loop is over
143 * full continuum and eats time
144 * >>chng 05 oct 20, following had tests on ion being within bounds, bug,
145 * changed to ion_lo and ion_hi */
146 ion_lo = 1;
147 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
148 ++ion_lo;
149 ion_hi = LIMELM;
150 while( sumion[ion_hi]==0 && ion_hi>0 )
151 --ion_hi;
152
153 const molezone *MH2p = findspecieslocal("H2+");
154 const molezone *MH3p = findspecieslocal("H3+");
155 for( i=0; i < rfield.nflux; i++ )
156 {
157 /* do hydrogen first, before main loop since want to add on H- brems */
158 nelem = ipHYDROGEN;
159 ion = 1;/* >>chng 02 nov 07, add charged ions as H+ */
160
161 BremsThisEner = bfac * rfield.gff[ion][i] * (dense.xIonDense[ipHYDROGEN][1]+MH2p->den+MH3p->den);
162 /* for case of hydrogen, add H- brems - OpacStack contains the ratio
163 * of the H- to H brems cross section - multiply by this and H(1s) population */
164 BremsThisEner += bfac * rfield.gff[ion][i] * opac.OpacStack[i-1+opac.iphmra] * iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
165
166 /* there is at least one standard test (grainlte) which has ZERO ionization -
167 * this assert will pass that test (the ==0) and also the usual case */
168 /* ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons[i] > 0.) );*/
169
170 /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
171 * using gaunt factors from rfield.gff. */
172 /* >>chng 05 jul 11 reorganize loop for speed */
173 for(ion=ion_lo; ion<=ion_hi; ++ion )
174 {
175 BremsThisEner += sumion[ion]*rfield.gff[ion][i];
176 }
177 ASSERT( !isnan( BremsThisEner ) );
178 TotBremsAllIons[i] = BremsThisEner;
179
180 /* >>>chng 05 jul 12, bring these two loops together */
181 if( rfield.ContBoltz[i] < 0.995 )
182 {
183 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
184 (1. - rfield.ContBoltz[i]);
185 }
186 else
187 {
188 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
189 rfield.anu[i]*TE1RYD/phycon.te;
190 }
191 opac.FreeFreeOpacity[i] = TotBremsAllIons[i];
192 /* >>chng 02 jul 23, from >0 to >=0, some models have no ionization,
193 * like grainlte.in */
194 /*ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );*/
195 opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
196 }
197
198
199 /* H minus absorption, with correction for stimulated emission */
200 if( hmi.hmidep > SMALLFLOAT )
201 {
202 DepartCoefInv = 1./hmi.hmidep;
203 }
204 else
205 {
206 /* the hmidep departure coef can become vastly small in totally
207 * neutral gas (no electrons) */
208 DepartCoefInv = 1.;
209 }
210
211 limit = iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon;
212 if(limit > rfield.nflux)
213 limit = rfield.nflux;
214
215 const molezone *MHm = findspecieslocal("H-");
216 for( i=hmi.iphmin-1; i < limit; i++ )
217 {
218 double factor;
219 factor = 1. - rfield.ContBoltz[i]*DepartCoefInv;
220 if(factor > 0)
221 opac.opacity_abs[i] += opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
222 MHm->den*factor;
223 }
224
225 /* H2P h2plus bound free opacity */
226 limit = opac.ih2pnt[1];
227 if(limit > rfield.nflux)
228 limit = rfield.nflux;
229 for( i=opac.ih2pnt[0]-1; i < limit; i++ )
230 {
231 opac.opacity_abs[i] += MH2p->den*opac.OpacStack[i-opac.ih2pnt[0]+
232 opac.ih2pof];
233 ASSERT( !isnan( opac.opacity_abs[i] ) );
234 }
235
236 /* H2 continuum dissociation opacity */
237 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
238 {
239 if( (*diatom)->lgEnabled && mole_global.lgStancil )
240 {
241 for( vector< diss_tran >::iterator tran = (*diatom)->Diss_Trans.begin(); tran != (*diatom)->Diss_Trans.end(); ++tran )
242 {
243 long lower_limit = ipoint(tran->energies[0]);
244 long upper_limit = ipoint(tran->energies.back());
245 upper_limit = MIN2( upper_limit, rfield.nflux-1 );
246 for(i = lower_limit; i <= upper_limit; ++i)
247 {
248 opac.opacity_abs[i] += (*diatom)->MolDissocOpacity( *tran, rfield.anu[i]);
249 }
250 }
251 }
252 }
253
254 /* get total population of hydrogen ground to do Rayleigh scattering */
255 if( dense.xIonDense[ipHYDROGEN][1] <= 0. )
256 {
257 fac = dense.xIonDense[ipHYDROGEN][0];
258 }
259 else
260 {
261 fac = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
262 }
263
264 /* Ly a damp wing opac (Rayleigh scattering) */
265 limit = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon;
266 if(limit > rfield.nflux)
267 limit = rfield.nflux;
268 for( i=0; i < limit; i++ )
269 {
270 opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]);
271 }
272
273 /* remember largest correction for stimulated emission */
274 if( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].DepartCoef > 1e-30 && !conv.lgSearch )
275 {
276 realnum factor;
277 factor = (realnum)(rfield.ContBoltz[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1]/iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].DepartCoef);
278 if(opac.stimax[0] < factor)
279 opac.stimax[0] = factor;
280 }
281
282 if( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].DepartCoef > 1e-30 && !conv.lgSearch )
283 {
284 realnum factor;
285 factor = (realnum)(rfield.ContBoltz[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1]/iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].DepartCoef);
286 if(opac.stimax[1] < factor)
287 opac.stimax[1] = factor;
288 }
289
290# if 0
291 /* check whether hydrogen or Helium singlets mased, if not in search mode */
292 if( !conv.lgSearch )
293 {
294 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).PopOpc() < 0. )
295 {
296 hydro.lgHLyaMased = true;
297 }
298 }
299# endif
300
301 /* >>chng 05 nov 25, use Yan et al. H2 photo cs
302 * following reference gives cross section for all energies
303 * >>refer H2 photo cs Yan, M., Sadeghpour, H.R., & Dalgarno, A., 1998, ApJ, 496, 1044
304 * Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 */
305 /* >>chng 02 jan 16, approximate inclusion of H_2 photoelectric opacity
306 * include H_2 in total photoelectric opacity as twice H0 cs */
307 /* set lower and upper limits to this range */
308 /*>>KEYWORD H2 photoionization opacity */
309 low = h2.ip_photo_opac_thresh;
310 ipHi = rfield.nupper;
311 ipop = h2.ip_photo_opac_offset;
312 /* OpacityAdd1Subshell just returns for static opacities if opac.lgRedoStatic not set*/
313 /* >>chng 05 nov 27, change on nov 25 had left 2*density from H0, so
314 * twice the H2 density was used - * also changed to static opacity
315 * this assumes that all v,J levels contribute the same opacity */
316 OpacityAdd1Subshell( ipop , low , ipHi , hmi.H2_total , 's' );
317
318 /*>>KEYWORD CO photoionization opacity */
319 /* include photoionization of CO - assume C and O in CO each have
320 * same photo cs as atom - this should only be significant in highly
321 * shielded regions where only very hard photons penetrate
322 * also H2O condensed onto grain surfaces - very important deep in cloud */
323 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
324 SaveCarbon1 = dense.xIonDense[ipCARBON][0];
325 /* atomic C and O will include CO during the heating sum loop */
326 fixit(); /* This breaks the invariant for total mass.
327 *
328 * In any case, is CO the only species which contributes?
329 * -- might expect all other molecular species to do so,
330 * i.e. the item added should perhaps be
331 * dense.xMolecules[nelem]) -- code duplicated in heatsum
332 */
333 dense.xIonDense[ipOXYGEN][0] += (realnum) (findspecieslocal("CO")->den + findspecieslocal("H2Ogrn")->den);
334 dense.xIonDense[ipCARBON][0] += (realnum) (findspecieslocal("CO")->den);
335
336 /* following loop adds standard opacities for first 30 elements
337 * most heavy element opacity added here */
338 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
339 {
340 /* this element may be turned off */
341 if( dense.lgElmtOn[nelem] )
342 {
343 OpacityAdd1Element(nelem);
344 }
345 }
346
347 /* now reset the abundances */
348 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
349 dense.xIonDense[ipCARBON][0] = SaveCarbon1;
350
351 /* following are opacities due to specific excited levels */
352
353 /* nitrogen opacity
354 * excited level of N+ */
355 OpacityAdd1Subshell(opac.in1[2],opac.in1[0],opac.in1[1],
356 dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' );
357
358 /* oxygen opacity
359 * excited level of Oo */
360 OpacityAdd1Subshell(opac.ipo1exc[2],opac.ipo1exc[0],opac.ipo1exc[1],
361 dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v');
362
363 /* O2+ excited states */
364 OpacityAdd1Subshell(opac.ipo3exc[2],opac.ipo3exc[0],opac.ipo3exc[1],
365 dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v');
366
367 OpacityAdd1Subshell(opac.ipo3exc3[2],opac.ipo3exc3[0],opac.ipo3exc3[1],
368 dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v');
369
370 /* magnesium opacity
371 * excited level of Mg+ */
372 OpacityAdd1Subshell(opac.ipOpMgEx,opac.ipmgex,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon,
373 dense.xIonDense[ipMAGNESIUM][1]* atoms.popmg2,'v');
374
375 /* calcium opacity
376 * photoionization of excited levels of Ca+ */
377 OpacityAdd1Subshell(opac.ica2op,opac.ica2ex[0],opac.ica2ex[1],
378 ca.popca2ex,'v');
379
380 /*******************************************************************
381 *
382 * complete evaluation of total opacity by adding in the static part and grains
383 *
384 *******************************************************************/
385
386 /* this loop defines the variable iso_sp[ipH_LIKE][nelem].fb[n].ConOpacRatio,
387 * the ratio of not H to Hydrogen opacity. for grain free environments
388 * at low densities this is nearly zero. The correction includes
389 * stimulated emission correction */
390 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
391 {
392 for( nelem=ipISO; nelem < LIMELM; nelem++ )
393 {
394 /* this element may be turned off */
395 if( dense.lgElmtOn[nelem] )
396 {
397 /* this branch is for startup only */
398 if( nzone < 1 )
399 {
400 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */
401 for( n=0; n < iso_sp[ipISO][nelem].numLevels_local; n++ )
402 {
403 if(iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon < rfield.nflux )
404 {
405 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
406 * greater than this - caused oscillations as opacity fell below
407 * and around this value - change to greater than 0 */
408 /*if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 1e-30 )*/
409 if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 0. )
410 {
411 /* >>chng 02 may 06, use general form of threshold cs */
412 /*double t1 = atmdat_sth(n)/(POW2(nelem+1.-ipISO));*/
413 long int ip = iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon;
414 double t2 = csphot(
415 ip ,
416 ip ,
417 iso_sp[ipISO][nelem].fb[n].ipOpac );
418
419 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 1.-
420 (iso_sp[ipISO][nelem].st[n].Pop()*t2)/
421 opac.opacity_abs[ip-1];
422 }
423 else
424 {
425 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 0.;
426 }
427 }
428 }
429 }
430 /* end branch for startup only, start branch for all zones including startup */
431 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */
432 for( n=0; n < iso_sp[ipISO][nelem].numLevels_local; n++ )
433 {
434 /* ratios of other to total opacity for continua of all atoms done with iso model */
435 if(iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon < rfield.nflux )
436 {
437 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
438 * greater than this - caused oscillations as opacity fell below
439 * and around this value - change to greater than 0 */
440 /*if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 1e-30 )*/
441 if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 0. )
442 {
443 /* first get departure coef */
444 if( iso_sp[ipISO][nelem].fb[n].DepartCoef > 1e-30 && (!conv.lgSearch ) )
445 {
446 /* this is the usual case, use inverse of departure coef */
447 fac = 1./iso_sp[ipISO][nelem].fb[n].DepartCoef;
448 }
449 else if( conv.lgSearch )
450 {
451 /* do not make correction for stim emission during search
452 * for initial temperature solution, since trys are very non-equil */
453 fac = 0.;
454 }
455 else
456 {
457 fac = 1.;
458 }
459
462 /* now get opaicty ratio with correction for stimulated emission */
463 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
464 * greater than this - caused oscillations as opacity fell below
465 * and around this value - change to greater than 0 */
466 /*if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 1e-30 )*/
467 if( opac.opacity_abs[iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon-1] > 0. )
468 {
469 /* >>chng 02 may 06, use general form of threshold cs */
470 long int ip = iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon;
471
472 double t2 = csphot(
473 ip ,
474 ip ,
475 iso_sp[ipISO][nelem].fb[n].ipOpac );
476
477 double opacity_this_species =
478 iso_sp[ipISO][nelem].st[n].Pop()*t2*
479 (1. - fac*rfield.ContBoltz[ip-1]);
480
481 double opacity_fraction = 1. - opacity_this_species / opac.opacity_abs[ip-1];
482 if(opacity_fraction < 0)
483 opacity_fraction = 0.;
484
485 /* use mean of old and new ratios */
486 iso_sp[ipISO][nelem].fb[n].ConOpacRatio =
487 iso_sp[ipISO][nelem].fb[n].ConOpacRatio* 0.75 + 0.25*opacity_fraction;
488
489 if(iso_sp[ipISO][nelem].fb[n].ConOpacRatio < 0.)
490 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 0.;
491 }
492 else
493 {
494 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 0.;
495 }
496 }
497 else
498 {
499 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 0.;
500 }
501 }
502 else
503 {
504 iso_sp[ipISO][nelem].fb[n].ConOpacRatio = 0.;
505 }
506 }
507 }
508 }
509 }
510
511 /* add dust grain opacity if dust present */
512 if( gv.lgDustOn() )
513 {
514 /* generate current grain opacities since may be function of depth */
515 /* >>chng 01 may 11, removed code to update grain opacities, already done by GrainChargeTemp */
516 for( i=0; i < rfield.nflux; i++ )
517 {
518 /* units cm-1 */
519 opac.opacity_sct[i] += gv.dstsc[i]*dense.gas_phase[ipHYDROGEN];
520 opac.opacity_abs[i] += gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
521 }
522 }
523
524 /* check that opacity is sane */
525 for( i=0; i < rfield.nflux; i++ )
526 {
527 /* OpacStatic was zeroed in OpacityZero, incremented in opacityadd1subshell */
528 opac.opacity_abs[i] += opac.OpacStatic[i];
529 /* make sure that opacity is positive */
530 /*ASSERT( opac.opacity_abs[i] > 0. );*/
531 }
532
533 /* compute gas albedo here */
534 for( i=0; i < rfield.nflux; i++ )
535 {
536 opac.albedo[i] = opac.opacity_sct[i]/
537 (opac.opacity_sct[i] + opac.opacity_abs[i]);
538 }
539
540 /* during search phase set old opacity array to current value */
541 if( conv.lgSearch )
543
544 if( trace.lgTrace )
545 {
546 fprintf( ioQQQ, " OpacityAddTotal returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e HeI,II:%10.2e%10.2e\n",
547 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ConOpacRatio,
548 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ConOpacRatio,
549 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ConOpacRatio,
550 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3s].ConOpacRatio,
551 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3p].ConOpacRatio,
552 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH3d].ConOpacRatio,
553 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH4s].ConOpacRatio,
554 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH4p].ConOpacRatio,
555 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH4d].ConOpacRatio,
556 iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH4f].ConOpacRatio,
557 iso_sp[ipHE_LIKE][ipHELIUM].fb[ipHe1s1S].ConOpacRatio,
558 iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ConOpacRatio );
559 }
560
561 {
562 /* following should be set true to print strongest ots contributors */
563 enum {DEBUG_LOC=false};
564 if( DEBUG_LOC && (nzone>=378)/**/ )
565 {
566 if( nzone > 380 )
568 for( i=0; i<rfield.nflux; ++i )
569 {
570 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
571 conv.nPres2Ioniz,
572 rfield.anu[i],
573 opac.opacity_abs[i],
574 rfield.otscon[i],
575 rfield.otslin[i]);
576 }
577 }
578 }
579 return;
580}
t_atoms atoms
Definition atoms.cpp:5
t_ca ca
Definition ca.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
double csphot(long int inu, long int ithr, long int iofset)
Definition service.cpp:1602
#define isnan
Definition cddefines.h:620
#define ASSERT(exp)
Definition cddefines.h:578
const int ipNITROGEN
Definition cddefines.h:311
const int ipCARBON
Definition cddefines.h:310
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
const int ipMAGNESIUM
Definition cddefines.h:316
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
double den
Definition mole.h:358
long ipoint(double energy_ryd)
t_conv conv
Definition conv.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_hmi hmi
Definition hmi.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipHe1s1S
Definition iso.h:41
const int ipH3p
Definition iso.h:31
const int ipH1s
Definition iso.h:27
const int ipH4p
Definition iso.h:34
const int ipHE_LIKE
Definition iso.h:63
const int ipH3s
Definition iso.h:30
const int ipH4d
Definition iso.h:35
const int ipH2p
Definition iso.h:29
const int ipH4s
Definition iso.h:33
const int ipH4f
Definition iso.h:36
const int ipH3d
Definition iso.h:32
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_opac opac
Definition opacity.cpp:5
void OpacityZero(void)
void OpacityAdd1Element(long int ipZ)
void OpacityAdd1Subshell(long int ipOpac, long int ipLowLim, long int ipUpLim, realnum abundance, char chStat)
void OpacityZeroOld(void)
void OpacityAddTotal(void)
t_oxy oxy
Definition oxy.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double TE1RYD
Definition physconst.h:183
t_rfield rfield
Definition rfield.cpp:8
t_trace trace
Definition trace.cpp:5