00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "optimize.h"
00007 #include "doppvel.h"
00008 #include "stopcalc.h"
00009 #include "abund.h"
00010 #include "geometry.h"
00011 #include "dense.h"
00012 #include "plot.h"
00013 #include "grid.h"
00014 #include "rfield.h"
00015 #include "grainvar.h"
00016 #include "dynamics.h"
00017 #include "magnetic.h"
00018 #include "trace.h"
00019 #include "h2.h"
00020 #include "mole.h"
00021 #include "hmi.h"
00022 #include "rt.h"
00023 #include "thermal.h"
00024 #include "opacity.h"
00025 #include "atomfeii.h"
00026 #include "called.h"
00027 #include "wind.h"
00028 #include "hextra.h"
00029 #include "iterations.h"
00030 #include "radius.h"
00031 #include "input.h"
00032 #include "monitor_results.h"
00033 #include "phycon.h"
00034 #include "fudgec.h"
00035 #include "version.h"
00036 #include "conv.h"
00037 #include "parse.h"
00038 #include "cosmology.h"
00039 #include "pressure.h"
00040 #include "parser.h"
00041 #include "dark_matter.h"
00042
00043 void ParseAperture(Parser &p);
00044 void ParseAtom(Parser &p);
00045 void ParseBremsstrahlung(Parser &p);
00046 void ParseCExtra(Parser &p);
00047 void ParseCMBOuter(Parser &p);
00048 void ParseCosm(Parser &p);
00049 void ParseCovering(Parser &p);
00050 void ParseCylinder(Parser &p);
00051 void ParseDarkMatter(Parser &p);
00052 void ParseDielectronic(Parser &);
00053 void ParseDiffuse(Parser &p);
00054 void ParseDistance(Parser &p);
00055 void ParseDoubleTau(Parser &);
00056 void ParseEden(Parser &p);
00057 void ParseEnergy(Parser &p);
00058 void ParseFail(Parser &p);
00059 void ParseFill(Parser &p);
00060 void ParseF_nuSpecific(Parser &p);
00061 void ParseForceTemperature(Parser &p);
00062 void ParseFudge(Parser &p);
00063 void ParsePGrains(Parser &);
00064 void ParseGravity(Parser &p);
00065 void ParseHeLike(Parser &);
00066 void ParseHelp(Parser &);
00067 void ParseHExtra(Parser &p);
00068 void ParseConvHighT(Parser &);
00069 void ParseHydrogen(Parser &);
00070 void ParseInitCount(Parser &p);
00071 void ParseIntensity(Parser &p);
00072 void ParseIterations(Parser &p);
00073 void ParseL_nu(Parser &p);
00074 void ParseLaser(Parser &p);
00075 void ParseLuminosity(Parser &p);
00076 void ParseNeutrons(Parser &p);
00077 void ParseNuF_nu(Parser &p);
00078 void ParseNuL_nu(Parser &p);
00079 void ParsePhi(Parser &p);
00080 void ParseQH(Parser &p);
00081 void ParseRoberto(Parser &);
00082 void ParseSpecial(Parser &);
00083 void ParseTauMin(Parser &p);
00084 void ParseTitle(Parser &);
00085 void ParseTolerance(Parser &);
00086 void ParseVLaw(Parser &p);
00087 void ParseTurbulence(Parser &p);
00088
00089 void ParseCommands(void)
00090 {
00091 bool lgStop ,
00092 lgStop_not_enough_info;
00093
00094 DEBUG_ENTRY( "ParseCommands()" );
00095
00096
00097 abund.lgAbnSolar = true;
00098
00099
00100 plotCom.lgPlotON = false;
00101 plotCom.nplot = 0;
00102
00103
00104
00105
00106
00107
00108
00109 radius.Radius = 0.;
00110 radius.rinner = 0.;
00111 rfield.nShape = 0;
00112
00113
00114 input.lgSetNoBuffering = false;
00115
00116
00117 InitMonitorResults();
00118
00119 for( long int i=0; i < LIMSPC; i++ )
00120 {
00121 strcpy( rfield.chRSpec[i], "UNKN" );
00122 }
00123 optimize.nparm = 0;
00124
00125
00126
00127 if( optimize.lgTrOpt )
00128 {
00129
00130
00131 if( optimize.nTrOpt == optimize.nOptimiz )
00132 {
00133 trace.lgTrace = true;
00134 called.lgTalk = cpu.lgMPI_talk();
00135 trace.lgTrOvrd = true;
00136 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
00137 }
00138 }
00139
00140
00141 if( t_version::Inst().nBetaVer > 0 && called.lgTalk )
00142 {
00143 fprintf( ioQQQ,
00144 "\n This is a beta release of Cloudy, and is intended for testing only.\n\n" );
00145 }
00146
00147 if( called.lgTalk )
00148 {
00149
00150 int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2);
00151 fprintf( ioQQQ, "%*cCloudy %s\n", indent, ' ', t_version::Inst().chVersion );
00152
00153 fprintf( ioQQQ, "%57cwww.nublado.org\n\n", ' ' );
00154
00155
00156 fprintf( ioQQQ, "%23c", ' ' );
00157 fprintf( ioQQQ, "**************************************");
00158 fprintf( ioQQQ, "%7.7s", t_version::Inst().chDate);
00159 fprintf( ioQQQ, "**************************************\n");
00160
00161 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
00162 }
00163
00164
00165
00166
00167
00168
00169 input.iReadWay = 1;
00170
00171
00172
00173 input.init();
00174
00175 static const CloudyCommand commands[] = {
00176 {"ABSO",ParseAbsMag},
00177
00178 {"AGE",ParseAge},
00179
00180 {"AGN ",ParseAgn},
00181
00182 {"ABUN",ParseAbundancesNonSolar},
00183 {"APER",ParseAperture},
00184 {"ATOM",ParseAtom},
00185 {"BACK", ParseBackgrd},
00186
00187 {"BLAC", ParseBlackbody},
00188
00189 {"BREM", ParseBremsstrahlung},
00190 {"CASE",ParseCaseB},
00191
00192 {"CEXT",ParseCExtra},
00193
00194 {"CMB", ParseCMBOuter},
00195
00196
00197 {"COMP", ParseCompile},
00198
00199 {"CONS",ParseConstant},
00200
00201
00202 {"CORO",ParseCoronal},
00203
00204
00205 {"COSMOLOGY",ParseCosmology},
00206 {"COSMI", ParseCosmicRays},
00207 {"COSM",ParseCosm},
00208 {"COVE",ParseCovering},
00209 {"CRAS",ParseCrashDo},
00210
00211 {"CYLI",ParseCylinder},
00212 {"DARK",ParseDarkMatter},
00213 {"DIEL",ParseDielectronic},
00214 {"DIFF",ParseDiffuse},
00215 {"DIST",ParseDistance},
00216 {"DLAW",ParseDLaw},
00217
00218 {"DOUB",ParseDoubleTau},
00219
00220 {"DRIV",ParseDrive},
00221
00222 {"EDEN",ParseEden},
00223
00224
00225 {"ELEM",ParseElement},
00226
00227
00228
00229
00230 {"ENER",ParseEnergy},
00231 {"EXTI",ParseExtinguish},
00232
00233
00234
00235
00236
00237
00238
00239 {"FAIL",ParseFail},
00240
00241 {"FILL",ParseFill},
00242 {"FLUC",ParseFluc},
00243
00244 {"F(NU",ParseF_nuSpecific},
00245
00246
00247 {"FORC",ParseForceTemperature},
00248
00249
00250
00251 {"FUDG",ParseFudge},
00252 {"GLOB",ParseGlobule},
00253
00254
00255 {"GRAI",ParseGrain},
00256
00257 {"PGRA",ParsePGrains},
00258 {"GRAV",ParseGravity},
00259
00260 {"GRID",ParseGrid},
00261
00262
00263 {"HDEN",ParseHDEN},
00264
00265 {"HELI",ParseHeLike},
00266 {"HELP",ParseHelp},
00267 {"HEXT",ParseHExtra},
00268
00269
00270
00271
00272 {"HIGH",ParseConvHighT},
00273
00274 {"HYDROGEN",ParseHydrogen},
00275 {"ILLU",ParseIlluminate},
00276
00277 {"INIT",ParseInitCount},
00278 {"INTEN",ParseIntensity},
00279 {"INTER",ParseInterp},
00280
00281
00282
00283
00284
00285 {"IONI",ParseIonParI},
00286
00287
00288
00289 {"ITER",ParseIterations},
00290 {"L(NU",ParseL_nu},
00291 {"LASE",ParseLaser},
00292 {"LUMI",ParseLuminosity},
00293 {"MAGN", ParseMagnet},
00294
00295 {"MAP ",ParseMap},
00296
00297 {"META",ParseMetal},
00298
00299
00300 {"MONI",ParseMonitorResults},
00301
00302 {"NEUT",ParseNeutrons},
00303 {"NO ", ParseDont},
00304
00305 {"NORM",ParseNorm},
00306
00307 {"NUF(",ParseNuF_nu},
00308 {"NUL(",ParseNuL_nu},
00309 {"OPTI",ParseOptimize},
00310
00311
00312 {"PHI(",ParsePhi},
00313 {"PLOT", ParsePlot},
00314
00315
00316 {"POWE", ParsePowerlawContinuum},
00317
00318 {"PRIN",ParsePrint},
00319
00320 {"PUNC",ParseSave},
00321
00322 {"Q(H)",ParseQH},
00323 {"RATI",ParseRatio},
00324
00325
00326
00327
00328
00329
00330
00331 {"RADI", ParseRadius},
00332
00333
00334
00335
00336 {"ROBE",ParseRoberto},
00337 {"SAVE",ParseSave},
00338 {"SET ",ParseSet},
00339
00340 {"SPEC", ParseSpecial},
00341 {"SPHE", ParseSphere},
00342
00343
00344 {"STAT",ParseState},
00345
00346 {"STOP",ParseStop},
00347
00348
00349 {"TABL",ParseTable},
00350
00351
00352
00353
00354 {"TAUM",ParseTauMin},
00355 {"TEST",ParseTest},
00356
00357 {"TIME",ParseDynaTime},
00358
00359 {"TITL",ParseTitle},
00360 {"TLAW", ParseTLaw},
00361
00362 {"TOLE", ParseTolerance},
00363 {"TRAC", ParseTrace},
00364
00365 {"VLAW",ParseVLaw},
00366 {"TURB",ParseTurbulence},
00367 {"WIND",ParseDynaWind},
00368
00369
00370 {"XI",ParseIonParX},
00371 {NULL,NULL}};
00372
00373 Parser p(commands);
00374
00375 p.m_nqh = 0;
00376 p.m_nInitFile=0;
00377 p.m_lgDSet = false;
00378 p.m_lgEOF = false;
00379
00380
00381
00382 while (p.getline())
00383 {
00384
00385 if ( p.last( ) )
00386 break;
00387
00388
00389 p.echo( );
00390
00391
00392 if( p.nMatch("VARY") )
00393 {
00394 optimize.lgVarOn = true;
00395 if( optimize.nparm + 1 > LIMPAR )
00396 {
00397 fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n",
00398 LIMPAR );
00399 cdEXIT(EXIT_FAILURE);
00400 }
00401 }
00402 else
00403 {
00404 optimize.lgVarOn = false;
00405 }
00406
00407 if( p.isCommandComment() )
00408 {
00409 ((void)0);
00410 }
00411 else
00412 {
00413 long int i;
00414 for (i=0; commands[i].name != NULL; ++i)
00415 {
00416 if (p.Command(commands[i].name,commands[i].action))
00417 break;
00418 }
00419 if (commands[i].name == NULL)
00420 {
00421 p.CommandError();
00422 }
00423 }
00424 }
00425
00426
00427
00428
00429 if( called.lgTalk )
00430 {
00431 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' );
00432 fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n\n", ' ' );
00433 }
00434
00435
00436
00437
00438 if( optimize.lgTrOpt )
00439 {
00440
00441
00442 if( optimize.nTrOpt != optimize.nOptimiz )
00443 {
00444 trace.lgTrace = false;
00445
00446 trace.lgTrOvrd = false;
00447 }
00448 else
00449 {
00450 trace.lgTrace = true;
00451 called.lgTalk = cpu.lgMPI_talk();
00452 trace.lgTrOvrd = true;
00453 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" );
00454 }
00455 }
00456
00457
00458 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00459 {
00460 dense.gas_phase[ipHYDROGEN] = (realnum)dense_fabden(radius.Radius,radius.depth);
00461
00462 if( dense.gas_phase[ipHYDROGEN] <= 0. )
00463 {
00464 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
00465 cdEXIT(EXIT_FAILURE);
00466 }
00467 }
00468 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00469 {
00470 dense.gas_phase[ipHYDROGEN] = (realnum)dense_tabden(radius.Radius,radius.depth);
00471
00472 if( dense.gas_phase[ipHYDROGEN] <= 0. )
00473 {
00474 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
00475 cdEXIT(EXIT_FAILURE);
00476 }
00477 }
00478 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
00479 {
00480 dense.gas_phase[ipHYDROGEN] = (realnum)dense_parametric_wind(radius.Radius);
00481
00482 if( dense.gas_phase[ipHYDROGEN] <= 0. )
00483 {
00484 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" );
00485 cdEXIT(EXIT_FAILURE);
00486 }
00487 }
00488
00489
00490
00491
00492
00493 lgStop_not_enough_info = false;
00494 lgStop = false;
00495
00496
00497 if( dense.gas_phase[ipHYDROGEN] <= 0. )
00498 {
00499 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density MUST be specified.\n" );
00500 lgStop_not_enough_info = true;
00501 lgStop = true;
00502
00503
00504 dense.gas_phase[ipHYDROGEN] = 1.;
00505 }
00506
00507 if( geometry.covaper < 0.f || geometry.iEmissPower == 2 )
00508 geometry.covaper = geometry.covgeo;
00509
00510 radius.rinner = radius.Radius;
00511
00512
00513 wind.emdot = dense.gas_phase[ipHYDROGEN]*wind.windv0;
00514
00515
00516 if( iterations.lgConverge_set)
00517 {
00518 iterations.itermx = MIN2( iterations.itermx , iterations.lim_iter );
00519 for( long int j=0; j < iterations.iter_malloc; j++ )
00520 {
00521 geometry.nend[j] = MIN2( geometry.nend[j] , iterations.lim_zone );
00522 }
00523 }
00524
00525
00526 if( input.nSave < 0 )
00527 {
00528 fprintf( ioQQQ, " PROBLEM DISASTER No commands were entered - whats up?\n" );
00529 cdEXIT(EXIT_FAILURE);
00530 }
00531
00532
00533 if( wind.lgBallistic() && conv.lgAutoIt )
00534 {
00535 if( called.lgTalk )
00536 {
00537 fprintf( ioQQQ, " NOTE PROBLEM Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" );
00538 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off\n\n\n" );
00539 }
00540 conv.lgAutoIt = false;
00541 iterations.itermx = 0;
00542 }
00543
00544
00545
00546 if( opac.lgCaseB && conv.lgAutoIt && (wind.lgBallistic() || wind.lgStatic()) )
00547 {
00548 if( called.lgTalk )
00549 {
00550 fprintf( ioQQQ, " NOTE Case B is an artificial test, it makes no sense to converge this model.\n" );
00551 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off.\n\n\n" );
00552 }
00553 conv.lgAutoIt = false;
00554 iterations.itermx = 0;
00555 }
00556
00557
00558 if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 )
00559 {
00560 if( called.lgTalk )
00561 {
00562 fprintf( ioQQQ, " NOTE Specifying both a density power law and constant pressure is impossible.\n" );
00563 }
00564 lgStop = true;
00565 }
00566
00567 if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
00568 {
00569 if( called.lgTalk )
00570 {
00571 fprintf( ioQQQ, " NOTE NO REEVALUATE IONIZATION can only be used with constant density.\n" );
00572 fprintf( ioQQQ, " NOTE Resetting to reevaluate ionization.\n\n" );
00573 }
00574 rfield.lgIonizReevaluate = true;
00575 }
00576
00577 if( !rfield.lgOpacityReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 )
00578 {
00579 if( called.lgTalk )
00580 {
00581 fprintf( ioQQQ, " NOTE NO REEVALUATE OPACITY can only be used with constant density.\n" );
00582 fprintf( ioQQQ, " NOTE Resetting to reevaluate opacity.\n\n" );
00583 }
00584 rfield.lgOpacityReevaluate = true;
00585 }
00586
00587
00588 if( pressure.external_mass[0].size()>0 && pressure.gravity_symmetry==-1 )
00589 {
00590 if( called.lgTalk )
00591 {
00592 fprintf( ioQQQ, " NOTE Gravity from an external mass has been added, but no symmetry (spherical/mid-plane) was specified.\n" );
00593 fprintf( ioQQQ, " NOTE It will be ignored.\n\n\n" );
00594 }
00595 }
00596
00597
00598 double thickness = min( MIN3( StopCalc.tauend, StopCalc.colpls, StopCalc.colnut ),
00599 MIN3( StopCalc.col_h2, StopCalc.col_h2_nut, StopCalc.HColStop ) );
00600 if( thickness < COLUMN_INIT )
00601 {
00602
00603 thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
00604
00605 if( thickness > 1e25 && radius.StopThickness[0] > 1e25 )
00606 {
00607 fprintf( ioQQQ,
00608 "NOTE The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n",
00609 thickness);
00610 fprintf( ioQQQ,
00611 "NOTE This seems large to me.\n");
00612 fprintf(ioQQQ,"NOTE a very large radius may cause overflow.\n\n");
00613 }
00614 }
00615
00616 if( gv.lgDColOn && thermal.ConstGrainTemp>0 && called.lgTalk )
00617 {
00618
00619
00620 fprintf( ioQQQ,
00621 "NOTE The grain temperatures are set to a constant value with the "
00622 "CONSTANT GRAIN TEMPERATURE command, but "
00623 "energy exchange \n");
00624 fprintf( ioQQQ,
00625 "NOTE is still included. The grain-gas heating-cooling will be incorrect. "
00626 "Consider turning off gas-grain collisional energy\n");
00627 fprintf( ioQQQ,
00628 "NOTE exchange with the NO GRAIN GAS COLLISIONAL ENERGY EXCHANGE command.\n\n\n");
00629 }
00630
00631 if( !rfield.lgDoLineTrans && rfield.lgOpacityFine )
00632 {
00633 if( called.lgTalk )
00634 {
00635 fprintf( ioQQQ, " NOTE NO LINE TRANSER set but fine opacities still computed.\n" );
00636 fprintf( ioQQQ, " NOTE Turning off fine opacities.\n\n" );
00637 }
00638 rfield.lgOpacityFine = false;
00639 }
00640
00641 if( h2.lgH2ON && (!rfield.lgDoLineTrans || !rfield.lgOpacityFine) )
00642 {
00643 if( called.lgTalk )
00644 {
00645 fprintf( ioQQQ, " NOTE Large H2 molecule turned on but line transfer and fine opacities are not.\n" );
00646 fprintf( ioQQQ, " NOTE Turning on line transfer and fine opacities.\n\n" );
00647 }
00648 rfield.lgOpacityFine = true;
00649 rfield.lgDoLineTrans = true;
00650 }
00651
00652 if( rfield.lgMustBlockHIon && !rfield.lgBlockHIon )
00653 {
00654
00655
00656 fprintf( ioQQQ, "\n NOTE\n"
00657 " NOTE One of the incident continuum is a form used when no H-ionizing radiation is produced.\n" );
00658 fprintf( ioQQQ, " NOTE You must also include the EXTINGUISH command to make sure this is done.\n" );
00659 fprintf( ioQQQ, " NOTE The EXTINGUISH command was not included.\n" );
00660 fprintf( ioQQQ, " NOTE YOU MAY BE MAKING A BIG MISTAKE!!\n NOTE\n\n\n\n" );
00661 }
00662
00663
00664
00665 if( called.lgTalk && (StopCalc.TempLoStopZone < phycon.TEMP_STOP_DEFAULT ||
00666
00667 (thermal.ConstTemp > 0. && thermal.ConstTemp < phycon.TEMP_STOP_DEFAULT ) ) )
00668 {
00669
00670
00671 if( (hextra.cryden == 0.) && !co.lgNoCOMole && !hmi.lgNoH2Mole )
00672 {
00673 fprintf( ioQQQ, "\n NOTE\n"
00674 " NOTE The simulation is going into neutral gas but cosmic rays are not included.\n" );
00675 fprintf( ioQQQ, " NOTE Ion-molecule chemistry will not occur without a source of ionization.\n" );
00676 fprintf( ioQQQ, " NOTE The chemistry network may collapse deep in molecular regions.\n" );
00677 fprintf( ioQQQ, " NOTE Consider adding galactic background cosmic rays with the COSMIC RAYS BACKGROUND command.\n" );
00678 fprintf( ioQQQ, " NOTE You may be making a BIG mistake.\n NOTE\n\n\n\n" );
00679 }
00680 }
00681
00682
00683
00684 if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 )
00685 {
00686 fprintf( ioQQQ, " NOTE Is the entered value of the hydrogen density (%.2e) reasonable?\n",
00687 dense.gas_phase[ipHYDROGEN]);
00688 fprintf( ioQQQ, " NOTE It seems pretty low to me.\n\n\n" );
00689 }
00690 else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 )
00691 {
00692 fprintf( ioQQQ, " NOTE Is this value of the hydrogen density reasonable?\n" );
00693 fprintf( ioQQQ, " NOTE It seems pretty high to me.\n\n\n" );
00694 }
00695
00696
00697 if( called.lgTalk && !lgStop && !lgStop_not_enough_info )
00698 {
00699 if( dense.gas_phase[ipHYDROGEN] < 1e-6 || dense.gas_phase[ipHYDROGEN] > 1e19 )
00700 {
00701 fprintf( ioQQQ, " NOTE Simulation may crash because of extreme "
00702 "density. The value was %.2e\n\n" ,
00703 dense.gas_phase[ipHYDROGEN] );
00704 }
00705 }
00706
00707 if( rfield.nShape == 0 )
00708 {
00709 fprintf( ioQQQ, " PROBLEM DISASTER No incident radiation field was specified - "
00710 "at least put in the CMB.\n" );
00711 lgStop = true;
00712 lgStop_not_enough_info = true;
00713 }
00714
00715 if( (p.m_nqh) == 0 )
00716 {
00717 fprintf( ioQQQ, " PROBLEM DISASTER Luminosity of continuum MUST be specified.\n" );
00718 lgStop = true;
00719 lgStop_not_enough_info = true;
00720 }
00721
00722
00723
00724
00725 if( radius.Radius == 0. && rfield.nShape> 0)
00726 {
00727 fprintf( ioQQQ, " PROBLEM DISASTER Starting radius MUST be specified.\n" );
00728 lgStop = true;
00729 lgStop_not_enough_info = true;
00730 }
00731
00732 if( rfield.nShape != (p.m_nqh) )
00733 {
00734 fprintf( ioQQQ, " PROBLEM DISASTER There were not the same number of continuum shapes and luminosities entered.\n" );
00735 lgStop = true;
00736 }
00737
00738
00739
00740
00741 static bool lgFirstPass = true;
00742
00743
00744
00745 if( optimize.lgNoVary && grid.lgGrid )
00746 {
00747
00748 grid.lgGrid = false;
00749 optimize.nparm = 0;
00750 grid.nGridCommands = 0;
00751 }
00752
00753 if( lgFirstPass && grid.lgGrid && (optimize.nparm!=grid.nGridCommands) )
00754 {
00755
00756 fprintf( ioQQQ, " PROBLEM DISASTER The GRID command was entered "
00757 "but there were %li GRID commands and %li commands with a VARY option.\n" ,
00758 grid.nGridCommands , optimize.nparm);
00759 fprintf( ioQQQ, " There must be the same number of GRIDs and VARY.\n" );
00760 lgStop = true;
00761 }
00762 lgFirstPass = false;
00763
00764 if( lgStop_not_enough_info )
00765 {
00766 fprintf( ioQQQ, " PROBLEM DISASTER I do not have enough information to do the simulation, I cannot go on.\n" );
00767 fprintf( ioQQQ, "\n\n Sorry.\n\n\n" );
00768 cdEXIT(EXIT_FAILURE);
00769 }
00770
00771 {
00772 bool lgParserTest = false;
00773 if ( lgParserTest )
00774 {
00775
00776
00777 fprintf(ioQQQ,"Parser phase PASSED\n");
00778 exit(0);
00779 }
00780 }
00781
00782 if( lgStop )
00783 {
00784 cdEXIT(EXIT_FAILURE);
00785 }
00786
00787
00788 return;
00789 }
00790
00791 void ParseAbundancesNonSolar(Parser &p)
00792 {
00793
00794
00795 ParseAbundances(p);
00796
00797 abund.lgAbnSolar = false;
00798 }
00799
00800 void ParseAperture(Parser &p)
00801 {
00802 DEBUG_ENTRY("ParseAperture()");
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814 if( p.nMatch("SLIT") )
00815 {
00816
00817
00818 geometry.iEmissPower = 1;
00819 }
00820 else if( p.nMatch("BEAM") )
00821 {
00822
00823
00824 geometry.iEmissPower = 0;
00825 }
00826 else if( p.nMatch("SIZE") )
00827 {
00828
00829
00830 geometry.size = realnum(p.FFmtRead());
00831 if( p.lgEOL() )
00832 {
00833 p.NoNumb("aperture size");
00834 }
00835 if( geometry.size <= 0.f )
00836 {
00837 fprintf( ioQQQ, " The aperture size must be positive. Sorry.\n" );
00838 cdEXIT(EXIT_FAILURE);
00839 }
00840 geometry.lgSizeSet = true;
00841 }
00842 else if( p.nMatch("COVE") )
00843 {
00844
00845
00846 geometry.covaper = realnum(p.FFmtRead());
00847 if( p.lgEOL() )
00848 {
00849 p.NoNumb("aperture covering factor");
00850 }
00851 if( geometry.covaper <= 0.f || geometry.covaper > 1.f )
00852 {
00853 fprintf( ioQQQ, " The aperture covering factor must be > 0 and <= 1. Sorry.\n" );
00854 cdEXIT(EXIT_FAILURE);
00855 }
00856 }
00857 else
00858 {
00859 fprintf( ioQQQ, " One of the keywords SLIT, BEAM, SIZE or COVEring factor must appear.\n" );
00860 fprintf( ioQQQ, " Sorry.\n" );
00861 cdEXIT(EXIT_FAILURE);
00862 }
00863 }
00864 void ParseAtom(Parser &p)
00865 {
00866 DEBUG_ENTRY("ParseAtom()");
00867
00868 if( p.nMatch("FEII") || p.nMatch("FE II") )
00869 {
00870
00871 ParseAtomFeII(p);
00872 }
00873
00874 else if( p.nMatch("H-LI") )
00875 {
00876
00877 ParseAtomISO(ipH_LIKE, p);
00878 }
00879
00880 else if( p.nMatch("HE-L") )
00881 {
00882
00883 ParseAtomISO(ipHE_LIKE, p);
00884 }
00885
00886 else if( p.nMatch("ROTO") || p.nMatch(" CO ") )
00887 {
00888
00889 fprintf(ioQQQ," The old CO models no longer exist, and this command is no longer supported.\n" );
00890 fprintf( ioQQQ, " Sorry.\n" );
00891 cdEXIT(EXIT_FAILURE);
00892 }
00893
00894 else if( p.nMatch(" H2 ") )
00895 {
00896
00897
00898 ParseAtomH2(p);
00899 }
00900
00901 else
00902 {
00903 fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n");
00904 fprintf( ioQQQ, " The available keys are FeII, H-Like, He-like, rotor and H2.\n");
00905 fprintf( ioQQQ, " Sorry.\n" );
00906 cdEXIT(EXIT_FAILURE);
00907 }
00908 }
00909 void ParseBremsstrahlung(Parser &p)
00910 {
00911 DEBUG_ENTRY("ParseBremsstrahlung()");
00912
00913 strcpy( rfield.chSpType[rfield.nShape], "BREMS" );
00914 rfield.slope[rfield.nShape] =
00915 (realnum)p.FFmtRead();
00916 if( p.lgEOL() )
00917 {
00918 p.NoNumb("temperature");
00919 }
00920
00921
00922 if( rfield.slope[rfield.nShape] <= 10. || p.nMatch(" LOG") )
00923 {
00924 rfield.slope[rfield.nShape] =
00925 pow(10.,rfield.slope[rfield.nShape]);
00926 }
00927 rfield.cutoff[rfield.nShape][0] = 0.;
00928
00929
00930 if( optimize.lgVarOn )
00931 {
00932
00933 optimize.nvarxt[optimize.nparm] = 1;
00934 strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f LOG" );
00935
00936 optimize.nvfpnt[optimize.nparm] = input.nRead;
00937
00938 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nShape]);
00939 optimize.vincr[optimize.nparm] = 0.5;
00940 ++optimize.nparm;
00941 }
00942 ++rfield.nShape;
00943 if( rfield.nShape >= LIMSPC )
00944 {
00945
00946 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
00947 cdEXIT(EXIT_FAILURE);
00948 }
00949 }
00950 void ParseCExtra(Parser &p)
00951 {
00952
00953 thermal.lgCExtraOn = true;
00954 thermal.CoolExtra = (realnum)pow(10.,p.FFmtRead());
00955 if( p.lgEOL() )
00956 {
00957 p.NoNumb("extra cooling");
00958 }
00959 thermal.cextpw = (realnum)p.FFmtRead();
00960 }
00961 void ParseCMBOuter(Parser &p)
00962 {
00963 cosmology.redshift_start = (realnum)p.FFmtRead();
00964 cosmology.redshift_current = cosmology.redshift_start;
00965
00966
00967 if( p.nMatch( "TIME" ) )
00968 rfield.lgTimeVary[(p.m_nqh)] = true;
00969
00970 ParseCMB(cosmology.redshift_start,&(p.m_nqh));
00971
00972
00973 if( p.nMatch("DENS") && !p.lgEOL() )
00974 {
00975 char chStuff[INPUT_LINE_LENGTH];
00976
00977
00978 sprintf( chStuff , "HDEN %.2e LINEAR", GetDensity( cosmology.redshift_start ) );
00979 p.setline(chStuff);
00980 p.set_point(4);
00981 ParseHDEN(p);
00982 }
00983 }
00984 void ParseCosm(Parser &)
00985 {
00986 DEBUG_ENTRY("ParseCosm()");
00987 fprintf(ioQQQ,
00988 "This command is now ambiguous -- please specify either COSMIC RAYS or COSMOLOGY.\nSorry.\n");
00989 cdEXIT(EXIT_FAILURE);
00990 }
00991 void ParseCovering(Parser &p)
00992 {
00993 DEBUG_ENTRY("ParseCovering()");
00994
00995
00996
00997 geometry.covgeo = (realnum)p.FFmtRead();
00998
00999 if( p.lgEOL() )
01000 {
01001 p.NoNumb("covering factor");
01002 }
01003
01004
01005 if( geometry.covgeo <= 0. )
01006 {
01007 geometry.covgeo = (realnum)pow((realnum)10.f,geometry.covgeo);
01008 }
01009
01010 if( geometry.covgeo > 1. )
01011 {
01012 fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense. Sorry.\n" );
01013 cdEXIT(EXIT_FAILURE);
01014 }
01015
01016
01017 geometry.covrt = geometry.covgeo;
01018 }
01019 void ParseCylinder(Parser &p)
01020 {
01021
01022 radius.lgCylnOn = true;
01023 radius.CylindHigh = pow(10.,p.FFmtRead());
01024 if( p.lgEOL() )
01025 {
01026 p.NoNumb("height");
01027 }
01028 }
01029 void ParseDarkMatter(Parser &p)
01030 {
01031 DEBUG_ENTRY( "ParseDarkMatter()" );
01032
01033 if( p.nMatch(" NFW") )
01034 {
01035
01036
01037
01038 dark.r_200 = (realnum)p.getNumberCheckAlwaysLog("NFW r_200");
01039
01040 dark.r_s = (realnum)p.getNumberDefaultAlwaysLog("NFW r_s", log10(dark.r_200)-1. );
01041 dark.lgNFW_Set = true;
01042
01043
01044 if( optimize.lgVarOn )
01045 {
01046
01047 optimize.nvarxt[optimize.nparm] = 1;
01048 strcpy( optimize.chVarFmt[optimize.nparm], "DARK NFW %f" );
01049
01050 optimize.nvfpnt[optimize.nparm] = input.nRead;
01051
01052 optimize.vparm[0][optimize.nparm] = (realnum)log10(dark.r_200);
01053 optimize.vincr[optimize.nparm] = 0.5;
01054 ++optimize.nparm;
01055 }
01056 }
01057 else
01058 {
01059 fprintf( ioQQQ, " Did not recognize a valid option for this DARK command.\nSorry.\n\n" );
01060 cdEXIT(EXIT_FAILURE);
01061 }
01062
01063 return;
01064 }
01065 void ParseDielectronic(Parser &)
01066 {
01067 DEBUG_ENTRY("ParseDielectronic()");
01068
01069 fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" );
01070 fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\n" );
01071 cdEXIT(EXIT_FAILURE);
01072 }
01073 void ParseDiffuse(Parser &p)
01074 {
01075 DEBUG_ENTRY("ParseDiffuse()");
01076
01077
01078 if( p.nMatch(" OTS") )
01079 {
01080 if( p.nMatch("SIMP") )
01081 {
01082
01083 strcpy( rfield.chDffTrns, "OSS" );
01084 }
01085 else
01086 {
01087
01088 strcpy( rfield.chDffTrns, "OTS" );
01089 }
01090 rfield.lgOutOnly = false;
01091 }
01092 else if( p.nMatch(" OUT") )
01093 {
01094 rfield.lgOutOnly = true;
01095 long int j = (long int)p.FFmtRead();
01096 if( p.lgEOL() )
01097 {
01098
01099 strcpy( rfield.chDffTrns, "OU2" );
01100 }
01101 else
01102 {
01103 if( j > 0 && j < 10 )
01104 {
01105 sprintf( rfield.chDffTrns, "OU%1ld", j );
01106 }
01107 else
01108 {
01109 fprintf( ioQQQ, " must be between 1 and 9 \n" );
01110 cdEXIT(EXIT_FAILURE);
01111 }
01112 }
01113 }
01114
01115 else
01116 {
01117 fprintf( ioQQQ, " There should have been OUTward or OTS on this line. Sorry.\n" );
01118 cdEXIT(EXIT_FAILURE);
01119 }
01120 }
01121 void ParseDistance(Parser &p)
01122 {
01123
01124 radius.distance = p.FFmtRead();
01125 if( p.lgEOL() )
01126 {
01127 p.NoNumb("distance");
01128 }
01129
01130
01131 if( !p.nMatch("LINE" ) )
01132 {
01133 radius.distance = pow(10., radius.distance );
01134 }
01135
01136
01137 if( p.nMatch("PARS") )
01138 {
01139 radius.distance *= PARSEC;
01140 }
01141
01142
01143 if( optimize.lgVarOn )
01144 {
01145 strcpy( optimize.chVarFmt[optimize.nparm], "DISTANCE %f LOG" );
01146 optimize.nvfpnt[optimize.nparm] = input.nRead;
01147 optimize.vparm[0][optimize.nparm] = realnum(log10(radius.distance));
01148 optimize.vincr[optimize.nparm] = 0.3f;
01149 optimize.nvarxt[optimize.nparm] = 1;
01150 ++optimize.nparm;
01151 }
01152 }
01153 void ParseDoubleTau(Parser &)
01154 {
01155 rt.DoubleTau = 2.;
01156 }
01157 void ParseEden(Parser &p)
01158 {
01159 dense.EdenExtra = (realnum)pow(10.,p.FFmtRead());
01160 if( p.lgEOL() )
01161 {
01162 p.NoNumb("electron density");
01163 }
01164
01165 phycon.lgPhysOK = false;
01166 }
01167 void ParseEnergy(Parser &p)
01168 {
01169 DEBUG_ENTRY("ParseEnergy()");
01170
01171 if( (p.m_nqh) >= LIMSPC )
01172 {
01173
01174 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01175 cdEXIT(EXIT_FAILURE);
01176 }
01177
01178 ASSERT( (p.m_nqh) < LIMSPC );
01179 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
01180 realnum teset = (realnum)p.FFmtRead();
01181 if( p.lgEOL() )
01182 {
01183 p.NoNumb("energy density");
01184 }
01185
01186
01187 if( !radius.lgRadiusKnown )
01188 {
01189
01190 radius.Radius = pow(10.,radius.rdfalt);
01191 }
01192
01193
01194 if( !p.nMatch(" LOG") && (p.nMatch("LINE") || teset > 10.) )
01195 {
01196
01197 teset = (realnum)log10(teset);
01198 }
01199
01200 if( teset > 5. )
01201 {
01202 fprintf( ioQQQ, " This intensity may be too large. The code may crash due to overflow. Was log intended?\n" );
01203 }
01204
01205
01206 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
01207
01208
01209 rfield.range[(p.m_nqh)][0] = rfield.emm;
01210 rfield.range[(p.m_nqh)][1] = rfield.egamry;
01211 rfield.totpow[(p.m_nqh)] = (4.*(teset) - 4.2464476 + 0.60206);
01212
01213
01214 if( p.nMatch( "TIME" ) )
01215 rfield.lgTimeVary[(p.m_nqh)] = true;
01216
01217
01218 if( optimize.lgVarOn )
01219 {
01220 strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f LOG" );
01221 if( rfield.lgTimeVary[(p.m_nqh)] )
01222 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01223
01224 optimize.nvfpnt[optimize.nparm] = input.nRead;
01225 optimize.vparm[0][optimize.nparm] = teset;
01226 optimize.vincr[optimize.nparm] = 0.1f;
01227 optimize.nvarxt[optimize.nparm] = 1;
01228 ++optimize.nparm;
01229 }
01230
01231
01232 ++(p.m_nqh);
01233 }
01234 void ParseFail(Parser &p)
01235 {
01236
01237 long int j = conv.LimFail;
01238 conv.LimFail = (long int)p.FFmtRead();
01239 if( p.lgEOL() )
01240 {
01241 p.NoNumb("limit");
01242 }
01243
01244
01245
01246
01247 if( p.nMatch(" MAP") && !p.nMatch(" NO ") )
01248 {
01249 conv.lgMap = true;
01250 }
01251
01252
01253 if( conv.LimFail > j )
01254 {
01255 fprintf( ioQQQ, " This command should not be necessary.\n" );
01256 fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\n" );
01257 }
01258 }
01259 void ParseFill(Parser &p)
01260 {
01261
01262 realnum a = (realnum)p.FFmtRead();
01263 if( p.lgEOL() )
01264 {
01265 p.NoNumb("filling factor");
01266 }
01267
01268 if( a <= 0. || p.nMatch(" LOG") )
01269 {
01270
01271 geometry.FillFac = (realnum)pow((realnum)10.f,a);
01272 }
01273 else
01274 {
01275
01276 geometry.FillFac = a;
01277 }
01278 if( geometry.FillFac > 1.0 )
01279 {
01280 if( called.lgTalk )
01281 fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" );
01282 geometry.FillFac = 1.;
01283 }
01284 geometry.fiscal = geometry.FillFac;
01285
01286
01287 geometry.filpow = (realnum)p.FFmtRead();
01288
01289
01290 if( optimize.lgVarOn )
01291 {
01292 strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f LOG power= %f" );
01293
01294
01295 optimize.nvfpnt[optimize.nparm] = input.nRead;
01296 optimize.vparm[0][optimize.nparm] = (realnum)log10(geometry.FillFac);
01297 optimize.vincr[optimize.nparm] = 0.5f;
01298
01299
01300 optimize.vparm[1][optimize.nparm] = geometry.filpow;
01301 optimize.nvarxt[optimize.nparm] = 2;
01302
01303
01304 optimize.varang[optimize.nparm][0] = -FLT_MAX;
01305 optimize.varang[optimize.nparm][1] = 0.f;
01306 ++optimize.nparm;
01307 }
01308 }
01309 void ParseF_nuSpecific(Parser &p)
01310 {
01311 bool lgNu2 = false;
01312
01313 ParseF_nu(p,"SQCM",lgNu2);
01314 }
01315 void ParseForceTemperature(Parser &p)
01316 {
01317 thermal.ConstTemp = (realnum)p.FFmtRead();
01318 if( p.lgEOL() )
01319 {
01320 p.NoNumb("temperature");
01321 }
01322
01323 if( p.nMatch(" LOG") || (thermal.ConstTemp <= 10. &&
01324 !p.nMatch("LINE")) )
01325 {
01326 thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp);
01327 }
01328
01329
01330 if( thermal.ConstTemp < 3. )
01331 {
01332 fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" );
01333 thermal.ConstTemp = 3.;
01334 }
01335 }
01336 void ParseFudge(Parser &p)
01337 {
01338
01339 fudgec.nfudge = 0;
01340 for( long int j=0; j < NFUDGC; j++ )
01341 {
01342 fudgec.fudgea[j] = (realnum)p.FFmtRead();
01343
01344 if( !p.lgEOL() )
01345 fudgec.nfudge = j+1;
01346 }
01347 if( fudgec.nfudge == 0 )
01348 p.NoNumb("fudge factor");
01349
01350
01351 if( optimize.lgVarOn )
01352 {
01353
01354 optimize.nvarxt[optimize.nparm] = 1;
01355 strcpy( optimize.chVarFmt[optimize.nparm], "FUDGE= %f" );
01356
01357 char chHold[1000];
01358 for( long int j=1; j<fudgec.nfudge; ++j )
01359 {
01360 sprintf( chHold , " %f" , fudgec.fudgea[j] );
01361 strcat( optimize.chVarFmt[optimize.nparm] , chHold );
01362 }
01363 optimize.lgOptimizeAsLinear[optimize.nparm] = true;
01364
01365
01366 optimize.nvfpnt[optimize.nparm] = input.nRead;
01367
01368 optimize.vparm[0][optimize.nparm] = fudgec.fudgea[0];
01369
01370 optimize.vincr[optimize.nparm] = abs(0.2f*fudgec.fudgea[0]);
01371
01372 if( optimize.vincr[optimize.nparm] == 0.f )
01373 optimize.vincr[optimize.nparm] = 1.f;
01374 ++optimize.nparm;
01375 }
01376 }
01377 void ParsePGrains(Parser &)
01378 {
01379 DEBUG_ENTRY("ParsePGrains()");
01380 fprintf(ioQQQ," Sorry, this command is obsolete, you can now use the normal GRAINS command.\n");
01381 cdEXIT(EXIT_FAILURE);
01382 }
01383 void ParseGravity(Parser &p)
01384 {
01385 DEBUG_ENTRY("ParseGravity()");
01386
01387 if( p.nMatch("EXTE") )
01388 {
01389
01390
01391
01392 double M_i = p.FFmtRead();
01393
01394 if( !p.lgEOL() && p.nMatch("LOG") )
01395 M_i = pow( 10., M_i );
01396 pressure.external_mass[0].push_back( M_i );
01397
01398
01399
01400 double x_i = p.FFmtRead();
01401
01402 if( !p.lgEOL() && p.nMatch("LOG") )
01403 x_i = pow( 10., x_i );
01404 pressure.external_mass[1].push_back( x_i * PARSEC );
01405
01406
01407 pressure.external_mass[2].push_back( p.FFmtRead() );
01408 }
01409 else
01410 {
01411
01412 if( p.nMatch("SPHE") )
01413 {
01414 pressure.gravity_symmetry = 0;
01415 }
01416 else if( p.nMatch("PLAN") )
01417 {
01418 pressure.gravity_symmetry = 1;
01419 }
01420 else
01421 {
01422 fprintf( ioQQQ, " The symmetry of the gravitational mass must be specified explicitly. Sorry.\n" );
01423 cdEXIT(EXIT_FAILURE);
01424 }
01425
01426
01427
01428 pressure.self_mass_factor = p.FFmtRead();
01429
01430 if( p.lgEOL() )
01431 pressure.self_mass_factor = 1.;
01432 else if( p.nMatch("LOG") )
01433 {
01434 pressure.self_mass_factor = pow( 10., pressure.self_mass_factor );
01435 }
01436 }
01437 }
01438 void ParseHeLike(Parser &)
01439 {
01440 DEBUG_ENTRY("ParseHeLike()");
01441 fprintf(ioQQQ,"Sorry, this command is replaced with ATOM HE-LIKE\n");
01442 cdEXIT(EXIT_FAILURE);
01443 }
01444 void ParseHelp(Parser &p)
01445 {
01446 DEBUG_ENTRY("ParseHelp()");
01447 p.help(ioQQQ);
01448 }
01449 void ParseHExtra(Parser &p)
01450 {
01451 hextra.TurbHeat = (realnum)pow(10.,p.FFmtRead());
01452 if( p.lgEOL() )
01453 p.NoNumb("extra heating" );
01454
01455 hextra.TurbHeatSave = hextra.TurbHeat;
01456
01457
01458 const char *chHextraScale;
01459
01460 realnum par1=0. , par2=0.;
01461 long int nPar;
01462 if( p.nMatch("DEPT") )
01463 {
01464
01465 hextra.lgHextraDepth = true;
01466 chHextraScale = "DEPTH";
01467
01468 realnum a = (realnum)p.FFmtRead();
01469 if( p.lgEOL() )
01470 p.NoNumb("depth");
01471 else
01472 hextra.turrad = (realnum)pow((realnum)10.f,a);
01473
01474
01475
01476 realnum b = (realnum)p.FFmtRead();
01477 if( p.lgEOL() )
01478 {
01479 hextra.turback = 0.;
01480 nPar = 2;
01481 }
01482 else
01483 {
01484 hextra.turback = (realnum)pow((realnum)10.f,b);
01485 nPar = 3;
01486 }
01487 par1 = a;
01488 par2 = b;
01489 }
01490 else if( p.nMatch("DENS") )
01491 {
01492
01493 chHextraScale = "DENSITY";
01494 hextra.lgHextraDensity = true;
01495
01496
01497 realnum a = (realnum)p.FFmtRead();
01498
01499
01500 hextra.HextraScaleDensity = (realnum)pow((realnum)10.f,a);
01501 par1 = a;
01502 par2 = 0;
01503 nPar = 2;
01504 }
01505 else if( p.nMatch("SS") )
01506 {
01507
01508 chHextraScale = "SS";
01509 hextra.lgHextraSS = true;
01510
01511
01512 hextra.HextraSSalpha = hextra.TurbHeat;
01513
01514
01515 realnum a = (realnum)p.FFmtRead();
01516 hextra.HextraSS_M = a*SOLAR_MASS;
01517
01518
01519 realnum b = (realnum)p.FFmtRead();
01520 hextra.HextraSSradius = b;
01521 par1 = a;
01522 par2 = 0;
01523 nPar = 2;
01524 }
01525 else
01526 {
01527
01528 chHextraScale = "";
01529 nPar = 1;
01530 }
01531
01532
01533 if( p.nMatch( "TIME" ) )
01534 hextra.lgTurbHeatVaryTime = true;
01535
01536
01537 if( optimize.lgVarOn )
01538 {
01539
01540 optimize.nvarxt[optimize.nparm] = nPar;
01541 optimize.vparm[0][optimize.nparm] = log10(hextra.TurbHeat);
01542 optimize.vparm[1][optimize.nparm] = par1;
01543 optimize.vparm[2][optimize.nparm] = par2;
01544
01545
01546 strcpy( optimize.chVarFmt[optimize.nparm], "HEXTra %f LOG " );
01547 strcat( optimize.chVarFmt[optimize.nparm], chHextraScale );
01548 while( nPar > 1 )
01549 {
01550
01551 --nPar;
01552 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
01553 }
01554 if( hextra.lgTurbHeatVaryTime )
01555 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01556
01557 optimize.nvfpnt[optimize.nparm] = input.nRead;
01558
01559 optimize.vincr[optimize.nparm] = 0.1f;
01560 ++optimize.nparm;
01561 }
01562 }
01563
01564 void ParseConvHighT(Parser &)
01565 {
01566 thermal.lgTeHigh = true;
01567 }
01568 void ParseHydrogen(Parser &)
01569 {
01570 DEBUG_ENTRY("ParseHydrogen()");
01571 fprintf(ioQQQ," Sorry, this command has been replaced with the ATOM H-LIKE command.\n");
01572 cdEXIT(EXIT_FAILURE);
01573 }
01574 void ParseInitCount(Parser &p)
01575 {
01576 DEBUG_ENTRY("ParseInitCount()");
01577
01578
01579
01580 ParseInit(p);
01581
01582
01583
01584 ++p.m_nInitFile;
01585 if( p.m_nInitFile > 1 )
01586 {
01587 fprintf( ioQQQ,
01588 " This is the second init file, I can only handle one.\nSorry.\n" );
01589 cdEXIT(EXIT_FAILURE);
01590 }
01591
01592
01593
01594 input.iReadWay = -1;
01595
01596
01597
01598 input.init();
01599 }
01600 void ParseIntensity(Parser &p)
01601 {
01602 DEBUG_ENTRY("ParseIntensity()");
01603
01604 if( (p.m_nqh) >= LIMSPC )
01605 {
01606
01607 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01608 cdEXIT(EXIT_FAILURE);
01609 }
01610
01611
01612 ASSERT( (p.m_nqh) < LIMSPC );
01613 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
01614 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
01615 if( p.lgEOL() )
01616 p.NoNumb("intensity");
01617
01618
01619
01620 if( !radius.lgRadiusKnown )
01621 {
01622
01623 radius.Radius = pow(10.,radius.rdfalt);
01624 }
01625 if( p.nMatch("LINE") )
01626 {
01627
01628 ASSERT( (p.m_nqh) < LIMSPC );
01629
01630 rfield.totpow[(p.m_nqh)] = log10(rfield.totpow[(p.m_nqh)]);
01631 }
01632 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
01633
01634 ParseRangeOption(p);
01635
01636
01637 if( p.nMatch( "TIME" ) )
01638 rfield.lgTimeVary[(p.m_nqh)] = true;
01639
01640
01641 if( optimize.lgVarOn )
01642 {
01643
01644 strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f LOG range %f %f" );
01645 if( rfield.lgTimeVary[(p.m_nqh)] )
01646 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01647
01648 optimize.nvfpnt[optimize.nparm] = input.nRead;
01649 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
01650 optimize.vincr[optimize.nparm] = 0.5;
01651
01652 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
01653 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
01654 optimize.nvarxt[optimize.nparm] = 3;
01655 ++optimize.nparm;
01656 }
01657 ++(p.m_nqh);
01658 }
01659 void ParseIterations(Parser &p)
01660 {
01661
01662 iterations.itermx = (long int)p.FFmtRead() - 1;
01663 iterations.itermx = MAX2(iterations.itermx,1);
01664
01665
01666
01667 if( iterations.itermx > iterations.iter_malloc - 1 )
01668 {
01669 long int iter_malloc_save = iterations.iter_malloc;
01670
01671
01672 iterations.iter_malloc = iterations.itermx+3;
01673 iterations.IterPrnt = (long int*)REALLOC(iterations.IterPrnt ,
01674 (size_t)iterations.iter_malloc*sizeof(long int) );
01675 geometry.nend = (long int*)REALLOC(geometry.nend ,
01676 (size_t)iterations.iter_malloc*sizeof(long int) );
01677 radius.StopThickness = (double*)REALLOC(radius.StopThickness ,
01678 (size_t)iterations.iter_malloc*sizeof(double) );
01679
01680 for(long int j=iter_malloc_save; j<iterations.iter_malloc; ++j )
01681 {
01682 iterations.IterPrnt[j] = iterations.IterPrnt[iter_malloc_save-1];
01683 geometry.nend[j] = geometry.nend[iter_malloc_save-1];
01684 radius.StopThickness[j] = radius.StopThickness[iter_malloc_save-1];
01685 }
01686 }
01687
01688 if( p.nMatch("CONV") )
01689 {
01690
01691
01692 conv.lgAutoIt = true;
01693
01694 if( p.lgEOL() )
01695 {
01696 iterations.itermx = 10 - 1;
01697 }
01698 realnum a = (realnum)p.FFmtRead();
01699
01700 if( !p.lgEOL() )
01701 {
01702 conv.autocv = a;
01703 }
01704 }
01705 }
01706 void ParseL_nu(Parser &p)
01707 {
01708
01709
01710 bool lgNu2 = false;
01711
01712 ParseF_nu(p,"4 PI",lgNu2);
01713 }
01714 void ParseLaser(Parser &p)
01715 {
01716 DEBUG_ENTRY("ParseLaser()");
01717
01718
01719 strcpy( rfield.chSpType[rfield.nShape], "LASER" );
01720
01721
01722 rfield.slope[rfield.nShape] = p.FFmtRead();
01723
01724
01725 if( rfield.slope[rfield.nShape] <= 0. )
01726 {
01727 rfield.slope[rfield.nShape] =
01728 pow(10.,rfield.slope[rfield.nShape]);
01729 }
01730 if( p.lgEOL() )
01731 {
01732 p.NoNumb("energy");
01733 }
01734
01735
01736 rfield.cutoff[rfield.nShape][0] = p.FFmtRead();
01737 if( p.lgEOL() )
01738 {
01739
01740 rfield.cutoff[rfield.nShape][0] = 0.05;
01741 }
01742
01743
01744 ++rfield.nShape;
01745 if( rfield.nShape >= LIMSPC )
01746 {
01747
01748 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01749 cdEXIT(EXIT_FAILURE);
01750 }
01751 }
01752 void ParseLuminosity(Parser &p)
01753 {
01754 DEBUG_ENTRY("ParseLuminosity()");
01755
01756 if( (p.m_nqh) >= LIMSPC )
01757 {
01758
01759 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01760 cdEXIT(EXIT_FAILURE);
01761 }
01762
01763 strcpy( rfield.chRSpec[(p.m_nqh)], "4 PI" );
01764 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
01765 if( p.lgEOL() )
01766 {
01767 p.NoNumb("luminosity");
01768 }
01769 if( p.nMatch("LINE") )
01770 {
01771
01772 rfield.totpow[(p.m_nqh)] = log10(rfield.totpow[(p.m_nqh)]);
01773 }
01774
01775 strcpy( rfield.chSpNorm[(p.m_nqh)], "LUMI" );
01776
01777
01778 if( p.nMatch("SOLA") )
01779 {
01780
01781 rfield.range[(p.m_nqh)][0] = rfield.emm;
01782 rfield.range[(p.m_nqh)][1] = rfield.egamry;
01783
01784 rfield.totpow[(p.m_nqh)] += 33.5827f;
01785 }
01786 else
01787 {
01788
01789
01790 ParseRangeOption(p);
01791 }
01792
01793
01794 if( p.nMatch( "TIME" ) )
01795 rfield.lgTimeVary[(p.m_nqh)] = true;
01796
01797
01798 if( optimize.lgVarOn )
01799 {
01800 strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f LOG range %f %f" );
01801 if( rfield.lgTimeVary[(p.m_nqh)] )
01802 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01803
01804 optimize.nvfpnt[optimize.nparm] = input.nRead;
01805 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
01806 optimize.vincr[optimize.nparm] = 0.5;
01807
01808 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
01809 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
01810 optimize.nvarxt[optimize.nparm] = 3;
01811 ++optimize.nparm;
01812 }
01813 ++(p.m_nqh);
01814 }
01815 void ParseNeutrons(Parser &p)
01816 {
01817
01818 hextra.lgNeutrnHeatOn = true;
01819
01820
01821
01822 hextra.frcneu = (realnum)p.FFmtRead();
01823 if( p.lgEOL() )
01824 p.NoNumb("neutron luminosity");
01825
01826
01827 if( hextra.frcneu > 0. )
01828 hextra.frcneu = (realnum)log10(hextra.frcneu);
01829
01830
01831 hextra.effneu = (realnum)p.FFmtRead();
01832 if( p.lgEOL() )
01833 {
01834 hextra.effneu = 1.0;
01835 }
01836 else
01837 {
01838 if( hextra.effneu <= 0. )
01839 hextra.effneu = (realnum)pow((realnum)10.f,hextra.effneu);
01840 }
01841 }
01842 void ParseNuF_nu(Parser &p)
01843 {
01844
01845
01846 bool lgNu2 = true;
01847
01848 ParseF_nu(p,"SQCM",lgNu2);
01849 }
01850 void ParseNuL_nu(Parser &p)
01851 {
01852
01853
01854 bool lgNu2 = true;
01855
01856 ParseF_nu(p,"4 PI",lgNu2);
01857 }
01858 void ParsePhi(Parser &p)
01859 {
01860 DEBUG_ENTRY("ParsePhi()");
01861
01862 if( (p.m_nqh) >= LIMSPC )
01863 {
01864
01865 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01866 cdEXIT(EXIT_FAILURE);
01867 }
01868
01869 ASSERT( (p.m_nqh) < LIMSPC );
01870 strcpy( rfield.chRSpec[(p.m_nqh)], "SQCM" );
01871 strcpy( rfield.chSpNorm[(p.m_nqh)], "PHI " );
01872 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
01873 if( p.lgEOL() )
01874 {
01875 p.NoNumb("number of h-ionizing photons");
01876 }
01877
01878
01879 if( !radius.lgRadiusKnown )
01880 {
01881
01882 radius.Radius = pow(10.,radius.rdfalt);
01883 }
01884
01885 if( rfield.totpow[(p.m_nqh)] > 29. )
01886 {
01887 fprintf( ioQQQ, " Is the flux for this continuum correct?\n" );
01888 fprintf( ioQQQ, " It appears too bright to me.\n" );
01889 }
01890
01891 ParseRangeOption(p);
01892
01893
01894 if( p.nMatch( "TIME" ) )
01895 rfield.lgTimeVary[(p.m_nqh)] = true;
01896
01897
01898 if( optimize.lgVarOn )
01899 {
01900 strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f LOG range %f %f" );
01901 if( rfield.lgTimeVary[(p.m_nqh)] )
01902 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01903
01904 optimize.nvfpnt[optimize.nparm] = input.nRead;
01905 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
01906 optimize.vincr[optimize.nparm] = 0.5;
01907
01908 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
01909 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
01910 optimize.nvarxt[optimize.nparm] = 3;
01911 ++optimize.nparm;
01912 }
01913 ++(p.m_nqh);
01914 }
01915 void ParseQH(Parser &p)
01916 {
01917 DEBUG_ENTRY("ParseQH()");
01918
01919 if( (p.m_nqh) >= LIMSPC )
01920 {
01921
01922 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
01923 cdEXIT(EXIT_FAILURE);
01924 }
01925
01926
01927 ASSERT( (p.m_nqh) < LIMSPC );
01928 strcpy( rfield.chRSpec[(p.m_nqh)], "4 PI" );
01929 strcpy( rfield.chSpNorm[(p.m_nqh)], "Q(H)" );
01930 rfield.totpow[(p.m_nqh)] = p.FFmtRead();
01931 if( rfield.totpow[(p.m_nqh)] > 100. && called.lgTalk )
01932 {
01933 fprintf( ioQQQ, " Is this reasonable?\n" );
01934 }
01935 if( p.lgEOL() )
01936 {
01937 p.NoNumb("number of ionizing photons");
01938 }
01939
01940 ParseRangeOption(p);
01941
01942
01943 if( p.nMatch( "TIME" ) )
01944 rfield.lgTimeVary[(p.m_nqh)] = true;
01945
01946
01947 if( optimize.lgVarOn )
01948 {
01949 strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f LOG range %f %f" );
01950 if( rfield.lgTimeVary[(p.m_nqh)] )
01951 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01952
01953 optimize.nvfpnt[optimize.nparm] = input.nRead;
01954 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[(p.m_nqh)];
01955 optimize.vincr[optimize.nparm] = 0.5;
01956
01957 optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][0]);
01958 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[(p.m_nqh)][1]);
01959 optimize.nvarxt[optimize.nparm] = 3;
01960 ++optimize.nparm;
01961 }
01962
01963 ++(p.m_nqh);
01964 }
01965 void ParseRoberto(Parser &)
01966 {
01967
01968 radius.dRadSign = -1.;
01969 }
01970 void ParseSpecial(Parser &)
01971 {
01972 DEBUG_ENTRY("ParseSpecial()");
01973
01974 cdEXIT(EXIT_FAILURE);
01975 }
01976 void ParseTauMin(Parser &p)
01977 {
01978
01979 opac.taumin = (realnum)pow(10.,p.FFmtRead());
01980 if( p.lgEOL() )
01981 p.NoNumb("minimum optical depth");
01982 }
01983 void ParseTitle(Parser &p)
01984 {
01985
01986
01987
01988 if (p.GetQuote(input.chTitle,false) != 0)
01989 strcpy( input.chTitle , p.getRawTail().c_str()+1 );
01990 }
01991 void ParseTolerance(Parser &)
01992 {
01993 DEBUG_ENTRY("ParseTolerance()");
01994 fprintf(ioQQQ,
01995 "Sorry, this command has been replaced with the SET TEMPERATURE TOLERANCE command.\n");
01996 cdEXIT(EXIT_FAILURE);
01997 }
01998 void ParseVLaw(Parser &p)
01999 {
02000
02001 DoppVel.TurbVelLaw = (realnum)p.FFmtRead();
02002
02003 DoppVel.lgTurbLawOn = true;
02005
02006
02007 ASSERT( DoppVel.TurbVelLaw <= 0.f );
02008 }
02009 void ParseTurbulence(Parser &p)
02010 {
02011 DEBUG_ENTRY("ParseTurbulence()");
02012 string ExtraPars;
02013
02014 if( p.nMatch("EQUIPART") )
02015 {
02016
02017 DoppVel.lgTurbEquiMag = true;
02018
02019
02020
02021
02022 DoppVel.Heiles_Troland_F = (realnum)p.FFmtRead();
02023 if( p.lgEOL() )
02024 {
02025
02026 DoppVel.Heiles_Troland_F = 3.f;
02027 }
02028 }
02029 else
02030 {
02031
02032 DoppVel.lgTurbEquiMag = false;
02033 DoppVel.TurbVel = (realnum)p.FFmtRead();
02034 if( p.lgEOL() )
02035 p.NoNumb("microturbulent velocity");
02036
02037 if( p.nMatch(" LOG") )
02038 {
02039 if( DoppVel.TurbVel > 32. )
02040 {
02041 fprintf( ioQQQ, "PROBLEM the log of the turbulence is "
02042 "%.2e - I cannot handle a number this big.\n",
02043 DoppVel.TurbVel );
02044 fprintf( ioQQQ, " The line image was\n" );
02045 p.PrintLine(ioQQQ);
02046 fprintf( ioQQQ, " Sorry.\n" );
02047 cdEXIT(EXIT_FAILURE);
02048 }
02049 DoppVel.TurbVel = (realnum)pow((realnum)10.f,DoppVel.TurbVel);
02050 }
02051
02052 DoppVel.TurbVel *= 1e5;
02053
02054 if( DoppVel.TurbVel <= 0. )
02055 {
02056 fprintf( ioQQQ, " PROBLEM: the turbulent velocity needs to be > 0, but this was entered: %e\n",
02057 DoppVel.TurbVel );
02058 fprintf( ioQQQ, " Bailing out. Sorry.\n" );
02059 cdEXIT(EXIT_FAILURE);
02060 }
02061 else if( DoppVel.TurbVel >= SPEEDLIGHT )
02062 {
02063 fprintf( ioQQQ, " PROBLEM: A turbulent velocity greater than speed of light is not allowed, this was entered: %e\n",
02064 DoppVel.TurbVel );
02065 fprintf( ioQQQ, " Bailing out. Sorry.\n" );
02066 cdEXIT(EXIT_FAILURE);
02067 }
02068
02069
02070
02071
02072 DoppVel.Heiles_Troland_F = (realnum)p.FFmtRead();
02073 if( p.lgEOL() )
02074 {
02075
02076 DoppVel.Heiles_Troland_F = 3.f;
02077 }
02078
02079
02080
02081 if( p.nMatch("DISS") )
02082 {
02083 DoppVel.DispScale = (realnum)pow(10., p.FFmtRead() );
02084 if( p.lgEOL() )
02085 p.NoNumb("turbulence dissipation scale");
02086 ExtraPars += " DISSIPATE %f";
02087 }
02088 }
02089
02090
02091
02092
02093 if( p.nMatch(" NO ") && p.nMatch("PRES") )
02094 {
02095 DoppVel.lgTurb_pressure = false;
02096 ExtraPars += " NO PRESSURE";
02097 }
02098 else
02099 {
02100
02101 DoppVel.lgTurb_pressure = true;
02102 }
02103
02104
02105 if( optimize.lgVarOn && !p.nMatch("EQUIPART") )
02106 {
02107
02108 optimize.nvarxt[optimize.nparm] = 2;
02109
02110 strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f LOG %f" );
02111 strcat( optimize.chVarFmt[optimize.nparm], ExtraPars.c_str() );
02112
02113
02114 optimize.nvfpnt[optimize.nparm] = input.nRead;
02115
02116 optimize.vparm[0][optimize.nparm] = log10(DoppVel.TurbVel/1e5f);
02117 optimize.vparm[1][optimize.nparm] = DoppVel.Heiles_Troland_F;
02118 if( p.nMatch("DISS") )
02119 {
02120 optimize.nvarxt[optimize.nparm] = 3;
02121 optimize.vparm[2][optimize.nparm] = log10(DoppVel.DispScale);
02122 }
02123 optimize.vincr[optimize.nparm] = 0.1f;
02124 ++optimize.nparm;
02125 }
02126
02127
02128 DoppVel.TurbVelZero = DoppVel.TurbVel;
02129 }