cloudy trunk
Loading...
Searching...
No Matches
parse_interp.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/*ParseInterp parse parameters on interpolate command */
4#include "cddefines.h"
5#include "called.h"
6#include "rfield.h"
7#include "trace.h"
8#include "input.h"
9#include "parser.h"
10
12{
13 DEBUG_ENTRY( "ParseInterp()" );
14
15 /*
16 * this sub reads in the "interpolate" command, and has
17 * logic for the "continue" lines as well.
18 * OUTPUT:
19 * TNU is vector of energies where the grid is defined
20 * TSLOP initially is vector of log fnu at each freq
21 * converted into slopes here too
22 */
23
24 if( rfield.nShape >= LIMSPC )
25 {
26 fprintf( ioQQQ, " Too many spectra entered. Increase LIMSPC\n" );
28 }
29
30 /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */
31 if( !rfield.lgContMalloc[rfield.nShape] )
32 {
33 rfield.tNu[rfield.nShape].resize( NCELL );
34 rfield.tslop[rfield.nShape].resize( NCELL );
35 rfield.tFluxLog[rfield.nShape].resize( NCELL );
36 rfield.lgContMalloc[rfield.nShape] = true;
37 }
38
39 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
40
41 /* read all of the numbers on a line */
42 long npairs = 0;
43
44 /* this is flag saying that all numbers are in */
45 bool lgDONE = false;
46
47 /* this flag says we hit end of command stream */
48 p.m_lgEOF = false;
49 while( !lgDONE && !p.m_lgEOF )
50 {
51 /* keep scanning numbers until we hit eol for current line image */
52 /* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per
53 * Ilse van Bemmel bug report */
54 while( !p.lgEOL() && npairs < NCELL )
55 {
56 rfield.tNu[rfield.nShape][npairs].set( p.FFmtRead() );
57 rfield.tFluxLog[rfield.nShape][npairs] = (realnum)p.FFmtRead();
58 ++npairs;
59 }
60 /* check if we fell through due to too many points, did not hit EOL */
61 if( !p.lgEOL() )
62 {
63 fprintf( ioQQQ, "Too many continuum points were entered.\n" );
64 fprintf( ioQQQ,
65 "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
66 fprintf( ioQQQ,
67 "Increase the value of NCELL in rfield.h.\nSorry.\n" );
69 }
70
71 /* added too many above, so back off */
72 --npairs;
73
74 do
75 {
76 /* read a new line, checking for EOF */
77 p.getline();
78
79 /* option to ignore all lines starting with #, *, or %.
80 * >>chng 06 sep 04 use routine to check for comments */
81 }
82 while( !p.m_lgEOF && p.isComment());
83
84 /* print the line, but only if it is a continue line */
85 if( called.lgTalk && (p.strcmp("CONT")==0) )
86 {
87 p.echo();
88 }
89
90 /* is this a continue line? */
91 if( p.strcmp("CONT") != 0 )
92 {
93 /* we have a line but next command, not continue */
94 lgDONE = true;
95 }
96
97 /* this is another way to hit end of input stream - blank lines */
98 if( p.last() )
99 p.m_lgEOF = true;
100 }
101
102 /* if valid next line, backup one line */
103 if( lgDONE )
104 {
105 input.nRead -= input.iReadWay;
106 }
107
108 /* done reading all of the possible lines */
109 if( npairs < 2 )
110 {
111 fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
113 }
114
115 if( rfield.tNu[rfield.nShape][0].Ryd() == 0. )
116 {
117 /* special case - if first energy is zero then it is low energy */
118 if( rfield.tNu[rfield.nShape][1].Ryd() > 0. )
119 {
120 /* second energy positive, assume linear Ryd */
121 rfield.tNu[rfield.nShape][0].set(rfield.emm);
122 }
123 else if( rfield.tNu[rfield.nShape][1].Ryd() < 0. )
124 {
125 /* second energy negative, assume log of Ryd */
126 rfield.tNu[rfield.nShape][0].set(log10(rfield.emm));
127 }
128 else
129 {
130 /* second energy zero, not allowed */
131 fprintf( ioQQQ,
132 "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
133 rfield.nShape );
135 }
136 }
137
138 /* convert from log(Hz) to Ryd if first nu>5 */
139 if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
140 {
141 for( long i=0; i < npairs; i++ )
142 {
143 rfield.tNu[rfield.nShape][i].set(
144 pow(10.,rfield.tNu[rfield.nShape][i].Ryd())/FR1RYD);
145 }
146 }
147 else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
148 {
149 /* energies entered as logs */
150 for( long i=0; i < npairs; i++ )
151 {
152 rfield.tNu[rfield.nShape][i].set(
153 pow(10.,(double)rfield.tNu[rfield.nShape][i].Ryd()));
154 }
155 }
156 else
157 {
158 /* numbers are linear Rydbergs */
159 for( long i=0; i < npairs; i++ )
160 {
161 if( rfield.tNu[rfield.nShape][i].Ryd() == 0. )
162 {
163 fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
164 i );
166 }
167 }
168 }
169
170 /* option to print debugging file then exit */
171 {
172 /* following should be set true to print information */
173 enum {DEBUG_LOC=false};
174 if( DEBUG_LOC )
175 {
176 for( long i=0; i < npairs; i++ )
177 {
178 fprintf(ioQQQ,"%.4e\t%.3e\n",
179 rfield.tNu[rfield.nShape][i].Ryd(),
180 pow(10.,(double)rfield.tFluxLog[rfield.nShape][i]) * rfield.tNu[rfield.nShape][i].Ryd());
181 }
183 }
184 }
185
186 for( long i=0; i < npairs-1; i++ )
187 {
188 /* check that frequencies are monotonically increasing */
189 if( rfield.tNu[rfield.nShape][i+1].Ryd() <= rfield.tNu[rfield.nShape][i].Ryd() )
190 {
191 fprintf( ioQQQ, "The energies MUST be in increasing order. Energy #%3ld=%10.2e Ryd was greater than or equal to the next one.\nSorry.\n",
192 i, rfield.tNu[rfield.nShape][i].Ryd() );
194 }
195
196 /* TFAC is energy, and TSLOP is slope in f_nu; not photons */
197 rfield.tslop[rfield.nShape][i] = (realnum)((rfield.tFluxLog[rfield.nShape][i+1] -
198 rfield.tFluxLog[rfield.nShape][i])/log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
199 rfield.tNu[rfield.nShape][i].Ryd()));
200 }
201
202 /* now check that array is defined over all energies */
203 if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
204 {
205 /* not defined over low energy part of array */
206 fprintf( ioQQQ,
207 "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
208 fprintf( ioQQQ,
209 " NOTE You may be making a BIG mistake.\n\n" );
210 }
211
212 /* check on IR */
213 if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm )
214 rfield.lgMMok = false;
215
216 if( rfield.tNu[rfield.nShape][0].Ryd() > 1/36. )
217 rfield.lgHPhtOK = false;
218
219 /* gamma ray, EGAMRY is roughly 100MeV */
220 if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.egamry )
221 rfield.lgGamrOK = false;
222
223 /* EnerGammaRay is roughly 100keV; high is gamma ray */
224 if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.EnerGammaRay )
225 rfield.lgXRayOK = false;
226
227 /* renormalize continuum so that flux we will interpolate upon passes through unity
228 * at near 1 Ryd. but first we must find 1 Ryd in the array.
229 * find 1 Ryd, npairs is one less than number of continuum pairs */
230 /* If no flux is defined at 1 Ryd, use the nearest endpoint of the supplied spectrum */
231 long n;
232 for( n=0; n < npairs; n++ )
233 {
234 if( rfield.tNu[rfield.nShape][n].Ryd() > 1. )
235 break;
236 }
237 /* if present, n is now the table point where rfield.tNuRyd[n-1] is <= 1 ryd,
238 * and rfield.tNuRyd[n] > 1 ryd; max(n-1,0) is the nearest endpoint otherwise */
239 realnum fac = rfield.tFluxLog[rfield.nShape][max(n-1,0)];
240
241 for( long i=0; i < npairs; i++ )
242 {
243 rfield.tFluxLog[rfield.nShape][i] -= fac;
244 /*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nShape][i] );*/
245 }
246
247 /* option to print out results at this stage - "trace continuum" */
248 if( trace.lgConBug && trace.lgTrace )
249 {
250 fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
251 npairs );
252 for( long i=0; i < npairs-1; i++ )
253 {
254 fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n",
255 i , rfield.tNu[rfield.nShape][i].Ryd(),
256 rfield.tFluxLog[rfield.nShape][i], rfield.tslop[rfield.nShape][i] );
257 }
258 fprintf( ioQQQ, "%li\t%.4e\t%.4e\n",
259 npairs , rfield.tNu[rfield.nShape][npairs].Ryd(),
260 rfield.tFluxLog[rfield.nShape][npairs]);
261 }
262
263 /* finally check that we are within dynamic range of this cpu */
264 double cmin = log10( FLT_MIN );
265 double cmax = log10( FLT_MAX );
266 bool lgHit = false;
267 for( long i=0; i < npairs; i++ )
268 {
269 if( rfield.tFluxLog[rfield.nShape][i] <= cmin ||
270 rfield.tFluxLog[rfield.nShape][i] >= cmax )
271 {
272 fprintf(ioQQQ,
273 " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
274 fprintf(ioQQQ,
275 " The frequency is %f and the log of the flux is %f.\n\n",
276 rfield.tNu[rfield.nShape][i].Ryd() ,
277 rfield.tFluxLog[rfield.nShape][i]);
278 lgHit = true;
279 }
280 }
281 if( lgHit )
282 {
283 fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
284 fprintf(ioQQQ," NOTE I will try to renormalize it to be within the range of this cpu, but if I crash, this is a likely reason.\n");
285 fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
286 fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
287 fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
288 }
289
290 rfield.ncont[rfield.nShape] = npairs;
291
292 /*>>chng 06 may 17, must insure that rNuRyd is sane through NCELL values, not just
293 * nupper values. This bit when stellar continuum had more cells than cloudy grid,
294 * so npairs is much greater than nupper */
295 /* zero out remainder of array */
296 rfield.tslop[rfield.nShape][npairs-1] = 0.;
297 for( long i=npairs; i < NCELL; i++ )
298 {
299 /* >>chng 05 jul 27, next three indices had been ipspec, changed to nShape
300 * bug caught by Ralf Quast */
301 rfield.tFluxLog[rfield.nShape][i] = 0.;
302 rfield.tslop[rfield.nShape][i] = 0.;
303 rfield.tNu[rfield.nShape][i].set(0.);
304 }
305
306 /* now increment number of continua */
307 ++rfield.nShape;
308
309 return;
310}
t_called called
Definition called.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
long max(int a, long b)
Definition cddefines.h:775
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
bool getline(void)
Definition parser.cpp:164
double FFmtRead(void)
Definition parser.cpp:353
int strcmp(const char *s2)
Definition parser.h:177
bool isComment(void) const
Definition parser.cpp:93
bool last(void) const
Definition parser.h:200
bool lgEOL(void) const
Definition parser.h:98
void echo(void) const
Definition parser.cpp:147
bool m_lgEOF
Definition parser.h:42
t_input input
Definition input.cpp:12
void ParseInterp(Parser &p)
UNUSED const double FR1RYD
Definition physconst.h:195
t_rfield rfield
Definition rfield.cpp:8
const int NCELL
Definition rfield.h:21
const int LIMSPC
Definition rfield.h:18
t_trace trace
Definition trace.cpp:5