cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_blackbody.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 "rfield.h"
7 #include "radius.h"
8 #include "parse.h"
9 
10 void ParseBlackbody(char *chCard, /* input command line, already changed to caps */
11  /* counter for which continuum source this is */
12  long int *nqh,
13  /* optional area that might be set here */
14  realnum *ar1)
15 {
16  bool lgEOL;
17  long int i;
18  double a,
19  dil,
20  rlogl;
21 
22  DEBUG_ENTRY( "ParseBlackbody()" );
23 
24  /* type is blackbody */
25  strcpy( rfield.chSpType[rfield.nspec], "BLACK" );
26 
27  /* these two are not used for this continuum shape */
28  rfield.cutoff[rfield.nspec][0] = 0.;
29  rfield.cutoff[rfield.nspec][1] = 0.;
30 
31  /* get the temperature */
32  i = 5;
33  rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
34 
35  /* there must be at least one number, the temperature */
36  if( lgEOL )
37  {
38  fprintf( ioQQQ, " There must be at least 1 number on a BLACKBODY command line. Sorry.\n" );
39  cdEXIT(EXIT_FAILURE);
40  }
41 
42  /* this is the temperature - make sure its linear in the end
43  * there are two keys, LINEAR and LOG, that could be here,
44  * else choose which is here by which side of 10 */
45  if( (rfield.slope[rfield.nspec] <= 10. && !nMatch("LINE",chCard)) ||
46  nMatch(" LOG",chCard) )
47  {
48  /* log option */
50  }
51 
52  /* check that temp is not too low - could happen if log misused */
53  if( rfield.slope[rfield.nspec] < 1e4 )
54  {
55  fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n",
57  }
58 
59  /* now check that temp not too low - would peak below low
60  * energy limit of the code
61  * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
63  {
64  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",
65  rfield.emm );
66  fprintf( ioQQQ, " Was this intended?\n" );
67  }
68 
69  /* now check that temp not too high - would extend beyond high
70  * energy limit of the code
71  * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
73  {
74  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",
75  rfield.egamry );
76  fprintf( ioQQQ, " Was this intended?\n" );
77  }
78 
79  /* increment continuum indices */
80  ++rfield.nspec;
81  if( rfield.nspec >= LIMSPC )
82  {
83  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
84  cdEXIT(EXIT_FAILURE);
85  }
86 
87  /* also possible to input log(total luminosity)=real log(l) */
88  a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
89  if( lgEOL )
90  {
91  /* there was not a second number on the line; check if LTE or STE */
92  if( nMatch(" LTE",chCard) || nMatch("LTE ",chCard) ||
93  nMatch(" STE",chCard) || nMatch("STE ",chCard) )
94  {
95  /* set energy density to the STE - strict thermodynamic equilibrium - value */
96  strcpy( rfield.chSpNorm[*nqh], "LUMI" );
97 
98  /* need nspec-1 since was incremented above */
99  a = log10(rfield.slope[rfield.nspec-1]);
100  rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
101 
102  /* set radius to very large value if not already set */
103  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
104  if( !radius.lgRadiusKnown )
105  {
106  *ar1 = (realnum)radius.rdfalt;
107  radius.Radius = pow(10.,radius.rdfalt);
108  }
109  strcpy( rfield.chRSpec[*nqh], "SQCM" );
110  rfield.range[*nqh][0] = rfield.emm;
111  rfield.range[*nqh][1] = rfield.egamry;
112  rfield.totpow[*nqh] = rlogl;
113  *nqh = MIN2(*nqh+1,LIMSPC);
114 
115  /* check that stack of shape and luminosity specifications
116  * is parallel, stop if not - this happens is background comes
117  * BETWEEN another set of shape and luminosity commands */
118  if( rfield.nspec != *nqh )
119  {
120  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" );
121  fprintf( ioQQQ, " Sorry.\n" );
122  cdEXIT(EXIT_FAILURE);
123  }
124  }
125 
126  return;
127  }
128 
129  strcpy( rfield.chSpNorm[*nqh], "LUMI" );
130 
131  /* a second number was entered, what was it? */
132  if( nMatch("LUMI",chCard) )
133  {
134  rlogl = a;
135  strcpy( rfield.chRSpec[*nqh], "4 PI" );
136  }
137 
138  else if( nMatch("RADI",chCard) )
139  {
140  /* radius was entered, convert to total lumin */
141  rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nspec-1]);
142  strcpy( rfield.chRSpec[*nqh], "4 PI" );
143  }
144 
145  else if( nMatch("DENS",chCard) )
146  {
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( nMatch("LINE",chCard) || 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  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
157  if( !radius.lgRadiusKnown )
158  {
159  *ar1 = (realnum)radius.rdfalt;
160  radius.Radius = pow(10.,radius.rdfalt);
161  }
162  strcpy( rfield.chRSpec[*nqh], "SQCM" );
163  }
164 
165  else if( nMatch("DILU",chCard) )
166  {
167  /* number is dilution factor, if negative then its log */
168  if( a < 0. )
169  {
170  dil = a;
171  }
172  else
173  {
174  dil = log10(a);
175  }
176  if( dil > 0. )
177  {
178  fprintf( ioQQQ, " Is a dilution factor greater than one physical?\n" );
179  }
180 
181  /* this is LTE energy density */
182  a = log10(rfield.slope[rfield.nspec-1]);
183  rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
184 
185  /* add on dilution factor */
186  rlogl += dil;
187 
188  /* set radius to very large value if not already set */
189  /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
190  if( !radius.lgRadiusKnown )
191  {
192  *ar1 = (realnum)radius.rdfalt;
193  radius.Radius = pow(10.,radius.rdfalt);
194  }
195  strcpy( rfield.chRSpec[*nqh], "SQCM" );
196 
197  }
198  else
199  {
200  fprintf( ioQQQ, " I do not understand what the second number was- keywords are LUMINOSITY, RADIUS, DILUTION, and DENSITY.\n" );
201  cdEXIT(EXIT_FAILURE);
202  }
203 
204  rfield.range[*nqh][0] = rfield.emm;
205  rfield.range[*nqh][1] = rfield.egamry;
206  rfield.totpow[*nqh] = rlogl;
207  /* the second number will be some sort of luminosity */
208  *nqh = MIN2(*nqh+1,LIMSPC);
209 
210  /* check that stack of shape and luminosity specifications
211  * is parallel, stop if not - this happens is background comes
212  * BETWEEN another set of shape and luminosity commands */
213  if( rfield.nspec != *nqh )
214  {
215  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" );
216  fprintf( ioQQQ, " Sorry.\n" );
217  cdEXIT(EXIT_FAILURE);
218  }
219 
220  return;
221 }

Generated for cloudy by doxygen 1.8.3.1