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