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 /*ParseSet scan parameters off SET command */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "geometry.h" 00007 #include "input.h" 00008 #include "prt.h" 00009 #include "mole.h" 00010 #include "rt.h" 00011 #include "phycon.h" 00012 #include "optimize.h" 00013 #include "hcmap.h" 00014 #include "hmi.h" 00015 #include "dense.h" 00016 #include "h2.h" 00017 #include "iterations.h" 00018 #include "conv.h" 00019 #include "secondaries.h" 00020 #include "rfield.h" 00021 #include "ionbal.h" 00022 #include "numderiv.h" 00023 #include "dynamics.h" 00024 #include "iso.h" 00025 #include "punch.h" 00026 #include "stopcalc.h" 00027 #include "opacity.h" 00028 #include "hydrogenic.h" 00029 #include "peimbt.h" 00030 #include "radius.h" 00031 #include "atmdat.h" 00032 #include "continuum.h" 00033 #include "grains.h" 00034 #include "grainvar.h" 00035 #include "parse.h" 00036 #include "lines.h" 00037 #include "assertresults.h" 00038 #include "thirdparty.h" 00039 00040 void ParseSet(char *chCard ) 00041 { 00042 bool lgEOL; 00043 long int i, 00044 ip; 00045 char chString_quotes_lowercase[INPUT_LINE_LENGTH]; 00046 bool lgQuotesFound; 00047 00048 DEBUG_ENTRY( "ParseSet()" ); 00049 00050 /* first check for any strings within quotes - this will get the string 00051 * and blank it out, so that these are not confused with keywords. if 00052 * double quotes not present routine returns unity, zero if found*/ 00053 lgQuotesFound = true; 00054 if( GetQuote( chString_quotes_lowercase , chCard , false ) ) 00055 lgQuotesFound = false; 00056 00057 /* commands to set certain variables, "SET XXX=" */ 00058 if( nMatch("ASSE",chCard) && nMatch("ABOR",chCard) ) 00059 { 00060 /* set assert abort command, to tell the code to raise SIGABRT when an 00061 * assert is thrown - this can be caught by the debugger. The keyword 00062 * _FPE is accepted for backward compatibility */ 00063 cpu.setAssertAbort( true ); 00064 } 00065 00066 else if( nMatch("ASSE",chCard) &&nMatch("SCIE",chCard) ) 00067 { 00068 /* print asserts using scientific notation. Useful when an asserted value is 00069 * less than 10^-4 of the normalization line. */ 00070 lgPrtSciNot = true; 00071 } 00072 00073 else if( nMatch("ATOM",chCard) && nMatch( "DATA", chCard ) ) 00074 { 00075 /* set atomic data */ 00076 long int nelem , ion; 00077 bool lgAs=false , lgCS=false , lgDR=false; 00078 const char *chMole; 00079 00080 /* As? */ 00081 if( nMatch( " AS " , chCard ) ) 00082 lgAs = true; 00083 /* DR - dielectronic recombination */ 00084 else if( nMatch( " DR " , chCard ) ) 00085 lgDR = true; 00086 /* lgCS collision strengths */ 00087 else if( nMatch( " CS " , chCard ) ) 00088 lgCS = true; 00089 else 00090 { 00091 /* no keyword */ 00092 fprintf( ioQQQ, " One of the keywords AS DR or CS must appear to change" 00093 " the transition probabilities, dielectronic recombination rate," 00094 " or collision strength.\n Sorry.\n" ); 00095 cdEXIT(EXIT_FAILURE); 00096 } 00097 00098 /* ion or molecule? */ 00099 if( nMatch(" ION",chCard) ) 00100 { 00101 /* set atomic data - must have an element name and possibly an 00102 * ionization stage - nelem is atomic number on C scale */ 00103 if( (nelem = GetElem(chCard))<0 ) 00104 { 00105 fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" ); 00106 cdEXIT(EXIT_FAILURE); 00107 } 00108 i = 5; 00109 /* make sure that the ionization stage is ok if it is specified 00110 * this is entered on physics scale, 1 is atom, 2 first ion */ 00111 ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00112 if( !lgEOL && (ion< 1 || ion>nelem+1) ) 00113 { 00114 fprintf( ioQQQ, " The ionization stage is not valid\n Sorry.\n" ); 00115 cdEXIT(EXIT_FAILURE); 00116 } 00117 00118 /* two possible sets of transition probabilities for OII */ 00119 if( nelem==ipOXYGEN && ion == 2 && lgAs ) 00120 { 00121 if( nMatch( " NEW" , chCard ) ) 00122 { 00123 dense.lgAsChoose[nelem][ion-1] = true; 00124 } 00125 else if( nMatch( " OLD" , chCard ) ) 00126 { 00127 /*>>chng 06 mar 30, this is the default */ 00128 dense.lgAsChoose[nelem][ion-1] = false; 00129 } 00130 else 00131 { 00132 fprintf( ioQQQ, " The keyword old or new must appear\n Sorry.\n" ); 00133 cdEXIT(EXIT_FAILURE); 00134 } 00135 } 00136 else if( nelem==ipSULPHUR && lgDR ) 00137 { 00138 /* change DR data set for S 00139 * three cases for S DR - ionbal.nDR_S_guess 00140 * 0, default larger of guess and Badnell 00141 * 1, pure Badnell 00142 * 3, scaled oxygen */ 00143 if( nMatch( " MIX" , chCard ) ) 00144 /* this is the default, larger of Badnell or guess */ 00145 ionbal.nDR_S_guess = 0; 00146 else if( nMatch( "PURE" , chCard ) ) 00147 /* use only Badnell 1991 */ 00148 ionbal.nDR_S_guess = 1; 00149 else if( nMatch( "SCAL" , chCard ) ) 00150 { 00151 /* scaled oxygen */ 00152 ionbal.nDR_S_guess = 2; 00153 i = 5; 00154 ionbal.DR_S_scale[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00155 if( lgEOL ) 00156 NoNumb(chCard); 00157 for( ion=1; ion<5; ++ion ) 00158 { 00159 /* optionally get rest - if too few specified then use 00160 * previously set value */ 00161 ionbal.DR_S_scale[ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00162 if( lgEOL ) 00163 ionbal.DR_S_scale[ion] = ionbal.DR_S_scale[ion-1]; 00164 } 00165 } 00166 else 00167 { 00168 /* disaster - no match */ 00169 fprintf( 00170 ioQQQ, " One of the keywords MIX, PURE, or SCALe must" 00171 "appear on the set atomic physics sulphur dr command.\n" 00172 "Sorry.\n" ); 00173 cdEXIT(EXIT_FAILURE); 00174 } 00175 } 00176 else 00177 { 00178 fprintf( ioQQQ, " None of the valid set atomic data atoms/ions were found\n Sorry.\n" ); 00179 cdEXIT(EXIT_FAILURE); 00180 } 00181 } 00182 else 00183 { 00184 /* a molecule */ 00185 if( nMatch(" H2 ",chCard) ) 00186 { 00187 /* H2 */ 00188 chMole = "H2"; 00189 } 00190 else 00191 { 00192 /* nothing recognized */ 00193 fprintf( ioQQQ, " No molecule was on this SET ATOMIC DATA command.\n Sorry.\n" ); 00194 fprintf( ioQQQ, " Use SET ATOMIC DATA ION to change an ion.\n Sorry.\n" ); 00195 cdEXIT(EXIT_FAILURE); 00196 } 00197 00198 if( strcmp( chMole , "H2" )==0 ) 00199 { 00200 if( nMatch(" H " , chCard ) && lgCS ) 00201 { 00202 /* change the H - H2 collision data set */ 00203 i = 5; 00204 ion = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00205 if( ion!=2 ) 00206 TotalInsanity(); 00207 long int nYear = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00208 if( lgEOL ) 00209 NoNumb( chCard ); 00210 if( nYear == 1999 ) 00211 /* use the coefficients from 00212 *>>refer H2 H collision Le Bourlot, J., Pineau des Forets, 00213 *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802 */ 00214 h2.lgH2_H_coll_07 = false; 00215 else if( nYear == 2007 ) 00216 /* use the coefficients from 00217 *>>refer H2 H collision Wrathmall, S. A., Gusdorf A., 00218 *>>refercon & Flower, D.R. 2007, MNRAS, 382, 133 */ 00219 h2.lgH2_H_coll_07 = true; 00220 else 00221 { 00222 /* not an option */ 00223 fprintf(ioQQQ," the SET ATOMIC DATA MOLECULE H2" 00224 " H CS command must have year 1999 or 2007.\n" ); 00225 cdEXIT(EXIT_FAILURE); 00226 } 00227 } 00228 else if( nMatch(" He" , chCard ) && lgCS ) 00229 { 00230 /* change the He - H2 collision data set */ 00231 if( nMatch("ORNL" , chCard ) ) 00232 { 00233 /* use the coefficients from 00234 *>>refer H2 He collision Lee et al in preparation */ 00235 mole.lgH2_He_ORNL = true; 00236 } 00237 else if( nMatch( "BOUR",chCard ) ) 00238 { 00239 /* use the coefficients from 00240 *>>refer H2 He collision Le Bourlot, J., Pineau des Forets, 00241 *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802*/ 00242 mole.lgH2_He_ORNL = false; 00243 } 00244 else 00245 { 00246 /* not an option */ 00247 fprintf(ioQQQ," the SET ATOMIC DATA MOLECULE H2" 00248 "He CS command must have ORNL or Le BOURlot.\n" ); 00249 cdEXIT(EXIT_FAILURE); 00250 } 00251 } 00252 } 00253 else 00254 TotalInsanity(); 00255 } 00256 } 00257 00258 else if( nMatch(" CHA",chCard) && !nMatch( "HO ", chCard ) ) 00259 { 00260 /* set log of minimum charge transfer rate for high ions and H 00261 * default of 1.92e-10 set in zero */ 00262 i = 5; 00263 atmdat.HCTAlex = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00264 if( lgEOL ) 00265 { 00266 NoNumb(chCard); 00267 } 00268 if( atmdat.HCTAlex < 0. ) 00269 { 00270 atmdat.HCTAlex = pow(10.,atmdat.HCTAlex); 00271 } 00272 } 00273 00274 else if( nMatch("CHEM",chCard ) && !nMatch( "HO ", chCard ) ) 00275 { 00276 /* turn on Steve Federman's chemistry */ 00277 if( nMatch("FEDE",chCard ) ) 00278 { 00279 if( nMatch( " ON " , chCard ) ) 00280 { 00281 /* This turns on the diffuse cloud chemical rates of 00282 * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/ 00283 co.lgFederman = true; 00284 } 00285 else if( nMatch( " OFF" , chCard ) ) 00286 { 00287 co.lgFederman = false; 00288 } 00289 else 00290 { 00291 /* this is the default when command used - true */ 00292 co.lgFederman = true; 00293 } 00294 } 00295 /* >>chng 06 may 30 --NPA. Turn on non-equilibrium chemistry */ 00296 else if( nMatch(" NON",chCard ) && nMatch( "EQUI", chCard )) 00297 { 00298 00299 /* option to use effective temperature as defined in 00300 * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319 00301 * By default, this is false - changed with set chemistry command */ 00302 00303 co.lgNonEquilChem = true; 00304 00305 /* >>chng 06 jul 21 -- NPA. Option to include non-equilibrium 00306 * effects for neutral reactions with a temperature dependent rate. 00307 * Reasoning is that non-equilibrium chemistry is caused by MHD, and 00308 * if magnetic field is only coupled to ions, then neutrals may not be 00309 * affected. 00310 * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319 00311 * By default, this is true - changed with set chemistry command */ 00312 00313 if( nMatch("NEUT",chCard ) ) 00314 { 00315 if( nMatch( " ON " , chCard ) ) 00316 { 00317 /* This turns on the diffuse cloud chemical rates of 00318 * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/ 00319 co.lgNeutrals = true; 00320 } 00321 else if( nMatch( " OFF" , chCard ) ) 00322 { 00323 co.lgNeutrals = false; 00324 } 00325 else 00326 { 00327 /* this is the default when command used - true */ 00328 co.lgNeutrals = true; 00329 } 00330 } 00331 } 00332 00333 /* turn off proton elimination rates, which are of the form 00334 * 00335 * 00336 * A + BH+ --> AB + H+ or 00337 * AH + B+ --> AB + H+ 00338 * 00339 * the following paper: 00340 * 00341 * >>refer CO chemistry Huntress, W. T., 1977, ApJS, 33, 495 00342 * says reactions of these types are much less likely than 00343 * identical reactions which leave the product AB ionized (AB+), 00344 * leaving an H instead of H+ (this is called H atom elimination 00345 * currently we only have one reaction of this type, it is 00346 * C+ + OH -> CO + H+ */ 00347 else if( nMatch("PROT",chCard ) && nMatch( "ELIM", chCard ) ) 00348 { 00349 if( nMatch( " ON " , chCard ) ) 00350 { 00351 /* This turns on the diffuse cloud chemical rates of 00352 * >>refer CO chemistry Zsargo, J. & Federman, S. R. 2003, ApJ, 589, 319*/ 00353 co.lgProtElim = true; 00354 } 00355 else if( nMatch( " OFF" , chCard ) ) 00356 { 00357 co.lgProtElim = false; 00358 } 00359 else 00360 { 00361 /* this is the default when command used - true */ 00362 co.lgProtElim = true; 00363 } 00364 } 00365 00366 else 00367 { 00368 /* should not have happened ... */ 00369 fprintf( ioQQQ, " There should have been an option on this SET CHEMISTRY command.\n" ); 00370 fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" ); 00371 cdEXIT(EXIT_FAILURE); 00372 } 00373 } 00374 00375 /* set collision strength averaging ON / OFF */ 00376 else if( nMatch("COLL",chCard) && nMatch("STRE",chCard) && nMatch("AVER",chCard) ) 00377 { 00378 if( nMatch(" OFF",chCard) ) 00379 { 00380 iso.lgCollStrenThermAver = false; 00381 } 00382 else 00383 { 00384 /* this is default behavior of this command */ 00385 iso.lgCollStrenThermAver = true; 00386 } 00387 } 00388 00389 else if( nMatch("COVE",chCard) ) 00390 { 00391 iterations.lgConverge_set = true; 00392 /* set coverage - limit number of iterations and zones */ 00393 if( nMatch("FAST",chCard) ) 00394 { 00395 iterations.lim_zone = 1; 00396 iterations.lim_iter = 0; 00397 } 00398 else 00399 { 00400 iterations.lim_zone = 10; 00401 iterations.lim_iter = 1; 00402 } 00403 } 00404 00405 else if( nMatch("CSUP",chCard) ) 00406 { 00407 /* force secondary ionization rate, log entered */ 00408 i = 5; 00409 secondaries.SetCsupra = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00410 secondaries.lgCSetOn = true; 00411 if( lgEOL ) 00412 { 00413 NoNumb(chCard); 00414 } 00415 secondaries.SetCsupra = (realnum)pow((realnum)10.f,secondaries.SetCsupra); 00416 } 00417 00418 else if( nMatch(" D/H",chCard) ) 00419 { 00420 /* set deuterium abundance, D to H ratio */ 00421 i = 5; 00422 hydro.D2H_ratio = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00423 if( hydro.D2H_ratio <= 0. || nMatch( " LOG", chCard ) ) 00424 { 00425 hydro.D2H_ratio = pow(10.,hydro.D2H_ratio); 00426 } 00427 if( lgEOL ) 00428 { 00429 NoNumb(chCard); 00430 } 00431 } 00432 00433 else if( nMatch( " HO " , chCard) && nMatch( "CHAR" , chCard) ) 00434 { 00435 /* which solver will include H + O charge transfer? done in 00436 * chemistry by default */ 00437 if( nMatch( "CHEM" , chCard ) ) 00438 { 00439 /* this is the default */ 00440 ionbal.lgHO_ct_chem = true; 00441 } 00442 else if( nMatch( "IONI" , chCard ) ) 00443 { 00444 /* do it in the ionization solver */ 00445 ionbal.lgHO_ct_chem = false; 00446 } 00447 else 00448 { 00449 fprintf( ioQQQ, " I did not recognize a subkey on this SET OH CHARge transfer line.\n" ); 00450 fprintf( ioQQQ, " Valid keys are CHEMistry and IONIzation.\n" ); 00451 cdEXIT(EXIT_FAILURE); 00452 } 00453 } 00454 00455 else if( nMatch("12C1",chCard) ) 00456 { 00457 /* set the 12C to 13C abundance ratio - default is 30 */ 00458 i = 5; 00459 /* first two numbers on the line are 12 and 13 - we don't want them */ 00460 co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00461 co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00462 00463 /* now we can get the ratio */ 00464 co.C12_C13_isotope_ratio = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00465 if( lgEOL ) 00466 NoNumb(chCard); 00467 00468 if( co.C12_C13_isotope_ratio <= 0. || nMatch( " LOG", chCard ) ) 00469 co.C12_C13_isotope_ratio = (realnum)pow((realnum)10.f,co.C12_C13_isotope_ratio); 00470 } 00471 00472 /* set dynamics ... */ 00473 else if( nMatch("DYNA",chCard) ) 00474 { 00475 /* set dynamics advection length */ 00476 if( nMatch("ADVE",chCard) && nMatch("LENG",chCard) ) 00477 { 00478 /* <0 => relative fraction of length, + val in cm */ 00479 i = 5; 00480 dynamics.AdvecLengthInit = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00481 if( lgEOL ) 00482 NoNumb( chCard ); 00483 /* if fraction is present, then number was linear fraction, if not present 00484 * then physical length in cm, log10 */ 00485 if( nMatch("FRAC",chCard) ) 00486 { 00487 /* we scanned in the number, if it is a negative then it is the log of the fraction */ 00488 if( dynamics.AdvecLengthInit <= 0. ) 00489 dynamics.AdvecLengthInit = pow(10.,dynamics.AdvecLengthInit); 00490 00491 /* make neg sign as flag for this in dynamics.c */ 00492 dynamics.AdvecLengthInit *= -1.; 00493 } 00494 else 00495 { 00496 /* fraction did not occur, the number is the log of the length in cm - 00497 * convert to linear cm */ 00498 dynamics.AdvecLengthInit = pow(10.,dynamics.AdvecLengthInit); 00499 } 00500 } 00501 else if( nMatch("PRES",chCard) && nMatch("MODE",chCard) ) 00502 { 00503 dynamics.lgSetPresMode = true; 00504 if( nMatch("SUBS",chCard) ) 00505 { 00506 /* subsonic */ 00507 strcpy( dynamics.chPresMode , "subsonic" ); 00508 } 00509 else if( nMatch("SUPE",chCard) ) 00510 { 00511 /* supersonic */ 00512 strcpy( dynamics.chPresMode , "supersonic" ); 00513 } 00514 else if( nMatch("STRO",chCard) ) 00515 { 00516 /* strong d */ 00517 strcpy( dynamics.chPresMode , "strongd" ); 00518 } 00519 else if( nMatch("ORIG",chCard) ) 00520 { 00521 /* original */ 00522 strcpy( dynamics.chPresMode , "original" ); 00523 } 00524 } 00525 else if( nMatch("ANTI",chCard) && nMatch("DEPT",chCard) ) 00526 { 00527 dynamics.lgSetPresMode = true; 00528 strcpy( dynamics.chPresMode , "antishock" ); 00529 /* shock depth */ 00530 i = 5; 00531 /* get log of shock depth in cm */ 00532 dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00533 if( lgEOL ) 00534 NoNumb( chCard ); 00535 dynamics.ShockDepth = pow( 10., dynamics.ShockDepth ); 00536 } 00537 else if( nMatch("ANTI",chCard) && nMatch("MACH",chCard) ) 00538 { 00539 dynamics.lgSetPresMode = true; 00540 strcpy( dynamics.chPresMode , "antishock-by-mach" ); 00541 /* Mach number */ 00542 i = 5; 00543 /* get (isothermal) Mach number where we want antishock to occur */ 00544 dynamics.ShockMach = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00545 if( lgEOL ) 00546 NoNumb( chCard ); 00547 } 00548 else if( nMatch("RELA",chCard) ) 00549 { 00550 /* set how many iterations we will start with, before allowing 00551 * changes. This allows the solution to relax to an equilibrium */ 00552 i = 5; 00553 dynamics.n_initial_relax = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00554 if( lgEOL ) 00555 NoNumb( chCard ); 00556 else if( dynamics.n_initial_relax < 2 ) 00557 { 00558 fprintf(ioQQQ," First iteration to relax dynamics must be > 1." 00559 "It was %li. Sorry.\n", 00560 dynamics.n_initial_relax ); 00561 cdEXIT(EXIT_FAILURE); 00562 } 00563 } 00564 else if( nMatch("SHOC",chCard) && nMatch("DEPT",chCard) ) 00565 { 00566 dynamics.lgSetPresMode = true; 00567 strcpy( dynamics.chPresMode , "shock" ); 00568 /* shock depth */ 00569 i = 5; 00570 /* get log of shock depth in cm */ 00571 dynamics.ShockDepth = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00572 if( lgEOL ) 00573 NoNumb( chCard ); 00574 dynamics.ShockDepth = pow( 10., dynamics.ShockDepth ); 00575 } 00576 else 00577 { 00578 /* should not have happened ... */ 00579 fprintf( ioQQQ, " There should have been an option on this SET DYNAMICS command.\n" ); 00580 fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" ); 00581 cdEXIT(EXIT_FAILURE); 00582 } 00583 } 00584 00585 else if( nMatch("DIDZ",chCard) ) 00586 { 00587 /* set parameter to do with choice of dr; 00588 * par is the largest optical depth to allow in the zone. 00589 * >>chng 96 jan 08 had been two numbers - dtau1 removed */ 00590 i = 5; 00591 radius.drChange = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00592 if( radius.drChange <= 0. ) 00593 { 00594 radius.drChange = (realnum)pow((realnum)10.f,radius.drChange); 00595 } 00596 if( lgEOL ) 00597 { 00598 NoNumb(chCard); 00599 } 00600 } 00601 00602 /* something to do with electron density */ 00603 else if( nMatch("EDEN",chCard) ) 00604 { 00605 /* >>chng 02 apr 20, from ERROR to TOLERANCE to be parallel to set temp equivalent, 00606 * >>chng 04 jun 03, but also accept error as keyword */ 00607 if( nMatch("CONV",chCard) || nMatch("ERRO",chCard)) 00608 { 00609 /* keyword is eden convergence 00610 * set tolerance in eden match */ 00611 i = 5; 00612 conv.EdenErrorAllowed = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00613 if( lgEOL ) 00614 { 00615 NoNumb(chCard); 00616 } 00617 00618 if( conv.EdenErrorAllowed < 0. ) 00619 { 00620 conv.EdenErrorAllowed = pow(10.,conv.EdenErrorAllowed); 00621 } 00622 } 00623 00624 else if( nMatch("SOLV",chCard) ) 00625 { 00626 /* the electron density solver */ 00627 if( nMatch("NEW",chCard) ) 00628 { 00629 /* new method */ 00630 strcpy( conv.chSolverEden , "new" ); 00631 } 00632 else 00633 { 00634 /* default method */ 00635 strcpy( conv.chSolverEden , "simple" ); 00636 } 00637 } 00638 else 00639 { 00640 /* no keyword, assume log of electron density */ 00641 i = 5; 00642 /* set the electron density */ 00643 dense.EdenSet = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00644 /* warn that this model is meaningless */ 00645 phycon.lgPhysOK = false; 00646 } 00647 } 00648 00649 else if( nMatch("FINE",chCard ) && nMatch("CONT",chCard ) ) 00650 { 00651 /* set fine continuum resolution - an element name, used to get 00652 * thermal width, and how many resolution elements to use to resolve 00653 * a line of this element at 1e4 K 00654 * first get an element name, nelem is atomic number on C scale 00655 * default is iron */ 00656 if( (rfield.fine_opac_nelem = GetElem(chCard))<0 ) 00657 { 00658 fprintf( ioQQQ, " An element name must appear on this line\n Sorry.\n" ); 00659 cdEXIT(EXIT_FAILURE); 00660 } 00661 i = 5; 00662 /* set the number of resolution elements in HWHM at 1e4 K for turbulent 00663 * velocity field, default is one element */ 00664 rfield.fine_opac_nresolv = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00665 if( rfield.fine_opac_nresolv< 1 ) 00666 { 00667 fprintf( ioQQQ, " The number of resolution elements within FWHM of line must appear\n Sorry.\n" ); 00668 cdEXIT(EXIT_FAILURE); 00669 } 00670 } 00671 00672 /* set grain command - but not set H2 grain command */ 00673 else if( nMatch("GRAI",chCard ) && !nMatch(" H2 ",chCard) ) 00674 { 00675 if( nMatch("HEAT",chCard ) ) 00676 { 00677 /* scale factor to change grain heating as per Allers et al. */ 00678 i = 5; 00679 gv.GrainHeatScaleFactor = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00680 /* warn that this model is not the best we can do */ 00681 phycon.lgPhysOK = false; 00682 if( lgEOL ) 00683 { 00684 NoNumb(chCard); 00685 } 00686 } 00687 else 00688 { 00689 fprintf( ioQQQ, " A keyword must appear on the SET GRAIN line - options are HEAT \n Sorry.\n" ); 00690 cdEXIT(EXIT_FAILURE); 00691 } 00692 } 00693 00694 /* these are the "set Leiden hack" commands, used to turn off physics and 00695 * sometimes replace with simple approximation */ 00696 else if( nMatch("LEID",chCard ) && nMatch("HACK",chCard ) ) 00697 { 00698 if( nMatch( "H2* " , chCard ) && nMatch( " OFF" , chCard ) ) 00699 { 00700 /* turn off reactions with H2* in the chemistry network */ 00701 hmi.lgLeiden_Keep_ipMH2s = false; 00702 /* warn that this model is not the best we can do */ 00703 phycon.lgPhysOK = false; 00704 } 00705 else if( nMatch( "CR " , chCard ) && nMatch( " OFF" , chCard ) ) 00706 { 00707 /* the CR Leiden hack option - turn off CR excitations of H2 */ 00708 hmi.lgLeidenCRHack = false; 00709 00710 } 00711 else if( nMatch("RATE",chCard ) && nMatch("UMIS",chCard )) 00712 { 00713 /* This command will use the rates given in the UMIST database, It 00714 * will set to zero many reactions in hmole_step.c that are not 00715 * in UMIST */ 00716 co.lgUMISTrates = false; 00717 } 00718 } 00719 00720 /* set H2 ... */ 00721 else if( nMatch(" H2 ",chCard) ) 00722 { 00723 i = 5; 00724 ip = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00725 if( ip != 2 ) 00726 { 00727 fprintf( ioQQQ, " The first number on this line must be the 2 in H2\n Sorry.\n" ); 00728 cdEXIT(EXIT_FAILURE); 00729 } 00730 00731 /* SET H2 SMALL MODEL or SOLOMON - which approximation to Solomon process */ 00732 if( nMatch("SOLO",chCard) || ( nMatch("SMAL",chCard)&&nMatch("MODE",chCard) ) ) 00733 { 00734 if( nMatch("SOLO",chCard) ) 00735 { 00736 /* warn that Solomon will not be used forever */ 00737 fprintf(ioQQQ,"PROBLEM - *set H2 Solomon* has been changed to *set H2 small model*."\ 00738 " This is OK for now but it may not work in a future version.\n"); 00739 } 00740 if( nMatch("TH85", chCard ) ) 00741 { 00742 /* rate from is eqn A8 of Tielens and Hollenbach 85a, */ 00743 hmi.chH2_small_model_type = 'T'; 00744 } 00745 else if( nMatch( " BHT" , chCard ) ) 00746 { 00747 /* the improved H2 formalism given by 00748 *>>refer H2 dissoc Burton, M.G., Hollenbach, D.J. & Tielens, A.G.G.M 00749 >>refcon 1990, ApJ, 365, 620 */ 00750 hmi.chH2_small_model_type = 'H'; 00751 } 00752 else if( nMatch( "BD96" , chCard ) ) 00753 { 00754 /* this rate is equation 23 of 00755 *>>refer H2 dissoc Bertoldi, F., & Draine, B.T., 1996, 458, 222 */ 00756 /* this is the default */ 00757 hmi.chH2_small_model_type = 'B'; 00758 } 00759 else if( nMatch( "ELWE" , chCard ) ) 00760 { 00761 /* this rate is equation 23 of 00762 *>>refer H2 dissoc Elwert et al., in preparation */ 00763 /* this is the default */ 00764 hmi.chH2_small_model_type = 'E'; 00765 } 00766 else 00767 { 00768 fprintf( ioQQQ, " One of the keywords TH85, _BHT, BD96 or ELWErt must appear.\n Sorry.\n" ); 00769 cdEXIT(EXIT_FAILURE); 00770 } 00771 } 00772 00773 /* series of commands that deal with grains */ 00774 /* which approximation to grain formation pumping */ 00775 if( nMatch("GRAI",chCard ) && nMatch("FORM",chCard ) && nMatch("PUMP",chCard ) ) 00776 { 00777 if( nMatch( "DB96" , chCard ) ) 00778 { 00779 /* Draine & Bertoldi 1996 */ 00780 hmi.chGrainFormPump = 'D'; 00781 } 00782 else if( nMatch( "TAKA" , chCard ) ) 00783 { 00784 /* Takahashi 2001 */ 00785 hmi.chGrainFormPump = 'T'; 00786 } 00787 else if( nMatch( "THER" , chCard ) ) 00788 { 00789 /* thermal distribution, upper right column of page 239 of 00790 *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235 */ 00791 hmi.chGrainFormPump = 't'; 00792 } 00793 else 00794 { 00795 fprintf( ioQQQ, " The grain form pump option is wrong.\n Sorry.\n" ); 00796 cdEXIT(EXIT_FAILURE); 00797 } 00798 } 00799 00800 /* which approximation to Jura rate */ 00801 else if( nMatch("JURA",chCard) ) 00802 { 00803 if( nMatch("TH85", chCard ) ) 00804 { 00805 /* rate from is eqn A8 of Tielens and Hollenbach 85a*/ 00806 hmi.chJura = 'T'; 00807 } 00808 else if( nMatch( "CT02" , chCard ) ) 00809 { 00810 /* this rate is equation Cazeux & Tielens */ 00811 hmi.chJura = 'C'; 00812 } 00813 else if( nMatch( "SN99" , chCard ) ) 00814 { 00815 /* this rate is from Sternberg & Neufeld 99 */ 00816 hmi.chJura = 'S'; 00817 } 00818 else if( nMatch( "RATE" , chCard ) ) 00819 { 00820 /* set H2 rate - enters log of Jura rate - F for fixed 00821 * no dependence on grain properties 00822 * had been C, a bug since triggered Cazeux & Tielens 00823 * >>chng 07 jan 10, bug caught by Robin Williams & Fixed by Nick Abel */ 00824 hmi.chJura = 'F'; 00825 hmi.rate_h2_form_grains_set = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) ); 00826 if( lgEOL ) 00827 { 00828 /* no number on the line so use Jura's value, 3e-17 00829 * >>refer H2 Jura Jura, M. 1975, ApJ, 197, 575 */ 00830 hmi.rate_h2_form_grains_set = 3e-17; 00831 } 00832 } 00833 else if( nMatch( "SCAL" , chCard ) ) 00834 { 00835 /* this is a scale factor to multiply the Jura rate */ 00836 hmi.ScaleJura = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00837 /* log or negative number means log was entered */ 00838 if( nMatch( " LOG" , chCard ) || hmi.ScaleJura<0. ) 00839 { 00840 hmi.ScaleJura = (realnum)pow( 10., (double)hmi.ScaleJura ); 00841 } 00842 if( lgEOL ) 00843 { 00844 NoNumb( chCard ); 00845 } 00846 00847 /* option to vary scale factor */ 00848 if( optimize.lgVarOn ) 00849 { 00850 optimize.nvarxt[optimize.nparm] = 1; 00851 strcpy( optimize.chVarFmt[optimize.nparm], "SET H2 JURA SCALE %f" ); 00852 00853 /* pointer to where to write */ 00854 optimize.nvfpnt[optimize.nparm] = input.nRead; 00855 00856 /* log of Jura rate scale factor will be parameter */ 00857 optimize.vparm[0][optimize.nparm] = (realnum)log10(hmi.ScaleJura); 00858 optimize.vincr[optimize.nparm] = 0.3f; 00859 00860 ++optimize.nparm; 00861 } 00862 } 00863 else 00864 { 00865 fprintf( ioQQQ, " The Jura rate option is wrong.\n Sorry.\n" ); 00866 cdEXIT(EXIT_FAILURE); 00867 } 00868 } 00869 00870 /* what temperature to use for binding energy, Tad in Le Bourlot, J., 2000, A&A, 360, 656-662 */ 00871 else if( nMatch(" TAD",chCard) ) 00872 { 00873 hmi.Tad = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00874 if( lgEOL ) 00875 NoNumb(chCard); 00876 /* log if <= 10. unless linear key appears too */ 00877 if( hmi.Tad <=10. && !nMatch("LINE",chCard) ) 00878 hmi.Tad = (realnum)pow((realnum)10.f,hmi.Tad); 00879 } 00880 00881 else if( nMatch("FRAC",chCard ) ) 00882 { 00883 /* this is special option to force H2 abundance to value for testing 00884 * this factor will multiply the hydrogen density to become the H2 density 00885 * no attempt to conserve particles, or do the rest of the molecular equilibrium 00886 * set consistently is made */ 00887 hmi.H2_frac_abund_set = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00888 if( lgEOL ) 00889 NoNumb( chCard ); 00890 00891 /* a number <= 0 is the log of the ratio */ 00892 if( hmi.H2_frac_abund_set <= 0. ) 00893 hmi.H2_frac_abund_set = pow(10., hmi.H2_frac_abund_set); 00894 /* don't let it exceed 0.5 */ 00895 /* >>chng 03 jul 19, from 0.5 to 0.4999, do not want atomic density exactly zero */ 00896 hmi.H2_frac_abund_set = MIN2(0.49999 , hmi.H2_frac_abund_set ); 00897 } 00898 } 00899 00900 /* this is a scale factor that changes the n(H0)*1.7e-4 that is added to the 00901 * electron density to account for collisions with atomic H. it is an order of 00902 * magnitude guess, so this command provides ability to see whether it affects results */ 00903 else if( nMatch("HCOR",chCard) ) 00904 { 00905 i = 5; 00906 dense.HCorrFac = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00907 if( lgEOL ) 00908 NoNumb( chCard ); 00909 } 00910 00911 else if( nMatch(" PAH",chCard) ) 00912 { 00913 /* set one of several possible PAH abundance distribution functions */ 00914 if( nMatch(" H0 " , chCard ) ) 00915 { 00916 /* the default, abundance depends on H0 / Htot */ 00917 strcpy( gv.chPAH_abundance_fcn , "H0" ); 00918 } 00919 else if( nMatch("CONS" , chCard ) ) 00920 { 00921 /* constant PAH abundance */ 00922 strcpy( gv.chPAH_abundance_fcn , "CON" ); 00923 } 00924 00925 else if(nMatch("BAKE",chCard ) ) 00926 { 00927 /* turn on simple PAH heating from Bakes & Tielens - this is very approximate */ 00928 /*>>>refer PAH heating Bakes, E.L.O., & Tielens, A.G.G.M. 1994, ApJ, 427, 822 */ 00929 gv.lgBakesPAH_heat = true; 00930 /* warn that this model is not the best we can do */ 00931 phycon.lgPhysOK = false; 00932 } 00933 else 00934 { 00935 fprintf( ioQQQ, " one of the keywords H0, BAKES, or CONStant must appear.\n Sorry." ); 00936 cdEXIT(EXIT_FAILURE); 00937 } 00938 } 00939 00940 else if( nMatch("PRES",chCard) && nMatch("IONI",chCard) ) 00941 { 00942 /* set limit to number of calls from pressure to ionize solver, 00943 * this set limit to conv.nPres2Ioniz */ 00944 i = 5; 00945 conv.limPres2Ioniz = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00946 if( lgEOL ) 00947 { 00948 NoNumb(chCard); 00949 } 00950 else if( conv.limPres2Ioniz <= 0 ) 00951 { 00952 fprintf( ioQQQ, " The limit must be greater than zero.\n Sorry." ); 00953 cdEXIT(EXIT_FAILURE); 00954 } 00955 } 00956 /* something to do with pressure */ 00957 else if( nMatch("PRES",chCard) ) 00958 { 00959 /* tolerance on pressure convergence */ 00960 if( nMatch("CONV",chCard) ) 00961 { 00962 /* keyword is tolerance 00963 * set tolerance or relative error in pressure match */ 00964 i = 5; 00965 conv.PressureErrorAllowed = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00966 if( lgEOL ) 00967 NoNumb(chCard); 00968 00969 if( conv.PressureErrorAllowed < 0. ) 00970 conv.PressureErrorAllowed = (realnum)pow((realnum)10.f,conv.PressureErrorAllowed); 00971 } 00972 00973 else 00974 { 00975 /* >>chng 04 mar 02, printout had been wrong, saying TOLErange 00976 * rather than CONVergence. Nick Abel */ 00977 fprintf( ioQQQ, " I didn\'t recognize a key on this SET PRESSURE line.\n" ); 00978 fprintf( ioQQQ, " The ones I know about are: CONVergence.\n" ); 00979 cdEXIT(EXIT_FAILURE); 00980 } 00981 } 00982 else if( nMatch("RECO",chCard) && nMatch("MBIN",chCard) ) 00983 { 00984 /* dielectronic recombination */ 00985 if( nMatch("DIEL",chCard) && nMatch("ECTR",chCard) ) 00986 { 00987 /* change various factors for dielectronic recombination */ 00988 if( nMatch("KLUD",chCard) ) 00989 { 00990 if( nMatch("BADN",chCard) ) 00991 { 00992 /* set true to use Badnell mean in place of existing hacks */ 00993 ionbal.lg_use_DR_Badnell_rate_coef_mean_ion = true; 00994 } 00995 else 00996 { 00997 int j; 00998 /* option to turn off guess of dielectronic rec coefficient for 3rd row elements */ 00999 i = 3; 01000 /* this is first call, no number, lgEOL true but zero returned, the intended effect*/ 01001 ionbal.GuessDiel[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01002 /* >>chng 06 jan 30 also turn off guess if GuessDiel set zero */ 01003 if( ionbal.GuessDiel[0]<=0. || nMatch(" OFF",chCard) ) 01004 ionbal.lg_guess_coef = false; 01005 01006 for( j=1; j<4; ++j ) 01007 { 01008 ionbal.GuessDiel[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01009 /* here lgEOL means to use previous values */ 01010 if( lgEOL ) 01011 ionbal.GuessDiel[j] = ionbal.GuessDiel[j-1]; 01012 } 01013 } 01014 /* option to add log normal noise */ 01015 if( nMatch("NOISE" , chCard ) ) 01016 { 01017 i = 3; 01018 ionbal.guess_noise = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01019 if( lgEOL ) 01020 ionbal.guess_noise = 2.; 01021 /* Seed the random-number generator with current time so that 01022 * the numbers will be different every time we run. */ 01023 init_genrand( (unsigned)time( NULL ) ); 01024 } 01025 else 01026 { 01027 ionbal.guess_noise = 0.; 01028 } 01029 } 01030 01031 else if( nMatch("BURG" , chCard) ) 01032 { 01033 /* modify suppression of burgess dielectronic recombinations */ 01034 if( nMatch(" ON ",chCard) ) 01035 { 01036 /* leave suppression on - this is default state */ 01037 ionbal.lgSupDie[0] = true; 01038 } 01039 01040 else if( nMatch(" OFF",chCard) ) 01041 { 01042 /* turn suppression off */ 01043 ionbal.lgSupDie[0] = false; 01044 } 01045 01046 else 01047 { 01048 fprintf( ioQQQ, " flag ON or OFF must appear.\n" ); 01049 cdEXIT(EXIT_FAILURE); 01050 } 01051 } 01052 01053 else if( nMatch("NUSS" , chCard) ) 01054 { 01055 /* modify suppression of Nussbaumer and Storey dielectronic recombination */ 01056 if( nMatch(" ON ",chCard) ) 01057 { 01058 /* turn suppression on */ 01059 ionbal.lgSupDie[1] = true; 01060 } 01061 else if( nMatch(" OFF",chCard) ) 01062 { 01063 /* leave suppression off - this is default state */ 01064 ionbal.lgSupDie[1] = false; 01065 } 01066 else 01067 { 01068 fprintf( ioQQQ, " flag ON or OFF must appear.\n" ); 01069 cdEXIT(EXIT_FAILURE); 01070 } 01071 } 01072 01073 else if( nMatch("BADN" , chCard ) ) 01074 { 01075 /* use the new (in early 2006) Badnell DR rates */ 01076 if( nMatch( " OFF" , chCard ) ) 01077 { 01078 ionbal.lgDR_recom_Badnell_use = false; 01079 } 01080 else 01081 { 01082 ionbal.lgDR_recom_Badnell_use = true; 01083 } 01084 01085 /* option to print rates then exit */ 01086 if( nMatch("PRIN" , chCard ) ) 01087 ionbal.lgRecom_Badnell_print = true; 01088 } 01089 else 01090 { 01091 fprintf( ioQQQ, " key KLUDge, BURGess, NUSSbaumer, or BADNell must appear.\n" ); 01092 cdEXIT(EXIT_FAILURE); 01093 } 01094 } 01095 /* radiative recombination */ 01096 else if( nMatch("RADI",chCard) && nMatch("ATIV",chCard) ) 01097 { 01098 if( nMatch("BADN" , chCard ) ) 01099 { 01100 /* use the new (in early 2006) Badnell DR rates */ 01101 if( nMatch( " OFF" , chCard ) ) 01102 { 01103 ionbal.lgRR_recom_Badnell_use = false; 01104 } 01105 else 01106 { 01107 ionbal.lgRR_recom_Badnell_use = true; 01108 } 01109 01110 /* option to print rates then exit */ 01111 if( nMatch("PRIN" , chCard ) ) 01112 ionbal.lgRecom_Badnell_print = true; 01113 } 01114 else 01115 { 01116 fprintf( ioQQQ, " key BADNell must appear.\n" ); 01117 cdEXIT(EXIT_FAILURE); 01118 } 01119 } 01120 else 01121 { 01122 fprintf( ioQQQ, " key RADIATIve or DIELECTRonic must appear on set recombination command.\n" ); 01123 cdEXIT(EXIT_FAILURE); 01124 } 01125 } 01126 01127 else if( nMatch("SPEC",chCard) ) 01128 { 01129 /* "set spectrum" for optional parameters for punch spectrum command */ 01130 if( nMatch("RANG" , chCard ) ) 01131 { 01132 /* energy range option - argument is energy in Rydbergs */ 01133 i = 5; 01134 /* >>chng 05 aug 25, rm pow(10., - was wrong */ 01135 punch.cp_range_min[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01136 punch.cp_range_max[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01137 if( lgEOL ) 01138 { 01139 NoNumb(chCard); 01140 } 01141 /* confirm that wavelengths are in correct order */ 01142 if( punch.cp_range_min[punch.cp_npun] >= punch.cp_range_max[punch.cp_npun] ) 01143 { 01144 fprintf( ioQQQ, " The limits must be in increasing order.\n" ); 01145 cdEXIT(EXIT_FAILURE); 01146 } 01147 } 01148 01149 else if( nMatch("RESO" , chCard ) ) 01150 { 01151 /* resolving power, default is zero, which means leave at native resolution */ 01152 /* >>chng 05 aug 25, following had pow(number), so crashed for high resolving power 01153 * bug caught by Yan Changshou */ 01154 i = 5; 01155 punch.cp_resolving_power[punch.cp_npun] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01156 if( lgEOL ) 01157 { 01158 NoNumb(chCard); 01159 } 01160 } 01161 } 01162 01163 else if( nMatch(" DR ",chCard) ) 01164 { 01165 /* set zone thickness by forcing drmax and drmin */ 01166 i = 5; 01167 /* at this stage linear, but default is log */ 01168 radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01169 if( !nMatch("LINE",chCard) ) 01170 { 01171 /* linear was not on command, so default is log */ 01172 radius.sdrmax = pow( 10. , radius.sdrmax ); 01173 } 01174 if( lgEOL ) 01175 { 01176 NoNumb(chCard); 01177 } 01178 01179 /* NB these being equal are tested in convinittemp to set dr */ 01180 radius.sdrmin = radius.sdrmax; 01181 if( radius.sdrmax < DEPTH_OFFSET*1e4 ) 01182 { 01183 fprintf( ioQQQ, "\n Thicknesses less than about %.0e will NOT give accurate results. If tricking the code\n", 01184 DEPTH_OFFSET*1e4 ); 01185 fprintf( ioQQQ, " into computing emissivities instead of intensities, try to instead use a thickness of unity,\n" ); 01186 fprintf( ioQQQ, " and then multiply (divide) the results by the necessary thickness (product of densities).\n\n" ); 01187 cdEXIT(EXIT_FAILURE); 01188 } 01189 } 01190 01191 else if( nMatch("DRMA",chCard) ) 01192 { 01193 /* set maximum zone thickness */ 01194 i = 5; 01195 radius.sdrmax = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01196 if( lgEOL ) 01197 NoNumb(chCard); 01198 01199 /* log if less than 38 or if log keyword occurs */ 01200 if( radius.sdrmax < log10(BIGFLOAT) || nMatch(" LOG",chCard) ) 01201 radius.sdrmax = pow(10.,radius.sdrmax); 01202 } 01203 01204 else if( nMatch("DRMI",chCard) ) 01205 { 01206 /* option to set minimum zone thickness */ 01207 i = 5; 01208 radius.sdrmin = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01209 if( lgEOL ) 01210 NoNumb(chCard); 01211 01212 /* log if less than 38 or if log keyword occurs */ 01213 if( radius.sdrmin < log10(BIGFLOAT ) || nMatch(" LOG",chCard) ) 01214 radius.sdrmin = pow(10.,radius.sdrmin); 01215 01216 radius.lgSMinON = true; 01217 } 01218 01219 else if( nMatch("FLXF",chCard) ) 01220 { 01221 /* faintest continuum flux to consider */ 01222 i = 5; 01223 rfield.FluxFaint = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01224 if( lgEOL ) 01225 { 01226 NoNumb(chCard); 01227 } 01228 if( rfield.FluxFaint < 0. ) 01229 { 01230 rfield.FluxFaint = (realnum)pow((realnum)10.f,rfield.FluxFaint); 01231 } 01232 } 01233 01234 else if( nMatch("LINE",chCard) && nMatch("PREC",chCard) ) 01235 { 01236 /* set line precision (number of significant figures) 01237 * this affects all aspects of reading and writing lines */ 01238 i = 5; 01239 LineSave.sig_figs = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01240 if( lgEOL ) 01241 { 01242 NoNumb(chCard); 01243 } 01244 else if( LineSave.sig_figs > 6 ) 01245 { 01246 fprintf( ioQQQ, " set line precision currently only works for up to 6 significant figures.\n" ); 01247 cdEXIT(EXIT_FAILURE); 01248 } 01249 /* option to print main line array as a single column */ 01250 /* prt.lgPrtLineArray = false; */ 01251 } 01252 01253 /* >>chng 00 dec 08, command added by Peter van Hoof */ 01254 else if( nMatch("NFNU",chCard) ) 01255 { 01256 /* set nFnu [incident_reflected] [incident_transmitted] [diffuse_inward] [diffuse_outward] 01257 * command for specifying what to include in the nFnu entries in LineSv */ 01258 /* >>chng 01 nov 12, also accept form with space */ 01259 /* "incident reflected" keyword */ 01260 prt.lgSourceReflected = nMatch("NT R",chCard) || nMatch("NT_R",chCard); 01261 /* "incident transmitted" keyword */ 01262 prt.lgSourceTransmitted = nMatch("NT_T",chCard) || nMatch("NT T",chCard); 01263 /* "diffuse inward" keyword */ 01264 prt.lgDiffuseInward = nMatch("SE_I",chCard) || nMatch("SE I",chCard); 01265 /* "diffuse outward" keyword */ 01266 prt.lgDiffuseOutward = nMatch("SE_O",chCard) || nMatch("SE O",chCard); 01267 01268 /* at least one of these needs to be set ! */ 01269 if( ! ( prt.lgSourceReflected || prt.lgSourceTransmitted || 01270 prt.lgDiffuseInward || prt.lgDiffuseOutward ) ) 01271 { 01272 fprintf( ioQQQ, " set nFnu expects one or more of the following keywords:\n" ); 01273 fprintf( ioQQQ, " INCIDENT_REFLECTED, INCIDENT_TRANSMITTED, DIFFUSE_INWARD, DIFFUSE_OUTWARD\n" ); 01274 cdEXIT(EXIT_FAILURE); 01275 } 01276 /* automatically print the nFnu items in the output - it is not necessary to also include 01277 * a print continua command if this is entered */ 01278 prt.lgPrnDiff = true; 01279 } 01280 01281 else if( nMatch("IND2",chCard) ) 01282 { 01283 if( nMatch(" ON ",chCard) ) 01284 { 01285 /* set flag saying to off or on induced two photon processes */ 01286 iso.lgInd2nu_On = true; 01287 } 01288 else if( nMatch(" OFF",chCard) ) 01289 { 01290 iso.lgInd2nu_On = false; 01291 } 01292 else 01293 { 01294 fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" ); 01295 cdEXIT(EXIT_FAILURE); 01296 } 01297 } 01298 01299 else if( nMatch("TEMP",chCard) ) 01300 { 01301 /* set something to do with temperature, currently solver, tolerance, floor */ 01302 if( nMatch("SOLV",chCard) ) 01303 { 01304 /* solver */ 01305 /* the electron density solver */ 01306 if( nMatch("NEW",chCard) ) 01307 { 01308 /* new method */ 01309 strcpy( conv.chSolverTemp , "new" ); 01310 } 01311 else 01312 { 01313 /* default method */ 01314 strcpy( conv.chSolverTemp , "simple" ); 01315 } 01316 } 01317 01318 else if( nMatch("FLOO",chCard) ) 01319 { 01320 i = 5; 01321 StopCalc.TeFloor = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01322 01323 /* linear option */ 01324 if( StopCalc.TeFloor <= 10. && !nMatch("LINE",chCard) ) 01325 { 01326 StopCalc.TeFloor = (realnum)pow(10.,StopCalc.TeFloor); 01327 } 01328 01329 if( lgEOL ) 01330 { 01331 NoNumb(chCard); 01332 } 01333 01334 if( StopCalc.TeFloor < 2.8 ) 01335 { 01336 fprintf( ioQQQ, " TE < 3K, reset to 2.8K.\n" ); 01337 StopCalc.TeFloor = 2.8f; 01338 } 01339 } 01340 01341 else if( nMatch("CONV",chCard) || nMatch("TOLE",chCard) ) 01342 { 01343 /* error tolerance in heating cooling match, number is error/total */ 01344 i = 5; 01345 conv.HeatCoolRelErrorAllowed = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01346 if( lgEOL ) 01347 { 01348 NoNumb(chCard); 01349 } 01350 if( conv.HeatCoolRelErrorAllowed <= 0. ) 01351 { 01352 conv.HeatCoolRelErrorAllowed = (realnum)pow((realnum)10.f,conv.HeatCoolRelErrorAllowed); 01353 } 01354 } 01355 01356 else 01357 { 01358 fprintf( ioQQQ, "\nI did not recognize a keyword on this SET TEPERATURE command.\n%s\n" , chCard); 01359 fprintf( ioQQQ, "The keywords are SOLVer, FLOOr, and CONVergence.\n" ); 01360 cdEXIT(EXIT_FAILURE); 01361 } 01362 } 01363 01364 else if( nMatch("TEST",chCard) ) 01365 { 01366 /* set flag saying to turn on some test - this is in cddefines.h in the global namespace */ 01367 lgTestCodeEnabled = true; 01368 } 01369 01370 else if( nMatch("TRIM",chCard) ) 01371 { 01372 /* set trim upper or lower, for ionization stage trimming 01373 * ion trimming ionization trimming 01374 * in routine TrimStage */ 01375 i = 5; 01376 if( nMatch("UPPE",chCard) ) 01377 { 01378 /* set trim upper - either set value or turn off */ 01379 double save = ionbal.trimhi; 01380 ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 01381 if( lgEOL && nMatch(" OFF",chCard) ) 01382 { 01383 /* turn off upward trimming */ 01384 lgEOL = false; 01385 ionbal.lgTrimhiOn = false; 01386 /* reset high limit to proper value */ 01387 ionbal.trimhi = save; 01388 } 01389 } 01390 01391 else if( nMatch("LOWE",chCard) ) 01392 { 01393 /* set trim lower */ 01394 ionbal.trimlo = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 01395 } 01396 01397 /* turn off ionization stage trimming */ 01398 else if( nMatch("SMAL",chCard) || nMatch(" OFF",chCard) ) 01399 { 01400 /* set small limits to both upper and lower limit*/ 01401 ionbal.trimlo = SMALLFLOAT; 01402 ionbal.trimhi = SMALLFLOAT; 01403 lgEOL = false; 01404 } 01405 01406 else 01407 { 01408 /* set trim upper */ 01409 ionbal.trimhi = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 01410 01411 /* set trim lower to same number */ 01412 ionbal.trimlo = ionbal.trimhi; 01413 } 01414 01415 if( lgEOL ) 01416 { 01417 NoNumb(chCard); 01418 } 01419 01420 if( ionbal.trimlo >= 1. || ionbal.trimhi >= 1. ) 01421 { 01422 fprintf( ioQQQ, " number must be negative since log\n" ); 01423 cdEXIT(EXIT_FAILURE); 01424 } 01425 } 01426 01427 else if( nMatch("SKIP",chCard) ) 01428 { 01429 /* skip punch command, for punching every n't point */ 01430 i = 5; 01431 punch.ncPunchSkip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01432 if( lgEOL ) 01433 { 01434 NoNumb(chCard); 01435 } 01436 } 01437 01438 else if( nMatch(" UTA",chCard) ) 01439 { 01440 /* set UTA command, to determine which UTA data set to use 01441 * default is to use Gu et al. 2006 data set */ 01442 if( nMatch("GU06" , chCard ) ) 01443 { 01444 if( nMatch(" OFF" , chCard ) ) 01445 { 01446 /* use Behar et al. 2001 */ 01447 ionbal.lgInnerShell_Gu06 = false; 01448 } 01449 else 01450 { 01451 /* the default, use 01452 *>>refer Fe UTA Gu, M.F., Holczer, T., Behar, E. & Kahn, S.M. 01453 *>>refercon 2006, ApJ, 641, 1227 */ 01454 ionbal.lgInnerShell_Gu06 = true; 01455 } 01456 } 01457 else if( nMatch("KISI" , chCard ) ) 01458 { 01459 if( nMatch(" OFF" , chCard ) ) 01460 { 01461 /* the default, do not use it */ 01462 ionbal.lgInnerShell_Kisielius = false; 01463 } 01464 else 01465 { 01466 /* use the data from Kisielius et al. 2003, MNRAS, 344, 696 */ 01467 ionbal.lgInnerShell_Kisielius = true; 01468 } 01469 } 01470 } 01471 01472 else if( nMatch("EAKH",chCard) ) 01473 { 01474 /* set WeakHeatCool, threshold on punch heating and cooling, default 0.05 */ 01475 i = 5; 01476 punch.WeakHeatCool = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01477 01478 if( lgEOL ) 01479 { 01480 NoNumb(chCard); 01481 } 01482 01483 if( punch.WeakHeatCool < 0. ) 01484 { 01485 punch.WeakHeatCool = (realnum)pow((realnum)10.f,punch.WeakHeatCool); 01486 } 01487 } 01488 01489 else if( nMatch("KSHE",chCard) ) 01490 { 01491 /* upper limit to opacities for k-shell ionization */ 01492 i = 5; 01493 continuum.EnergyKshell = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01494 if( lgEOL ) 01495 { 01496 NoNumb(chCard); 01497 } 01498 01499 if( continuum.EnergyKshell == 0. ) 01500 { 01501 /* arg of 0 causes upper limit to energy array */ 01502 continuum.EnergyKshell = rfield.egamry; 01503 } 01504 01505 else if( continuum.EnergyKshell < 194. ) 01506 { 01507 fprintf( ioQQQ, " k-shell energy must be greater than 194 Ryd\n" ); 01508 cdEXIT(EXIT_FAILURE); 01509 } 01510 } 01511 01512 else if( nMatch("NCHR",chCard) ) 01513 { 01514 /* option to set the number of charge states for grain model */ 01515 long ii = 5; 01516 double val = FFmtRead(chCard,&ii,INPUT_LINE_LENGTH,&lgEOL); 01517 01518 if( lgEOL ) 01519 { 01520 NoNumb(chCard); 01521 } 01522 else 01523 { 01524 long nChrg = nint(val); 01525 if( nChrg < 2 || nChrg > NCHS ) 01526 { 01527 fprintf( ioQQQ, " illegal value for number of charge states: %ld\n", nChrg ); 01528 fprintf( ioQQQ, " choose a value between 2 and %d\n", NCHS ); 01529 fprintf( ioQQQ, " or increase NCHS in grainvar.h and recompile\n" ); 01530 cdEXIT(EXIT_FAILURE); 01531 } 01532 else 01533 { 01534 SetNChrgStates(nChrg); 01535 } 01536 } 01537 } 01538 01539 else if( nMatch("NEGO",chCard) ) 01540 { 01541 /* punch negative opacities if they occur, set negopac */ 01542 opac.lgNegOpacIO = true; 01543 } 01544 01545 else if( nMatch("NEND",chCard) ) 01546 { 01547 /* default limit to number of zones to be computed 01548 * only do this if nend is NOT currently left at its default 01549 * nend is set to nEndDflt in routine zero 01550 * this command only has effect if stop zone not entered */ 01551 if( geometry.lgEndDflt ) 01552 { 01553 i = 5; 01554 /* this is default limit to number of zones, change it to this value */ 01555 geometry.nEndDflt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01556 geometry.lgEndDflt = false; 01557 if( lgEOL ) 01558 { 01559 NoNumb(chCard); 01560 } 01561 01562 /* now change all limits, for all iterations, to this value */ 01563 for( i=0; i < iterations.iter_malloc; i++ ) 01564 { 01565 geometry.nend[i] = geometry.nEndDflt; 01566 } 01567 } 01568 } 01569 01570 else if( nMatch("TSQD",chCard) ) 01571 { 01572 /* upper limit for highest density considered in the 01573 * Peimbert-style t^2 section of the printout */ 01574 i = 5; 01575 peimbt.tsqden = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01576 01577 if( lgEOL ) 01578 { 01579 NoNumb(chCard); 01580 } 01581 peimbt.tsqden = (realnum)pow((realnum)10.f,peimbt.tsqden); 01582 } 01583 01584 else if( nMatch("NMAP",chCard) ) 01585 { 01586 /* how many steps in plot or punch of heating-cooling map */ 01587 i = 5; 01588 hcmap.nMapStep = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01589 if( lgEOL ) 01590 { 01591 NoNumb(chCard); 01592 } 01593 } 01594 01595 else if( nMatch("NUME",chCard) && nMatch("DERI",chCard) ) 01596 { 01597 /* this is an option to use numerical derivatives for heating and cooling */ 01598 NumDeriv.lgNumDeriv = true; 01599 } 01600 01601 else if( nMatch("PATH",chCard) ) 01602 { 01603 fprintf( ioQQQ, " The SET PATH command is no longer supported.\n" ); 01604 fprintf( ioQQQ, " Please set the correct path in the file path.h.\n" ); 01605 fprintf( ioQQQ, " Or set the environment variable CLOUDY_DATA_PATH.\n Sorry.\n" ); 01606 cdEXIT(EXIT_FAILURE); 01607 } 01608 01609 else if( nMatch("PHFI",chCard) ) 01610 { 01611 /* which version of PHFIT to use, 1995 or 1996 */ 01612 i = 5; 01613 ip = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01614 if( lgEOL ) 01615 { 01616 NoNumb(chCard); 01617 } 01618 01619 if( ip == 1995 ) 01620 { 01621 /* option to go back to old results, pre op */ 01622 t_ADfA::Inst().set_version( PHFIT95 ); 01623 } 01624 else if( ip == 1996 ) 01625 { 01626 /* default is to use newer results, including opacity project */ 01627 t_ADfA::Inst().set_version( PHFIT96 ); 01628 } 01629 else 01630 { 01631 fprintf( ioQQQ, " Two possible values are 1995 and 1996.\n" ); 01632 cdEXIT(EXIT_FAILURE); 01633 } 01634 } 01635 01636 /* set punch xxx command */ 01637 else if( nMatch("PUNC",chCard) ) 01638 { 01639 if( nMatch("HASH",chCard) ) 01640 { 01641 /* to use the hash option there must be double quotes on line - were there? */ 01642 if( !lgQuotesFound ) 01643 { 01644 fprintf( ioQQQ, " I didn\'t recognize a key on this SET PUNCH HASH line.\n" ); 01645 cdEXIT(EXIT_FAILURE); 01646 } 01647 /* specify the terminator between punch output sets - normally a set of hash marks */ 01648 /* 01649 * get any string within double quotes, and return it as 01650 * null terminated string 01651 * also sets name in OrgCard and chCard to blanks so that 01652 * do not trigger off it later */ 01653 if( strcmp( chString_quotes_lowercase , "return" ) == 0 ) 01654 { 01655 /* special case - return becomes new line */ 01656 sprintf(punch.chHashString , "\n" ); 01657 } 01658 else if( strcmp( chString_quotes_lowercase , "time" ) == 0 ) 01659 { 01660 /* special case where output time between iterations */ 01661 sprintf(punch.chHashString , "TIME_DEP" ); 01662 } 01663 else 01664 { 01665 /* usual case, simply copy what is in quotes */ 01666 sprintf(punch.chHashString , chString_quotes_lowercase ); 01667 } 01668 } 01669 01670 else if( nMatch("WIDT",chCard) ) 01671 { 01672 /* set punch line width for contrast in continuum plots */ 01673 i = 5; 01674 /* units are km/sec */ 01675 punch.PunchLWidth = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01676 01677 /* lgEOL if no number hit, could have been c */ 01678 if( lgEOL ) 01679 { 01680 if( nMatch(" C " , chCard ) ) 01681 { 01682 /* no number but _c_ entered, so enter speed of light in km/s */ 01683 punch.PunchLWidth = (realnum)(SPEEDLIGHT/1e5); 01684 } 01685 else 01686 { 01687 NoNumb(chCard); 01688 } 01689 } 01690 01691 if( punch.PunchLWidth <= 0. ) 01692 { 01693 fprintf( ioQQQ, " line width must be greater than zero.\n" ); 01694 cdEXIT(EXIT_FAILURE); 01695 } 01696 01697 /* >>chng 06 dec 06, value of PunchLWidth is km/s but SPEEDLIGHT 01698 * is cm/s - add 1e5 before compare with SPEEDLIGHT */ 01699 else if( punch.PunchLWidth*0.9999e5 > SPEEDLIGHT ) 01700 { 01701 fprintf( ioQQQ, " line width must be entered in km/s and be less than c.\n" ); 01702 cdEXIT(EXIT_FAILURE); 01703 } 01704 /* this is the factor that is used to add the lines to the continuum, 01705 * 15 converts to cm/s */ 01706 punch.PunchLWidth = (realnum)(SPEEDLIGHT/(punch.PunchLWidth*1e5)); 01707 01708 } 01709 01710 else if( nMatch("PREF",chCard) ) 01711 { 01712 /* specify a prefix before all punch filenames */ 01713 /* 01714 * get any string within double quotes, and return it as 01715 * null terminated string 01716 * also sets name in OrgCard and chCard to blanks so that 01717 * do not trigger off it later */ 01718 01719 strcpy( punch.chFilenamePrefix , chString_quotes_lowercase ); 01720 } 01721 01722 else if( nMatch("FLUS",chCard) ) 01723 { 01724 /* flus the output after every iteration */ 01725 punch.lgFLUSH = true; 01726 } 01727 01728 else 01729 { 01730 fprintf( ioQQQ, " There should have been an option on this command.\n" ); 01731 fprintf( ioQQQ, " Valid options are HASH and PREFIX.\n" ); 01732 cdEXIT(EXIT_FAILURE); 01733 } 01734 } 01735 01736 /* set continuum .... options */ 01737 else if( nMatch("CONT" , chCard ) ) 01738 { 01739 if( nMatch("RESO" , chCard ) ) 01740 { 01741 /* set resolution, get factor that will multiply continuum resolution that 01742 * is contained in the file continuum_mesh.ini */ 01743 i = 5; 01744 continuum.ResolutionScaleFactor = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01745 if( lgEOL ) 01746 { 01747 NoNumb(chCard); 01748 } 01749 /* negative numbers were logs */ 01750 if( continuum.ResolutionScaleFactor <= 0. ) 01751 continuum.ResolutionScaleFactor = pow(10.,continuum.ResolutionScaleFactor); 01752 } 01753 01754 else if( nMatch("SHIE" , chCard ) ) 01755 { 01756 /* set continuum shielding function */ 01757 /* these are all possible values of rt.nLineContShield, 01758 * first is default, these are set with set continuum shielding */ 01759 /*#define LINE_CONT_SHIELD_PESC 1*/ 01760 /*#define LINE_CONT_SHIELD_FEDERMAN 2*/ 01761 /*#define LINE_CONT_SHIELD_FERLAND 3*/ 01762 if( nMatch("PESC" , chCard ) ) 01763 { 01764 /* this uses an inward looking escape probability */ 01765 rt.nLineContShield = LINE_CONT_SHIELD_PESC; 01766 } 01767 else if( nMatch("FEDE" , chCard ) ) 01768 { 01769 /* set continuum shielding Federman, 01770 * this is the default, and uses the appendix of 01771 * >>refer RT continuum shielding Federman, S.R., Glassgold, A.E., & 01772 * >>refercon Kwan, J. 1979, ApJ, 227, 466*/ 01773 rt.nLineContShield = LINE_CONT_SHIELD_FEDERMAN; 01774 } 01775 else if( nMatch("FERL" , chCard ) ) 01776 { 01777 rt.nLineContShield = LINE_CONT_SHIELD_FERLAND; 01778 } 01779 else if( nMatch("NONE" , chCard ) ) 01780 { 01781 /* turn off continuum shielding */ 01782 rt.nLineContShield = 0; 01783 } 01784 else 01785 { 01786 fprintf( ioQQQ, " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n" ); 01787 fprintf( ioQQQ, " The ones I know about are: PESC, FEDErman, & FERLand.\n" ); 01788 cdEXIT(EXIT_FAILURE); 01789 } 01790 } 01791 01792 else 01793 { 01794 fprintf( ioQQQ, " I didn\'t recognize a key on this SET CONTINUUM line.\n" ); 01795 fprintf( ioQQQ, " The ones I know about are: RESOlution and SHIEld.\n" ); 01796 cdEXIT(EXIT_FAILURE); 01797 } 01798 } 01799 01800 else 01801 { 01802 fprintf( ioQQQ, " I didn\'t recognize a key on this SET command.\n" ); 01803 cdEXIT(EXIT_FAILURE); 01804 } 01805 01806 return; 01807 }