cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_interp.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 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  // we reserve space for 1000 table points, no worries if that is too small though...
31  ASSERT( rfield.tNu[rfield.nShape].empty() );
32  rfield.tNu[rfield.nShape].reserve( 1000 );
33  ASSERT( rfield.tslop[rfield.nShape].empty() );
34  rfield.tslop[rfield.nShape].reserve( 1000 );
35  ASSERT( rfield.tFluxLog[rfield.nShape].empty() );
36  rfield.tFluxLog[rfield.nShape].reserve( 1000 );
37 
38  strcpy( rfield.chSpType[rfield.nShape], "INTER" );
39 
40  /* read all of the numbers on a line */
41  long npairs = 0;
42 
43  /* this is flag saying that all numbers are in */
44  bool lgDONE = false;
45 
46  /* this flag says we hit end of command stream */
47  p.m_lgEOF = false;
48  while( !lgDONE && !p.m_lgEOF )
49  {
50  /* keep scanning numbers until we hit eol for current line image */
51  while( true )
52  {
53  Energy E( p.FFmtRead() );
54  realnum FluxLog = (realnum)p.FFmtRead();
55  if( p.lgEOL() )
56  break;
57  rfield.tNu[rfield.nShape].push_back( E );
58  rfield.tFluxLog[rfield.nShape].push_back( FluxLog );
59  ++npairs;
60  }
61 
62  do
63  {
64  /* read a new line, checking for EOF */
65  p.getline();
66 
67  /* option to ignore all lines starting with #, *, or %.
68  * >>chng 06 sep 04 use routine to check for comments */
69  }
70  while( !p.m_lgEOF && p.isComment());
71 
72  /* print the line, but only if it is a continue line */
73  if( called.lgTalk && p.hasCommand("CONT") )
74  {
75  p.echo();
76  }
77 
78  /* is this a continue line? */
79  if( ! p.hasCommand("CONT") )
80  {
81  /* we have a line but next command, not continue */
82  lgDONE = true;
83  }
84 
85  /* this is another way to hit end of input stream - blank lines */
86  if( p.last() )
87  p.m_lgEOF = true;
88  }
89 
90  /* if valid next line, backup one line */
91  if( lgDONE )
92  {
94  }
95 
96  /* done reading all of the possible lines */
97  if( npairs < 2 )
98  {
99  fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
101  }
102 
103  if( rfield.tNu[rfield.nShape][0].Ryd() == 0. )
104  {
105  /* special case - if first energy is zero then it is low energy */
106  if( rfield.tNu[rfield.nShape][1].Ryd() > 0. )
107  {
108  /* second energy positive, assume linear Ryd */
109  rfield.tNu[rfield.nShape][0].set(rfield.emm());
110  }
111  else if( rfield.tNu[rfield.nShape][1].Ryd() < 0. )
112  {
113  /* second energy negative, assume log of Ryd */
114  rfield.tNu[rfield.nShape][0].set(log10(rfield.emm()));
115  }
116  else
117  {
118  /* second energy zero, not allowed */
119  fprintf( ioQQQ,
120  "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
121  rfield.nShape );
123  }
124  }
125 
126  /* convert from log(Hz) to Ryd if first nu>5 */
127  if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
128  {
129  for( long i=0; i < npairs; i++ )
130  {
131  if( rfield.tNu[rfield.nShape][i].Ryd()>300.)
132  {
133  fprintf(ioQQQ,"DISASTER in INTERPPLATE command syntax\n");
134  fprintf(ioQQQ,"log of photon energy in Hz was entered, but the value, %g is too large.\nSorry.\n",
135  rfield.tNu[rfield.nShape][i].Ryd());
137  }
138  rfield.tNu[rfield.nShape][i].set(
139  exp10(rfield.tNu[rfield.nShape][i].Ryd())/FR1RYD);
140  }
141  }
142  else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
143  {
144  /* energies entered as logs */
145  for( long i=0; i < npairs; i++ )
146  {
147  rfield.tNu[rfield.nShape][i].set(
148  exp10((double)rfield.tNu[rfield.nShape][i].Ryd()));
149  }
150  }
151  else
152  {
153  /* numbers are linear Rydbergs */
154  for( long i=0; i < npairs; i++ )
155  {
156  if( rfield.tNu[rfield.nShape][i].Ryd() == 0. )
157  {
158  fprintf( ioQQQ, "An energy of zero was entered for pair %3ld in INTERPOLATE and is not allowed.\nSorry\n",
159  i );
161  }
162  }
163  }
164 
165  /* option to print debugging file then exit */
166  {
167  /* following should be set true to print information */
168  enum {DEBUG_LOC=false};
169  if( DEBUG_LOC )
170  {
171  for( long i=0; i < npairs; i++ )
172  {
173  fprintf(ioQQQ,"%.4e\t%.3e\n",
174  rfield.tNu[rfield.nShape][i].Ryd(),
175  exp10((double)rfield.tFluxLog[rfield.nShape][i]) * rfield.tNu[rfield.nShape][i].Ryd());
176  }
178  }
179  }
180 
181  for( long i=0; i < npairs-1; i++ )
182  {
183  /* check that frequencies are monotonically increasing */
184  if( rfield.tNu[rfield.nShape][i+1].Ryd() <= rfield.tNu[rfield.nShape][i].Ryd() )
185  {
186  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",
187  i, rfield.tNu[rfield.nShape][i].Ryd() );
189  }
190 
191  /* TFAC is energy, and TSLOP is slope in f_nu; not photons */
192  rfield.tslop[rfield.nShape].push_back(
194  log10(rfield.tNu[rfield.nShape][i+1].Ryd()/rfield.tNu[rfield.nShape][i].Ryd()))
195  );
196  }
197  rfield.tslop[rfield.nShape].push_back( 0. );
198 
199  /* now check that array is defined over all energies */
200  if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm() )
201  {
202  /* not defined over low energy part of array */
203  fprintf( ioQQQ,
204  "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
205  fprintf( ioQQQ,
206  " NOTE You may be making a BIG mistake.\n\n" );
207  }
208 
209  /* check on IR */
210  if( rfield.tNu[rfield.nShape][0].Ryd() > rfield.emm() )
211  rfield.lgMMok = false;
212 
213  if( rfield.tNu[rfield.nShape][0].Ryd() > 1/36. )
214  rfield.lgHPhtOK = false;
215 
216  /* gamma ray, EGAMRY is roughly 100MeV */
217  if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.egamry() )
218  rfield.lgGamrOK = false;
219 
220  /* EnerGammaRay is roughly 100keV; high is gamma ray */
221  if( rfield.tNu[rfield.nShape][npairs-1].Ryd() < rfield.EnerGammaRay )
222  rfield.lgXRayOK = false;
223 
224  /* renormalize continuum so that flux we will interpolate upon passes through unity
225  * at near 1 Ryd. but first we must find 1 Ryd in the array.
226  * find 1 Ryd, npairs is one less than number of continuum pairs */
227  /* If no flux is defined at 1 Ryd, use the nearest endpoint of the supplied spectrum */
228  long n;
229  for( n=0; n < npairs; n++ )
230  {
231  if( rfield.tNu[rfield.nShape][n].Ryd() > 1. )
232  break;
233  }
234  /* if present, n is now the table point where rfield.tNuRyd[n-1] is <= 1 ryd,
235  * and rfield.tNuRyd[n] > 1 ryd; max(n-1,0) is the nearest endpoint otherwise */
236  realnum fac = rfield.tFluxLog[rfield.nShape][max(n-1,0)];
237 
238  for( long i=0; i < npairs; i++ )
239  {
240  rfield.tFluxLog[rfield.nShape][i] -= fac;
241  /*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nShape][i] );*/
242  }
243 
244  /* option to print out results at this stage - "trace continuum" */
245  if( trace.lgConBug && trace.lgTrace )
246  {
247  fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
248  npairs );
249  for( long i=0; i < npairs-1; i++ )
250  {
251  fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n",
252  i , rfield.tNu[rfield.nShape][i].Ryd(),
254  }
255  fprintf( ioQQQ, "%li\t%.4e\t%.4e\n",
256  npairs , rfield.tNu[rfield.nShape][npairs].Ryd(),
257  rfield.tFluxLog[rfield.nShape][npairs]);
258  }
259 
260  /* finally check that we are within dynamic range of this cpu */
261  double cmin = log10( FLT_MIN );
262  double cmax = log10( FLT_MAX );
263  bool lgHit = false;
264  for( long i=0; i < npairs; i++ )
265  {
266  if( rfield.tFluxLog[rfield.nShape][i] <= cmin ||
267  rfield.tFluxLog[rfield.nShape][i] >= cmax )
268  {
269  fprintf(ioQQQ,
270  " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
271  fprintf(ioQQQ,
272  " The photon energy is %f Rydberg and the log of the flux is %f.\n\n",
273  rfield.tNu[rfield.nShape][i].Ryd() ,
275  lgHit = true;
276  }
277  }
278  if( lgHit )
279  {
280  fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
281  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");
282  fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
283  fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
284  fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
285  }
286 
287  rfield.ncont[rfield.nShape] = npairs;
288 
289  ASSERT( rfield.tNu[rfield.nShape].size() == (size_t)npairs );
290  ASSERT( rfield.tFluxLog[rfield.nShape].size() == (size_t)npairs );
291  ASSERT( rfield.tslop[rfield.nShape].size() == (size_t)npairs );
292 
293  /* now increment number of continua */
294  ++rfield.nShape;
295 
296  return;
297 }
double emm() const
Definition: mesh.h:84
bool hasCommand(const char *s2)
Definition: parser.cpp:705
bool lgGamrOK
Definition: rfield.h:440
void echo(void) const
Definition: parser.cpp:167
double FFmtRead(void)
Definition: parser.cpp:438
double exp10(double x)
Definition: cddefines.h:1383
realnum EnerGammaRay
Definition: rfield.h:446
t_input input
Definition: input.cpp:12
long int nRead
Definition: input.h:62
vector< Energy > tNu[LIMSPC]
Definition: rfield.h:316
vector< realnum > tFluxLog[LIMSPC]
Definition: rfield.h:319
FILE * ioQQQ
Definition: cddefines.cpp:7
long ncont[LIMSPC]
Definition: rfield.h:321
vector< realnum > tslop[LIMSPC]
Definition: rfield.h:317
bool lgTalk
Definition: called.h:12
Definition: parser.h:42
const int LIMSPC
Definition: rfield.h:20
t_trace trace
Definition: trace.cpp:5
bool lgConBug
Definition: trace.h:97
bool lgHPhtOK
Definition: rfield.h:440
bool lgTrace
Definition: trace.h:12
t_rfield rfield
Definition: rfield.cpp:9
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
bool last(void) const
Definition: parser.h:196
long max(int a, long b)
Definition: cddefines.h:821
#define cdEXIT(FAIL)
Definition: cddefines.h:484
bool lgXRayOK
Definition: rfield.h:440
#define ASSERT(exp)
Definition: cddefines.h:617
long int iReadWay
Definition: input.h:62
bool getline()
Definition: parser.cpp:249
Definition: energy.h:9
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
bool lgMMok
Definition: rfield.h:440
double egamry() const
Definition: mesh.h:88
bool lgEOL(void) const
Definition: parser.h:103
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
long int nShape
Definition: rfield.h:308
bool m_lgEOF
Definition: parser.h:53
char chSpType[LIMSPC][6]
Definition: rfield.h:333
t_called called
Definition: called.cpp:4
bool isComment(void) const
Definition: parser.cpp:113
void ParseInterp(Parser &p)
#define EXIT_SUCCESS
Definition: cddefines.h:166