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