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

Generated for cloudy by doxygen 1.8.4