cloudy  trunk
 All Data Structures 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-2008 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 /* faintest pressure we will bother with */
13 #define THRESH 0.05
14 
15 void PrtLinePres(void)
16 {
17  long int i,
18  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 # define NLINE 100
27  /* labels, wavelengths, and fraction of total pressure, for important lines */
28  char chLab[NLINE][5];
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(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.numLevels_local[ipISO][nelem] - iso.nCollapsed_local[ipISO][nelem]; ipHi++ )
64  {
65  for( ipLo=0; ipLo < ipHi; ipLo++ )
66  {
67  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
68  continue;
69 
70  ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > iso.SmallA );
71 
72  /*NB - this code must be kept parallel with code in pressure_total */
73  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc > smallfloat &&
74  /* >>chng 01 nov 01, add test that have not overrun optical depth scale */
75  ( (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn) > smallfloat ) )
76  {
77  RadPres1 = PressureRadiationLine( &Transitions[ipISO][nelem][ipHi][ipLo], dense.xIonDense[nelem][nelem+1-ipISO] );
78 
80  {
81  wl[ip] = Transitions[ipISO][nelem][ipHi][ipLo].WLAng;
82 
83  /* put null terminated line label into chLab using line array*/
84  chIonLbl(chLab[ip], &Transitions[ipISO][nelem][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  /* line radiation pressure from large set of level 1 lines */
97  for( i=1; i <= nLevel1; i++ )
98  {
99  if( TauLines[i].Hi->Pop > 1e-30 )
100  {
101  RadPres1 = PressureRadiationLine( &TauLines[i], 1. );
102  }
103  else
104  {
105  RadPres1 = 0.;
106  }
107 
108  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
109  {
110  wl[ip] = TauLines[i].WLAng;
111 
112  /* put null terminated line label into chLab using line array*/
113  chIonLbl(chLab[ip], &TauLines[i]);
114  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
115 
116  ip = MIN2((long)NLINE-1,ip+1);
117  }
118  }
119 
120  for( i=0; i < nWindLine; i++ )
121  {
122  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
123  {
124  if( TauLine2[i].Hi->Pop > 1e-30 )
125  {
126  RadPres1 = PressureRadiationLine( &TauLine2[i], 1. );
127  }
128  else
129  {
130  RadPres1 = 0.;
131  }
132 
133  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
134  {
135  wl[ip] = TauLine2[i].WLAng;
136 
137  /* put null terminated line label into chLab using line array*/
138  chIonLbl(chLab[ip], &TauLine2[i]);
139  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
140 
141  ip = MIN2((long)NLINE-1,ip+1);
142  }
143  }
144  }
145 
146  for( i=0; i < nHFLines; i++ )
147  {
148  if( HFLines[i].Hi->Pop > 1e-30 )
149  {
150  RadPres1 = PressureRadiationLine( &HFLines[i], 1. );
151  }
152  else
153  {
154  RadPres1 = 0.;
155  }
156 
157  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
158  {
159  wl[ip] = HFLines[i].WLAng;
160 
161  /* put null terminated line label into chLab using line array*/
162  chIonLbl(chLab[ip], &HFLines[i]);
163  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
164 
165  ip = MIN2((long)NLINE-1,ip+1);
166  }
167  }
168 
169  /*Chiani & Leiden Atoms & Molecules*/
170  for( i=0; i <linesAdded2; i++ )
171  {
172  if( atmolEmis[i].tran->Hi->Pop > 1e-30 )
173  {
174  RadPres1 = PressureRadiationLine( atmolEmis[i].tran, 1. );
175  }
176  else
177  {
178  RadPres1 = 0.;
179  }
180 
181  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
182  {
183  wl[ip] = atmolEmis[i].tran->WLAng;
184 
185  /* put null terminated line label into chLab using line array*/
186  chIonLbl(chLab[ip], atmolEmis[i].tran);
187  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
188 
189  ip = MIN2((long)NLINE-1,ip+1);
190  }
191  }
192 
193  /* radiation pressure due to H2 */
194  RadPres1 = H2_RadPress();
195  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
196  {
197  wl[ip] = 0;
198 
199  /* put null terminated 4 char line label into chLab using line array*/
200  strcpy(chLab[ip], "H2 ");
201  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
202 
203  ip = MIN2((long)NLINE-1,ip+1);
204  }
205 
206  for( i=0; i < nCORotate; i++ )
207  {
208  if( C12O16Rotate[i].Hi->Pop > 1e-30 )
209  {
210  RadPres1 = PressureRadiationLine( &C12O16Rotate[i], 1. );
211  }
212  else
213  {
214  RadPres1 = 0.;
215  }
216 
217  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
218  {
219  wl[ip] = C12O16Rotate[i].WLAng;
220 
221  /* put null terminated line label into chLab using line array*/
222  chIonLbl(chLab[ip], &C12O16Rotate[i]);
223  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
224 
225  ip = MIN2((long)NLINE-1,ip+1);
226  }
227  }
228 
229  for( i=0; i < nCORotate; i++ )
230  {
231  if( C13O16Rotate[i].Hi->Pop > 1e-30 )
232  {
233  RadPres1 = PressureRadiationLine( &C13O16Rotate[i], 1. );
234  }
235  else
236  {
237  RadPres1 = 0.;
238  }
239 
240  if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
241  {
242  wl[ip] = C13O16Rotate[i].WLAng;
243 
244  /* put null terminated line label into chLab using line array*/
245  chIonLbl(chLab[ip], &C13O16Rotate[i]);
246  frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
247 
248  ip = MIN2((long)NLINE-1,ip+1);
249  }
250  }
251 
252  /* return if no significant contributors to radiation pressure found */
253  if( ip<= 0 )
254  return;
255 
256  /* sort from strong to weak, then only print strongest */
257  spsort(
258  /* input array to be sorted */
259  frac,
260  /* number of values in x */
261  ip,
262  /* permutation output array */
263  iperm,
264  /* flag saying what to do - 1 sorts into increasing order, not changing
265  * the original routine */
266  -1,
267  /* error condition, should be 0 */
268  &ier);
269 
270  /* now print up to ten strongest contributors to radiation pressure */
271  fprintf( ioQQQ, " P(Lines):" );
272  for( i=0; i < MIN2(10,ip); i++ )
273  {
274  int ipline = iperm[i];
275  fprintf( ioQQQ, "(%4.4s ", chLab[ipline]);
276  prt_wl(ioQQQ, wl[ipline] );
277  fprintf(ioQQQ," %.2f) ",frac[ipline]);
278  }
279 
280  /* finally end the line */
281  fprintf( ioQQQ, "\n" );
282  }
283  return;
284 }

Generated for cloudy by doxygen 1.8.1.2