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 /*ParseAbundances parse and read in composition as set by abundances command */ 00004 #include "cddefines.h" 00005 #include "grains.h" 00006 #include "abund.h" 00007 #include "phycon.h" 00008 #include "called.h" 00009 #include "elementnames.h" 00010 #include "input.h" 00011 #include "parse.h" 00012 00013 void ParseAbundances(char *chCard, 00014 /* following set true by grains command, 00015 * so this will not set more grains if grains already on. */ 00016 bool lgDSet) 00017 { 00018 bool lgEOF, 00019 lgEOL, 00020 lgLog; 00021 long int i, 00022 j, 00023 k; 00024 double absav[LIMELM], 00025 chk; 00026 GrainPar gp; 00027 00028 DEBUG_ENTRY( "ParseAbundances()" ); 00029 00030 j = 5; 00031 absav[0] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL); 00032 /* check whether it scanned off end of line on first try 00033 * at least one number on the line if it passes this test 00034 * no numbers means a keyword */ 00035 if( !lgEOL ) 00036 { 00037 absav[1] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL); 00038 if( lgEOL ) 00039 { 00040 /* this branch, we found one, but not two, numbers - 00041 * must be one of a few special cases */ 00042 if( nMatch(" ALL",chCard) ) 00043 { 00044 /* special option, all abundances will be this number */ 00045 if( absav[0] <= 0. ) 00046 { 00047 absav[0] = pow(10.,absav[0]); 00048 } 00049 for( i=1; i < LIMELM; i++ ) 00050 { 00051 abund.solar[i] = (realnum)absav[0]; 00052 } 00053 00054 } 00055 else if( nMatch("OLD ",chCard) && nMatch("SOLA",chCard) ) 00056 { 00057 i = (int)absav[0]; 00058 /* old solar - better be the number "84" on the line */ 00059 if( i!=84 ) 00060 { 00061 fprintf( ioQQQ, 00062 " The only old abundance set I have is for version 84 - %3ld was requested. Sorry.\n", 00063 i ); 00064 cdEXIT(EXIT_FAILURE); 00065 } 00066 for( i=1; i < LIMELM; i++ ) 00067 { 00068 /* these are the old abundances */ 00069 abund.SolarSave[i] = abund.OldSolar84[i]; 00070 abund.solar[i] = abund.OldSolar84[i]; 00071 } 00072 } 00073 00074 else if( nMatch("STAR",chCard) ) 00075 { 00076 /* Fred Hamonn's star burst galaxy mixture */ 00077 abund_starburst(chCard); 00078 } 00079 else 00080 { 00081 fprintf( ioQQQ, 00082 " I did not recognize a sub-keyword - options are ALL and OLD SOLAR 84. Sorry.\n"); 00083 cdEXIT(EXIT_FAILURE); 00084 } 00085 00086 /* normal return */ 00087 return; 00088 } 00089 00090 /* we get here if there is a second number - read in all abundances */ 00091 for( i=2; i < abund.npSolar; i++ ) 00092 { 00093 absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL); 00094 if( lgEOL ) 00095 { 00096 /* read CONTINUE line if J scanned off end before reading all abund */ 00097 input_readarray(chCard,&lgEOF); 00098 if( lgEOF ) 00099 { 00100 fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld. Sorry.\n", 00101 abund.npSolar, i ); 00102 cdEXIT(EXIT_FAILURE); 00103 } 00104 00105 /* option to ignore all lines starting with #, *, or % */ 00106 while( lgInputComment(chCard) ) 00107 { 00108 if( lgEOF ) 00109 { 00110 fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld. Sorry.\n", 00111 abund.npSolar, i ); 00112 cdEXIT(EXIT_FAILURE); 00113 } 00114 /* read in another line */ 00115 input_readarray(chCard,&lgEOF); 00116 } 00117 00118 /* we have the line image, print it */ 00119 if( called.lgTalk ) 00120 { 00121 fprintf( ioQQQ, " * "); 00122 k=0; 00123 while( chCard[k]!='\0' ) 00124 { 00125 fprintf(ioQQQ,"%c",chCard[k]); 00126 ++k; 00127 } 00128 while( k<80 ) 00129 { 00130 fprintf(ioQQQ,"%c",' '); 00131 ++k; 00132 } 00133 fprintf( ioQQQ,"*\n"); 00134 } 00135 00136 /* now convert to caps */ 00137 caps(chCard); 00138 if( strncmp(chCard,"CONT",4) != 0 ) 00139 { 00140 fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld. Sorry.\n", 00141 abund.npSolar, i ); 00142 cdEXIT(EXIT_FAILURE); 00143 } 00144 else 00145 { 00146 j = 1; 00147 absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL); 00148 if( lgEOF ) 00149 { 00150 fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld. Sorry.\n", 00151 abund.npSolar, i); 00152 cdEXIT(EXIT_FAILURE); 00153 } 00154 } 00155 } 00156 } 00157 00158 /* fell through to here after reading in N abundances for N elements 00159 * check that there are no more abundances on the line - that would 00160 * be an error - a typo */ 00161 chk = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL); 00162 if( !lgEOL || (chk!=0.) ) 00163 { 00164 /* got another number, not lgEOL, so too many numbers */ 00165 fprintf( ioQQQ, " There were more than %ld abundances entered\n", 00166 abund.npSolar ); 00167 fprintf( ioQQQ, " Could there have been a typo somewhere?\n" ); 00168 } 00169 00170 /* are numbers scale factors, or log of abund rel to H?? */ 00171 lgLog = false; 00172 for( i=0; i < abund.npSolar; i++ ) 00173 { 00174 if( absav[i] < 0. ) 00175 lgLog = true; 00176 } 00177 00178 if( lgLog ) 00179 { 00180 /* entered as log of number rel to hydrogen */ 00181 for( i=0; i < abund.npSolar; i++ ) 00182 { 00183 abund.solar[abund.ipSolar[i]-1] = (realnum)pow(10.,absav[i]); 00184 } 00185 } 00186 else 00187 { 00188 /* scale factors relative to solar */ 00189 for( i=0; i < abund.npSolar; i++ ) 00190 { 00191 abund.solar[abund.ipSolar[i]-1] *= (realnum)absav[i]; 00192 } 00193 } 00194 00195 /* check whether the abundances are reasonable */ 00196 if( abund.solar[1] > 0.2 ) 00197 { 00198 fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 00199 abund.solar[1], elementnames.chElementSym[1] ); 00200 } 00201 00202 for( i=2; i < LIMELM; i++ ) 00203 { 00204 if( abund.solar[i] > 0.1 ) 00205 { 00206 fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 00207 abund.solar[i], elementnames.chElementSym[i] ); 00208 } 00209 } 00210 return; 00211 } 00212 00213 /* this branch if no numbers on the line - option to use one of several 00214 * special stored abundance sets */ 00215 if( nMatch(" AGB",chCard) || nMatch("AGB ",chCard) || nMatch("PLAN",chCard) ) 00216 { 00217 /* AGB/planetary nebula abundances */ 00218 /* only turn on grains if "no grains" is not present */ 00219 if( !nMatch("O GR",chCard) ) 00220 { 00221 /* turn on grains if not already done */ 00222 if( !lgDSet ) 00223 { 00224 /* return bins allocated by previous abundances ... commands */ 00225 ReturnGrainBins(); 00226 /* now allocate new grain bins */ 00227 gp.dep = 1.; 00228 gp.lgAbunVsDepth = false; 00229 gp.lgRequestQHeating = true; 00230 gp.lgForbidQHeating = false; 00231 gp.lgGreyGrain = false; 00232 00233 /* NO QHEAT option to turn off quantum heating for grains */ 00234 if( nMatch("O QH",chCard) ) 00235 { 00236 gp.lgForbidQHeating = true; 00237 gp.lgRequestQHeating = false; 00238 phycon.lgPhysOK = false; 00239 } 00240 00241 /* actually set up the grains */ 00242 mie_read_opc("graphite_ism_10.opc",gp); 00243 mie_read_opc("silicate_ism_10.opc",gp); 00244 } 00245 } 00246 00247 for( i=0; i < LIMELM; i++ ) 00248 { 00249 abund.solar[i] = abund.apn[i]; 00250 } 00251 } 00252 00253 else if( nMatch("HII ",chCard) || nMatch("H II",chCard) || nMatch("ORIO",chCard) ) 00254 { 00255 /* H II region abundances */ 00256 /* only turn on grains if "no grains" is not present */ 00257 if( !nMatch("O GR",chCard) ) 00258 { 00259 /* option to turn on grains */ 00260 if( !lgDSet ) 00261 { 00262 /* return bins allocated by previous abundances ... commands */ 00263 ReturnGrainBins(); 00264 /* now allocate new grain bins */ 00265 gp.dep = 1.; 00266 gp.lgAbunVsDepth = false; 00267 gp.lgRequestQHeating = true; 00268 gp.lgForbidQHeating = false; 00269 gp.lgGreyGrain = false; 00270 00271 /* NO QHEAT option to turn off quantum heating for grains */ 00272 if( nMatch("O QH",chCard) ) 00273 { 00274 gp.lgForbidQHeating = true; 00275 gp.lgRequestQHeating = false; 00276 phycon.lgPhysOK = false; 00277 } 00278 /* This scales the Orion grain abundance so that the observed 00279 * dust to gas ratio that Cloudy predicts is in agreement with 00280 * that observed in the Veil, 00281 *>>refer grain Abel, N., Brogan, C., Ferland, G., O'Dell, C.R., 00282 *>>refercon Shaw, G., Troland, T., 2004, ApJ, submitted */ 00283 gp.dep *= 0.85; 00284 00285 mie_read_opc("graphite_orion_10.opc",gp); 00286 mie_read_opc("silicate_orion_10.opc",gp); 00287 } 00288 } 00289 00290 for( i=0; i < LIMELM; i++ ) 00291 { 00292 abund.solar[i] = abund.ahii[i]; 00293 } 00294 } 00295 00296 else if( nMatch("ISM ",chCard) || nMatch(" ISM",chCard) ) 00297 { 00298 /* ISM abundances from Cowie and Songaila Ann Rev '86 */ 00299 /* only turn on grains if "no grains" is not present */ 00300 if( !nMatch("O GR",chCard) ) 00301 { 00302 if( !lgDSet ) 00303 { 00304 /* return bins allocated by previous abundances ... commands */ 00305 ReturnGrainBins(); 00306 /* now allocate new grain bins */ 00307 gp.dep = 1.; 00308 gp.lgAbunVsDepth = false; 00309 gp.lgRequestQHeating = true; 00310 gp.lgForbidQHeating = false; 00311 gp.lgGreyGrain = false; 00312 00313 /* NO QHEAT option to turn off quantum heating */ 00314 if( nMatch("O QH",chCard) ) 00315 { 00316 gp.lgForbidQHeating = true; 00317 gp.lgRequestQHeating = false; 00318 phycon.lgPhysOK = false; 00319 } 00320 00321 /* actually set up the grains */ 00322 mie_read_opc("graphite_ism_10.opc",gp); 00323 mie_read_opc("silicate_ism_10.opc",gp); 00324 } 00325 } 00326 00327 for( i=0; i < LIMELM; i++ ) 00328 { 00329 abund.solar[i] = abund.aism[i]; 00330 } 00331 } 00332 00333 else if( nMatch("NOVA",chCard) ) 00334 { 00335 /* Nova Cyg abundances */ 00336 for( i=0; i < LIMELM; i++ ) 00337 { 00338 abund.solar[i] = abund.anova[i]; 00339 } 00340 } 00341 00342 else if( nMatch("PRIM",chCard) ) 00343 { 00344 /* roughly primordial abundances: He/H=.072 */ 00345 for( i=0; i < LIMELM; i++ ) 00346 { 00347 abund.solar[i] = abund.aprim[i]; 00348 } 00349 00350 /* now turn off the heavy elements */ 00351 for( i=4; i < LIMELM; i++ ) 00352 { 00353 /* turn off heavy elements - do this way to make sure, 00354 * that H-like and He-like level limits are handled properly */ 00355 char chDUMMY[LIMELM]; 00356 sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] ); 00357 /* now convert to caps */ 00358 caps(chDUMMY); 00359 ParseElement( chDUMMY ); 00360 } 00361 } 00362 00363 else if( nMatch("CAME",chCard) ) 00364 { 00365 /* mix from Cameron 1982, "Essays on Nuclear Astrophysics" */ 00366 for( i=0; i < LIMELM; i++ ) 00367 { 00368 abund.solar[i] = abund.camern[i]; 00369 } 00370 } 00371 00372 else 00373 { 00374 fprintf( ioQQQ, 00375 " ABUND must have PLAN, H II, CAMERON, NOVA, ALL, STARBURST, OLD SOLAR 84 or PRIMORDIAL. Sorry.\n" ); 00376 cdEXIT(EXIT_FAILURE); 00377 } 00378 return; 00379 }