cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ParsePunch parse the punch command */ 00004 /*PunchFilesInit initialize punch file pointers, called from cdInit */ 00005 /*ClosePunchFiles close all punch files */ 00006 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present */ 00007 #include "cddefines.h" 00008 #include "cddrive.h" 00009 #include "physconst.h" 00010 #include "elementnames.h" 00011 #include "input.h" 00012 #include "geometry.h" 00013 #include "prt.h" 00014 #include "optimize.h" 00015 #include "rfield.h" 00016 #include "hcmap.h" 00017 #include "atomfeii.h" 00018 #include "h2.h" 00019 #include "mole.h" 00020 #include "hmi.h" 00021 #include "version.h" 00022 #include "grainvar.h" 00023 #include "parse.h" 00024 #include "grid.h" 00025 #include "punch.h" 00026 /* check for keyword UNITS on line, then scan wavelength or energy units if present */ 00027 STATIC void ChkUnits( char *chCard); 00028 00029 /* option to append instead of overwrite - default was to overwrite, 00030 * >>chng 07 feb 02 invert logic, default not to overwrite, but clobber 00031 * keyword says to overwrite */ 00032 static bool lgNoClobber[LIMPUN]; 00033 00034 /* these are for some special cases, same purpose as previous no clobber */ 00035 static bool lgPunConv_noclobber , lgDROn_noclobber , 00036 lgPunPoint_noclobber , lgioRecom_noclobber , lgQHPunchFile_noclobber, 00037 lgTraceConvergeBase_noclobber; 00038 00039 /* NB NB NB NB NB NB NB NB NB NB 00040 * 00041 * if any "special" punch commands are added (commands that copy the file pointer 00042 * into a separate variable, e.g. like PUNCH _DR_), be sure to add that file pointer 00043 * to PunchFilesInit and ClosePunchFiles !!! 00044 * 00045 * PUNCH FILE POINTERS SHOULD NEVER BE ALTERED BY ROUTINES OUTSIDE THIS MODULE !!! 00046 * 00047 * hence initializations of punch file pointers should never be included in zero() !! 00048 * the pointer might be lost without the file being closed 00049 * 00050 * NB NB NB NB NB NB NB NB NB NB */ 00051 00052 /* punch file header headers - these are written into the string chHeader when 00053 * the command is parsed 00054 * punch.lgPunHeader determines whether header is punched 00055 */ 00056 00057 00058 void ParsePunch(char *chCard ) 00059 { 00060 char chLabel[INPUT_LINE_LENGTH] , 00061 chFilename[INPUT_LINE_LENGTH] , 00062 chSecondFilename[INPUT_LINE_LENGTH]; 00063 bool lgEOL, 00064 lgSecondFilename; 00065 /* pointer to column across line image for free format number reader*/ 00066 long int ipFFmt, 00067 i, 00068 nelem; 00069 00070 #define MAX_HEADER_SIZE 2000 00071 00072 char chHeader[MAX_HEADER_SIZE], 00073 chTemp[MAX_HEADER_SIZE]; 00074 00075 00076 DEBUG_ENTRY( "ParsePunch()" ); 00077 00078 /* initialize chHeader string with nonsense, compare at end to see if we have any actual headers. */ 00079 sprintf( chHeader, "ArNdY38dZ9us4N4e12SEcuQ" ); 00080 00081 /* check that limit not exceeded */ 00082 if( punch.npunch >= LIMPUN ) 00083 { 00084 fprintf( ioQQQ, 00085 "The limit to the number of PUNCH options is %ld. Increase " 00086 "LIMPUN in punch.h if more are needed.\nSorry.\n", 00087 LIMPUN ); 00088 cdEXIT(EXIT_FAILURE); 00089 } 00090 00091 /* if this is first call to this routine we will need to init some vars */ 00092 { 00093 static bool lgMustInit=true; 00094 if( lgMustInit ) 00095 { 00096 lgMustInit = false; 00097 00098 /* following are for punch LineList option 00099 * nLineList is number of em lines, -1 if not defined */ 00100 punch.nLineList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) ); 00101 /* lgLineListRatio flag saying to take ratios of pairs of lines */ 00102 punch.lgLineListRatio = (bool*)MALLOC((size_t)(LIMPUN*sizeof(bool)) ); 00103 /* wlLineList is set of emission lines for LineList */ 00104 punch.wlLineList = (realnum **)MALLOC((size_t)(LIMPUN*sizeof(realnum *)) ); 00105 /* chLineListLabel is label for line list */ 00106 punch.chLineListLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) ); 00107 /* remaining dimensions will be created if punch LineList command entered 00108 * but init nLineList to -1 to signal not set */ 00109 for( i=0; i<LIMPUN; ++i ) 00110 { 00111 punch.nLineList[i] = -1; 00112 punch.lgFITS[i] = false; 00113 } 00114 00115 /* following are parallel but for punch averages option 00116 * nAverageList is number of averages, -1 if not defined */ 00117 punch.nAverageList = (long int*)MALLOC((size_t)(LIMPUN*sizeof(long int)) ); 00118 /* chAverageSpeciesLabel is label for species of each average */ 00119 punch.chAverageSpeciesLabel = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) ); 00120 /* chAverageType is label for species of type of average */ 00121 punch.chAverageType = (char***)MALLOC((size_t)(LIMPUN*sizeof(char**)) ); 00122 /* nAverageIonList is set of ion stages for average */ 00123 punch.nAverageIonList = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) ); 00124 /* nAverage2ndPar is second parameters of each average */ 00125 punch.nAverage2ndPar = (int **)MALLOC((size_t)(LIMPUN*sizeof(int *)) ); 00126 /* remaining dimensions will be created if punch averages command entered 00127 * but init nLineList to -1 to signal not set */ 00128 for( i=0; i<LIMPUN; ++i ) 00129 { 00130 punch.nAverageList[i] = -1; 00131 } 00132 } 00133 } 00134 00135 /* initialize this flag to false, set true for special cases below */ 00136 punch.lgPunchToSeparateFiles[punch.npunch] = false; 00137 00138 /* LAST keyword is an option to punch only on last iteration */ 00139 if( nMatch("LAST",chCard) ) 00140 punch.lgPunLstIter[punch.npunch] = true; 00141 else 00142 punch.lgPunLstIter[punch.npunch] = false; 00143 00144 /* get file name for this punch output. 00145 * GetQuote does the following - 00146 * first copy original version of file name into chLabel, 00147 * string does include null termination. 00148 * set filename in OrgCard and second parameter to spaces so 00149 * that not picked up below as keyword 00150 * last parameter says whether to abort if no quote found */ 00151 if( GetQuote( chLabel , chCard , true ) ) 00152 /* this can't happen since routine would not return at all if no double quotes found */ 00153 TotalInsanity(); 00154 00155 /* check that name is not same as opacity.opc, a special file */ 00156 if( strcmp(chLabel , "opacity.opc") == 0 ) 00157 { 00158 fprintf( ioQQQ, "ParsePunch will not allow punch file name %s, please choose another.\nSorry.\n", 00159 chLabel); 00160 cdEXIT(EXIT_FAILURE); 00161 } 00162 else if( chLabel[0]==0 ) 00163 { 00164 fprintf( ioQQQ, "ParsePunch found a null file name between double quotes, please check command line.\nSorry.\n"); 00165 cdEXIT(EXIT_FAILURE); 00166 } 00167 00168 /* now copy to chFilename, with optional prefix first */ 00169 /* this is optional prefix, normally a null string, set with set punch prefix command */ 00170 strcpy( chFilename , punch.chFilenamePrefix ); 00171 strcat( chFilename , chLabel ); 00172 00173 /* there may be a second file name, and we need to get it off the line 00174 * before we parse options, last false parameter says not to abort if 00175 * missing - this is not a problem at this stage */ 00176 if( GetQuote( chSecondFilename , chCard , false ) ) 00177 lgSecondFilename = false; 00178 else 00179 lgSecondFilename = true; 00180 00181 /* CLOBBER clobber keyword is an option to overwrite rather than 00182 * append to a given file */ 00183 if( nMatch("CLOB",chCard) ) 00184 { 00185 if( nMatch(" NO ",chCard) ) 00186 { 00187 /* do not clobber files */ 00188 lgNoClobber[punch.npunch] = true; 00189 } 00190 else 00191 { 00192 /* clobber files */ 00193 lgNoClobber[punch.npunch] = false; 00194 } 00195 } 00196 00197 /* put version number and title of model on output file, but only if 00198 * this is requested with a "title" on the line*/ 00199 /* >>chng 02 may 10, invert logic from before - default had been title */ 00200 /* put version number and title of model on output file, but only if 00201 * there is not a "no title" on the line*/ 00202 if( !nMatch("O TI",chCard) && nMatch("TITL",chCard)) 00203 { 00204 sprintf( chHeader, 00205 "# %s %s\n", 00206 t_version::Inst().chVersion, input.chTitle ); 00207 } 00208 00209 /* usually results for each iteration are followed by a series of 00210 * hash marks, ####, which fool excel. This is an option to not print 00211 * the line. If the keyword NO HASH no hash appears the hash marks 00212 * will not occur */ 00213 if( nMatch("O HA",chCard) ) 00214 punch.lgHashEndIter[punch.npunch] = false; 00215 00216 /* punch opacity must come first since elements appear as sub-keywords */ 00217 if( nMatch("OPAC",chCard) ) 00218 { 00219 /* for grid run, punch results of different models to different files. */ 00220 punch.lgPunchToSeparateFiles[punch.npunch] = true; 00221 00222 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 00223 * units are copied into punch.chConPunEnr */ 00224 ChkUnits(chCard); 00225 00226 strcpy( punch.chPunch[punch.npunch], "OPAC" ); 00227 if( nMatch("TOTA",chCard) ) 00228 { 00229 /* DoPunch will call punch_opacity to parse the subcommands 00230 * punch total opacity */ 00231 strcpy( punch.chOpcTyp[punch.npunch], "TOTL" ); 00232 sprintf( chHeader, 00233 "#nu\tTot opac\tAbs opac\tScat opac\tAlbedo\telem\n" ); 00234 } 00235 00236 else if( nMatch("FIGU",chCard) ) 00237 { 00238 /* do figure for hazy */ 00239 strcpy( punch.chOpcTyp[punch.npunch], "FIGU" ); 00240 sprintf( chHeader, 00241 "#nu, H, He, tot opac\n" ); 00242 } 00243 00244 else if( nMatch("FINE",chCard) ) 00245 { 00246 /* punch the fine opacity array */ 00247 rfield.lgPunchOpacityFine = true; 00248 strcpy( punch.chOpcTyp[punch.npunch], "FINE" ); 00249 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 00250 * units are copied into punch.chConPunEnr */ 00251 ChkUnits(chCard); 00252 00253 sprintf( chHeader, 00254 "#nu\topac\n" ); 00255 ipFFmt = 5; 00256 /* range option - important since so much data - usually want to 00257 * only give portion of the continuum */ 00258 if( nMatch("RANGE",chCard) ) 00259 { 00260 /* get lower and upper range, must be in Ryd */ 00261 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 00262 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 00263 if( lgEOL ) 00264 { 00265 fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n"); 00266 cdEXIT(EXIT_FAILURE); 00267 } 00268 if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] ) 00269 { 00270 fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n"); 00271 cdEXIT(EXIT_FAILURE); 00272 } 00273 } 00274 else 00275 { 00276 /* these mean full energy range */ 00277 punch.punarg[punch.npunch][0] = 0.; 00278 punch.punarg[punch.npunch][1] = 0.; 00279 } 00280 /* optional last parameter - how many points to bring together */ 00281 punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt, 00282 INPUT_LINE_LENGTH,&lgEOL); 00283 00284 /* default is to average together ten */ 00285 if( lgEOL ) 00286 punch.punarg[punch.npunch][2] = 10; 00287 00288 if( punch.punarg[punch.npunch][2] < 1 ) 00289 { 00290 fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n"); 00291 cdEXIT(EXIT_FAILURE); 00292 } 00293 } 00294 00295 else if( nMatch("GRAI",chCard) ) 00296 { 00297 /* punch grain opacity command, give optical properties of gains in calculation */ 00298 strcpy( punch.chPunch[punch.npunch], "DUSO" ); 00299 /* punch grain opacity command in twice, here and above in opacity */ 00300 sprintf( chHeader, 00301 "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" ); 00302 } 00303 00304 else if( nMatch("BREM",chCard) ) 00305 { 00306 /* punch bremsstrahlung opacity */ 00307 strcpy( punch.chOpcTyp[punch.npunch], "BREM" ); 00308 sprintf( chHeader, 00309 "#nu\tbrem opac\n" ); 00310 } 00311 00312 else if( nMatch("SHEL",chCard) ) 00313 { 00314 /* punch shells, a form of the punch opacity command for showing subshell crossections*/ 00315 strcpy( punch.chPunch[punch.npunch], "OPAC" ); 00316 00317 /* punch subshell cross sections */ 00318 strcpy( punch.chOpcTyp[punch.npunch], "SHEL" ); 00319 00320 /* this is element */ 00321 ipFFmt = 3; 00322 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt, 00323 INPUT_LINE_LENGTH,&lgEOL); 00324 00325 /* this is ion */ 00326 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt, 00327 INPUT_LINE_LENGTH,&lgEOL); 00328 00329 /* this is shell */ 00330 punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt, 00331 INPUT_LINE_LENGTH,&lgEOL); 00332 00333 if( lgEOL ) 00334 { 00335 fprintf( ioQQQ, "There must be IO unit, atom weight, ion, shell\nSorry.\n" ); 00336 cdEXIT(EXIT_FAILURE); 00337 } 00338 sprintf( chHeader, 00339 "#sub shell cross section\n" ); 00340 } 00341 00342 else if( nMatch("ELEM",chCard) ) 00343 { 00344 /* punch element opacity, produces n name.n files, one for each stage of 00345 * ionization. the name is the 4-char version of the element's name, and 00346 * n is the stage of ionization. the file name on the card is ignored. 00347 * The code stops in punch_opacity after these files are produced. */ 00348 00349 /* this will be used as check that we did find an element on the command lines */ 00350 /* nelem is -1 if an element was not found */ 00351 if( (nelem = GetElem( chCard ) ) < 0 ) 00352 { 00353 fprintf( ioQQQ, "There must be a second keyword on the opacity command. Sorry.\n" ); 00354 fprintf( ioQQQ, "Options are total, figure, grain, or element name.\n" ); 00355 cdEXIT(EXIT_FAILURE); 00356 } 00357 00358 /* copy string over */ 00359 strcpy( punch.chOpcTyp[punch.npunch], elementnames.chElementNameShort[nelem] ); 00360 } 00361 else 00362 { 00363 fprintf( ioQQQ, " I did not recognize a keyword on this punch opacity command.\n" ); 00364 fprintf( ioQQQ, " Sorry.\n" ); 00365 cdEXIT(EXIT_FAILURE); 00366 } 00367 } 00368 00369 /* punch H2 has to come early since it has many suboptions */ 00370 else if( nMatch(" H2 ",chCard) ) 00371 { 00372 /* this is in mole_h2_io.c */ 00373 H2_ParsePunch( chCard , chHeader ); 00374 } 00375 00376 /* punch grain abundance will be handled later */ 00377 else if( nMatch("ABUN",chCard) && !nMatch("GRAI",chCard) ) 00378 { 00379 /* punch abundances */ 00380 strcpy( punch.chPunch[punch.npunch], "ABUN" ); 00381 sprintf( chHeader, 00382 "#abund H" ); 00383 for(nelem=ipHELIUM;nelem<LIMELM; ++nelem ) 00384 { 00385 sprintf( chTemp, "\t%s", 00386 elementnames.chElementNameShort[nelem] ); 00387 strcat( chHeader, chTemp ); 00388 } 00389 strcat( chHeader, "\n"); 00390 } 00391 00392 else if( nMatch(" AGE",chCard) ) 00393 { 00394 /* punch ages */ 00395 strcpy( punch.chPunch[punch.npunch], "AGES" ); 00396 sprintf( chHeader, 00397 "#ages depth\tt(cool)\tt(H2 dest)\tt(CO dest)\tt(OH dest)\tt(H rec)\n" ); 00398 } 00399 00400 else if( nMatch(" AGN",chCard) ) 00401 { 00402 /* punch tables needed for AGN3 */ 00403 strcpy( punch.chPunch[punch.npunch], " AGN" ); 00404 /* this is the AGN option, to produce a table for AGN */ 00405 00406 /* charge exchange rate coefficients */ 00407 if( nMatch("CHAR",chCard) ) 00408 { 00409 strcpy( punch.chPunch[punch.npunch], "CHAG" ); 00410 sprintf( chHeader, 00411 "#charge exchange rate coefficnt\n" ); 00412 } 00413 00414 else if( nMatch("RECO",chCard) ) 00415 { 00416 /* punch recombination rates for AGN3 table */ 00417 strcpy( punch.chPunch[punch.npunch], "RECA" ); 00418 sprintf( chHeader, 00419 "#Recom rates for AGN3 table\n" ); 00420 } 00421 00422 else if( nMatch("OPAC",chCard) ) 00423 { 00424 /* create table for appendix in AGN */ 00425 strcpy( punch.chOpcTyp[punch.npunch], " AGN" ); 00426 strcpy( punch.chPunch[punch.npunch], "OPAC" ); 00427 } 00428 00429 else if( nMatch("HECS",chCard) ) 00430 { 00431 /* create table for appendix in AGN */ 00432 strcpy( punch.chPunchArgs[punch.npunch], "HECS" ); 00433 sprintf( chHeader, 00434 "#AGN3 he cs \n" ); 00435 } 00436 00437 else if( nMatch("HEMI",chCard) ) 00438 { 00439 /* HEMIS - continuum emission needed for chap 4 of AGN3 */ 00440 strcpy( punch.chPunchArgs[punch.npunch], "HEMI" ); 00441 00442 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 00443 * units are copied into punch.chConPunEnr */ 00444 ChkUnits(chCard); 00445 } 00446 else if( nMatch("RECC",chCard) ) 00447 { 00448 /* recombination cooling, for AGN */ 00449 strcpy( punch.chPunch[punch.npunch], "HYDr" ); 00450 sprintf( chHeader, 00451 "#T\tbAS\tb1\tbB\n" ); 00452 } 00453 else 00454 { 00455 fprintf( ioQQQ, " I did not recognize this option on the PUNCH HYDROGEN command.\n" ); 00456 fprintf( ioQQQ, " Sorry.\n" ); 00457 cdEXIT(EXIT_FAILURE); 00458 } 00459 } 00460 00461 else if( nMatch("ASSE",chCard) ) 00462 { 00463 /* punch asserts */ 00464 strcpy( punch.chPunch[punch.npunch], "ASSE" ); 00465 /* no need to print this standard line of explanation*/ 00466 /*sprintf( chHeader, " asserts\n" );*/ 00467 } 00468 00469 else if( nMatch("AVER",chCard) ) 00470 { 00471 /* punch averages */ 00472 strcpy( punch.chPunch[punch.npunch], "AVER" ); 00473 /* no need to print this standard line of explanation*/ 00474 /*sprintf( chHeader, " asserts\n" );*/ 00475 00476 /* actually get the averages from the input stream, and malloc the 00477 * space in the arrays 00478 * punch io unit not used in read */ 00479 punch_average( punch.npunch, "READ", chHeader ); 00480 } 00481 00482 /* punch charge transfer */ 00483 else if( nMatch("CHAR",chCard) && nMatch("TRAN",chCard) ) 00484 { 00485 /* NB in PunchDo only the first three characters are compared to find this option, 00486 * search for "CHA" */ 00487 /* punch charge transfer */ 00488 strcpy( punch.chPunch[punch.npunch], "CHAR" ); 00489 sprintf( chHeader, 00490 "#charge exchange rate coefficient\n" ); 00491 } 00492 00493 else if( nMatch("CHEM",chCard) ) 00494 { 00495 00496 if( nMatch( "RATE" , chCard ) ) 00497 { 00498 /* >>chng 06 May 30, NPA. Punch reaction rates for selected species */ 00499 00500 if( nMatch( " CO " , chCard ) ) 00501 { 00502 00503 strcpy( punch.chPunch[punch.npunch], "RCO " ); 00504 sprintf( chHeader, 00505 "#Depth\tH2O_C3HP_CO_C2H3P\tC2H2_HCOP_CO_C2H3P\tC2H_HCOP_CO_C2H2P\tO_C3H_CO_C2H\t" 00506 "O_C2H2_CO_CH2\tOH_C2H2_CO_CH3\tHCOP_C3_C3HP_CO\tO2_C3_CO_C2_O\tO_C3_CO_C2\t" 00507 "H_COP_CO_HP\tHminus_HCOP_CO_H2\tH2P_CO_HCOP_H\tH2P_CO_COP_H2\tH3P_CO_HCOP_H2\t" 00508 "HeP_CO_OP_C_He\tHeP_CO_O_CP_He\tcrnu_CO_O_C\tCRP_CO_COP_e\tnu_CO_O_C\te_HCOP_CO_H\t" 00509 "C_COP_CO_CP\tC_HCOP_CO_CHP\tC_O_CO_nu\tC_O2_CO_O\tC_OH_CO_H\tC_SiOP_SiP_CO\tCP_OH_CO_HP\t" 00510 "CP_SiO_SiP_CO\tCP_O2_CO_OP\tO_CH_CO_H\tO_CH2_CO_H_H\tO_CH2_CO_H2\tO_COP_CO_OP\tOP_CO_COP_O\t" 00511 "CH_COP_CO_CHP\tCH_HCOP_CO_CH2P\tCH_O2_CO_OH\tCH2_COP_CO_CH2P\tCH2_HCOP_CO_CH3P\t" 00512 "CH2_O2_CO_H2O\tCOP_O2_O2P_CO\tH2O_COP_CO_H2OP\tH2O_HCOP_CO_H3OP\tH2OP_CO_HCOP_OH\t" 00513 "HCOP_SiH_SiH2P_CO\tHCOP_SiO_SiOHP_CO\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tOHP_CO_HCOP_O\t" 00514 "COP_CH4_CO_CH4P\tCO_CH4P_HCOP_CH3\tCO_CH5P_HCOP_CH4\tHP_OCS_HSP_CO\tHeP_OCS_SP_CO_He\t" 00515 "OCNP_e_CO_N\tOCSP_e_S_CO\tOCS_NU_S_CO\tOCS_CRP_S_CO\tC_NO_CO_N\tC_OCN_CO_CN\t" 00516 "C_SO_S_CO\tO_CN_CO_N\tO_HCN_CO_NH\tO_OCN_NO_CO\tO_CS_S_CO\tO_OCS_SO_CO\tOH_HCN_CO_NH2\t" 00517 "CN_NO_N2_CO\tCN_O2_NO_CO\tCO_HS_OCS_H\tCP_SO_SP_CO\tCP_OCS_CSP_CO\tCHP_OCS_HCSP_CO\t" 00518 "NP_CO_NOP_C\tNP_OCS_SP_CO_N\tNHP_CO_HCOP_N\tNHP_CO_OCNP_H\tNH_HCOP_CO_NH2P\t" 00519 "NH2_HCOP_CO_NH3P\tNH3_HCOP_CO_NH4P\tCNP_O2_NOP_CO\tHCNP_CO_HCOP_CN\tCO_HNOP_NO_HCOP\t" 00520 "N2P_OCS_SP_N2_CO\tHCOP_S_HSP_CO\tHCOP_CS_HCSP_CO\tNH_COP_CO_NHP\tNH2_COP_CO_NH2P\t" 00521 "NH3_COP_CO_NH3P\tCNP_CO_COP_CN\tHCN_COP_CO_HCNP\tCO_N2P_N2_COP\tCOP_NO_NOP_CO\t" 00522 "CO_S_OCS_NU\tNP_CO_COP_N\tCOP_S_SP_CO\tC_CO_C2_O\tO_C2_CO_C\tC2_O2_CO_CO\t" 00523 "C2_O2P_COP_CO\tC2_COP_CO_C2P\tC2P_O2_COP_CO\tO_C2H_CO_CH\tO_CCl_Cl_CO\t" 00524 "CO_H2ClP_HCl_HCOP\tHNC_HCOP_HCNHP_CO\tHCN_HCOP_HCNHP_CO\tC2H_COP_CO_C2HP\tC2_HCOP_CO_C2HP\n"); 00525 } 00526 00527 else if( nMatch( " OH " , chCard ) ) 00528 { 00529 strcpy( punch.chPunch[punch.npunch], "ROH " ); 00530 sprintf( chHeader, 00531 "#Depth\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tnu_H2O_OH_H\tnu_OH_OHP_e\tnu_OH_O_H\tcrnu_H2O_OH_H\tcrnu_OH_O_H\tCP_OH_CO_HP\tCP_OH_COP_H\t" 00532 "CP_OH_CO_HP\tCP_OH_COP_H\tCH2_OHP_OH_CH2P\tCH2P_O2_HCOP_OH\tH2O_COP_HCOP_OH\tH2OP_H2O_H3OP_OH\tC_H2OP_OH_CHP\tOP_OH_O2P_H\tOP_OH_OHP_O\tSiP_OH_SiOP_H\t" 00533 "CH_H2OP_OH_CH2P\tCH_OHP_OH_CHP\tCHP_OH_COP_H2\tCHP_O2_COP_OH\tCH2_H2OP_OH_CH3P\tOH_COP_HCOP_O\tOH_H2OP_H3OP_O\tOHP_H2O_H2OP_OH\tOHP_O2_O2P_OH\tOHP_OH_H2OP_O\t" 00534 "O_CH4P_OH_CH3P\tOP_CH4_OH_CH3P\tCH5P_OH_H2OP_CH4\tOHP_C2_C2P_OH\tOH_C2P_C2_OHP\tCH_NO_CN_OH\tNH_O_OH_N\tNH_OH_HNO_H\tO_HNO_NO_OH\tNH2_OH_H2O_NH\t" 00535 "NH2_NO_N2_OH_H\tOH_CN_OCN_H\tOH_S_HS_O\tOH_S_SO_H\tNHP_OH_H2OP_N\tNHP_H2O_OH_NH2P\tNHP_O2_NOP_OH\tO_HSP_SP_OH\tNH2P_OH_H2OP_NH\tNH2P_H2O_NH3P_OH\t" 00536 "NH2_H2OP_NH3P_OH\tNH2P_O2_HNOP_OH\tOH_NH3P_NH4P_O\tOH_HCNP_CN_H2OP\tOH_HNOP_NO_H2OP\tOH_SP_SOP_H\tNH3P_H2O_NH4P_OH\tNH3_H2OP_NH4P_OH\tH2O_CNP_HCNP_OH\tH2OP_S_HSP_OH\t" 00537 "NH_OHP_OH_NHP\tNH2_OHP_OH_NH2P\tOHP_NH3_NH3P_OH\tOH_CNP_CN_OHP\tOH_N2P_N2_OHP\tOHP_NO_NOP_OH\tNP_OH_OHP_N\tOHP_S_SP_OH\tH2OP_HNC_HCNHP_OH\tH2OP_HCN_HCNHP_OH\t" 00538 "H2O_C2P_C2HP_OH\tH2OP_C2_C2HP_OH\tH2OP_C2H_C2H2P_OH\tOHP_C2H_C2HP_OH\tHminus_OH_H2O_e\tHminus_O_OH_e\tHP_OH_OHP_H\tH2s_OH_O_H2_H\tH2s_H2O_OH_H2_H\tH2P_OH_H2OP_H\t" 00539 "H2P_OH_OHP_H2\tH3P_OH_H2OP_H2\tHeP_OH_OP_He_H\tHeP_H2O_OH_He_HP\tH3P_NO2_NOP_OH_H2\tCH_O2_CO_OH\tH2OP_CO_HCOP_OH\tOH_COP_CO_OHP\tOH_HCOP_CO_H2OP\tCH2_OH_H2O_CH\t" 00540 "C_OH_O_CH\tO_CH_OH_C\tO_CH2_OH_CH\tO_H2O_OH_OH\tO_OH_O2_H\tSi_OH_SiO_H\tOH_OH_H2O_O\tO_CH4_OH_CH3\tOH_CH2_O_CH3\tOH_CH3_CH4_O\t" 00541 "OH_CH3_H2O_CH2\tH2O_CH3_OH_CH4\tOH_CH4_H2O_CH3\tN_OH_O_NH\tN_OH_NO_H\tCH2_NO_HCN_OH\tNH_OH_NH2_O\tNH_OH_H2O_N\tNH_H2O_OH_NH2\tNH_NO_N2_OH\t" 00542 "NH_NO2_N2O_OH\tO_NH2_OH_NH\tO_NH3_OH_NH2\tO_HCN_CN_OH\tO_HS_S_OH\tNH2_OH_NH3_O\tOH_NH3_H2O_NH2\tOH_CN_HCN_O\tOH_HCN_CN_H2O\tOH_NO_NO2_H\t" 00543 "OH_N2O_HNO_NO\tOH_CS_OCS_H\tNO_HNO_N2O_OH\tO_C2H2_C2H_OH\tOH_C2H2_C2H_H2O\te_H2OP_OH_H\te_H3OP_OH_H_H\te_H3OP_OH_H2\te_SiOHP_Si_OH\tH_OH_O_H_H\t" 00544 "H_H2O_OH_H_H\tH_OH_O_H2\tH_H2O_OH_H2\tH_OH_H2O_nu\tHminus_H3OP_OH_H2_H\tH_O_OH_nu\tH2_OH_H2O_H\tH2_O_OH_H\tH2_OH_O_H2_H\tH2_H2O_OH_H2_H\t" 00545 "H2_O2_OH_OH\tH_O2_OH_O\tC_OH_CO_H\tOH_HCN_CO_NH2\tOH_C2H2_CO_CH3\tOH_C2H2_CO_CH3\n"); 00546 } 00547 00548 else 00549 { 00550 fprintf(ioQQQ," The keyword CO or OH must appear on the PUNCH CHEMISTRY RATES command.\n"); 00551 fprintf( ioQQQ, " Sorry.\n" ); 00552 cdEXIT(EXIT_FAILURE); 00553 } 00554 00555 } 00556 } 00557 00558 else if( nMatch("COMP",chCard) ) 00559 { 00560 /* punch Compton, details of the energy exchange problem */ 00561 strcpy( punch.chPunch[punch.npunch], "COMP" ); 00562 sprintf( chHeader, 00563 "#nu, comup, comdn\n" ); 00564 } 00565 00566 else if( nMatch("COOL",chCard) ) 00567 { 00568 /* punch cooling, actually done by routine cool_punch */ 00569 strcpy( punch.chPunch[punch.npunch], "COOL" ); 00570 /*>>chng 06 jun 06, revise to be same as punch cooling */ 00571 sprintf( chHeader, 00572 "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\tcool fracs\n" ); 00573 } 00574 00575 else if( nMatch("DYNA",chCard) ) 00576 { 00577 /* punch something dealing with dynamics 00578 * in PunchDo the DYN part of key is used to call DunaPunch, 00579 * with the 4th char as its second argument. DynaPunch uses that 00580 * 4th letter to decide the job */ 00581 if( nMatch( "ADVE" , chCard ) ) 00582 { 00583 /* punch information relating to advection */ 00584 strcpy( punch.chPunch[punch.npunch], "DYNa"); 00585 sprintf( chHeader, 00586 "#advection depth\tHtot\tadCool\tadHeat\tdCoolHeatdT\t" 00587 "Source[hyd][hyd]\tRate\tEnthalph\tadSpecEnthal\n" ); 00588 } 00589 else 00590 { 00591 fprintf( ioQQQ, " I did not recognize a sub keyword on this PUNCH DYNAMICS command.\n" ); 00592 fprintf( ioQQQ, " Sorry.\n" ); 00593 cdEXIT(EXIT_FAILURE); 00594 } 00595 } 00596 00597 else if( nMatch("ENTH",chCard) ) 00598 { 00599 /* contributors to the total enthalpy */ 00600 strcpy( punch.chPunch[punch.npunch], "ENTH" ); 00601 sprintf( chHeader, 00602 "#depth\tTotal\tExcit\tIoniz\tBind\tKE\tther+PdV\tmag \n" ); 00603 } 00604 00605 else if( nMatch("EXEC",chCard) && nMatch("TIME",chCard) ) 00606 { 00607 /* output the execution time per zone */ 00608 strcpy( punch.chPunch[punch.npunch], "XTIM" ); 00609 sprintf( chHeader, 00610 "#zone\tdTime\tElapsed t\n" ); 00611 } 00612 00613 else if( nMatch("FEII",chCard) || nMatch("FE II",chCard) ) 00614 { 00615 /* something to do with FeII atom - several options 00616 * FeII column densities */ 00617 if( nMatch("COLU",chCard) && nMatch("DENS",chCard) ) 00618 { 00619 /* punch FeII column density */ 00620 strcpy( punch.chPunch[punch.npunch], "FEcl" ); 00621 00622 /* file will give energy, statistical weight, and column density [cm-2] */ 00623 sprintf( chHeader, 00624 "#FeII: energy\tstat wght\tcol den\n" ); 00625 } 00626 00627 /* FeII continuum, only valid if large atom is set */ 00628 else if( nMatch("CONT",chCard) ) 00629 { 00630 /* punch FeII continuum */ 00631 strcpy( punch.chPunch[punch.npunch], "FEco" ); 00632 00633 sprintf( chHeader, 00634 "#FeII: Wl(A)\tInt[erg cm-2 s-1]\n" ); 00635 } 00636 00637 else if( nMatch("DEPA",chCard) ) 00638 { 00639 /* punch out departure coefficient for the large FeII atom */ 00640 sprintf( chHeader, 00641 "#FeII departure coefficient \n" ); 00642 /* optional keyword all means do all levels, if not then just do subset */ 00643 if( nMatch(" ALL",chCard) ) 00644 { 00645 /* punch all levels, calls routine FeIIPunDepart */ 00646 strcpy( punch.chPunch[punch.npunch], "FE2D" ); 00647 } 00648 else 00649 { 00650 /* punch a very few selected levels, calls routine FeIIPunDepart */ 00651 strcpy( punch.chPunch[punch.npunch], "FE2d" ); 00652 } 00653 } 00654 00655 else if( nMatch("LEVE",chCard) ) 00656 { 00657 /* punch level energies and statistical weights for large FeII atom */ 00658 sprintf( chHeader, 00659 "#FeII energy(wn)\tstat weight\n" ); 00660 strcpy( punch.chPunch[punch.npunch], "FE2l" ); 00661 } 00662 00663 else if( nMatch("LINE",chCard) ) 00664 { 00665 /* punch FeII lines command 00666 * three optional parameters, threshold for faintest 00667 * line to print, lower and upper energy bounds */ 00668 00669 /* punch intensities from large FeII atom */ 00670 strcpy( punch.chPunch[punch.npunch], "FEli" ); 00671 00672 /* short keyword makes punch half as big */ 00673 if( nMatch("SHOR",chCard) ) 00674 { 00675 FeII.lgShortFe2 = true; 00676 } 00677 else 00678 { 00679 FeII.lgShortFe2 = false; 00680 } 00681 00682 /* first optional number changes the threshold of weakest line to print*/ 00683 i = 5; 00684 /* fe2thresh is intensity relative to normalization line, 00685 * normally Hbeta, and is set to zero in zero.c */ 00686 FeII.fe2thresh = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00687 if( lgEOL ) 00688 { 00689 FeII.fe2thresh = 0.; 00690 } 00691 00692 /* it is a log if negative */ 00693 if( FeII.fe2thresh < 0. ) 00694 { 00695 FeII.fe2thresh = (realnum)pow((realnum)10.f,FeII.fe2thresh); 00696 } 00697 00698 /* check for energy range (Rydberg) of lines to be punched, 00699 * this is to limit size of output file */ 00700 FeII.fe2ener[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00701 if( lgEOL ) 00702 { 00703 FeII.fe2ener[0] = 0.; 00704 } 00705 00706 FeII.fe2ener[1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00707 if( lgEOL ) 00708 { 00709 FeII.fe2ener[1] = 1e8; 00710 } 00711 /* if either is negative then both are logs */ 00712 if( FeII.fe2ener[0] < 0. || FeII.fe2ener[1] < 0. ) 00713 { 00714 FeII.fe2ener[0] = (realnum)pow((realnum)10.f,FeII.fe2ener[0]); 00715 FeII.fe2ener[1] = (realnum)pow((realnum)10.f,FeII.fe2ener[1]); 00716 } 00717 00718 /* entered in Ryd in above, convert to wavenumbers */ 00719 FeII.fe2ener[0] /= (realnum)WAVNRYD; 00720 FeII.fe2ener[1] /= (realnum)WAVNRYD; 00721 00722 /* these results are actually created by the FeIIPunchLines routine 00723 * that lives in the FeIILevelPops file */ 00724 sprintf( chHeader, 00725 "#FeII ipHi\tipLo\tWL(A)\tlog I\tI/Inorm\t\tTau\n" ); 00726 } 00727 00728 else if( nMatch("OPTI",chCard) && nMatch("DEPT",chCard) ) 00729 { 00730 /* punch optical depths for large FeII atom */ 00731 sprintf( chHeader, 00732 "#FeII hi\tlow\twl(A)\ttau\n" ); 00733 strcpy( punch.chPunch[punch.npunch], "FE2o" ); 00734 } 00735 00736 else if( nMatch("POPU",chCard) ) 00737 { 00738 /* punch out level populations for the large FeII atom */ 00739 sprintf( chHeader, 00740 "#FeII level populations [cm^-3]\n" ); 00741 00742 /* this is keyword RELATIVE that says to punch relative to total Fe+, 00743 * default is actual density (cm-3) */ 00744 if( nMatch("RELA",chCard) ) 00745 { 00746 punch.punarg[punch.npunch][2] = 0.; 00747 } 00748 else 00749 { 00750 /* default is to punch density (cm-3) */ 00751 punch.punarg[punch.npunch][2] = 1.; 00752 } 00753 00754 /* optional keyword all means do all levels, if not then just do subset */ 00755 if( nMatch(" ALL",chCard) ) 00756 { 00757 /* punch all levels, calls routine FeIIPunPop */ 00758 strcpy( punch.chPunch[punch.npunch], "FE2P" ); 00759 punch.punarg[punch.npunch][0] = 0.; 00760 punch.punarg[punch.npunch][1] = (realnum)FeII.nFeIILevel; 00761 } 00762 00763 /* optional keyword range means read lower and upper bound, do these */ 00764 else if( nMatch("RANG",chCard) ) 00765 { 00766 /* punch range of levels, calls routine FeIIPunPop */ 00767 strcpy( punch.chPunch[punch.npunch], "FE2P" ); 00768 ipFFmt = 5; 00769 punch.punarg[punch.npunch][0] = 00770 (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL); 00771 punch.punarg[punch.npunch][1] = 00772 (realnum)FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL); 00773 if( lgEOL || punch.punarg[punch.npunch][0] <0 || 00774 punch.punarg[punch.npunch][0]> punch.punarg[punch.npunch][1] || 00775 punch.punarg[punch.npunch][1]> FeII.nFeIILevel) 00776 { 00777 fprintf( ioQQQ, "There must be two numbers on this punch FeII populations range command.\n" ); 00778 fprintf( ioQQQ, "These give the lower and upper levels for the range of FeII levels.\n" ); 00779 fprintf( ioQQQ, "The first must be less than the second and the second must be <= the number of FeII levels.\n" ); 00780 fprintf( ioQQQ, "Sorry.\n" ); 00781 cdEXIT(EXIT_FAILURE); 00782 } 00783 } 00784 00785 else 00786 { 00787 /* punch a very few selected levels, calls routine FeIIPunPop */ 00788 strcpy( punch.chPunch[punch.npunch], "FE2p" ); 00789 } 00790 } 00791 else 00792 { 00793 fprintf( ioQQQ, "There must be a second keyword on this PUNCH FEII command.\n" ); 00794 fprintf( ioQQQ, "The ones I know about are COLUmn, CONTinuum, " 00795 "DEPArture, LEVEls, LINE, OPTIcal DEPTh, and POPUlations.\n" ); 00796 fprintf( ioQQQ, "Sorry.\n" ); 00797 cdEXIT(EXIT_FAILURE); 00798 } 00799 } 00800 00801 /* the punch continuum command, with many options, 00802 * the first 3 char of the chPunch flag will always be "CON" 00803 * with the last indicating which one */ 00804 else if( nMatch("CONT",chCard) && !nMatch("XSPE",chCard) ) 00805 { 00806 /* this flag is checked in PrtComment to generate a caution 00807 * if continuum is punched but iterations not performed */ 00808 punch.lgPunContinuum = true; 00809 00810 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 00811 * units are copied into punch.chConPunEnr */ 00812 ChkUnits(chCard); 00813 00814 if( nMatch("BINS",chCard) ) 00815 { 00816 /* continuum binning */ 00817 strcpy( punch.chPunch[punch.npunch], "CONB" ); 00818 00819 sprintf( chHeader, 00820 "#Continuum binning, enrOrg, Energy, width of cells\n" ); 00821 } 00822 00823 else if( nMatch("DIFF",chCard) ) 00824 { 00825 /* diffuse continuum, the locally emitted continuum stored in rfield.ConEmitLocal */ 00826 strcpy( punch.chPunch[punch.npunch], "COND" ); 00827 00828 sprintf( chHeader, 00829 "#energy\t ConEmitLocal \n" ); 00830 } 00831 00832 else if( nMatch("EMIT",chCard) ) 00833 { 00834 /* continuum emitted by cloud */ 00835 strcpy( punch.chPunch[punch.npunch], "CONE" ); 00836 00837 sprintf( chHeader, 00838 "#Energy\treflec\toutward\ttotal\tline\tcont\n" ); 00839 } 00840 00841 else if( nMatch("FINE" , chCard ) ) 00842 { 00843 rfield.lgPunchOpacityFine = true; 00844 /* fine transmitted continuum cloud */ 00845 strcpy( punch.chPunch[punch.npunch], "CONf" ); 00846 00847 sprintf( chHeader, 00848 "#Energy\tTransmitted\n" ); 00849 00850 ipFFmt = 5; 00851 /* range option - important since so much data */ 00852 if( nMatch("RANGE",chCard) ) 00853 { 00854 /* get lower and upper range, must be in Ryd */ 00855 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 00856 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 00857 if( lgEOL ) 00858 { 00859 fprintf(ioQQQ,"There must be two numbers, the lower and upper energies in Ryd.\nSorry.\n"); 00860 cdEXIT(EXIT_FAILURE); 00861 } 00862 if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] ) 00863 { 00864 fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n"); 00865 cdEXIT(EXIT_FAILURE); 00866 } 00867 } 00868 else 00869 { 00870 /* these mean full energy range */ 00871 punch.punarg[punch.npunch][0] = 0.; 00872 punch.punarg[punch.npunch][1] = 0.; 00873 } 00874 /* optional last parameter - how many points to bring together */ 00875 punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt, 00876 INPUT_LINE_LENGTH,&lgEOL); 00877 00878 /* default is to bring together ten */ 00879 if( lgEOL ) 00880 punch.punarg[punch.npunch][2] = 10; 00881 00882 if( punch.punarg[punch.npunch][2] < 1 ) 00883 { 00884 fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n"); 00885 cdEXIT(EXIT_FAILURE); 00886 } 00887 } 00888 00889 else if( nMatch("GRAI",chCard) ) 00890 { 00891 /* punch grain continuum in optically thin limit */ 00892 strcpy( punch.chPunch[punch.npunch], "CONG" ); 00893 00894 sprintf( chHeader, 00895 "#energy\tgraphite\trest\ttotal\n" ); 00896 } 00897 00898 else if( nMatch("INCI",chCard) ) 00899 { 00900 /* incident continuum */ 00901 strcpy( punch.chPunch[punch.npunch], "CONC" ); 00902 00903 sprintf( chHeader, 00904 "#Incident Continuum, Enr\tnFn \n" ); 00905 } 00906 00907 else if( nMatch("INTE",chCard) ) 00908 { 00909 /* continuum interactions */ 00910 strcpy( punch.chPunch[punch.npunch], "CONi" ); 00911 00912 sprintf( chHeader, 00913 "#Continuum interactions, inc, otslin. otscon, ConInterOut, outlin \n" ); 00914 /* this is option for lowest energy, if nothing then zero */ 00915 ipFFmt = 3; 00916 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt, 00917 INPUT_LINE_LENGTH,&lgEOL); 00918 } 00919 00920 else if( nMatch("IONI",chCard) ) 00921 { 00922 /* punch ionizing continuum*/ 00923 strcpy( punch.chPunch[punch.npunch], "CONI" ); 00924 00925 /* this is option for lowest energy, if nothing then zero */ 00926 ipFFmt = 3; 00927 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt, 00928 INPUT_LINE_LENGTH,&lgEOL); 00929 00930 /* this is option for smallest interaction to punch, def is 1 percent */ 00931 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 00932 if( lgEOL ) 00933 { 00934 punch.punarg[punch.npunch][1] = 0.01f; 00935 } 00936 00937 /* >>chng 03 may 03, add "every" option to punch this on every zone - 00938 * if every is not present then only last zone is punched */ 00939 if( nMatch("EVER", chCard ) ) 00940 { 00941 /* punch every zone */ 00942 punch.lgPunchEveryZone[punch.npunch] = true; 00943 punch.nPunchEveryZone[punch.npunch] = 1; 00944 } 00945 else 00946 { 00947 /* only punch last zone */ 00948 punch.lgPunchEveryZone[punch.npunch] = false; 00949 punch.nPunchEveryZone[punch.npunch] = 1; 00950 } 00951 00952 /* put the header at the top of the file */ 00953 sprintf( chHeader, 00954 "#cell(on C scale)\tnu\tflux\tflx*cs\tFinc\totsl\totsc\toutlin\toutcon\trate/tot\tintegral\tline\tcont\n" ); 00955 } 00956 00957 else if( nMatch("OUTW",chCard) ) 00958 { 00959 /* outward only continuum */ 00960 if( nMatch("LOCA",chCard) ) 00961 { 00962 strcpy( punch.chPunch[punch.npunch], "CONo" ); 00963 sprintf( chHeader, 00964 "#Local Out ConInterOut+line SvOt*opc pass*opc\n" ); 00965 } 00966 else 00967 { 00968 strcpy( punch.chPunch[punch.npunch], "CONO" ); 00969 sprintf( chHeader, 00970 "#Out Con OutIncid OutConD OutLinD OutConS\n" ); 00971 } 00972 } 00973 00974 else if( nMatch("TRAN",chCard) ) 00975 { 00976 /* transmitted continuum */ 00977 strcpy( punch.chPunch[punch.npunch], "CONT" ); 00978 00979 sprintf( chHeader, 00980 "#ener\tTran Contin\ttrn coef\n" ); 00981 } 00982 00983 else if( nMatch(" TWO",chCard) ) 00984 { 00985 /* total two photon continua rfield.TotDiff2Pht */ 00986 strcpy( punch.chPunch[punch.npunch], "CON2" ); 00987 00988 sprintf( chHeader, 00989 "#energy\t n_nu\tnuF_nu \n" ); 00990 } 00991 00992 else if( nMatch(" RAW",chCard) ) 00993 { 00994 /* "raw" continua */ 00995 strcpy( punch.chPunch[punch.npunch], "CORA" ); 00996 00997 sprintf( chHeader, 00998 "#Raw Con anu\tflux\totslin\totscon\tConRefIncid\tConEmitReflec\tConInterOut\toutlin\tConEmitOut\tline\tcont\tnLines\n" ); 00999 } 01000 01001 else if( nMatch("REFL",chCard) ) 01002 { 01003 /* reflected continuum */ 01004 strcpy( punch.chPunch[punch.npunch], "CONR" ); 01005 01006 sprintf( chHeader, 01007 "#Reflected\tcont\tline\ttotal\talbedo\tConID\n" ); 01008 } 01009 01010 else 01011 { 01012 /* this is the usual punch continuum command */ 01013 strcpy( punch.chPunch[punch.npunch], "CON " ); 01014 sprintf( chHeader, 01015 "#Cont nu\tincident\ttrans\tDiffOut\tnet trans\treflc\ttotal\tline\tcont\tnLine\n" ); 01016 01017 /* >>chng 06 apr 03, add "every" option to punch this on every zone - 01018 * if every is not present then only last zone is punched */ 01019 if( nMatch("EVER", chCard ) ) 01020 { 01021 /* punch every zone */ 01022 punch.lgPunchEveryZone[punch.npunch] = true; 01023 /* option to say how many to skip */ 01024 ipFFmt = 5; 01025 punch.nPunchEveryZone[punch.npunch] = (long)FFmtRead(chCard,&ipFFmt, 01026 INPUT_LINE_LENGTH,&lgEOL); 01027 if( lgEOL ) 01028 punch.nPunchEveryZone[punch.npunch] = 1; 01029 } 01030 else 01031 { 01032 /* only punch last zone */ 01033 punch.lgPunchEveryZone[punch.npunch] = false; 01034 punch.nPunchEveryZone[punch.npunch] = 1; 01035 } 01036 } 01037 } 01038 01039 /* punch information about convergence of this model 01040 * reason - why it did not converge an iteration 01041 * error - zone by zone display of various convergence errors */ 01042 else if( nMatch("CONV",chCard) ) 01043 { 01044 if( nMatch("REAS",chCard) ) 01045 { 01046 /* this does not count as a punch option (really) */ 01047 punch.lgPunConv = true; 01048 /* this is done below */ 01049 strcpy( punch.chPunch[punch.npunch], "" ); 01050 punch.lgRealPunch[punch.npunch] = false; 01051 } 01052 else if( nMatch("ERRO",chCard) ) 01053 { 01054 /* punch zone by zone errors in pressure, electron density, and heating-cooling */ 01055 /* convergence error */ 01056 strcpy( punch.chPunch[punch.npunch], "CNVE" ); 01057 sprintf( chHeader, 01058 "#depth\tnPres2Ioniz\tP(cor)\tP(cur)\tP%%error\tNE(cor)\tNE(cur)\tNE%%error\tHeat\tCool\tHC%%error\n" ); 01059 } 01060 else if( nMatch("BASE",chCard) ) 01061 { 01062 /* punch converged quantities in Converge base for each pass through 01063 * solvers - not one pass per zone */ 01064 strcpy( punch.chPunch[punch.npunch], "CNVB" ); 01065 strcpy( punch.chPunch[punch.npunch], "" ); 01066 punch.lgRealPunch[punch.npunch] = false; 01067 } 01068 else 01069 { 01070 fprintf( ioQQQ, "There must be a second keyword on this command.\n" ); 01071 fprintf( ioQQQ, "The ones I know about are REASON and ERROR.\n" ); 01072 fprintf( ioQQQ, "Sorry.\n" ); 01073 cdEXIT(EXIT_FAILURE); 01074 } 01075 } 01076 01077 else if( nMatch(" DR ",chCard) ) 01078 { 01079 /* first occurrence of punch dr to follow choice in change of zoning */ 01080 punch.lgDROn = true; 01081 strcpy( punch.chPunch[punch.npunch], "" ); 01082 punch.lgRealPunch[punch.npunch] = false; 01083 } 01084 01085 else if( nMatch("ELEM",chCard) && !nMatch("GAMMA" , chCard) ) 01086 { 01087 /* option to punch ionization structure of some element 01088 * will give each stage of ionization, vs depth */ 01089 strcpy( punch.chPunch[punch.npunch], "ELEM" ); 01090 01091 /* this returns element number on c scale */ 01092 /* >>chng 04 nov 23, had converted to f scale, leave on c */ 01093 nelem = GetElem(chCard); 01094 01095 /* this is the atomic number on the c scale */ 01096 punch.punarg[punch.npunch][0] = (realnum)nelem; 01097 ASSERT( nelem>=ipHYDROGEN); 01098 01099 /* >>chng 04 nov 24, add DENSE option to print density rather than fraction */ 01100 punch.punarg[punch.npunch][1] = 0; 01101 if( nMatch("DENS",chCard) ) 01102 punch.punarg[punch.npunch][1] = 1.; 01103 01104 /* start printing header line - first will be the depth in cm */ 01105 sprintf( chHeader, "#depth"); 01106 01107 /* next come the nelem+1 ion stages */ 01108 for(i=0; i<=nelem+1;++i ) 01109 { 01110 sprintf( chTemp, 01111 "\t%2s%2li", elementnames.chElementSym[nelem],i+1); 01112 strcat( chHeader, chTemp ); 01113 } 01114 01115 /* finally some fine structure or molecular populations */ 01116 /* >>chng 04 nov 23, add fs pops of C, O 01117 * >>chng 04 nov 25, add molecules */ 01118 if( nelem==ipHYDROGEN ) 01119 { 01120 sprintf( chTemp, "\tH2"); 01121 strcat( chHeader, chTemp ); 01122 } 01123 if( nelem==ipCARBON ) 01124 { 01125 sprintf( chTemp, "\tC1\tC1*\tC1**\tC2\tC2*\tCO"); 01126 strcat( chHeader, chTemp ); 01127 } 01128 else if( nelem==ipOXYGEN ) 01129 { 01130 sprintf( chTemp, "\tO1\tO1*\tO1**"); 01131 strcat( chHeader, chTemp ); 01132 } 01133 01134 /* finally the new line */ 01135 strcat( chHeader, "\n"); 01136 } 01137 01138 else if( nMatch("FITS",chCard) ) 01139 { 01140 01141 #ifdef FLT_IS_DBL 01142 fprintf( ioQQQ, "Punching FITS files is not currently supported in double precision.\n" ); 01143 fprintf( ioQQQ, "Please recompile without the FLT_IS_DBL option.\n" ); 01144 fprintf( ioQQQ, "Sorry.\n" ); 01145 cdEXIT(EXIT_FAILURE); 01146 #else 01147 /* say that this is a FITS file output */ 01148 punch.lgFITS[punch.npunch] = true; 01149 01150 strcpy( punch.chPunch[punch.npunch], "FITS" ); 01151 #endif 01152 01153 } 01154 01155 else if( nMatch("FRED",chCard) ) 01156 { 01157 /* punch out some stuff for Fred's dynamics project */ 01158 sprintf( chHeader, 01159 "#Radius\tDepth\tVelocity\thden\teden\tTemperature\tRadAccel line\tRadAccel con\t" 01160 "Force multiplier\t" 01161 "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t" 01162 "O2\tO3\tO4\tO5\tO6\tO7\tO8\t" 01163 "HI\tHII\tHeI\tHeII\tHeIII\tC2\tC3\tC4\tO1\t" 01164 "O2\tO3\tO4\tO5\tO6\tO7\tO8\tMg2\tMg2\n"); 01165 strcpy( punch.chPunch[punch.npunch], "FRED" ); 01166 } 01167 01168 else if( nMatch("GAMM",chCard) ) 01169 { 01170 /* punch all photoionization rates for all subshells */ 01171 sprintf( chHeader, 01172 "#Photoionization rates \n" ); 01173 if( nMatch("ELEMENT" , chCard ) ) 01174 { 01175 /* element keyword, find element name and stage of ionization, 01176 * will print photoionization rates for valence of that element */ 01177 strcpy( punch.chPunch[punch.npunch], "GAMe" ); 01178 01179 /* this returns element number on c scale */ 01180 nelem = GetElem(chCard); 01181 /* this is the atomic number on the C scale */ 01182 punch.punarg[punch.npunch][0] = (realnum)nelem; 01183 01184 /* this will become the ionization stage on C scale */ 01185 ipFFmt = 5; 01186 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt, 01187 INPUT_LINE_LENGTH,&lgEOL) - 1; 01188 if( lgEOL ) 01189 NoNumb( chCard ); 01190 if( punch.punarg[punch.npunch][1]<0 || punch.punarg[punch.npunch][1]> nelem+1 ) 01191 { 01192 fprintf(ioQQQ,"Bad ionization stage - please check Hazy.\nSorry.\n"); 01193 cdEXIT(EXIT_FAILURE); 01194 } 01195 } 01196 else 01197 { 01198 /* no element - so make table of all rates */ 01199 strcpy( punch.chPunch[punch.npunch], "GAMt" ); 01200 } 01201 01202 } 01203 else if( nMatch("GRAI",chCard) ) 01204 { 01205 /* punch grain ... options */ 01206 if( nMatch("OPAC",chCard) ) 01207 { 01208 /* check for keyword UNITS on line, then scan wavelength or energy units, 01209 * sets punch.chConPunEnr*/ 01210 ChkUnits(chCard); 01211 01212 strcpy( punch.chPunch[punch.npunch], "DUSO" ); 01213 /* punch grain opacity command in twice, here and above in opacity */ 01214 sprintf( chHeader, 01215 "#grain\tnu\tabs+scat*(1-g)\tabs\tscat*(1-g)\tscat\tscat*(1-g)/[abs+scat*(1-g)]\n" ); 01216 } 01217 else if( nMatch("ABUN",chCard) ) 01218 { 01219 /* punch grain abundance */ 01220 strcpy( punch.chPunch[punch.npunch], "DUSA" ); 01221 sprintf( chHeader, 01222 "#grain\tdepth\tabundance (g/cm^3)\n" ); 01223 } 01224 else if( nMatch("D/G ",chCard) ) 01225 { 01226 /* punch grain dust/gas mass ratio */ 01227 strcpy( punch.chPunch[punch.npunch], "DUSD" ); 01228 sprintf( chHeader, 01229 "#grain\tdepth\tdust/gas mass ratio\n" ); 01230 } 01231 else if( nMatch("PHYS",chCard) ) 01232 { 01233 /* punch grain physical conditions */ 01234 strcpy( punch.chPunch[punch.npunch], "DUSP" ); 01235 sprintf( chHeader, 01236 "#grain\tdepth\tpotential\n" ); 01237 } 01238 else if( nMatch(" QS ",chCard) ) 01239 { 01240 strcpy( punch.chPunch[punch.npunch], "DUSQ" ); 01241 sprintf( chHeader, 01242 "#grain\tnu\tQ_abs\tQ_scat\n" ); 01243 } 01244 else if( nMatch("TEMP",chCard) ) 01245 { 01246 /* punch temperatures of each grain species */ 01247 strcpy( punch.chPunch[punch.npunch], "DUST" ); 01248 /* cannot punch grain labels since they are not known yet */ 01249 sprintf( chHeader, 01250 "#grain temperature\n" ); 01251 } 01252 else if( nMatch("DRIF",chCard) ) 01253 { 01254 /* punch drift velocity of each grain species */ 01255 strcpy( punch.chPunch[punch.npunch], "DUSV" ); 01256 /* cannot punch grain labels since they are not known yet */ 01257 sprintf( chHeader, 01258 "#grain drift velocity\n" ); 01259 } 01260 else if( nMatch("EXTI",chCard) ) 01261 { 01262 /* punch grain extinction */ 01263 strcpy( punch.chPunch[punch.npunch], "DUSE" ); 01264 /* cannot punch grain labels since they are not known yet */ 01265 sprintf( chHeader, 01266 "#depth\tA_V(extended)\tA_V(point)\n" ); 01267 } 01268 else if( nMatch("CHAR",chCard) ) 01269 { 01270 /* punch charge per grain (# elec/grain) for each grain species */ 01271 strcpy( punch.chPunch[punch.npunch], "DUSC" ); 01272 /* cannot punch grain labels since they are not known yet */ 01273 sprintf( chHeader, 01274 "#grain charge\n" ); 01275 } 01276 else if( nMatch("HEAT",chCard) ) 01277 { 01278 /* punch heating due to each grain species */ 01279 strcpy( punch.chPunch[punch.npunch], "DUSH" ); 01280 /* cannot punch grain labels since they are not known yet */ 01281 sprintf( chHeader, 01282 "#grain heating\n" ); 01283 } 01284 else if( nMatch("POTE",chCard) ) 01285 { 01286 /* punch floating potential of each grain species */ 01287 strcpy( punch.chPunch[punch.npunch], "DUSP" ); 01288 /* cannot punch grain labels since they are not known yet */ 01289 sprintf( chHeader, 01290 "#grain\tdepth\tpotential\n" ); 01291 } 01292 else if( nMatch("H2RA",chCard) ) 01293 { 01294 /* punch grain H2rate - H2 formation rate for each type of grains */ 01295 strcpy( punch.chPunch[punch.npunch], "DUSR" ); 01296 /* cannot punch grain labels since they are not known yet */ 01297 sprintf( chHeader, 01298 "#grain H2 formation rates\n" ); 01299 } 01300 else 01301 { 01302 fprintf( ioQQQ, "There must be a second key on this GRAIN command; The options I know about follow (required key in CAPS):\n"); 01303 fprintf( ioQQQ, "OPACity, ABUNdance, D/G mass ratio, PHYSical conditions, QS , TEMPerature, DRIFt velocity, EXTInction, CHARge, HEATing, POTEntial, H2RAtes\nSorry.\n" ); 01304 cdEXIT(EXIT_FAILURE); 01305 } 01306 } 01307 01308 else if( nMatch("GAUN",chCard) ) 01309 { 01310 strcpy( punch.chPunch[punch.npunch], "GAUN" ); 01311 sprintf( chHeader, 01312 "#Gaunt factors.\n" ); 01313 } 01314 else if( nMatch("GRID",chCard) ) 01315 { 01316 strcpy( punch.chPunch[punch.npunch], "GRID" ); 01317 /* automatically generate no hash option */ 01318 punch.lgHashEndIter[punch.npunch] = false; 01319 } 01320 else if( nMatch( "HIST" , chCard ) ) 01321 { 01322 /* punch pressure history of current zone */ 01323 if( nMatch( "PRES" , chCard ) ) 01324 { 01325 /* punch pressure history - density - pressure for this zone */ 01326 strcpy( punch.chPunch[punch.npunch], "HISp" ); 01327 sprintf( chHeader, 01328 "#iter zon\tdensity\tpres cur\tpres corr\n" ); 01329 } 01330 /* punch temperature history of current zone */ 01331 else if( nMatch( "TEMP" , chCard ) ) 01332 { 01333 /* punch pressure history - density - pressure for this zone */ 01334 strcpy( punch.chPunch[punch.npunch], "HISt" ); 01335 sprintf( chHeader, 01336 "#iter zon\ttemperature\theating\tcooling\n" ); 01337 } 01338 } 01339 01340 else if( nMatch("HTWO",chCard) ) 01341 { 01342 fprintf(ioQQQ," Sorry, this command has been replaced with the " 01343 "PUNCH H2 CREATION and PUNCH H2 DESTRUCTION commands.\n"); 01344 cdEXIT(EXIT_FAILURE); 01345 } 01346 01347 /* QHEAT has to come before HEAT... */ 01348 else if( nMatch("QHEA",chCard) ) 01349 { 01350 /* this is just a dummy clause, do the work below after parsing is over. 01351 * this is a no-nothing, picked up to stop optimizer */ 01352 ((void)0); 01353 } 01354 01355 else if( nMatch("HEAT",chCard) ) 01356 { 01357 /* punch heating */ 01358 strcpy( punch.chPunch[punch.npunch], "HEAT" ); 01359 /*>>chng 06 jun 06, revise to be same as punch cooling */ 01360 sprintf( chHeader, 01361 "#depth cm\tTemp K\tHtot erg/cm3/s\tCtot erg/cm3/s\theat fracs\n" ); 01362 } 01363 01364 else if( nMatch("HELI",chCard) &&!( nMatch("IONI",chCard))) 01365 { 01366 /* punch helium & helium-like iso sequence, but not punch helium ionization rate 01367 * punch helium line wavelengths */ 01368 if( nMatch("LINE",chCard) && nMatch("WAVE",chCard) ) 01369 { 01370 strcpy( punch.chPunch[punch.npunch], "HELW" ); 01371 sprintf( chHeader, 01372 "#wavelengths of lines from He-like ions\n" ); 01373 } 01374 else 01375 { 01376 fprintf( ioQQQ, "punch helium has options: LINE WAVElength.\nSorry.\n" ); 01377 cdEXIT(EXIT_FAILURE); 01378 /* no key */ 01379 } 01380 } 01381 01382 else if( nMatch("HUMM",chCard) ) 01383 { 01384 strcpy( punch.chPunch[punch.npunch], "HUMM" ); 01385 sprintf( chHeader, 01386 "#input to DHs routine.\n" ); 01387 } 01388 01389 else if( nMatch("HYDR",chCard) ) 01390 { 01391 /* punch hydrogen physical conditions */ 01392 if( nMatch("COND",chCard) ) 01393 { 01394 strcpy( punch.chPunch[punch.npunch], "HYDc" ); 01395 sprintf( chHeader, 01396 "#depth\tTe\tHDEN\tEDEN\tHI/H\tHII/H\tH2/H\tH2+/H\tH3+/H\tH-/H\n" ); 01397 /* punch hydrogen ionization */ 01398 } 01399 01400 /* punch information on 21 cm excitation processes - accept either keyword 21cm or 21 cm */ 01401 else if( nMatch("21 CM",chCard) ||nMatch("21CM",chCard)) 01402 { 01403 /* punch information about 21 cm line */ 01404 strcpy( punch.chPunch[punch.npunch], "21CM" ); 01405 sprintf( chHeader, 01406 "#depth\tT(spin)\tT(kin)\tT(Lya/21cm)\tnLo\tnHi\tOccLya\ttau(21cm)" 01407 "\ttau(Lya)\topac(21 cm)\tn/Ts\ttau(21)\tTex(Lya)\tN(H0)/Tspin" 01408 "\tSum_F0\tSum_F1\tSum_T21\n" ); 01409 } 01410 01411 else if( nMatch("IONI",chCard) ) 01412 { 01413 /* punch hydrogen ionization */ 01414 strcpy( punch.chPunch[punch.npunch], "HYDi" ); 01415 sprintf( chHeader, 01416 "#hion\tzn\tgam1\tcoll ion1\tRecTot\tHRecCaB\thii/hi\tSim hii/hi" 01417 "\time_Hrecom_long(esc)\tdec2grd\texc pht\texc col\trec eff\tsec ion\n" ); 01418 } 01419 else if( nMatch("POPU",chCard) ) 01420 { 01421 /* punch hydrogen populations */ 01422 strcpy( punch.chPunch[punch.npunch], "HYDp" ); 01423 sprintf( chHeader, 01424 "#depth\tn(H0)\tn(H+)\tn(1s)\tn(2s)\tn(2p)\tetc\n" ); 01425 } 01426 else if( nMatch("LINE",chCard) ) 01427 { 01428 /* punch hydrogen lines 01429 * hydrogen line intensities and optical depths */ 01430 strcpy( punch.chPunch[punch.npunch], "HYDl" ); 01431 sprintf( chHeader, 01432 "#nup\tnlo\tE(ryd)\ttau\n" ); 01433 } 01434 else if( nMatch(" LYA",chCard) ) 01435 { 01436 /* punch hydrogen Lya some details about Lyman alpha */ 01437 strcpy( punch.chPunch[punch.npunch], "HYDL" ); 01438 sprintf( chHeader, 01439 "#depth\tTauIn\tTauTot\tn(2p)/n(1s)\tTexc\tTe\tTex/T\tPesc\tPdes\tpump\topacity\talbedo\n" ); 01440 } 01441 else 01442 { 01443 fprintf( ioQQQ, "Punch hydrogen has options: CONDitions, 21 CM, LINE, POPUlations, and IONIzation.\nSorry.\n" ); 01444 cdEXIT(EXIT_FAILURE); 01445 } 01446 } 01447 01448 else if( nMatch("IONI",chCard) ) 01449 { 01450 if( nMatch("RATE",chCard) ) 01451 { 01452 /* punch ionization rates, search for the name of an element */ 01453 if( (nelem = GetElem( chCard ) ) < 0 ) 01454 { 01455 fprintf( ioQQQ, "There must be an element name on the ionization rates command. Sorry.\n" ); 01456 cdEXIT(EXIT_FAILURE); 01457 } 01458 punch.punarg[punch.npunch][0] = (realnum)nelem; 01459 strcpy( punch.chPunch[punch.npunch], "IONR" ); 01460 sprintf( chHeader, 01461 "#%s depth\teden\tdynamics.Rate\tabund\tTotIonize\tTotRecom\tSource\t ... \n", 01462 elementnames.chElementSym[nelem]); 01463 } 01464 else 01465 { 01466 /* punch table giving ionization means */ 01467 strcpy( punch.chPunch[punch.npunch], "IONI" ); 01468 sprintf( chHeader, 01469 "#Mean ionization distribution\n" ); 01470 } 01471 } 01472 01473 else if( nMatch(" IP ",chCard) ) 01474 { 01475 strcpy( punch.chPunch[punch.npunch], " IP " ); 01476 sprintf( chHeader, 01477 "#ionization potentials, valence shell\n" ); 01478 } 01479 01480 else if( nMatch("LEID",chCard) ) 01481 { 01482 if( nMatch( "LINE" , chCard ) ) 01483 { 01484 /* punch Leiden lines 01485 * final intensities of the Leiden PDR models */ 01486 strcpy( punch.chPunch[punch.npunch], "LEIL" ); 01487 sprintf( chHeader, "#ion\twl\tInt\trel int\n"); 01488 } 01489 else 01490 { 01491 /* punch Leiden structure 01492 * structure of the Leiden PDR models */ 01493 strcpy( punch.chPunch[punch.npunch], "LEIS" ); 01494 sprintf( chHeader, 01495 /* 1-17 */ 01496 "#Leid depth\tA_V(extentd)\tA_V(point)\tTe\tH0\tH2\tCo\tC+\tOo\tCO\tO2\tCH\tOH\te\tHe+\tH+\tH3+\t" 01497 /* 18 - 30 */ 01498 "N(H0)\tN(H2)\tN(Co)\tN(C+)\tN(Oo)\tN(CO)\tN(O2)\tN(CH)\t(OH)\tN(e)\tN(He+)\tN(H+)\tN(H3+)\t" 01499 /* 31 - 32 */ 01500 "H2(Sol)\tH2(FrmGrn)\t" 01501 /* 33 - 46*/ 01502 "G0(DB96)\trate(CO)\trate(C)\theat\tcool\tGrnP\tGr-Gas-Cool\tGr-Gas-Heat\tCOds\tH2dH\tH2vH\tChaT\tCR H\tMgI\tSI\t" 01503 "Si\tFe\tNa\tAl\tC\tC610\tC370\tC157\tC63\tC146\n" ); 01504 } 01505 } 01506 01507 else if( nMatch("LINE",chCard) && nMatch("LIST",chCard) ) 01508 { 01509 /* punch line list "output file" "Line List file" */ 01510 long int j; 01511 /* 01512 * we parsed off the second file name at start of this routine 01513 * check if file was found, use it if it was, else abort 01514 */ 01515 if( !lgSecondFilename ) 01516 { 01517 fprintf(ioQQQ , "There must be a second file name between " 01518 "double quotes on the PUNCH LINE LIST command. This second" 01519 " file contains the input line list. I did not find it.\nSorry.\n"); 01520 cdEXIT(EXIT_FAILURE); 01521 } 01522 01523 /* actually get the lines, and malloc the space in the arrays 01524 * cdGetLineList will look on path */ 01525 if( punch.ipPnunit[punch.npunch] == NULL && 01526 ( punch.nLineList[punch.npunch] = cdGetLineList( chSecondFilename , 01527 &punch.chLineListLabel[punch.npunch] , 01528 &punch.wlLineList[punch.npunch] ) ) < 0 ) 01529 { 01530 fprintf(ioQQQ,"DISASTER could not open PUNCH LINE LIST file %s \n", 01531 chSecondFilename ); 01532 cdEXIT(EXIT_FAILURE); 01533 } 01534 01535 /* ratio option, in which pairs of lines form ratios, first over 01536 * second */ 01537 if( nMatch("RATI",chCard) ) 01538 { 01539 punch.lgLineListRatio[punch.npunch] = true; 01540 if( punch.nLineList[punch.npunch]%2 ) 01541 { 01542 /* odd number of lines - cannot take ratio */ 01543 fprintf(ioQQQ , "There must be an even number of lines to" 01544 " take ratios of lines. There were %li, an odd number." 01545 "\nSorry.\n", punch.nLineList[punch.npunch]); 01546 cdEXIT(EXIT_FAILURE); 01547 } 01548 } 01549 else 01550 { 01551 /* no ratio */ 01552 punch.lgLineListRatio[punch.npunch] = false; 01553 } 01554 01555 /* do special string saying our job, and then print labels */ 01556 strcpy( punch.chPunch[punch.npunch], "LLST" ); 01557 01558 /* keyword absolute says to do absolute rather than relative intensities 01559 * relative intensities are the default */ 01560 if( nMatch("ABSO",chCard) ) 01561 { 01562 punch.punarg[punch.npunch][0] = 1; 01563 } 01564 else 01565 { 01566 punch.punarg[punch.npunch][0] = 0; 01567 } 01568 01569 /* give header line */ 01570 sprintf( chHeader, "#lineslist" ); 01571 for( j=0; j<punch.nLineList[punch.npunch]; ++j ) 01572 { 01573 /* if taking ratio then put div sign between pairs */ 01574 if( punch.lgLineListRatio[punch.npunch] && is_odd(j) ) 01575 strcat( chHeader , "/" ); 01576 else 01577 strcat( chHeader , "\t" ); 01578 sprintf( chTemp, "%s ", punch.chLineListLabel[punch.npunch][j] ); 01579 strcat( chHeader, chTemp ); 01580 sprt_wl( chTemp, punch.wlLineList[punch.npunch][j] ); 01581 strcat( chHeader, chTemp ); 01582 } 01583 strcat( chHeader, "\n" ); 01584 } 01585 01586 else if( nMatch("LINE",chCard) && !nMatch("XSPE",chCard) && !nMatch("NEAR",chCard)) 01587 { 01588 /* punch line emissivity - 01589 * this is not punch xspec lines and not linear option 01590 * check for keyword UNITS on line, then scan wavelength or energy units, 01591 * sets punch.chConPunEnr*/ 01592 ChkUnits(chCard); 01593 01594 /* punch line emissivity, line intensity, line array, 01595 * and line data */ 01596 if( nMatch("STRU",chCard) ) 01597 { 01598 fprintf(ioQQQ," The PUNCH LINES STRUCTURE command is now PUNCH LINES " 01599 "EMISSIVITY.\n Sorry.\n\n"); 01600 cdEXIT(EXIT_FAILURE); 01601 } 01602 01603 else if( nMatch("EMIS",chCard) ) 01604 { 01605 /* this used to be the punch lines structure command, is now 01606 * the punch lines emissivity command 01607 * give line emissivity vs depth */ 01608 strcpy( punch.chPunch[punch.npunch], "LINS" ); 01609 sprintf( chHeader, 01610 "#Emission line structure:"); 01611 /* read in the list of lines to examine */ 01612 punch_line(punch.ipPnunit[punch.npunch],"READ",false, chTemp); 01613 strcat( chHeader, chTemp ); 01614 } 01615 01616 else if( nMatch(" RT ", chCard ) ) 01617 { 01618 /* punch line RT */ 01619 strcpy( punch.chPunch[punch.npunch], "LINR" ); 01620 /* punch some details needed for line radiative transfer 01621 * routine in punch_line.c */ 01622 Punch_Line_RT(punch.ipPnunit[punch.npunch],"READ"); 01623 } 01624 01625 else if( nMatch("CUMU",chCard) ) 01626 { 01627 /* punch lines cumulative 01628 * this will be integrated line intensity, function of depth */ 01629 strcpy( punch.chPunch[punch.npunch], "LINC" ); 01630 /* option for either relative intensity or abs luminosity */ 01631 if( nMatch("RELA",chCard) ) 01632 { 01633 lgEOL = true; 01634 sprintf( chHeader, 01635 "#Cumulative emission line relative intensity.\n" ); 01636 } 01637 else 01638 { 01639 sprintf( chHeader, 01640 "#Cumulative emission line absolute intensity.\n" ); 01641 lgEOL = false; 01642 } 01643 /* read in the list of lines to examine */ 01644 punch_line(punch.ipPnunit[punch.npunch],"READ",lgEOL,chTemp); 01645 strcat( chHeader, chTemp ); 01646 } 01647 01648 else if( nMatch("DATA",chCard) ) 01649 { 01650 /* punch line data, done in PunchLineData */ 01651 strcpy( punch.chPunch[punch.npunch], "LIND" ); 01652 sprintf( chHeader, 01653 "#Emission line data.\n" ); 01654 } 01655 01656 else if( nMatch("ARRA",chCard) ) 01657 { 01658 /* punch line array - 01659 * output energies and luminosities of predicted lines */ 01660 strcpy( punch.chPunch[punch.npunch], "LINA" ); 01661 sprintf( chHeader, 01662 "#enr\tID\tI(intrinsic)\tI(emergent)\ttype\n" ); 01663 } 01664 01665 else if( nMatch("LABE",chCard) ) 01666 { 01667 /* punch line labels */ 01668 strcpy( punch.chPunch[punch.npunch], "LINL" ); 01669 sprintf( chHeader, 01670 "#index\tlabel\twavelength\tcomment\n" ); 01671 /* this controls whether we will print lots of redundant 01672 * info labels for transferred lines - if keyword LONG appears 01673 * then do so, if does not appear then do not - this is default */ 01674 if( nMatch("LONG",chCard) ) 01675 punch.punarg[punch.npunch][0] = 1; 01676 else 01677 punch.punarg[punch.npunch][0] = 0; 01678 } 01679 01680 else if( nMatch("OPTI",chCard) ) 01681 { 01682 /* punch line optical depths, done in PunchLineStuff, located in punchdo.c */ 01683 strcpy( punch.chPunch[punch.npunch], "LINO" ); 01684 sprintf( chHeader, 01685 "#species\tenergy\topt depth\tdamp\n" ); 01686 01687 /* the default will be to make wavelengths line in the printout, called labels, 01688 * if units appears then other units will be used instead */ 01689 strcpy( punch.chConPunEnr[punch.npunch], 01690 "labl" ); 01691 01692 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 01693 * units are copied into punch.chConPunEnr */ 01694 if( nMatch("UNIT",chCard) ) 01695 ChkUnits(chCard); 01696 01697 /* this is optional limit to smallest optical depths */ 01698 ipFFmt = 5; 01699 punch.punarg[punch.npunch][0] = (realnum)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL)); 01700 /* this is default of 0.1 napier */ 01701 if( lgEOL ) 01702 { 01703 punch.punarg[punch.npunch][0] = 0.1f; 01704 } 01705 } 01706 01707 else if( nMatch("POPU",chCard) ) 01708 { 01709 /* punch line populations command - first give index and inforamtion 01710 * for all lines, then populations for lines as a function of 01711 * depth, using this index */ 01712 strcpy( punch.chPunch[punch.npunch], "LINP" ); 01713 sprintf( chHeader, 01714 "#population information\n" ); 01715 /* this is optional limit to smallest population to punch - always 01716 * interpreted as a log */ 01717 ipFFmt = 5; 01718 punch.punarg[punch.npunch][0] = (realnum)pow(10.,FFmtRead(chCard,&ipFFmt, INPUT_LINE_LENGTH,&lgEOL)); 01719 01720 /* this is default - all positive populations */ 01721 if( lgEOL ) 01722 punch.punarg[punch.npunch][0] = 0.f; 01723 01724 if( nMatch(" OFF",chCard) ) 01725 { 01726 /* no lower limit - print all lines */ 01727 punch.punarg[punch.npunch][0] = -1.f; 01728 } 01729 } 01730 01731 else if( nMatch("INTE",chCard) ) 01732 { 01733 /* this will be full set of line intensities */ 01734 strcpy( punch.chPunch[punch.npunch], "LINI" ); 01735 sprintf( chHeader, 01736 "#Emission line intensities per unit inner area\n" ); 01737 if( nMatch("COLU",chCard) ) 01738 { 01739 /* column is key to punch single column */ 01740 strcpy( punch.chPunRltType, "column" ); 01741 } 01742 else 01743 { 01744 /* array is key to punch large array */ 01745 strcpy( punch.chPunRltType, "array " ); 01746 } 01747 if( nMatch("EVER",chCard) ) 01748 { 01749 ipFFmt = 3; 01750 punch.LinEvery = (long int)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 01751 punch.lgLinEvery = true; 01752 if( lgEOL ) 01753 { 01754 fprintf( ioQQQ, 01755 "There must be a second number, the number of zones to print.\nSorry.\n" ); 01756 cdEXIT(EXIT_FAILURE); 01757 } 01758 } 01759 else 01760 { 01761 punch.LinEvery = geometry.nend[0]; 01762 punch.lgLinEvery = false; 01763 } 01764 } 01765 else 01766 { 01767 fprintf( ioQQQ, 01768 "This option for PUNCH LINE is something that I do not understand. Sorry.\n" ); 01769 cdEXIT(EXIT_FAILURE); 01770 } 01771 } 01772 01773 else if( nMatch(" MAP",chCard) ) 01774 { 01775 strcpy( punch.chPunch[punch.npunch], "MAP " ); 01776 sprintf( chHeader, 01777 "#te, heating, cooling.\n" ); 01778 /* do cooling space map for specified zones 01779 * if no number, or <0, do map and punch out without doing first zone 01780 * does map by calling punt(" map") 01781 */ 01782 i = 5; 01783 hcmap.MapZone = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01784 if( lgEOL ) 01785 { 01786 hcmap.MapZone = 1; 01787 } 01788 01789 if( nMatch("RANG",chCard) ) 01790 { 01791 bool lgLogOn; 01792 hcmap.RangeMap[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01793 if( hcmap.RangeMap[0] <= 10. && !nMatch("LINE",chCard) ) 01794 { 01795 hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,hcmap.RangeMap[0]); 01796 lgLogOn = true; 01797 } 01798 else 01799 { 01800 lgLogOn = false; 01801 } 01802 01803 hcmap.RangeMap[1] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01804 if( lgLogOn ) 01805 hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,hcmap.RangeMap[1]); 01806 01807 if( lgEOL ) 01808 { 01809 fprintf( ioQQQ, "There must be a zone number, followed by two temperatures, on this line. Sorry.\n" ); 01810 cdEXIT(EXIT_FAILURE); 01811 } 01812 } 01813 } 01814 01815 else if( nMatch("MOLE",chCard) ) 01816 { 01817 /* molecules, especially for PDR calculations */ 01818 strcpy( punch.chPunch[punch.npunch], "MOLE" ); 01819 sprintf( chHeader, 01820 "#molecular species will follow:\n"); 01821 } 01822 01823 else if( nMatch("OPTICAL",chCard) && nMatch("DEPTH",chCard) ) 01824 { 01825 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 01826 * units are copied into punch.chConPunEnr */ 01827 ChkUnits(chCard); 01828 if( nMatch("FINE",chCard) ) 01829 { 01830 /* punch fine continuum optical depths */ 01831 rfield.lgPunchOpacityFine = true; 01832 strcpy( punch.chPunch[punch.npunch], "OPTf" ); 01833 sprintf( chHeader, "#energy\tTau tot\topacity\n" ); 01834 ipFFmt = 5; 01835 /* range option - important since so much data */ 01836 if( nMatch("RANGE",chCard) ) 01837 { 01838 /* get lower and upper range, must be in Ryd */ 01839 punch.punarg[punch.npunch][0] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 01840 punch.punarg[punch.npunch][1] = (realnum)FFmtRead(chCard,&ipFFmt,INPUT_LINE_LENGTH,&lgEOL); 01841 if( lgEOL ) 01842 { 01843 fprintf(ioQQQ,"There must be two numbers, the lower and upper energy range in Ryd.\nSorry.\n"); 01844 cdEXIT(EXIT_FAILURE); 01845 } 01846 if( punch.punarg[punch.npunch][0] >=punch.punarg[punch.npunch][1] ) 01847 { 01848 fprintf(ioQQQ,"The two energies for the range must be in increasing order.\nSorry.\n"); 01849 cdEXIT(EXIT_FAILURE); 01850 } 01851 } 01852 else 01853 { 01854 /* these mean full energy range */ 01855 punch.punarg[punch.npunch][0] = 0.; 01856 punch.punarg[punch.npunch][1] = 0.; 01857 } 01858 /* optional last parameter - how many points to bring together */ 01859 punch.punarg[punch.npunch][2] = (realnum)FFmtRead(chCard,&ipFFmt, 01860 INPUT_LINE_LENGTH,&lgEOL); 01861 /* default is to bring together ten */ 01862 if( lgEOL ) 01863 punch.punarg[punch.npunch][2] = 10; 01864 if( punch.punarg[punch.npunch][2] < 1 ) 01865 { 01866 fprintf(ioQQQ,"The number of fine opacities to skip must be > 0 \nSorry.\n"); 01867 cdEXIT(EXIT_FAILURE); 01868 } 01869 } 01870 else 01871 { 01872 /* punch coarse continuum optical depths */ 01873 strcpy( punch.chPunch[punch.npunch], "OPTc" ); 01874 sprintf( chHeader, 01875 "#energy\ttotal\tabsorp\tscat\n" ); 01876 } 01877 01878 } 01879 else if( nMatch(" OTS",chCard) ) 01880 { 01881 strcpy( punch.chPunch[punch.npunch], " OTS" ); 01882 sprintf( chHeader, 01883 "#otscon, lin, conOpac LinOpc\n" ); 01884 } 01885 01886 else if( nMatch("OVER",chCard) && nMatch(" OVE",chCard) ) 01887 { 01888 /* punch overview of model results */ 01889 strcpy( punch.chPunch[punch.npunch], "OVER" ); 01890 sprintf( chHeader, 01891 "#depth\tTe\tHtot\thden\teden\t2H_2/H\tHI\tHII\tHeI\tHeII\tHeIII\tCO/C\tC1\tC2\tC3\tC4\tO1\tO2\tO3\tO4\tO5\tO6\tAV(point)\tAV(extend)\n" ); 01892 } 01893 01894 else if( nMatch(" PDR",chCard) ) 01895 { 01896 strcpy( punch.chPunch[punch.npunch], " PDR" ); 01897 sprintf( chHeader, 01898 "#depth\tH colden\tTe\tHI/HDEN\tH2/HDEN\tH2*/HDEN\tCI/C\tCO/C\tH2O/O\tG0\tAV(point)\tAV(extend)\tTauV(point)\n" ); 01899 } 01900 01901 else if( nMatch("PHYS",chCard) ) 01902 { 01903 /* punch physical conditions */ 01904 strcpy( punch.chPunch[punch.npunch], "PHYS" ); 01905 sprintf( chHeader, 01906 "#PhyC depth\tTe\tn(H)\tn(e)\tHtot\taccel\tfillfac\n" ); 01907 } 01908 01909 else if( nMatch("POIN",chCard) ) 01910 { 01911 /* punch out the pointers */ 01912 punch.lgPunPoint = true; 01913 /* this does not count as a punch option (really) */ 01914 strcpy( punch.chPunch[punch.npunch], "" ); 01915 punch.lgRealPunch[punch.npunch] = false; 01916 } 01917 01918 else if( nMatch("PRES",chCard) ) 01919 { 01920 /* the punch pressure command */ 01921 strcpy( punch.chPunch[punch.npunch], "PRES" ); 01922 sprintf( chHeader, 01923 "#P depth\tPcorrect\tPcurrent\tPIn+Pinteg\tPgas(r0)\tPgas\tPram\tPrad(line)\tPinteg\tV(wind km/s)\tcad(wind km/s)\tP(mag)\tV(turb km/s)\tP(turb)\tconv?\n" ); 01924 } 01925 01926 else if( nMatch("RADI",chCard) ) 01927 { 01928 /* the punch radius command */ 01929 sprintf( chHeader, "#NZONE\tradius\tdepth\tdr\n" ); 01930 /* option to only punch the outer radius */ 01931 if( nMatch( "OUTE" , chCard ) ) 01932 { 01933 /* only outer radius */ 01934 strcpy( punch.chPunch[punch.npunch], "RADO" ); 01935 } 01936 else 01937 { 01938 /* all radii */ 01939 strcpy( punch.chPunch[punch.npunch], "RADI" ); 01940 } 01941 } 01942 01943 else if( nMatch("RECO",chCard) ) 01944 { 01945 if( nMatch("COEF",chCard) ) 01946 { 01947 /* recombination coefficients for everything */ 01948 01949 /* this is logical flag used in routine ion_recom to create the punch output */ 01950 punch.lgioRecom = true; 01951 /* this does not count as a punch option (really) */ 01952 strcpy( punch.chPunch[punch.npunch], "" ); 01953 punch.lgRealPunch[punch.npunch] = false; 01954 } 01955 01956 else if( nMatch("EFFI",chCard) ) 01957 { 01958 /* punch recombination efficiency */ 01959 strcpy( punch.chPunch[punch.npunch], "RECE" ); 01960 sprintf( chHeader, 01961 "#Recom effic H, Heo, He+\n" ); 01962 } 01963 01964 else 01965 { 01966 fprintf( ioQQQ, "No option recognized on this punch recombination command\n" ); 01967 fprintf( ioQQQ, "Valid options are COEFFICIENTS, AGN, and EFFICIENCY\nSorry.\n" ); 01968 cdEXIT(EXIT_FAILURE); 01969 } 01970 } 01971 01972 /* punch results command, either as single column or wide array */ 01973 else if( nMatch("RESU",chCard) ) 01974 { 01975 strcpy( punch.chPunch[punch.npunch], "RESU" ); 01976 if( nMatch("COLU",chCard) ) 01977 { 01978 /* column is key to punch single column */ 01979 strcpy( punch.chPunRltType, "column" ); 01980 } 01981 else 01982 { 01983 /* array is key to punch large array */ 01984 strcpy( punch.chPunRltType, "array " ); 01985 } 01986 01987 /* do not change following, is used as flag in getlines */ 01988 sprintf( chHeader, 01989 "#results of calculation\n" ); 01990 } 01991 01992 else if( nMatch("SECO",chCard) ) 01993 { 01994 /* punch secondary ionization rate */ 01995 strcpy( punch.chPunch[punch.npunch], "SECO" ); 01996 sprintf( chHeader, 01997 "#depth\tIon(H^0)\tDiss(H_2)\tExcit(Lya)\n" ); 01998 } 01999 02000 else if( nMatch("SOUR",chCard) ) 02001 { 02002 if( nMatch("DEPT",chCard) ) 02003 { 02004 /* print continuum source function as function of depth */ 02005 strcpy( punch.chPunch[punch.npunch], "SOUD" ); 02006 sprintf( chHeader, 02007 "#continuum source function vs depth\n" ); 02008 } 02009 else if( nMatch("SPEC",chCard) ) 02010 { 02011 /* print spectrum continuum source function at 1 depth */ 02012 strcpy( punch.chPunch[punch.npunch], "SOUS" ); 02013 sprintf( chHeader, 02014 "#continuum source function\n" ); 02015 } 02016 else 02017 { 02018 fprintf( ioQQQ, "A second keyword must appear on this line.\n" ); 02019 fprintf( ioQQQ, "They are DEPTH and SPECTRUM.\n" ); 02020 fprintf( ioQQQ, "Sorry.\n" ); 02021 cdEXIT(EXIT_FAILURE); 02022 } 02023 } 02024 02025 02026 /* punch spectrum the new form of the punch continuum, will eventually replace the standard 02027 * punch continuum command */ 02028 else if( nMatch("SPEC",chCard) && nMatch("TRUM",chCard) && !nMatch("XSPE",chCard) ) 02029 { 02030 /* this flag is checked in PrtComment to generate a caution 02031 * if continuum is punched but iterations not performed */ 02032 punch.lgPunContinuum = true; 02033 02034 /* set flag for spectrum */ 02035 strcpy( punch.chPunch[punch.npunch], "CONN" ); 02036 02037 sprintf( chHeader, 02038 "#Cont Enr\tincid nFn\ttrans\tdiff\tlines \n" ); 02039 02040 /* check for keyword UNITS on line, then scan wavelength or energy units if present, 02041 * units are copied into punch.chConPunEnr */ 02042 ChkUnits(chCard); 02043 02044 /* this points to which punch new continuum this is - 02045 * parameters were stored in previous set spectrum commands */ 02046 punch.punarg[punch.npunch][0] = (realnum)punch.cp_npun; 02047 02048 ++punch.cp_npun; 02049 /* check that limit not exceeded */ 02050 if( punch.cp_npun > LIMPUN ) 02051 { 02052 fprintf( ioQQQ, 02053 "The limit to the number of PUNCH options is %ld. Increase LIMPUN in punch.h if more are needed.\nSorry.\n", 02054 LIMPUN ); 02055 cdEXIT(EXIT_FAILURE); 02056 } 02057 02058 } 02059 02060 else if( nMatch("SPEC",chCard) && nMatch("CIAL",chCard) ) 02061 { 02062 /* punch special, will call routine PunchSpecial */ 02063 strcpy( punch.chPunch[punch.npunch], "SPEC" ); 02064 sprintf( chHeader, "#Special.\n" ); 02065 } 02066 02067 else if( nMatch("TEGR",chCard) ) 02068 { 02069 /* punch tegrid */ 02070 strcpy( punch.chPunch[punch.npunch], "TEGR" ); 02071 sprintf( chHeader, 02072 "#zone, te, heat, cool.\n" ); 02073 } 02074 02075 else if( nMatch("TEMP",chCard) ) 02076 { 02077 /* punch temperature command */ 02078 strcpy( punch.chPunch[punch.npunch], "TEMP" ); 02079 sprintf( chHeader, 02080 "#depth\tTe\tcC/dT\tdt/dr\td^2T/dr^2\n" ); 02081 } 02082 02083 else if( nMatch("TIME",chCard) && nMatch("DEPE",chCard) ) 02084 { 02085 /* information about time dependent solutions */ 02086 strcpy( punch.chPunch[punch.npunch], "TIMD" ); 02087 /* do not want to separate iterations with special character */ 02088 punch.lg_separate_iterations[punch.npunch] = false; 02089 /* write header */ 02090 sprintf( chHeader , 02091 "#elapsed time\ttime step \tscale cont\tn(H)\t<T>\t<H+/H rad>\t<H0/H rad>\t<H2/H rad>\t<He+/H rad>\t<CO/H>\n" ); 02092 } 02093 02094 else if( nMatch("TPRE",chCard) ) 02095 { 02096 /* debug output from the temperature predictor in zonestart, 02097 * set with punch tpred command */ 02098 strcpy( punch.chPunch[punch.npunch], "TPRE" ); 02099 sprintf( chHeader, 02100 "#zone old temp, guess Tnew, new temp delta \n" ); 02101 } 02102 02103 else if( nMatch("WIND",chCard) ) 02104 { 02105 strcpy( punch.chPunch[punch.npunch], "WIND" ); 02106 sprintf( chHeader, 02107 "#radius\tdepth\tvel [cm/s]\tTot accel [cm s-2]\tLin accel [cm s-2]" 02108 "\tCon accel [cm s-2]\tforce multiplier\ta_gravity\n" ); 02109 if( nMatch( "TERM" , chCard ) ) 02110 { 02111 /* only punch for last zone, the terminal velocity, for grids */ 02112 punch.punarg[punch.npunch][0] = 0.; 02113 } 02114 else 02115 { 02116 /* one means punch every zone */ 02117 punch.punarg[punch.npunch][0] = 1.; 02118 } 02119 } 02120 02121 else if( nMatch("XSPE",chCard) ) 02122 { 02123 /* say that this is a FITS file output */ 02124 punch.lgFITS[punch.npunch] = true; 02125 02126 /* the punch xspec commands */ 02127 if( !punch.lgPunLstIter[punch.npunch] ) 02128 { 02129 punch.lgPunLstIter[punch.npunch] = true; 02130 } 02131 02132 if( nMatch("ATAB",chCard) ) 02133 { 02134 /* punch xspec atable command */ 02135 if( nMatch("INCI",chCard) ) 02136 { 02137 if( nMatch("ATTE",chCard) ) 02138 { 02139 /* attenuated incident continuum */ 02140 strcpy( punch.chPunch[punch.npunch], "XATT" ); 02141 grid.lgOutputTypeOn[2] = true; 02142 } 02143 else if( nMatch("REFL",chCard) ) 02144 { 02145 /* reflected incident continuum */ 02146 strcpy( punch.chPunch[punch.npunch], "XRFI" ); 02147 grid.lgOutputTypeOn[3] = true; 02148 } 02149 else 02150 { 02151 /* incident continuum */ 02152 strcpy( punch.chPunch[punch.npunch], "XINC" ); 02153 grid.lgOutputTypeOn[1] = true; 02154 } 02155 } 02156 else if( nMatch("DIFF",chCard) ) 02157 { 02158 if( nMatch("REFL",chCard) ) 02159 { 02160 /* reflected diffuse continuous emission */ 02161 strcpy( punch.chPunch[punch.npunch], "XDFR" ); 02162 grid.lgOutputTypeOn[5] = true; 02163 } 02164 else 02165 { 02166 /* diffuse continuous emission outward */ 02167 strcpy( punch.chPunch[punch.npunch], "XDFO" ); 02168 grid.lgOutputTypeOn[4] = true; 02169 } 02170 } 02171 else if( nMatch("LINE",chCard) ) 02172 { 02173 if( nMatch("REFL",chCard) ) 02174 { 02175 /* reflected lines */ 02176 strcpy( punch.chPunch[punch.npunch], "XLNR" ); 02177 grid.lgOutputTypeOn[7] = true; 02178 } 02179 else 02180 { 02181 /* outward lines */ 02182 strcpy( punch.chPunch[punch.npunch], "XLNO" ); 02183 grid.lgOutputTypeOn[6] = true; 02184 } 02185 } 02186 else if( nMatch("SPEC",chCard) ) 02187 { 02188 if( nMatch("REFL",chCard) ) 02189 { 02190 /* reflected spectrum */ 02191 strcpy( punch.chPunch[punch.npunch], "XREF" ); 02192 grid.lgOutputTypeOn[9] = true; 02193 } 02194 else 02195 { 02196 /* transmitted spectrum */ 02197 strcpy( punch.chPunch[punch.npunch], "XTRN" ); 02198 grid.lgOutputTypeOn[8] = true; 02199 } 02200 } 02201 else 02202 { 02203 /* transmitted spectrum */ 02204 strcpy( punch.chPunch[punch.npunch], "XTRN" ); 02205 grid.lgOutputTypeOn[8] = true; 02206 } 02207 } 02208 else if( nMatch("MTAB",chCard) ) 02209 { 02210 /* punch xspec mtable */ 02211 strcpy( punch.chPunch[punch.npunch], "XSPM" ); 02212 grid.lgOutputTypeOn[10] = true; 02213 } 02214 else 02215 { 02216 fprintf( ioQQQ, "Support only for xspec atable and xspec mtable.\n" ); 02217 cdEXIT( EXIT_FAILURE ); 02218 } 02219 } 02220 02221 /* punch column density has to come last so do not trigger specific column 02222 * densities, H2, FeII, etc. 02223 * Need both keywords since column is also the keyword for one line per line */ 02224 else if( nMatch("COLU",chCard) && nMatch("DENS",chCard) ) 02225 { 02226 if( nMatch("SOME" , chCard )) 02227 { 02228 /* flag saying punch some column densities */ 02229 strcpy( punch.chPunch[punch.npunch], "COLS" ); 02230 punch_colden( chHeader , punch.ipPnunit[punch.npunch] , "READ" ); 02231 } 02232 else 02233 { 02234 /* punch them all as in table */ 02235 /* punch column density */ 02236 strcpy( punch.chPunch[punch.npunch], "COLU" ); 02237 sprintf( chHeader, 02238 "#column densities\n" ); 02239 } 02240 } 02241 else 02242 { 02243 fprintf( ioQQQ, 02244 "ParsePunch cannot find a recognized keyword on this PUNCH command line.\nSorry.\n" ); 02245 cdEXIT(EXIT_FAILURE); 02246 } 02247 02248 /* only open if file has not already been opened during a previous call */ 02249 if( punch.ipPnunit[punch.npunch] == NULL ) 02250 { 02251 if( nMatch("XSPE",chCard) || nMatch("FITS",chCard) ) 02252 { 02253 /* open the file with the name generated above */ 02254 punch.ipPnunit[punch.npunch] = open_data( chFilename, "wb", AS_LOCAL_ONLY ); 02255 } 02256 else if( punch.lgPunchToSeparateFiles[punch.npunch] && optimize.lgOptimr ) 02257 { 02258 char chIndex[4]; 02259 02260 /* take max because optimize.nOptimiz has not yet been incremented on first pass. */ 02261 sprintf( chIndex, "%li", MAX2( 1, optimize.nOptimiz ) ); 02262 02263 strcat( chFilename, chIndex ); 02264 02265 fixit(); 02266 /* must create unique name here. considerations include: 02267 * 1. must zero fill to the left 02268 * 2. must make sure that we have enough digits (don't hard-wire 3 or 4 digits, 02269 count the number of terms. 02270 * 3. should we look for a file ending (i.e., ".txt" and insert number before that? */ 02271 02272 /* open the file with the name generated above */ 02273 punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY ); 02274 } 02275 else 02276 { 02277 /* open the file with the name generated above */ 02278 punch.ipPnunit[punch.npunch] = open_data( chFilename, "w", AS_LOCAL_ONLY ); 02279 } 02280 02281 /* option to set no buffering for this file. The setbuf command may 02282 * ONLY be issued right after the open of the file. Giving it after 02283 * i/o has been done may result in loss of the contents of the buffer, PvH */ 02284 if( nMatch("NO BUFFER",chCard) ) 02285 setbuf( punch.ipPnunit[punch.npunch] , NULL ); 02286 } 02287 02288 /***************************************************************/ 02289 /* */ 02290 /* The following are special punch options and must be done */ 02291 /* after the parsing and file opening above. */ 02292 /* */ 02293 /* NB: these are ALSO parsed above. Here we DO something. */ 02294 /* */ 02295 /***************************************************************/ 02296 02297 if( nMatch("CONV",chCard) && nMatch("REAS",chCard) ) 02298 { 02299 /* punch reason model declared not converged 02300 * not a true punch command, since done elsewhere */ 02301 punch.ipPunConv = punch.ipPnunit[punch.npunch]; 02302 lgPunConv_noclobber = lgNoClobber[punch.npunch]; 02303 punch.lgPunConv = true; 02304 fprintf( punch.ipPunConv, 02305 "# reason for continued iterations\n" ); 02306 strcpy( punch.chPunch[punch.npunch], "" ); 02307 punch.lgRealPunch[punch.npunch] = false; 02308 } 02309 02310 else if( nMatch("CONV",chCard) && nMatch("BASE",chCard) ) 02311 { 02312 /* punch some quantities we are converging */ 02313 punch.lgTraceConvergeBase = true; 02314 /* the second punch occurrence - file has been opened, 02315 * copy handle, also pass on special no hash option */ 02316 if( nMatch("O HA",chCard) ) 02317 punch.lgTraceConvergeBaseHash = false; 02318 punch.ipTraceConvergeBase = punch.ipPnunit[punch.npunch]; 02319 /* set punch last flag to whatever it was above */ 02320 lgTraceConvergeBase_noclobber = lgNoClobber[punch.npunch]; 02321 static bool lgPrtHeader = true; 02322 if( lgPrtHeader ) 02323 fprintf( punch.ipTraceConvergeBase, 02324 "#zone\theat\tcool\teden\n" ); 02325 lgPrtHeader = false; 02326 } 02327 02328 else if( nMatch(" DR ",chCard) ) 02329 { 02330 static bool lgPrtHeader = true; 02331 /* the second punch dr occurrence - file has been opened, 02332 * copy handle to ipDRout, also pass on special no hash option */ 02333 if( nMatch("O HA",chCard) ) 02334 punch.lgDRHash = false; 02335 punch.ipDRout = punch.ipPnunit[punch.npunch]; 02336 /* set punch last flag to whatever it was above */ 02337 punch.lgDRPLst = punch.lgPunLstIter[punch.npunch]; 02338 lgDROn_noclobber = lgNoClobber[punch.npunch]; 02339 if( lgPrtHeader ) 02340 fprintf( punch.ipDRout, 02341 "#zone\tdepth\tdr\tdr 2 go\treason \n" ); 02342 lgPrtHeader = false; 02343 strcpy( punch.chPunch[punch.npunch], "" ); 02344 punch.lgRealPunch[punch.npunch] = false; 02345 } 02346 02347 else if( nMatch("QHEA",chCard) ) 02348 { 02349 gv.QHPunchFile = punch.ipPnunit[punch.npunch]; 02350 gv.lgQHPunLast = punch.lgPunLstIter[punch.npunch]; 02351 lgQHPunchFile_noclobber = lgNoClobber[punch.npunch]; 02352 fprintf( gv.QHPunchFile, 02353 "#Probability distributions from quantum heating routine.\n" ); 02354 } 02355 02356 else if( nMatch("POIN",chCard) ) 02357 { 02358 /* punch out the pointers */ 02359 punch.ipPoint = punch.ipPnunit[punch.npunch]; 02360 lgPunPoint_noclobber = lgNoClobber[punch.npunch]; 02361 punch.lgPunPoint = true; 02362 fprintf( punch.ipPoint, 02363 "#pointers. \n" ); 02364 strcpy( punch.chPunch[punch.npunch], "" ); 02365 punch.lgRealPunch[punch.npunch] = false; 02366 } 02367 02368 else if( nMatch("RECO",chCard) && nMatch("COEF",chCard) ) 02369 { 02370 /* recombination coefficients for everything 02371 * punch.lgioRecom set to false in routine zero, non-zero value 02372 * is flag to punch recombination coefficients. the output is actually 02373 * produced by a series of routines, as they generate the recombination 02374 * coefficients. these include 02375 * diel supres, helium, hydrorecom, iibod, and makerecomb*/ 02376 punch.ioRecom = punch.ipPnunit[punch.npunch]; 02377 lgioRecom_noclobber = lgNoClobber[punch.npunch]; 02378 /* this is logical flag used in routine ion_recom to create the punch output */ 02379 punch.lgioRecom = true; 02380 fprintf( punch.ioRecom, 02381 "#recombination coefficients cm3 s-1 for current density and temperature\n" ); 02382 strcpy( punch.chPunch[punch.npunch], "" ); 02383 punch.lgRealPunch[punch.npunch] = false; 02384 } 02385 02386 else if( nMatch(" MAP",chCard) ) 02387 { 02388 /* say output goes to special punch */ 02389 ioMAP = punch.ipPnunit[punch.npunch]; 02390 } 02391 02392 /* check that string written into chHeader can actually fit there 02393 * we may have overrun this buffer, an internal error */ 02394 /* check that there are less than nChar characters in the line */ 02395 char *chEOL = strchr(chHeader , '\0' ); 02396 02397 /* return null if input string longer than nChar, the longest we can read. 02398 * Print and return null but chLine still has as much of the line as 02399 * could be placed in cdLine */ 02400 if( (chEOL==NULL) || (chEOL - chHeader)>=MAX_HEADER_SIZE-1 ) 02401 { 02402 fprintf( ioQQQ, "DISASTER chHeader has been overwritten " 02403 "with a line too long to be read.\n" ); 02404 cdEXIT(EXIT_FAILURE); 02405 } 02406 02407 /* if flag set and cdHeader not still equal to initialized string, print header 02408 * punch.lgPunHeader normally true, is false during grid calculation */ 02409 if( nSimThisCoreload > 1 ) 02410 punch.lgPunHeader = false; 02411 if( punch.lgPunHeader && !nMatch("ArNdY38dZ9us4N4e12SEcuQ",chHeader) ) 02412 { 02413 fprintf( punch.ipPnunit[punch.npunch], "%s", chHeader ); 02414 } 02415 02416 /* increment total number of punch commands, */ 02417 ++punch.npunch; 02418 return; 02419 } 02420 02421 /*PunchFilesInit initialize punch file pointers, called from InitCoreload 02422 * called one time per core load 02423 * NB KEEP THIS ROUTINE SYNCHED UP WITH THE NEXT ONE, ClosePunchFiles */ 02424 void PunchFilesInit(void) 02425 { 02426 long int i; 02427 static bool lgFIRST = true; 02428 02429 DEBUG_ENTRY( "PunchFilesInit()" ); 02430 02431 ASSERT( lgFIRST ); 02432 lgFIRST = false; 02433 02434 /* set lgNoClobber to not overwrite files, reset with clobber on punch line 02435 * if we are running a grid (grid command entered in cdRead) grid.lgGrid 02436 * true, is false if single sim. For grid we want to not clobber files 02437 * by default, do clobber for optimizer since this was behavior before */ 02438 bool lgNoClobberDefault = false; 02439 if( grid.lgGrid ) 02440 { 02441 /* cdRead encountered grid command - do not want to clobber files */ 02442 lgNoClobberDefault = true; 02443 } 02444 02445 for( i=0; i < LIMPUN; i++ ) 02446 { 02447 lgNoClobber[i] = lgNoClobberDefault; 02448 } 02449 lgPunConv_noclobber = lgNoClobberDefault; 02450 lgDROn_noclobber = lgNoClobberDefault; 02451 lgTraceConvergeBase_noclobber = lgNoClobberDefault; 02452 lgPunPoint_noclobber = lgNoClobberDefault; 02453 lgioRecom_noclobber = lgNoClobberDefault; 02454 lgQHPunchFile_noclobber = lgNoClobberDefault; 02455 02456 /* done parsing, now turn flag to punch headers, the default behavior 02457 * this is set false during a grid calculation to do only one header per grid*/ 02458 punch.lgPunHeader = true; 02459 02460 for( i=0; i < LIMPUN; i++ ) 02461 { 02462 if( !lgNoClobber[i] ) 02463 { 02464 punch.ipPnunit[i] = NULL; 02465 } 02466 /* set false with the dummy punch commands like punch dr */ 02467 punch.lgRealPunch[i] = true; 02468 /* this sets energy range of continuum, zero says to do full continuum */ 02469 punch.cp_range_min[i] = 0.; 02470 punch.cp_range_max[i] = 0.; 02471 /* this means to keep native resolution of code, reset to another 02472 * resolution with set punch resolution command */ 02473 punch.cp_resolving_power[i] = 0.; 02474 } 02475 02476 punch.lgTraceConvergeBase = false; 02477 02478 if( !lgDROn_noclobber ) 02479 { 02480 punch.ipDRout = NULL; 02481 punch.lgDROn = false; 02482 } 02483 02484 if( !lgTraceConvergeBase_noclobber ) 02485 { 02486 punch.ipTraceConvergeBase = NULL; 02487 punch.lgTraceConvergeBase = false; 02488 } 02489 02490 if( !lgPunConv_noclobber ) 02491 { 02492 punch.ipPunConv = NULL; 02493 punch.lgPunConv = false; 02494 } 02495 02496 if( !lgPunPoint_noclobber ) 02497 { 02498 punch.ipPoint = NULL; 02499 punch.lgPunPoint = false; 02500 } 02501 02502 if( !lgQHPunchFile_noclobber ) 02503 gv.QHPunchFile = NULL; 02504 02505 if( !lgioRecom_noclobber ) 02506 { 02507 punch.ioRecom = NULL; 02508 punch.lgioRecom = false; 02509 } 02510 02511 ioMAP = NULL; 02512 return; 02513 } 02514 02515 /*ClosePunchFiles close punch files called from cdEXIT upon termination, 02516 * from cloudy before returning 02517 * NB - KEEP THIS ROUTINE SYNCHED UP WITH THE PREVIOUS ONE, PunchFilesInit */ 02518 void ClosePunchFiles( bool lgFinal ) 02519 { 02520 long int i; 02521 02522 DEBUG_ENTRY( "ClosePunchFiles()" ); 02523 02524 /* close all punch units cloudy opened with punch command, 02525 * lgNoClobber is set false with CLOBBER option on punch, says to 02526 * overwrite the files */ 02527 for( i=0; i < punch.npunch; i++ ) 02528 { 02529 /* if lgFinal is true, we close everything, no matter what. 02530 * this means ignoring "no clobber" options */ 02531 if( punch.ipPnunit[i] != NULL && ( !lgNoClobber[i] || lgFinal ) ) 02532 { 02533 /* Test that any FITS files are the right size! */ 02534 if( punch.lgFITS[i] ) 02535 { 02536 /* \todo 2 This overflows for file sizes larger (in bytes) than 02537 * a long int can represent (about 2GB on most 2007 systems) */ 02538 fseek(punch.ipPnunit[i], 0, SEEK_END); 02539 long file_size = ftell(punch.ipPnunit[i]); 02540 if( file_size%2880 ) 02541 { 02542 fprintf( ioQQQ, " PROBLEM FITS file is wrong size!\n" ); 02543 } 02544 } 02545 02546 fclose( punch.ipPnunit[i] ); 02547 punch.ipPnunit[i] = NULL; 02548 } 02549 } 02550 02551 /* following file handles are aliased to ipPnunit which 02552 * was closed above */ 02553 if( punch.ipDRout != NULL && !lgDROn_noclobber ) 02554 { 02555 /*fclose( punch.ipDRout );*/ 02556 punch.ipDRout = NULL; 02557 punch.lgDROn = false; 02558 } 02559 02560 if( punch.ipTraceConvergeBase != NULL && !lgTraceConvergeBase_noclobber ) 02561 { 02562 /*fclose( punch.ipDRout );*/ 02563 punch.ipTraceConvergeBase = NULL; 02564 punch.lgTraceConvergeBase = false; 02565 } 02566 02567 if( punch.ipPunConv != NULL && !lgPunConv_noclobber ) 02568 { 02569 /*fclose( punch.ipPunConv );*/ 02570 punch.ipPunConv = NULL; 02571 punch.lgPunConv = false; 02572 } 02573 if( punch.ipPoint != NULL && !lgPunPoint_noclobber ) 02574 { 02575 /*fclose( punch.ipPoint );*/ 02576 punch.ipPoint = NULL; 02577 punch.lgPunPoint = false; 02578 } 02579 if( gv.QHPunchFile != NULL && !lgQHPunchFile_noclobber ) 02580 { 02581 /*fclose( gv.QHPunchFile );*/ 02582 gv.QHPunchFile = NULL; 02583 } 02584 if( punch.ioRecom != NULL && !lgioRecom_noclobber ) 02585 { 02586 /*fclose( punch.ioRecom );*/ 02587 punch.ioRecom = NULL; 02588 punch.lgioRecom = false; 02589 } 02590 /* this file was already closed as a punch.ipPnunit, erase copy of pointer as well */ 02591 ioMAP = NULL; 02592 return; 02593 } 02594 02595 /*ChkUnits check for keyword UNITS on line, then scan wavelength or energy units if present, 02596 * units are copied into punch.chConPunEnr - when doing output, the routine call 02597 * AnuUnit( energy ) will automatically return the energy in the right units, 02598 * when called to do punch output */ 02599 STATIC void ChkUnits( char *chCard) 02600 { 02601 02602 DEBUG_ENTRY( "ChkUnits()" ); 02603 02604 /* option to set units for continuum energy in punch output */ 02605 if( nMatch("NITS",chCard) ) 02606 { 02607 if( nMatch("MICR",chCard) ) 02608 { 02609 /* microns */ 02610 strcpy( punch.chConPunEnr[punch.npunch], "micr" ); 02611 } 02612 else if( nMatch(" KEV",chCard) ) 02613 { 02614 /* keV */ 02615 strcpy( punch.chConPunEnr[punch.npunch], " kev" ); 02616 } 02617 else if( nMatch("CENT",chCard) || nMatch(" CM ",chCard) ) 02618 { 02619 /* centimeters*/ 02620 strcpy( punch.chConPunEnr[punch.npunch], "cent" ); 02621 } 02622 else if( nMatch(" EV ",chCard) ) 02623 { 02624 /* eV */ 02625 strcpy( punch.chConPunEnr[punch.npunch], " ev " ); 02626 } 02627 else if( nMatch("ANGS",chCard) ) 02628 { 02629 /* angstroms */ 02630 strcpy( punch.chConPunEnr[punch.npunch], "angs" ); 02631 } 02632 else if( nMatch("WAVE",chCard) ) 02633 { 02634 /* wavenumbers */ 02635 strcpy( punch.chConPunEnr[punch.npunch], "wave" ); 02636 } 02637 else if( nMatch(" MHZ",chCard) ) 02638 { 02639 /* megaherz */ 02640 strcpy( punch.chConPunEnr[punch.npunch], " mhz" ); 02641 } 02642 else if( nMatch(" RYD",chCard) ) 02643 { 02644 /* the noble Rydberg */ 02645 /* >>chng 01 sep 02, had been rdyb unlike proper ryd */ 02646 strcpy( punch.chConPunEnr[punch.npunch], "ryd " ); 02647 } 02648 else 02649 { 02650 fprintf( ioQQQ, "I did not recognize units on this line.\n" ); 02651 fprintf( ioQQQ, "Units are keV, eV, Angstroms, Rydbergs, centimeters, and microns.\nSorry.\n" ); 02652 cdEXIT(EXIT_FAILURE); 02653 } 02654 } 02655 else 02656 { 02657 strcpy( punch.chConPunEnr[punch.npunch], "ryd " ); 02658 } 02659 return; 02660 }