cloudy  trunk
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 
16 void 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 }
dense
t_dense dense
Definition: dense.cpp:24
nLevel1
long int nLevel1
Definition: taulines.cpp:28
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
diatom_iter
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
realnum
float realnum
Definition: cddefines.h:103
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
nWindLine
long nWindLine
Definition: cdinit.cpp:19
PressureRadiationLine
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition: pressure.h:18
mole.h
spsort
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition: service.cpp:1100
diatoms
vector< diatomics * > diatoms
Definition: h2.cpp:8
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
ProxyIterator
Definition: proxy_iterator.h:58
lines_service.h
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
MIN2
#define MIN2
Definition: cddefines.h:761
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
dense.h
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
cddefines.h
chIonLbl
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)
Definition: transition.cpp:195
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
pressure.h
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:101
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
TauLines
TransitionList TauLines("TauLines", &AnonStates)
prt.h
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
THRESH
#define THRESH
Definition: prt_linepres.cpp:14
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:499
nLine
static long int nLine
Definition: save_line.cpp:288
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
nHFLines
long int nHFLines
Definition: taulines.cpp:31
pressure
t_pressure pressure
Definition: pressure.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
h2.h
nSpecies
long int nSpecies
Definition: taulines.cpp:21
NISO
const int NISO
Definition: cddefines.h:261
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
HFLines
TransitionList HFLines("HFLines", &AnonStates)
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
Generated by   doxygen 1.8.17