cloudy trunk
Loading...
Searching...
No Matches
parse_blackbody.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*ParseBlackbody parse parameters off black body command */
4#include "cddefines.h"
5#include "physconst.h"
6#include "optimize.h"
7#include "input.h"
8#include "rfield.h"
9#include "radius.h"
10#include "parse.h"
11#include "parser.h"
12
14 /* input command line, already changed to caps */
15 Parser &p)
16{
17 bool
18 lgIntensitySet=false;
19 double a,
20 dil,
21 rlogl;
22
23 char chParamType[20] = "";
24 int nParam = 0;
25
26 DEBUG_ENTRY( "ParseBlackbody()" );
27
28 set_NaN( rlogl );
29
30 /* type is blackbody */
31 strcpy( rfield.chSpType[rfield.nShape], "BLACK" );
32 strcpy( rfield.chSpNorm[p.m_nqh], "LUMI" );
33
34 /* these two are not used for this continuum shape */
35 rfield.cutoff[rfield.nShape][0] = 0.;
36 rfield.cutoff[rfield.nShape][1] = 0.;
37
38 /* get the blackbody temperature */
39 rfield.slope[rfield.nShape] = p.FFmtRead();
40 if( p.lgEOL() )
41 p.NoNumb("blackbody temperature");
42
43 /* this is the temperature - make sure its linear in the end
44 * there are two keys, LINEAR and LOG, that could be here,
45 * else choose which is here by which side of 10 */
46 if( (rfield.slope[rfield.nShape] <= 10. && !p.nMatch("LINE")) ||
47 p.nMatch(" LOG") )
48 {
49 /* log option */
50 if( rfield.slope[rfield.nShape]>log10(BIGFLOAT) )
51 {
52 fprintf(ioQQQ,"PROBLEM The specified log of the temperature, %.3e, is too large.\nSorry.\n",
53 rfield.slope[rfield.nShape]);
55 }
56 rfield.slope[rfield.nShape] = pow(10.,rfield.slope[rfield.nShape]);
57 }
58
59 /* check that temp is not too low - could happen if log misused */
60 if( rfield.slope[rfield.nShape] < 1e4 )
61 {
62 fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n",
63 rfield.slope[rfield.nShape] );
64 }
65
66 /* now check that temp not too low - would peak below low
67 * energy limit of the code
68 * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
69 if( rfield.slope[rfield.nShape]/TE1RYD < rfield.emm )
70 {
71 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",
72 rfield.emm );
73 fprintf( ioQQQ, " Was this intended?\n" );
74 }
75
76 /* now check that temp not too high - would extend beyond high
77 * energy limit of the code
78 * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
79 if( rfield.slope[rfield.nShape]/TE1RYD*2. > rfield.egamry )
80 {
81 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",
82 rfield.egamry );
83 fprintf( ioQQQ, " Was this intended?\n" );
84 }
85
86 /* also possible to input log(total luminosity)=real log(l) */
87 a = p.FFmtRead();
88
89 /* there was not a second number on the line; check if LTE or STE */
90 if( p.nMatch(" LTE") || p.nMatch("LTE ") ||
91 p.nMatch(" STE") || p.nMatch("STE ") )
92 {
93 /* set energy density to the STE - strict thermodynamic equilibrium - value */
94 strcpy( chParamType , "STE" );
95 nParam = 1;
96
97 if( !p.lgEOL() )
98 {
99 fprintf(ioQQQ,"PROBLEM the luminosity was specified on "
100 "the BLACKBODY K STE command.\n");
101 fprintf(ioQQQ,"Do not specify the luminosity since STE does this.\n");
102 fprintf( ioQQQ, " Sorry.\n" );
104 }
105
106 /* use blackbody relations to get intensity from temperature */
107 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);
108
109 /* set radius to very large value if not already set */
110 if( !radius.lgRadiusKnown )
111 {
112 radius.Radius = pow(10.,radius.rdfalt);
113 }
114
115 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
116 lgIntensitySet = true;
117 }
118
119 /* a second number was entered, what was it? */
120 else if( p.nMatch("LUMI") )
121 {
122 strcpy( chParamType , "LUMINOSITY" );
123 nParam = 2;
124 rlogl = a;
125 strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
126 if( p.lgEOL() )
127 p.NoNumb("luminosity" );
128 lgIntensitySet = true;
129 }
130
131 else if( p.nMatch("RADI") )
132 {
133 strcpy( chParamType , "RADIUS" );
134 nParam = 2;
135 /* radius was entered, convert to total luminosity */
136 rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nShape]);
137 strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
138 if( p.lgEOL() )
139 p.NoNumb("radius" );
140 lgIntensitySet = true;
141 }
142
143 else if( p.nMatch("DENS") )
144 {
145 strcpy( chParamType , "ENERGY DENSITY" );
146 nParam = 2;
147 /* number was temperature to deduce energy density
148 * number is linear if greater than 10, or if LINEAR appears on line
149 * want number to be log of temperature at end of this */
150 if( !p.nMatch(" LOG") && (p.nMatch("LINE") || a > 10.) )
151 {
152 a = log10(a);
153 }
154 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
155 /* set radius to very large value if not already set */
156 if( !radius.lgRadiusKnown )
157 {
158 radius.Radius = pow(10.,radius.rdfalt);
159 }
160 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
161 if( p.lgEOL() )
162 p.NoNumb("energy density");
163 lgIntensitySet = true;
164 }
165
166 else if( p.nMatch("DILU") )
167 {
168 strcpy( chParamType , "DILUTION FACTOR" );
169 nParam = 2;
170 /* number is dilution factor, if negative then its log */
171 if( a <= 0. )
172 dil = a;
173 else
174 dil = log10(a);
175
176 if( dil > 0. )
177 fprintf( ioQQQ, "PROBLEM Is the dilution factor > 1 on this "
178 "blackbody command physical?\n" );
179
180 /* intensity from black body relations and temperature */
181 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);
182
183 /* add on dilution factor */
184 rlogl += dil;
185
186 /* set radius to very large value if not already set */
187 if( !radius.lgRadiusKnown )
188 {
189 radius.Radius = pow(10.,radius.rdfalt);
190 }
191 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
192 if( p.lgEOL() )
193 p.NoNumb("dilution factor" );
194 lgIntensitySet = true;
195 }
196
197 else if( p.nMatch("DISK") )
198 {
199 if( p.lgEOL() )
200 p.NoNumb("disk Te" );
201
202 strcpy( chParamType , "DISK" );
203 nParam = 2;
204
205 rfield.cutoff[rfield.nShape][0] = a;
206 /* this is the temperature - make sure its linear in the end
207 * there are two keys, LINEAR and LOG, that could be here,
208 * else choose which is here by which side of 10 */
209 if( (rfield.cutoff[rfield.nShape][0] <= 10. && !p.nMatch("LINE")) ||
210 p.nMatch(" LOG") )
211 {
212 /* log option */
213 rfield.cutoff[rfield.nShape][0] = pow(10.,rfield.cutoff[rfield.nShape][0]);
214 }
215 a = log10( rfield.cutoff[rfield.nShape][0] );
216
217 strcpy( rfield.chSpType[rfield.nShape], "DISKB" );
218 lgIntensitySet = false;
219 }
220
221 if( lgIntensitySet )
222 {
223 /* a luminosity option was specified
224 * check that stack of shape and luminosity specifications
225 * is parallel, stop if not - this happens is background comes
226 * BETWEEN another set of shape and luminosity commands */
227 if( rfield.nShape != p.m_nqh )
228 {
229 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" );
230 fprintf( ioQQQ, " Sorry.\n" );
232 }
233
234 rfield.range[p.m_nqh][0] = rfield.emm;
235 rfield.range[p.m_nqh][1] = rfield.egamry;
236 rfield.totpow[p.m_nqh] = rlogl;
237 ++p.m_nqh;
238 }
239 /* vary option */
240 if( optimize.lgVarOn )
241 {
242 /* this test no option on blackbody command */
243 if( strcmp(chParamType,"") == 0 )
244 {
245 /* no luminosity options on vary */
246 optimize.nvarxt[optimize.nparm] = 1;
247 strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody= %f LOG" );
248 }
249 else
250 {
251 char chHold[100];
252 /* there was an option - honor it */
253 if( nParam==1 )
254 {
255 optimize.nvarxt[optimize.nparm] = 1;
256 strcpy( chHold , "BLACKbody= %f LOG ");
257 strcat( chHold , chParamType );
258 }
259 else if( nParam==2 )
260 {
261 optimize.nvarxt[optimize.nparm] = 2;
262 optimize.vparm[1][optimize.nparm] = (realnum)a;
263 strcpy( chHold , "BLACKbody= %f LOG %f ");
264 strcat( chHold , chParamType );
265 }
266 else
268 strcpy( optimize.chVarFmt[optimize.nparm], chHold );
269 }
270
271 /* pointer to where to write */
272 optimize.nvfpnt[optimize.nparm] = input.nRead;
273 /* log of temp stored here */
274 optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nShape]);
275 /* the increment in the first steps away from the original value */
276 optimize.vincr[optimize.nparm] = 0.5f;
277 ++optimize.nparm;
278 }
279
280 /* increment SED indices */
281 ++rfield.nShape;
282 if( rfield.nShape >= LIMSPC )
283 {
284 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
286 }
287
288 return;
289}
FILE * ioQQQ
Definition cddefines.cpp:7
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
bool lgEOL(void) const
Definition parser.h:98
NORETURN void NoNumb(const char *chDesc) const
Definition parser.cpp:233
long int m_nqh
Definition parser.h:41
void set_NaN(sys_float &x)
Definition cpu.cpp:682
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
t_input input
Definition input.cpp:12
t_optimize optimize
Definition optimize.cpp:5
void ParseBlackbody(Parser &p)
UNUSED const double STEFAN_BOLTZ
Definition physconst.h:210
UNUSED const double TE1RYD
Definition physconst.h:183
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
const int LIMSPC
Definition rfield.h:18