cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_punch.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 /*ParsePunch parse the punch command */
4 /*PunchFilesInit initialize punch file pointers, called from cdInit */
5 /*ClosePunchFiles close all punch files */
6 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present */
7 #include "cddefines.h"
8 #include "cddrive.h"
9 #include "physconst.h"
10 #include "elementnames.h"
11 #include "input.h"
12 #include "geometry.h"
13 #include "prt.h"
14 #include "optimize.h"
15 #include "rfield.h"
16 #include "hcmap.h"
17 #include "atomfeii.h"
18 #include "h2.h"
19 #include "mole.h"
20 #include "hmi.h"
21 #include "version.h"
22 #include "grainvar.h"
23 #include "parse.h"
24 #include "grid.h"
25 #include "punch.h"
26 /* check for keyword UNITS on line, then scan wavelength or energy units if present */
27 STATIC void ChkUnits( char *chCard);
28 
29 /* option to append instead of overwrite - default was to overwrite,
30  * >>chng 07 feb 02 invert logic, default not to overwrite, but clobber
31  * keyword says to overwrite */
32 static bool lgNoClobber[LIMPUN];
33 
34 /* these are for some special cases, same purpose as previous no clobber */
38 
39 /* NB NB NB NB NB NB NB NB NB NB
40  *
41  * if any "special" punch commands are added (commands that copy the file pointer
42  * into a separate variable, e.g. like PUNCH _DR_), be sure to add that file pointer
43  * to PunchFilesInit and ClosePunchFiles !!!
44  *
45  * PUNCH FILE POINTERS SHOULD NEVER BE ALTERED BY ROUTINES OUTSIDE THIS MODULE !!!
46  *
47  * hence initializations of punch file pointers should never be included in zero() !!
48  * the pointer might be lost without the file being closed
49  *
50  * NB NB NB NB NB NB NB NB NB NB */
51 
52 /* punch file header headers - these are written into the string chHeader when
53  * the command is parsed
54  * punch.lgPunHeader determines whether header is punched
55  */
56 
57 
58 void ParsePunch(char *chCard )
59 {
61  chFilename[INPUT_LINE_LENGTH] ,
62  chSecondFilename[INPUT_LINE_LENGTH];
63  bool lgEOL,
64  lgSecondFilename;
65  /* pointer to column across line image for free format number reader*/
66  long int ipFFmt,
67  i,
68  nelem;
69 
70 #define MAX_HEADER_SIZE 2000
71 
72  char chHeader[MAX_HEADER_SIZE],
73  chTemp[MAX_HEADER_SIZE];
74 
75 
76  DEBUG_ENTRY( "ParsePunch()" );
77 
78  /* initialize chHeader string with nonsense, compare at end to see if we have any actual headers. */
79  sprintf( chHeader, "ArNdY38dZ9us4N4e12SEcuQ" );
80 
81  /* check that limit not exceeded */
82  if( punch.npunch >= LIMPUN )
83  {
84  fprintf( ioQQQ,
85  "The limit to the number of PUNCH options is %ld. Increase "
86  "LIMPUN in punch.h if more are needed.\nSorry.\n",
87  LIMPUN );
88  cdEXIT(EXIT_FAILURE);
89  }
90 
91  /* if this is first call to this routine we will need to init some vars */
92  {
93  static bool lgMustInit=true;
94  if( lgMustInit )
95  {
96  lgMustInit = false;
97 
98  /* following are for punch LineList option
99  * nLineList is number of em lines, -1 if not defined */
100  punch.nLineList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) );
101  /* lgLineListRatio flag saying to take ratios of pairs of lines */
102  punch.lgLineListRatio = (bool*)MALLOC((size_t)(LIMPUN*sizeof(bool)) );
103  /* wlLineList is set of emission lines for LineList */
104  punch.wlLineList = (realnum **)MALLOC((size_t)(LIMPUN*sizeof(realnum *)) );
105  /* chLineListLabel is label for line list */
106  punch.chLineListLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
107  /* remaining dimensions will be created if punch LineList command entered
108  * but init nLineList to -1 to signal not set */
109  for( i=0; i<LIMPUN; ++i )
110  {
111  punch.nLineList[i] = -1;
112  punch.lgFITS[i] = false;
113  }
114 
115  /* following are parallel but for punch averages option
116  * nAverageList is number of averages, -1 if not defined */
117  punch.nAverageList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) );
118  /* chAverageSpeciesLabel is label for species of each average */
119  punch.chAverageSpeciesLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
120  /* chAverageType is label for species of type of average */
121  punch.chAverageType = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) );
122  /* nAverageIonList is set of ion stages for average */
123  punch.nAverageIonList = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) );
124  /* nAverage2ndPar is second parameters of each average */
125  punch.nAverage2ndPar = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) );
126  /* remaining dimensions will be created if punch averages command entered
127  * but init nLineList to -1 to signal not set */
128  for( i=0; i<LIMPUN; ++i )
129  {
130  punch.nAverageList[i] = -1;
131  }
132  }
133  }
134 
135  /* initialize this flag to false, set true for special cases below */
137 
138  /* LAST keyword is an option to punch only on last iteration */
139  if( nMatch("LAST",chCard) )
140  punch.lgPunLstIter[punch.npunch] = true;
141  else
142  punch.lgPunLstIter[punch.npunch] = false;
143 
144  /* get file name for this punch output.
145  * GetQuote does the following -
146  * first copy original version of file name into chLabel,
147  * string does include null termination.
148  * set filename in OrgCard and second parameter to spaces so
149  * that not picked up below as keyword
150  * last parameter says whether to abort if no quote found */
151  if( GetQuote( chLabel , chCard , true ) )
152  /* this can't happen since routine would not return at all if no double quotes found */
153  TotalInsanity();
154 
155  /* check that name is not same as opacity.opc, a special file */
156  if( strcmp(chLabel , "opacity.opc") == 0 )
157  {
158  fprintf( ioQQQ, "ParsePunch will not allow punch file name %s, please choose another.\nSorry.\n",
159  chLabel);
160  cdEXIT(EXIT_FAILURE);
161  }
162  else if( chLabel[0]==0 )
163  {
164  fprintf( ioQQQ, "ParsePunch found a null file name between double quotes, please check command line.\nSorry.\n");
165  cdEXIT(EXIT_FAILURE);
166  }
167 
168  /* now copy to chFilename, with optional prefix first */
169  /* this is optional prefix, normally a null string, set with set punch prefix command */
170  strcpy( chFilename , punch.chFilenamePrefix );
171  strcat( chFilename , chLabel );
172 
173  /* there may be a second file name, and we need to get it off the line
174  * before we parse options, last false parameter says not to abort if
175  * missing - this is not a problem at this stage */
176  if( GetQuote( chSecondFilename , chCard , false ) )
177  lgSecondFilename = false;
178  else
179  lgSecondFilename = true;
180 
181  /* CLOBBER clobber keyword is an option to overwrite rather than
182  * append to a given file */
183  if( nMatch("CLOB",chCard) )
184  {
185  if( nMatch(" NO ",chCard) )
186  {
187  /* do not clobber files */
188  lgNoClobber[punch.npunch] = true;
189  }
190  else
191  {
192  /* clobber files */
193  lgNoClobber[punch.npunch] = false;
194  }
195  }
196 
197  /* put version number and title of model on output file, but only if
198  * this is requested with a "title" on the line*/
199  /* >>chng 02 may 10, invert logic from before - default had been title */
200  /* put version number and title of model on output file, but only if
201  * there is not a "no title" on the line*/
202  if( !nMatch("O TI",chCard) && nMatch("TITL",chCard))
203  {
204  sprintf( chHeader,
205  "# %s %s\n",
206  t_version::Inst().chVersion, input.chTitle );
207  }
208 
209  /* usually results for each iteration are followed by a series of
210  * hash marks, ####, which fool excel. This is an option to not print
211  * the line. If the keyword NO HASH no hash appears the hash marks
212  * will not occur */
213  if( nMatch("O HA",chCard) )
214  punch.lgHashEndIter[punch.npunch] = false;
215 
216  /* punch opacity must come first since elements appear as sub-keywords */
217  if( nMatch("OPAC",chCard) )
218  {
219  /* for grid run, punch results of different models to different files. */
221 
222  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
223  * units are copied into punch.chConPunEnr */
224  ChkUnits(chCard);
225 
226  strcpy( punch.chPunch[punch.npunch], "OPAC" );
227  if( nMatch("TOTA",chCard) )
228  {
229  /* DoPunch will call punch_opacity to parse the subcommands
230  * punch total opacity */
231  strcpy( punch.chOpcTyp[punch.npunch], "TOTL" );
232  sprintf( chHeader,
233  "#nu\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n" );
234  }
235 
236  else if( nMatch("FIGU",chCard) )
237  {
238  /* do figure for hazy */
239  strcpy( punch.chOpcTyp[punch.npunch], "FIGU" );
240  sprintf( chHeader,
241  "#nu, H, He, tot opac\n" );
242  }
243 
244  else if( nMatch("FINE",chCard) )
245  {
246  /* punch the fine opacity array */
247  rfield.lgPunchOpacityFine = true;
248  strcpy( punch.chOpcTyp[punch.npunch], "FINE" );
249  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
250  * units are copied into punch.chConPunEnr */
251  ChkUnits(chCard);
252 
253  sprintf( chHeader,
254  "#nu\topac\n" );
255  ipFFmt = 5;
256  /* range option - important since so much data - usually want to
257  * only give portion of the continuum */
258  if( nMatch("RANGE",chCard) )
259  {
260  /* get lower and upper range, must be in Ryd */
261  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
262  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
263  if( lgEOL )
264  {
265  fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
266  cdEXIT(EXIT_FAILURE);
267  }
268  if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
269  {
270  fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
271  cdEXIT(EXIT_FAILURE);
272  }
273  }
274  else
275  {
276  /* these mean full energy range */
277  punch.punarg[punch.npunch][0] = 0.;
278  punch.punarg[punch.npunch][1] = 0.;
279  }
280  /* optional last parameter - how many points to bring together */
281  punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
282  INPUT_LINE_LENGTH,&lgEOL);
283 
284  /* default is to average together ten */
285  if( lgEOL )
286  punch.punarg[punch.npunch][2] = 10;
287 
288  if( punch.punarg[punch.npunch][2] < 1 )
289  {
290  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
291  cdEXIT(EXIT_FAILURE);
292  }
293  }
294 
295  else if( nMatch("GRAI",chCard) )
296  {
297  /* punch grain opacity command, give optical properties of gains in calculation */
298  strcpy( punch.chPunch[punch.npunch], "DUSO" );
299  /* punch grain opacity command in twice, here and above in opacity */
300  sprintf( chHeader,
301  "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
302  }
303 
304  else if( nMatch("BREM",chCard) )
305  {
306  /* punch bremsstrahlung opacity */
307  strcpy( punch.chOpcTyp[punch.npunch], "BREM" );
308  sprintf( chHeader,
309  "#nu\tbrem opac\n" );
310  }
311 
312  else if( nMatch("SHEL",chCard) )
313  {
314  /* punch shells, a form of the punch opacity command for showing subshell crossections*/
315  strcpy( punch.chPunch[punch.npunch], "OPAC" );
316 
317  /* punch subshell cross sections */
318  strcpy( punch.chOpcTyp[punch.npunch], "SHEL" );
319 
320  /* this is element */
321  ipFFmt = 3;
322  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,
323  INPUT_LINE_LENGTH,&lgEOL);
324 
325  /* this is ion */
326  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,
327  INPUT_LINE_LENGTH,&lgEOL);
328 
329  /* this is shell */
330  punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
331  INPUT_LINE_LENGTH,&lgEOL);
332 
333  if( lgEOL )
334  {
335  fprintf( ioQQQ, "There must be IO unit, atom weight, ion, shell\nSorry.\n" );
336  cdEXIT(EXIT_FAILURE);
337  }
338  sprintf( chHeader,
339  "#sub shell cross section\n" );
340  }
341 
342  else if( nMatch("ELEM",chCard) )
343  {
344  /* punch element opacity, produces n name.n files, one for each stage of
345  * ionization. the name is the 4-char version of the element's name, and
346  * n is the stage of ionization. the file name on the card is ignored.
347  * The code stops in punch_opacity after these files are produced. */
348 
349  /* this will be used as check that we did find an element on the command lines */
350  /* nelem is -1 if an element was not found */
351  if( (nelem = GetElem( chCard ) ) < 0 )
352  {
353  fprintf( ioQQQ, "There must be a second keyword on the opacity command. Sorry.\n" );
354  fprintf( ioQQQ, "Options are total, figure, grain, or element name.\n" );
355  cdEXIT(EXIT_FAILURE);
356  }
357 
358  /* copy string over */
360  }
361  else
362  {
363  fprintf( ioQQQ, " I did not recognize a keyword on this punch opacity command.\n" );
364  fprintf( ioQQQ, " Sorry.\n" );
365  cdEXIT(EXIT_FAILURE);
366  }
367  }
368 
369  /* punch H2 has to come early since it has many suboptions */
370  else if( nMatch(" H2 ",chCard) )
371  {
372  /* this is in mole_h2_io.c */
373  H2_ParsePunch( chCard , chHeader );
374  }
375 
376  /* punch grain abundance will be handled later */
377  else if( nMatch("ABUN",chCard) && !nMatch("GRAI",chCard) )
378  {
379  /* punch abundances */
380  strcpy( punch.chPunch[punch.npunch], "ABUN" );
381  sprintf( chHeader,
382  "#abund H" );
383  for(nelem=ipHELIUM;nelem<LIMELM; ++nelem )
384  {
385  sprintf( chTemp, "\t%s",
387  strcat( chHeader, chTemp );
388  }
389  strcat( chHeader, "\n");
390  }
391 
392  else if( nMatch(" AGE",chCard) )
393  {
394  /* punch ages */
395  strcpy( punch.chPunch[punch.npunch], "AGES" );
396  sprintf( chHeader,
397  "#ages depth\tt(cool)\tt(H2 dest)\tt(CO dest)\tt(OH dest)\tt(H rec)\n" );
398  }
399 
400  else if( nMatch(" AGN",chCard) )
401  {
402  /* punch tables needed for AGN3 */
403  strcpy( punch.chPunch[punch.npunch], " AGN" );
404  /* this is the AGN option, to produce a table for AGN */
405 
406  /* charge exchange rate coefficients */
407  if( nMatch("CHAR",chCard) )
408  {
409  strcpy( punch.chPunch[punch.npunch], "CHAG" );
410  sprintf( chHeader,
411  "#charge exchange rate coefficnt\n" );
412  }
413 
414  else if( nMatch("RECO",chCard) )
415  {
416  /* punch recombination rates for AGN3 table */
417  strcpy( punch.chPunch[punch.npunch], "RECA" );
418  sprintf( chHeader,
419  "#Recom rates for AGN3 table\n" );
420  }
421 
422  else if( nMatch("OPAC",chCard) )
423  {
424  /* create table for appendix in AGN */
425  strcpy( punch.chOpcTyp[punch.npunch], " AGN" );
426  strcpy( punch.chPunch[punch.npunch], "OPAC" );
427  }
428 
429  else if( nMatch("HECS",chCard) )
430  {
431  /* create table for appendix in AGN */
432  strcpy( punch.chPunchArgs[punch.npunch], "HECS" );
433  sprintf( chHeader,
434  "#AGN3 he cs \n" );
435  }
436 
437  else if( nMatch("HEMI",chCard) )
438  {
439  /* HEMIS - continuum emission needed for chap 4 of AGN3 */
440  strcpy( punch.chPunchArgs[punch.npunch], "HEMI" );
441 
442  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
443  * units are copied into punch.chConPunEnr */
444  ChkUnits(chCard);
445  }
446  else if( nMatch("RECC",chCard) )
447  {
448  /* recombination cooling, for AGN */
449  strcpy( punch.chPunch[punch.npunch], "HYDr" );
450  sprintf( chHeader,
451  "#T\tbAS\tb1\tbB\n" );
452  }
453  else
454  {
455  fprintf( ioQQQ, " I did not recognize this option on the PUNCH HYDROGEN command.\n" );
456  fprintf( ioQQQ, " Sorry.\n" );
457  cdEXIT(EXIT_FAILURE);
458  }
459  }
460 
461  else if( nMatch("ASSE",chCard) )
462  {
463  /* punch asserts */
464  strcpy( punch.chPunch[punch.npunch], "ASSE" );
465  /* no need to print this standard line of explanation*/
466  /*sprintf( chHeader, " asserts\n" );*/
467  }
468 
469  else if( nMatch("AVER",chCard) )
470  {
471  /* punch averages */
472  strcpy( punch.chPunch[punch.npunch], "AVER" );
473  /* no need to print this standard line of explanation*/
474  /*sprintf( chHeader, " asserts\n" );*/
475 
476  /* actually get the averages from the input stream, and malloc the
477  * space in the arrays
478  * punch io unit not used in read */
479  punch_average( punch.npunch, "READ", chHeader );
480  }
481 
482  /* punch charge transfer */
483  else if( nMatch("CHAR",chCard) && nMatch("TRAN",chCard) )
484  {
485  /* NB in PunchDo only the first three characters are compared to find this option,
486  * search for "CHA" */
487  /* punch charge transfer */
488  strcpy( punch.chPunch[punch.npunch], "CHAR" );
489  sprintf( chHeader,
490  "#charge exchange rate coefficient\n" );
491  }
492 
493  else if( nMatch("CHEM",chCard) )
494  {
495 
496  if( nMatch( "RATE" , chCard ) )
497  {
498  /* >>chng 06 May 30, NPA. Punch reaction rates for selected species */
499 
500  if( nMatch( " CO " , chCard ) )
501  {
502 
503  strcpy( punch.chPunch[punch.npunch], "RCO " );
504  sprintf( chHeader,
505  "#Depth\tH2O_C3HP_CO_C2H3P\tC2H2_HCOP_CO_C2H3P\tC2H_HCOP_CO_C2H2P\tO_C3H_CO_C2H\t"
506  "O_C2H2_CO_CH2\tOH_C2H2_CO_CH3\tHCOP_C3_C3HP_CO\tO2_C3_CO_C2_O\tO_C3_CO_C2\t"
507  "H_COP_CO_HP\tHminus_HCOP_CO_H2\tH2P_CO_HCOP_H\tH2P_CO_COP_H2\tH3P_CO_HCOP_H2\t"
508  "HeP_CO_OP_C_He\tHeP_CO_O_CP_He\tcrnu_CO_O_C\tCRP_CO_COP_e\tnu_CO_O_C\te_HCOP_CO_H\t"
509  "C_COP_CO_CP\tC_HCOP_CO_CHP\tC_O_CO_nu\tC_O2_CO_O\tC_OH_CO_H\tC_SiOP_SiP_CO\tCP_OH_CO_HP\t"
510  "CP_SiO_SiP_CO\tCP_O2_CO_OP\tO_CH_CO_H\tO_CH2_CO_H_H\tO_CH2_CO_H2\tO_COP_CO_OP\tOP_CO_COP_O\t"
511  "CH_COP_CO_CHP\tCH_HCOP_CO_CH2P\tCH_O2_CO_OH\tCH2_COP_CO_CH2P\tCH2_HCOP_CO_CH3P\t"
512  "CH2_O2_CO_H2O\tCOP_O2_O2P_CO\tH2O_COP_CO_H2OP\tH2O_HCOP_CO_H3OP\tH2OP_CO_HCOP_OH\t"
513  "HCOP_SiH_SiH2P_CO\tHCOP_SiO_SiOHP_CO\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tOHP_CO_HCOP_O\t"
514  "COP_CH4_CO_CH4P\tCO_CH4P_HCOP_CH3\tCO_CH5P_HCOP_CH4\tHP_OCS_HSP_CO\tHeP_OCS_SP_CO_He\t"
515  "OCNP_e_CO_N\tOCSP_e_S_CO\tOCS_NU_S_CO\tOCS_CRP_S_CO\tC_NO_CO_N\tC_OCN_CO_CN\t"
516  "C_SO_S_CO\tO_CN_CO_N\tO_HCN_CO_NH\tO_OCN_NO_CO\tO_CS_S_CO\tO_OCS_SO_CO\tOH_HCN_CO_NH2\t"
517  "CN_NO_N2_CO\tCN_O2_NO_CO\tCO_HS_OCS_H\tCP_SO_SP_CO\tCP_OCS_CSP_CO\tCHP_OCS_HCSP_CO\t"
518  "NP_CO_NOP_C\tNP_OCS_SP_CO_N\tNHP_CO_HCOP_N\tNHP_CO_OCNP_H\tNH_HCOP_CO_NH2P\t"
519  "NH2_HCOP_CO_NH3P\tNH3_HCOP_CO_NH4P\tCNP_O2_NOP_CO\tHCNP_CO_HCOP_CN\tCO_HNOP_NO_HCOP\t"
520  "N2P_OCS_SP_N2_CO\tHCOP_S_HSP_CO\tHCOP_CS_HCSP_CO\tNH_COP_CO_NHP\tNH2_COP_CO_NH2P\t"
521  "NH3_COP_CO_NH3P\tCNP_CO_COP_CN\tHCN_COP_CO_HCNP\tCO_N2P_N2_COP\tCOP_NO_NOP_CO\t"
522  "CO_S_OCS_NU\tNP_CO_COP_N\tCOP_S_SP_CO\tC_CO_C2_O\tO_C2_CO_C\tC2_O2_CO_CO\t"
523  "C2_O2P_COP_CO\tC2_COP_CO_C2P\tC2P_O2_COP_CO\tO_C2H_CO_CH\tO_CCl_Cl_CO\t"
524  "CO_H2ClP_HCl_HCOP\tHNC_HCOP_HCNHP_CO\tHCN_HCOP_HCNHP_CO\tC2H_COP_CO_C2HP\tC2_HCOP_CO_C2HP\n");
525  }
526 
527  else if( nMatch( " OH " , chCard ) )
528  {
529  strcpy( punch.chPunch[punch.npunch], "ROH " );
530  sprintf( chHeader,
531  "#Depth\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tcrnu_H2O_OH_H\tcrnu_OH_O_H\tCP_OH_CO_HP\tCP_OH_COP_H\t"
532  "CP_OH_CO_HP\tCP_OH_COP_H\tCH2_OHP_OH_CH2P\tCH2P_O2_HCOP_OH\tH2O_COP_HCOP_OH\tH2OP_H2O_H3OP_OH\tC_H2OP_OH_CHP\tOP_OH_O2P_H\tOP_OH_OHP_O\tSiP_OH_SiOP_H\t"
533  "CH_H2OP_OH_CH2P\tCH_OHP_OH_CHP\tCHP_OH_COP_H2\tCHP_O2_COP_OH\tCH2_H2OP_OH_CH3P\tOH_COP_HCOP_O\tOH_H2OP_H3OP_O\tOHP_H2O_H2OP_OH\tOHP_O2_O2P_OH\tOHP_OH_H2OP_O\t"
534  "O_CH4P_OH_CH3P\tOP_CH4_OH_CH3P\tCH5P_OH_H2OP_CH4\tOHP_C2_C2P_OH\tOH_C2P_C2_OHP\tCH_NO_CN_OH\tNH_O_OH_N\tNH_OH_HNO_H\tO_HNO_NO_OH\tNH2_OH_H2O_NH\t"
535  "NH2_NO_N2_OH_H\tOH_CN_OCN_H\tOH_S_HS_O\tOH_S_SO_H\tNHP_OH_H2OP_N\tNHP_H2O_OH_NH2P\tNHP_O2_NOP_OH\tO_HSP_SP_OH\tNH2P_OH_H2OP_NH\tNH2P_H2O_NH3P_OH\t"
536  "NH2_H2OP_NH3P_OH\tNH2P_O2_HNOP_OH\tOH_NH3P_NH4P_O\tOH_HCNP_CN_H2OP\tOH_HNOP_NO_H2OP\tOH_SP_SOP_H\tNH3P_H2O_NH4P_OH\tNH3_H2OP_NH4P_OH\tH2O_CNP_HCNP_OH\tH2OP_S_HSP_OH\t"
537  "NH_OHP_OH_NHP\tNH2_OHP_OH_NH2P\tOHP_NH3_NH3P_OH\tOH_CNP_CN_OHP\tOH_N2P_N2_OHP\tOHP_NO_NOP_OH\tNP_OH_OHP_N\tOHP_S_SP_OH\tH2OP_HNC_HCNHP_OH\tH2OP_HCN_HCNHP_OH\t"
538  "H2O_C2P_C2HP_OH\tH2OP_C2_C2HP_OH\tH2OP_C2H_C2H2P_OH\tOHP_C2H_C2HP_OH\tHminus_OH_H2O_e\tHminus_O_OH_e\tHP_OH_OHP_H\tH2s_OH_O_H2_H\tH2s_H2O_OH_H2_H\tH2P_OH_H2OP_H\t"
539  "H2P_OH_OHP_H2\tH3P_OH_H2OP_H2\tHeP_OH_OP_He_H\tHeP_H2O_OH_He_HP\tH3P_NO2_NOP_OH_H2\tCH_O2_CO_OH\tH2OP_CO_HCOP_OH\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tCH2_OH_H2O_CH\t"
540  "C_OH_O_CH\tO_CH_OH_C\tO_CH2_OH_CH\tO_H2O_OH_OH\tO_OH_O2_H\tSi_OH_SiO_H\tOH_OH_H2O_O\tO_CH4_OH_CH3\tOH_CH2_O_CH3\tOH_CH3_CH4_O\t"
541  "OH_CH3_H2O_CH2\tH2O_CH3_OH_CH4\tOH_CH4_H2O_CH3\tN_OH_O_NH\tN_OH_NO_H\tCH2_NO_HCN_OH\tNH_OH_NH2_O\tNH_OH_H2O_N\tNH_H2O_OH_NH2\tNH_NO_N2_OH\t"
542  "NH_NO2_N2O_OH\tO_NH2_OH_NH\tO_NH3_OH_NH2\tO_HCN_CN_OH\tO_HS_S_OH\tNH2_OH_NH3_O\tOH_NH3_H2O_NH2\tOH_CN_HCN_O\tOH_HCN_CN_H2O\tOH_NO_NO2_H\t"
543  "OH_N2O_HNO_NO\tOH_CS_OCS_H\tNO_HNO_N2O_OH\tO_C2H2_C2H_OH\tOH_C2H2_C2H_H2O\te_H2OP_OH_H\te_H3OP_OH_H_H\te_H3OP_OH_H2\te_SiOHP_Si_OH\tH_OH_O_H_H\t"
544  "H_H2O_OH_H_H\tH_OH_O_H2\tH_H2O_OH_H2\tH_OH_H2O_nu\tHminus_H3OP_OH_H2_H\tH_O_OH_nu\tH2_OH_H2O_H\tH2_O_OH_H\tH2_OH_O_H2_H\tH2_H2O_OH_H2_H\t"
545  "H2_O2_OH_OH\tH_O2_OH_O\tC_OH_CO_H\tOH_HCN_CO_NH2\tOH_C2H2_CO_CH3\tOH_C2H2_CO_CH3\n");
546  }
547 
548  else
549  {
550  fprintf(ioQQQ," The keyword CO or OH must appear on the PUNCH CHEMISTRY RATES command.\n");
551  fprintf( ioQQQ, " Sorry.\n" );
552  cdEXIT(EXIT_FAILURE);
553  }
554 
555  }
556  }
557 
558  else if( nMatch("COMP",chCard) )
559  {
560  /* punch Compton, details of the energy exchange problem */
561  strcpy( punch.chPunch[punch.npunch], "COMP" );
562  sprintf( chHeader,
563  "#nu, comup, comdn\n" );
564  }
565 
566  else if( nMatch("COOL",chCard) )
567  {
568  /* punch cooling, actually done by routine cool_punch */
569  strcpy( punch.chPunch[punch.npunch], "COOL" );
570  /*>>chng 06 jun 06, revise to be same as punch cooling */
571  sprintf( chHeader,
572  "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\tcool fracs\n" );
573  }
574 
575  else if( nMatch("DYNA",chCard) )
576  {
577  /* punch something dealing with dynamics
578  * in PunchDo the DYN part of key is used to call DunaPunch,
579  * with the 4th char as its second argument. DynaPunch uses that
580  * 4th letter to decide the job */
581  if( nMatch( "ADVE" , chCard ) )
582  {
583  /* punch information relating to advection */
584  strcpy( punch.chPunch[punch.npunch], "DYNa");
585  sprintf( chHeader,
586  "#advection depth\tHtot\tadCool\tadHeat\tdCoolHeatdT\t"
587  "Source[hyd][hyd]\tRate\tEnthalph\tadSpecEnthal\n" );
588  }
589  else
590  {
591  fprintf( ioQQQ, " I did not recognize a sub keyword on this PUNCH DYNAMICS command.\n" );
592  fprintf( ioQQQ, " Sorry.\n" );
593  cdEXIT(EXIT_FAILURE);
594  }
595  }
596 
597  else if( nMatch("ENTH",chCard) )
598  {
599  /* contributors to the total enthalpy */
600  strcpy( punch.chPunch[punch.npunch], "ENTH" );
601  sprintf( chHeader,
602  "#depth\tTotal\tExcit\tIoniz\tBind\tKE\tther+PdV\tmag \n" );
603  }
604 
605  else if( nMatch("EXEC",chCard) && nMatch("TIME",chCard) )
606  {
607  /* output the execution time per zone */
608  strcpy( punch.chPunch[punch.npunch], "XTIM" );
609  sprintf( chHeader,
610  "#zone\tdTime\tElapsed t\n" );
611  }
612 
613  else if( nMatch("FEII",chCard) || nMatch("FE II",chCard) )
614  {
615  /* something to do with FeII atom - several options
616  * FeII column densities */
617  if( nMatch("COLU",chCard) && nMatch("DENS",chCard) )
618  {
619  /* punch FeII column density */
620  strcpy( punch.chPunch[punch.npunch], "FEcl" );
621 
622  /* file will give energy, statistical weight, and column density [cm-2] */
623  sprintf( chHeader,
624  "#FeII: energy\tstat wght\tcol den\n" );
625  }
626 
627  /* FeII continuum, only valid if large atom is set */
628  else if( nMatch("CONT",chCard) )
629  {
630  /* punch FeII continuum */
631  strcpy( punch.chPunch[punch.npunch], "FEco" );
632 
633  sprintf( chHeader,
634  "#FeII: Wl(A)\tInt[erg cm-2 s-1]\n" );
635  }
636 
637  else if( nMatch("DEPA",chCard) )
638  {
639  /* punch out departure coefficient for the large FeII atom */
640  sprintf( chHeader,
641  "#FeII departure coefficient \n" );
642  /* optional keyword all means do all levels, if not then just do subset */
643  if( nMatch(" ALL",chCard) )
644  {
645  /* punch all levels, calls routine FeIIPunDepart */
646  strcpy( punch.chPunch[punch.npunch], "FE2D" );
647  }
648  else
649  {
650  /* punch a very few selected levels, calls routine FeIIPunDepart */
651  strcpy( punch.chPunch[punch.npunch], "FE2d" );
652  }
653  }
654 
655  else if( nMatch("LEVE",chCard) )
656  {
657  /* punch level energies and statistical weights for large FeII atom */
658  sprintf( chHeader,
659  "#FeII energy(wn)\tstat weight\n" );
660  strcpy( punch.chPunch[punch.npunch], "FE2l" );
661  }
662 
663  else if( nMatch("LINE",chCard) )
664  {
665  /* punch FeII lines command
666  * three optional parameters, threshold for faintest
667  * line to print, lower and upper energy bounds */
668 
669  /* punch intensities from large FeII atom */
670  strcpy( punch.chPunch[punch.npunch], "FEli" );
671 
672  /* short keyword makes punch half as big */
673  if( nMatch("SHOR",chCard) )
674  {
675  FeII.lgShortFe2 = true;
676  }
677  else
678  {
679  FeII.lgShortFe2 = false;
680  }
681 
682  /* first optional number changes the threshold of weakest line to print*/
683  i = 5;
684  /* fe2thresh is intensity relative to normalization line,
685  * normally Hbeta, and is set to zero in zero.c */
686  FeII.fe2thresh = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
687  if( lgEOL )
688  {
689  FeII.fe2thresh = 0.;
690  }
691 
692  /* it is a log if negative */
693  if( FeII.fe2thresh < 0. )
694  {
695  FeII.fe2thresh = (realnum)pow((realnum)10.f,FeII.fe2thresh);
696  }
697 
698  /* check for energy range (Rydberg) of lines to be punched,
699  * this is to limit size of output file */
700  FeII.fe2ener[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
701  if( lgEOL )
702  {
703  FeII.fe2ener[0] = 0.;
704  }
705 
706  FeII.fe2ener[1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
707  if( lgEOL )
708  {
709  FeII.fe2ener[1] = 1e8;
710  }
711  /* if either is negative then both are logs */
712  if( FeII.fe2ener[0] < 0. || FeII.fe2ener[1] < 0. )
713  {
714  FeII.fe2ener[0] = (realnum)pow((realnum)10.f,FeII.fe2ener[0]);
715  FeII.fe2ener[1] = (realnum)pow((realnum)10.f,FeII.fe2ener[1]);
716  }
717 
718  /* entered in Ryd in above, convert to wavenumbers */
719  FeII.fe2ener[0] /= (realnum)WAVNRYD;
720  FeII.fe2ener[1] /= (realnum)WAVNRYD;
721 
722  /* these results are actually created by the FeIIPunchLines routine
723  * that lives in the FeIILevelPops file */
724  sprintf( chHeader,
725  "#FeII ipHi\tipLo\tWL(A)\tlog I\tI/Inorm\t\tTau\n" );
726  }
727 
728  else if( nMatch("OPTI",chCard) && nMatch("DEPT",chCard) )
729  {
730  /* punch optical depths for large FeII atom */
731  sprintf( chHeader,
732  "#FeII hi\tlow\twl(A)\ttau\n" );
733  strcpy( punch.chPunch[punch.npunch], "FE2o" );
734  }
735 
736  else if( nMatch("POPU",chCard) )
737  {
738  /* punch out level populations for the large FeII atom */
739  sprintf( chHeader,
740  "#FeII level populations [cm^-3]\n" );
741 
742  /* this is keyword RELATIVE that says to punch relative to total Fe+,
743  * default is actual density (cm-3) */
744  if( nMatch("RELA",chCard) )
745  {
746  punch.punarg[punch.npunch][2] = 0.;
747  }
748  else
749  {
750  /* default is to punch density (cm-3) */
751  punch.punarg[punch.npunch][2] = 1.;
752  }
753 
754  /* optional keyword all means do all levels, if not then just do subset */
755  if( nMatch(" ALL",chCard) )
756  {
757  /* punch all levels, calls routine FeIIPunPop */
758  strcpy( punch.chPunch[punch.npunch], "FE2P" );
759  punch.punarg[punch.npunch][0] = 0.;
761  }
762 
763  /* optional keyword range means read lower and upper bound, do these */
764  else if( nMatch("RANG",chCard) )
765  {
766  /* punch range of levels, calls routine FeIIPunPop */
767  strcpy( punch.chPunch[punch.npunch], "FE2P" );
768  ipFFmt = 5;
769  punch.punarg[punch.npunch][0] =
770  (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL);
771  punch.punarg[punch.npunch][1] =
772  (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL);
773  if( lgEOL || punch.punarg[punch.npunch][0] <0 ||
776  {
777  fprintf( ioQQQ, "There must be two numbers on this punch FeII populations range command.\n" );
778  fprintf( ioQQQ, "These give the lower and upper levels for the range of FeII levels.\n" );
779  fprintf( ioQQQ, "The first must be less than the second and the second must be <= the number of FeII levels.\n" );
780  fprintf( ioQQQ, "Sorry.\n" );
781  cdEXIT(EXIT_FAILURE);
782  }
783  }
784 
785  else
786  {
787  /* punch a very few selected levels, calls routine FeIIPunPop */
788  strcpy( punch.chPunch[punch.npunch], "FE2p" );
789  }
790  }
791  else
792  {
793  fprintf( ioQQQ, "There must be a second keyword on this PUNCH FEII command.\n" );
794  fprintf( ioQQQ, "The ones I know about are COLUmn, CONTinuum, "
795  "DEPArture, LEVEls, LINE, OPTIcal DEPTh, and POPUlations.\n" );
796  fprintf( ioQQQ, "Sorry.\n" );
797  cdEXIT(EXIT_FAILURE);
798  }
799  }
800 
801  /* the punch continuum command, with many options,
802  * the first 3 char of the chPunch flag will always be "CON"
803  * with the last indicating which one */
804  else if( nMatch("CONT",chCard) && !nMatch("XSPE",chCard) )
805  {
806  /* this flag is checked in PrtComment to generate a caution
807  * if continuum is punched but iterations not performed */
808  punch.lgPunContinuum = true;
809 
810  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
811  * units are copied into punch.chConPunEnr */
812  ChkUnits(chCard);
813 
814  if( nMatch("BINS",chCard) )
815  {
816  /* continuum binning */
817  strcpy( punch.chPunch[punch.npunch], "CONB" );
818 
819  sprintf( chHeader,
820  "#Continuum binning, enrOrg, Energy, width of cells\n" );
821  }
822 
823  else if( nMatch("DIFF",chCard) )
824  {
825  /* diffuse continuum, the locally emitted continuum stored in rfield.ConEmitLocal */
826  strcpy( punch.chPunch[punch.npunch], "COND" );
827 
828  sprintf( chHeader,
829  "#energy\t ConEmitLocal \n" );
830  }
831 
832  else if( nMatch("EMIT",chCard) )
833  {
834  /* continuum emitted by cloud */
835  strcpy( punch.chPunch[punch.npunch], "CONE" );
836 
837  sprintf( chHeader,
838  "#Energy\treflec\toutward\ttotal\tline\tcont\n" );
839  }
840 
841  else if( nMatch("FINE" , chCard ) )
842  {
843  rfield.lgPunchOpacityFine = true;
844  /* fine transmitted continuum cloud */
845  strcpy( punch.chPunch[punch.npunch], "CONf" );
846 
847  sprintf( chHeader,
848  "#Energy\tTransmitted\n" );
849 
850  ipFFmt = 5;
851  /* range option - important since so much data */
852  if( nMatch("RANGE",chCard) )
853  {
854  /* get lower and upper range, must be in Ryd */
855  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
856  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
857  if( lgEOL )
858  {
859  fprintf(ioQQQ,"There must be two numbers, the lower and upper energies in Ryd.\nSorry.\n");
860  cdEXIT(EXIT_FAILURE);
861  }
862  if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
863  {
864  fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
865  cdEXIT(EXIT_FAILURE);
866  }
867  }
868  else
869  {
870  /* these mean full energy range */
871  punch.punarg[punch.npunch][0] = 0.;
872  punch.punarg[punch.npunch][1] = 0.;
873  }
874  /* optional last parameter - how many points to bring together */
875  punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
876  INPUT_LINE_LENGTH,&lgEOL);
877 
878  /* default is to bring together ten */
879  if( lgEOL )
880  punch.punarg[punch.npunch][2] = 10;
881 
882  if( punch.punarg[punch.npunch][2] < 1 )
883  {
884  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
885  cdEXIT(EXIT_FAILURE);
886  }
887  }
888 
889  else if( nMatch("GRAI",chCard) )
890  {
891  /* punch grain continuum in optically thin limit */
892  strcpy( punch.chPunch[punch.npunch], "CONG" );
893 
894  sprintf( chHeader,
895  "#energy\tgraphite\trest\ttotal\n" );
896  }
897 
898  else if( nMatch("INCI",chCard) )
899  {
900  /* incident continuum */
901  strcpy( punch.chPunch[punch.npunch], "CONC" );
902 
903  sprintf( chHeader,
904  "#Incident Continuum, Enr\tnFn \n" );
905  }
906 
907  else if( nMatch("INTE",chCard) )
908  {
909  /* continuum interactions */
910  strcpy( punch.chPunch[punch.npunch], "CONi" );
911 
912  sprintf( chHeader,
913  "#Continuum interactions, inc, otslin. otscon, ConInterOut, outlin \n" );
914  /* this is option for lowest energy, if nothing then zero */
915  ipFFmt = 3;
916  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,
917  INPUT_LINE_LENGTH,&lgEOL);
918  }
919 
920  else if( nMatch("IONI",chCard) )
921  {
922  /* punch ionizing continuum*/
923  strcpy( punch.chPunch[punch.npunch], "CONI" );
924 
925  /* this is option for lowest energy, if nothing then zero */
926  ipFFmt = 3;
927  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,
928  INPUT_LINE_LENGTH,&lgEOL);
929 
930  /* this is option for smallest interaction to punch, def is 1 percent */
931  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
932  if( lgEOL )
933  {
934  punch.punarg[punch.npunch][1] = 0.01f;
935  }
936 
937  /* >>chng 03 may 03, add "every" option to punch this on every zone -
938  * if every is not present then only last zone is punched */
939  if( nMatch("EVER", chCard ) )
940  {
941  /* punch every zone */
944  }
945  else
946  {
947  /* only punch last zone */
950  }
951 
952  /* put the header at the top of the file */
953  sprintf( chHeader,
954  "#cell(on C scale)\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" );
955  }
956 
957  else if( nMatch("OUTW",chCard) )
958  {
959  /* outward only continuum */
960  if( nMatch("LOCA",chCard) )
961  {
962  strcpy( punch.chPunch[punch.npunch], "CONo" );
963  sprintf( chHeader,
964  "#Local Out ConInterOut+line SvOt*opc pass*opc\n" );
965  }
966  else
967  {
968  strcpy( punch.chPunch[punch.npunch], "CONO" );
969  sprintf( chHeader,
970  "#Out Con OutIncid OutConD OutLinD OutConS\n" );
971  }
972  }
973 
974  else if( nMatch("TRAN",chCard) )
975  {
976  /* transmitted continuum */
977  strcpy( punch.chPunch[punch.npunch], "CONT" );
978 
979  sprintf( chHeader,
980  "#ener\tTran Contin\ttrn coef\n" );
981  }
982 
983  else if( nMatch(" TWO",chCard) )
984  {
985  /* total two photon continua rfield.TotDiff2Pht */
986  strcpy( punch.chPunch[punch.npunch], "CON2" );
987 
988  sprintf( chHeader,
989  "#energy\t n_nu\tnuF_nu \n" );
990  }
991 
992  else if( nMatch(" RAW",chCard) )
993  {
994  /* "raw" continua */
995  strcpy( punch.chPunch[punch.npunch], "CORA" );
996 
997  sprintf( chHeader,
998  "#Raw Con anu\tflux\totslin\totscon\tConRefIncid\tConEmitReflec\tConInterOut\toutlin\tConEmitOut\tline\tcont\tnLines\n" );
999  }
1000 
1001  else if( nMatch("REFL",chCard) )
1002  {
1003  /* reflected continuum */
1004  strcpy( punch.chPunch[punch.npunch], "CONR" );
1005 
1006  sprintf( chHeader,
1007  "#Reflected\tcont\tline\ttotal\talbedo\tConID\n" );
1008  }
1009 
1010  else
1011  {
1012  /* this is the usual punch continuum command */
1013  strcpy( punch.chPunch[punch.npunch], "CON " );
1014  sprintf( chHeader,
1015  "#Cont nu\tincident\ttrans\tDiffOut\tnet trans\treflc\ttotal\tline\tcont\tnLine\n" );
1016 
1017  /* >>chng 06 apr 03, add "every" option to punch this on every zone -
1018  * if every is not present then only last zone is punched */
1019  if( nMatch("EVER", chCard ) )
1020  {
1021  /* punch every zone */
1023  /* option to say how many to skip */
1024  ipFFmt = 5;
1025  punch.nPunchEveryZone[punch.npunch] = (long)FFmtRead(chCard,&ipFFmt,
1026  INPUT_LINE_LENGTH,&lgEOL);
1027  if( lgEOL )
1029  }
1030  else
1031  {
1032  /* only punch last zone */
1035  }
1036  }
1037  }
1038 
1039  /* punch information about convergence of this model
1040  * reason - why it did not converge an iteration
1041  * error - zone by zone display of various convergence errors */
1042  else if( nMatch("CONV",chCard) )
1043  {
1044  if( nMatch("REAS",chCard) )
1045  {
1046  /* this does not count as a punch option (really) */
1047  punch.lgPunConv = true;
1048  /* this is done below */
1049  strcpy( punch.chPunch[punch.npunch], "" );
1050  punch.lgRealPunch[punch.npunch] = false;
1051  }
1052  else if( nMatch("ERRO",chCard) )
1053  {
1054  /* punch zone by zone errors in pressure, electron density, and heating-cooling */
1055  /* convergence error */
1056  strcpy( punch.chPunch[punch.npunch], "CNVE" );
1057  sprintf( chHeader,
1058  "#depth\tnPres2Ioniz\tP(cor)\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" );
1059  }
1060  else if( nMatch("BASE",chCard) )
1061  {
1062  /* punch converged quantities in Converge base for each pass through
1063  * solvers - not one pass per zone */
1064  strcpy( punch.chPunch[punch.npunch], "CNVB" );
1065  strcpy( punch.chPunch[punch.npunch], "" );
1066  punch.lgRealPunch[punch.npunch] = false;
1067  }
1068  else
1069  {
1070  fprintf( ioQQQ, "There must be a second keyword on this command.\n" );
1071  fprintf( ioQQQ, "The ones I know about are REASON and ERROR.\n" );
1072  fprintf( ioQQQ, "Sorry.\n" );
1073  cdEXIT(EXIT_FAILURE);
1074  }
1075  }
1076 
1077  else if( nMatch(" DR ",chCard) )
1078  {
1079  /* first occurrence of punch dr to follow choice in change of zoning */
1080  punch.lgDROn = true;
1081  strcpy( punch.chPunch[punch.npunch], "" );
1082  punch.lgRealPunch[punch.npunch] = false;
1083  }
1084 
1085  else if( nMatch("ELEM",chCard) && !nMatch("GAMMA" , chCard) )
1086  {
1087  /* option to punch ionization structure of some element
1088  * will give each stage of ionization, vs depth */
1089  strcpy( punch.chPunch[punch.npunch], "ELEM" );
1090 
1091  /* this returns element number on c scale */
1092  /* >>chng 04 nov 23, had converted to f scale, leave on c */
1093  nelem = GetElem(chCard);
1094 
1095  /* this is the atomic number on the c scale */
1096  punch.punarg[punch.npunch][0] = (realnum)nelem;
1097  ASSERT( nelem>=ipHYDROGEN);
1098 
1099  /* >>chng 04 nov 24, add DENSE option to print density rather than fraction */
1100  punch.punarg[punch.npunch][1] = 0;
1101  if( nMatch("DENS",chCard) )
1102  punch.punarg[punch.npunch][1] = 1.;
1103 
1104  /* start printing header line - first will be the depth in cm */
1105  sprintf( chHeader, "#depth");
1106 
1107  /* next come the nelem+1 ion stages */
1108  for(i=0; i<=nelem+1;++i )
1109  {
1110  sprintf( chTemp,
1111  "\t%2s%2li", elementnames.chElementSym[nelem],i+1);
1112  strcat( chHeader, chTemp );
1113  }
1114 
1115  /* finally some fine structure or molecular populations */
1116  /* >>chng 04 nov 23, add fs pops of C, O
1117  * >>chng 04 nov 25, add molecules */
1118  if( nelem==ipHYDROGEN )
1119  {
1120  sprintf( chTemp, "\tH2");
1121  strcat( chHeader, chTemp );
1122  }
1123  if( nelem==ipCARBON )
1124  {
1125  sprintf( chTemp, "\tC1\tC1*\tC1**\tC2\tC2*\tCO");
1126  strcat( chHeader, chTemp );
1127  }
1128  else if( nelem==ipOXYGEN )
1129  {
1130  sprintf( chTemp, "\tO1\tO1*\tO1**");
1131  strcat( chHeader, chTemp );
1132  }
1133 
1134  /* finally the new line */
1135  strcat( chHeader, "\n");
1136  }
1137 
1138  else if( nMatch("FITS",chCard) )
1139  {
1140 
1141 #ifdef FLT_IS_DBL
1142  fprintf( ioQQQ, "Punching FITS files is not currently supported in double precision.\n" );
1143  fprintf( ioQQQ, "Please recompile without the FLT_IS_DBL option.\n" );
1144  fprintf( ioQQQ, "Sorry.\n" );
1145  cdEXIT(EXIT_FAILURE);
1146 #else
1147  /* say that this is a FITS file output */
1148  punch.lgFITS[punch.npunch] = true;
1149 
1150  strcpy( punch.chPunch[punch.npunch], "FITS" );
1151 #endif
1152 
1153  }
1154 
1155  else if( nMatch("FRED",chCard) )
1156  {
1157  /* punch out some stuff for Fred's dynamics project */
1158  sprintf( chHeader,
1159  "#Radius\tDepth\tVelocity\thden\teden\tTemperature\tRadAccel line\tRadAccel con\t"
1160  "Force multiplier\t"
1161  "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
1162  "O2\tO3\tO4\tO5\tO6\tO7\tO8\t"
1163  "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t"
1164  "O2\tO3\tO4\tO5\tO6\tO7\tO8\tMg2\tMg2\n");
1165  strcpy( punch.chPunch[punch.npunch], "FRED" );
1166  }
1167 
1168  else if( nMatch("GAMM",chCard) )
1169  {
1170  /* punch all photoionization rates for all subshells */
1171  sprintf( chHeader,
1172  "#Photoionization rates \n" );
1173  if( nMatch("ELEMENT" , chCard ) )
1174  {
1175  /* element keyword, find element name and stage of ionization,
1176  * will print photoionization rates for valence of that element */
1177  strcpy( punch.chPunch[punch.npunch], "GAMe" );
1178 
1179  /* this returns element number on c scale */
1180  nelem = GetElem(chCard);
1181  /* this is the atomic number on the C scale */
1182  punch.punarg[punch.npunch][0] = (realnum)nelem;
1183 
1184  /* this will become the ionization stage on C scale */
1185  ipFFmt = 5;
1186  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,
1187  INPUT_LINE_LENGTH,&lgEOL) - 1;
1188  if( lgEOL )
1189  NoNumb( chCard );
1190  if( punch.punarg[punch.npunch][1]<0 || punch.punarg[punch.npunch][1]> nelem+1 )
1191  {
1192  fprintf(ioQQQ,"Bad ionization stage - please check Hazy.\nSorry.\n");
1193  cdEXIT(EXIT_FAILURE);
1194  }
1195  }
1196  else
1197  {
1198  /* no element - so make table of all rates */
1199  strcpy( punch.chPunch[punch.npunch], "GAMt" );
1200  }
1201 
1202  }
1203  else if( nMatch("GRAI",chCard) )
1204  {
1205  /* punch grain ... options */
1206  if( nMatch("OPAC",chCard) )
1207  {
1208  /* check for keyword UNITS on line, then scan wavelength or energy units,
1209  * sets punch.chConPunEnr*/
1210  ChkUnits(chCard);
1211 
1212  strcpy( punch.chPunch[punch.npunch], "DUSO" );
1213  /* punch grain opacity command in twice, here and above in opacity */
1214  sprintf( chHeader,
1215  "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" );
1216  }
1217  else if( nMatch("ABUN",chCard) )
1218  {
1219  /* punch grain abundance */
1220  strcpy( punch.chPunch[punch.npunch], "DUSA" );
1221  sprintf( chHeader,
1222  "#grain\tdepth\tabundance (g/cm^3)\n" );
1223  }
1224  else if( nMatch("D/G ",chCard) )
1225  {
1226  /* punch grain dust/gas mass ratio */
1227  strcpy( punch.chPunch[punch.npunch], "DUSD" );
1228  sprintf( chHeader,
1229  "#grain\tdepth\tdust/gas mass ratio\n" );
1230  }
1231  else if( nMatch("PHYS",chCard) )
1232  {
1233  /* punch grain physical conditions */
1234  strcpy( punch.chPunch[punch.npunch], "DUSP" );
1235  sprintf( chHeader,
1236  "#grain\tdepth\tpotential\n" );
1237  }
1238  else if( nMatch(" QS ",chCard) )
1239  {
1240  strcpy( punch.chPunch[punch.npunch], "DUSQ" );
1241  sprintf( chHeader,
1242  "#grain\tnu\tQ_abs\tQ_scat\n" );
1243  }
1244  else if( nMatch("TEMP",chCard) )
1245  {
1246  /* punch temperatures of each grain species */
1247  strcpy( punch.chPunch[punch.npunch], "DUST" );
1248  /* cannot punch grain labels since they are not known yet */
1249  sprintf( chHeader,
1250  "#grain temperature\n" );
1251  }
1252  else if( nMatch("DRIF",chCard) )
1253  {
1254  /* punch drift velocity of each grain species */
1255  strcpy( punch.chPunch[punch.npunch], "DUSV" );
1256  /* cannot punch grain labels since they are not known yet */
1257  sprintf( chHeader,
1258  "#grain drift velocity\n" );
1259  }
1260  else if( nMatch("EXTI",chCard) )
1261  {
1262  /* punch grain extinction */
1263  strcpy( punch.chPunch[punch.npunch], "DUSE" );
1264  /* cannot punch grain labels since they are not known yet */
1265  sprintf( chHeader,
1266  "#depth\tA_V(extended)\tA_V(point)\n" );
1267  }
1268  else if( nMatch("CHAR",chCard) )
1269  {
1270  /* punch charge per grain (# elec/grain) for each grain species */
1271  strcpy( punch.chPunch[punch.npunch], "DUSC" );
1272  /* cannot punch grain labels since they are not known yet */
1273  sprintf( chHeader,
1274  "#grain charge\n" );
1275  }
1276  else if( nMatch("HEAT",chCard) )
1277  {
1278  /* punch heating due to each grain species */
1279  strcpy( punch.chPunch[punch.npunch], "DUSH" );
1280  /* cannot punch grain labels since they are not known yet */
1281  sprintf( chHeader,
1282  "#grain heating\n" );
1283  }
1284  else if( nMatch("POTE",chCard) )
1285  {
1286  /* punch floating potential of each grain species */
1287  strcpy( punch.chPunch[punch.npunch], "DUSP" );
1288  /* cannot punch grain labels since they are not known yet */
1289  sprintf( chHeader,
1290  "#grain\tdepth\tpotential\n" );
1291  }
1292  else if( nMatch("H2RA",chCard) )
1293  {
1294  /* punch grain H2rate - H2 formation rate for each type of grains */
1295  strcpy( punch.chPunch[punch.npunch], "DUSR" );
1296  /* cannot punch grain labels since they are not known yet */
1297  sprintf( chHeader,
1298  "#grain H2 formation rates\n" );
1299  }
1300  else
1301  {
1302  fprintf( ioQQQ, "There must be a second key on this GRAIN command; The options I know about follow (required key in CAPS):\n");
1303  fprintf( ioQQQ, "OPACity, ABUNdance, D/G mass ratio, PHYSical conditions, QS , TEMPerature, DRIFt velocity, EXTInction, CHARge, HEATing, POTEntial, H2RAtes\nSorry.\n" );
1304  cdEXIT(EXIT_FAILURE);
1305  }
1306  }
1307 
1308  else if( nMatch("GAUN",chCard) )
1309  {
1310  strcpy( punch.chPunch[punch.npunch], "GAUN" );
1311  sprintf( chHeader,
1312  "#Gaunt factors.\n" );
1313  }
1314  else if( nMatch("GRID",chCard) )
1315  {
1316  strcpy( punch.chPunch[punch.npunch], "GRID" );
1317  /* automatically generate no hash option */
1318  punch.lgHashEndIter[punch.npunch] = false;
1319  }
1320  else if( nMatch( "HIST" , chCard ) )
1321  {
1322  /* punch pressure history of current zone */
1323  if( nMatch( "PRES" , chCard ) )
1324  {
1325  /* punch pressure history - density - pressure for this zone */
1326  strcpy( punch.chPunch[punch.npunch], "HISp" );
1327  sprintf( chHeader,
1328  "#iter zon\tdensity\tpres cur\tpres corr\n" );
1329  }
1330  /* punch temperature history of current zone */
1331  else if( nMatch( "TEMP" , chCard ) )
1332  {
1333  /* punch pressure history - density - pressure for this zone */
1334  strcpy( punch.chPunch[punch.npunch], "HISt" );
1335  sprintf( chHeader,
1336  "#iter zon\ttemperature\theating\tcooling\n" );
1337  }
1338  }
1339 
1340  else if( nMatch("HTWO",chCard) )
1341  {
1342  fprintf(ioQQQ," Sorry, this command has been replaced with the "
1343  "PUNCH H2 CREATION and PUNCH H2 DESTRUCTION commands.\n");
1344  cdEXIT(EXIT_FAILURE);
1345  }
1346 
1347  /* QHEAT has to come before HEAT... */
1348  else if( nMatch("QHEA",chCard) )
1349  {
1350  /* this is just a dummy clause, do the work below after parsing is over.
1351  * this is a no-nothing, picked up to stop optimizer */
1352  ((void)0);
1353  }
1354 
1355  else if( nMatch("HEAT",chCard) )
1356  {
1357  /* punch heating */
1358  strcpy( punch.chPunch[punch.npunch], "HEAT" );
1359  /*>>chng 06 jun 06, revise to be same as punch cooling */
1360  sprintf( chHeader,
1361  "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\theat fracs\n" );
1362  }
1363 
1364  else if( nMatch("HELI",chCard) &&!( nMatch("IONI",chCard)))
1365  {
1366  /* punch helium & helium-like iso sequence, but not punch helium ionization rate
1367  * punch helium line wavelengths */
1368  if( nMatch("LINE",chCard) && nMatch("WAVE",chCard) )
1369  {
1370  strcpy( punch.chPunch[punch.npunch], "HELW" );
1371  sprintf( chHeader,
1372  "#wavelengths of lines from He-like ions\n" );
1373  }
1374  else
1375  {
1376  fprintf( ioQQQ, "punch helium has options: LINE WAVElength.\nSorry.\n" );
1377  cdEXIT(EXIT_FAILURE);
1378  /* no key */
1379  }
1380  }
1381 
1382  else if( nMatch("HUMM",chCard) )
1383  {
1384  strcpy( punch.chPunch[punch.npunch], "HUMM" );
1385  sprintf( chHeader,
1386  "#input to DHs routine.\n" );
1387  }
1388 
1389  else if( nMatch("HYDR",chCard) )
1390  {
1391  /* punch hydrogen physical conditions */
1392  if( nMatch("COND",chCard) )
1393  {
1394  strcpy( punch.chPunch[punch.npunch], "HYDc" );
1395  sprintf( chHeader,
1396  "#depth\tTe\tHDEN\tEDEN\tHI/H\tHII/H\tH2/H\tH2+/H\tH3+/H\tH-/H\n" );
1397  /* punch hydrogen ionization */
1398  }
1399 
1400  /* punch information on 21 cm excitation processes - accept either keyword 21cm or 21 cm */
1401  else if( nMatch("21 CM",chCard) ||nMatch("21CM",chCard))
1402  {
1403  /* punch information about 21 cm line */
1404  strcpy( punch.chPunch[punch.npunch], "21CM" );
1405  sprintf( chHeader,
1406  "#depth\tT(spin)\tT(kin)\tT(Lya/21cm)\tnLo\tnHi\tOccLya\ttau(21cm)"
1407  "\ttau(Lya)\topac(21 cm)\tn/Ts\ttau(21)\tTex(Lya)\tN(H0)/Tspin"
1408  "\tSum_F0\tSum_F1\tSum_T21\n" );
1409  }
1410 
1411  else if( nMatch("IONI",chCard) )
1412  {
1413  /* punch hydrogen ionization */
1414  strcpy( punch.chPunch[punch.npunch], "HYDi" );
1415  sprintf( chHeader,
1416  "#hion\tzn\tgam1\tcoll ion1\tRecTot\tHRecCaB\thii/hi\tSim hii/hi"
1417  "\time_Hrecom_long(esc)\tdec2grd\texc pht\texc col\trec eff\tsec ion\n" );
1418  }
1419  else if( nMatch("POPU",chCard) )
1420  {
1421  /* punch hydrogen populations */
1422  strcpy( punch.chPunch[punch.npunch], "HYDp" );
1423  sprintf( chHeader,
1424  "#depth\tn(H0)\tn(H+)\tn(1s)\tn(2s)\tn(2p)\tetc\n" );
1425  }
1426  else if( nMatch("LINE",chCard) )
1427  {
1428  /* punch hydrogen lines
1429  * hydrogen line intensities and optical depths */
1430  strcpy( punch.chPunch[punch.npunch], "HYDl" );
1431  sprintf( chHeader,
1432  "#nup\tnlo\tE(ryd)\ttau\n" );
1433  }
1434  else if( nMatch(" LYA",chCard) )
1435  {
1436  /* punch hydrogen Lya some details about Lyman alpha */
1437  strcpy( punch.chPunch[punch.npunch], "HYDL" );
1438  sprintf( chHeader,
1439  "#depth\tTauIn\tTauTot\tn(2p)/n(1s)\tTexc\tTe\tTex/T\tPesc\tPdes\tpump\topacity\talbedo\n" );
1440  }
1441  else
1442  {
1443  fprintf( ioQQQ, "Punch hydrogen has options: CONDitions, 21 CM, LINE, POPUlations, and IONIzation.\nSorry.\n" );
1444  cdEXIT(EXIT_FAILURE);
1445  }
1446  }
1447 
1448  else if( nMatch("IONI",chCard) )
1449  {
1450  if( nMatch("RATE",chCard) )
1451  {
1452  /* punch ionization rates, search for the name of an element */
1453  if( (nelem = GetElem( chCard ) ) < 0 )
1454  {
1455  fprintf( ioQQQ, "There must be an element name on the ionization rates command. Sorry.\n" );
1456  cdEXIT(EXIT_FAILURE);
1457  }
1458  punch.punarg[punch.npunch][0] = (realnum)nelem;
1459  strcpy( punch.chPunch[punch.npunch], "IONR" );
1460  sprintf( chHeader,
1461  "#%s depth\teden\tdynamics.Rate\tabund\tTotIonize\tTotRecom\tSource\t ... \n",
1462  elementnames.chElementSym[nelem]);
1463  }
1464  else
1465  {
1466  /* punch table giving ionization means */
1467  strcpy( punch.chPunch[punch.npunch], "IONI" );
1468  sprintf( chHeader,
1469  "#Mean ionization distribution\n" );
1470  }
1471  }
1472 
1473  else if( nMatch(" IP ",chCard) )
1474  {
1475  strcpy( punch.chPunch[punch.npunch], " IP " );
1476  sprintf( chHeader,
1477  "#ionization potentials, valence shell\n" );
1478  }
1479 
1480  else if( nMatch("LEID",chCard) )
1481  {
1482  if( nMatch( "LINE" , chCard ) )
1483  {
1484  /* punch Leiden lines
1485  * final intensities of the Leiden PDR models */
1486  strcpy( punch.chPunch[punch.npunch], "LEIL" );
1487  sprintf( chHeader, "#ion\twl\tInt\trel int\n");
1488  }
1489  else
1490  {
1491  /* punch Leiden structure
1492  * structure of the Leiden PDR models */
1493  strcpy( punch.chPunch[punch.npunch], "LEIS" );
1494  sprintf( chHeader,
1495  /* 1-17 */
1496  "#Leid depth\tA_V(extentd)\tA_V(point)\tTe\tH0\tH2\tCo\tC+\tOo\tCO\tO2\tCH\tOH\te\tHe+\tH+\tH3+\t"
1497  /* 18 - 30 */
1498  "N(H0)\tN(H2)\tN(Co)\tN(C+)\tN(Oo)\tN(CO)\tN(O2)\tN(CH)\t(OH)\tN(e)\tN(He+)\tN(H+)\tN(H3+)\t"
1499  /* 31 - 32 */
1500  "H2(Sol)\tH2(FrmGrn)\t"
1501  /* 33 - 46*/
1502  "G0(DB96)\trate(CO)\trate(C)\theat\tcool\tGrnP\tGr-Gas-Cool\tGr-Gas-Heat\tCOds\tH2dH\tH2vH\tChaT\tCR H\tMgI\tSI\t"
1503  "Si\tFe\tNa\tAl\tC\tC610\tC370\tC157\tC63\tC146\n" );
1504  }
1505  }
1506 
1507  else if( nMatch("LINE",chCard) && nMatch("LIST",chCard) )
1508  {
1509  /* punch line list "output file" "Line List file" */
1510  long int j;
1511  /*
1512  * we parsed off the second file name at start of this routine
1513  * check if file was found, use it if it was, else abort
1514  */
1515  if( !lgSecondFilename )
1516  {
1517  fprintf(ioQQQ , "There must be a second file name between "
1518  "double quotes on the PUNCH LINE LIST command. This second"
1519  " file contains the input line list. I did not find it.\nSorry.\n");
1520  cdEXIT(EXIT_FAILURE);
1521  }
1522 
1523  /* actually get the lines, and malloc the space in the arrays
1524  * cdGetLineList will look on path */
1525  if( punch.ipPnunit[punch.npunch] == NULL &&
1526  ( punch.nLineList[punch.npunch] = cdGetLineList( chSecondFilename ,
1528  &punch.wlLineList[punch.npunch] ) ) < 0 )
1529  {
1530  fprintf(ioQQQ,"DISASTER could not open PUNCH LINE LIST file %s \n",
1531  chSecondFilename );
1532  cdEXIT(EXIT_FAILURE);
1533  }
1534 
1535  /* ratio option, in which pairs of lines form ratios, first over
1536  * second */
1537  if( nMatch("RATI",chCard) )
1538  {
1540  if( punch.nLineList[punch.npunch]%2 )
1541  {
1542  /* odd number of lines - cannot take ratio */
1543  fprintf(ioQQQ , "There must be an even number of lines to"
1544  " take ratios of lines. There were %li, an odd number."
1545  "\nSorry.\n", punch.nLineList[punch.npunch]);
1546  cdEXIT(EXIT_FAILURE);
1547  }
1548  }
1549  else
1550  {
1551  /* no ratio */
1552  punch.lgLineListRatio[punch.npunch] = false;
1553  }
1554 
1555  /* do special string saying our job, and then print labels */
1556  strcpy( punch.chPunch[punch.npunch], "LLST" );
1557 
1558  /* keyword absolute says to do absolute rather than relative intensities
1559  * relative intensities are the default */
1560  if( nMatch("ABSO",chCard) )
1561  {
1562  punch.punarg[punch.npunch][0] = 1;
1563  }
1564  else
1565  {
1566  punch.punarg[punch.npunch][0] = 0;
1567  }
1568 
1569  /* give header line */
1570  sprintf( chHeader, "#lineslist" );
1571  for( j=0; j<punch.nLineList[punch.npunch]; ++j )
1572  {
1573  /* if taking ratio then put div sign between pairs */
1574  if( punch.lgLineListRatio[punch.npunch] && is_odd(j) )
1575  strcat( chHeader , "/" );
1576  else
1577  strcat( chHeader , "\t" );
1578  sprintf( chTemp, "%s ", punch.chLineListLabel[punch.npunch][j] );
1579  strcat( chHeader, chTemp );
1580  sprt_wl( chTemp, punch.wlLineList[punch.npunch][j] );
1581  strcat( chHeader, chTemp );
1582  }
1583  strcat( chHeader, "\n" );
1584  }
1585 
1586  else if( nMatch("LINE",chCard) && !nMatch("XSPE",chCard) && !nMatch("NEAR",chCard))
1587  {
1588  /* punch line emissivity -
1589  * this is not punch xspec lines and not linear option
1590  * check for keyword UNITS on line, then scan wavelength or energy units,
1591  * sets punch.chConPunEnr*/
1592  ChkUnits(chCard);
1593 
1594  /* punch line emissivity, line intensity, line array,
1595  * and line data */
1596  if( nMatch("STRU",chCard) )
1597  {
1598  fprintf(ioQQQ," The PUNCH LINES STRUCTURE command is now PUNCH LINES "
1599  "EMISSIVITY.\n Sorry.\n\n");
1600  cdEXIT(EXIT_FAILURE);
1601  }
1602 
1603  else if( nMatch("EMIS",chCard) )
1604  {
1605  /* this used to be the punch lines structure command, is now
1606  * the punch lines emissivity command
1607  * give line emissivity vs depth */
1608  strcpy( punch.chPunch[punch.npunch], "LINS" );
1609  sprintf( chHeader,
1610  "#Emission line structure:");
1611  /* read in the list of lines to examine */
1612  punch_line(punch.ipPnunit[punch.npunch],"READ",false, chTemp);
1613  strcat( chHeader, chTemp );
1614  }
1615 
1616  else if( nMatch(" RT ", chCard ) )
1617  {
1618  /* punch line RT */
1619  strcpy( punch.chPunch[punch.npunch], "LINR" );
1620  /* punch some details needed for line radiative transfer
1621  * routine in punch_line.c */
1623  }
1624 
1625  else if( nMatch("CUMU",chCard) )
1626  {
1627  /* punch lines cumulative
1628  * this will be integrated line intensity, function of depth */
1629  strcpy( punch.chPunch[punch.npunch], "LINC" );
1630  /* option for either relative intensity or abs luminosity */
1631  if( nMatch("RELA",chCard) )
1632  {
1633  lgEOL = true;
1634  sprintf( chHeader,
1635  "#Cumulative emission line relative intensity.\n" );
1636  }
1637  else
1638  {
1639  sprintf( chHeader,
1640  "#Cumulative emission line absolute intensity.\n" );
1641  lgEOL = false;
1642  }
1643  /* read in the list of lines to examine */
1644  punch_line(punch.ipPnunit[punch.npunch],"READ",lgEOL,chTemp);
1645  strcat( chHeader, chTemp );
1646  }
1647 
1648  else if( nMatch("DATA",chCard) )
1649  {
1650  /* punch line data, done in PunchLineData */
1651  strcpy( punch.chPunch[punch.npunch], "LIND" );
1652  sprintf( chHeader,
1653  "#Emission line data.\n" );
1654  }
1655 
1656  else if( nMatch("ARRA",chCard) )
1657  {
1658  /* punch line array -
1659  * output energies and luminosities of predicted lines */
1660  strcpy( punch.chPunch[punch.npunch], "LINA" );
1661  sprintf( chHeader,
1662  "#enr\tID\tI(intrinsic)\tI(emergent)\ttype\n" );
1663  }
1664 
1665  else if( nMatch("LABE",chCard) )
1666  {
1667  /* punch line labels */
1668  strcpy( punch.chPunch[punch.npunch], "LINL" );
1669  sprintf( chHeader,
1670  "#index\tlabel\twavelength\tcomment\n" );
1671  /* this controls whether we will print lots of redundant
1672  * info labels for transferred lines - if keyword LONG appears
1673  * then do so, if does not appear then do not - this is default */
1674  if( nMatch("LONG",chCard) )
1675  punch.punarg[punch.npunch][0] = 1;
1676  else
1677  punch.punarg[punch.npunch][0] = 0;
1678  }
1679 
1680  else if( nMatch("OPTI",chCard) )
1681  {
1682  /* punch line optical depths, done in PunchLineStuff, located in punchdo.c */
1683  strcpy( punch.chPunch[punch.npunch], "LINO" );
1684  sprintf( chHeader,
1685  "#species\tenergy\topt depth\tdamp\n" );
1686 
1687  /* the default will be to make wavelengths line in the printout, called labels,
1688  * if units appears then other units will be used instead */
1689  strcpy( punch.chConPunEnr[punch.npunch],
1690  "labl" );
1691 
1692  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
1693  * units are copied into punch.chConPunEnr */
1694  if( nMatch("UNIT",chCard) )
1695  ChkUnits(chCard);
1696 
1697  /* this is optional limit to smallest optical depths */
1698  ipFFmt = 5;
1699  punch.punarg[punch.npunch][0] = (realnum)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL));
1700  /* this is default of 0.1 napier */
1701  if( lgEOL )
1702  {
1703  punch.punarg[punch.npunch][0] = 0.1f;
1704  }
1705  }
1706 
1707  else if( nMatch("POPU",chCard) )
1708  {
1709  /* punch line populations command - first give index and inforamtion
1710  * for all lines, then populations for lines as a function of
1711  * depth, using this index */
1712  strcpy( punch.chPunch[punch.npunch], "LINP" );
1713  sprintf( chHeader,
1714  "#population information\n" );
1715  /* this is optional limit to smallest population to punch - always
1716  * interpreted as a log */
1717  ipFFmt = 5;
1718  punch.punarg[punch.npunch][0] = (realnum)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL));
1719 
1720  /* this is default - all positive populations */
1721  if( lgEOL )
1722  punch.punarg[punch.npunch][0] = 0.f;
1723 
1724  if( nMatch(" OFF",chCard) )
1725  {
1726  /* no lower limit - print all lines */
1727  punch.punarg[punch.npunch][0] = -1.f;
1728  }
1729  }
1730 
1731  else if( nMatch("INTE",chCard) )
1732  {
1733  /* this will be full set of line intensities */
1734  strcpy( punch.chPunch[punch.npunch], "LINI" );
1735  sprintf( chHeader,
1736  "#Emission line intensities per unit inner area\n" );
1737  if( nMatch("COLU",chCard) )
1738  {
1739  /* column is key to punch single column */
1740  strcpy( punch.chPunRltType, "column" );
1741  }
1742  else
1743  {
1744  /* array is key to punch large array */
1745  strcpy( punch.chPunRltType, "array " );
1746  }
1747  if( nMatch("EVER",chCard) )
1748  {
1749  ipFFmt = 3;
1750  punch.LinEvery = (long int)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
1751  punch.lgLinEvery = true;
1752  if( lgEOL )
1753  {
1754  fprintf( ioQQQ,
1755  "There must be a second number, the number of zones to print.\nSorry.\n" );
1756  cdEXIT(EXIT_FAILURE);
1757  }
1758  }
1759  else
1760  {
1761  punch.LinEvery = geometry.nend[0];
1762  punch.lgLinEvery = false;
1763  }
1764  }
1765  else
1766  {
1767  fprintf( ioQQQ,
1768  "This option for PUNCH LINE is something that I do not understand. Sorry.\n" );
1769  cdEXIT(EXIT_FAILURE);
1770  }
1771  }
1772 
1773  else if( nMatch(" MAP",chCard) )
1774  {
1775  strcpy( punch.chPunch[punch.npunch], "MAP " );
1776  sprintf( chHeader,
1777  "#te, heating, cooling.\n" );
1778  /* do cooling space map for specified zones
1779  * if no number, or <0, do map and punch out without doing first zone
1780  * does map by calling punt(" map")
1781  */
1782  i = 5;
1783  hcmap.MapZone = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1784  if( lgEOL )
1785  {
1786  hcmap.MapZone = 1;
1787  }
1788 
1789  if( nMatch("RANG",chCard) )
1790  {
1791  bool lgLogOn;
1792  hcmap.RangeMap[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1793  if( hcmap.RangeMap[0] <= 10. && !nMatch("LINE",chCard) )
1794  {
1795  hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,hcmap.RangeMap[0]);
1796  lgLogOn = true;
1797  }
1798  else
1799  {
1800  lgLogOn = false;
1801  }
1802 
1803  hcmap.RangeMap[1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
1804  if( lgLogOn )
1805  hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,hcmap.RangeMap[1]);
1806 
1807  if( lgEOL )
1808  {
1809  fprintf( ioQQQ, "There must be a zone number, followed by two temperatures, on this line. Sorry.\n" );
1810  cdEXIT(EXIT_FAILURE);
1811  }
1812  }
1813  }
1814 
1815  else if( nMatch("MOLE",chCard) )
1816  {
1817  /* molecules, especially for PDR calculations */
1818  strcpy( punch.chPunch[punch.npunch], "MOLE" );
1819  sprintf( chHeader,
1820  "#molecular species will follow:\n");
1821  }
1822 
1823  else if( nMatch("OPTICAL",chCard) && nMatch("DEPTH",chCard) )
1824  {
1825  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
1826  * units are copied into punch.chConPunEnr */
1827  ChkUnits(chCard);
1828  if( nMatch("FINE",chCard) )
1829  {
1830  /* punch fine continuum optical depths */
1831  rfield.lgPunchOpacityFine = true;
1832  strcpy( punch.chPunch[punch.npunch], "OPTf" );
1833  sprintf( chHeader, "#energy\tTau tot\topacity\n" );
1834  ipFFmt = 5;
1835  /* range option - important since so much data */
1836  if( nMatch("RANGE",chCard) )
1837  {
1838  /* get lower and upper range, must be in Ryd */
1839  punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
1840  punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL);
1841  if( lgEOL )
1842  {
1843  fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n");
1844  cdEXIT(EXIT_FAILURE);
1845  }
1846  if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] )
1847  {
1848  fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n");
1849  cdEXIT(EXIT_FAILURE);
1850  }
1851  }
1852  else
1853  {
1854  /* these mean full energy range */
1855  punch.punarg[punch.npunch][0] = 0.;
1856  punch.punarg[punch.npunch][1] = 0.;
1857  }
1858  /* optional last parameter - how many points to bring together */
1859  punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt,
1860  INPUT_LINE_LENGTH,&lgEOL);
1861  /* default is to bring together ten */
1862  if( lgEOL )
1863  punch.punarg[punch.npunch][2] = 10;
1864  if( punch.punarg[punch.npunch][2] < 1 )
1865  {
1866  fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n");
1867  cdEXIT(EXIT_FAILURE);
1868  }
1869  }
1870  else
1871  {
1872  /* punch coarse continuum optical depths */
1873  strcpy( punch.chPunch[punch.npunch], "OPTc" );
1874  sprintf( chHeader,
1875  "#energy\ttotal\tabsorp\tscat\n" );
1876  }
1877 
1878  }
1879  else if( nMatch(" OTS",chCard) )
1880  {
1881  strcpy( punch.chPunch[punch.npunch], " OTS" );
1882  sprintf( chHeader,
1883  "#otscon, lin, conOpac LinOpc\n" );
1884  }
1885 
1886  else if( nMatch("OVER",chCard) && nMatch(" OVE",chCard) )
1887  {
1888  /* punch overview of model results */
1889  strcpy( punch.chPunch[punch.npunch], "OVER" );
1890  sprintf( chHeader,
1891  "#depth\tTe\tHtot\thden\teden\t2H_2/H\tHI\tHII\tHeI\tHeII\tHeIII\tCO/C\tC1\tC2\tC3\tC4\tO1\tO2\tO3\tO4\tO5\tO6\tAV(point)\tAV(extend)\n" );
1892  }
1893 
1894  else if( nMatch(" PDR",chCard) )
1895  {
1896  strcpy( punch.chPunch[punch.npunch], " PDR" );
1897  sprintf( chHeader,
1898  "#depth\tH colden\tTe\tHI/HDEN\tH2/HDEN\tH2*/HDEN\tCI/C\tCO/C\tH2O/O\tG0\tAV(point)\tAV(extend)\tTauV(point)\n" );
1899  }
1900 
1901  else if( nMatch("PHYS",chCard) )
1902  {
1903  /* punch physical conditions */
1904  strcpy( punch.chPunch[punch.npunch], "PHYS" );
1905  sprintf( chHeader,
1906  "#PhyC depth\tTe\tn(H)\tn(e)\tHtot\taccel\tfillfac\n" );
1907  }
1908 
1909  else if( nMatch("POIN",chCard) )
1910  {
1911  /* punch out the pointers */
1912  punch.lgPunPoint = true;
1913  /* this does not count as a punch option (really) */
1914  strcpy( punch.chPunch[punch.npunch], "" );
1915  punch.lgRealPunch[punch.npunch] = false;
1916  }
1917 
1918  else if( nMatch("PRES",chCard) )
1919  {
1920  /* the punch pressure command */
1921  strcpy( punch.chPunch[punch.npunch], "PRES" );
1922  sprintf( chHeader,
1923  "#P depth\tPcorrect\tPcurrent\tPIn+Pinteg\tPgas(r0)\tPgas\tPram\tPrad(line)\tPinteg\tV(wind km/s)\tcad(wind km/s)\tP(mag)\tV(turb km/s)\tP(turb)\tconv?\n" );
1924  }
1925 
1926  else if( nMatch("RADI",chCard) )
1927  {
1928  /* the punch radius command */
1929  sprintf( chHeader, "#NZONE\tradius\tdepth\tdr\n" );
1930  /* option to only punch the outer radius */
1931  if( nMatch( "OUTE" , chCard ) )
1932  {
1933  /* only outer radius */
1934  strcpy( punch.chPunch[punch.npunch], "RADO" );
1935  }
1936  else
1937  {
1938  /* all radii */
1939  strcpy( punch.chPunch[punch.npunch], "RADI" );
1940  }
1941  }
1942 
1943  else if( nMatch("RECO",chCard) )
1944  {
1945  if( nMatch("COEF",chCard) )
1946  {
1947  /* recombination coefficients for everything */
1948 
1949  /* this is logical flag used in routine ion_recom to create the punch output */
1950  punch.lgioRecom = true;
1951  /* this does not count as a punch option (really) */
1952  strcpy( punch.chPunch[punch.npunch], "" );
1953  punch.lgRealPunch[punch.npunch] = false;
1954  }
1955 
1956  else if( nMatch("EFFI",chCard) )
1957  {
1958  /* punch recombination efficiency */
1959  strcpy( punch.chPunch[punch.npunch], "RECE" );
1960  sprintf( chHeader,
1961  "#Recom effic H, Heo, He+\n" );
1962  }
1963 
1964  else
1965  {
1966  fprintf( ioQQQ, "No option recognized on this punch recombination command\n" );
1967  fprintf( ioQQQ, "Valid options are COEFFICIENTS, AGN, and EFFICIENCY\nSorry.\n" );
1968  cdEXIT(EXIT_FAILURE);
1969  }
1970  }
1971 
1972  /* punch results command, either as single column or wide array */
1973  else if( nMatch("RESU",chCard) )
1974  {
1975  strcpy( punch.chPunch[punch.npunch], "RESU" );
1976  if( nMatch("COLU",chCard) )
1977  {
1978  /* column is key to punch single column */
1979  strcpy( punch.chPunRltType, "column" );
1980  }
1981  else
1982  {
1983  /* array is key to punch large array */
1984  strcpy( punch.chPunRltType, "array " );
1985  }
1986 
1987  /* do not change following, is used as flag in getlines */
1988  sprintf( chHeader,
1989  "#results of calculation\n" );
1990  }
1991 
1992  else if( nMatch("SECO",chCard) )
1993  {
1994  /* punch secondary ionization rate */
1995  strcpy( punch.chPunch[punch.npunch], "SECO" );
1996  sprintf( chHeader,
1997  "#depth\tIon(H^0)\tDiss(H_2)\tExcit(Lya)\n" );
1998  }
1999 
2000  else if( nMatch("SOUR",chCard) )
2001  {
2002  if( nMatch("DEPT",chCard) )
2003  {
2004  /* print continuum source function as function of depth */
2005  strcpy( punch.chPunch[punch.npunch], "SOUD" );
2006  sprintf( chHeader,
2007  "#continuum source function vs depth\n" );
2008  }
2009  else if( nMatch("SPEC",chCard) )
2010  {
2011  /* print spectrum continuum source function at 1 depth */
2012  strcpy( punch.chPunch[punch.npunch], "SOUS" );
2013  sprintf( chHeader,
2014  "#continuum source function\n" );
2015  }
2016  else
2017  {
2018  fprintf( ioQQQ, "A second keyword must appear on this line.\n" );
2019  fprintf( ioQQQ, "They are DEPTH and SPECTRUM.\n" );
2020  fprintf( ioQQQ, "Sorry.\n" );
2021  cdEXIT(EXIT_FAILURE);
2022  }
2023  }
2024 
2025 
2026  /* punch spectrum the new form of the punch continuum, will eventually replace the standard
2027  * punch continuum command */
2028  else if( nMatch("SPEC",chCard) && nMatch("TRUM",chCard) && !nMatch("XSPE",chCard) )
2029  {
2030  /* this flag is checked in PrtComment to generate a caution
2031  * if continuum is punched but iterations not performed */
2032  punch.lgPunContinuum = true;
2033 
2034  /* set flag for spectrum */
2035  strcpy( punch.chPunch[punch.npunch], "CONN" );
2036 
2037  sprintf( chHeader,
2038  "#Cont Enr\tincid nFn\ttrans\tdiff\tlines \n" );
2039 
2040  /* check for keyword UNITS on line, then scan wavelength or energy units if present,
2041  * units are copied into punch.chConPunEnr */
2042  ChkUnits(chCard);
2043 
2044  /* this points to which punch new continuum this is -
2045  * parameters were stored in previous set spectrum commands */
2047 
2048  ++punch.cp_npun;
2049  /* check that limit not exceeded */
2050  if( punch.cp_npun > LIMPUN )
2051  {
2052  fprintf( ioQQQ,
2053  "The limit to the number of PUNCH options is %ld. Increase LIMPUN in punch.h if more are needed.\nSorry.\n",
2054  LIMPUN );
2055  cdEXIT(EXIT_FAILURE);
2056  }
2057 
2058  }
2059 
2060  else if( nMatch("SPEC",chCard) && nMatch("CIAL",chCard) )
2061  {
2062  /* punch special, will call routine PunchSpecial */
2063  strcpy( punch.chPunch[punch.npunch], "SPEC" );
2064  sprintf( chHeader, "#Special.\n" );
2065  }
2066 
2067  else if( nMatch("TEGR",chCard) )
2068  {
2069  /* punch tegrid */
2070  strcpy( punch.chPunch[punch.npunch], "TEGR" );
2071  sprintf( chHeader,
2072  "#zone, te, heat, cool.\n" );
2073  }
2074 
2075  else if( nMatch("TEMP",chCard) )
2076  {
2077  /* punch temperature command */
2078  strcpy( punch.chPunch[punch.npunch], "TEMP" );
2079  sprintf( chHeader,
2080  "#depth\tTe\tcC/dT\tdt/dr\td^2T/dr^2\n" );
2081  }
2082 
2083  else if( nMatch("TIME",chCard) && nMatch("DEPE",chCard) )
2084  {
2085  /* information about time dependent solutions */
2086  strcpy( punch.chPunch[punch.npunch], "TIMD" );
2087  /* do not want to separate iterations with special character */
2089  /* write header */
2090  sprintf( chHeader ,
2091  "#elapsed time\ttime step \tscale cont\tn(H)\t<T>\t<H+/H rad>\t<H0/H rad>\t<H2/H rad>\t<He+/H rad>\t<CO/H>\n" );
2092  }
2093 
2094  else if( nMatch("TPRE",chCard) )
2095  {
2096  /* debug output from the temperature predictor in zonestart,
2097  * set with punch tpred command */
2098  strcpy( punch.chPunch[punch.npunch], "TPRE" );
2099  sprintf( chHeader,
2100  "#zone old temp, guess Tnew, new temp delta \n" );
2101  }
2102 
2103  else if( nMatch("WIND",chCard) )
2104  {
2105  strcpy( punch.chPunch[punch.npunch], "WIND" );
2106  sprintf( chHeader,
2107  "#radius\tdepth\tvel [cm/s]\tTot accel [cm s-2]\tLin accel [cm s-2]"
2108  "\tCon accel [cm s-2]\tforce multiplier\ta_gravity\n" );
2109  if( nMatch( "TERM" , chCard ) )
2110  {
2111  /* only punch for last zone, the terminal velocity, for grids */
2112  punch.punarg[punch.npunch][0] = 0.;
2113  }
2114  else
2115  {
2116  /* one means punch every zone */
2117  punch.punarg[punch.npunch][0] = 1.;
2118  }
2119  }
2120 
2121  else if( nMatch("XSPE",chCard) )
2122  {
2123  /* say that this is a FITS file output */
2124  punch.lgFITS[punch.npunch] = true;
2125 
2126  /* the punch xspec commands */
2127  if( !punch.lgPunLstIter[punch.npunch] )
2128  {
2129  punch.lgPunLstIter[punch.npunch] = true;
2130  }
2131 
2132  if( nMatch("ATAB",chCard) )
2133  {
2134  /* punch xspec atable command */
2135  if( nMatch("INCI",chCard) )
2136  {
2137  if( nMatch("ATTE",chCard) )
2138  {
2139  /* attenuated incident continuum */
2140  strcpy( punch.chPunch[punch.npunch], "XATT" );
2141  grid.lgOutputTypeOn[2] = true;
2142  }
2143  else if( nMatch("REFL",chCard) )
2144  {
2145  /* reflected incident continuum */
2146  strcpy( punch.chPunch[punch.npunch], "XRFI" );
2147  grid.lgOutputTypeOn[3] = true;
2148  }
2149  else
2150  {
2151  /* incident continuum */
2152  strcpy( punch.chPunch[punch.npunch], "XINC" );
2153  grid.lgOutputTypeOn[1] = true;
2154  }
2155  }
2156  else if( nMatch("DIFF",chCard) )
2157  {
2158  if( nMatch("REFL",chCard) )
2159  {
2160  /* reflected diffuse continuous emission */
2161  strcpy( punch.chPunch[punch.npunch], "XDFR" );
2162  grid.lgOutputTypeOn[5] = true;
2163  }
2164  else
2165  {
2166  /* diffuse continuous emission outward */
2167  strcpy( punch.chPunch[punch.npunch], "XDFO" );
2168  grid.lgOutputTypeOn[4] = true;
2169  }
2170  }
2171  else if( nMatch("LINE",chCard) )
2172  {
2173  if( nMatch("REFL",chCard) )
2174  {
2175  /* reflected lines */
2176  strcpy( punch.chPunch[punch.npunch], "XLNR" );
2177  grid.lgOutputTypeOn[7] = true;
2178  }
2179  else
2180  {
2181  /* outward lines */
2182  strcpy( punch.chPunch[punch.npunch], "XLNO" );
2183  grid.lgOutputTypeOn[6] = true;
2184  }
2185  }
2186  else if( nMatch("SPEC",chCard) )
2187  {
2188  if( nMatch("REFL",chCard) )
2189  {
2190  /* reflected spectrum */
2191  strcpy( punch.chPunch[punch.npunch], "XREF" );
2192  grid.lgOutputTypeOn[9] = true;
2193  }
2194  else
2195  {
2196  /* transmitted spectrum */
2197  strcpy( punch.chPunch[punch.npunch], "XTRN" );
2198  grid.lgOutputTypeOn[8] = true;
2199  }
2200  }
2201  else
2202  {
2203  /* transmitted spectrum */
2204  strcpy( punch.chPunch[punch.npunch], "XTRN" );
2205  grid.lgOutputTypeOn[8] = true;
2206  }
2207  }
2208  else if( nMatch("MTAB",chCard) )
2209  {
2210  /* punch xspec mtable */
2211  strcpy( punch.chPunch[punch.npunch], "XSPM" );
2212  grid.lgOutputTypeOn[10] = true;
2213  }
2214  else
2215  {
2216  fprintf( ioQQQ, "Support only for xspec atable and xspec mtable.\n" );
2217  cdEXIT( EXIT_FAILURE );
2218  }
2219  }
2220 
2221  /* punch column density has to come last so do not trigger specific column
2222  * densities, H2, FeII, etc.
2223  * Need both keywords since column is also the keyword for one line per line */
2224  else if( nMatch("COLU",chCard) && nMatch("DENS",chCard) )
2225  {
2226  if( nMatch("SOME" , chCard ))
2227  {
2228  /* flag saying punch some column densities */
2229  strcpy( punch.chPunch[punch.npunch], "COLS" );
2230  punch_colden( chHeader , punch.ipPnunit[punch.npunch] , "READ" );
2231  }
2232  else
2233  {
2234  /* punch them all as in table */
2235  /* punch column density */
2236  strcpy( punch.chPunch[punch.npunch], "COLU" );
2237  sprintf( chHeader,
2238  "#column densities\n" );
2239  }
2240  }
2241  else
2242  {
2243  fprintf( ioQQQ,
2244  "ParsePunch cannot find a recognized keyword on this PUNCH command line.\nSorry.\n" );
2245  cdEXIT(EXIT_FAILURE);
2246  }
2247 
2248  /* only open if file has not already been opened during a previous call */
2249  if( punch.ipPnunit[punch.npunch] == NULL )
2250  {
2251  if( nMatch("XSPE",chCard) || nMatch("FITS",chCard) )
2252  {
2253  /* open the file with the name generated above */
2254  punch.ipPnunit[punch.npunch] = open_data( chFilename, "wb", AS_LOCAL_ONLY );
2255  }
2257  {
2258  char chIndex[4];
2259 
2260  /* take max because optimize.nOptimiz has not yet been incremented on first pass. */
2261  sprintf( chIndex, "%li", MAX2( 1, optimize.nOptimiz ) );
2262 
2263  strcat( chFilename, chIndex );
2264 
2265  fixit();
2266  /* must create unique name here. considerations include:
2267  * 1. must zero fill to the left
2268  * 2. must make sure that we have enough digits (don't hard-wire 3 or 4 digits,
2269  count the number of terms.
2270  * 3. should we look for a file ending (i.e., ".txt" and insert number before that? */
2271 
2272  /* open the file with the name generated above */
2273  punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY );
2274  }
2275  else
2276  {
2277  /* open the file with the name generated above */
2278  punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY );
2279  }
2280 
2281  /* option to set no buffering for this file. The setbuf command may
2282  * ONLY be issued right after the open of the file. Giving it after
2283  * i/o has been done may result in loss of the contents of the buffer, PvH */
2284  if( nMatch("NO BUFFER",chCard) )
2285  setbuf( punch.ipPnunit[punch.npunch] , NULL );
2286  }
2287 
2288  /***************************************************************/
2289  /* */
2290  /* The following are special punch options and must be done */
2291  /* after the parsing and file opening above. */
2292  /* */
2293  /* NB: these are ALSO parsed above. Here we DO something. */
2294  /* */
2295  /***************************************************************/
2296 
2297  if( nMatch("CONV",chCard) && nMatch("REAS",chCard) )
2298  {
2299  /* punch reason model declared not converged
2300  * not a true punch command, since done elsewhere */
2303  punch.lgPunConv = true;
2304  fprintf( punch.ipPunConv,
2305  "# reason for continued iterations\n" );
2306  strcpy( punch.chPunch[punch.npunch], "" );
2307  punch.lgRealPunch[punch.npunch] = false;
2308  }
2309 
2310  else if( nMatch("CONV",chCard) && nMatch("BASE",chCard) )
2311  {
2312  /* punch some quantities we are converging */
2313  punch.lgTraceConvergeBase = true;
2314  /* the second punch occurrence - file has been opened,
2315  * copy handle, also pass on special no hash option */
2316  if( nMatch("O HA",chCard) )
2319  /* set punch last flag to whatever it was above */
2321  static bool lgPrtHeader = true;
2322  if( lgPrtHeader )
2323  fprintf( punch.ipTraceConvergeBase,
2324  "#zone\theat\tcool\teden\n" );
2325  lgPrtHeader = false;
2326  }
2327 
2328  else if( nMatch(" DR ",chCard) )
2329  {
2330  static bool lgPrtHeader = true;
2331  /* the second punch dr occurrence - file has been opened,
2332  * copy handle to ipDRout, also pass on special no hash option */
2333  if( nMatch("O HA",chCard) )
2334  punch.lgDRHash = false;
2336  /* set punch last flag to whatever it was above */
2339  if( lgPrtHeader )
2340  fprintf( punch.ipDRout,
2341  "#zone\tdepth\tdr\tdr 2 go\treason \n" );
2342  lgPrtHeader = false;
2343  strcpy( punch.chPunch[punch.npunch], "" );
2344  punch.lgRealPunch[punch.npunch] = false;
2345  }
2346 
2347  else if( nMatch("QHEA",chCard) )
2348  {
2352  fprintf( gv.QHPunchFile,
2353  "#Probability distributions from quantum heating routine.\n" );
2354  }
2355 
2356  else if( nMatch("POIN",chCard) )
2357  {
2358  /* punch out the pointers */
2361  punch.lgPunPoint = true;
2362  fprintf( punch.ipPoint,
2363  "#pointers. \n" );
2364  strcpy( punch.chPunch[punch.npunch], "" );
2365  punch.lgRealPunch[punch.npunch] = false;
2366  }
2367 
2368  else if( nMatch("RECO",chCard) && nMatch("COEF",chCard) )
2369  {
2370  /* recombination coefficients for everything
2371  * punch.lgioRecom set to false in routine zero, non-zero value
2372  * is flag to punch recombination coefficients. the output is actually
2373  * produced by a series of routines, as they generate the recombination
2374  * coefficients. these include
2375  * diel supres, helium, hydrorecom, iibod, and makerecomb*/
2378  /* this is logical flag used in routine ion_recom to create the punch output */
2379  punch.lgioRecom = true;
2380  fprintf( punch.ioRecom,
2381  "#recombination coefficients cm3 s-1 for current density and temperature\n" );
2382  strcpy( punch.chPunch[punch.npunch], "" );
2383  punch.lgRealPunch[punch.npunch] = false;
2384  }
2385 
2386  else if( nMatch(" MAP",chCard) )
2387  {
2388  /* say output goes to special punch */
2390  }
2391 
2392  /* check that string written into chHeader can actually fit there
2393  * we may have overrun this buffer, an internal error */
2394  /* check that there are less than nChar characters in the line */
2395  char *chEOL = strchr(chHeader , '\0' );
2396 
2397  /* return null if input string longer than nChar, the longest we can read.
2398  * Print and return null but chLine still has as much of the line as
2399  * could be placed in cdLine */
2400  if( (chEOL==NULL) || (chEOL - chHeader)>=MAX_HEADER_SIZE-1 )
2401  {
2402  fprintf( ioQQQ, "DISASTER chHeader has been overwritten "
2403  "with a line too long to be read.\n" );
2404  cdEXIT(EXIT_FAILURE);
2405  }
2406 
2407  /* if flag set and cdHeader not still equal to initialized string, print header
2408  * punch.lgPunHeader normally true, is false during grid calculation */
2409  if( nSimThisCoreload > 1 )
2410  punch.lgPunHeader = false;
2411  if( punch.lgPunHeader && !nMatch("ArNdY38dZ9us4N4e12SEcuQ",chHeader) )
2412  {
2413  fprintf( punch.ipPnunit[punch.npunch], "%s", chHeader );
2414  }
2415 
2416  /* increment total number of punch commands, */
2417  ++punch.npunch;
2418  return;
2419 }
2420 
2421 /*PunchFilesInit initialize punch file pointers, called from InitCoreload
2422  * called one time per core load
2423  * NB KEEP THIS ROUTINE SYNCHED UP WITH THE NEXT ONE, ClosePunchFiles */
2424 void PunchFilesInit(void)
2425 {
2426  long int i;
2427  static bool lgFIRST = true;
2428 
2429  DEBUG_ENTRY( "PunchFilesInit()" );
2430 
2431  ASSERT( lgFIRST );
2432  lgFIRST = false;
2433 
2434  /* set lgNoClobber to not overwrite files, reset with clobber on punch line
2435  * if we are running a grid (grid command entered in cdRead) grid.lgGrid
2436  * true, is false if single sim. For grid we want to not clobber files
2437  * by default, do clobber for optimizer since this was behavior before */
2438  bool lgNoClobberDefault = false;
2439  if( grid.lgGrid )
2440  {
2441  /* cdRead encountered grid command - do not want to clobber files */
2442  lgNoClobberDefault = true;
2443  }
2444 
2445  for( i=0; i < LIMPUN; i++ )
2446  {
2447  lgNoClobber[i] = lgNoClobberDefault;
2448  }
2449  lgPunConv_noclobber = lgNoClobberDefault;
2450  lgDROn_noclobber = lgNoClobberDefault;
2451  lgTraceConvergeBase_noclobber = lgNoClobberDefault;
2452  lgPunPoint_noclobber = lgNoClobberDefault;
2453  lgioRecom_noclobber = lgNoClobberDefault;
2454  lgQHPunchFile_noclobber = lgNoClobberDefault;
2455 
2456  /* done parsing, now turn flag to punch headers, the default behavior
2457  * this is set false during a grid calculation to do only one header per grid*/
2458  punch.lgPunHeader = true;
2459 
2460  for( i=0; i < LIMPUN; i++ )
2461  {
2462  if( !lgNoClobber[i] )
2463  {
2464  punch.ipPnunit[i] = NULL;
2465  }
2466  /* set false with the dummy punch commands like punch dr */
2467  punch.lgRealPunch[i] = true;
2468  /* this sets energy range of continuum, zero says to do full continuum */
2469  punch.cp_range_min[i] = 0.;
2470  punch.cp_range_max[i] = 0.;
2471  /* this means to keep native resolution of code, reset to another
2472  * resolution with set punch resolution command */
2473  punch.cp_resolving_power[i] = 0.;
2474  }
2475 
2476  punch.lgTraceConvergeBase = false;
2477 
2478  if( !lgDROn_noclobber )
2479  {
2480  punch.ipDRout = NULL;
2481  punch.lgDROn = false;
2482  }
2483 
2485  {
2486  punch.ipTraceConvergeBase = NULL;
2487  punch.lgTraceConvergeBase = false;
2488  }
2489 
2490  if( !lgPunConv_noclobber )
2491  {
2492  punch.ipPunConv = NULL;
2493  punch.lgPunConv = false;
2494  }
2495 
2496  if( !lgPunPoint_noclobber )
2497  {
2498  punch.ipPoint = NULL;
2499  punch.lgPunPoint = false;
2500  }
2501 
2503  gv.QHPunchFile = NULL;
2504 
2505  if( !lgioRecom_noclobber )
2506  {
2507  punch.ioRecom = NULL;
2508  punch.lgioRecom = false;
2509  }
2510 
2511  ioMAP = NULL;
2512  return;
2513 }
2514 
2515 /*ClosePunchFiles close punch files called from cdEXIT upon termination,
2516  * from cloudy before returning
2517  * NB - KEEP THIS ROUTINE SYNCHED UP WITH THE PREVIOUS ONE, PunchFilesInit */
2518 void ClosePunchFiles( bool lgFinal )
2519 {
2520  long int i;
2521 
2522  DEBUG_ENTRY( "ClosePunchFiles()" );
2523 
2524  /* close all punch units cloudy opened with punch command,
2525  * lgNoClobber is set false with CLOBBER option on punch, says to
2526  * overwrite the files */
2527  for( i=0; i < punch.npunch; i++ )
2528  {
2529  /* if lgFinal is true, we close everything, no matter what.
2530  * this means ignoring "no clobber" options */
2531  if( punch.ipPnunit[i] != NULL && ( !lgNoClobber[i] || lgFinal ) )
2532  {
2533  /* Test that any FITS files are the right size! */
2534  if( punch.lgFITS[i] )
2535  {
2536  /* \todo 2 This overflows for file sizes larger (in bytes) than
2537  * a long int can represent (about 2GB on most 2007 systems) */
2538  fseek(punch.ipPnunit[i], 0, SEEK_END);
2539  long file_size = ftell(punch.ipPnunit[i]);
2540  if( file_size%2880 )
2541  {
2542  fprintf( ioQQQ, " PROBLEM FITS file is wrong size!\n" );
2543  }
2544  }
2545 
2546  fclose( punch.ipPnunit[i] );
2547  punch.ipPnunit[i] = NULL;
2548  }
2549  }
2550 
2551  /* following file handles are aliased to ipPnunit which
2552  * was closed above */
2553  if( punch.ipDRout != NULL && !lgDROn_noclobber )
2554  {
2555  /*fclose( punch.ipDRout );*/
2556  punch.ipDRout = NULL;
2557  punch.lgDROn = false;
2558  }
2559 
2561  {
2562  /*fclose( punch.ipDRout );*/
2563  punch.ipTraceConvergeBase = NULL;
2564  punch.lgTraceConvergeBase = false;
2565  }
2566 
2567  if( punch.ipPunConv != NULL && !lgPunConv_noclobber )
2568  {
2569  /*fclose( punch.ipPunConv );*/
2570  punch.ipPunConv = NULL;
2571  punch.lgPunConv = false;
2572  }
2573  if( punch.ipPoint != NULL && !lgPunPoint_noclobber )
2574  {
2575  /*fclose( punch.ipPoint );*/
2576  punch.ipPoint = NULL;
2577  punch.lgPunPoint = false;
2578  }
2579  if( gv.QHPunchFile != NULL && !lgQHPunchFile_noclobber )
2580  {
2581  /*fclose( gv.QHPunchFile );*/
2582  gv.QHPunchFile = NULL;
2583  }
2584  if( punch.ioRecom != NULL && !lgioRecom_noclobber )
2585  {
2586  /*fclose( punch.ioRecom );*/
2587  punch.ioRecom = NULL;
2588  punch.lgioRecom = false;
2589  }
2590  /* this file was already closed as a punch.ipPnunit, erase copy of pointer as well */
2591  ioMAP = NULL;
2592  return;
2593 }
2594 
2595 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present,
2596  * units are copied into punch.chConPunEnr - when doing output, the routine call
2597  * AnuUnit( energy ) will automatically return the energy in the right units,
2598  * when called to do punch output */
2599 STATIC void ChkUnits( char *chCard)
2600 {
2601 
2602  DEBUG_ENTRY( "ChkUnits()" );
2603 
2604  /* option to set units for continuum energy in punch output */
2605  if( nMatch("NITS",chCard) )
2606  {
2607  if( nMatch("MICR",chCard) )
2608  {
2609  /* microns */
2610  strcpy( punch.chConPunEnr[punch.npunch], "micr" );
2611  }
2612  else if( nMatch(" KEV",chCard) )
2613  {
2614  /* keV */
2615  strcpy( punch.chConPunEnr[punch.npunch], " kev" );
2616  }
2617  else if( nMatch("CENT",chCard) || nMatch(" CM ",chCard) )
2618  {
2619  /* centimeters*/
2620  strcpy( punch.chConPunEnr[punch.npunch], "cent" );
2621  }
2622  else if( nMatch(" EV ",chCard) )
2623  {
2624  /* eV */
2625  strcpy( punch.chConPunEnr[punch.npunch], " ev " );
2626  }
2627  else if( nMatch("ANGS",chCard) )
2628  {
2629  /* angstroms */
2630  strcpy( punch.chConPunEnr[punch.npunch], "angs" );
2631  }
2632  else if( nMatch("WAVE",chCard) )
2633  {
2634  /* wavenumbers */
2635  strcpy( punch.chConPunEnr[punch.npunch], "wave" );
2636  }
2637  else if( nMatch(" MHZ",chCard) )
2638  {
2639  /* megaherz */
2640  strcpy( punch.chConPunEnr[punch.npunch], " mhz" );
2641  }
2642  else if( nMatch(" RYD",chCard) )
2643  {
2644  /* the noble Rydberg */
2645  /* >>chng 01 sep 02, had been rdyb unlike proper ryd */
2646  strcpy( punch.chConPunEnr[punch.npunch], "ryd " );
2647  }
2648  else
2649  {
2650  fprintf( ioQQQ, "I did not recognize units on this line.\n" );
2651  fprintf( ioQQQ, "Units are keV, eV, Angstroms, Rydbergs, centimeters, and microns.\nSorry.\n" );
2652  cdEXIT(EXIT_FAILURE);
2653  }
2654  }
2655  else
2656  {
2657  strcpy( punch.chConPunEnr[punch.npunch], "ryd " );
2658  }
2659  return;
2660 }

Generated for cloudy by doxygen 1.8.1.2