00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "optimize.h"
00007 #include "input.h"
00008 #include "rfield.h"
00009 #include "radius.h"
00010 #include "parse.h"
00011 #include "parser.h"
00012
00013 void ParseBlackbody(
00014
00015 Parser &p)
00016 {
00017 bool
00018 lgIntensitySet=false;
00019 double a,
00020 dil,
00021 rlogl;
00022
00023 char chParamType[20] = "";
00024 int nParam = 0;
00025
00026 DEBUG_ENTRY( "ParseBlackbody()" );
00027
00028 set_NaN( rlogl );
00029
00030
00031 strcpy( rfield.chSpType[rfield.nShape], "BLACK" );
00032 strcpy( rfield.chSpNorm[p.m_nqh], "LUMI" );
00033
00034
00035 rfield.cutoff[rfield.nShape][0] = 0.;
00036 rfield.cutoff[rfield.nShape][1] = 0.;
00037
00038
00039 rfield.slope[rfield.nShape] = p.FFmtRead();
00040 if( p.lgEOL() )
00041 p.NoNumb("blackbody temperature");
00042
00043
00044
00045
00046 if( (rfield.slope[rfield.nShape] <= 10. && !p.nMatch("LINE")) ||
00047 p.nMatch(" LOG") )
00048 {
00049
00050 rfield.slope[rfield.nShape] = pow(10.,rfield.slope[rfield.nShape]);
00051 }
00052
00053
00054 if( rfield.slope[rfield.nShape] < 1e4 )
00055 {
00056 fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n",
00057 rfield.slope[rfield.nShape] );
00058 }
00059
00060
00061
00062
00063 if( rfield.slope[rfield.nShape]/TE1RYD < rfield.emm )
00064 {
00065 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",
00066 rfield.emm );
00067 fprintf( ioQQQ, " Was this intended?\n" );
00068 }
00069
00070
00071
00072
00073 if( rfield.slope[rfield.nShape]/TE1RYD*2. > rfield.egamry )
00074 {
00075 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",
00076 rfield.egamry );
00077 fprintf( ioQQQ, " Was this intended?\n" );
00078 }
00079
00080
00081 a = p.FFmtRead();
00082
00083
00084 if( p.nMatch(" LTE") || p.nMatch("LTE ") ||
00085 p.nMatch(" STE") || p.nMatch("STE ") )
00086 {
00087
00088 strcpy( chParamType , "STE" );
00089 nParam = 1;
00090
00091 if( !p.lgEOL() )
00092 {
00093 fprintf(ioQQQ,"PROBLEM the luminosity was specified on "
00094 "the BLACKBODY K STE command.\n");
00095 fprintf(ioQQQ,"Do not specify the luminosity since STE does this.\n");
00096 fprintf( ioQQQ, " Sorry.\n" );
00097 cdEXIT(EXIT_FAILURE);
00098 }
00099
00100
00101 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);
00102
00103
00104 if( !radius.lgRadiusKnown )
00105 {
00106 radius.Radius = pow(10.,radius.rdfalt);
00107 }
00108
00109 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00110 lgIntensitySet = true;
00111 }
00112
00113
00114 else if( p.nMatch("LUMI") )
00115 {
00116 strcpy( chParamType , "LUMINOSITY" );
00117 nParam = 2;
00118 rlogl = a;
00119 strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
00120 if( p.lgEOL() )
00121 p.NoNumb("luminosity" );
00122 lgIntensitySet = true;
00123 }
00124
00125 else if( p.nMatch("RADI") )
00126 {
00127 strcpy( chParamType , "RADIUS" );
00128 nParam = 2;
00129
00130 rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nShape]);
00131 strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
00132 if( p.lgEOL() )
00133 p.NoNumb("radius" );
00134 lgIntensitySet = true;
00135 }
00136
00137 else if( p.nMatch("DENS") )
00138 {
00139 strcpy( chParamType , "ENERGY DENSITY" );
00140 nParam = 2;
00141
00142
00143
00144 if( !p.nMatch(" LOG") && (p.nMatch("LINE") || a > 10.) )
00145 {
00146 a = log10(a);
00147 }
00148 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
00149
00150 if( !radius.lgRadiusKnown )
00151 {
00152 radius.Radius = pow(10.,radius.rdfalt);
00153 }
00154 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00155 if( p.lgEOL() )
00156 p.NoNumb("energy density");
00157 lgIntensitySet = true;
00158 }
00159
00160 else if( p.nMatch("DILU") )
00161 {
00162 strcpy( chParamType , "DILUTION FACTOR" );
00163 nParam = 2;
00164
00165 if( a <= 0. )
00166 dil = a;
00167 else
00168 dil = log10(a);
00169
00170 if( dil > 0. )
00171 fprintf( ioQQQ, "PROBLEM Is the dilution factor > 1 on this "
00172 "blackbody command physical?\n" );
00173
00174
00175 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);
00176
00177
00178 rlogl += dil;
00179
00180
00181 if( !radius.lgRadiusKnown )
00182 {
00183 radius.Radius = pow(10.,radius.rdfalt);
00184 }
00185 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00186 if( p.lgEOL() )
00187 p.NoNumb("dilution factor" );
00188 lgIntensitySet = true;
00189 }
00190
00191 else if( p.nMatch("DISK") )
00192 {
00193 if( p.lgEOL() )
00194 p.NoNumb("disk Te" );
00195
00196 strcpy( chParamType , "DISK" );
00197 nParam = 2;
00198
00199 rfield.cutoff[rfield.nShape][0] = a;
00200
00201
00202
00203 if( (rfield.cutoff[rfield.nShape][0] <= 10. && !p.nMatch("LINE")) ||
00204 p.nMatch(" LOG") )
00205 {
00206
00207 rfield.cutoff[rfield.nShape][0] = pow(10.,rfield.cutoff[rfield.nShape][0]);
00208 }
00209 a = log10( rfield.cutoff[rfield.nShape][0] );
00210
00211 strcpy( rfield.chSpType[rfield.nShape], "DISKB" );
00212 lgIntensitySet = false;
00213 }
00214
00215 if( lgIntensitySet )
00216 {
00217
00218
00219
00220
00221 if( rfield.nShape != p.m_nqh )
00222 {
00223 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" );
00224 fprintf( ioQQQ, " Sorry.\n" );
00225 cdEXIT(EXIT_FAILURE);
00226 }
00227
00228 rfield.range[p.m_nqh][0] = rfield.emm;
00229 rfield.range[p.m_nqh][1] = rfield.egamry;
00230 rfield.totpow[p.m_nqh] = rlogl;
00231 ++p.m_nqh;
00232 }
00233
00234 if( optimize.lgVarOn )
00235 {
00236
00237 if( strcmp(chParamType,"") == 0 )
00238 {
00239
00240 optimize.nvarxt[optimize.nparm] = 1;
00241 strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody= %f LOG" );
00242 }
00243 else
00244 {
00245 char chHold[100];
00246
00247 if( nParam==1 )
00248 {
00249 optimize.nvarxt[optimize.nparm] = 1;
00250 strcpy( chHold , "BLACKbody= %f LOG ");
00251 strcat( chHold , chParamType );
00252 }
00253 else if( nParam==2 )
00254 {
00255 optimize.nvarxt[optimize.nparm] = 2;
00256 optimize.vparm[1][optimize.nparm] = (realnum)a;
00257 strcpy( chHold , "BLACKbody= %f LOG %f ");
00258 strcat( chHold , chParamType );
00259 }
00260 else
00261 TotalInsanity();
00262 strcpy( optimize.chVarFmt[optimize.nparm], chHold );
00263 }
00264
00265
00266 optimize.nvfpnt[optimize.nparm] = input.nRead;
00267
00268 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nShape]);
00269
00270 optimize.vincr[optimize.nparm] = 0.5f;
00271 ++optimize.nparm;
00272 }
00273
00274
00275 ++rfield.nShape;
00276 if( rfield.nShape >= LIMSPC )
00277 {
00278 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
00279 cdEXIT(EXIT_FAILURE);
00280 }
00281
00282 return;
00283 }