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