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 /*ParseGrain parse parameters on grains command */ 00004 #include "cddefines.h" 00005 #include "grainvar.h" 00006 #include "phycon.h" 00007 #include "input.h" 00008 #include "optimize.h" 00009 #include "parse.h" 00010 #include "grains.h" 00011 00012 void ParseGrain(char *chCard, 00013 bool *lgDSet) 00014 { 00015 bool lgC15 = false, 00016 lgC120 = false, 00017 lgEOL, 00018 lgLinSet, 00019 lgLogLinSet, 00020 lgQuoteFound, 00021 lgSizeDistribution; 00022 char *ptr; 00023 GrainPar gp; 00024 00025 /*possible name of input file with opacities */ 00026 char chFile[FILENAME_PATH_LENGTH_2]; 00027 const char *chOption = NULL; 00028 long int i; 00029 00030 DEBUG_ENTRY( "ParseGrain()" ); 00031 00032 *lgDSet = true; 00033 00034 /* >>chng 00 dec 20, initialize chFile to empty string..., PvH */ 00035 chFile[0] = '\0'; 00036 /* >>chng 01 jan 17, read filename first, it may contain digits that would upset FFmtRead, 00037 * as well as other tests that will follow. GetQuote will erase the filename from chCard, PvH */ 00038 lgQuoteFound = strchr( input.chOrgCard, '\"' ) != NULL; 00039 if( lgQuoteFound ) 00040 { 00041 /* this will both scan in whatever label is inside the quotes in OrgCard, 00042 * but also remove the contents there and in chCard, 00043 * so that following keywords will not trigger off it */ 00044 GetQuote( chFile, chCard, true ); 00045 } 00046 00047 if( nMatch("GREY",chCard) || nMatch("GRAY",chCard) || nMatch("grey_",chFile) ) 00048 gp.lgGreyGrain = true; 00049 else 00050 gp.lgGreyGrain = false; 00051 00052 /* now check for the keyword "function" - this is a special case, 00053 * where grain abundance varies with depth 00054 * if set true then will call routine GrnVryDpth */ 00055 if( nMatch("FUNC",chCard) ) 00056 gp.lgAbunVsDepth = true; 00057 else 00058 gp.lgAbunVsDepth = false; 00059 00060 /* check for the keyword "single bin" - 00061 * use resolved grain size distributions if not present 00062 * NB - logic is backwards here */ 00063 if( !nMatch("SING",chCard) ) 00064 lgSizeDistribution = true; 00065 else 00066 lgSizeDistribution = false; 00067 00068 /* check for the keyword "qheating" - 00069 * quantum heating should be turned on */ 00070 gp.lgForbidQHeating = false; 00071 00072 if( nMatch("QHEA",chCard) ) 00073 gp.lgRequestQHeating = true; 00074 else 00075 gp.lgRequestQHeating = false; 00076 00077 /* the keyword "no qheat" always takes precedence */ 00078 if( nMatch("O QH",chCard) ) 00079 { 00080 gp.lgForbidQHeating = true; 00081 gp.lgRequestQHeating = false; 00082 phycon.lgPhysOK = false; 00083 } 00084 00085 /* option to force constant reevaluation of grain physics - 00086 * usually reevaluate grains at all times, but NO REEVALUATE will 00087 * save some time but may affect stability */ 00088 gv.lgReevaluate = !nMatch(" NO REEV",chCard); 00089 00090 /* option to turn off photoelectric heating by grain, NO HEATING */ 00091 if( nMatch("O HE",chCard) ) 00092 { 00093 phycon.lgPhysOK = false; 00094 gv.lgDHetOn = false; 00095 } 00096 00097 /* option to turn off gas cooling by grain, NO COOLING */ 00098 if( nMatch("O CO",chCard) ) 00099 { 00100 phycon.lgPhysOK = false; 00101 gv.lgDColOn = false; 00102 } 00103 00104 /* these are keywords for PAH's, they need to be read before the depletion factor */ 00105 if( (ptr = strstr(chCard,"C120 ")) != NULL ) 00106 { 00107 lgC120 = true; 00108 /* erase this keyword, it upsets FFmtRead */ 00109 strncpy(ptr," ",5); 00110 } 00111 else if( (ptr = strstr(chCard,"C15 ")) != NULL ) 00112 { 00113 lgC15 = true; 00114 /* erase this keyword, it upsets FFmtRead */ 00115 strncpy(ptr," ",4); 00116 } 00117 00118 /* log - linear option for grain abundance */ 00119 lgLogLinSet = false; 00120 lgLinSet = false; 00121 if( nMatch(" LOG",chCard) ) 00122 { 00123 lgLogLinSet = true; 00124 lgLinSet = false; 00125 } 00126 else if( nMatch("LINE",chCard) ) 00127 { 00128 lgLogLinSet = true; 00129 lgLinSet = true; 00130 } 00131 00132 /* get the grain abundance as the first parameter, 00133 * returns 0 if no number, ok since interpreted as log, or unity*/ 00134 i = 5; 00135 gp.dep = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00136 00137 /* was keyword log or linear on the line? */ 00138 if( lgLogLinSet ) 00139 { 00140 /* log or linear was specified, which was it */ 00141 if( lgLinSet ) 00142 { 00143 /* linear quantity entered, check it */ 00144 if( gp.dep <= 0. ) 00145 { 00146 fprintf( ioQQQ, " Impossible value for linear abundance.\n" ); 00147 fprintf( ioQQQ, " Abundance entered was%10.2e\n", gp.dep ); 00148 fprintf( ioQQQ, " Sorry.\n" ); 00149 cdEXIT(EXIT_FAILURE); 00150 } 00151 } 00152 else 00153 { 00154 gp.dep = pow(10.,gp.dep); 00155 } 00156 } 00157 else 00158 { 00159 /* neither log nor linear specified, check sign 00160 * force it to be a log - linear if greater than 0 */ 00161 if( gp.dep <= 0. ) 00162 { 00163 gp.dep = pow(10.,gp.dep); 00164 } 00165 } 00166 00167 if( gp.dep < FLT_MIN ) 00168 { 00169 fprintf( ioQQQ, " Grain abundance entered here (%f) is impossible.\n", gp.dep ); 00170 cdEXIT(EXIT_FAILURE); 00171 } 00172 00173 /* it is possible that there is a file name on the command line - 00174 * if so then we want to call correct reoutine, and not look for keywords 00175 * the signature of a keyword is a pair of quotes - is one present? */ 00176 if( lgQuoteFound ) 00177 { 00178 /* read the file name that was specified */ 00179 chOption = ""; 00180 mie_read_opc(chFile,gp); 00181 } 00182 else 00183 { 00184 if( nMatch("ORIO",chCard) ) 00185 { 00186 /* This scales the Orion grain abundance so that the observed 00187 * dust to gas ratio that Cloudy predicts is in agreement with 00188 * that observed in the Veil, 00189 *>>refer grain Abel, N., Brogan, C., Ferland, G., O'Dell, C.R., 00190 *>>refercon Shaw, G., Troland, T., 2004, ApJ, submitted */ 00191 gp.dep *= 0.85; 00192 00193 /* optional keyword ORION to use orion curves for large R grains */ 00194 /* only turn on one if graphite or silicate is on line, both if not */ 00195 if( nMatch("GRAP",chCard) ) 00196 { 00197 /* only turn on orion graphite */ 00198 chOption = "ORION GRAPHITE "; 00199 if( lgSizeDistribution ) 00200 { 00201 mie_read_opc("graphite_orion_10.opc",gp); 00202 } 00203 else 00204 { 00205 mie_read_opc("graphite_orion_01.opc",gp); 00206 } 00207 } 00208 else if( nMatch("SILI",chCard) ) 00209 { 00210 /* only turn on orion silicate */ 00211 chOption = "ORION SILICATE "; 00212 if( lgSizeDistribution ) 00213 { 00214 mie_read_opc("silicate_orion_10.opc",gp); 00215 } 00216 else 00217 { 00218 mie_read_opc("silicate_orion_01.opc",gp); 00219 } 00220 } 00221 else 00222 { 00223 /* turn both on */ 00224 chOption = "ORION "; 00225 if( lgSizeDistribution ) 00226 { 00227 mie_read_opc("graphite_orion_10.opc",gp); 00228 mie_read_opc("silicate_orion_10.opc",gp); 00229 } 00230 else 00231 { 00232 mie_read_opc("graphite_orion_01.opc",gp); 00233 mie_read_opc("silicate_orion_01.opc",gp); 00234 } 00235 } 00236 } 00237 00238 else if( nMatch(" PAH",chCard) ) 00239 { 00240 /* only turn on the large PAH */ 00241 if( lgC120 ) 00242 { 00243 chOption = "PAH C120 "; 00244 mie_read_opc("pah1_0n682.opc",gp); 00245 } 00246 /* only turn on the small PAH */ 00247 else if( lgC15 ) 00248 { 00249 chOption = "PAH C15 "; 00250 mie_read_opc("pah1_0n341.opc",gp); 00251 } 00252 /* turn on size-distributed PAHs */ 00253 else 00254 { 00255 chOption = "PAH "; 00256 if( lgSizeDistribution ) 00257 { 00258 mie_read_opc("pah1_ab08_10.opc",gp); 00259 } 00260 else 00261 { 00262 mie_read_opc("pah1_ab08_01.opc",gp); 00263 } 00264 } 00265 } 00266 00267 else if( nMatch("GREY",chCard) || nMatch("GRAY",chCard) ) 00268 { 00269 /* grey grains */ 00270 chOption = "GREY "; 00271 if( lgSizeDistribution ) 00272 { 00273 mie_read_opc("grey_ism_10.opc",gp); 00274 } 00275 else 00276 { 00277 mie_read_opc("grey_ism_01.opc",gp); 00278 } 00279 } 00280 00281 else if( nMatch(" ISM",chCard) ) 00282 { 00283 if( nMatch("GRAP",chCard) ) 00284 { 00285 /* only turn on ism graphite */ 00286 chOption = "ISM GRAPHITE "; 00287 if( lgSizeDistribution ) 00288 { 00289 mie_read_opc("graphite_ism_10.opc",gp); 00290 } 00291 else 00292 { 00293 mie_read_opc("graphite_ism_01.opc",gp); 00294 } 00295 } 00296 else if( nMatch("SILI",chCard) ) 00297 { 00298 /* only turn on orion silicate */ 00299 chOption = "ISM SILICATE "; 00300 if( lgSizeDistribution ) 00301 { 00302 mie_read_opc("silicate_ism_10.opc",gp); 00303 } 00304 else 00305 { 00306 mie_read_opc("silicate_ism_01.opc",gp); 00307 } 00308 } 00309 else 00310 { 00311 /* turn both ISM graphite and silicate on */ 00312 chOption = "ISM "; 00313 if( lgSizeDistribution ) 00314 { 00315 mie_read_opc("graphite_ism_10.opc",gp); 00316 mie_read_opc("silicate_ism_10.opc",gp); 00317 } 00318 else 00319 { 00320 mie_read_opc("graphite_ism_01.opc",gp); 00321 mie_read_opc("silicate_ism_01.opc",gp); 00322 } 00323 } 00324 } 00325 00326 /* default case */ 00327 else 00328 { 00329 /* turn both ISM graphite and silicate on */ 00330 chOption = ""; 00331 if( lgSizeDistribution ) 00332 { 00333 mie_read_opc("graphite_ism_10.opc",gp); 00334 mie_read_opc("silicate_ism_10.opc",gp); 00335 } 00336 else 00337 { 00338 mie_read_opc("graphite_ism_01.opc",gp); 00339 mie_read_opc("silicate_ism_01.opc",gp); 00340 } 00341 } 00342 } 00343 00344 /* vary option */ 00345 if( optimize.lgVarOn ) 00346 { 00347 optimize.nvfpnt[optimize.nparm] = input.nRead; 00348 optimize.vparm[0][optimize.nparm] = (realnum)log10(gp.dep); 00349 optimize.vincr[optimize.nparm] = 1.; 00350 00351 // now build the input command string with all the necessary options... 00352 string command( "GRAIN ABUND=%f LOG " ); 00353 if( chFile[0] != '\0' ) 00354 { 00355 command += "\""; 00356 command += chFile; 00357 command += "\" "; 00358 } 00359 00360 // every branch of the if-statement above must set chOption. 00361 // this sentinel is here in case somebody forgets... 00362 if( chOption == NULL ) 00363 TotalInsanity(); 00364 00365 command += chOption; 00366 00367 if( gp.lgAbunVsDepth ) 00368 command += "FUNCTION "; 00369 if( !lgSizeDistribution ) 00370 command += "SINGLE "; 00371 if( gp.lgForbidQHeating ) 00372 command += "NO QHEAT "; 00373 else if( gp.lgRequestQHeating ) 00374 command += "QHEAT "; 00375 if( !gv.lgReevaluate ) 00376 command += "NO REEVALUATE "; 00377 if( !gv.lgDHetOn ) 00378 command += "NO HEATING "; 00379 if( !gv.lgDColOn ) 00380 command += "NO COOLING "; 00381 00382 if( command.length() < static_cast<string::size_type>(FILENAME_PATH_LENGTH_2) ) 00383 strcpy( optimize.chVarFmt[optimize.nparm], command.c_str() ); 00384 else 00385 { 00386 fprintf(ioQQQ," grain command string is too long. This is parse_grain\n"); 00387 TotalInsanity(); 00388 } 00389 00390 optimize.nvarxt[optimize.nparm] = 1; 00391 ++optimize.nparm; 00392 } 00393 return; 00394 }