cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ParseCommands main command line parser, decode command, then call other routines to read */ 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 "assertresults.h" 00033 #include "phycon.h" 00034 #include "fudgec.h" 00035 #include "version.h" 00036 #include "conv.h" 00037 #include "parse.h" 00038 00039 void ParseCommands(void) 00040 { 00041 char chCard[INPUT_LINE_LENGTH], 00042 chKey2[3], 00043 chKey3[4], 00044 chKey4[5], 00045 chKey5[6]; 00046 00047 bool lgDSet, 00048 lgEOF, 00049 lgEOL, 00050 lgNu2; 00051 bool lgStop , 00052 lgStop_not_enough_info; 00053 00054 long int i, 00055 j, 00056 nqh; 00057 long nInitFile=0;/* used to count number of init files read in */ 00058 00059 realnum a, 00060 /* \todo 1 ar1 saves initial radius but does not appear 00061 * to be used for anything. all uses are on left side of = 00062 * this could be removed across the code. But first check 00063 * that it is indeed not ever used */ 00064 ar1, 00065 teset, 00066 z; 00067 double thickness; 00068 00069 DEBUG_ENTRY( "ParseCommands()" ); 00070 00071 /* following says abundances are solar */ 00072 abund.lgAbnSolar = true; 00073 00074 /* there are no plots desired yet */ 00075 plotCom.lgPlotON = false; 00076 plotCom.nplot = 0; 00077 00078 /* this flag remembers whether grains have ever been turned on. It is passed 00079 * to routine ParseAbundances - there grains will be turned on with commands 00080 * such as abundances ism, unless grains were previously set 00081 * with a grains command. this way abundances will not clobber explicitly set 00082 * grains. */ 00083 lgDSet = false; 00084 00085 radius.Radius = 0.; 00086 radius.rinner = 0.; 00087 nqh = 0; 00088 rfield.nspec = 0; 00089 00090 /* will be set true if no buffering command entered */ 00091 input.lgSetNoBuffering = false; 00092 00093 /* initialize some code variables in case assert command used in input stream */ 00094 InitAssertResults(); 00095 00096 for( i=0; i < LIMSPC; i++ ) 00097 { 00098 strcpy( rfield.chRSpec[i], "UNKN" ); 00099 } 00100 optimize.nparm = 0; 00101 00102 /* this is an option to turn on trace printout on the nth 00103 * call from the optimizer */ 00104 if( optimize.lgTrOpt ) 00105 { 00106 /* nTrOpt was set with "optimize trace" command, 00107 * is iteration to turn on trace */ 00108 if( optimize.nTrOpt == optimize.nOptimiz ) 00109 { 00110 trace.lgTrace = true; 00111 called.lgTalk = true; 00112 trace.lgTrOvrd = true; 00113 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" ); 00114 } 00115 } 00116 00117 /* say this is a beta version if we are talking and it is the truth */ 00118 if( t_version::Inst().nBetaVer > 0 && called.lgTalk ) 00119 { 00120 fprintf( ioQQQ, 00121 "\n This is a beta release of Cloudy, and is intended for testing only.\n\n" ); 00122 } 00123 00124 if( called.lgTalk ) 00125 { 00126 /* this code prints pretty lines at top of output box */ 00127 int indent = (int)((122 - strlen(t_version::Inst().chVersion))/2); 00128 fprintf( ioQQQ, "%*cCloudy %s\n", indent, ' ', t_version::Inst().chVersion ); 00129 00130 fprintf( ioQQQ, "%57cwww.nublado.org\n\n", ' ' ); 00131 00132 /* now print box and date of version, before printing commands */ 00133 fprintf( ioQQQ, "%23c", ' ' ); 00134 fprintf( ioQQQ, "**************************************"); 00135 fprintf( ioQQQ, "%7.7s", t_version::Inst().chDate); 00136 fprintf( ioQQQ, "**************************************\n"); 00137 00138 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' ); 00139 } 00140 00141 /* read in commands and print them */ 00142 00143 /* following signals which read is in progress, 00144 * 1 is forward read, of input commands 00145 * -1 is reverse read from bottom of line array, of cloudy.ini file */ 00146 input.iReadWay = 1; 00147 00148 /* initialize array reader, this sub does nothing but set 00149 * initial value of a variable, depending on value of iReadWay */ 00150 input_init(); 00151 00152 /* read until eof or blank line, then return control back to main program */ 00153 00154 /* 00155 * input_readarray puts the next stored line image into chCard 00156 * it will be <=INPUT_LINE_LENGTH char long, with eol at end 00157 * it will be in mixed upper and lower case 00158 * lgEOF is set true if we hit eof and no more lines present 00159 * code moving over to info contained in structures, 00160 * 00161 * input.chOrgCard == original image of command line 00162 * input.chCARDCAPS == original image converted to caps 00163 */ 00164 input_readarray(chCard,&lgEOF); 00165 00166 /* line starting with blank is eof, this checks we are not at eof */ 00167 while( !lgEOF && chCard[0] != ' ' ) 00168 { 00169 /* echo the line before we cap it, but only if it does 00170 * not contain the keyword HIDE */ 00171 /* >>chng 04 jan 21, add HIDE option, mostly for print quiet command */ 00172 if( called.lgTalk && !nMatch("HIDE",input.chCARDCAPS) ) 00173 fprintf( ioQQQ, "%23c* %-80s*\n", ' ', chCard ); 00174 00175 /* change chCard to all caps */ 00176 caps( chCard ); 00177 00178 /* now set up several partial keys */ 00179 /* first two characters into chKey2 */ 00180 strncpy( chKey2 , chCard , 2 ); 00181 chKey2[2] = '\0'; 00182 00183 /* first three characters into chKey3 */ 00184 strncpy( chKey3 , chCard , 3 ); 00185 chKey3[3] = '\0'; 00186 00187 /* first four characters into chKey4 */ 00188 strncpy( chKey4 , chCard , 4 ); 00189 chKey4[4] = '\0'; 00190 00191 /* first four characters into chKey4 */ 00192 strncpy( chKey5 , chCard , 5 ); 00193 chKey5[5] = '\0'; 00194 00195 /* check whether VARY is on line */ 00196 if( nMatch("VARY",chCard) ) 00197 { 00198 optimize.lgVarOn = true; 00199 if( optimize.nparm + 1 > LIMPAR ) 00200 { 00201 fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n", 00202 LIMPAR ); 00203 cdEXIT(EXIT_FAILURE); 00204 } 00205 } 00206 00207 else 00208 { 00209 optimize.lgVarOn = false; 00210 } 00211 00212 /* look for "C " comment - this form of test should not be necessary since 00213 * input line is padded with extra space - but works and would continue to 00214 * work if padding were removed */ 00215 if( chCard[0]=='C' && (chCard[1]==' ' || chCard[1]== 0) ) 00216 { 00217 /* this is null test since lines starting with "C " are comments, 00218 * the char after the c can be a space or a newline */ 00219 i = 1; 00220 } 00221 00222 /* start to look for keywords */ 00223 else if( strcmp(chKey4,"ABSO") == 0 ) 00224 { 00225 /* enter luminosity in absolute magnitudes, in reads2 */ 00226 ParseAbsMag(chCard,&nqh); 00227 } 00228 00229 else if( strcmp(chKey3,"AGE") == 0 ) 00230 { 00231 /* enter age of cloud so we can check for time-steady reads2 */ 00232 ParseAge(chCard); 00233 } 00234 00235 else if( strcmp(chKey4,"AGN ") == 0 ) 00236 { 00237 /* enter generic style AGN continuum, in reads2 */ 00238 ParseAgn(chCard); 00239 } 00240 00241 else if( strcmp(chKey4,"ABUN") == 0 ) 00242 { 00243 /* chemical abundances, readsun, lgDSet is set true if grains command 00244 * ever entered, so that abundances command does not override grains command */ 00245 ParseAbundances(chCard,lgDSet); 00246 /* abundances no longer solar */ 00247 abund.lgAbnSolar = false; 00248 } 00249 00250 else if( strcmp(chKey4,"APER") == 0 ) 00251 { 00252 /* aperture command to simulate pencil beam or long slit */ 00253 00254 /* if the "BEAM" or "SLIT" option is specified then only part 00255 * of the geometry is observed, and intensities 00256 * should not be weighted by r^2. There are two limiting cases, SLIT in which 00257 * the slit is longer than the diameter of the nebula and the contribution to the 00258 * detected luminosity goes as r^1, and BEAM when the contribution is r^0, 00259 * the same as plane parallel 00260 * 00261 * default value of geometry.iEmissPower is 2 (set in zero.c) for full geometry 00262 */ 00263 if( nMatch("SLIT",chCard) ) 00264 { 00265 /* long slit is case where slit is longer than diameter, so emissivity contributes 00266 * r^1 to the observed luminosity */ 00267 geometry.iEmissPower = 1; 00268 } 00269 else if( nMatch("BEAM",chCard) ) 00270 { 00271 /* long slit is case where slit is longer than diameter, so emissivity contributes 00272 * r^1 to the observed luminosity */ 00273 /* this is the default or SHORT case, where r^0 is dependence */ 00274 geometry.iEmissPower = 0; 00275 } 00276 else 00277 { 00278 fprintf( ioQQQ, " One of the keywords SLIT or BEAM must appear.\n" ); 00279 fprintf( ioQQQ, " Sorry.\n" ); 00280 cdEXIT(EXIT_FAILURE); 00281 } 00282 } 00283 00284 else if( strcmp(chKey4,"ASSE") == 0 ) 00285 { 00286 /* assert that code predicts certain results, in assertresults.h */ 00287 ParseAssertResults(); 00288 } 00289 00290 else if( strcmp(chKey4,"ATOM") == 0 ) 00291 { 00292 /* accept both forms of feii */ 00293 if( nMatch("FEII",chCard) || nMatch("FE II",chCard) ) 00294 { 00295 /* parse the atom FeII command - routine in atom_feii.c */ 00296 ParseAtomFeII(chCard); 00297 } 00298 00299 else if( nMatch("H-LI",chCard) ) 00300 { 00301 /* parse the atom h-like command */ 00302 ParseAtomISO(ipH_LIKE, chCard); 00303 } 00304 00305 else if( nMatch("HE-L",chCard) ) 00306 { 00307 /* parse the atom he-like command */ 00308 ParseAtomISO(ipHE_LIKE, chCard); 00309 } 00310 00311 else if( nMatch("ROTO",chCard) ) 00312 { 00313 /* ATOM ROTOR same as ATOM CO */ 00314 fprintf(ioQQQ,"PROBLEM - the atom rotor command is now the ATOM CO command. " 00315 "Please use that instead. I will accept this for now but may not in future versions.\n"); 00316 ParseAtomCO(chCard); 00317 } 00318 00319 else if( nMatch(" CO ",chCard) ) 00320 { 00321 /* parameters for the rigid rotor CO molecules (yes, its not an atom) 00322 * command is ATOM CO */ 00323 ParseAtomCO(chCard); 00324 } 00325 00326 else if( nMatch(" H2 ",chCard) ) 00327 { 00328 /* parameters for the rigid rotor CO molecules (yes, its not an atom) 00329 * command is ATOM ROTOR, routine is stored in parse_atomh2 */ 00330 ParseAtomH2(chCard); 00331 } 00332 00333 else 00334 { 00335 fprintf( ioQQQ, " I could not recognize a keyword on this atom command.\n"); 00336 fprintf( ioQQQ, " The available keys are FeII, H-Like, He-like, rotor and H2.\n"); 00337 fprintf( ioQQQ, " Sorry.\n" ); 00338 cdEXIT(EXIT_FAILURE); 00339 } 00340 } 00341 00342 else if( strcmp(chKey4,"BACK") == 0 ) 00343 { 00344 /* cosmic background, in parse_backgrd */ 00345 ParseBackgrd(&nqh,chCard,&ar1); 00346 } 00347 00348 else if( strcmp(chKey4,"BLAC") == 0 ) 00349 { 00350 /* black body, in reads2 */ 00351 ParseBlackbody(chCard,&nqh,&ar1); 00352 00353 /* vary option */ 00354 if( optimize.lgVarOn ) 00355 { 00356 /* no luminosity options on vary */ 00357 optimize.nvarxt[optimize.nparm] = 1; 00358 strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody=%f" ); 00359 /* pointer to where to write */ 00360 optimize.nvfpnt[optimize.nparm] = input.nRead; 00361 /* log of temp stored here */ 00362 /* >>chng 02 feb 11, log was not here, don't know what happened to it */ 00363 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nspec-1]); 00364 /* the increment in the first steps away from the original value */ 00365 optimize.vincr[optimize.nparm] = 0.5; 00366 optimize.nparm += 1; 00367 } 00368 } 00369 00370 else if( strcmp(chKey4,"BREM") == 0 ) 00371 { 00372 /* bremsstrahlung continuum from central object */ 00373 strcpy( rfield.chSpType[rfield.nspec], "BREMS" ); 00374 i = 5; 00375 rfield.slope[rfield.nspec] = 00376 (realnum)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL); 00377 if( lgEOL ) 00378 { 00379 NoNumb(chCard); 00380 } 00381 00382 /* temperature is interpreted as log if <=10, linear otherwise*/ 00383 if( rfield.slope[rfield.nspec] <= 10. ) 00384 { 00385 rfield.slope[rfield.nspec] = 00386 pow(10.,rfield.slope[rfield.nspec]); 00387 } 00388 rfield.cutoff[rfield.nspec][0] = 0.; 00389 00390 /* option for vary keyword */ 00391 if( optimize.lgVarOn ) 00392 { 00393 /* only one parameter */ 00394 optimize.nvarxt[optimize.nparm] = 1; 00395 strcpy( optimize.chVarFmt[optimize.nparm], "BREMS, T=%f" ); 00396 /* pointer to where to write */ 00397 optimize.nvfpnt[optimize.nparm] = input.nRead; 00398 /* log of temp will be pointer */ 00399 optimize.vparm[0][optimize.nparm] = (realnum)rfield.slope[rfield.nspec]; 00400 optimize.vincr[optimize.nparm] = 0.5; 00401 ++optimize.nparm; 00402 } 00403 ++rfield.nspec; 00404 if( rfield.nspec >= LIMSPC ) 00405 { 00406 /* too many continua were entered */ 00407 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00408 cdEXIT(EXIT_FAILURE); 00409 } 00410 } 00411 00412 else if( strcmp(chKey4,"CASE") == 0 ) 00413 { 00414 /* do hydrogen case b */ 00415 ParseCaseB( chCard ); 00416 } 00417 00418 else if( strcmp(chKey4,"CEXT") == 0 ) 00419 { 00420 /* add "extra" cooling, power law temp dependence */ 00421 thermal.lgCExtraOn = true; 00422 i = 5; 00423 thermal.CoolExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00424 if( lgEOL ) 00425 { 00426 NoNumb(chCard); 00427 } 00428 thermal.cextpw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00429 } 00430 00431 /* replace with CMB command, both will continue to be parsed */ 00432 else if( (strcmp(chKey3,"CMB") == 0) || (strcmp(chKey4,"FIRE") == 0) ) 00433 { 00434 /* cosmic thermal background radiation, argument is redshift */ 00435 i = 5; 00436 /* if no number on line then (realnum)FFmtRead returns z=0; i.e., now */ 00437 z = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00438 /* in readsun */ 00439 ParseCMB(z,&nqh,&ar1); 00440 } 00441 00442 else if( strcmp(chKey4,"COMP") == 0 ) 00443 { 00444 /* compile ascii version of stellar atmosphere continua in volk */ 00445 ParseCompile(chCard); 00446 } 00447 00448 else if( strcmp(chKey4,"CONS") == 0 ) 00449 { 00450 /* constant temperature, pressure, density, or gas pressure 00451 * in readsun */ 00452 ParseConstant(chCard); 00453 } 00454 00455 else if( strcmp(chKey4,"CORO") == 0 ) 00456 { 00457 /* coronal equilibrium; set constant temperature to number on line 00458 * in readsun */ 00459 ParseCoronal( chCard,&nqh,&ar1); 00460 } 00461 00462 else if( strcmp(chKey4,"COSM") == 0 ) 00463 { 00464 /* cosmic ray density, log of rate relative to background, or density of rel. electrons, 00465 * quantity is log unless keyword linear appears */ 00466 ParseCosmicRays( chCard ); 00467 } 00468 00469 else if( strcmp(chKey4,"COVE") == 0 ) 00470 { 00471 /* covering factor for gas */ 00472 i = 5; 00473 /* The geometric covering factor accounts for how much of 4\pi is 00474 * covered by gas, and so linearly multiplies the predicted intensities */ 00475 geometry.covgeo = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00476 00477 if( lgEOL ) 00478 { 00479 NoNumb(chCard); 00480 } 00481 00482 /* if negative then log, convert to linear */ 00483 if( geometry.covgeo <= 0. ) 00484 { 00485 geometry.covgeo = (realnum)pow((realnum)10.f,geometry.covgeo); 00486 } 00487 00488 if( geometry.covgeo > 1. ) 00489 { 00490 fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense. Sorry.\n" ); 00491 cdEXIT(EXIT_FAILURE); 00492 } 00493 00494 /* radiative transfer covering factor will be equal to the geometric one */ 00495 geometry.covrt = geometry.covgeo; 00496 } 00497 00498 else if( strcmp(chKey4,"CRAS") == 0 ) 00499 { 00500 /* any of several tests to cause the code to crash */ 00501 ParseCrashDo(chCard); 00502 } 00503 00504 else if( strcmp(chKey4,"CYLI") == 0 ) 00505 { 00506 /* height of cylinder in cm */ 00507 i = 5; 00508 radius.lgCylnOn = true; 00509 radius.CylindHigh = pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00510 if( lgEOL ) 00511 { 00512 NoNumb(chCard); 00513 } 00514 } 00515 00516 else if( strcmp(chKey4,"DIEL") == 0 ) 00517 { 00518 /* >>chng 05 dec 21, change dielectronic command to set dielectronic recombination command */ 00519 fprintf( ioQQQ, " The DIELectronic command has been replaced with the SET DIELectronic recombination command.\n" ); 00520 fprintf( ioQQQ, " Please have a look at Hazy.\n Sorry.\n\n" ); 00521 cdEXIT(EXIT_FAILURE); 00522 } 00523 00524 else if( strcmp(chKey4,"DIFF") == 0 ) 00525 { 00526 /* set method of transferring diffuse fields, 00527 * default is outward only, cdDffTrns label "OU2", set in zero.c */ 00528 if( nMatch(" OTS",chCard) ) 00529 { 00530 if( nMatch("SIMP",chCard) ) 00531 { 00532 /* this is simple ots, which never allows photons to escape */ 00533 strcpy( rfield.chDffTrns, "OSS" ); 00534 } 00535 else 00536 { 00537 /* this is regular ots, which allows photons to escape */ 00538 strcpy( rfield.chDffTrns, "OTS" ); 00539 } 00540 rfield.lgOutOnly = false; 00541 } 00542 else if( nMatch(" OUT",chCard) ) 00543 { 00544 rfield.lgOutOnly = true; 00545 i = 4; 00546 j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00547 if( lgEOL ) 00548 { 00549 /* this is the default set in zero */ 00550 strcpy( rfield.chDffTrns, "OU2" ); 00551 } 00552 else 00553 { 00554 if( j > 0 && j < 10 ) 00555 { 00556 sprintf( rfield.chDffTrns, "OU%1ld", j ); 00557 } 00558 else 00559 { 00560 fprintf( ioQQQ, " must be between 1 and 9 \n" ); 00561 cdEXIT(EXIT_FAILURE); 00562 } 00563 } 00564 } 00565 00566 else 00567 { 00568 fprintf( ioQQQ, " There should have been OUTward or OTS on this line. Sorry.\n" ); 00569 cdEXIT(EXIT_FAILURE); 00570 } 00571 } 00572 00573 else if( strcmp(chKey4,"DIST") == 0 ) 00574 { 00575 /* distance to the object */ 00576 i = 5; 00577 radius.distance = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00578 if( lgEOL ) 00579 { 00580 NoNumb(chCard); 00581 } 00582 /* default is for quantity to be log of distance, linear keyword 00583 * overrides this - is LINEAR is not on line then exp */ 00584 if( !nMatch("LINE",chCard ) ) 00585 { 00586 radius.distance = pow(10., radius.distance ); 00587 } 00588 /* default is radius in cm - if parsecs appears then must 00589 * convert to cm */ 00590 if( nMatch("PARS",chCard ) ) 00591 { 00592 radius.distance *= PARSEC; 00593 } 00594 } 00595 00596 else if( strcmp(chKey4,"DLAW") == 0 ) 00597 { 00598 /* either use dense_fabden routine, or read in table of points */ 00599 ParseDLaw(chCard); 00600 } 00601 00602 else if( strcmp(chKey4,"DOUB") == 0 ) 00603 { 00604 /* double optical depth scale after each iteration */ 00605 rt.DoubleTau = 2.; 00606 } 00607 00608 else if( strcmp(chKey4,"DRIV") == 0 ) 00609 { 00610 /* call one of several drivers, readsun */ 00611 ParseDrive(chCard); 00612 } 00613 00614 else if( strcmp(chKey4,"EDEN") == 0 ) 00615 { 00616 /* option to add "extra" electrons, to test Compton limit 00617 * for very low T(star) - option is log of eden */ 00618 i = 5; 00619 dense.EdenExtra = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00620 if( lgEOL ) 00621 { 00622 NoNumb(chCard); 00623 } 00624 /* warn that this model is meaningless */ 00625 phycon.lgPhysOK = false; 00626 } 00627 00628 else if( strcmp(chKey4,"ELEM") == 0 ) 00629 { 00630 /* element command; 00631 * scale or abundance options, to change abundance of specific element 00632 * read option to change order of elements 00633 * in reads2.f */ 00634 ParseElement(chCard); 00635 } 00636 00637 else if( strcmp(chKey4,"ENER") == 0 ) 00638 { 00639 /* energy density (degrees K) of this continuum source */ 00640 i = 5; 00641 if( nqh >= LIMSPC ) 00642 { 00643 /* too many continua were entered */ 00644 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00645 cdEXIT(EXIT_FAILURE); 00646 } 00647 /* silly, but calms down the lint */ 00648 ASSERT( nqh < LIMSPC ); 00649 strcpy( rfield.chRSpec[nqh], "SQCM" ); 00650 teset = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00651 if( lgEOL ) 00652 { 00653 NoNumb(chCard); 00654 } 00655 /* set radius to very large value if not already set */ 00656 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00657 if( !radius.lgRadiusKnown ) 00658 { 00659 /* set radius to large value */ 00660 ar1 = (realnum)radius.rdfalt; 00661 radius.Radius = pow(10.,radius.rdfalt); 00662 } 00663 00664 /* convert temp to log, recognize linear option */ 00665 if( nMatch("LINE",chCard) || teset > 10. ) 00666 { 00667 /* option for linear temperature, must store log */ 00668 teset = (realnum)log10(teset); 00669 } 00670 00671 if( teset > 5. ) 00672 { 00673 fprintf( ioQQQ, " This intensity may be too large. The code may crash due to overflow. Was log intended?\n" ); 00674 } 00675 00676 /* teset is not log of temp, now get log of total luminosity */ 00677 strcpy( rfield.chSpNorm[nqh], "LUMI" ); 00678 00679 /* full range of continuum will be used */ 00680 rfield.range[nqh][0] = rfield.emm; 00681 rfield.range[nqh][1] = rfield.egamry; 00682 rfield.totpow[nqh] = (4.*teset - 4.2464476 + 0.60206); 00683 00684 /* >>chng 06 mar 22, add time option to vary only some continua with time */ 00685 if( nMatch( "TIME" , chCard ) ) 00686 rfield.lgTimeVary[nqh] = true; 00687 00688 /* vary option */ 00689 if( optimize.lgVarOn ) 00690 { 00691 strcpy( optimize.chVarFmt[optimize.nparm], "ENERGY DENSITY %f log " ); 00692 /* pointer to where to write */ 00693 optimize.nvfpnt[optimize.nparm] = input.nRead; 00694 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.totpow[nqh]); 00695 optimize.vincr[optimize.nparm] = 0.1f; 00696 optimize.nvarxt[optimize.nparm] = 1; 00697 optimize.nparm += 1; 00698 } 00699 00700 /* finally increment number of continua */ 00701 ++nqh; 00702 } 00703 00704 else if( strcmp(chKey4,"EXTI") == 0 ) 00705 { 00706 /* extinguish ionizing continuum by absorbing column AFTER 00707 * setting luminosity or Q(H). First number is the column 00708 * density (log), second number is leakage (def=0%) 00709 * last number is lowest energy (ryd), last two may be omitted 00710 * from right to left 00711 * 00712 * extinction is actually done in extin, which is called by ContSetIntensity */ 00713 ParseExtinguish( chCard ); 00714 } 00715 00716 else if( strcmp(chKey4,"FAIL") == 0 ) 00717 { 00718 /* reset number of temp failures allowed, default=20 */ 00719 i = 5; 00720 00721 /* save previous value */ 00722 j = conv.LimFail; 00723 conv.LimFail = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00724 if( lgEOL ) 00725 { 00726 NoNumb(chCard); 00727 } 00728 00729 /* >>chng 01 mar 14, switch logic on maps */ 00730 /* ' map' flag, make map when failure, default is no map, 00731 * second check is so that no map does not trigger a map */ 00732 if( nMatch(" MAP",chCard) && !nMatch(" NO ",chCard) ) 00733 { 00734 conv.lgMap = true; 00735 } 00736 00737 /* complain if failures was increased above default */ 00738 if( conv.LimFail > j ) 00739 { 00740 fprintf( ioQQQ, " This command should not be necessary.\n" ); 00741 fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this simulation.\n" ); 00742 } 00743 } 00744 00745 else if( strcmp(chKey4,"FILL") == 0 ) 00746 { 00747 /* filling factor, power law exponent (default=1., 0.) */ 00748 i = 5; 00749 a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00750 if( lgEOL ) 00751 { 00752 NoNumb(chCard); 00753 } 00754 00755 if( a <= 0. ) 00756 { 00757 /* number less than or equal to 0, must have been entered as log */ 00758 geometry.FillFac = (realnum)pow((realnum)10.f,a); 00759 } 00760 else 00761 { 00762 /* number greater than zero, must have been the real thing */ 00763 geometry.FillFac = a; 00764 if( geometry.FillFac > 1.0 ) 00765 { 00766 if( called.lgTalk ) 00767 { 00768 fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" ); 00769 } 00770 geometry.FillFac = 1.; 00771 } 00772 } 00773 geometry.fiscal = geometry.FillFac; 00774 00775 /* option to have power law dependence, filpow set to 0 in zerologic */ 00776 geometry.filpow = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00777 00778 /* vary option */ 00779 if( optimize.lgVarOn ) 00780 { 00781 strcpy( optimize.chVarFmt[optimize.nparm], "FILLING FACTOR= %f power=%f" ); 00782 00783 /* pointer to where to write */ 00784 optimize.nvfpnt[optimize.nparm] = input.nRead; 00785 optimize.vparm[0][optimize.nparm] = (realnum)log10(geometry.FillFac); 00786 optimize.vincr[optimize.nparm] = 0.5; 00787 00788 /* power law dependence here, but cannot be varied */ 00789 optimize.vparm[1][optimize.nparm] = geometry.filpow; 00790 optimize.nvarxt[optimize.nparm] = 2; 00791 00792 /* do not allow filling factor to go positive */ 00793 optimize.varang[optimize.nparm][0] = -1e38f; 00794 optimize.varang[optimize.nparm][1] = 0.f; 00795 optimize.nparm += 1; 00796 } 00797 } 00798 00799 else if( strcmp(chKey4,"FLUC") == 0 ) 00800 { 00801 /* rapid density fluctuations, in readsun */ 00802 ParseFluc(chCard); 00803 } 00804 00805 else if( strcmp(chKey4,"F(NU") == 0 ) 00806 { 00807 /* this is the specific flux at nu 00808 * following says F(nu) not nuF(nu) */ 00809 lgNu2 = false; 00810 /* in reads2 */ 00811 ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2); 00812 } 00813 00814 else if( strcmp(chKey4,"FORC") == 0 ) 00815 { 00816 /* force temperature of first zone, but don't keep constant 00817 * allow to then go to nearest equilibrium 00818 * log if < 10 */ 00819 i = 5; 00820 thermal.ConstTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00821 if( lgEOL ) 00822 { 00823 NoNumb(chCard); 00824 } 00825 00826 if( nMatch(" LOG",chCard) || (thermal.ConstTemp <= 10. && 00827 !nMatch("LINE",chCard)) ) 00828 { 00829 thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp); 00830 } 00831 00832 /* low energy bounds of continuum array do not permit too-low a temp */ 00833 if( thermal.ConstTemp < 3. ) 00834 { 00835 fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" ); 00836 thermal.ConstTemp = 3.; 00837 } 00838 } 00839 00840 else if( strcmp(chKey4,"FUDG") == 0 ) 00841 { 00842 /* enter a fudge factor */ 00843 i = 5; 00844 for( j=0; j < NFUDGC; j++ ) 00845 { 00846 fudgec.fudgea[j] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00847 /* fudgec.nfudge is number of factors on the line */ 00848 if( !lgEOL ) 00849 fudgec.nfudge = j+1; 00850 } 00851 00852 /* vary option */ 00853 if( optimize.lgVarOn ) 00854 { 00855 /* no luminosity options on vary */ 00856 optimize.nvarxt[optimize.nparm] = 1; 00857 strcpy( optimize.chVarFmt[optimize.nparm], "FUDGE=%f" ); 00858 /* enter remaining parameters as constants */ 00859 char chHold[1000]; 00860 for( j=1; j<fudgec.nfudge; ++j ) 00861 { 00862 sprintf( chHold , " %f" , fudgec.fudgea[j] ); 00863 strcat( optimize.chVarFmt[optimize.nparm] , chHold ); 00864 } 00865 00866 /* pointer to where to write */ 00867 optimize.nvfpnt[optimize.nparm] = input.nRead; 00868 /* fudge factor stored here */ 00869 optimize.vparm[0][optimize.nparm] = fudgec.fudgea[0]; 00870 /* the increment in the first steps away from the original value */ 00871 optimize.vincr[optimize.nparm] = 0.5; 00872 optimize.nparm += 1; 00873 } 00874 } 00875 00876 else if( strcmp(chKey4,"GLOB") == 0 ) 00877 { 00878 /* globule with density increasing inward 00879 * parameters are outer density, radius of globule, and density power */ 00880 ParseGlobule(chCard); 00881 } 00882 00883 else if( (strcmp(chKey4,"GRAI") == 0 )||( strcmp(chKey4,"PGRA") == 0 ) ) 00884 { 00885 /* read parameters dealing with grains, in reads2 */ 00886 ParseGrain(chCard,&lgDSet); 00887 } 00888 00889 else if( strcmp(chKey4,"GRID") == 0 ) 00890 { 00891 /* option to run grid of models by varying certain parameters 00892 * in reads2 */ 00893 ParseGrid(chCard); 00894 } 00895 00896 else if( strcmp(chKey4,"HDEN") == 0 ) 00897 { 00898 /* parse the hden command to set the hydrogen density, in reads2 */ 00899 ParseHDEN(chCard); 00900 } 00901 00902 else if( strcmp(chKey4,"HELI") == 0 ) 00903 { 00904 fprintf(ioQQQ,"Sorry, this command is replaced with ATOM HE-LIKE\n"); 00905 cdEXIT(EXIT_FAILURE); 00906 } 00907 00908 else if( strcmp(chKey4,"HEXT") == 0 ) 00909 { 00910 /* "extra" heating rate, so that first= log(erg(cm-3, s-1)) 00911 * second optional number is scale radius, so that HXTOT = TurbHeat*SEXP(DEPTH/SCALE) 00912 * if missing then constant heating. 00913 * third option is depth from shielded face, to mimic irradiation from both sides*/ 00914 i = 5; 00915 hextra.TurbHeat = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 00916 if( lgEOL ) 00917 NoNumb( chCard ); 00918 /* save initial value in case reset with time dependent command */ 00919 hextra.TurbHeatSave = hextra.TurbHeat; 00920 00921 /* remember which type of scale dependence we will use */ 00922 const char *chHextraScale; 00923 /* these are optional numbers on depth or density option */ 00924 realnum par1=0. , par2=0.; 00925 long int nPar; 00926 if( nMatch("DEPT" , chCard ) ) 00927 { 00928 /* depend on depth */ 00929 hextra.lgHextraDepth = true; 00930 chHextraScale = "DEPTH"; 00931 /* optional scale radius */ 00932 a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00933 if( lgEOL ) 00934 NoNumb(chCard); 00935 else 00936 hextra.turrad = (realnum)pow((realnum)10.f,a); 00937 00938 /* depth from shielded face, to mimic illumination from both sides 00939 * may not be specified */ 00940 realnum b = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00941 if( lgEOL ) 00942 { 00943 hextra.turback = 0.; 00944 nPar = 2; 00945 } 00946 else 00947 { 00948 hextra.turback = (realnum)pow((realnum)10.f,b); 00949 nPar = 3; 00950 } 00951 par1 = a; 00952 par2 = b; 00953 } 00954 else if( nMatch("DENS" , chCard ) ) 00955 { 00956 /* depends on density */ 00957 chHextraScale = "DENSITY"; 00958 hextra.lgHextraDensity = true; 00959 00960 /* the log of the density */ 00961 a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00962 /* if number not entered then unity is used, heating 00963 * proportional to density */ 00964 hextra.HextraScaleDensity = (realnum)pow((realnum)10.f,a); 00965 par1 = a; 00966 par2 = 0; 00967 nPar = 2; 00968 } 00969 else 00970 { 00971 /* constant heating */ 00972 chHextraScale = ""; 00973 nPar = 1; 00974 } 00975 00976 /* option to have heating vary with time */ 00977 if( nMatch( "TIME" , chCard ) ) 00978 hextra.lgTurbHeatVaryTime = true; 00979 00980 /* vary option */ 00981 if( optimize.lgVarOn ) 00982 { 00983 /* 1 to 3 options on hextra vary */ 00984 optimize.nvarxt[optimize.nparm] = nPar; 00985 optimize.vparm[0][optimize.nparm] = log10(hextra.TurbHeat); 00986 optimize.vparm[1][optimize.nparm] = par1; 00987 optimize.vparm[2][optimize.nparm] = par2; 00988 00989 strcpy( optimize.chVarFmt[optimize.nparm], "HEXTra %f " ); 00990 strcat( optimize.chVarFmt[optimize.nparm], chHextraScale ); 00991 while( nPar > 1 ) 00992 { 00993 /* add extra spots to write parameters */ 00994 --nPar; 00995 strcat( optimize.chVarFmt[optimize.nparm], " %f" ); 00996 } 00997 if( hextra.lgTurbHeatVaryTime ) 00998 strcat( optimize.chVarFmt[optimize.nparm], " TIME" ); 00999 /* pointer to where to write */ 01000 optimize.nvfpnt[optimize.nparm] = input.nRead; 01001 /* the increment in the first steps away from the original value */ 01002 optimize.vincr[optimize.nparm] = 0.1f; 01003 optimize.nparm += 1; 01004 } 01005 } 01006 01007 else if( strcmp(chKey4,"HIGH") == 0 ) 01008 { 01009 /* approach equilibrium from high te */ 01010 thermal.lgTeHigh = true; 01011 } 01012 01013 else if( strncmp( chCard ,"HYDROGEN",8) == 0 ) 01014 { 01015 fprintf(ioQQQ," Sorry, this command has been replaced with the ATOM H-LIKE command.\n"); 01016 cdEXIT(EXIT_FAILURE); 01017 } 01018 01019 else if( strcmp(chKey4,"ILLU") == 0 ) 01020 { 01021 /* illumination angle */ 01022 i = 5; 01023 geometry.AngleIllumRadian = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01024 if( lgEOL ) 01025 { 01026 NoNumb(chCard); 01027 } 01028 /* default is degrees, but radian key also there */ 01029 if( nMatch("RADI",chCard) ) 01030 { 01031 /* entered as radians, convert to degrees first */ 01032 geometry.AngleIllumRadian /= (realnum)(PI/180.); 01033 } 01034 /* we now have degrees, make sure between 0 and 90 */ 01035 if( geometry.AngleIllumRadian < 0. || geometry.AngleIllumRadian >= 90. ) 01036 { 01037 fprintf( ioQQQ, " Angle of illumination must be between 0 and 90 degrees " 01038 "or 0 and pi/2 radians.\n" ); 01039 cdEXIT(EXIT_FAILURE); 01040 } 01041 /* we really want 1. / COS( theta ) with theta in radians 01042 * first convert angle to radians */ 01043 geometry.AngleIllumRadian = (realnum)(geometry.AngleIllumRadian*PI/180.); 01044 01045 /* this is effective path - dTau_eff = dTau_normal * geometry.DirectionalCosin */ 01046 geometry.DirectionalCosin = (realnum)(1./cos(geometry.AngleIllumRadian)); 01047 01048 /* vary option */ 01049 if( optimize.lgVarOn ) 01050 { 01051 /* no luminosity options on vary */ 01052 optimize.nvarxt[optimize.nparm] = 1; 01053 strcpy( optimize.chVarFmt[optimize.nparm], "ILLUminate %f radians " ); 01054 /* pointer to where to write */ 01055 optimize.nvfpnt[optimize.nparm] = input.nRead; 01056 optimize.vparm[0][optimize.nparm] = geometry.AngleIllumRadian; 01057 /* the increment in the first steps away from the original value */ 01058 optimize.vincr[optimize.nparm] = 0.1f; 01059 optimize.nparm += 1; 01060 } 01061 } 01062 01063 else if( strcmp(chKey4,"INIT") == 0 ) 01064 { 01065 /* read cloudy.ini initialization file 01066 * need to read in file every time, since some vars 01067 * get reset in zero - would be unsafe not to read in again */ 01068 ParseInit(chCard); 01069 01070 /* check that only one init file was in the input stream - 01071 * we cannot now read more than one */ 01072 ++nInitFile; 01073 if( nInitFile > 1 ) 01074 { 01075 fprintf( ioQQQ, 01076 " This is the second init file, I can only handle one.\nSorry.\n" ); 01077 cdEXIT(EXIT_FAILURE); 01078 } 01079 01080 /* we will put the data for the ini file at the end of the array of 01081 * line images, this is the increment for stuffing the lines in - negative */ 01082 input.iReadWay = -1; 01083 01084 /* initialize array reader, this sub does nothing but set 01085 * initial value of a variable, depending on value of iReadWay */ 01086 input_init(); 01087 } 01088 01089 else if( strcmp(chKey5,"INTEN") == 0 ) 01090 { 01091 /* intensity of this continuum source */ 01092 i = 5; 01093 if( nqh >= LIMSPC ) 01094 { 01095 /* too many continua were entered */ 01096 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 01097 cdEXIT(EXIT_FAILURE); 01098 } 01099 01100 /* silly, but calms down the lint */ 01101 ASSERT( nqh < LIMSPC ); 01102 strcpy( rfield.chRSpec[nqh], "SQCM" ); 01103 rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01104 if( lgEOL ) 01105 { 01106 NoNumb(chCard); 01107 } 01108 01109 /* set radius to very large value if not already set */ 01110 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 01111 if( !radius.lgRadiusKnown ) 01112 { 01113 /* set radius to large value */ 01114 ar1 = (realnum)radius.rdfalt; 01115 radius.Radius = pow(10.,radius.rdfalt); 01116 } 01117 if( nMatch("LINE",chCard) ) 01118 { 01119 /* silly, but calms down the lint */ 01120 ASSERT( nqh < LIMSPC ); 01121 /* option for linear input parameter */ 01122 rfield.totpow[nqh] = log10(rfield.totpow[nqh]); 01123 } 01124 strcpy( rfield.chSpNorm[nqh], "LUMI" ); 01125 /* ParseRangeOption in readsun */ 01126 ParseRangeOption(nqh,chCard); 01127 01128 /* >>chng 06 mar 22, add time option to vary only some continua with time */ 01129 if( nMatch( "TIME" , chCard ) ) 01130 rfield.lgTimeVary[nqh] = true; 01131 01132 /* vary option */ 01133 if( optimize.lgVarOn ) 01134 { 01135 /* create the format we will use to write the command */ 01136 /* >>chng 03 apr 30, had used log for range option, but if one or other 01137 * was 1 ryd, then became 0 in the log, which was misinterpreted as 01138 * one of the continuum bounds */ 01139 /*strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY%f log range %f %f" );*/ 01140 strcpy( optimize.chVarFmt[optimize.nparm], "INTENSITY %f range %f %f " ); 01141 /* array index for this command */ 01142 optimize.nvfpnt[optimize.nparm] = input.nRead; 01143 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh]; 01144 optimize.vincr[optimize.nparm] = 0.5; 01145 /* range option, but cannot be varied */ 01146 /*optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.range[nqh][0]); 01147 optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.range[nqh][1]);*/ 01148 /* >>change to linear option */ 01149 optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0]; 01150 optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1]; 01151 optimize.nvarxt[optimize.nparm] = 3; 01152 ++optimize.nparm; 01153 } 01154 ++nqh; 01155 } 01156 01157 else if( strcmp(chKey5,"INTER") == 0 ) 01158 { 01159 /* interpolate on input tables for continuum, set of power laws used 01160 * input ordered pairs nu( ryd or log(Hz) ), log(fnu) 01161 * additional lines begin CONTINUE 01162 * first check that this is the one and only INTERP command 01163 * in readsun */ 01164 ParseInterp(chCard,&lgEOF); 01165 } 01166 01167 else if( strcmp(chKey4,"IONI") == 0 ) 01168 { 01169 /* inter ionization parameter U=Q/12 R*R N C; 01170 * defined per hydrogen, not per electron (as before) 01171 * radius must also be entered if spherical, but not needed if plane */ 01172 ParseIonPar(&nqh,chCard,&ar1); 01173 } 01174 01175 else if( strcmp(chKey4,"ITER") == 0 ) 01176 { 01177 /* iterate to get optical depths and diffuse field properly */ 01178 i = 5; 01179 iterations.itermx = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) - 1; 01180 iterations.itermx = MAX2(iterations.itermx,1); 01181 /* >>chng 06 mar 19, malloc itrdim arrays 01182 * iterations.iter_malloc is -1 before space allocated and set to 01183 * itermx if not previously set */ 01184 if( iterations.itermx > iterations.iter_malloc - 1 ) 01185 { 01186 long int iter_malloc_save = iterations.iter_malloc; 01187 /* add one for mismatch between array dim and iterations defn, 01188 * another few for safety */ 01189 iterations.iter_malloc = iterations.itermx+3; 01190 iterations.IterPrnt = (long int*)REALLOC(iterations.IterPrnt , 01191 (size_t)iterations.iter_malloc*sizeof(long int) ); 01192 geometry.nend = (long int*)REALLOC(geometry.nend , 01193 (size_t)iterations.iter_malloc*sizeof(long int) ); 01194 radius.router = (double*)REALLOC(radius.router , 01195 (size_t)iterations.iter_malloc*sizeof(double) ); 01196 /* >>chng 06 jun 07, need to set IterPrnt, router, and nend */ 01197 for(j=iter_malloc_save; j<iterations.iter_malloc; ++j ) 01198 { 01199 iterations.IterPrnt[j] = iterations.IterPrnt[iter_malloc_save-1]; 01200 geometry.nend[j] = geometry.nend[iter_malloc_save-1]; 01201 radius.router[j] = radius.router[iter_malloc_save-1]; 01202 } 01203 } 01204 01205 if( nMatch("CONV",chCard) ) 01206 { 01207 /* option to keep iterating until it converges, checks on convergence 01208 * are done in update, and checked again in prtcomment */ 01209 conv.lgAutoIt = true; 01210 /* above would have been legitimate setting of ITERMX, set to default 10 */ 01211 if( lgEOL ) 01212 { 01213 iterations.itermx = 10 - 1; 01214 } 01215 a = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01216 /* change convergence criteria, preset in zero */ 01217 if( !lgEOL ) 01218 { 01219 conv.autocv = a; 01220 } 01221 } 01222 } 01223 01224 else if( strcmp(chKey4,"L(NU") == 0 ) 01225 { 01226 /* this is the specific luminosity at nu 01227 * following says n(nu) not nuF(nu) */ 01228 lgNu2 = false; 01229 /* in reads2 */ 01230 ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2); 01231 } 01232 01233 else if( strcmp(chKey4,"LASE") == 0 ) 01234 { 01235 /* mostly one frequency (a laser) to check gamma's,*/ 01236 /* say the continuum type */ 01237 strcpy( rfield.chSpType[rfield.nspec], "LASER" ); 01238 01239 /* scan off the laser's energy */ 01240 i = 5; 01241 rfield.slope[rfield.nspec] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL); 01242 01243 /* negative energies are logs */ 01244 if( rfield.slope[rfield.nspec] <= 0. ) 01245 { 01246 rfield.slope[rfield.nspec] = 01247 pow(10.,rfield.slope[rfield.nspec]); 01248 } 01249 if( lgEOL ) 01250 { 01251 NoNumb(chCard); 01252 } 01253 01254 /* next number is optional relative width of laser */ 01255 rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL); 01256 if( lgEOL ) 01257 { 01258 /* default width is 0.05 */ 01259 rfield.cutoff[rfield.nspec][0] = 0.05; 01260 } 01261 01262 /* increment counter making sure we still have space enough */ 01263 ++rfield.nspec; 01264 if( rfield.nspec >= LIMSPC ) 01265 { 01266 /* too many continua were entered */ 01267 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 01268 cdEXIT(EXIT_FAILURE); 01269 } 01270 } 01271 01272 else if( strcmp(chKey4,"LUMI") == 0 ) 01273 { 01274 /* luminosity of this continuum source */ 01275 i = 5; 01276 if( nqh >= LIMSPC ) 01277 { 01278 /* too many continua were entered */ 01279 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 01280 cdEXIT(EXIT_FAILURE); 01281 } 01282 01283 strcpy( rfield.chRSpec[nqh], "4 PI" ); 01284 rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01285 if( lgEOL ) 01286 { 01287 NoNumb(chCard); 01288 } 01289 if( nMatch("LINE",chCard) ) 01290 { 01291 /* option for linear input parameter */ 01292 rfield.totpow[nqh] = log10(rfield.totpow[nqh]); 01293 } 01294 01295 strcpy( rfield.chSpNorm[nqh], "LUMI" ); 01296 01297 /* the solar option - number is log total solar luminosity */ 01298 if( nMatch("SOLA",chCard) ) 01299 { 01300 /* option to use log of total luminosity in solar units */ 01301 rfield.range[nqh][0] = rfield.emm; 01302 rfield.range[nqh][1] = rfield.egamry; 01303 /* log of solar luminosity in erg s-1 */ 01304 rfield.totpow[nqh] += 33.5827f; 01305 } 01306 else 01307 { 01308 /* check if range is present and parse it if it is - ParseRangeOption in readsun 01309 * this includes TOTAL keyword for total luminosity */ 01310 ParseRangeOption(nqh,chCard); 01311 } 01312 01313 /* >>chng 06 mar 22, add time option to vary only some continua with time */ 01314 if( nMatch( "TIME" , chCard ) ) 01315 rfield.lgTimeVary[nqh] = true; 01316 01317 /* vary option */ 01318 if( optimize.lgVarOn ) 01319 { 01320 /* >>chng 03 apr 30, had used log for range option, but if one or other 01321 * was 1 ryd, then became 0 in the log, which was misinterpreted as 01322 * one of the continuum bounds */ 01323 strcpy( optimize.chVarFmt[optimize.nparm], "LUMINOSITY %f range %f %f " ); 01324 /* pointer to where to write */ 01325 optimize.nvfpnt[optimize.nparm] = input.nRead; 01326 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh]; 01327 optimize.vincr[optimize.nparm] = 0.5; 01328 /* range option in, but cannot be varied */ 01329 optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0]; 01330 optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1]; 01331 optimize.nvarxt[optimize.nparm] = 3; 01332 optimize.nparm += 1; 01333 } 01334 ++nqh; 01335 } 01336 01337 else if( strcmp(chKey4,"MAGN") == 0 ) 01338 { 01339 /* parse the magnetic field command, routine in magnetic.c */ 01340 ParseMagnet( chCard ); 01341 } 01342 01343 else if( strcmp(chKey4,"MAP ") == 0 ) 01344 { 01345 /* do cooling space map for specified zones, in reads2 */ 01346 ParseMap(chCard); 01347 } 01348 01349 else if( strcmp(chKey4,"META") == 0 ) 01350 { 01351 /* read depletion for metals, all elements heavier than He 01352 * in reads2 */ 01353 ParseMetal(chCard); 01354 /* abundances no longer solar */ 01355 abund.lgAbnSolar = false; 01356 } 01357 01358 else if( strcmp(chKey4,"NEUT") == 0 ) 01359 { 01360 /* heating and ionization due to fast neutrons */ 01361 hextra.lgNeutrnHeatOn = true; 01362 01363 /* first number is neutron luminosity 01364 * relative to bolometric luminosity */ 01365 i = 5; 01366 hextra.frcneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01367 if( lgEOL ) 01368 NoNumb(chCard); 01369 01370 /* save as log of fraction */ 01371 if( hextra.frcneu > 0. ) 01372 hextra.frcneu = (realnum)log10(hextra.frcneu); 01373 01374 /* second number is efficiency */ 01375 hextra.effneu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01376 if( lgEOL ) 01377 { 01378 hextra.effneu = 1.0; 01379 } 01380 else 01381 { 01382 if( hextra.effneu <= 0. ) 01383 hextra.effneu = (realnum)pow((realnum)10.f,hextra.effneu); 01384 } 01385 } 01386 01387 else if( strcmp(chKey3,"NO ") == 0 ) 01388 { 01389 /* don't do something, in readsun */ 01390 ParseDont(chCard); 01391 } 01392 01393 else if( strcmp(chKey4,"NORM") == 0 ) 01394 { 01395 /* normalize lines to this rather than h-b, sec number is scale factor */ 01396 ParseNorm(chCard); 01397 } 01398 01399 else if( strcmp(chKey4,"NUF(") == 0 ) 01400 { 01401 /* flux density of this continuum source, at optional frequency 01402 * expressed as product nu*f_nu */ 01403 lgNu2 = true; 01404 /* in reads2 */ 01405 ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2); 01406 } 01407 01408 else if( strcmp(chKey4,"NUL(") == 0 ) 01409 { 01410 /* specific luminosity density of this continuum source, at opt freq 01411 * expressed as product nu*f_nu */ 01412 lgNu2 = true; 01413 /* in reads2 */ 01414 ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2); 01415 } 01416 01417 else if( strcmp(chKey4,"OPTI") == 0 ) 01418 { 01419 /* option to optimize the model by varying certain parameters 01420 * in reads2 */ 01421 ParseOptimize(chCard); 01422 } 01423 01424 else if( strcmp(chKey4,"PHI(") == 0 ) 01425 { 01426 /* enter phi(h), the number of h-ionizing photons/cm2 */ 01427 i = 5; 01428 if( nqh >= LIMSPC ) 01429 { 01430 /* too many continua were entered */ 01431 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 01432 cdEXIT(EXIT_FAILURE); 01433 } 01434 /* silly, but calms down the lint */ 01435 ASSERT( nqh < LIMSPC ); 01436 strcpy( rfield.chRSpec[nqh], "SQCM" ); 01437 strcpy( rfield.chSpNorm[nqh], "PHI " ); 01438 rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01439 if( lgEOL ) 01440 { 01441 NoNumb(chCard); 01442 } 01443 /* set radius to very large value if not already set */ 01444 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 01445 if( !radius.lgRadiusKnown ) 01446 { 01447 /* set radius to large value */ 01448 ar1 = (realnum)radius.rdfalt; 01449 radius.Radius = pow(10.,radius.rdfalt); 01450 } 01451 /* check if continuum so intense that crash is likely */ 01452 if( rfield.totpow[nqh] > 29. ) 01453 { 01454 fprintf( ioQQQ, " Is the flux for this continuum correct?\n" ); 01455 fprintf( ioQQQ, " It appears too bright to me.\n" ); 01456 } 01457 /* ParseRangeOption in readsun xx parse_rangeoption*/ 01458 ParseRangeOption(nqh,chCard); 01459 01460 /* >>chng 06 mar 22, add time option to vary only some continua with time */ 01461 if( nMatch( "TIME" , chCard ) ) 01462 rfield.lgTimeVary[nqh] = true; 01463 01464 /* vary option */ 01465 if( optimize.lgVarOn ) 01466 { 01467 /* >>chng 03 apr 30, had used log for range option, but if one or other 01468 * was 1 ryd, then became 0 in the log, which was misinterpreted as 01469 * one of the continuum bounds */ 01470 strcpy( optimize.chVarFmt[optimize.nparm], "phi(h) %f range %f %f" ); 01471 /* pointer to where to write */ 01472 optimize.nvfpnt[optimize.nparm] = input.nRead; 01473 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh]; 01474 optimize.vincr[optimize.nparm] = 0.5; 01475 /* range option in, but cannot be varied */ 01476 optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0]; 01477 optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1]; 01478 optimize.nvarxt[optimize.nparm] = 3; 01479 optimize.nparm += 1; 01480 } 01481 ++nqh; 01482 } 01483 01484 else if( strcmp(chKey4,"PLOT") == 0 ) 01485 { 01486 /* make plot of log(nu.f(nu)) vs log(nu) or opacity 01487 * in file plot */ 01488 ParsePlot(chCard); 01489 } 01490 01491 else if( strcmp(chKey4,"POWE") == 0 ) 01492 { 01493 /* power law with cutoff, in reads2 */ 01494 ParsePowerlawContinuum(chCard); 01495 } 01496 01497 else if( strcmp(chKey4,"PRIN") == 0 ) 01498 { 01499 /* adjust print schedule, in readsun */ 01500 ParsePrint(chCard); 01501 } 01502 01503 else if( strcmp(chKey4,"PUNC") == 0 ) 01504 { 01505 /* punch something, in punch */ 01506 ParsePunch(chCard); 01507 } 01508 01509 else if( strcmp(chKey4,"Q(H)") == 0 ) 01510 { 01511 /* log of number of ionizing photons */ 01512 i = 5; 01513 if( nqh >= LIMSPC ) 01514 { 01515 /* too many continua were entered */ 01516 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 01517 cdEXIT(EXIT_FAILURE); 01518 } 01519 01520 /* silly, but calms down the lint */ 01521 ASSERT( nqh < LIMSPC ); 01522 strcpy( rfield.chRSpec[nqh], "4 PI" ); 01523 strcpy( rfield.chSpNorm[nqh], "Q(H)" ); 01524 rfield.totpow[nqh] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01525 if( rfield.totpow[nqh] > 100. && called.lgTalk ) 01526 { 01527 fprintf( ioQQQ, " Is this reasonable?\n" ); 01528 } 01529 if( lgEOL ) 01530 { 01531 NoNumb(chCard); 01532 } 01533 /* ParseRangeOption in readsun */ 01534 ParseRangeOption(nqh,chCard); 01535 01536 /* >>chng 06 mar 22, add time option to vary only some continua with time */ 01537 if( nMatch( "TIME" , chCard ) ) 01538 rfield.lgTimeVary[nqh] = true; 01539 01540 /* vary option */ 01541 if( optimize.lgVarOn ) 01542 { 01543 /* >>chng 03 apr 30, had used log for range option, but if one or other 01544 * was 1 ryd, then became 0 in the log, which was misinterpreted as 01545 * one of the continuum bounds */ 01546 strcpy( optimize.chVarFmt[optimize.nparm], "Q(H) %f range %f %f" ); 01547 /* pointer to where to write */ 01548 optimize.nvfpnt[optimize.nparm] = input.nRead; 01549 optimize.vparm[0][optimize.nparm] = (realnum)rfield.totpow[nqh]; 01550 optimize.vincr[optimize.nparm] = 0.5; 01551 /* range option in, but cannot be varied */ 01552 optimize.vparm[1][optimize.nparm] = (realnum)rfield.range[nqh][0]; 01553 optimize.vparm[2][optimize.nparm] = (realnum)rfield.range[nqh][1]; 01554 optimize.nvarxt[optimize.nparm] = 3; 01555 ++optimize.nparm; 01556 } 01557 /* increment number of luminosity sources specified */ 01558 ++nqh; 01559 } 01560 01561 else if( strcmp(chKey4,"RATI") == 0 ) 01562 { 01563 /* enter a continuum luminosity as a ratio of 01564 * nuFnu for this continuum relative to a previous continuum 01565 * format; first number is ratio of second to first continuum 01566 * second number is energy for this ratio 01567 * if third number on line, then 2nd number is energy of 01568 * first continuum, while 3rd number is energy of second continuum 01569 * in reads2 */ 01570 ParseRatio(chCard,&nqh); 01571 } 01572 01573 else if( strcmp(chKey4,"RADI") == 0 ) 01574 { 01575 /* log of inner and outer radii, default second=infinity, 01576 * if R2<R1 then R2=R1+R2 01577 * there is an optional keyword, "PARSEC" on the line 01578 * to use PC as units, reads2 */ 01579 ParseRadius(chCard,&ar1); 01580 } 01581 01582 else if( strcmp(chKey4,"ROBE") == 0 ) 01583 { 01584 /* this is the Roberto Terlivich command, no telling if it still works */ 01585 radius.dRadSign = -1.; 01586 } 01587 01588 else if( strcmp(chKey4,"SET ") == 0 ) 01589 { 01590 /* set something, in reads2 */ 01591 ParseSet(chCard); 01592 } 01593 01594 else if( strcmp(chKey4,"SPEC") == 0 ) 01595 { 01596 /* special key, can do anything */ 01597 cdEXIT(EXIT_FAILURE); 01598 } 01599 01600 else if( strcmp(chKey4,"SPHE") == 0 ) 01601 { 01602 /* compute a spherical model, diffuse field from other side in 01603 * in reads2 */ 01604 ParseSphere(chCard); 01605 } 01606 01607 else if( strcmp(chKey4,"STAT") == 0 ) 01608 { 01609 /* either get or put the code's state as a file on disk */ 01610 ParseState(chCard); 01611 } 01612 01613 else if( strcmp(chKey4,"STOP") == 0 ) 01614 { 01615 /* stop model at desired zone, temperature, column density or tau-912 01616 * in readsun */ 01617 ParseStop(chCard); 01618 } 01619 01620 else if( strcmp(chKey4,"TABL") == 0 ) 01621 { 01622 /* interpolate on input tables for continuum, set of power laws used 01623 * input stored in big BLOCK data 01624 * first check that this is the one and only INTERP command 01625 * in readsun */ 01626 ParseTable(&nqh,chCard,&ar1); 01627 } 01628 01629 else if( strcmp(chKey4,"TAUM") == 0 ) 01630 { 01631 /* minimum optical depths for lines (normally 1e-20) 01632 * (used in STARTER to set lines up) */ 01633 i = 5; 01634 opac.taumin = (realnum)pow(10.,FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL)); 01635 if( lgEOL ) 01636 { 01637 NoNumb(chCard); 01638 } 01639 } 01640 01641 else if( strcmp(chKey4 , "TEST" ) == 0 ) 01642 { 01643 /* parse the test command and its options */ 01644 ParseTest(chCard,&nqh , &ar1, lgDSet); 01645 } 01646 01647 else if( strcmp(chKey4,"TIME") == 0 ) 01648 { 01649 /* parse the time dependent command, in dynamics.c */ 01650 ParseDynaTime(chCard); 01651 } 01652 01653 else if( strcmp(chKey4,"TITL") == 0 ) 01654 { 01655 /* read in title of model starting in col 5 */ 01656 strcpy( input.chTitle , &input.chOrgCard[5] ); 01657 } 01658 01659 else if( strcmp(chKey4,"TLAW") == 0 ) 01660 { 01661 /* some temperature vs depth law */ 01662 ParseTLaw(chCard); 01663 } 01664 01665 else if( strcmp(chKey4,"TOLE") == 0 ) 01666 { 01667 fprintf(ioQQQ, 01668 "Sorry, this command has been replaced with the SET TEMPERATURE TOLERANCE command.\n"); 01669 cdEXIT(EXIT_FAILURE); 01670 } 01671 01672 else if( strcmp(chKey4,"TRAC") == 0 ) 01673 { 01674 /* turn on trace output, in reads2 */ 01675 ParseTrace(chCard); 01676 } 01677 01678 else if( strcmp(chKey4,"VLAW") == 0 ) 01679 { 01680 01681 /* velocity power law as a function of radius. */ 01682 i = 5; 01683 DoppVel.TurbVelLaw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01684 01685 DoppVel.lgTurbLawOn = true; 01687 /* for now, insist upon non-positive power, 01688 * so that velocity decreases with increasing radius. */ 01689 ASSERT( DoppVel.TurbVelLaw <= 0.f ); 01690 } 01691 01692 else if( strcmp(chKey4,"TURB") == 0 ) 01693 { 01694 /* line has turbulent velocity in km/s */ 01695 i = 5; 01696 DoppVel.TurbVel = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01697 01698 /* include turbulent pressure in equation of state? 01699 * >>chng 06 mar 24, include turbulent pressure in gas equation of state by default, 01700 * but not if NO PRESSURE occurs */ 01701 if( nMatch(" NO ",chCard) && nMatch("PRES",chCard) ) 01702 { 01703 DoppVel.lgTurb_pressure = false; 01704 } 01705 else 01706 { 01707 /* default is to include turbulent pressure in equation of state */ 01708 DoppVel.lgTurb_pressure = true; 01709 } 01710 01711 /* true if no number on line - check for equipartition */ 01712 DoppVel.lgTurbEquiMag = false; 01713 01714 if( nMatch("EQUI",chCard) && nMatch("PART",chCard) ) 01715 { 01716 /* turbulence equipartition option */ 01717 DoppVel.lgTurbEquiMag = true; 01718 } 01719 01720 if( nMatch(" LOG",chCard) ) 01721 { 01722 DoppVel.TurbVel = (realnum)pow((realnum)10.f,DoppVel.TurbVel); 01723 } 01724 /* now convert from km/s to cm/s */ 01725 DoppVel.TurbVel *= 1e5; 01726 01727 /* this is optional F parameter from equation 34 of 01728 *>>refer pressure turb Heiles, C. & Troland, T.H. 2005, ApJ, 624, 773 01729 * turbulent energy density will be rho F v^2 / 2 */ 01730 DoppVel.Heiles_Troland_F = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01731 if( lgEOL ) 01732 { 01733 /* this is the default value of 3 for isotropic turbulence */ 01734 DoppVel.Heiles_Troland_F = 3.f; 01735 } 01736 01737 /* option to have turbulence be dissipative, keyword is dissipate, 01738 * argument is log of scale length in cm */ 01739 if( nMatch("DISS",chCard) ) 01740 { 01741 DoppVel.DispScale = (realnum)pow(10., FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL) ); 01742 if( lgEOL ) 01743 { 01744 NoNumb(chCard); 01745 } 01746 } 01747 01748 /* vary option */ 01749 if( optimize.lgVarOn ) 01750 { 01751 /* only one parameter to vary */ 01752 optimize.nvarxt[optimize.nparm] = 1; 01753 strcpy( optimize.chVarFmt[optimize.nparm], "TURBULENCE= %f " ); 01754 /* pointer to where to write */ 01755 optimize.nvfpnt[optimize.nparm] = input.nRead; 01756 /* turbulent velocity */ 01757 optimize.vparm[0][optimize.nparm] = DoppVel.TurbVel/1e5f; 01758 optimize.vincr[optimize.nparm] = 0.2f*optimize.vparm[0][optimize.nparm]; 01759 ++optimize.nparm; 01760 } 01761 01762 /* set initial turbulence */ 01763 DoppVel.TurbVelZero = DoppVel.TurbVel; 01764 } 01765 01766 else if( strcmp(chKey4,"WIND") == 0 ) 01767 { 01768 /* NB - advection and wind commands are now a single command */ 01769 /* parse the wind command, in dynamics.c */ 01770 ParseDynaWind(chCard); 01771 } 01772 01773 else if( strcmp(chKey2,"XI") == 0 ) 01774 { 01775 /* inter x-ray ionization parameter xi=L/r^2 n; 01776 * defined per hydrogen 01777 * radius must also be entered if spherical, but not needed if plane */ 01778 ParseIonPar(&nqh,chCard,&ar1); 01779 } 01780 01781 /* a line starting with these characters is a comment and so ok 01782 * >>chng 06 sep 04 use routine to check for comments */ 01783 else if( !lgInputComment(chCard) /*(chK1 != '#' && chK1 != '*') && chK1 != '%'*/ ) 01784 { 01785 /* end of keys - if we get here key is unrecognized---------- */ 01786 fprintf( ioQQQ, " Unrecognized command. Key=\"%4.4s\". This is routine ParseCommands.\n", 01787 chKey4 ); 01788 fprintf( ioQQQ, " The line image was==%s==\n", 01789 chCard ); 01790 fprintf( ioQQQ, " Sorry.\n" ); 01791 cdEXIT(EXIT_FAILURE); 01792 } 01793 01794 /* get next line image, this is tail of while loop that will 01795 * look for lgEOF or blank card */ 01796 input_readarray(chCard,&lgEOF); 01797 } 01798 01799 /*------------------------------------------------------------------- */ 01800 /* fall through - hit lgEOF or blank line */ 01801 01802 /* >>chng 00 mar 31, removed dummy call to ParseQuantHeat, PvH */ 01803 01804 for( i=0; i < INPUT_LINE_LENGTH; i++ ) 01805 { 01806 chCard[i] = ' '; 01807 } 01808 chCard[INPUT_LINE_LENGTH-1] = '\0'; 01809 01810 if( called.lgTalk ) 01811 { 01812 fprintf( ioQQQ, "%23c*%81c*\n", ' ', ' ' ); 01813 fprintf( ioQQQ, "%23c***********************************************************************************\n\n\n\n", ' ' ); 01814 } 01815 01816 /* this is an option to turn on trace printout on the nth 01817 * call from the optimizer - only allow trace if 01818 * this is the case and nOptimiz 1 below nTrOpt */ 01819 if( optimize.lgTrOpt ) 01820 { 01821 /* nTrOpt was set with "optimize trace" command, 01822 * is iteration to turn on trace */ 01823 if( optimize.nTrOpt != optimize.nOptimiz ) 01824 { 01825 trace.lgTrace = false; 01826 /* following overrides turning on trace elsewhere */ 01827 trace.lgTrOvrd = false; 01828 } 01829 else 01830 { 01831 trace.lgTrace = true; 01832 called.lgTalk = true; 01833 trace.lgTrOvrd = true; 01834 fprintf( ioQQQ, " READR turns on trace from optimize option.\n" ); 01835 } 01836 } 01837 01838 /* set density from DLAW command, must be done here since it may depend on later input commands */ 01839 if( strcmp(dense.chDenseLaw,"DLW1") == 0 ) 01840 { 01841 dense.gas_phase[ipHYDROGEN] = (realnum)dense_fabden(radius.Radius,radius.depth); 01842 01843 if( dense.gas_phase[ipHYDROGEN] <= 0. ) 01844 { 01845 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" ); 01846 cdEXIT(EXIT_FAILURE); 01847 } 01848 } 01849 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 ) 01850 { 01851 dense.gas_phase[ipHYDROGEN] = (realnum)dense_tabden(radius.Radius,radius.depth); 01852 01853 if( dense.gas_phase[ipHYDROGEN] <= 0. ) 01854 { 01855 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density set by DLAW must be > 0.\n" ); 01856 cdEXIT(EXIT_FAILURE); 01857 } 01858 } 01859 01860 /* start checks on parameters set properly - this begins with same line saying start .. */ 01861 01862 /* lgStop_not_enough_info says that not enough info for model, so stop 01863 * set true in following tests if anything missing */ 01864 lgStop_not_enough_info = false; 01865 lgStop = false; 01866 01867 /* check whether hydrogen density has been set - this value was set to 0 in zero */ 01868 if( dense.gas_phase[ipHYDROGEN] <= 0. ) 01869 { 01870 fprintf( ioQQQ, " PROBLEM DISASTER Hydrogen density MUST be specified.\n" ); 01871 lgStop_not_enough_info = true; 01872 lgStop = true; 01873 /* need to set something since used below - will abort 01874 * since lgStop is set */ 01875 dense.gas_phase[ipHYDROGEN] = 1.; 01876 } 01877 01878 radius.rinner = radius.Radius; 01879 01880 /* mass flux for wind model - used for mass conservation */ 01881 wind.emdot = dense.gas_phase[ipHYDROGEN]*wind.windv0; 01882 01883 /* set converge criteria - limit number of iterations and zones */ 01884 if( iterations.lgConverge_set) 01885 { 01886 iterations.itermx = MIN2( iterations.itermx , iterations.lim_iter ); 01887 for( j=0; j < iterations.iter_malloc; j++ ) 01888 { 01889 geometry.nend[j] = MIN2( geometry.nend[j] , iterations.lim_zone ); 01890 } 01891 } 01892 01893 /* NSAVE is number of lines saved, =0 if no commands entered */ 01894 if( input.nSave < 0 ) 01895 { 01896 fprintf( ioQQQ, " PROBLEM DISASTER No commands were entered - whats up?\n" ); 01897 cdEXIT(EXIT_FAILURE); 01898 } 01899 01900 /* iterate to convergence and wind models are mutually exclusive */ 01901 if( wind.windv0 > 0. && conv.lgAutoIt ) 01902 { 01903 if( called.lgTalk ) 01904 { 01905 fprintf( ioQQQ, " NOTE PROBLEM Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" ); 01906 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off\n\n\n" ); 01907 } 01908 conv.lgAutoIt = false; 01909 iterations.itermx = 0; 01910 } 01911 01912 /* iterate to convergence and case b do not make sense together */ 01913 /* WJH 22 May 2004: unless we are doing i-front dynamics (negative wind.windv0) */ 01914 if( opac.lgCaseB && conv.lgAutoIt && wind.windv0 >= 0. ) 01915 { 01916 if( called.lgTalk ) 01917 { 01918 fprintf( ioQQQ, " NOTE Case B is an artificial test, it makes no sense to converge this model.\n" ); 01919 fprintf( ioQQQ, " NOTE Iterate to convergence is turned off.\n\n\n" ); 01920 } 01921 conv.lgAutoIt = false; 01922 iterations.itermx = 0; 01923 } 01924 01925 /* specifying a density power and constant pressure makes no sense */ 01926 if( dense.DensityPower!=0. && strcmp( dense.chDenseLaw, "CPRE" )==0 ) 01927 { 01928 if( called.lgTalk ) 01929 { 01930 fprintf( ioQQQ, " NOTE Specifying both a density power law and constant pressure is impossible.\n" ); 01931 } 01932 lgStop = true; 01933 } 01934 01935 if( !rfield.lgIonizReevaluate && strcmp( dense.chDenseLaw, "CDEN" )!=0 ) 01936 { 01937 if( called.lgTalk ) 01938 { 01939 fprintf( ioQQQ, " NOTE NO REEVALUATE IONIZATION can only be used with constant density.\n" ); 01940 fprintf( ioQQQ, " NOTE Resetting to reevaluate ionization.\n\n" ); 01941 } 01942 rfield.lgIonizReevaluate = true; 01943 } 01944 01945 /* check if the combination of stopping column density and H density are physically plausible */ 01946 thickness = min( MIN3( StopCalc.tauend, StopCalc.colpls, StopCalc.colnut ), 01947 MIN3( StopCalc.col_h2, StopCalc.col_h2_nut, StopCalc.HColStop ) ); 01948 if( thickness < COLUMN_INIT ) 01949 { 01950 /* a stop column density was specified - check on physical thickness this corresponds to */ 01951 thickness /= (dense.gas_phase[ipHYDROGEN]*geometry.FillFac); 01952 /* >> chng 06 dec 21, don't complain if outer radius set small with `stop thickness' command. */ 01953 if( thickness > 1e25 && radius.router[0] > 1e25 ) 01954 { 01955 fprintf( ioQQQ, 01956 "NOTE The specified column density and hydrogen density correspond to a thickness of %.2e cm.\n", 01957 thickness); 01958 fprintf( ioQQQ, 01959 "NOTE This seems large to me.\n"); 01960 fprintf(ioQQQ,"NOTE a very large radius may cause overflow.\n\n"); 01961 } 01962 } 01963 01964 if( gv.lgDColOn && thermal.ConstGrainTemp>0 && called.lgTalk ) 01965 { 01966 /* warn if constant grain temperature but gas-grain thermal effects 01967 * are still included */ 01968 fprintf( ioQQQ, 01969 "NOTE The grain temperatures are set to a constant value with the " 01970 "CONSTANT GRAIN TEMPERATURE command, but " 01971 "energy exchange \n"); 01972 fprintf( ioQQQ, 01973 "NOTE is still included. The grain-gas heating-cooling will be incorrect. " 01974 "Consider turning off gas-grain collisional energy\n"); 01975 fprintf( ioQQQ, 01976 "NOTE exchange with the NO GRAIN GAS COLLISIONAL ENERGY EXCHANGE command.\n\n\n"); 01977 } 01978 01979 /* >>chng 04 nov 27, test on these two */ 01980 if( !rfield.lgDoLineTrans && rfield.lgOpacityFine ) 01981 { 01982 if( called.lgTalk ) 01983 { 01984 fprintf( ioQQQ, " NOTE NO LINE TRANSER set but fine opacities still computed.\n" ); 01985 fprintf( ioQQQ, " NOTE Turning off fine opacities.\n\n" ); 01986 } 01987 rfield.lgOpacityFine = false; 01988 } 01989 01990 /* >>chng 04 nov 27, test of these two */ 01991 if( h2.lgH2ON && (!rfield.lgDoLineTrans || !rfield.lgOpacityFine) ) 01992 { 01993 if( called.lgTalk ) 01994 { 01995 fprintf( ioQQQ, " NOTE Large H2 molecule turned on but line transfer and fine opacities are not.\n" ); 01996 fprintf( ioQQQ, " NOTE Turning on line transfer and fine opacities.\n\n" ); 01997 } 01998 rfield.lgOpacityFine = true; 01999 rfield.lgDoLineTrans = true; 02000 } 02001 02002 if( rfield.lgMustBlockHIon && !rfield.lgBlockHIon ) 02003 { 02004 /* one of the input continua needs to have H-ionizing radiation 02005 * blocked with extinguish command, but it was not done */ 02006 fprintf( ioQQQ, "\n NOTE\n" 02007 " NOTE One of the incident continuum is a form used when no H-ionizing radiation is produced.\n" ); 02008 fprintf( ioQQQ, " NOTE You must also include the EXTINGUISH command to make sure this is done.\n" ); 02009 fprintf( ioQQQ, " NOTE The EXTINGUISH command was not included.\n" ); 02010 fprintf( ioQQQ, " NOTE YOU MAY BE MAKING A BIG MISTAKE!!\n NOTE\n\n\n\n" ); 02011 } 02012 02013 /* if stop temp set below default then we are going into cold and possibly molecular 02014 * gas - check some parameters in this case */ 02015 if( called.lgTalk && (StopCalc.tend < phycon.TEMP_STOP_DEFAULT || 02016 /* thermal.ConstTemp def is zero, set pos when constant temperature is set */ 02017 (thermal.ConstTemp > 0. && thermal.ConstTemp < phycon.TEMP_STOP_DEFAULT ) ) ) 02018 { 02019 /* print warning if temperature set below default but cosmic rays not turned on */ 02020 /* >>chng 06 mar 02, do not print if molecules are off */ 02021 if( (hextra.cryden == 0.) && !co.lgNoCOMole && !hmi.lgNoH2Mole ) 02022 { 02023 fprintf( ioQQQ, "\n NOTE\n" 02024 " NOTE The simulation is going into neutral gas but cosmic rays are not included.\n" ); 02025 fprintf( ioQQQ, " NOTE Ion-molecule chemistry will not occur without a source of ionization.\n" ); 02026 fprintf( ioQQQ, " NOTE The chemistry network may collapse deep in molecular regions.\n" ); 02027 fprintf( ioQQQ, " NOTE Consider adding galactic background cosmic rays with the COSMIC RAYS BACKGROUND command.\n" ); 02028 fprintf( ioQQQ, " NOTE You may be making a BIG mistake.\n NOTE\n\n\n\n" ); 02029 } 02030 } 02031 02032 /* dense.gas_phase[ipHYDROGEN] is linear hydrogen density (cm-3) */ 02033 /* test for hydrogen density properly set has already been done above */ 02034 if( called.lgTalk && dense.gas_phase[ipHYDROGEN] < 1e-4 ) 02035 { 02036 fprintf( ioQQQ, " NOTE Is the entered value of the hydrogen density (%.2e) reasonable?\n", 02037 dense.gas_phase[ipHYDROGEN]); 02038 fprintf( ioQQQ, " NOTE It seems pretty low to me.\n\n\n" ); 02039 } 02040 else if( called.lgTalk && dense.gas_phase[ipHYDROGEN] > 1e15 ) 02041 { 02042 fprintf( ioQQQ, " NOTE Is this value of the hydrogen density reasonable?\n" ); 02043 fprintf( ioQQQ, " NOTE It seems pretty high to me.\n\n\n" ); 02044 } 02045 02046 /* is the model going to crash because of extreme density? */ 02047 if( called.lgTalk && !lgStop && !lgStop_not_enough_info ) 02048 { 02049 if( dense.gas_phase[ipHYDROGEN] < 1e-6 || dense.gas_phase[ipHYDROGEN] > 1e19 ) 02050 { 02051 fprintf( ioQQQ, " NOTE Simulation may crash because of extreme " 02052 "density. The value was %.2e\n\n" , 02053 dense.gas_phase[ipHYDROGEN] ); 02054 } 02055 } 02056 02057 if( rfield.nspec == 0 ) 02058 { 02059 fprintf( ioQQQ, " PROBLEM DISASTER No incident radiation field was specified - " 02060 "at least put in the CMB.\n" ); 02061 lgStop = true; 02062 lgStop_not_enough_info = true; 02063 } 02064 02065 if( nqh == 0 ) 02066 { 02067 fprintf( ioQQQ, " PROBLEM DISASTER Luminosity of continuum MUST be specified.\n" ); 02068 lgStop = true; 02069 lgStop_not_enough_info = true; 02070 } 02071 02072 /* Testing on 0 is safe since this it is initialized that way 02073 * only print comment if continuum has been specified, 02074 * if continuum not given then we are aborting anyway */ 02075 if( radius.Radius == 0. && rfield.nspec> 0) 02076 { 02077 fprintf( ioQQQ, " PROBLEM DISASTER Starting radius MUST be specified.\n" ); 02078 lgStop = true; 02079 lgStop_not_enough_info = true; 02080 } 02081 02082 if( rfield.nspec != nqh ) 02083 { 02084 fprintf( ioQQQ, " PROBLEM DISASTER There were not the same number of continuum shapes and luminosities entered.\n" ); 02085 lgStop = true; 02086 } 02087 02088 /* we only want to do this test on the first call to the command 02089 * parser - it will be called many more times but with no grid command 02090 * during the actual grid calculation */ 02091 static bool lgFirstPass = true; 02092 02093 /* the NO VARY option sets this flag, and can be used to turn off 02094 * the grid command, as well as the optimizer */ 02095 if( optimize.lgNoVary && grid.lgGrid ) 02096 { 02097 /* ignore grids */ 02098 grid.lgGrid = false; 02099 optimize.nparm = 0; 02100 grid.nGridCommands = 0; 02101 } 02102 02103 if( lgFirstPass && grid.lgGrid && (optimize.nparm!=grid.nGridCommands) ) 02104 { 02105 /* number of grid vary options do match */ 02106 fprintf( ioQQQ, " PROBLEM DISASTER The GRID command was entered " 02107 "but there were %li GRID commands and %li commands with a VARY option.\n" , 02108 grid.nGridCommands , optimize.nparm); 02109 fprintf( ioQQQ, " There must be the same number of GRIDs and VARY.\n" ); 02110 lgStop = true; 02111 } 02112 lgFirstPass = false; 02113 02114 if( lgStop_not_enough_info ) 02115 { 02116 fprintf( ioQQQ, " PROBLEM DISASTER I do not have enough information to do the simulation, I cannot go on.\n" ); 02117 fprintf( ioQQQ, "\n\n Sorry.\n\n\n" ); 02118 cdEXIT(EXIT_FAILURE); 02119 } 02120 02121 if( lgStop ) 02122 { 02123 cdEXIT(EXIT_FAILURE); 02124 } 02125 02126 /* end checks on parameters set properly - this begins with same line saying start .. */ 02127 return; 02128 }