cloudy trunk
Loading...
Searching...
No Matches
parse_powerlawcontinuum.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/*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 "parser.h"
9#include "physconst.h"
10
12{
13 DEBUG_ENTRY( "ParsePowerlawContinuum()" );
14
15 /* power law with cutoff and X-ray continuum */
16 strcpy( rfield.chSpType[rfield.nShape], "POWER" );
17
18 /* first parameter is slope of continuum, probably should be negative */
19 rfield.slope[rfield.nShape] = p.FFmtRead();
20 if( p.lgEOL() )
21 {
22 fprintf( ioQQQ, " There should have been a number on this line. Sorry.\n" );
24 }
25
26 if( rfield.slope[rfield.nShape] >= 0. )
27 {
28 fprintf( ioQQQ, " Is the slope of this power law correct?\n" );
29 }
30
31 /* second optional parameter is high energy cut off */
32 rfield.cutoff[rfield.nShape][0] = p.FFmtRead();
33
34 /* no cutoff if eof hit */
35 if( p.lgEOL() )
36 {
37 /* no extra parameters at all, so put in extreme cutoffs */
38 rfield.cutoff[rfield.nShape][0] = 1e4;
39 rfield.cutoff[rfield.nShape][1] = 1e-4;
40 }
41 else
42 {
43 /* first cutoff was present, check for second */
44 rfield.cutoff[rfield.nShape][1] = p.FFmtRead();
45 if( p.lgEOL() )
46 rfield.cutoff[rfield.nShape][1] = 1e-4;
47 }
48
49 /* check that energies were entered in the correct order */
50 if( rfield.cutoff[rfield.nShape][1] > rfield.cutoff[rfield.nShape][0] )
51 {
52 fprintf( ioQQQ, " The optional cutoff energies do not appear to be in the correct order. Check Hazy.\n" );
54 }
55
56 /* check for keyword KELVIN to interpret cutoff energies as degrees */
57 if( p.nMatch("KELV") )
58 {
59 /* temps are log if first le 10 */
60 if( rfield.cutoff[rfield.nShape][0] <= 10. || p.nMatch(" LOG") )
61 rfield.cutoff[rfield.nShape][0] = pow(10.,rfield.cutoff[rfield.nShape][0]);
62 if( rfield.cutoff[rfield.nShape][1] <= 10. || p.nMatch(" LOG") )
63 rfield.cutoff[rfield.nShape][1] = pow(10.,rfield.cutoff[rfield.nShape][1]);
64 rfield.cutoff[rfield.nShape][0] /= TE1RYD;
65 rfield.cutoff[rfield.nShape][1] /= TE1RYD;
66 }
67
68 if( rfield.cutoff[rfield.nShape][0] < 0. ||
69 rfield.cutoff[rfield.nShape][1] < 0. )
70 {
71 fprintf( ioQQQ, " A negative cutoff energy is not physical. Sorry.\n" );
73 }
74
75 if( rfield.cutoff[rfield.nShape][1] == 0. &&
76 rfield.slope[rfield.nShape] <= -1. )
77 {
78 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" );
79 }
80
81 /* vary option */
82 if( optimize.lgVarOn )
83 {
84 /* pointer to where to write */
85 optimize.nvfpnt[optimize.nparm] = input.nRead;
86 if( p.nMatch("VARYB") )
87 {
88 /* this test is for key "varyb", meaning to vary second parameter
89 * the cutoff temperature
90 * this is the number of parameters to feed onto the input line */
91 optimize.nvarxt[optimize.nparm] = 1;
92 sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW %f KELVIN %%f %f LOG",
93 rfield.slope[rfield.nShape],
94 log10(rfield.cutoff[rfield.nShape][1]*TE1RYD) );
95 optimize.vparm[0][optimize.nparm] =
96 (realnum)log10(rfield.cutoff[rfield.nShape][0]*TE1RYD);
97 optimize.varang[optimize.nparm][0] = (realnum)log10(rfield.cutoff[rfield.nShape][1]*TE1RYD);
98 optimize.varang[optimize.nparm][1] = FLT_MAX;
99 optimize.vincr[optimize.nparm] = 0.2f;
100 }
101 else if( p.nMatch("VARYC") )
102 {
103 /* the keyword was "varyc"
104 * this is the number of parameters to feed onto the input line */
105 optimize.nvarxt[optimize.nparm] = 1;
106 sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW %f KELVIN %f %%f LOG",
107 rfield.slope[rfield.nShape],
108 log10(rfield.cutoff[rfield.nShape][0]*TE1RYD) );
109 optimize.vparm[0][optimize.nparm] =
110 (realnum)log10(rfield.cutoff[rfield.nShape][1]*TE1RYD);
111 optimize.varang[optimize.nparm][0] = -FLT_MAX;
112 optimize.varang[optimize.nparm][1] = (realnum)log10(rfield.cutoff[rfield.nShape][0]*TE1RYD);
113 optimize.vincr[optimize.nparm] = 0.2f;
114 }
115 else
116 {
117 /* vary the first parameter only, but still are two more
118 * this is the number of parameters to feed onto the input line */
119 optimize.nvarxt[optimize.nparm] = 1;
120 sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW %%f KELVIN %f %f LOG",
121 log10(rfield.cutoff[rfield.nShape][0]*TE1RYD),
122 log10(rfield.cutoff[rfield.nShape][1]*TE1RYD) );
123 optimize.vparm[0][optimize.nparm] = (realnum)rfield.slope[rfield.nShape];
124 optimize.vincr[optimize.nparm] = 0.2f;
125 }
126 ++optimize.nparm;
127 }
128
129 /*>>chng 06 nov 10, BUGFIX, nShape was incremented before previous branch
130 * and so crashed with log10 0 since nSpage was beyond set values
131 * caused fpe domain function error with log 0
132 * fpe caught by Pavel Abolmasov */
133 ++rfield.nShape;
134 if( rfield.nShape >= LIMSPC )
135 {
136 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
138 }
139
140 return;
141}
FILE * ioQQQ
Definition cddefines.cpp:7
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double FFmtRead(void)
Definition parser.cpp:353
bool nMatch(const char *chKey) const
Definition parser.h:135
bool lgEOL(void) const
Definition parser.h:98
t_input input
Definition input.cpp:12
t_optimize optimize
Definition optimize.cpp:5
void ParsePowerlawContinuum(Parser &p)
UNUSED const double TE1RYD
Definition physconst.h:183
t_rfield rfield
Definition rfield.cpp:8
const int LIMSPC
Definition rfield.h:18