cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_powerlawcontinuum.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 /*ParsePowerlawContinuum parse the power law continuum command */
4 #include "cddefines.h"
5 #include "rfield.h"
6 #include "optimize.h"
7 #include "input.h"
8 #include "parse.h"
9 #include "physconst.h"
10 
11 void ParsePowerlawContinuum(char *chCard)
12 {
13  bool lgEOL;
14  long int i;
15 
16  DEBUG_ENTRY( "ParsePowerlawContinuum()" );
17 
18  /* power law with cutoff and X-ray continuum */
19  strcpy( rfield.chSpType[rfield.nspec], "POWER" );
20 
21  /* first parameter is slope of continuum, probably should be negative */
22  i = 5;
23  rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
24  if( lgEOL )
25  {
26  fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" );
27  cdEXIT(EXIT_FAILURE);
28  }
29 
30  if( rfield.slope[rfield.nspec] >= 0. )
31  {
32  fprintf( ioQQQ, " Is the slope of this power law correct?\n" );
33  }
34 
35  /* second optional parameter is high energy cut off */
37  &lgEOL);
38 
39  /* no cutoff if eof hit */
40  if( lgEOL )
41  {
42  /* no extra parameters at all, so put in extreme cutoffs */
43  rfield.cutoff[rfield.nspec][0] = 1e4;
44  rfield.cutoff[rfield.nspec][1] = 1e-4;
45  }
46  else
47  {
48  /* first cutoff was present, check for second */
49  rfield.cutoff[rfield.nspec][1] = FFmtRead(chCard,&i,
50  INPUT_LINE_LENGTH,&lgEOL);
51  if( lgEOL )
52  rfield.cutoff[rfield.nspec][1] = 1e-4;
53  }
54 
55  /* check that energies were entered in the correct order */
57  {
58  fprintf( ioQQQ, " The optional cutoff energies do not appear to be in the correct order. Check Hazy.\n" );
59  cdEXIT(EXIT_FAILURE);
60  }
61 
62  /* check for keyword KELVIN to interpret cutoff energies as degrees */
63  if( nMatch("KELV",chCard) )
64  {
65  /* temps are log if first le 10 */
66  if( rfield.cutoff[rfield.nspec][0] <= 10. )
67  rfield.cutoff[rfield.nspec][0] = pow(10.,rfield.cutoff[rfield.nspec][0]);
68  if( rfield.cutoff[rfield.nspec][1] <= 10. )
69  rfield.cutoff[rfield.nspec][1] = pow(10.,rfield.cutoff[rfield.nspec][1]);
72  }
73 
74  if( rfield.cutoff[rfield.nspec][0] < 0. ||
75  rfield.cutoff[rfield.nspec][1] < 0. )
76  {
77  fprintf( ioQQQ, " A negative cutoff energy is not physical. Sorry.\n" );
78  cdEXIT(EXIT_FAILURE);
79  }
80 
81  if( rfield.cutoff[rfield.nspec][1] == 0. &&
82  rfield.slope[rfield.nspec] <= -1. )
83  {
84  fprintf( ioQQQ, " A power-law with this slope, and no low energy cutoff, may have an unphysically large\n brightness temperature in the radio.\n" );
85  }
86 
87  /* vary option */
88  if( optimize.lgVarOn )
89  {
90  /* pointer to where to write */
92  if( nMatch("ARYB",chCard) )
93  {
94  /* this test is for key "varyb", meaning to vary second parameter
95  * the cutoff temperature
96  * this is the number of parameters to feed onto the input line */
98  /* vary ?? */
99  sprintf( optimize.chVarFmt[optimize.nparm],
100  "POWER LAW %f KELVIN%%f %%f",
102  /* param is linear scale factor */
104  TE1RYD);
106  TE1RYD);
107  optimize.vincr[optimize.nparm] = 0.2f;
108  }
109  else if( nMatch("ARYC",chCard) )
110  {
111  /* the keyword was "varyc"
112  * this is the number of parameters to feed onto the input line */
114 
115  /* vary ?? */
116  sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW%f %f KELVIN %%f )",
118  log10(rfield.cutoff[rfield.nspec][0]* TE1RYD) );
119 
121  TE1RYD);
122  optimize.vincr[optimize.nparm] = 0.2f;
123  }
124  else
125  {
126  /* vary the first parameter only, but still are two more
127  * this is the number of parameters to feed onto the input line */
129  strcpy( optimize.chVarFmt[optimize.nparm],
130  "POWER LAW KELVIN%f %f %f" );
131  /* param is slope of the power law */
134  TE1RYD);
136  TE1RYD);
137  optimize.vincr[optimize.nparm] = 0.2f;
138  }
139  ++optimize.nparm;
140  }
141 
142  /*>>chng 06 nov 10, BUGFIX, nspec was incremented before previous branch
143  * and so crashed with log10 0 since nspec was beyond set values
144  * caused fpe domain function error with log 0
145  * fpe caught by Pavel Abolmasov */
146  ++rfield.nspec;
147  if( rfield.nspec >= LIMSPC )
148  {
149  fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
150  cdEXIT(EXIT_FAILURE);
151  }
152 
153  return;
154 }

Generated for cloudy by doxygen 1.8.4