/*ParseCommands main command line parser, decode command, then call other routines to read */ #include #include #include "cddefines.h" #include "physconst.h" #include "varypar.h" #include "grainvar.h" #include "incidspec.h" #include "abnslr.h" #include "zeronum.h" #include "supdie.h" #include "elecden.h" #include "sphere.h" #include "plot.h" #include "trovrd.h" #include "cludon.h" #include "bounds.h" #include "physok.h" #include "double.h" #include "trace.h" #include "forcte.h" #include "casbhs.h" #include "called.h" #include "extinc.h" #include "toler.h" #include "taumin.h" #include "wind.h" #include "turhet.h" #include "hextra.h" #include "cyclot.h" #include "limfal.h" #include "itercnt.h" #include "autoit.h" #include "turb.h" #include "timed.h" #include "caseb.h" #include "thigh.h" #include "radius.h" #include "input.h" #include "pressure.h" #include "assertresults.h" #include "phycon.h" #include "filfac.h" #include "fudgec.h" #include "redis.h" #include "date.h" #include "neutrn.h" #include "dtrans.h" #include "cextra.h" #include "vrrfit.h" #include "betaver.h" #include "angllum.h" #include "caps.h" #include "ffmtread.h" #include "lgmatch.h" #include "readar.h" #include "rdinit.h" #include "nonumb.h" #include "toomancon.h" #include "tabden.h" #include "fabden.h" #include "readck.h" #include "parse.h" void ParseCommands(void) { char chCard[81], chKey3[4], chKey4[5], chKey5[6]; int lgDSet, lgEOF, lgEOL, lgNu2; char chK1; long int i, j, nqh; long nInitFile=0;/* used to count number of init files read in */ double var; float a, ar1, teset, z; static int lgNoInit = TRUE; # ifdef DEBUG_FUN fputs( "<+>ParseCommands()\n", debug_fp ); # endif /* following says abundances are solar */ abnslr.lgAbnSolar = TRUE; /* there are no plots desired yet */ plotCom.lgPlotON = FALSE; plotCom.nplot = 0; /* this flag remembers whether grains have ever been turned on. It is passed * to routine ParseAbundances - there grains will be turned on with commands * such as abundances ism, unless grains were previously set * with a grains command. this way abundances will not clobber explictly set * grains. */ lgDSet = FALSE; radius.Radius = 0.; radius.rinner = 0.; nqh = 0; IncidSpec.nspec = 0; /* initialize some code variables in case assert command used in input stream */ InitAssertResults(); for( i=0; i < LIMSPC; i++ ) { strcpy( spnorm.chRSpec[i], "UNKN" ); } VaryPar.nparm = 0; /* this is an option to turn on trace printout on the nth * call from the optimizer */ if( VaryPar.lgTrOpt ) { /* nTrOpt was set with "optimize trace" command, * is iteration to turn on trace */ if( VaryPar.nTrOpt == VaryPar.nOptimiz ) { trace.lgTrace = TRUE; called.lgTalk = TRUE; trovrd.lgTrOvrd = TRUE; fprintf( ioQQQ, " READR turns on trace from optimize option.\n" ); } } /* say this is a beta version if we are talking and it is the truth */ if( BetaVer.nBetaVer > 0 && called.lgTalk ) { fprintf( ioQQQ, "\n This is %s (beta %ld) of Cloudy, and is intended for testing only.\n\n", date.chVersion, BetaVer.nBetaVer ); } if( called.lgTalk ) { /* this code prints pretty lines at top of output box */ for( i=0; i<57 ; i++ ) { fprintf(ioQQQ,"%c",' '); } fprintf( ioQQQ, "cCloudy %7.7s\n\n", date.chVersion); for( i=0; i<23 ; i++ ) { fprintf(ioQQQ,"%c",' '); } /* now print box and date of version, before printing commands */ fprintf( ioQQQ, "**************************************"); fprintf( ioQQQ, "%7.7s", date.chDate); fprintf( ioQQQ, "**************************************\n"); for( i=0; i<23 ; i++ ) { fprintf(ioQQQ,"%c",' '); } fprintf( ioQQQ, "*"); for( i=0; i<81 ; i++ ) { fprintf(ioQQQ,"%c",' '); } fprintf( ioQQQ, "*\n"); } /* read in commands and print them */ /* following signals which read is in progress, * 1 is forward read, of input commands * -1 is reverse read from bottom of line array, of cloudy.ini file */ input.iReadWay = 1; /* initialize array reader, this sub does nothing but set * initial value of a variable, depending on value of iReadWay */ rdinit(); /* read until eof or blank line, then exit */ /* * readar puts the next stored line image into chCard * it will be <=80 char long, with eol at end * it will be in mixed upper and lower case * lgEOF is set true if we hit eof and no more lines present * code moving over to info contained in structures, * * input.chOrgCard == original image of command line * input.chCARDCAPS == original image converted to caps */ readar(chCard,&lgEOF); /* line starting with blank is eof, this checks we are not at eof */ while( !lgEOF && chCard[0] != ' ' ) { /* echo the line before we cap it */ if( called.lgTalk ) { fprintf( ioQQQ, " * "); i=0; while( chCard[i]!='\0' ) { fprintf(ioQQQ,"%c",chCard[i]); ++i; } while( i<80 ) { fprintf(ioQQQ,"%c",' '); ++i; } fprintf(ioQQQ,"*\n"); } /* change chCard to all caps */ caps( chCard ); /* now set up several partial keys */ chK1 = chCard[0]; /* first three characters into chKey3 */ strncpy( chKey3 , chCard , 3 ); chKey3[3] = '\0'; /* first four characters into chKey4 */ strncpy( chKey4 , chCard , 4 ); chKey4[4] = '\0'; /* first four characters into chKey4 */ strncpy( chKey5 , chCard , 5 ); chKey5[5] = '\0'; /* check whether VARY is on line */ if( lgMatch("VARY",chCard) ) { VaryPar.lgVarOn = TRUE; if( VaryPar.nparm + 1 > LIMPAR ) { fprintf( ioQQQ, " Too many VARY lines entered; the limit is%4ld\n", LIMPAR ); puts( "[Stop in ParseCommands]" ); exit(1); } } else { VaryPar.lgVarOn = FALSE; } /*if( strncmp(chCard,"C ",2) == 0 )*/ if( chCard[0]=='C' && (chCard[1]==' ' || chCard[1]== 0) ) { /* this is null test since lines starting with "C " are comments, * the char after the c can be a space or a newline */ i = 1; } /* start to look for keywords */ else if( strcmp(chKey4,"ABSO") == 0 ) { /* enter luminosity in absolute magnitudes, in reads2 */ ParseAbsMag(chCard,&nqh); } else if( strcmp(chKey3,"AGE") == 0 ) { /* enter age of cloud so we can check for time-steady reads2 */ ParseAge(chCard); } else if( strcmp(chKey4,"AGN ") == 0 ) { /* enter generic style agn continuum, in reads2 */ ParseAgn(chCard); } else if( strcmp(chKey4,"ABUN") == 0 ) { /* chemical abundances, readsun */ ParseAbundances(chCard,lgDSet); /* abundances no longer solar */ abnslr.lgAbnSolar = FALSE; } else if( strcmp(chKey4,"ASSE") == 0 ) { /* assert that code predicts certain results, in assertresults.h */ ParseAssertResults(); } else if( strcmp(chKey4,"ATOM") == 0 ) { /* parse the atom command, only option now is feii */ ParseAtom(chCard); } else if( strcmp(chKey4,"BACK") == 0 ) { /* cosmic background, in readsun */ ParseBackgrd(&nqh,chCard,&ar1); } else if( strcmp(chKey4,"BLAC") == 0 ) { /* black body, in reads2 */ ParseBlackbody(chCard,&nqh,&ar1); /* vary option */ if( VaryPar.lgVarOn ) { /* no luminosity options on vary */ VaryPar.nvarxt[VaryPar.nparm] = 1; strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "BLACKbody=%f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; /* log of temp will be pointer */ VaryPar.vparm[0][VaryPar.nparm] = (float)log10(IncidSpec.slope[IncidSpec.nspec-1]); VaryPar.vincr[VaryPar.nparm] = 0.5; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"BREM") == 0 ) { /* brems continuum from central object */ IncidSpec.nspec += 1; if( IncidSpec.nspec > LIMSPC ) TooManCon(); strcpy( spnorm.chSpType[IncidSpec.nspec-1], "BREMS" ); i = 5; IncidSpec.slope[IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i, 76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } /* temperature is interpreted as log if <=10, linear otherwise*/ if( IncidSpec.slope[IncidSpec.nspec-1] <= 10. ) { IncidSpec.slope[IncidSpec.nspec-1] = (float)pow(10.,IncidSpec.slope[IncidSpec.nspec-1]); } IncidSpec.cutoff[0][IncidSpec.nspec-1] = 0.; /* option for vary keyword */ if( VaryPar.lgVarOn ) { /* only one parameter */ VaryPar.nvarxt[VaryPar.nparm] = 1; strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "BREMS, T=%f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; /* log of temp will be pointer */ VaryPar.vparm[0][VaryPar.nparm] = (float)log10(IncidSpec.slope[IncidSpec.nspec-1]); VaryPar.vincr[VaryPar.nparm] = 0.5; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"CASE") == 0 ) /* do hydrogen case b */ { /* set flag saying we are doing case b */ caseb.lgCaseb = TRUE; /* scan in an optional optical depth in lya */ i = 5; tauminCom.tlamin = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { /* set tau-La to 10**9 if not specified */ tauminCom.tlamin = 1e9; } else { tauminCom.tlamin = (float)pow(10.,tauminCom.tlamin); } /* Hummer and Storey case B; no collisions from n=1, 2 (usually in) */ if( lgMatch("HUMM",chCard) ) { casbhs.lgCasBhs = TRUE; } } else if( strcmp(chKey4,"CEXT") == 0 ) { /* add "extra" cooling, power law temp dependence */ cextra.lgCExtraOn = TRUE; i = 5; cextra.CoolExtra = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); if( lgEOL ) { NoNumb(chCard); } cextra.cextpw = (float)FFmtRead(chCard,&i,76,&lgEOL); } else if( strcmp(chKey4,"COMP") == 0 ) { /* compile ascii version of stellar atmosphere continua in volk */ ParseCompile(chCard); } else if( strcmp(chKey4,"CONS") == 0 ) { /* constant temperature, pressure, density, or gas pressure * in readsun */ ParseConstant(chCard); } else if( strcmp(chKey4,"CORO") == 0 ) { /* coronal equilibrium; set constant temperature to number on line * in readsun */ ParseCoronal(chCard,&nqh,&ar1); } else if( strcmp(chKey4,"COSM") == 0 ) { /* C.R. heating, log of density of rel. electrons */ i = 5; hextra.cryden = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); if( lgEOL ) { if( lgMatch("BACK",chCard) ) { /* galactic background cosmic ray density to produce * secondary ionization rate quoted by Tielens and Hollenbach */ /* hextra.cryden = 2e-9f;*/ /* >>chng 99 Jun 24, slight change to value * quoted by McKee astro-ph 9901370, this will produce a total * secondary ionization rate of 2.5e-17 s^-1, as tested in * tsuite secondary.in. If each ionization produces 2.4 eV of heat, * the background heating rate should be 9.6e-29 * n*/ hextra.cryden = 7.07e-9f; hextra.crtemp = 2.6e9f; } else { NoNumb(chCard); } } else { /* optional power law density */ hextra.crpowr = (float)FFmtRead(chCard,&i,76,&lgEOL); /* option to specify a temp for non-rel electrons */ hextra.crtemp = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { /* relavistic limit (Balbus and McKee) */ hextra.crtemp = 2.6e9; } else { var = pow(10.,hextra.crtemp); hextra.crtemp = (float)MIN2(var,2.6e9); } } } else if( strcmp(chKey4,"COVE") == 0 ) { /* covering factor for gas */ i = 5; /* The geometric covering factor accounts for how much of 4\pi is * covered by gas, and so linearly multiplies the predictd intensities */ sphere.covgeo = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } /* if negative then log, convert to linear */ if( sphere.covgeo <= 0. ) { sphere.covgeo = (float)pow(10.,sphere.covgeo); } if( sphere.covgeo > 1. ) { fprintf( ioQQQ, " A covering factor greater than 1 makes no physical sense. Sorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } /* radiative transfer covering factor will be equal to the geometric one */ sphere.covrt = sphere.covgeo; } else if( strcmp(chKey4,"CRAS") == 0 ) { /* div by 0 to get crash as check on FP environment */ if( lgMatch("ZERO",chCard) ) { fprintf(ioQQQ," I will now div by 0 to get crash. Hold on.\n"); fflush(ioQQQ); sphere.covgeo = (float)(1. / ZeroNum ); fprintf(ioQQQ," I am still alive - something is wrong, result is %e\n", ZeroNum); puts( "[Stop in ParseCommands]" ); exit(1); } /* make overflow to get crash as check on FP environment */ else if( lgMatch("OVER",chCard) ) { ar1 = 1e-20f; fprintf(ioQQQ," I will now make overflow to get crash. Hold on.\n"); fflush(ioQQQ); sphere.covgeo = (float)(DBL_MAX / ar1 ); fprintf(ioQQQ," I am still alive - something is wrong, the result was %e\n", sphere.covgeo); puts( "[Stop in ParseCommands]" ); exit(1); } /* assert impossible to get crash as check on environment */ else if( lgMatch("ASSE",chCard) ) { fprintf(ioQQQ," I will now make assert to get crash. Hold on.\n"); fflush(ioQQQ); assert( DBL_MAX < ZeroNum ); fprintf(ioQQQ," I am still alive - something is wrong.\n"); puts( "[Stop in ParseCommands]" ); exit(1); } /* assert ratios of zeros (NaN) to get crash as check on environment */ else if( lgMatch(" NAN",chCard) ) { ar1 = 0.; fprintf(ioQQQ," I will now div 0 by 0 to get crash. Hold on.\n"); fflush(ioQQQ); sphere.covgeo = (float)(ar1 / ZeroNum ); fprintf(ioQQQ," I am still alive - something is wrong, the result was %e\n", sphere.covgeo); puts( "[Stop in ParseCommands]" ); exit(1); } else { fprintf(ioQQQ," crash options are ZERO, OVERflow, _NaN, and ASSErt.\n"); puts( "[Stop in ParseCommands]" ); exit(1); } } else if( strcmp(chKey4,"CYLI") == 0 ) { /* height of cylinder in cm */ i = 5; radius.lgCylnOn = TRUE; radius.CylindHigh = pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); if( lgEOL ) { NoNumb(chCard); } } else if( strcmp(chKey4,"DIEL") == 0 ) { /* change various factors for dielectronic recombination */ if( lgMatch("KLUD",chCard) ) { /* option to turn off guess of diel rec coef for 3rd row elements */ i = 3; cludon.GuessDiel = (float)FFmtRead(chCard,&i,76,&lgEOL); } else if( lgMatch("BURG",chCard) ) { /* modify suppression of burgess dielectronic recombinations */ if( lgMatch(" ON ",chCard) ) { /* leave suppression on - this is default state */ supdie.lgSupDie[0] = TRUE; } else if( lgMatch(" OFF",chCard) ) { /* turn suppression off */ supdie.lgSupDie[0] = FALSE; } else { fprintf( ioQQQ, " flag ON or OFF must appear.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } } else if( lgMatch("NUSS",chCard) ) { /* modify suppression of nussbaumer and storey dielectronic recomb */ if( lgMatch(" ON ",chCard) ) { /* turn suppression on */ supdie.lgSupDie[1] = TRUE; } else if( lgMatch(" OFF",chCard) ) { /* leave suppression off - this is default state */ supdie.lgSupDie[1] = FALSE; } else { fprintf( ioQQQ, " flag ON or OFF must appear.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } } else { fprintf( ioQQQ, " key KLUDge, BURGess, or NUSSbaumer must appear.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } } else if( strcmp(chKey4,"DIFF") == 0 ) { /* set method of transferring diffuse fields */ if( lgMatch(" OTS",chCard) ) { strcpy( dtrans.chDffTrns, "OTS" ); } else if( lgMatch(" OUT",chCard) ) { i = 4; j = (long int)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { /* this is the default set in block data */ strcpy( dtrans.chDffTrns, "OU2" ); } else { if( j > 0 && j < 10 ) { sprintf( dtrans.chDffTrns, "OU%1ld", j ); } else { fprintf( ioQQQ, " must be between 1 and 9 \n" ); puts( "[Stop in ParseCommands]" ); exit(1); } } } else { fprintf( ioQQQ, " There should have been OUTward or OTS on this line. Sorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } } else if( strcmp(chKey4,"DLAW") == 0 ) { /* either use fabden routine, or read in table of points */ ParseDLaw(chCard); } else if( strcmp(chKey4,"DOUB") == 0 ) { /* double optical depth scale after each iteration */ double_.DoubleTau = 2.; } else if( strcmp(chKey4,"DRIV") == 0 ) { /* call one of several drivers, readsun */ ParseDriveCmnd(chCard); } else if( strcmp(chKey4,"GRAI") == 0 ) { /* read parameters dealing with grains, in reads2 */ ParseGrain(chCard,&lgDSet); } else if( strcmp(chKey4,"EDEN") == 0 ) { /* option to add "extra" electrons, to test compton limit * for very low T(star) - option is log of eden */ i = 5; ElecDen.EdenExtra = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); if( lgEOL ) { NoNumb(chCard); } /* warn that this model is meaningless */ physok.lgPhysOK = FALSE; } else if( strcmp(chKey4,"ELEM") == 0 ) { /* element command; * scale or abundance options, to change abundance of specific element * read option to change order of elements * in reads2.f */ ParseElement(chCard); } else if( strcmp(chKey4,"ENER") == 0 ) { /* energy density (degrees K) of this continuum source */ i = 5; nqh += 1; if( nqh > LIMSPC ) { TooManCon(); exit(1); } strcpy( spnorm.chRSpec[nqh-1], "SQCM" ); teset = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( radius.Radius == 0. ) { /* set radius to large value */ ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,ar1); } /* convert temp to log, recognize linear option */ if( lgMatch("LINE",chCard) || teset > 10. ) { /* option for linear temperature, must store log */ teset = (float)log10(teset); } if( teset > 5. ) { fprintf( ioQQQ, " This intensity may be too large. The code may crash due to overflow. Was log intended?\n" ); } /* teset is not log of temp, now get log of total luminosity */ strcpy( spnorm.chSpNorm[nqh-1], "LUMI" ); /* full range of continuum will be used */ IncidSpec.range[0][nqh-1] = bounds.emm; IncidSpec.range[1][nqh-1] = bounds.egamry; IncidSpec.totpow[nqh-1] = (float)(4.*teset - 4.2464476 + 0.60206); /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "ENERGY DENSITY %f log " ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = (float)log10(IncidSpec.totpow[nqh-1]); VaryPar.vincr[VaryPar.nparm] = 0.1f; VaryPar.nvarxt[VaryPar.nparm] = 1; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"ESCA") == 0 ) { /* option to change escape probabilities, in reads2 */ ParseEscP(chCard); } else if( strcmp(chKey4,"EXTI") == 0 ) { /* extinguish ionizing continuum by absorbing column AFTER * setting luminosity or Q(H). First number is the column * density (log), second number is leakage (def=0%) * last number is lowest energy (ryd), last two may be omitted * from right to left */ i = 5; extinc.excolm = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); if( lgEOL ) { NoNumb(chCard); } /* option to set leakage - default is 0. */ extinc.exleak = (float)FFmtRead(chCard,&i,76,&lgEOL); /* optional leakage is zero */ if( lgEOL ) { extinc.exleak = 0.; } /* negative leaks are logs */ if( extinc.exleak < 0. ) { extinc.exleak = (float)pow(10.,extinc.exleak); } if( extinc.exleak > 1. ) { /* but leaks greater than 1 are not allowed */ if( called.lgTalk ) { fprintf( ioQQQ, " A leakage of%9.0f%% was entered - this must be less than 100%%\n", extinc.exleak*100. ); } puts( "[Stop in ParseCommands]" ); exit(1); } /* option to set lowest energy for absorber */ extinc.exlow = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { extinc.exlow = 0.99946f; } else { if( extinc.exlow <= 0. ) extinc.exlow = (float)pow(10.,extinc.exlow); if( extinc.exlow < 0.99946 ) { fprintf( ioQQQ, " Energy less than 1 Ryd!!\n" ); } } } else if( strcmp(chKey4,"FAIL") == 0 ) { /* reset number of temp failures allowed, default=20 */ i = 5; /* save previous value */ j = limfal.LimFail; limfal.LimFail = (long int)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } /* 'no map' flag, fail without mapping */ if( lgMatch("NO M",chCard) ) { limfal.lgNoMap = TRUE; } /* complain if failures was increased above default */ if( limfal.LimFail > j ) { fprintf( ioQQQ, " Cloudy should not have problems making this command necessary.\n" ); fprintf( ioQQQ, " Please show this input stream to Gary Ferland if this command is really needed for this model.\n" ); } } else if( strcmp(chKey4,"FILL") == 0 ) { /* filling factor, power law exponent (default=1., 0.) */ i = 5; a = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( a <= 0. ) { /* number less than or equal to 0, must have been entered as log */ filfac.FillFac = (float)pow(10.,a); } else { /* number greater than zero, must have been the real thing */ filfac.FillFac = a; if( filfac.FillFac > 1.0 ) { if( called.lgTalk ) { fprintf( ioQQQ, " Filling factor > 1, reset to 1\n" ); } filfac.FillFac = 1.; } } filfac.fiscal = filfac.FillFac; /* option to have power law dependence, filpow set to 0 in zerologic */ filfac.filpow = (float)FFmtRead(chCard,&i,76,&lgEOL); /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "FILLING FACTOR= %f power=%f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = (float)log10(filfac.FillFac); VaryPar.vincr[VaryPar.nparm] = 0.5; /* power law dependence here, but cannot be varied */ VaryPar.vparm[1][VaryPar.nparm] = filfac.filpow; VaryPar.nvarxt[VaryPar.nparm] = 2; /* do not allow filling factor to go positive */ VaryPar.varang[VaryPar.nparm][0] = -1e38f; VaryPar.varang[VaryPar.nparm][1] = 0.f; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"FIRE") == 0 ) { /* cosmic thermal background radiation, argument is redshift */ i = 5; /* if no number on line then (float)FFmtRead returns z=0; i.e., now */ z = (float)FFmtRead(chCard,&i,76,&lgEOL); /* in readsun */ ParseFireBall(z,&nqh,&ar1); } else if( strcmp(chKey4,"FLUC") == 0 ) { /* rapid density fluctuations, in readsun */ ParseFluc(chCard); } else if( strcmp(chKey4,"F(NU") == 0 ) { /* this is the specific flux at nu * following says n(nu) not nuF(nu) */ lgNu2 = FALSE; /* in reads2 */ ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2); } else if( strcmp(chKey4,"FORC") == 0 ) { /* set temperature of first zone, but don't keep constant * log if < 10 */ i = 5; forcte.ForceTemp = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( lgMatch(" LOG",chCard) || (forcte.ForceTemp <= 10. && !lgMatch("LINE",chCard)) ) { forcte.ForceTemp = (float)pow(10.,forcte.ForceTemp); } /* low energy bounds of continuum array do not permit too-low a temp */ if( forcte.ForceTemp < 3. ) { fprintf( ioQQQ, " TE reset to 3K: entered number too small.\n" ); forcte.ForceTemp = 3.; } } else if( strcmp(chKey4,"FUDG") == 0 ) { /* enter a fudge factor */ i = 5; for( j=0; j < NFUDGC; j++ ) { fudgec.fudgea[j] = (float)FFmtRead(chCard,&i,76,&lgEOL); if( !lgEOL ) fudgec.nfudge = j+1; /* this will be fortran counting for now */ } } else if( strcmp(chKey4,"GLOB") == 0 ) { /* globule with density increasing inward * parameters are outer density, radius of globule, and density power */ ParseGlobule(chCard); } else if( strcmp(chKey4,"HDEN") == 0 ) { /* parse the hden command to set the hydrogen density, in reads2 */ ParseHDEN(chCard); } else if( strcmp(chKey4,"HELI") == 0 ) { /* stuff to do with helium, in reads2 */ ParseHelium(chCard); } else if( strcmp(chKey4,"HEXT") == 0 ) { /* turbulent heating rate, so that first= log(erg(cm-3, s-1)) * second is scale radius, so that HXTOT = TurbHeat*SEXP(TURRAD/ * 1 DEPTH) */ i = 5; turhet.TurbHeat = (float)pow(10.,FFmtRead(chCard,&i,76,&lgEOL)); a = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { turhet.turrad = 0.; } else { turhet.turrad = (float)pow(10.,a); } } else if( strcmp(chKey4,"HIGH") == 0 ) { /* approach equilibrium from high te */ thigh.lgTeHigh = TRUE; } else if( strncmp( chCard ,"HYDROGEN",8) == 0 ) { /* stuff to do with hydrogen atom, in reads2 */ ParseHydrogen(chCard); } else if( strcmp(chKey4,"ILLU") == 0 ) { /* illumination angle */ i = 5; Angllum.AngleIllum = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } /* default is degrees, but radian key also there */ if( lgMatch("RADI",chCard) ) { /* entered as radians, convert to degrees first */ Angllum.AngleIllum /= (float)(PI/180.); } /* we now have degrees, make sure between 0 and 90 */ if( Angllum.AngleIllum < 0. || Angllum.AngleIllum >= 90. ) { fprintf( ioQQQ, " Angle of illumination must be between 0 and 90.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } /* we really want 1. / COS( theta ) with theta in radians */ Angllum.AngleIllum = (float)(Angllum.AngleIllum*PI/180.); Angllum.AngleIllum = (float)(1./cos(Angllum.AngleIllum)); } else if( strcmp(chKey4,"INIT") == 0 ) { /* read cloudy.ini initialization file * following is set to false after first pass through this sub * so that init file only read one time in multi run jobs */ if( lgNoInit ) { ParseInit(chCard); } /* check that only one init file was in the input stream - * we cannot now read more than one */ ++nInitFile; if( nInitFile > 1 ) { fprintf( ioQQQ, " This is the second init file, I can only handle one.\nSorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } /* we will put the data for the ini file at the end of the array of * line images, this is the increment for stuffing the lines in - negative */ input.iReadWay = -1; /* initialize array reader, this sub does nothing but set * initial value of a variable, depending on value of iReadWay */ rdinit(); } else if( strcmp(chKey5,"INTEN") == 0 ) { /* intensity of this continuum source */ i = 5; nqh += 1; if( nqh > LIMSPC ) { TooManCon(); exit(1); } strcpy( spnorm.chRSpec[nqh-1], "SQCM" ); IncidSpec.totpow[nqh-1] = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( radius.Radius == 0. ) { /* set radius to large value */ ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,ar1); } if( lgMatch("LINE",chCard) ) { /* option for linear input parameter */ IncidSpec.totpow[nqh-1] = (float)log10(IncidSpec.totpow[nqh-1]); } strcpy( spnorm.chSpNorm[nqh-1], "LUMI" ); /* ParseRangeOption in readsun */ ParseRangeOption(nqh,chCard); /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "INTENSITY%f log range %f %f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = IncidSpec.totpow[nqh-1]; VaryPar.vincr[VaryPar.nparm] = 0.5; /* range option, but cannot be varied */ VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.range[0][nqh-1]); VaryPar.vparm[2][VaryPar.nparm] = (float)log10(IncidSpec.range[1][nqh-1]); VaryPar.nvarxt[VaryPar.nparm] = 3; VaryPar.nparm += 1; } } else if( strcmp(chKey5,"INTER") == 0 ) { /* interpolate on input tables for continuum, set of power laws used * input ordered pairs nu( ryd or log(hz) ), log(fnu) * additional lines begin CONTINUE * first check that this is the one and only INTERP command * in readsun */ ParseInterp(chCard,&lgEOF); } else if( strcmp(chKey4,"IONI") == 0 ) { /* inter ionization parameter U=Q/12 R*R N C; * defined per hydrogen, not per electron (as before) * radius must also be entered if spherical, but not needed if plane */ ParseIonPar(&nqh,chCard,&ar1); } else if( strcmp(chKey4,"ITER") == 0 ) { /* iterate to get optical depths and diffuse field properly */ i = 5; IterCnt.itermx = (long int)FFmtRead(chCard,&i,76,&lgEOL) - 1; IterCnt.itermx = MAX2(IterCnt.itermx,1); if( IterCnt.itermx > ITRDIM - 1 ) { fprintf( ioQQQ, " No more than%3ld can be performed because of vector sizes; reduce ITER or increase ITRDIM in code.\n", ITRDIM ); puts( "[Stop in ParseCommands]" ); exit(1); } if( lgMatch("CONV",chCard) ) { /* option to keep iterating until it converges, checks on convergence * are done in update, and checked again in prtcomment */ autoit.lgAutoIt = TRUE; /* above would have been legitimite setting of ITERMX, set to default 10 */ if( lgEOL ) { IterCnt.itermx = 10 - 1; } a = (float)FFmtRead(chCard,&i,76,&lgEOL); /* change convergence criteria, preset in zero */ if( !lgEOL ) { autoit.autocv = a; } } } else if( strcmp(chKey4,"L(NU") == 0 ) { /* this is the specific luminositiy at nu * following says n(nu) not nuF(nu) */ lgNu2 = FALSE; /* in reads2 */ ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2); } else if( strcmp(chKey4,"LASE") == 0 ) { /* mostly one frequency (a laser) to check gamma's */ IncidSpec.nspec += 1; if( IncidSpec.nspec > LIMSPC ) { TooManCon(); } strcpy( spnorm.chSpType[IncidSpec.nspec-1], "LASER" ); i = 5; IncidSpec.slope[IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i, 76,&lgEOL); if( IncidSpec.slope[IncidSpec.nspec-1] <= 0. ) { IncidSpec.slope[IncidSpec.nspec-1] = (float)pow(10.,IncidSpec.slope[IncidSpec.nspec-1]); } if( lgEOL ) { NoNumb(chCard); } } else if( strcmp(chKey4,"LUMI") == 0 ) { /* luminosity of this continuum source */ i = 5; nqh += 1; if( nqh > LIMSPC ) { TooManCon(); exit(1); } strcpy( spnorm.chRSpec[nqh-1], "4 PI" ); IncidSpec.totpow[nqh-1] = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( lgMatch("LINE",chCard) ) { /* option for linear input parameter */ IncidSpec.totpow[nqh-1] = (float)log10(IncidSpec.totpow[nqh-1]); } strcpy( spnorm.chSpNorm[nqh-1], "LUMI" ); if( lgMatch("SOLA",chCard) ) { /* option to use log of total lumin in solar units */ IncidSpec.range[0][nqh-1] = bounds.emm; IncidSpec.range[1][nqh-1] = bounds.egamry; IncidSpec.totpow[nqh-1] += 33.5827f; } else { /* ParseRangeOption in readsun */ ParseRangeOption(nqh,chCard); } /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "LUMINOSITY %f log range %f %f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = IncidSpec.totpow[nqh-1]; VaryPar.vincr[VaryPar.nparm] = 0.5; /* range option in, but cannot be varied */ VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.range[0][nqh-1]); VaryPar.vparm[2][VaryPar.nparm] = (float)log10(IncidSpec.range[1][nqh-1]); VaryPar.nvarxt[VaryPar.nparm] = 3; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"MAGN") == 0 ) { /* log of magnetic field (gauss) to turn on cyclotron cooling */ i = 5; cyclot.CycloCoolCoef = (float)pow(10.,2.*FFmtRead(chCard,&i,76, &lgEOL)); if( lgEOL ) { NoNumb(chCard); } /* coef is 4/3 /8pi /c * vtherm(elec) */ cyclot.CycloCoolCoef *= 4.5433e-25f; } else if( strcmp(chKey4,"MAP ") == 0 ) { /* do cooling space map for specified zones, in reads2 */ ParseMap(chCard); } else if( strcmp(chKey4,"META") == 0 ) { /* read depletion for metals, all elements heavier than He * in reads2 */ ParseMetal(chCard); /* abundances no longer solar */ abnslr.lgAbnSolar = FALSE; } else if( strcmp(chKey4,"NEUT") == 0 ) { /* heating and ionization due to fast neutrons, arg is total luminosity * in neutrons rel to boblmetric lumin; sec number is efficiency */ i = 5; neutrn.frcneu = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( neutrn.frcneu > 0. ) { neutrn.frcneu = (float)log10(neutrn.frcneu); } neutrn.lgNeutrnHeatOn = TRUE; neutrn.effneu = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { neutrn.effneu = 1.0; } else { if( neutrn.effneu <= 0. ) neutrn.effneu = (float)pow(10.,neutrn.effneu); } } else if( strcmp(chKey3,"NO ") == 0 ) { /* don't do something, in readsun */ ParseDont(chCard); } else if( strcmp(chKey4,"NORM") == 0 ) { /* normalise lines to this rather than h-b, sec number is scale factor */ ParseNorm(chCard); } else if( strcmp(chKey4,"NUF(") == 0 ) { /* flux density of this continuum source, at optional frequency * expressed as product nu*f_nu */ lgNu2 = TRUE; /* in reads2 */ ParseF_nu(chCard,&nqh,&ar1,"SQCM",lgNu2); } else if( strcmp(chKey4,"NUL(") == 0 ) { /* specific luminosity density of this continuum source, at opt freq * expressed as product nu*f_nu */ lgNu2 = TRUE; /* in reads2 */ ParseF_nu(chCard,&nqh,&ar1,"4 PI",lgNu2); } else if( strcmp(chKey4,"OPTI") == 0 ) { /* option to optimize the model by varying certain parameters * in reads2 */ ParseOptimize(chCard); } else if( strcmp(chKey4,"PHI(") == 0 ) { /* enter phi(h), the number of h-ionizing photons/cm2 */ i = 5; nqh += 1; if( nqh > LIMSPC ) { TooManCon(); exit(1); } strcpy( spnorm.chRSpec[nqh-1], "SQCM" ); strcpy( spnorm.chSpNorm[nqh-1], "PHI " ); IncidSpec.totpow[nqh-1] = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { NoNumb(chCard); } if( radius.Radius == 0. ) { /* set radius to large value */ ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,ar1); } /* check if continuum so intense that crash is likely */ if( IncidSpec.totpow[nqh-1] > 29. ) { fprintf( ioQQQ, " Is the flux for this continuum correct?\n" ); fprintf( ioQQQ, " It appears too bright to me.\n" ); } /* ParseRangeOption in readsun */ ParseRangeOption(nqh,chCard); /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "phi(h) %f log range %f %f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = IncidSpec.totpow[nqh-1]; VaryPar.vincr[VaryPar.nparm] = 0.5; /* range option in, but cannot be varied */ VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.range[0][nqh-1]); VaryPar.vparm[2][VaryPar.nparm] = (float)log10(IncidSpec.range[1][nqh-1]); VaryPar.nvarxt[VaryPar.nparm] = 3; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"PLOT") == 0 ) { /* make plot of log(nu.f(nu)) vs log(nu) or opacity * in file plot */ ParsePlot(chCard); } else if( strcmp(chKey4,"POWE") == 0 ) { /* power law with cutoff, in reads2 */ ParsePowerlawContinuum(chCard); } else if( strcmp(chKey4,"PRIN") == 0 ) { /* adjust print schedule, in readsun */ ParsePrint(chCard); } else if( strcmp(chKey4,"PUNC") == 0 ) { /* punch something, in punch */ ParsePunch(chCard); } else if( strcmp(chKey4,"Q(H)") == 0 ) { /* log of number of ionizing photons */ i = 5; nqh += 1; if( nqh > LIMSPC ) { TooManCon(); exit(1); } strcpy( spnorm.chRSpec[nqh-1], "4 PI" ); strcpy( spnorm.chSpNorm[nqh-1], "Q(H)" ); IncidSpec.totpow[nqh-1] = (float)FFmtRead(chCard,&i,76,&lgEOL); if( IncidSpec.totpow[nqh-1] > 100. && called.lgTalk ) { fprintf( ioQQQ, " Is this reasonable?\n" ); } if( lgEOL ) { NoNumb(chCard); } /* ParseRangeOption in readsun */ ParseRangeOption(nqh,chCard); /* vary option */ if( VaryPar.lgVarOn ) { strcpy( chVaryPar.chVarFmt[VaryPar.nparm], "Q(H) %f log range %f %f" ); /* pointer to where to write */ VaryPar.nvfpnt[VaryPar.nparm] = input.nRead; VaryPar.vparm[0][VaryPar.nparm] = IncidSpec.totpow[nqh-1]; VaryPar.vincr[VaryPar.nparm] = 0.5; /* range option in, but cannot be varied */ VaryPar.vparm[1][VaryPar.nparm] = (float)log10(IncidSpec.range[0][nqh-1]); VaryPar.vparm[2][VaryPar.nparm] = (float)log10(IncidSpec.range[1][nqh-1]); VaryPar.nvarxt[VaryPar.nparm] = 3; VaryPar.nparm += 1; } } else if( strcmp(chKey4,"QHEA") == 0 ) { /* (K Volk) added QHEAT here */ /* Turn on quantum heating heating (ignored if grains are not PAHs unless * they have the qheat keyword on the card) */ ParseQuantHeat(chCard,0); } else if( strcmp(chKey4,"RATI") == 0 ) { /* enter a continuum luminosity as a ratio of * nuFnu for this continuum relative to a previous continuum * format; first number is ratio of second to first continuum * second number is energy for this ratio * if third numbewr on line, then 2nd number is energy of * first continuum, while 3rd number is energy of second continuum * in reads2 */ ParseRatio(chCard,&nqh); } else if( strcmp(chKey4,"RADI") == 0 ) { /* log of inner and outer radii, default second=infinity, * if R2= 0. ) { fprintf( ioQQQ, " The initial wind velocity must be negative for D-critical flow. Sorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } /* gravity is not included in hydro flow*/ wind.comass = 0.; } /* this is usual hypersonic outflow */ else { if( wind.windv0 <= 0. ) { fprintf( ioQQQ, " The initial wind velocity must be positive. Sorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } else if( wind.windv0 <= 1.e6 ) { /* speed of sound roughly 10 km/s */ fprintf( ioQQQ, " >>>>Initial wind velocity should be greater than speed of sound; calculation only valid above sonic point.\n" ); wind.lgWindOK = FALSE; } if( lgEOL ) { NoNumb(chCard); } wind.comass = (float)FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { wind.comass = 1.; } /* used for line width, this is not reset with d-critical flow */ strcpy( redis.chRedis, "WIND" ); } wind.windv = wind.windv0; /* this is used in convpres to say wind solution - both cases use this*/ strcpy( pressure.chCPres, "WIND" ); /* end end end end end end end end end end end end end end end end * these are first col characters that mean line totally * ignored. usually stripped out in fillar, but could * be input here is code called as sub */ } else if( (chK1 != '#' && chK1 != '*') && chK1 != '%' ) { /* end of keys - if we get here key is unrecognized---------- */ fprintf( ioQQQ, " Unrecognized command. Key=\"%4.4s\". This is routine ParseCommands.\n", chKey4 ); fprintf( ioQQQ, " The line image was==%s==\n", chCard ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in ParseCommands]" ); exit(1); } /* get next line image, this is tail of while loop that will * look for lgEOF or blank card */ readar(chCard,&lgEOF); } /*------------------------------------------------------------------- */ /* fall through - hit lgEOL or blank line */ /* qtmax was set to zero in zero. * if qheat command ever appeared then qtmax is non-zero */ if( GrainVar.qtmax == 0. ) { /* there were no qheat commands, * but qheat may have been turned on with a grains line * generate dummy command line and call ParseQuantHeat to get defaults */ strcpy( chCard, "QHEAT " ); for( j=5; j < 80; j++ ) { chCard[j] = ' '; } chCard[80] = '\0'; /* flag 1 indicates which job to do */ ParseQuantHeat(chCard,1); } /* end of changes (K Volk) */ for( i=0; i < 80; i++ ) { chCard[i] = ' '; } chCard[80] = '\0'; if( called.lgTalk ) { fprintf( ioQQQ, " * %80.80s*\n", chCard ); fprintf( ioQQQ, " ***********************************************************************************\n\n\n\n" ); } /* set flag so that we will not actually read cards in again */ lgNoInit = FALSE; if( strcmp(pressure.chCPres,"DLW1") == 0 ) { phycon.hden = (float)log10(fabden(radius.Radius,radius.depth)); } else if( strcmp(pressure.chCPres,"DLW2") == 0 ) { /* >>chng 96 nov 29, added tabden, based on Kevin Volk code */ phycon.hden = (float)log10(tabden(radius.Radius,radius.depth)); } /* perform sanity checks for some of the things read in, in reads2 */ readck(nqh); /* fixup some variables */ phycon.hden = (float)pow(10.,phycon.hden); radius.rinner = radius.Radius; /* mass loss rate for wind model, normally WINDV0=0 */ wind.emdot = (float)(phycon.hden*wind.windv0); /* iterate to convergence and wind models are mutually exclusive */ if( wind.windv0 > 0. && autoit.lgAutoIt ) { fprintf( ioQQQ, " Due to the nature of the Sobolev approximation, it makes no sense to converge a windy model.\n" ); fprintf( ioQQQ, " Iterate to convergence is turned off\n" ); autoit.lgAutoIt = FALSE; IterCnt.itermx = 0; } /* this is an option to turn on trace printout on the nth * call from the optimizer - only allow trace if * this is the case and nOptimiz 1 below nTrOpt */ if( VaryPar.lgTrOpt ) { /* nTrOpt was set with "optimize trace" command, * is iteration to turn on trace */ if( VaryPar.nTrOpt != VaryPar.nOptimiz ) { trace.lgTrace = FALSE; /* following overrides turning on trace elsewhere */ trovrd.lgTrOvrd = FALSE; } else { trace.lgTrace = TRUE; called.lgTalk = TRUE; trovrd.lgTrOvrd = TRUE; fprintf( ioQQQ, " READR turns on trace from optimize option.\n" ); } } # ifdef DEBUG_FUN fputs( " <->ParseCommands()\n", debug_fp ); # endif return; }