cloudy trunk
Loading...
Searching...
No Matches
opacity_add1element.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/*OpacityAdd1Element enter total photo cross section for all subshells into opacity array */
4#include "cddefines.h"
5#include "iso.h"
6#include "rfield.h"
7#include "dense.h"
8#include "heavy.h"
9#include "opacity.h"
10#include "taulines.h"
11
13 /* nelem is 0 for H, 1 for He, etc */
14 long int nelem)
15{
16 long int ipHi,
17 ipop,
18 limit,
19 low,
20 n,
21 ion,
22 nshell;
23 char chStat;
24 double abundance;
25
26 DEBUG_ENTRY( "OpacityAdd1Element()" );
27
28 /* this routine drives OpacityAdd1Subshell to put in total opacities for all shells*/
29
30 /*begin sanity check */
31 ASSERT( (nelem >=0 ) && (nelem < LIMELM) );
32
33 /* first do simple two-level systems -
34 * this is number of species that are not treated on common iso-electronic series */
35 limit = nelem + 1 - NISO;
36 /* this can be called with hydrogen itself, in which case nelem is 0, and limit is
37 * -1 - do not do any of the simple ions */
38 limit = MAX2( 0 , limit );
39
40 /* do not include the ion stages that have complete atoms,
41 * currently H and He like iso sequences */
42 for( ion=0; ion < limit; ion++ )
43 {
44 if( dense.xIonDense[nelem][ion] > 0. )
45 {
46 /*start with static opacities, then do volatile*/
47
48 chStat = 's';
49 /* number of bound electrons */
50 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
51 {
52 /* highest shell will be volatile*/
53 if( nshell== Heavy.nsShells[nelem][ion]-1 )
54 chStat = 'v';
55 /* set lower and upper limits to this range */
56 low = opac.ipElement[nelem][ion][nshell][0];
57 ipHi = opac.ipElement[nelem][ion][nshell][1];
58 ipop = opac.ipElement[nelem][ion][nshell][2];
59 /* OpacityAdd1Subshell will not do anything if static opacities do not need to be reset*/
60 OpacityAdd1Subshell(ipop,low,ipHi,dense.xIonDense[nelem][ion] , chStat );
61 }
62 }
63 }
64
65 /* now loop over all species done as large multi-level systems */
66 /* >>chng 02 jan 17, add loop over H and He like */
67 /* ion is on the c scale, =0 for HI, =1 for HeII */
68 for( ion=limit; ion<nelem+1; ++ion )
69 {
70 /* ipISO is 0 for H-like, 1 for He-like */
71 long int ipISO = nelem-ion;
72
73 /* do multi level systems, but only if present
74 * test for nelem+1 in case atom present but not ion, test is whether the
75 * abundance of the recombined species is present */
76 /* >>chng 02 jan 17, sec dim had been nelem+1, change to ion+1 */
77 /*if( dense.xIonDense[nelem][nelem] > 0. )*/
78 if( dense.xIonDense[nelem][ion] > 0. )
79 {
80 ASSERT( ipISO < NISO );
81
82 /* do ground first, then all excited states */
83 n = 0;
84 /* abundance of recombined species, which can be zero if no ion present */
85 abundance = iso_sp[ipISO][nelem].st[n].Pop();
86
87 /* >>chng 02 may 06, add second test, had been just the chck on helium,
88 * with no option to use new soln */
89 if( abundance == 0. )
90 {
91 /* no ionized species, assume everything in ground */
92 abundance = dense.xIonDense[nelem][ion];
93 }
94
95 /* >>chng 02 jan 17, to arbitrary iso sequence */
96 /* use computed opacities and departure coef for level */
98 iso_sp[ipISO][nelem].fb[n].ipOpac,
99 iso_sp[ipISO][nelem].fb[n].ipIsoLevNIonCon,
100 /* the upper limit to the integration,
101 * ground opacity goes up to the high-energy limit of code*/
102 rfield.nflux,
103 /* the abundance of the ion */
104 abundance,
105 /* departure coef, volatile opac, always reevaluate */
106 iso_sp[ipISO][nelem].fb[n].DepartCoef , 'v' );
107
108 /* do excited levvels,
109 * this loop only if upper levels have finite population*/
110 if( iso_sp[ipISO][nelem].st[3].Pop() > 0. )
111 {
112 char chType = 'v';
113 /* always want to evaluate all opacities for n=3, 4, use static opacities for higher levels */
114 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */
115 for( long level =1; level < iso_sp[ipISO][nelem].numLevels_local; level++ )
116 {
117 if( level==iso_sp[ipISO][nelem].numLevels_max-1 )
118 chType = 'v';
119 /* above 4 is static */
120 else if( iso_sp[ipISO][nelem].st[level].n() >= 5 )
121 chType = 's';
122
123 //fixit(); should third parameter go to rfield.nflux?
124 // gjf: don't think so - excited state opacities will be 4-8 dex
125 // smaller than ground state opacity due to population differences
126 // between ground and excited states. H and He like species
127 // have the n=2 level at ~3/4 of the ionization potential, so the
128 // level populations are small compared with ground for recombination
129 // or collision dominated plasmas. This is not true for non-H-like species
130 // so should be revisited if the iso sequence treatment is ever extended
131 // beyond He-like. This is a significant time sink due to large number
132 // of opacity cells across the Lyman coninuum.
133 /* include correction for stimulated emission */
135 iso_sp[ipISO][nelem].fb[level].ipOpac,
136 iso_sp[ipISO][nelem].fb[level].ipIsoLevNIonCon,
137 /* the high energy bound of excited states is the
138 * edge of the Lyman continuum */
139 iso_sp[ipISO][nelem].fb[0].ipIsoLevNIonCon,
140 iso_sp[ipISO][nelem].st[level].Pop(),
141 /* departure coef, volitile opacities */
142 iso_sp[ipISO][nelem].fb[level].DepartCoef , chType );
143 }
144 }
145 }
146 }
147 return;
148}
#define ASSERT(exp)
Definition cddefines.h:578
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_dense dense
Definition dense.cpp:24
t_Heavy Heavy
Definition heavy.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_opac opac
Definition opacity.cpp:5
void OpacityAdd1SubshellInduc(long int ipOpac, long int low, long int ihi, double a, double b, char chStat)
void OpacityAdd1Subshell(long int ipOpac, long int ipLowLim, long int ipUpLim, realnum abundance, char chStat)
void OpacityAdd1Element(long int nelem)
t_rfield rfield
Definition rfield.cpp:8