cloudy trunk
Loading...
Searching...
No Matches
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
27static char chPLab[NPUNLM][5];
28static long int nLinesEntered;
30static long int ipLine[NPUNLM];
32
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" );
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",
74 }
75
77
78 /* this is total number stored so far */
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" );
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
102void 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 {
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);
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] );
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 );
267 }
268 return;
269}
270
271#define LIMLINE 10
272static 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 },
276 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
277 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
279 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
280 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
282 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
283 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
285 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
286 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN };
287static bool lgMustPrintHeader=true;
288static long int nLine=-1;
289
290/*Save_Line_RT parse the save line rt command - read in a set of lines */
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" );
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");
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);
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);
348 }
349}
350
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;
362 {
363 fprintf( ioPUN , "Line\tP(con,inc)\tAul\tgl\tgu\n");
364 for( n=0; n<nLine; ++n )
365 {
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");
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 {
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
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define ASSERT(exp)
Definition cddefines.h:578
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long int cdEmis(char *chLabel, realnum wavelength, double *emiss)
Definition cddrive.cpp:533
void cdEmis_ip(long int ipLine, double *emiss, bool lgEmergent)
Definition cddrive.cpp:613
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition cddrive.cpp:1228
realnum ColUL(const ColliderList &colls) const
Definition collision.h:99
double & PopOpc() const
Definition emission.h:603
realnum & Aul() const
Definition emission.h:613
realnum & TauIn() const
Definition emission.h:423
double & pump() const
Definition emission.h:473
bool getline(void)
Definition parser.cpp:164
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
int strcmp(const char *s2)
Definition parser.h:177
void getLineID(char *LabelBuf, realnum *wave)
Definition parser.cpp:446
bool lgEOL(void) const
Definition parser.h:98
bool m_lgEOF
Definition parser.h:42
int PrintLine(FILE *fp) const
Definition parser.h:204
CollisionProxy Coll() const
Definition transition.h:424
qList::iterator Lo() const
Definition transition.h:392
long & ipCont() const
Definition transition.h:450
qList::iterator Hi() const
Definition transition.h:396
EmissionList::reference Emis() const
Definition transition.h:408
ColliderList colliders
Definition collision.cpp:57
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_LineSave LineSave
Definition lines.cpp:5
static realnum * wavelength
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
void sprt_wl(char *chString, realnum wl)
Definition prt.cpp:25
static long int * ipLine
t_radius radius
Definition radius.cpp:5
#define NPUNLM
void Parse_Save_Line_RT(Parser &p)
static long int line_RT_ipHi[LIMLINE]
static long int line_RT_ipLo[LIMLINE]
void parse_save_line(Parser &p, bool lgLog3, char *chHeader)
Definition save_line.cpp:33
static long int nLinesEntered
Definition save_line.cpp:28
static bool lgMustPrintHeader
static char chPLab[NPUNLM][5]
Definition save_line.cpp:27
void Save_Line_RT(FILE *ioPUN)
static long int line_RT_nelem[LIMLINE]
static long int line_RT_ipISO[LIMLINE]
#define LIMLINE
static long int line_RT_type[LIMLINE]
void save_line(FILE *ioPUN, const char *chDo, bool lgEmergent)
static bool lgRelativeIntensity
Definition save_line.cpp:31
static long int nLine
char * chLineLbl(const TransitionProxy &t)