cloudy  trunk
save_line.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 /*save_line parse save lines command, or actually do the save lines output */
4 /*Save_Line_RT parse the save line rt command - read in a set of lines */
5 #include "cddefines.h"
6 #include "cddrive.h"
7 #include "radius.h"
8 #include "taulines.h"
9 #include "opacity.h"
10 #include "phycon.h"
11 #include "dense.h"
12 #include "lines.h"
13 #include "h2.h"
14 #include "lines_service.h"
15 #include "input.h"
16 #include "prt.h"
17 #include "save.h"
18 #include "iso.h"
19 #include "parser.h"
20 /* this is the limit to the number of emission lines we can store */
21 #define NPUNLM 200L
22 
23 /* implement the save line xxx command. cumulative, structure, and
24  * emissivity all use same code base and variables, so only one can be used
25  * at present */
26 
27 static char chPLab[NPUNLM][5];
28 static long int nLinesEntered;
29 static realnum wavelength[NPUNLM];
30 static long int ipLine[NPUNLM];
31 static bool lgRelativeIntensity;
32 
33 void parse_save_line(Parser &p,
34  /* true, return rel intensity, false, log of luminosity or intensity I */
35  bool lgLog3,
36  char *chHeader)
37 {
38  char chTemp[INPUT_LINE_LENGTH];
39 
40  // save return value of cdLine, 0 for success, -number of lines for fail
41  long int i;
42 
43  DEBUG_ENTRY( "parse_save_line()" );
44 
45  /* very first time this routine is called, chDo is "READ" and we read
46  * in lines from the input stream. The line labels and wavelengths
47  * are store locally, and output in later calls to this routine
48  * following is flag saying whether to do relative intensity or
49  * absolute emissivity */
50  lgRelativeIntensity = lgLog3;
51 
52  /* number of lines we will save */
53  nLinesEntered = 0;
54 
55  /* get the next line, and check for eof */
56  p.getline();
57  if( p.m_lgEOF )
58  {
59  fprintf( ioQQQ,
60  " Hit EOF while reading line list; use END to end list.\n" );
61  cdEXIT(EXIT_FAILURE);
62  }
63 
64  /* convert line to caps */
65 
66  while( p.strcmp("END") != 0 )
67  {
68  if( nLinesEntered >= NPUNLM )
69  {
70  fprintf( ioQQQ,
71  " Too many lines have been entered; the limit is %ld. Increase variable NPUNLM in routine save_line.\n",
72  nLinesEntered );
73  cdEXIT(EXIT_FAILURE);
74  }
75 
76  p.getLineID(chPLab[nLinesEntered], &wavelength[nLinesEntered]);
77 
78  /* this is total number stored so far */
79  ++nLinesEntered;
80 
81  /* get next line and check for eof */
82  p.getline();
83  if( p.m_lgEOF )
84  {
85  fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
86  cdEXIT(EXIT_FAILURE);
87  }
88 
89  }
90 
91  sprintf( chHeader, "depth");
92  for( i=0; i < nLinesEntered; i++ )
93  {
94  sprintf( chTemp, "\t%s ", chPLab[i] );
95  strcat( chHeader, chTemp );
96  sprt_wl( chTemp, wavelength[i] );
97  strcat( chHeader, chTemp );
98  }
99  strcat( chHeader, "\n" );
100 }
101 
102 void save_line(FILE * ioPUN, /* the file we will write to */
103  const char *chDo,
104  // intrinsic or emergent line emission?
105  bool lgEmergent)
106 {
107  // save return value of cdLine, 0 for success, -number of lines for fail
108  long int nCdLineReturn;
109  long int i;
110  double a[32],
111  absint,
112  emiss,
113  relint;
114  double dlum[NPUNLM];
115 
116  DEBUG_ENTRY( "save_line()" );
117 
118  if( strcmp(chDo,"PUNS") == 0 )
119  {
120  /* save lines emissivity command */
121  static bool lgMustGetLines=true,
122  lgBadLine=false;
123  lgBadLine = false;
124  static bool lgBadH2Line;
125  /* it is possible that we will get here after an initial temperature
126  * too high abort, and the line arrays will not have been defined.
127  * do no lines in this case. must still do save so that there
128  * is not a missing line in the grid save output */
129  if( LineSave.nsum >0 )
130  {
131  lgBadH2Line = false;
132  lgBadLine = false;
133  /* save lines structure command */
134  for( i=0; i < nLinesEntered; i++ )
135  {
136  if( nzone <= 1 && lgMustGetLines )
137  {
138  if( (ipLine[i] = cdEmis((char*)chPLab[i],wavelength[i],
139  &emiss,lgEmergent)) <= 0 )
140  {
141  // missed line - ignore if H2 line since large model may be off
142  if( !h2.lgEnabled && strncmp( chPLab[i] , "H2 " , 4 )==0 )
143  {
144  static bool lgMustPrintFirstTime = true;
145  if( lgMustPrintFirstTime )
146  {
147  /* it's an H2 line and H2 is not being done - ignore it */
148  fprintf( ioQQQ,"\nPROBLEM Did not find an H2 line, the large model is not "
149  "included, so I will ignore it. Log intensity set to -30.\n" );
150  fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n\n");
151  lgMustPrintFirstTime = false;
152  }
153  /* flag saying to ignore this line */
154  ipLine[i] = -2;
155  lgBadH2Line = true;
156  }
157  else
158  {
159  fprintf( ioQQQ, " save_line could not find line: %s %f\n",
160  chPLab[i], wavelength[i] );
161  ipLine[i] = -1;
162  lgBadLine = true;
163  }
164  }
165  }
166  /* 0th line is dummy, can't be used, so this is safe */
167  ASSERT( ipLine[i] > 0 || lgBadLine || lgBadH2Line ||
168  /* this is case where we did not find line on previous iteration,
169  * perhaps because H2 is not turned on, and -1 or -2 was
170  * stored */
171  (ipLine[i]<0&&!lgMustGetLines) );
172  /* this version of cdEmis uses index, does not search, do not call if line could not be found */
173  /* test on case where we abort before first zone is done
174  * this happens in grid when temperature bounds of code
175  * are exceeded. In this case return small float */
176  if( lgAbort && nzone <=1 )
177  dlum[i] = SMALLFLOAT;
178  else if( ipLine[i]>0 )
179  cdEmis_ip(ipLine[i],&dlum[i],lgEmergent);
180  else
181  dlum[i] = 1e-30f;
182  }
183  if( lgBadLine )
184  {
185  cdEXIT(EXIT_FAILURE);
186  }
187  }
188  lgMustGetLines = false;
189 
190  /* print depth */
191  fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
192 
193  /* then print all line emissivity */
194  for( i=0; i < nLinesEntered; i++ )
195  {
196  /*lint -e644 dlum not initialized */
197  fprintf( ioPUN, "\t%.4f", log10( MAX2( SMALLFLOAT , dlum[i] ) ) );
198  /*lint +e644 dlum not initialized */
199  }
200  fprintf( ioPUN, "\n" );
201  }
202 
203  else if( strcmp(chDo,"PUNC") == 0 )
204  {
205  /* save lines cumulative command */
206  fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
207 
208  /* it is possible that we will get here after an initial temperature
209  * too high abort, and the line arrays will not have been defined.
210  * do no lines in this case. must still do save so that there
211  * is not a missing line in the grid save output */
212  if( LineSave.nsum >0 )
213  {
214  for( i=0; i < nLinesEntered; i++ )
215  {
216  nCdLineReturn = cdLine((char*)chPLab[i],wavelength[i],
217  &relint,&absint,lgEmergent);
218  if( lgRelativeIntensity )
219  /* relative intensity case */
220  a[i] = relint;
221  else
222  /* emissivity or luminosity case */
223  a[i] = absint;
224 
225  if( nCdLineReturn<=0 )
226  {
227  /* missed line - ignore if H2 line */
228  if( !h2.lgEnabled && strncmp( chPLab[i] , "H2 " , 4 )==0 )
229  {
230  static bool lgMustPrintFirstTime = true;
231  if( lgMustPrintFirstTime )
232  {
233  /* it's an H2 line and H2 is not being done - ignore it */
234  fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
235  "included, so I will ignore it. Log intensity set to -30.\n" );
236  fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
237  lgMustPrintFirstTime = false;
238  }
239  /* flag saying to ignore this line */
240  a[i] = -30.;
241  absint = -30.;
242  relint = -30.;
243  }
244  else
245  {
246  fprintf( ioQQQ, " save_line could not fine line: %s %f\n",
247  chPLab[i], wavelength[i] );
248  cdEXIT(EXIT_FAILURE);
249  }
250  }
251  }
252 
253  for( i=0; i < nLinesEntered; i++ )
254  {
255  fprintf( ioPUN, "\t%.4e", a[i] );
256  }
257  }
258  fprintf( ioPUN, "\n" );
259  }
260 
261  else
262  {
263  fprintf( ioQQQ,
264  " unrecognized key for save_line=%4.4s\n",
265  chDo );
266  cdEXIT(EXIT_FAILURE);
267  }
268  return;
269 }
270 
271 #define LIMLINE 10
272 static long int line_RT_type[LIMLINE] =
273  {LONG_MIN , LONG_MIN ,LONG_MIN , LONG_MIN ,LONG_MIN ,
274  LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
275  line_RT_ipISO[LIMLINE] =
276  {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
277  LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
278  line_RT_nelem[LIMLINE] =
279  {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
280  LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
281  line_RT_ipHi[LIMLINE] =
282  {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
283  LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
284  line_RT_ipLo[LIMLINE] =
285  {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
286  LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN };
287 static bool lgMustPrintHeader=true;
288 static long int nLine=-1;
289 
290 /*Save_Line_RT parse the save line rt command - read in a set of lines */
291 void Parse_Save_Line_RT(Parser &p)
292 {
293  /* save line rt parameters */
294  DEBUG_ENTRY( "Parse_Save_Line_RT()" );
295 
296  /* very first time this routine is called, chDo is "READ" and we read
297  * in lines from the input stream. The line labels and wavelengths
298  * are store locally, and output in later calls to this routine */
299 
300  /* say that we must print the header */
301  lgMustPrintHeader = true;
302 
303  /* get the next line, and check for eof */
304  nLine = 0;
305  p.getline();
306  if( p.m_lgEOF )
307  {
308  fprintf( ioQQQ,
309  " Hit EOF while reading line list; use END to end list.\n" );
310  cdEXIT(EXIT_FAILURE);
311  }
312 
313  do
314  {
315  if(nLine>=LIMLINE )
316  {
317  fprintf(ioQQQ," PUNCH RT has too many lines - increase LIMLINE in save_line.cpp\n");
318  cdEXIT(EXIT_FAILURE);
319  }
320 
321  /* right now it just does lines in the iso sequences */
322  line_RT_type[nLine] = (long int)p.FFmtRead();
323  line_RT_ipISO[nLine] = (long int)p.FFmtRead();
324  line_RT_nelem[nLine] = (long int)p.FFmtRead();
325  line_RT_ipHi[nLine] = (long int)p.FFmtRead();
326  line_RT_ipLo[nLine] = (long int)p.FFmtRead();
327 
328  if( p.lgEOL() )
329  {
330  fprintf( ioQQQ,
331  " there must be five numbers on this line\n");
332  p.PrintLine(ioQQQ);
333  cdEXIT(EXIT_FAILURE);
334  }
335 
336  /* now increment number of lines */
337  ++nLine;
338 
339  /* now get next line until we hit eof or the word END */
340  p.getline();
341  } while( !p.m_lgEOF && !p.nMatch( "END") );
342  if( p.m_lgEOF )
343  {
344  fprintf( ioQQQ,
345  " Save_Line_RT hit end of file looking for END of RT lines\n");
346  p.PrintLine(ioQQQ);
347  cdEXIT(EXIT_FAILURE);
348  }
349 }
350 
351 void Save_Line_RT(
352  FILE * ioPUN )
353 {
354  /* save line rt parameters */
355 
356  DEBUG_ENTRY( "Save_Line_RT()" );
357 
358 
359  static char chLabel[LIMLINE][30];
360  long int n;
361  if( lgMustPrintHeader )
362  {
363  fprintf( ioPUN , "Line\tP(con,inc)\tAul\tgl\tgu\n");
364  for( n=0; n<nLine; ++n )
365  {
366  TransitionProxy tr = iso_sp[line_RT_ipISO[n]][line_RT_nelem[n]].trans(line_RT_ipHi[n],line_RT_ipLo[n]);
367  /* print info for header of file, line id and pump rate */
368  sprintf( chLabel[n] , "%s ",
369  chLineLbl(tr) );
370  fprintf( ioPUN , "%s ", chLabel[n] );
371  fprintf( ioPUN , "%.4e ",
372  tr.Emis().pump());
373  fprintf( ioPUN , "%.4e ",
374  tr.Emis().Aul());
375  fprintf( ioPUN , "%.0f ",
376  (*tr.Lo()).g());
377  fprintf( ioPUN , "%.0f ",
378  (*tr.Hi()).g());
379  fprintf( ioPUN , "\n");
380 
381  if( line_RT_type[n]!=0. )
382  {
383  /* for now code only exists for H He like iso - assert this */
384  fprintf( ioQQQ,
385  " PunchLine_RT only H, He like allowed for now\n");
386  cdEXIT(EXIT_FAILURE);
387  }
388  }
389  fprintf( ioPUN , "Line\tTauIn\tPopLo\tPopHi\tCul\tk(line)\tk(con,abs)\tk(con,scat)\n");
390  lgMustPrintHeader = false;
391  }
392 
393  fprintf(ioPUN, "RADIUS\t%e\tDEPTH\t%e\tTe\t%e\tNe\t%e\n",
394  radius.Radius_mid_zone ,
395  radius.depth_mid_zone ,
396  phycon.te ,
397  dense.eden );
398  for( n=0; n<nLine; ++n )
399  {
400  TransitionProxy tr = iso_sp[line_RT_ipISO[n]][line_RT_nelem[n]].trans(line_RT_ipHi[n],line_RT_ipLo[n]);
401 
402  /* index for line within continuum array */
403  long int ipCont = tr.ipCont();
404  fprintf( ioPUN , "%s ", chLabel[n] );
405  fprintf( ioPUN , "\t%e\t%e\t%e",
406  tr.Emis().TauIn() ,
407  (*tr.Lo()).Pop(),
408  (*tr.Hi()).Pop()
409  );
410  fprintf( ioPUN , "\t%e",
411  tr.Coll().ColUL( colliders ) / dense.EdenHCorr
412  );
413 
414  fprintf( ioPUN , "\t%e\t%e\t%e\n",
415  tr.Emis().PopOpc(),
416  opac.opacity_abs[ipCont-1] ,
417  opac.opacity_sct[ipCont-1]
418  );
419  }
420 }
421 
422 # undef LIMELM
423 
chLineLbl
char * chLineLbl(const TransitionProxy &t)
Definition: transition.cpp:237
parse_save_line
void parse_save_line(Parser &p, bool lgLog3, char *chHeader)
Definition: save_line.cpp:33
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
sprt_wl
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:25
lines.h
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:353
t_dense::eden
double eden
Definition: dense.h:190
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
colliders
ColliderList colliders
Definition: collision.cpp:57
dense
t_dense dense
Definition: dense.cpp:24
Parse_Save_Line_RT
void Parse_Save_Line_RT(Parser &p)
Definition: save_line.cpp:291
lgMustPrintHeader
static bool lgMustPrintHeader
Definition: save_line.cpp:287
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:345
realnum
float realnum
Definition: cddefines.h:103
lgRelativeIntensity
static bool lgRelativeIntensity
Definition: save_line.cpp:31
Save_Line_RT
void Save_Line_RT(FILE *ioPUN)
Definition: save_line.cpp:351
EmissionProxy::PopOpc
double & PopOpc() const
Definition: emission.h:603
line_RT_ipHi
static long int line_RT_ipHi[LIMLINE]
Definition: save_line.cpp:281
Parser::m_lgEOF
bool m_lgEOF
Definition: parser.h:42
phycon
t_phycon phycon
Definition: phycon.cpp:6
lines_service.h
t_LineSave::nsum
long int nsum
Definition: lines.h:62
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:424
opac
t_opac opac
Definition: opacity.cpp:5
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
TransitionProxy
Definition: transition.h:23
cddrive.h
nzone
long int nzone
Definition: cddefines.cpp:14
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
radius
t_radius radius
Definition: radius.cpp:5
nLinesEntered
static long int nLinesEntered
Definition: save_line.cpp:28
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
chPLab
static char chPLab[NPUNLM][5]
Definition: save_line.cpp:27
Parser
Definition: parser.h:31
dense.h
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
Parser::getline
bool getline(void)
Definition: parser.cpp:164
LIMLINE
#define LIMLINE
Definition: save_line.cpp:271
cddefines.h
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:98
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
cdEmis_ip
void cdEmis_ip(long int ipLine, double *emiss, bool lgEmergent)
Definition: cddrive.cpp:613
radius.h
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
EmissionProxy::pump
double & pump() const
Definition: emission.h:473
MAX2
#define MAX2
Definition: cddefines.h:782
CollisionProxy::ColUL
realnum ColUL(const ColliderList &colls) const
Definition: collision.h:99
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:98
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
save.h
prt.h
line_RT_nelem
static long int line_RT_nelem[LIMLINE]
Definition: save_line.cpp:278
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
wavelength
static realnum wavelength[NPUNLM]
Definition: save_line.cpp:29
cdEmis
long int cdEmis(char *chLabel, realnum wavelength, double *emiss)
Definition: cddrive.cpp:533
parser.h
nLine
static long int nLine
Definition: save_line.cpp:288
Parser::getLineID
void getLineID(char *LabelBuf, realnum *wave)
Definition: parser.cpp:446
line_RT_ipLo
static long int line_RT_ipLo[LIMLINE]
Definition: save_line.cpp:284
line_RT_ipISO
static long int line_RT_ipISO[LIMLINE]
Definition: save_line.cpp:275
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
NPUNLM
#define NPUNLM
Definition: save_line.cpp:21
line_RT_type
static long int line_RT_type[LIMLINE]
Definition: save_line.cpp:272
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:41
phycon.h
ipLine
static long int ipLine[NPUNLM]
Definition: save_line.cpp:30
save_line
void save_line(FILE *ioPUN, const char *chDo, bool lgEmergent)
Definition: save_line.cpp:102
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
Parser::strcmp
int strcmp(const char *s2)
Definition: parser.h:177
h2.h
Parser::PrintLine
int PrintLine(FILE *fp) const
Definition: parser.h:204
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_phycon::te
double te
Definition: phycon.h:11
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:28
input.h
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
Generated by   doxygen 1.8.17