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