cloudy  trunk
 All Data Structures 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-2008 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 "parse.h"
10 
11 void ParseInterp(char *chCard,
12  bool *lgEOF)
13 {
14  char chLab4[5];
15  bool lgDONE,
16  lgEOL,
17  lgHit;
18  long int i,
19  k,
20  npairs;
21  double cmax,
22  cmin,
23  fac;
24 
25  DEBUG_ENTRY( "ParseInterp()" );
26 
27  /*
28  * this sub reads in the "interpolate" command, and has
29  * logic for the "continue" lines as well.
30  * OUTPUT:
31  * TNU is vector of energies where the grid is defined
32  * TSLOP initially is vector of log fnu at each freq
33  * converted into slopes here too
34  */
35 
36  if( rfield.nspec >= LIMSPC-1 )
37  {
38  fprintf( ioQQQ, " Too many spectra entered. Increase LIMSPC\n" );
39  cdEXIT(EXIT_FAILURE);
40  }
41  /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */
43  {
44  rfield.tNuRyd[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
45  rfield.tslop[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
46  rfield.tFluxLog[rfield.nspec] = (realnum*)MALLOC((size_t)(NCELL*sizeof(realnum)) );
48  }
49 
50  strcpy( rfield.chSpType[rfield.nspec], "INTER" );
51 
52  /* read all of the numbers on a line */
53  npairs = 0;
54 
55  /* this is flag saying that all numbers are in */
56  lgDONE = false;
57 
58  /* this flag says we hit end of command stream */
59  *lgEOF = false;
60  while( !lgDONE && !*lgEOF )
61  {
62  i = 5;
63  lgEOL = false;
64 
65  /* keep scanning numbers until we hit eol for current line image */
66  /* >>chng 01 aug 10, add check on npairs not exceeding NC_ELL as per
67  * Ilse van Bemmel bug report */
68  /*while( !lgEOL && npairs<rfield.nupper )*/
69  while( !lgEOL && npairs<NCELL )
70  {
71  rfield.tNuRyd[rfield.nspec][npairs] =
72  (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
73  rfield.tFluxLog[rfield.nspec][npairs] =
74  (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
75  ++npairs;
76  }
77  /* check if we fell through due to too many points, did not hit EOL */
78  if( !lgEOL )
79  {
80  fprintf( ioQQQ, "Too many continuum points were entered.\n" );
81  fprintf( ioQQQ,
82  "The current logic limits the number of possible points to the value of NCELL, which is %i.\n",NCELL );
83  fprintf( ioQQQ,
84  "Increase the value of NCELL in rfield.h.\nSorry.\n" );
85  cdEXIT(EXIT_FAILURE);
86  }
87 
88  /* added too many above, so back off */
89  --npairs;
90 
91  /* read a new line, checking for EOF */
92  input_readarray(chCard,lgEOF);
93 
94  /* option to ignore all lines starting with #, *, or %.
95  * >>chng 06 sep 04 use routine to check for comments */
96  while( !*lgEOF && lgInputComment(chCard)/*((chCard[0] == '#' || chCard[0] == '*') || chCard[0] == '%')*/ )
97  {
98  input_readarray(chCard,lgEOF);
99  }
100 
101  /* get capital first four char of chCard */
102  cap4( chLab4 , chCard );
103 
104  /* print the line, but only if it is a continue line */
105  if( called.lgTalk && (strncmp(chLab4,"CONT",4)==0) )
106  {
107  fprintf( ioQQQ, " * ");
108  /* next logic will print command plus spaces to right justified * */
109  k=0;
110  while( chCard[k]!='\0' )
111  {
112  fprintf(ioQQQ,"%c",chCard[k]);
113  ++k;
114  }
115  while( k<80 )
116  {
117  fprintf(ioQQQ,"%c",' ');
118  ++k;
119  }
120  fprintf( ioQQQ,"*\n");
121  }
122 
123  /* now convert entire line to caps */
124  caps(chCard);
125 
126  /* is this a continue line? */
127  if( strncmp(chCard,"CONT",4) != 0 )
128  {
129  /* we have a line but next command, not continue */
130  lgDONE = true;
131  }
132 
133  /* this is another way to hit end of input stream - blank lines */
134  if( chCard[0] == ' ' )
135  *lgEOF = true;
136  }
137 
138  /* if valid next line, backup one line */
139  if( lgDONE )
140  {
142  }
143 
144  /* done reading all of the possible lines */
145  --npairs;
146  if( npairs < 1 )
147  {
148  fprintf( ioQQQ, "There must be at least 2 pairs to interpolate,\nSorry\n" );
149  cdEXIT(EXIT_FAILURE);
150  }
151 
152  /* >>chng 02 apr 19, from > to >=, test did not trigger due to overrun protection
153  * in main loop */
154  else if( npairs >= NCELL - 2 )
155  {
156  fprintf( ioQQQ, " Too many continuum points entered.\n" );
157  fprintf( ioQQQ,
158  "The current logic limits the number of possible points to the value of NCELL, which is %i.\nSorry.\n",NCELL );
159  fprintf( ioQQQ, " Increase NCELL (in cddefines.h) to more than the number of continuum points.\n" );
160  cdEXIT(EXIT_FAILURE);
161  }
162 
163  if( rfield.tNuRyd[rfield.nspec][0] == 0. )
164  {
165  /* special case - if first energy is zero then it is low energy */
166  if( rfield.tNuRyd[rfield.nspec][1] > 0. )
167  {
168  /* second energy positive, assume linear Ryd */
170  }
171  else if( rfield.tNuRyd[rfield.nspec][1] < 0. )
172  {
173  /* second energy negative, assume log of Ryd */
174  rfield.tNuRyd[rfield.nspec][0] = (realnum)log10(rfield.emm);
175  }
176  else
177  {
178  /* second energy zero, not allowed */
179  fprintf( ioQQQ,
180  "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
181  rfield.nspec );
182  cdEXIT(EXIT_FAILURE);
183  }
184  }
185 
186  /* convert from log(Hz) to Ryd if first nu>5 */
187  if( rfield.tNuRyd[rfield.nspec][0] >= 5. )
188  {
189  for( i=0; i < (npairs + 1); i++ )
190  {
191  rfield.tNuRyd[rfield.nspec][i] =
192  (realnum)pow(10.,rfield.tNuRyd[rfield.nspec][i] - 15.517);
193  }
194  }
195 
196  else if( rfield.tNuRyd[rfield.nspec][0] < 0. )
197  {
198  /* energies entered as logs */
199  for( i=0; i < (npairs + 1); i++ )
200  {
201  rfield.tNuRyd[rfield.nspec][i] = (realnum)pow(10.,(double)rfield.tNuRyd[rfield.nspec][i]);
202  }
203  }
204 
205  else
206  {
207  /* numbers are linear Rydbergs */
208  for( i=0; i < (npairs + 1); i++ )
209  {
210  if( rfield.tNuRyd[rfield.nspec][i] == 0. )
211  {
212  fprintf( ioQQQ, "An energy of zero was entered for element%3ld in INTERPOLATE and is not allowed.\nSorry\n",
213  i );
214  cdEXIT(EXIT_FAILURE);
215  }
216  }
217  }
218 
219  /* option to print debugging file then exit */
220  {
221  /* following should be set true to print information */
222  enum {DEBUG_LOC=false};
223  if( DEBUG_LOC )
224  {
225  for( i=0; i < npairs; i++ )
226  {
227  fprintf(ioQQQ,"%.4e\t%.3e\n",
229  pow(10.,(double)rfield.tFluxLog[rfield.nspec][i]) * rfield.tNuRyd[rfield.nspec][i]);
230  }
231  cdEXIT(EXIT_SUCCESS);
232  }
233  }
234 
235  for( i=0; i < npairs; i++ )
236  {
237  /* check that frequencies are monotonically increasing */
238  if( rfield.tNuRyd[rfield.nspec][i+1] <= rfield.tNuRyd[rfield.nspec][i] )
239  {
240  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",
241  i, rfield.tNuRyd[rfield.nspec][i] );
242  cdEXIT(EXIT_FAILURE);
243  }
244 
245  /* TFAC is energy, and TSLOP is slope in f_nu; not photons */
247  rfield.tFluxLog[rfield.nspec][i])/log10(rfield.tNuRyd[rfield.nspec][i+1]/
248  rfield.tNuRyd[rfield.nspec][i]));
249  }
250 
251  /*>>chng 06 may 17, must insure that rNuRyd is sane through NCELL values, not just
252  * nupper values. This bit when stellar continuum had more cells than cloudy grid,
253  * so npairs is much greater than nupper */
254  /*if( npairs + 2 < rfield.nupper )*/
255  if( npairs + 2 < NCELL )
256  {
257  /* zero out remainder of array */
258  /*>>chng 06 may 17, same as above */
259  /*for( i=npairs + 1; i < rfield.nupper; i++ )*/
260  for( i=npairs + 1; i < NCELL; i++ )
261  {
262  rfield.tNuRyd[rfield.nspec][i] = 0.;
263  }
264  }
265 
266  /* now check that array is defined over all energies */
267  if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
268  {
269  /* not defined over low energy part of array */
270  fprintf( ioQQQ,
271  "\n NOTE The incident continuum was not defined over the entire energy range. Some energies are set to zero.\n" );
272  fprintf( ioQQQ,
273  " NOTE You may be making a BIG mistake.\n\n" );
274  }
275 
276  /* check on IR */
277  if( rfield.tNuRyd[rfield.nspec][0] > rfield.emm )
278  rfield.lgMMok = false;
279 
280  if( rfield.tNuRyd[rfield.nspec][0] > 1/36. )
281  rfield.lgHPhtOK = false;
282 
283  /* gamma ray, EGAMRY is roughly 100MeV */
284  if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.egamry )
285  rfield.lgGamrOK = false;
286 
287  /* EnerGammaRay is roughly 100keV; high is gamma ray */
288  if( rfield.tNuRyd[rfield.nspec][npairs] < rfield.EnerGammaRay )
289  rfield.lgXRayOK = false;
290 
291  /* find min and max of continuum */
292  cmax = -38.;
293  cmin = 28;
294 
295  /* rfield.tFluxLog can be very large or small */
296  /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
297  for( i=0; i <= npairs; i++ )
298  {
299  cmax = MAX2(cmax,rfield.tFluxLog[rfield.nspec][i]);
300  cmin = MIN2(cmin,rfield.tFluxLog[rfield.nspec][i]);
301  }
302 
303  /* check on dynamic range of input continuum */
304  if( cmax - cmin > 74. )
305  {
306  fprintf( ioQQQ, "The dynamic range of the specified continuum is too large.\nSorry.\n" );
307  cdEXIT(EXIT_FAILURE);
308  }
309 
310  if( cmin < -37. )
311  {
312  fac = -cmin - 37.;
313  /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
314  for( i=0; i <= npairs; i++ )
315  {
316  rfield.tFluxLog[rfield.nspec][i] += (realnum)fac;
317  }
318  }
319 
320  else if( cmax > 37. )
321  {
322  fac = cmax - 37.;
323  /* >>chng 01 jul 11, from <npairs to <=npairs, [npairs] is highest point in table */
324  for( i=0; i <= npairs; i++ )
325  {
326  rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac;
327  }
328  }
329 
330  /* option to print out results at this stage - "trace continuum" */
331  if( trace.lgConBug && trace.lgTrace )
332  {
333  fprintf( ioQQQ, " Table for this continuum;\ni\tTNU\tTFAC\tTSLOP, npairs=%li\n",
334  npairs );
335  for( i=0; i < npairs; i++ )
336  {
337  fprintf( ioQQQ, "%li\t%.4e\t%.4e\t%.4e\n",
338  i , rfield.tNuRyd[rfield.nspec][i],
340  }
341  i = npairs;
342  fprintf( ioQQQ, "%li\t%.4e\t%.4e\n",
343  i , rfield.tNuRyd[rfield.nspec][i],
345  }
346 
347  /* renormalize continuum so that flux we will interpolate upon passes through unity
348  * at near 1 Ryd. but first we must find 1 Ryd in the array.
349  * find 1 Ryd, npairs is one less than number of continuum pairs */
350  i = 0;
351  while( rfield.tNuRyd[rfield.nspec][i] <= 1. && i < npairs-1 )
352  {
353  i += 1;
354  }
355  /* i is now the table point where rfield.tNuRyd[i-1] is <= 1 ryd,
356  * and rfield.tNuRyd[i] > 1 ryd */
357 
358  /* following is impossible but sanity check */
359  if( i < 1 )
360  {
361  fprintf( ioQQQ, "ParseInput finds insane i after rfield.tNuRyd loop\n");
362  ShowMe();
363  cdEXIT(EXIT_FAILURE);
364  }
365 
366  /* do the renormalization so 1 ryd is about unity */
367  fac = rfield.tFluxLog[rfield.nspec][i-1];
368  for( i=0; i <= npairs; i++ )
369  {
370  rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac;
371  /*fprintf(ioQQQ,"DEBUG parse %li %e \n", i , rfield.tFluxLog[rfield.nspec][i] );*/
372  }
373 
374  /* finally check that we are within dynamic range of this cpu */
375  cmin = log10( FLT_MIN );
376  cmax = log10( FLT_MAX );
377  lgHit = false;
378  for( i=0; i <= npairs; i++ )
379  {
380  if( rfield.tFluxLog[rfield.nspec][i] <= cmin ||
381  rfield.tFluxLog[rfield.nspec][i] >= cmax )
382  {
383  fprintf(ioQQQ,
384  " The log of the flux specified in interpolate pair %li is not within dynamic range of this CPU - please rescale.\n",i);
385  fprintf(ioQQQ,
386  " The frequency is %f and the log of the flux is %f.\n\n",
387  rfield.tNuRyd[rfield.nspec][i] ,
389  lgHit = true;
390  }
391  }
392  if( lgHit )
393  {
394  fprintf(ioQQQ,"\n NOTE The log of the flux given in an interpolate command is outside the range of this cpu.\n");
395  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");
396  fprintf(ioQQQ," NOTE Note that the interpolate command only is used for the shape of the continuum.\n");
397  fprintf(ioQQQ," NOTE The order of magnitude of the flux is not used in any way.\n");
398  fprintf(ioQQQ," NOTE For safety this could be of order unity.\n\n");
399  }
400 
401  /* zero out remainder of array */
402  for( i=npairs+1; i < rfield.nupper; i++ )
403  {
404  /* >>chng 05 jul 27, next three indices had been ipspec, changed to nspec
405  * bug caught by Ralf Quast */
406  rfield.tFluxLog[rfield.nspec][i] = 0.;
407  rfield.tslop[rfield.nspec][i] = 0.;
408  rfield.tNuRyd[rfield.nspec][i] = 0.;
409  }
410 
411  /* now increment number of continua */
412  ++rfield.nspec;
413  if( rfield.nspec >= LIMSPC )
414  {
415  /* too many continua were entered */
416  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
417  cdEXIT(EXIT_FAILURE);
418  }
419 
420  return;
421 }

Generated for cloudy by doxygen 1.8.1.1