cloudy trunk
Loading...
Searching...
No Matches
prt_linepres.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/*PrtLinePres print line radiation pressures for current conditions */
4#include "cddefines.h"
5#include "pressure.h"
6#include "taulines.h"
7#include "iso.h"
8#include "dense.h"
9#include "lines_service.h"
10#include "h2.h"
11#include "prt.h"
12#include "mole.h"
13/* faintest pressure we will bother with */
14#define THRESH 0.05
15
16void PrtLinePres( FILE *ioPRESSURE )
17{
18 long int i,
19 ip,
20 ipLo,
21 ipHi,
22 nelem;
23 int ier;
24 double RadPres1;
25
26 // this is limit to number of lines we can store at one time
27 const int nLine = 100;
28 /* labels, wavelengths, and fraction of total pressure, for important lines */
29 char chLab[nLine][5];
30 realnum wl[nLine] , frac[nLine];
31 long int iperm[nLine];
32
33 /* will be used to check on size of opacity, was capped at this value */
34 realnum smallfloat=SMALLFLOAT*10.f;
35
36 DEBUG_ENTRY( "PrtLinePres()" );
37
38 /* this routine only called if printout of contributors to line
39 * radiation pressure is desired */
40
41 /* compute radiation pressure due to iso lines */
42 ip = 0;
43
44 for(i=0; i<nLine; ++i)
45 {
46 frac[i] = FLT_MAX;
47 wl[i] = FLT_MAX;
48 }
49
50 if( pressure.pres_radiation_lines_curr > 1e-30 )
51 {
55 for( long ipISO = 0; ipISO<NISO; ipISO++ )
56 {
57 for( nelem=ipISO; nelem < LIMELM; nelem++ )
58 {
59 /* does this ion stage exist? */
60 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
61 {
62 /* do not include highest levels since maser can occur due to topoff,
63 * and pops are set to small number in this case */
64 for( ipHi=1; ipHi <iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
65 {
66 for( ipLo=0; ipLo < ipHi; ipLo++ )
67 {
68 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
69 continue;
70
71 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > iso_ctrl.SmallA );
72
73 /*NB - this code must be kept parallel with code in pressure_total */
74 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() > smallfloat &&
75 /* >>chng 01 nov 01, add test that have not overrun optical depth scale */
76 ( (iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() - iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()) > smallfloat ) )
77 {
78 RadPres1 = PressureRadiationLine( iso_sp[ipISO][nelem].trans(ipHi,ipLo), GetDopplerWidth(dense.AtomicWeight[nelem]) );
79
80 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
81 {
82 wl[ip] = iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng();
83
84 /* put null terminated line label into chLab using line array*/
85 chIonLbl(chLab[ip], iso_sp[ipISO][nelem].trans(ipHi,ipLo));
86 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
87
88 ip = MIN2((long)nLine-1,ip+1);
89 }
90 }
91 }
92 }
93 }
94 }
95 }
96
97 /* line radiation pressure from large set of level 1 lines */
98 for( i=1; i <= nLevel1; i++ )
99 {
100 if( (*TauLines[i].Hi()).Pop() > 1e-30 )
101 {
102 RadPres1 = PressureRadiationLine( TauLines[i], GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
103 }
104 else
105 {
106 RadPres1 = 0.;
107 }
108
109 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
110 {
111 wl[ip] = TauLines[i].WLAng();
112
113 /* put null terminated line label into chLab using line array*/
114 chIonLbl(chLab[ip], TauLines[i]);
115 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
116
117 ip = MIN2((long)nLine-1,ip+1);
118 }
119 }
120
121 for( i=0; i < nWindLine; i++ )
122 {
123 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
124 {
125 if( (*TauLine2[i].Hi()).Pop() > 1e-30 )
126 {
127 RadPres1 = PressureRadiationLine( TauLine2[i], GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
128 }
129 else
130 {
131 RadPres1 = 0.;
132 }
133
134 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
135 {
136 wl[ip] = TauLine2[i].WLAng();
137
138 /* put null terminated line label into chLab using line array*/
139 chIonLbl(chLab[ip], TauLine2[i]);
140 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
141
142 ip = MIN2((long)nLine-1,ip+1);
143 }
144 }
145 }
146
147 for( i=0; i < nHFLines; i++ )
148 {
149 if( (*HFLines[i].Hi()).Pop() > 1e-30 )
150 {
151 RadPres1 = PressureRadiationLine( HFLines[i], GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
152 }
153 else
154 {
155 RadPres1 = 0.;
156 }
157
158 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
159 {
160 wl[ip] = HFLines[i].WLAng();
161
162 /* put null terminated line label into chLab using line array*/
163 chIonLbl(chLab[ip], HFLines[i]);
164 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
165
166 ip = MIN2((long)nLine-1,ip+1);
167 }
168 }
169
170 /* lines from external databases */
171 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
172 {
173 if( dBaseSpecies[ipSpecies].lgActive )
174 {
175 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
176 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
177 tr != dBaseTrans[ipSpecies].end(); ++tr)
178 {
179 int ipHi = (*tr).ipHi();
180 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
181 continue;
182 if( (*(*tr).Hi()).Pop() > 1e-30 )
183 {
184 RadPres1 = PressureRadiationLine( *tr, DopplerWidth );
185 }
186 else
187 {
188 RadPres1 = 0.;
189 }
190
191 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
192 {
193 wl[ip] = (*tr).WLAng();
194
195 /* put null terminated line label into chLab using line array*/
196 chIonLbl(chLab[ip], *tr );
197 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
198
199 ip = MIN2((long)nLine-1,ip+1);
200 }
201 }
202 }
203 }
204
205 /* radiation pressure due to H2 */
206 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
207 {
208 RadPres1 = (*diatom)->H2_RadPress();
209 if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
210 {
211 wl[ip] = 0;
212
213 /* put null terminated 4 char line label into chLab using line array*/
214 strcpy(chLab[ip], "H2 ");
215 frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
216
217 ip = MIN2((long)nLine-1,ip+1);
218 }
219 }
220
221 /* return if no significant contributors to radiation pressure found */
222 if( ip<= 0 )
223 {
224 fprintf( ioPRESSURE, "\n" );
225 return;
226 }
227
228 /* sort from strong to weak, then only print strongest */
229 spsort(
230 /* input array to be sorted */
231 frac,
232 /* number of values in x */
233 ip,
234 /* permutation output array */
235 iperm,
236 /* flag saying what to do - 1 sorts into increasing order, not changing
237 * the original routine */
238 -1,
239 /* error condition, should be 0 */
240 &ier);
241
242 /* now print up to ten strongest contributors to radiation pressure */
243 fprintf( ioPRESSURE, " P(Lines):" );
244 for( i=0; i < MIN2(10,ip); i++ )
245 {
246 int ipline = iperm[i];
247 fprintf( ioPRESSURE, "(%4.4s ", chLab[ipline]);
248 prt_wl(ioPRESSURE, wl[ipline] );
249 fprintf(ioPRESSURE," %.2f) ",frac[ipline]);
250 }
251
252 /* finally end the line */
253 fprintf( ioPRESSURE, "\n" );
254 }
255 return;
256}
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
TransitionProxy::iterator iterator
Definition transition.h:280
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
t_pressure pressure
Definition pressure.cpp:5
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition pressure.h:18
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
#define THRESH
void PrtLinePres(FILE *ioPRESSURE)
static long int nLine
long int nSpecies
Definition taulines.cpp:21
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
species * dBaseSpecies
Definition taulines.cpp:14
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)