00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "geometry.h"
00007 #include "input.h"
00008 #include "prt.h"
00009 #include "atomfeii.h"
00010 #include "mole.h"
00011 #include "rt.h"
00012 #include "phycon.h"
00013 #include "optimize.h"
00014 #include "hcmap.h"
00015 #include "hmi.h"
00016 #include "dense.h"
00017 #include "h2.h"
00018 #include "co.h"
00019 #include "iterations.h"
00020 #include "conv.h"
00021 #include "secondaries.h"
00022 #include "rfield.h"
00023 #include "ionbal.h"
00024 #include "numderiv.h"
00025 #include "dynamics.h"
00026 #include "iso.h"
00027 #include "mole.h"
00028 #include "predcont.h"
00029 #include "save.h"
00030 #include "stopcalc.h"
00031 #include "opacity.h"
00032 #include "hydrogenic.h"
00033 #include "peimbt.h"
00034 #include "radius.h"
00035 #include "atmdat.h"
00036 #include "continuum.h"
00037 #include "grains.h"
00038 #include "grainvar.h"
00039 #include "parser.h"
00040 #include "lines.h"
00041 #include "monitor_results.h"
00042 #include "thirdparty.h"
00043
00044 void ParseSet(Parser &p)
00045 {
00046 long int ip;
00047 char chString_quotes_lowercase[INPUT_LINE_LENGTH];
00048 DEBUG_ENTRY( "ParseSet()" );
00049
00050
00051
00052
00053 bool lgQuotesFound = true;
00054 if (p.GetQuote(chString_quotes_lowercase, false))
00055 lgQuotesFound = false;
00056
00057
00058 if (p.nMatch("ASSE") && p.nMatch("ABOR"))
00059 {
00060
00061
00062 cpu.i().setAssertAbort(true);
00063 }
00064
00065 else if (p.nMatch("MONI") && p.nMatch("SCIE"))
00066 {
00067
00068
00069 lgPrtSciNot = true;
00070 }
00071
00072 else if (p.nMatch("ATOM") && p.nMatch("DATA"))
00073 {
00074
00075 long int ion;
00076 bool UNUSED lgAs = false;
00077 bool lgCS = false;
00078 bool UNUSED lgDR = false;
00079 const char *chMole;
00080
00081
00082 if (p.nMatch(" AS "))
00083 lgAs = true;
00084
00085 else if (p.nMatch(" DR "))
00086 lgDR = true;
00087
00088 else if (p.nMatch(" CS "))
00089 lgCS = true;
00090 else
00091 {
00092
00093 fprintf(ioQQQ,
00094 " One of the keywords AS DR or CS must appear to change"
00095 " the transition probabilities, dielectronic recombination rate,"
00096 " or collision strength.\n Sorry.\n");
00097 cdEXIT(EXIT_FAILURE);
00098 }
00099
00100
00101 if (p.nMatch(" H2 "))
00102 {
00103
00104 chMole = "H2";
00105 }
00106 else
00107 {
00108
00109 fprintf(ioQQQ, " No molecule was on this SET ATOMIC DATA command.\n Sorry.\n");
00110 fprintf(ioQQQ, " Use SET ATOMIC DATA ION to change an ion.\n Sorry.\n");
00111 cdEXIT(EXIT_FAILURE);
00112 }
00113
00114 if (strcmp(chMole, "H2") == 0)
00115 {
00116 if (p.nMatch(" H ") && lgCS)
00117 {
00118
00119 ion = (long int) p.FFmtRead();
00120 if (ion != 2)
00121 TotalInsanity();
00122 long int nYear = (long) p.FFmtRead();
00123 if (p.lgEOL())
00124 p.NoNumb("collison data set (year)");
00125 if (nYear == 1999)
00126 {
00127
00128
00129
00130 h2.lgH2_H_coll_07 = false;
00131 h2.coll_source[0].filename = "coll_rates_H_99.dat";
00132 }
00133 else if (nYear == 2007)
00134 {
00135
00136
00137
00138 h2.lgH2_H_coll_07 = true;
00139 h2.coll_source[0].filename = "coll_rates_H_07.dat";
00140 }
00141 else
00142 {
00143
00144 fprintf(ioQQQ, " the SET ATOMIC DATA MOLECULE H2"
00145 " H CS command must have year 1999 or 2007.\n");
00146 cdEXIT(EXIT_FAILURE);
00147 }
00148 }
00149 else if (p.nMatch(" HE ") && lgCS)
00150 {
00151
00152 if (p.nMatch("ORNL"))
00153 {
00154
00155
00156 h2.lgH2_He_ORNL = true;
00157 h2.coll_source[1].filename = "coll_rates_He_ORNL.dat";
00158 }
00159 else if( p.nMatch( "BOUR") )
00160 {
00161
00162
00163
00164 h2.lgH2_He_ORNL = false;
00165 h2.coll_source[1].filename = "coll_rates_He_LeBourlot.dat";
00166 }
00167 else
00168 {
00169
00170 fprintf(ioQQQ, " the SET ATOMIC DATA MOLECULE H2"
00171 "He CS command must have ORNL or Le BOURlot.\n");
00172 cdEXIT(EXIT_FAILURE);
00173 }
00174 }
00175 else
00176 {
00177 fprintf( ioQQQ, " the SET ATOMIC DATA H2 CS command requires a collider keyword.\n" );
00178 cdEXIT(EXIT_FAILURE);
00179 }
00180 }
00181 }
00182
00183
00184 else if (p.nMatch("CHECK") && p.nMatch("ENERGY") && p.nMatch("EVERY")
00185 && p.nMatch("ZONE"))
00186 {
00187 continuum.lgCheckEnergyEveryZone = true;
00188 }
00189
00190 else if (p.nMatch("UPDA") && p.nMatch("COUP") &&
00191 p.nMatch("EVER") && p.nMatch("ION") )
00192 {
00193 conv.lgUpdateCouplings = true;
00194 }
00195
00196
00197
00198 else if (p.nMatch(" COLL") && p.nMatch(" IONIZ"))
00199 {
00200
00201 if (p.nMatch(" HYBRID"))
00202 {
00203 atmdat.CIRCData = t_atmdat::HYBRID;
00204 }
00205 else if (p.nMatch(" DIMA"))
00206 {
00207 atmdat.CIRCData = t_atmdat::DIMA;
00208 }
00209 else
00210 {
00211 fprintf(ioQQQ,
00212 " \nPROBLEM Unrecognized set collisional ionization data option.\n");
00213 fprintf(ioQQQ, " Valid options are Dima or Hybrid.\n");
00214 fprintf(ioQQQ, " See Hazy 1 for details.\n");
00215 cdEXIT( EXIT_FAILURE );
00216 }
00217 }
00218
00219 else if( p.nMatch(" H2 " ) && p.nMatch( "CONT" ) && p.nMatch( "DISS" ))
00220 {
00221 if( p.nMatch( "STAN" ) )
00222 {
00223
00224
00225 mole_global.lgStancil = true;
00226 }
00227
00228 else if( p.nMatch( "AD69" ) )
00229 {
00230
00231
00232
00233
00234 mole_global.lgStancil = false;
00235 }
00236 else
00237 {
00238
00239 fprintf( ioQQQ, " There should have been an option on this SET H2 CONTinuum DISSociation command.\n" );
00240 fprintf( ioQQQ, " consult Hazy to find valid options.\n Sorry.\n" );
00241 cdEXIT(EXIT_FAILURE);
00242 }
00243
00244 }
00245
00246 else if( p.nMatch(" CHA") && !p.nMatch( "HO ") && !p.nMatch( "HHE") )
00247 {
00248
00249
00250 atmdat.HCTAlex = p.FFmtRead();
00251 if (p.lgEOL())
00252 {
00253 p.NoNumb("minimum charge transfer rate");
00254 }
00255 if (atmdat.HCTAlex < 0.)
00256 {
00257 atmdat.HCTAlex = pow(10., atmdat.HCTAlex);
00258 }
00259 }
00260 else if( p.nMatch("CHEM") && !p.nMatch( "HO ") && !p.nMatch(" HHE") )
00261 {
00262
00263 if (p.nMatch("FEDE"))
00264 {
00265 if (p.nMatch(" ON "))
00266 {
00267
00268
00269 mole_global.lgFederman = true;
00270 }
00271 else if( p.nMatch( " OFF" ) )
00272 {
00273 mole_global.lgFederman = false;
00274 }
00275 else
00276 {
00277
00278 mole_global.lgFederman = true;
00279 }
00280 }
00281
00282 else if (p.nMatch(" NON") && p.nMatch("EQUI"))
00283 {
00284
00285
00286
00287
00288
00289 mole_global.lgNonEquilChem = true;
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299 if (p.nMatch("NEUT"))
00300 {
00301 if (p.nMatch(" ON "))
00302 {
00303
00304
00305 mole_global.lgNeutrals = true;
00306 }
00307 else if( p.nMatch( " OFF" ) )
00308 {
00309 mole_global.lgNeutrals = false;
00310 }
00311 else
00312 {
00313
00314 mole_global.lgNeutrals = true;
00315 }
00316 }
00317 }
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333 else if (p.nMatch("PROT") && p.nMatch("ELIM"))
00334 {
00335 if (p.nMatch(" ON "))
00336 {
00337
00338
00339 mole_global.lgProtElim = true;
00340 }
00341 else if( p.nMatch( " OFF" ) )
00342 {
00343 mole_global.lgProtElim = false;
00344 }
00345 else
00346 {
00347
00348 mole_global.lgProtElim = true;
00349 }
00350 }
00351
00352 else
00353 {
00354
00355 fprintf(ioQQQ,
00356 " There should have been an option on this SET CHEMISTRY command.\n");
00357 fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
00358 cdEXIT(EXIT_FAILURE);
00359 }
00360 }
00361
00362
00363 else if( p.nMatch("COLL") && p.nMatch("STRE") && p.nMatch("AVER") )
00364 {
00365 if( p.nMatch(" OFF") )
00366 {
00367 iso_ctrl.lgCollStrenThermAver = false;
00368 }
00369 else
00370 {
00371
00372 iso_ctrl.lgCollStrenThermAver = true;
00373 }
00374 }
00375
00376 else if (p.nMatch("COVE"))
00377 {
00378 iterations.lgConverge_set = true;
00379
00380 if (p.nMatch("FAST"))
00381 {
00382 iterations.lim_zone = 1;
00383 iterations.lim_iter = 0;
00384 }
00385 else
00386 {
00387 iterations.lim_zone = 10;
00388 iterations.lim_iter = 1;
00389 }
00390 }
00391
00392 else if (p.nMatch("CSUP"))
00393 {
00394
00395 secondaries.SetCsupra = (realnum) p.FFmtRead();
00396 secondaries.lgCSetOn = true;
00397 if (p.lgEOL())
00398 {
00399 p.NoNumb("secondary ionization rate");
00400 }
00401 secondaries.SetCsupra = (realnum) pow((realnum) 10.f,
00402 secondaries.SetCsupra);
00403 }
00404
00405 else if (p.nMatch(" D/H"))
00406 {
00407
00408 hydro.D2H_ratio = p.FFmtRead();
00409 if (hydro.D2H_ratio <= 0. || p.nMatch(" LOG"))
00410 {
00411 hydro.D2H_ratio = pow(10., hydro.D2H_ratio);
00412 }
00413 if (p.lgEOL())
00414 {
00415 p.NoNumb("D to H ratio");
00416 }
00417 }
00418
00419 else if( p.nMatch("DENS") && p.nMatch("TOLE") )
00420 {
00421
00422 conv.GasPhaseAbundErrorAllowed = (realnum)p.FFmtRead();
00423 if( p.lgEOL() )
00424 {
00425 p.NoNumb("density tolerance");
00426 }
00427 if( conv.GasPhaseAbundErrorAllowed <= 0. )
00428 {
00429 conv.GasPhaseAbundErrorAllowed = (realnum)pow((realnum)10.f,conv.GasPhaseAbundErrorAllowed);
00430 }
00431 }
00432
00433 else if( p.nMatch( " HO " ) && p.nMatch( "CHAR" ) )
00434 {
00435 fprintf(ioQQQ, " The SET HO CHAR command is no longer supported.\n" );
00436 cdEXIT(EXIT_FAILURE);
00437 }
00438
00439 else if( p.nMatch( " HHE" ) && p.nMatch( "CHAR" ) )
00440 {
00441 fprintf(ioQQQ, " The SET HHE CHAR command is no longer supported.\n" );
00442 cdEXIT(EXIT_FAILURE);
00443 }
00444
00445 else if( p.nMatch("12C1") )
00446 {
00447
00448
00449 co.C12_C13_isotope_ratio = (realnum) p.FFmtRead();
00450 co.C12_C13_isotope_ratio = (realnum) p.FFmtRead();
00451
00452
00453 co.C12_C13_isotope_ratio = (realnum) p.FFmtRead();
00454 if (p.lgEOL())
00455 p.NoNumb("12C to 13C abundance ratio");
00456
00457 if (co.C12_C13_isotope_ratio <= 0. || p.nMatch(" LOG"))
00458 co.C12_C13_isotope_ratio = (realnum) pow((realnum) 10.f,
00459 co.C12_C13_isotope_ratio);
00460 }
00461
00462
00463 else if (p.nMatch("DYNA"))
00464 {
00465
00466 if (p.nMatch("ADVE") && p.nMatch("LENG"))
00467 {
00468
00469 dynamics.AdvecLengthInit = p.FFmtRead();
00470 if (p.lgEOL())
00471 p.NoNumb("advection length");
00472
00473
00474 if (p.nMatch("FRAC"))
00475 {
00476
00477 if (dynamics.AdvecLengthInit <= 0.)
00478 dynamics.AdvecLengthInit = pow(10.,
00479 dynamics.AdvecLengthInit);
00480
00481
00482 dynamics.AdvecLengthInit *= -1.;
00483 }
00484 else
00485 {
00486
00487
00488 dynamics.AdvecLengthInit = pow(10., dynamics.AdvecLengthInit);
00489 }
00490 }
00491 else if (p.nMatch("PRES") && p.nMatch("MODE"))
00492 {
00493 dynamics.lgSetPresMode = true;
00494 if (p.nMatch("SUBS"))
00495 {
00496
00497 strcpy(dynamics.chPresMode, "subsonic");
00498 }
00499 else if (p.nMatch("SUPE"))
00500 {
00501
00502 strcpy(dynamics.chPresMode, "supersonic");
00503 }
00504 else if (p.nMatch("STRO"))
00505 {
00506
00507 strcpy(dynamics.chPresMode, "strongd");
00508 }
00509 else if (p.nMatch("ORIG"))
00510 {
00511
00512 strcpy(dynamics.chPresMode, "original");
00513 }
00514 }
00515 else if (p.nMatch("ANTI") && p.nMatch("DEPT"))
00516 {
00517 dynamics.lgSetPresMode = true;
00518 strcpy(dynamics.chPresMode, "antishock");
00519
00520
00521 dynamics.ShockDepth = p.FFmtRead();
00522 if (p.lgEOL())
00523 p.NoNumb("antishock depth");
00524 dynamics.ShockDepth = pow(10., dynamics.ShockDepth);
00525 }
00526 else if (p.nMatch("ANTI") && p.nMatch("MACH"))
00527 {
00528 dynamics.lgSetPresMode = true;
00529 strcpy(dynamics.chPresMode, "antishock-by-mach");
00530
00531
00532 dynamics.ShockMach = p.FFmtRead();
00533 if (p.lgEOL())
00534 p.NoNumb("antishock-by-mach");
00535 }
00536 else if (p.nMatch("RELA"))
00537 {
00538
00539
00540 dynamics.n_initial_relax = (long int) p.FFmtRead();
00541 if (p.lgEOL())
00542 p.NoNumb("relaxation cycles before start of dynamics");
00543 else if (dynamics.n_initial_relax < 2)
00544 {
00545 fprintf(ioQQQ,
00546 " First iteration to relax dynamics must be > 1."
00547 "It was %li. Sorry.\n", dynamics.n_initial_relax);
00548 cdEXIT(EXIT_FAILURE);
00549 }
00550 }
00551 else if (p.nMatch("SHOC") && p.nMatch("DEPT"))
00552 {
00553 dynamics.lgSetPresMode = true;
00554 strcpy(dynamics.chPresMode, "shock");
00555
00556
00557 dynamics.ShockDepth = p.FFmtRead();
00558 if (p.lgEOL())
00559 p.NoNumb("shock depth");
00560 dynamics.ShockDepth = pow(10., dynamics.ShockDepth);
00561 }
00562 else if (p.nMatch("POPU") && p.nMatch("EQUI"))
00563 {
00564
00565
00566 dynamics.lgEquilibrium = true;
00567 }
00568 else
00569 {
00570
00571 fprintf(ioQQQ,
00572 " There should have been an option on this SET DYNAMICS command.\n");
00573 fprintf(ioQQQ, " consult Hazy to find valid options.\n Sorry.\n");
00574 cdEXIT(EXIT_FAILURE);
00575 }
00576 }
00577
00578 else if (p.nMatch("DIDZ"))
00579 {
00580
00581
00582
00583 radius.drChange = (realnum) p.FFmtRead();
00584 if (radius.drChange <= 0.)
00585 {
00586 radius.drChange = (realnum) pow((realnum) 10.f, radius.drChange);
00587 }
00588 if (p.lgEOL())
00589 {
00590 p.NoNumb("largest optical depth allowed in zone");
00591 }
00592 }
00593
00594
00595 else if (p.nMatch("EDEN"))
00596 {
00597
00598
00599 if (p.nMatch("CONV") || p.nMatch("ERRO"))
00600 {
00601
00602
00603 conv.EdenErrorAllowed = p.FFmtRead();
00604 if (p.lgEOL())
00605 {
00606 p.NoNumb("electron density error allowed");
00607 }
00608
00609 if (conv.EdenErrorAllowed < 0.)
00610 {
00611 conv.EdenErrorAllowed = pow(10., conv.EdenErrorAllowed);
00612 }
00613 }
00614 else if(p.nMatch("FRACTION"))
00615 {
00616
00617 dense.EdenFraction = p.FFmtRead();
00618 if (p.lgEOL())
00619 {
00620 p.NoNumb("electron density fraction");
00621 }
00622
00623 if (dense.EdenFraction <= 0. )
00624 {
00625 dense.EdenFraction = pow(10., dense.EdenFraction);
00626 }
00627
00628 phycon.lgPhysOK = false;
00629 }
00630 else
00631 {
00632
00633
00634 dense.EdenSet = (realnum) pow(10., p.FFmtRead());
00635 if (p.lgEOL())
00636 {
00637 p.NoNumb("electron density");
00638 }
00639
00640
00641 phycon.lgPhysOK = false;
00642 }
00643 }
00644
00645
00646 else if( p.nMatch("GBAR"))
00647 {
00648 atmdat.lgGbarOn = false;
00649 }
00650
00651
00652 else if (p.nMatch("IONI") && p.nMatch("TOLE"))
00653 {
00654
00655 conv.IonizErrorAllowed = p.FFmtRead();
00656 if (p.lgEOL())
00657 {
00658 p.NoNumb("ionization error allowed");
00659 }
00660
00661 if (conv.IonizErrorAllowed < 0.)
00662 {
00663 conv.IonizErrorAllowed = pow(10., (double)conv.IonizErrorAllowed);
00664 }
00665 }
00666
00667
00668 else if (p.nMatch("ISOTOPE") || p.nMatch("ISOTOPO") )
00669 {
00670 if( p.nMatch(" ALL ") )
00671 {
00672
00673 for( long nelem = ipHELIUM ; nelem < LIMELM; ++nelem )
00674 {
00675 mole_global.lgTreatIsotopes[nelem] = true;
00676 }
00677 }
00678 }
00679
00680 else if (p.nMatch("FINE") && p.nMatch("CONT"))
00681 {
00682
00683
00684
00685
00686
00687 if ((rfield.fine_opac_nelem = p.GetElem()) < 0)
00688 {
00689 fprintf(ioQQQ,
00690 " An element name must appear on this line\n Sorry.\n");
00691 cdEXIT(EXIT_FAILURE);
00692 }
00693
00694
00695 rfield.fine_opac_nresolv = (long int) p.FFmtRead();
00696 if (rfield.fine_opac_nresolv < 1)
00697 {
00698 fprintf(ioQQQ,
00699 " The number of resolution elements within FWHM of line must appear\n Sorry.\n");
00700 cdEXIT(EXIT_FAILURE);
00701 }
00702
00703
00704
00705 if (!p.lgEOL())
00706 {
00707 realnum lower_limit = (realnum) p.FFmtRead();
00708 if (lower_limit < 0)
00709 lower_limit = pow((realnum) 10.f, lower_limit);
00710
00711 if (lower_limit > 0.2f)
00712 fprintf(
00713 ioQQQ,
00714 " The fine continuum lower limit is quite high (%f Ryd). Please check.\n",
00715 lower_limit);
00716
00717
00718 rfield.fine_ener_lo = MAX2(rfield.emm, lower_limit);
00719
00720 if (!p.lgEOL())
00721 {
00722 realnum upper_limit = (realnum) p.FFmtRead();
00723 if (upper_limit < 0)
00724 upper_limit = pow((realnum) 10.f, upper_limit);
00725
00726 if (upper_limit < 10.f)
00727 fprintf(
00728 ioQQQ,
00729 " The fine continuum upper limit is quite low (%f Ryd). Please check.\n",
00730 upper_limit);
00731
00732
00733 rfield.fine_ener_hi = MIN2(rfield.egamry, upper_limit);
00734 }
00735 }
00736 }
00737
00738
00739 else if (p.nMatch("GRAI") && !p.nMatch(" H2 "))
00740 {
00741 if (p.nMatch("HEAT"))
00742 {
00743
00744 gv.GrainHeatScaleFactor = (realnum) p.FFmtRead();
00745
00746 phycon.lgPhysOK = false;
00747 if (p.lgEOL())
00748 {
00749 p.NoNumb("grain heating");
00750 }
00751 }
00752 else
00753 {
00754 fprintf(ioQQQ,
00755 " A keyword must appear on the SET GRAIN line - options are HEAT \n Sorry.\n");
00756 cdEXIT(EXIT_FAILURE);
00757 }
00758 }
00759
00760
00761
00762 else if (p.nMatch("LEID") && p.nMatch("HACK"))
00763 {
00764 if (p.nMatch("H2* ") && p.nMatch(" OFF"))
00765 {
00766
00767 hmi.lgLeiden_Keep_ipMH2s = false;
00768
00769 phycon.lgPhysOK = false;
00770 }
00771 else if (p.nMatch("CR ") && p.nMatch(" OFF"))
00772 {
00773
00774 hmi.lgLeidenCRHack = false;
00775 }
00776 else if( p.nMatch("RATE") && p.nMatch("UMIS"))
00777 {
00778
00779
00780 mole_global.lgLeidenHack = true;
00781 }
00782 }
00783
00784
00785 else if (p.nMatch(" H2 "))
00786 {
00787 ip = (long int) p.FFmtRead();
00788 if (ip != 2)
00789 {
00790 fprintf(ioQQQ,
00791 " The first number on this line must be the 2 in H2\n Sorry.\n");
00792 cdEXIT(EXIT_FAILURE);
00793 }
00794
00795
00796 if (p.nMatch("SOLO") || (p.nMatch("SMAL") && p.nMatch("MODE")))
00797 {
00798 if (p.nMatch("SOLO"))
00799 {
00800
00801 fprintf(ioQQQ,
00802 "PROBLEM - *set H2 Solomon* has been changed to *set H2 small model*."
00803 " This is OK for now but it may not work in a future version.\n");
00804 }
00805 if (p.nMatch("TH85"))
00806 {
00807
00808 hmi.chH2_small_model_type = 'T';
00809 }
00810 else if (p.nMatch(" BHT"))
00811 {
00812
00813
00814
00815 hmi.chH2_small_model_type = 'H';
00816 }
00817 else if (p.nMatch("BD96"))
00818 {
00819
00820
00821
00822 hmi.chH2_small_model_type = 'B';
00823 }
00824 else if (p.nMatch("ELWE"))
00825 {
00826
00827
00828
00829 hmi.chH2_small_model_type = 'E';
00830 }
00831 else
00832 {
00833 fprintf(ioQQQ,
00834 " One of the keywords TH85, _BHT, BD96 or ELWErt must appear.\n Sorry.\n");
00835 cdEXIT(EXIT_FAILURE);
00836 }
00837 }
00838
00839
00840
00841 if (p.nMatch("GRAI") && p.nMatch("FORM") && p.nMatch("PUMP"))
00842 {
00843 if (p.nMatch("DB96"))
00844 {
00845
00846 hmi.chGrainFormPump = 'D';
00847 }
00848 else if (p.nMatch("TAKA"))
00849 {
00850
00851
00852 hmi.chGrainFormPump = 'T';
00853 }
00854 else if (p.nMatch("THER"))
00855 {
00856
00857
00858 hmi.chGrainFormPump = 't';
00859 }
00860 else
00861 {
00862 fprintf(ioQQQ,
00863 " The grain form pump option is wrong.\n Sorry.\n");
00864 cdEXIT(EXIT_FAILURE);
00865 }
00866 }
00867
00868
00869 else if (p.nMatch("JURA"))
00870 {
00871 if (p.nMatch("TH85"))
00872 {
00873
00874 hmi.chJura = 'T';
00875 }
00876 else if (p.nMatch("CT02"))
00877 {
00878
00879 hmi.chJura = 'C';
00880 }
00881 else if (p.nMatch("SN99"))
00882 {
00883
00884 hmi.chJura = 'S';
00885 }
00886 else if (p.nMatch("RATE"))
00887 {
00888
00889
00890
00891
00892 hmi.chJura = 'F';
00893 hmi.rate_h2_form_grains_set = pow(10., p.FFmtRead());
00894 if (p.lgEOL())
00895 {
00896
00897
00898 hmi.rate_h2_form_grains_set = 3e-17;
00899 }
00900 }
00901 else if (p.nMatch("SCAL"))
00902 {
00903
00904 hmi.ScaleJura = (realnum) p.FFmtRead();
00905
00906 if (p.nMatch(" LOG") || hmi.ScaleJura < 0.)
00907 {
00908 hmi.ScaleJura = (realnum) pow(10., (double) hmi.ScaleJura);
00909 }
00910 if (p.lgEOL())
00911 p.NoNumb("scale for Jura rate");
00912
00913
00914 if (optimize.lgVarOn)
00915 {
00916 optimize.nvarxt[optimize.nparm] = 1;
00917 strcpy(optimize.chVarFmt[optimize.nparm],
00918 "SET H2 JURA SCALE %f LOG");
00919
00920
00921 optimize.nvfpnt[optimize.nparm] = input.nRead;
00922
00923
00924 optimize.vparm[0][optimize.nparm] = (realnum) log10(
00925 hmi.ScaleJura);
00926 optimize.vincr[optimize.nparm] = 0.3f;
00927
00928 ++optimize.nparm;
00929 }
00930 }
00931 else
00932 {
00933 fprintf(ioQQQ, " The Jura rate option is wrong.\n Sorry.\n");
00934 cdEXIT(EXIT_FAILURE);
00935 }
00936 }
00937
00938
00939 else if (p.nMatch(" TAD"))
00940 {
00941 hmi.Tad = (realnum) p.FFmtRead();
00942 if (p.lgEOL())
00943 p.NoNumb("temperature for binding energy");
00944
00945 if (hmi.Tad <= 10. && !p.nMatch("LINE"))
00946 hmi.Tad = (realnum) pow((realnum) 10.f, hmi.Tad);
00947 }
00948
00949 else if (p.nMatch("FRAC"))
00950 {
00951
00952
00953
00954
00955 hmi.H2_frac_abund_set = p.FFmtRead();
00956 if (p.lgEOL())
00957 p.NoNumb("H2 fractional abundance");
00958
00959
00960 if (hmi.H2_frac_abund_set <= 0.)
00961 hmi.H2_frac_abund_set = pow(10., hmi.H2_frac_abund_set);
00962
00963
00964 hmi.H2_frac_abund_set = MIN2(0.49999, hmi.H2_frac_abund_set);
00965 }
00966 #if 0
00967 else if( p.nMatch("FORM") && p.nMatch("SCAL") )
00968 {
00969
00970
00971
00972
00973
00974 hmi.H2_formation_scale = p.FFmtRead();
00975 if (p.lgEOL())
00976 p.NoNumb("H2 formation scale");
00977
00978
00979 if (hmi.H2_formation_scale <= 0.)
00980 hmi.H2_formation_scale = pow(10., hmi.H2_frac_abund_set);
00981 }
00982 #endif
00983 }
00984
00985
00986
00987
00988 else if (p.nMatch("HCOR"))
00989 {
00990 dense.HCorrFac = (realnum) p.FFmtRead();
00991 if (p.lgEOL())
00992 p.NoNumb("scale for H0 correction to e- collision rate");
00993 }
00994
00995 else if (p.nMatch(" PAH"))
00996 {
00997
00998 if (lgQuotesFound)
00999 {
01000
01001 gv.chPAH_abundance = chString_quotes_lowercase;
01002 }
01003 else if (p.nMatch("CONS"))
01004 {
01005
01006 gv.chPAH_abundance = "CON";
01007 }
01008 else if (p.nMatch("BAKE"))
01009 {
01010
01011
01012 gv.lgBakesPAH_heat = true;
01013
01014 phycon.lgPhysOK = false;
01015 }
01016 else
01017 {
01018 fprintf(ioQQQ,
01019 " a string, or one of the keywords BAKES, or CONStant must appear.\n Sorry.");
01020 cdEXIT(EXIT_FAILURE);
01021 }
01022 }
01023
01024 else if (p.nMatch("PRES") && p.nMatch("IONI"))
01025 {
01026
01027
01028 conv.limPres2Ioniz = (long) p.FFmtRead();
01029 if (p.lgEOL())
01030 {
01031 p.NoNumb("number of calls from pressure to ion solver");
01032 }
01033 else if (conv.limPres2Ioniz <= 0)
01034 {
01035 fprintf(ioQQQ, " The limit must be greater than zero.\n Sorry.");
01036 cdEXIT(EXIT_FAILURE);
01037 }
01038 }
01039
01040 else if (p.nMatch("PRES"))
01041 {
01042
01043 if (p.nMatch("CONV"))
01044 {
01045
01046
01047 conv.PressureErrorAllowed = (realnum) p.FFmtRead();
01048 if (p.lgEOL())
01049 p.NoNumb("pressure convergence tolerance");
01050
01051 if (conv.PressureErrorAllowed < 0.)
01052 conv.PressureErrorAllowed = (realnum) pow((realnum) 10.f,
01053 conv.PressureErrorAllowed);
01054 }
01055
01056 else
01057 {
01058
01059
01060 fprintf(ioQQQ,
01061 " I didn\'t recognize a key on this SET PRESSURE line.\n");
01062 fprintf(ioQQQ, " The ones I know about are: CONVergence.\n");
01063 cdEXIT(EXIT_FAILURE);
01064 }
01065 }
01066 else if (p.nMatch("RECOMBIN"))
01067 {
01068
01069 if (p.nMatch("DIELECTR"))
01070 {
01071 if (p.nMatch("MEAN"))
01072 {
01073
01074 if (p.nMatch(" OFF"))
01075 {
01076 for (int ion = 0; ion < LIMELM; ++ion)
01077 ionbal.DR_mean_scale[ion] = 0.;
01078 }
01079
01080 else if (p.nMatch("SCALE"))
01081 {
01082
01083 ionbal.DR_mean_scale[0] = (realnum) p.FFmtRead();
01084 if (p.lgEOL())
01085 {
01086 fprintf(
01087 ioQQQ,
01088 " There must be at least one scale factor on the SET RECOMBIANTION MEAN command.\n");
01089 cdEXIT(EXIT_FAILURE);
01090 }
01091
01092 for (int ion = 1; ion < LIMELM; ++ion)
01093 {
01094 ionbal.DR_mean_scale[ion] = (realnum) p.FFmtRead();
01095 if (p.lgEOL())
01096 ionbal.DR_mean_scale[ion]
01097 = ionbal.DR_mean_scale[ion - 1];
01098 }
01099 for (int ion = 0; ion < LIMELM; ++ion)
01100 {
01101 if (ionbal.DR_mean_scale[ion] < 0)
01102 {
01103 fprintf(ioQQQ,
01104 " All scale factors on the SET RECOMBIANTION MEAN command must be >=0.\n");
01105 cdEXIT(EXIT_FAILURE);
01106 }
01107 }
01108 }
01109
01110 else if (p.nMatch("NOISE"))
01111 {
01112 ionbal.guess_noise = (realnum) p.FFmtRead();
01113 if (p.lgEOL())
01114 ionbal.guess_noise = 2.;
01115
01116
01117 init_genrand((unsigned) time(NULL));
01118 }
01119 else
01120 {
01121 fprintf(ioQQQ, " key OFF or NOISE must appear.\n");
01122 cdEXIT(EXIT_FAILURE);
01123 }
01124 }
01125 else
01126 {
01127 fprintf(ioQQQ, " key MEAN must appear.\n");
01128 cdEXIT(EXIT_FAILURE);
01129 }
01130 }
01131 else
01132 {
01133 fprintf(ioQQQ,
01134 " key DIELECTRonic must appear on set recombination command.\n");
01135 cdEXIT(EXIT_FAILURE);
01136 }
01137 }
01138
01139 else if (p.nMatch(" DR "))
01140 {
01141
01142
01143 radius.sdrmax = p.FFmtRead();
01144 if (!p.nMatch("LINE"))
01145 {
01146
01147 radius.sdrmax = pow(10., radius.sdrmax);
01148 }
01149 if (p.lgEOL())
01150 {
01151 p.NoNumb("zone thickness");
01152 }
01153
01154 radius.lgSdrmaxRel = p.nMatch("RELA");
01155
01156
01157 radius.lgFixed = true;
01158 radius.sdrmin = radius.sdrmax;
01159 radius.lgSdrminRel = radius.lgSdrmaxRel;
01160 if (!radius.lgSdrmaxRel && radius.sdrmax < DEPTH_OFFSET * 1e4)
01161 {
01162 fprintf(
01163 ioQQQ,
01164 "\n Thicknesses less than about %.0e will NOT give accurate results. If tricking the code\n",
01165 DEPTH_OFFSET * 1e4);
01166 fprintf(
01167 ioQQQ,
01168 " into computing emissivities instead of intensities, try to instead use a thickness of unity,\n");
01169 fprintf(
01170 ioQQQ,
01171 " and then multiply (divide) the results by the necessary thickness (product of densities).\n\n");
01172 cdEXIT(EXIT_FAILURE);
01173 }
01174 if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
01175 {
01176 fprintf(
01177 ioQQQ,
01178 " When using a relative dr, a fraction between 0 and 1 must be entered. Found: %g\n",
01179 radius.sdrmax);
01180 cdEXIT(EXIT_FAILURE);
01181 }
01182 }
01183
01184 else if (p.nMatch("DRMA"))
01185 {
01186
01187 radius.sdrmax = p.FFmtRead();
01188 if (p.lgEOL())
01189 p.NoNumb("maximum zone thickness");
01190
01191 if ((radius.sdrmax < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
01192 radius.sdrmax = pow(10., radius.sdrmax);
01193
01194
01195 radius.lgSdrmaxRel = p.nMatch("RELA");
01196 if (radius.lgSdrmaxRel && (radius.sdrmax <= 0. || radius.sdrmax >= 1.))
01197 {
01198 fprintf(
01199 ioQQQ,
01200 " When using a relative drmax, a fraction between 0 and 1 must be entered. Found: %g\n",
01201 radius.sdrmax);
01202 cdEXIT(EXIT_FAILURE);
01203 }
01204 }
01205
01206 else if (p.nMatch("DRMI") && p.nMatch("DEPT"))
01207 {
01208
01209 radius.sdrmin_rel_depth = p.FFmtRead();
01210 if (p.lgEOL())
01211 p.NoNumb("minimum zone thickness rel to depth");
01212
01213 if ((radius.sdrmin_rel_depth < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
01214 radius.sdrmin_rel_depth = pow(10., radius.sdrmin_rel_depth );
01215
01216 if( radius.sdrmin_rel_depth <= 0. || radius.sdrmin_rel_depth >= 1.)
01217 {
01218 fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
01219 radius.sdrmin_rel_depth );
01220 cdEXIT(EXIT_FAILURE);
01221 }
01222 }
01223
01224 else if (p.nMatch("DRMI"))
01225 {
01226
01227 radius.sdrmin = p.FFmtRead();
01228 if (p.lgEOL())
01229 p.NoNumb("minimum zone thickness");
01230
01231 if ((radius.sdrmin < 38. || p.nMatch(" LOG")) && !p.nMatch("LINE"))
01232 radius.sdrmin = pow(10., radius.sdrmin);
01233
01234
01235 radius.lgSdrminRel = p.nMatch("RELA");
01236 radius.lgSMinON = true;
01237 if (radius.lgSdrminRel && (radius.sdrmin <= 0. || radius.sdrmin >= 1.))
01238 {
01239 fprintf( ioQQQ, " When using a relative drmin, a fraction between 0 and 1 must be entered. Found: %g\n",
01240 radius.sdrmin);
01241 cdEXIT(EXIT_FAILURE);
01242 }
01243 }
01244
01245 else if (p.nMatch("FLXF"))
01246 {
01247
01248 rfield.FluxFaint = (realnum) p.FFmtRead();
01249 if (p.lgEOL())
01250 {
01251 p.NoNumb("faintest continuum flux to consider");
01252 }
01253 if (rfield.FluxFaint < 0.)
01254 {
01255 rfield.FluxFaint = (realnum) pow((realnum) 10.f, rfield.FluxFaint);
01256 }
01257 }
01258
01259 else if (p.nMatch("LINE") && p.nMatch("PREC"))
01260 {
01261
01262
01263 LineSave.sig_figs = (int) p.FFmtRead();
01264 if (p.lgEOL())
01265 {
01266 p.NoNumb("line precision");
01267 }
01268 else if (LineSave.sig_figs > 6)
01269 {
01270 fprintf(ioQQQ,
01271 " set line precision currently only works for up to 6 significant figures.\n");
01272 cdEXIT(EXIT_FAILURE);
01273 }
01274
01275
01276 }
01277
01278
01279 else if (p.nMatch("NFNU"))
01280 {
01281 if (p.nMatch(" ADD"))
01282 {
01283
01284 double energy = p.FFmtRead();
01285 if (p.lgEOL())
01286 p.NoNumb("energy");
01287 t_PredCont::Inst().add(energy, p.StandardEnergyUnit());
01288 }
01289 else
01290 {
01291
01292
01293
01294
01295 prt.lgSourceReflected = p.nMatch("INCIDENT R") || p.nMatch(
01296 "INCIDENT_R");
01297
01298 prt.lgSourceTransmitted = p.nMatch("INCIDENT_T") || p.nMatch(
01299 "INCIDENT T");
01300
01301 prt.lgDiffuseInward = p.nMatch("DIFFUSE_I")
01302 || p.nMatch("DIFFUSE I");
01303
01304 prt.lgDiffuseOutward = p.nMatch("DIFFUSE_O") || p.nMatch(
01305 "DIFFUSE O");
01306
01307
01308 if (!(prt.lgSourceReflected || prt.lgSourceTransmitted
01309 || prt.lgDiffuseInward || prt.lgDiffuseOutward))
01310 {
01311 fprintf(ioQQQ,
01312 " set nFnu expects one or more of the following keywords:\n");
01313 fprintf(ioQQQ, " INCIDENT_REFLECTED, INCIDENT_TRANSMITTED, "
01314 "DIFFUSE_INWARD, DIFFUSE_OUTWARD\n");
01315 cdEXIT(EXIT_FAILURE);
01316 }
01317 }
01318 }
01319
01320 else if (p.nMatch("IND2"))
01321 {
01322 if (p.nMatch(" ON "))
01323 {
01324
01325 iso_ctrl.lgInd2nu_On = true;
01326 }
01327 else if( p.nMatch(" OFF") )
01328 {
01329 iso_ctrl.lgInd2nu_On = false;
01330 }
01331 else
01332 {
01333 fprintf( ioQQQ, " set ind2 needs either ON or OFF.\n" );
01334 cdEXIT(EXIT_FAILURE);
01335 }
01336 }
01337
01338
01339
01340 else if (p.nMatch("SPECIES"))
01341 {
01342 if (p.nMatch(" GBAR"))
01343 {
01344 atmdat.collstrDefault = (long)p.FFmtRead();
01345 if (p.lgEOL())
01346 {
01347 p.NoNumb("the default collision strengths when no collision or radiative data are available");
01348 }
01349 }
01350 else
01351 {
01352 fprintf(ioQQQ, " SET SPECIES takes option GBAR.\n");
01353 cdEXIT(EXIT_FAILURE);
01354 }
01355 }
01356
01357 else if (p.nMatch("TEMP"))
01358 {
01359
01360 if (p.nMatch("FLOO"))
01361 {
01362 StopCalc.TeFloor = (realnum) p.FFmtRead();
01363 if (p.lgEOL())
01364 p.NoNumb("temperature floor");
01365
01366
01367 if (p.nMatch(" LOG") || (StopCalc.TeFloor <= 10. && !p.nMatch(
01368 "LINE")))
01369 StopCalc.TeFloor = (realnum) pow(10., StopCalc.TeFloor);
01370
01371 if (StopCalc.TeFloor < phycon.TEMP_LIMIT_LOW)
01372 {
01373 fprintf(ioQQQ, " TE < %gK, reset to %gK.\n",
01374 phycon.TEMP_LIMIT_LOW, 1.0001 * phycon.TEMP_LIMIT_LOW);
01375 StopCalc.TeFloor = realnum(1.0001 * phycon.TEMP_LIMIT_LOW);
01376 }
01377
01378 if (StopCalc.TeFloor > phycon.TEMP_LIMIT_HIGH)
01379 {
01380 fprintf(ioQQQ,
01381 " TE > %gK. Cloudy cannot handle this. Bailing out.\n",
01382 phycon.TEMP_LIMIT_HIGH);
01383 cdEXIT(EXIT_FAILURE);
01384 }
01385
01386
01387 if (optimize.lgVarOn)
01388 {
01389 optimize.nvarxt[optimize.nparm] = 1;
01390 strcpy(optimize.chVarFmt[optimize.nparm],
01391 "SET TEMPERATURE FLOOR %f LOG");
01392
01393
01394 optimize.nvfpnt[optimize.nparm] = input.nRead;
01395
01396
01397 optimize.vparm[0][optimize.nparm] = (realnum) log10(
01398 StopCalc.TeFloor);
01399 optimize.varang[optimize.nparm][0] = (realnum) log10(
01400 1.00001 * phycon.TEMP_LIMIT_LOW);
01401 optimize.varang[optimize.nparm][1] = (realnum) log10(
01402 0.99999 * phycon.TEMP_LIMIT_HIGH);
01403 optimize.vincr[optimize.nparm] = 0.3f;
01404
01405 ++optimize.nparm;
01406 }
01407 }
01408
01409 else if (p.nMatch("CONV") || p.nMatch("TOLE"))
01410 {
01411
01412 conv.HeatCoolRelErrorAllowed = (realnum) p.FFmtRead();
01413 if (p.lgEOL())
01414 {
01415 p.NoNumb("heating cooling tolerance");
01416 }
01417 if (conv.HeatCoolRelErrorAllowed <= 0.)
01418 {
01419 conv.HeatCoolRelErrorAllowed = (realnum) pow((realnum) 10.f,
01420 conv.HeatCoolRelErrorAllowed);
01421 }
01422 }
01423
01424 else
01425 {
01426 fprintf(ioQQQ,
01427 "\nI did not recognize a keyword on this SET TEMPERATURE command.\n");
01428 p.PrintLine(ioQQQ);
01429 fprintf(ioQQQ, "The keywords are FLOOr and CONVergence.\n");
01430 cdEXIT(EXIT_FAILURE);
01431 }
01432 }
01433
01434 else if (p.nMatch("TEST"))
01435 {
01436
01437 lgTestCodeEnabled = true;
01438 }
01439
01440 else if (p.nMatch("TRIM"))
01441 {
01442
01443
01444
01445 if (p.nMatch("UPPE"))
01446 {
01447
01448 double save = ionbal.trimhi;
01449 ionbal.trimhi = pow(10., p.FFmtRead());
01450 if (p.lgEOL() && p.nMatch(" OFF"))
01451 {
01452
01453 p.setEOL(false);
01454 ionbal.lgTrimhiOn = false;
01455
01456 ionbal.trimhi = save;
01457 }
01458 }
01459
01460 else if (p.nMatch("LOWE"))
01461 {
01462
01463 ionbal.trimlo = pow(10., p.FFmtRead());
01464 }
01465
01466
01467 else if (p.nMatch("SMAL") || p.nMatch(" OFF"))
01468 {
01469
01470 ionbal.trimlo = SMALLFLOAT;
01471 ionbal.trimhi = SMALLFLOAT;
01472 p.setEOL(false);
01473 }
01474
01475 else
01476 {
01477
01478 ionbal.trimhi = pow(10., p.FFmtRead());
01479
01480
01481 ionbal.trimlo = ionbal.trimhi;
01482 }
01483
01484 if (p.lgEOL())
01485 {
01486 p.NoNumb("trimming parameter");
01487 }
01488
01489 if (ionbal.trimlo >= 1. || ionbal.trimhi >= 1.)
01490 {
01491 fprintf(ioQQQ, " number must be negative since log\n");
01492 cdEXIT(EXIT_FAILURE);
01493 }
01494 }
01495
01496 else if (p.nMatch("SKIP"))
01497 {
01498
01499 save.ncSaveSkip = (long) p.FFmtRead();
01500 if (p.lgEOL())
01501 {
01502 p.NoNumb("number of points to skip in save");
01503 }
01504 }
01505
01506 else if (p.nMatch(" UTA"))
01507 {
01508
01509
01510 if (p.nMatch("GU06"))
01511 {
01512 if (p.nMatch(" OFF"))
01513 {
01514
01515 ionbal.lgInnerShell_Gu06 = false;
01516 }
01517 else
01518 {
01519
01520
01521
01522 ionbal.lgInnerShell_Gu06 = true;
01523 }
01524 }
01525 else if (p.nMatch("KISI"))
01526 {
01527 if (p.nMatch(" OFF"))
01528 {
01529
01530 ionbal.lgInnerShell_Kisielius = false;
01531 }
01532 else
01533 {
01534
01535 ionbal.lgInnerShell_Kisielius = true;
01536 }
01537 }
01538 }
01539
01540 else if (p.nMatch("WEAKH"))
01541 {
01542
01543 save.WeakHeatCool = (realnum) p.FFmtRead();
01544
01545 if (p.lgEOL())
01546 {
01547 p.NoNumb("threshold on save heating and cooling");
01548 }
01549
01550 if (save.WeakHeatCool < 0.)
01551 {
01552 save.WeakHeatCool
01553 = (realnum) pow((realnum) 10.f, save.WeakHeatCool);
01554 }
01555 }
01556
01557 else if (p.nMatch("KSHE"))
01558 {
01559
01560 continuum.EnergyKshell = (realnum) p.FFmtRead();
01561 if (p.lgEOL())
01562 {
01563 p.NoNumb("k-shell ionization opacity limit");
01564 }
01565
01566 if (continuum.EnergyKshell == 0.)
01567 {
01568
01569 continuum.EnergyKshell = rfield.egamry;
01570 }
01571
01572 else if (continuum.EnergyKshell < 194.)
01573 {
01574 fprintf(ioQQQ, " k-shell energy must be greater than 194 Ryd\n");
01575 cdEXIT(EXIT_FAILURE);
01576 }
01577 }
01578
01579 else if (p.nMatch("NCHR"))
01580 {
01581
01582 double val = p.FFmtRead();
01583
01584 if (p.lgEOL())
01585 {
01586 p.NoNumb("number of charge states");
01587 }
01588 else
01589 {
01590 long nChrg = nint(val);
01591 if (nChrg < 2 || nChrg > NCHS)
01592 {
01593 fprintf(ioQQQ,
01594 " illegal value for number of charge states: %ld\n",
01595 nChrg);
01596 fprintf(ioQQQ, " choose a value between 2 and %d\n", NCHS);
01597 fprintf(ioQQQ,
01598 " or increase NCHS in grainvar.h and recompile\n");
01599 cdEXIT(EXIT_FAILURE);
01600 }
01601 else
01602 {
01603 SetNChrgStates(nChrg);
01604 }
01605 }
01606 }
01607
01608 else if (p.nMatch("NEGO"))
01609 {
01610
01611 opac.lgNegOpacIO = true;
01612 }
01613
01614 else if (p.nMatch("NEND"))
01615 {
01616
01617
01618
01619
01620 if (geometry.lgEndDflt)
01621 {
01622
01623 geometry.nEndDflt = (long) p.FFmtRead();
01624 geometry.lgEndDflt = false;
01625 if (p.lgEOL())
01626 p.NoNumb("limit to zone number");
01627
01628
01629 for (long i = 0; i < iterations.iter_malloc; i++)
01630 geometry.nend[i] = geometry.nEndDflt;
01631
01632 if (geometry.nEndDflt > 2000)
01633 fprintf(
01634 ioQQQ,
01635 "CAUTION - it will take a lot of memory to save"
01636 " results for %li zones. Is this many zones really necessary?\n",
01637 geometry.nEndDflt);
01638 }
01639 }
01640
01641 else if (p.nMatch("TSQD"))
01642 {
01643
01644
01645 peimbt.tsqden = (realnum) p.FFmtRead();
01646
01647 if (p.lgEOL())
01648 {
01649 p.NoNumb("highest density in t^2 section of printout");
01650 }
01651 peimbt.tsqden = (realnum) pow((realnum) 10.f, peimbt.tsqden);
01652 }
01653
01654 else if (p.nMatch("NMAP"))
01655 {
01656
01657 hcmap.nMapStep = (long) p.FFmtRead();
01658 if (p.lgEOL())
01659 {
01660 p.NoNumb("steps in heating-cooling map");
01661 }
01662 }
01663
01664 else if (p.nMatch("NUME") && p.nMatch("DERI"))
01665 {
01666
01667 NumDeriv.lgNumDeriv = true;
01668 }
01669
01670 else if (p.nMatch("PATH"))
01671 {
01672 fprintf(ioQQQ, " The SET PATH command is no longer supported.\n");
01673 fprintf(ioQQQ, " Please set the correct path in the file path.h.\n");
01674 fprintf(ioQQQ,
01675 " Or set the environment variable CLOUDY_DATA_PATH.\n Sorry.\n");
01676 cdEXIT(EXIT_FAILURE);
01677 }
01678
01679 else if (p.nMatch("PHFI"))
01680 {
01681
01682 ip = (long) p.FFmtRead();
01683 if (p.lgEOL())
01684 {
01685 p.NoNumb("version of PHFIT");
01686 }
01687
01688 if (ip == 1995)
01689 {
01690
01691 t_ADfA::Inst().set_version(PHFIT95);
01692 }
01693 else if (ip == 1996)
01694 {
01695
01696 t_ADfA::Inst().set_version(PHFIT96);
01697 }
01698 else
01699 {
01700 fprintf(ioQQQ, " Two possible values are 1995 and 1996.\n");
01701 cdEXIT(EXIT_FAILURE);
01702 }
01703 }
01704
01705
01706 else if (p.nMatch("PUNC") || p.nMatch("SAVE"))
01707 {
01708 if (p.nMatch("HASH"))
01709 {
01710
01711 if (!lgQuotesFound)
01712 {
01713 fprintf(ioQQQ,
01714 " I didn\'t recognize a key on this SET SAVE HASH line.\n");
01715 cdEXIT(EXIT_FAILURE);
01716 }
01717
01718
01719
01720
01721
01722
01723 if (strcmp(chString_quotes_lowercase, "return") == 0)
01724 {
01725
01726 sprintf(save.chHashString, "\n");
01727 }
01728 else if (strcmp(chString_quotes_lowercase, "time") == 0)
01729 {
01730
01731 sprintf(save.chHashString, "TIME_DEP");
01732 }
01733 else
01734 {
01735
01736 strcpy(save.chHashString, chString_quotes_lowercase);
01737 }
01738 }
01739
01740 else if ((p.nMatch("LINE") && p.nMatch("WIDT")) || p.nMatch("RESO")
01741 || p.nMatch("LWID"))
01742 {
01743
01744 if (p.nMatch(" C "))
01745 {
01746 fprintf(ioQQQ, " the keyword C is no longer necessary since"
01747 " energy conservation is now supported by default.\n");
01748 cdEXIT(EXIT_FAILURE);
01749 }
01750 else if (p.nMatch("SUPP"))
01751
01752 save.Resolution = realnum(0.);
01753 else
01754 {
01755 double number = p.FFmtRead();
01756 if (p.lgEOL())
01757 p.NoNumb("line width or resolution");
01758 if (number <= 0.)
01759 {
01760 fprintf(ioQQQ,
01761 " line width or resolution must be greater than zero.\n");
01762 cdEXIT(EXIT_FAILURE);
01763 }
01764
01765 if (p.nMatch("RESO"))
01766
01767 save.Resolution = realnum(number);
01768 else
01769
01770 save.Resolution = realnum(SPEEDLIGHT / (number * 1.e5));
01771 }
01772
01773
01774
01775 if (p.nMatch("ABSO"))
01776 save.ResolutionAbs = save.Resolution;
01777 }
01778
01779 else if (p.nMatch("PREF"))
01780 {
01781
01782 save.chFilenamePrefix = chString_quotes_lowercase;
01783 }
01784
01785 else if (p.nMatch("FLUS"))
01786 {
01787
01788 save.lgFLUSH = true;
01789 }
01790
01791 else
01792 {
01793 fprintf(ioQQQ,
01794 " There should have been an option on this command.\n");
01795 fprintf(ioQQQ,
01796 " Valid options for SET SAVE are summarized in Hazy 1 "
01797 "Miscellaneous commands.\n");
01798 fprintf(ioQQQ,
01799 " The SET PUNCHLWIDTH command is now SET SAVE LINE WIDTH.\n");
01800 cdEXIT(EXIT_FAILURE);
01801 }
01802 }
01803
01804
01805 else if (p.nMatch("CONT"))
01806 {
01807 if (p.nMatch("FEII") || p.nMatch("FE I"))
01808 {
01809
01810
01811 FeII.feconwlLo = (realnum) p.FFmtRead();
01812 FeII.feconwlHi = (realnum) p.FFmtRead();
01813 FeII.nfe2con = (long) p.FFmtRead();
01814 if (p.lgEOL())
01815 {
01816 fprintf(ioQQQ,
01817 "PROBLEM The set FeII continuum command must have"
01818 " three numbers, the lower and upper wavelength range in Angstroms"
01819 " and the number of bins to divide this into.\n");
01820 cdEXIT(EXIT_FAILURE);
01821 }
01822
01823 if (FeII.feconwlLo >= FeII.feconwlHi)
01824 {
01825 fprintf(ioQQQ, "PROBLEM The first two numbers on the set "
01826 "FeII continuum command must be the lower and upper "
01827 "wavelength range in Angstroms and the first must be less "
01828 "than the second.\n");
01829 cdEXIT(EXIT_FAILURE);
01830 }
01831 if (FeII.nfe2con < 2)
01832 {
01833 fprintf(ioQQQ,
01834 "PROBLEM The third number on the set FeII continuum "
01835 "command must be the number of bins to divide the range into and"
01836 " it must be greater than 1.\n");
01837 cdEXIT(EXIT_FAILURE);
01838 }
01839
01840 }
01841
01842 else if (p.nMatch("RESO"))
01843 {
01844
01845
01846 continuum.ResolutionScaleFactor = p.FFmtRead();
01847 if (p.lgEOL())
01848 {
01849 p.NoNumb("continuum resolution scale");
01850 }
01851
01852 if (continuum.ResolutionScaleFactor <= 0.)
01853 continuum.ResolutionScaleFactor = pow(10.,
01854 continuum.ResolutionScaleFactor);
01855 }
01856
01857 else if (p.nMatch("SHIE"))
01858 {
01859
01860
01861
01862
01863
01864
01865 if (p.nMatch("PESC"))
01866 {
01867
01868 rt.nLineContShield = LINE_CONT_SHIELD_PESC;
01869 }
01870 else if (p.nMatch("FEDE"))
01871 {
01872
01873
01874
01875
01876 rt.nLineContShield = LINE_CONT_SHIELD_FEDERMAN;
01877 }
01878 else if (p.nMatch("FERL"))
01879 {
01880 rt.nLineContShield = LINE_CONT_SHIELD_FERLAND;
01881 }
01882 else if (p.nMatch("NONE"))
01883 {
01884
01885 rt.nLineContShield = 0;
01886 }
01887 else
01888 {
01889 fprintf(ioQQQ,
01890 " I didn\'t recognize a key on this SET CONTINUUM SHIELDing line.\n");
01891 fprintf(ioQQQ,
01892 " The ones I know about are: PESC, FEDErman, & FERLand.\n");
01893 cdEXIT(EXIT_FAILURE);
01894 }
01895 }
01896
01897 else
01898 {
01899 fprintf(ioQQQ,
01900 " I didn\'t recognize a key on this SET CONTINUUM line.\n");
01901 fprintf(ioQQQ,
01902 " The ones I know about are: RESOlution and SHIEld.\n");
01903 cdEXIT(EXIT_FAILURE);
01904 }
01905 }
01906
01907 else
01908 {
01909 fprintf(ioQQQ, " I didn\'t recognize a key on this SET command.\n");
01910 cdEXIT(EXIT_FAILURE);
01911 }
01912
01913 return;
01914 }