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 /*ParseBlackbody parse parameters off black body command */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "rfield.h" 00007 #include "radius.h" 00008 #include "parse.h" 00009 00010 void ParseBlackbody(char *chCard, /* input command line, already changed to caps */ 00011 /* counter for which continuum source this is */ 00012 long int *nqh, 00013 /* optional area that might be set here */ 00014 realnum *ar1) 00015 { 00016 bool lgEOL; 00017 long int i; 00018 double a, 00019 dil, 00020 rlogl; 00021 00022 DEBUG_ENTRY( "ParseBlackbody()" ); 00023 00024 /* type is blackbody */ 00025 strcpy( rfield.chSpType[rfield.nspec], "BLACK" ); 00026 00027 /* these two are not used for this continuum shape */ 00028 rfield.cutoff[rfield.nspec][0] = 0.; 00029 rfield.cutoff[rfield.nspec][1] = 0.; 00030 00031 /* get the temperature */ 00032 i = 5; 00033 rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00034 00035 /* there must be at least one number, the temperature */ 00036 if( lgEOL ) 00037 { 00038 fprintf( ioQQQ, " There must be at least 1 number on a BLACKBODY command line. Sorry.\n" ); 00039 cdEXIT(EXIT_FAILURE); 00040 } 00041 00042 /* this is the temperature - make sure its linear in the end 00043 * there are two keys, LINEAR and LOG, that could be here, 00044 * else choose which is here by which side of 10 */ 00045 if( (rfield.slope[rfield.nspec] <= 10. && !nMatch("LINE",chCard)) || 00046 nMatch(" LOG",chCard) ) 00047 { 00048 /* log option */ 00049 rfield.slope[rfield.nspec] = pow(10.,rfield.slope[rfield.nspec]); 00050 } 00051 00052 /* check that temp is not too low - could happen if log misused */ 00053 if( rfield.slope[rfield.nspec] < 1e4 ) 00054 { 00055 fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n", 00056 rfield.slope[rfield.nspec] ); 00057 } 00058 00059 /* now check that temp not too low - would peak below low 00060 * energy limit of the code 00061 * factor is temperature of 1 Ryd, egamry is high-energy limit of code */ 00062 if( rfield.slope[rfield.nspec]/TE1RYD < rfield.emm ) 00063 { 00064 fprintf( ioQQQ, " This temperature is very low - the blackbody will have significant flux low the low energy limit of the code, presently %10.2e Ryd.\n", 00065 rfield.emm ); 00066 fprintf( ioQQQ, " Was this intended?\n" ); 00067 } 00068 00069 /* now check that temp not too high - would extend beyond high 00070 * energy limit of the code 00071 * factor is temperature of 1 Ryd, egamry is high-energy limit of code */ 00072 if( rfield.slope[rfield.nspec]/TE1RYD*2. > rfield.egamry ) 00073 { 00074 fprintf( ioQQQ, " This temperature is very high - the blackbody will have significant flux above the high-energy limit of the code,%10.2e Ryd.\n", 00075 rfield.egamry ); 00076 fprintf( ioQQQ, " Was this intended?\n" ); 00077 } 00078 00079 /* increment continuum indices */ 00080 ++rfield.nspec; 00081 if( rfield.nspec >= LIMSPC ) 00082 { 00083 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00084 cdEXIT(EXIT_FAILURE); 00085 } 00086 00087 /* also possible to input log(total luminosity)=real log(l) */ 00088 a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00089 if( lgEOL ) 00090 { 00091 /* there was not a second number on the line; check if LTE or STE */ 00092 if( nMatch(" LTE",chCard) || nMatch("LTE ",chCard) || 00093 nMatch(" STE",chCard) || nMatch("STE ",chCard) ) 00094 { 00095 /* set energy density to the STE - strict thermodynamic equilibrium - value */ 00096 strcpy( rfield.chSpNorm[*nqh], "LUMI" ); 00097 00098 /* need nspec-1 since was incremented above */ 00099 a = log10(rfield.slope[rfield.nspec-1]); 00100 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a; 00101 00102 /* set radius to very large value if not already set */ 00103 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00104 if( !radius.lgRadiusKnown ) 00105 { 00106 *ar1 = (realnum)radius.rdfalt; 00107 radius.Radius = pow(10.,radius.rdfalt); 00108 } 00109 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00110 rfield.range[*nqh][0] = rfield.emm; 00111 rfield.range[*nqh][1] = rfield.egamry; 00112 rfield.totpow[*nqh] = rlogl; 00113 *nqh = MIN2(*nqh+1,LIMSPC); 00114 00115 /* check that stack of shape and luminosity specifications 00116 * is parallel, stop if not - this happens is background comes 00117 * BETWEEN another set of shape and luminosity commands */ 00118 if( rfield.nspec != *nqh ) 00119 { 00120 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00121 fprintf( ioQQQ, " Sorry.\n" ); 00122 cdEXIT(EXIT_FAILURE); 00123 } 00124 } 00125 00126 return; 00127 } 00128 00129 strcpy( rfield.chSpNorm[*nqh], "LUMI" ); 00130 00131 /* a second number was entered, what was it? */ 00132 if( nMatch("LUMI",chCard) ) 00133 { 00134 rlogl = a; 00135 strcpy( rfield.chRSpec[*nqh], "4 PI" ); 00136 } 00137 00138 else if( nMatch("RADI",chCard) ) 00139 { 00140 /* radius was entered, convert to total lumin */ 00141 rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nspec-1]); 00142 strcpy( rfield.chRSpec[*nqh], "4 PI" ); 00143 } 00144 00145 else if( nMatch("DENS",chCard) ) 00146 { 00147 /* number was temperature to deduce energy density 00148 * number is linear if greater than 10, or if LINEAR appears on line 00149 * want number to be log of temperature at end of this */ 00150 if( nMatch("LINE",chCard) || a > 10. ) 00151 { 00152 a = log10(a); 00153 } 00154 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a; 00155 /* set radius to very large value if not already set */ 00156 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00157 if( !radius.lgRadiusKnown ) 00158 { 00159 *ar1 = (realnum)radius.rdfalt; 00160 radius.Radius = pow(10.,radius.rdfalt); 00161 } 00162 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00163 } 00164 00165 else if( nMatch("DILU",chCard) ) 00166 { 00167 /* number is dilution factor, if negative then its log */ 00168 if( a < 0. ) 00169 { 00170 dil = a; 00171 } 00172 else 00173 { 00174 dil = log10(a); 00175 } 00176 if( dil > 0. ) 00177 { 00178 fprintf( ioQQQ, " Is a dilution factor greater than one physical?\n" ); 00179 } 00180 00181 /* this is LTE energy density */ 00182 a = log10(rfield.slope[rfield.nspec-1]); 00183 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a; 00184 00185 /* add on dilution factor */ 00186 rlogl += dil; 00187 00188 /* set radius to very large value if not already set */ 00189 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00190 if( !radius.lgRadiusKnown ) 00191 { 00192 *ar1 = (realnum)radius.rdfalt; 00193 radius.Radius = pow(10.,radius.rdfalt); 00194 } 00195 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00196 00197 } 00198 else 00199 { 00200 fprintf( ioQQQ, " I do not understand what the second number was- keywords are LUMINOSITY, RADIUS, DILUTION, and DENSITY.\n" ); 00201 cdEXIT(EXIT_FAILURE); 00202 } 00203 00204 rfield.range[*nqh][0] = rfield.emm; 00205 rfield.range[*nqh][1] = rfield.egamry; 00206 rfield.totpow[*nqh] = rlogl; 00207 /* the second number will be some sort of luminosity */ 00208 *nqh = MIN2(*nqh+1,LIMSPC); 00209 00210 /* check that stack of shape and luminosity specifications 00211 * is parallel, stop if not - this happens is background comes 00212 * BETWEEN another set of shape and luminosity commands */ 00213 if( rfield.nspec != *nqh ) 00214 { 00215 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00216 fprintf( ioQQQ, " Sorry.\n" ); 00217 cdEXIT(EXIT_FAILURE); 00218 } 00219 00220 return; 00221 }