1 | //--------------------------------------------------------------------------
|
---|
2 | // File and Version Information:
|
---|
3 | // $Id: radspec.cc,v 1.4 2000-04-13 14:10:44 ansari Exp $
|
---|
4 | //
|
---|
5 | // Description:
|
---|
6 | // Aim of the class: To give the energy density
|
---|
7 | // The unity used here is W/m^2/Hz/sr
|
---|
8 | //
|
---|
9 | // History (add to end):
|
---|
10 | // Sophie Oct, 1999 - creation
|
---|
11 | //
|
---|
12 | //------------------------------------------------------------------------
|
---|
13 |
|
---|
14 | //---------------
|
---|
15 | // C++ Headers --
|
---|
16 | //---------------
|
---|
17 | #include "machdefs.h"
|
---|
18 | #include <iostream.h>
|
---|
19 | #include <typeinfo>
|
---|
20 | #include <math.h>
|
---|
21 |
|
---|
22 | #include "radspec.h"
|
---|
23 | #include "integ.h"
|
---|
24 |
|
---|
25 | /*!
|
---|
26 | * \class SOPHYA::RadSpectra <BR>
|
---|
27 | * This class is an abstract base class for radiation emission spectra. The flux() function returns the value of the flux (the spectral <BR>
|
---|
28 | * energy distribution) as a function of the frequency. As in the SpectralResponse class, the () operator has been redefined <BR>
|
---|
29 | * at this level, so that the user can access the flux value, either by calling the function or directly by using this operator. <BR>
|
---|
30 | * For all the sub-classes, \nu is given in units of Hz and
|
---|
31 | * the flux is returned in units of W/m^2/sr/Hz.
|
---|
32 | */
|
---|
33 |
|
---|
34 |
|
---|
35 | //----------------
|
---|
36 | // Constructor --
|
---|
37 | //----------------
|
---|
38 | /*! Default constructor */
|
---|
39 | /*!
|
---|
40 | The constructor takes as an argument the minimum
|
---|
41 | and the maximum frequency of the spectrum, if any. <BR>
|
---|
42 | In the case the user does not want to specify these
|
---|
43 | values, there are set respectively to 0. and 9.E49
|
---|
44 | by default.
|
---|
45 | */
|
---|
46 | RadSpectra::RadSpectra(double numin, double numax)
|
---|
47 | {
|
---|
48 | _numin = numin;
|
---|
49 | _numax = numax;
|
---|
50 | }
|
---|
51 |
|
---|
52 |
|
---|
53 | //--------------
|
---|
54 | // Destructor --
|
---|
55 | //--------------
|
---|
56 | RadSpectra::~RadSpectra()
|
---|
57 | {
|
---|
58 | }
|
---|
59 |
|
---|
60 | // ---------------------------
|
---|
61 | // -- Function Definitions --
|
---|
62 | // ---------------------------
|
---|
63 |
|
---|
64 | double
|
---|
65 | RadSpectra::minFreq() const
|
---|
66 | {
|
---|
67 | return _numin;
|
---|
68 | }
|
---|
69 |
|
---|
70 | double
|
---|
71 | RadSpectra::maxFreq() const
|
---|
72 | {
|
---|
73 | return _numax;
|
---|
74 | }
|
---|
75 |
|
---|
76 | double
|
---|
77 | RadSpectra::meanFreq() const
|
---|
78 | {
|
---|
79 | double result = (_numax+_numin)/2.;
|
---|
80 | return result;
|
---|
81 | }
|
---|
82 |
|
---|
83 |
|
---|
84 | /*
|
---|
85 | The peakFreq() function returns the value of the
|
---|
86 | frequency for the maximum value of the flux
|
---|
87 | */
|
---|
88 | double
|
---|
89 | RadSpectra::peakFreq() const
|
---|
90 | {
|
---|
91 | double maxAnswer = -1.e99;
|
---|
92 | double maxNu = -10;
|
---|
93 | double nu;
|
---|
94 | for (int i=1; i<1000;i++)
|
---|
95 | {
|
---|
96 | nu=(_numax-_numin)*i/1000.+_numin;
|
---|
97 | double lookForMax = flux(nu);
|
---|
98 | if(maxAnswer <= lookForMax) {
|
---|
99 | maxAnswer= lookForMax;
|
---|
100 | maxNu = nu;
|
---|
101 | }
|
---|
102 | }
|
---|
103 | return maxNu;
|
---|
104 | }
|
---|
105 |
|
---|
106 | void
|
---|
107 | RadSpectra::setMinMaxFreq(double numin, double numax)
|
---|
108 | {
|
---|
109 | _numin = numin;
|
---|
110 | _numax = numax;
|
---|
111 | }
|
---|
112 |
|
---|
113 | // the RadSpectra_fluxFunction function is used to call TrpzInteg double(double)
|
---|
114 | // (integration over a range of frequencies)
|
---|
115 | static RadSpectra* _raypourfinteg = NULL;
|
---|
116 | static double RadSpectra_fluxFunction(double nu)
|
---|
117 | {
|
---|
118 | return(_raypourfinteg->flux(nu));
|
---|
119 | }
|
---|
120 | /*!
|
---|
121 | The integratedFlux() function performs the integration
|
---|
122 | of the flux function in a frequency range <BR> defined by
|
---|
123 | f1 and f2.
|
---|
124 | */
|
---|
125 | double
|
---|
126 | RadSpectra::integratedFlux(double f1, double f2) const
|
---|
127 | {
|
---|
128 | if(f1 < this->minFreq()) f1 = this->minFreq();
|
---|
129 | if(f2 > this->maxFreq()) f2 = this->maxFreq();
|
---|
130 | _raypourfinteg = const_cast<RadSpectra *>(this);
|
---|
131 | TrpzInteg I(RadSpectra_fluxFunction , f1, f2);
|
---|
132 | double val = (double)I;
|
---|
133 | _raypourfinteg = NULL; // On ne peut pas faire ca avant la destruction de I
|
---|
134 | return(val);
|
---|
135 | }
|
---|
136 |
|
---|
137 | /*!
|
---|
138 | Same than integratedFlux() over the frequency range
|
---|
139 | of definition of the flux function
|
---|
140 | */
|
---|
141 | double
|
---|
142 | RadSpectra::integratedFlux() const
|
---|
143 | {
|
---|
144 | return integratedFlux(this->minFreq(),this->maxFreq());
|
---|
145 | }
|
---|
146 |
|
---|
147 | // integration using the logarithm !!
|
---|
148 | // Carefull!! Base 10....
|
---|
149 | static RadSpectra* _rayIntLog = NULL;
|
---|
150 |
|
---|
151 | static double RadSpectra_logFluxFunction(double tau)
|
---|
152 | {
|
---|
153 | double value = _rayIntLog->flux(pow(10,tau))*pow(10,tau);
|
---|
154 | return(value);
|
---|
155 | }
|
---|
156 |
|
---|
157 | /*!
|
---|
158 | The logIntegratedFlux() function performs the integration
|
---|
159 | of the flux function in a frequency range <BR> defined by
|
---|
160 | f1 and f2. The integration is here performed
|
---|
161 | on the logarithm of the flux function.
|
---|
162 | */
|
---|
163 | double
|
---|
164 | RadSpectra::logIntegratedFlux(double f1, double f2) const
|
---|
165 | {
|
---|
166 | if(f1 < this->minFreq()) f1 = this->minFreq();
|
---|
167 | if(f2 > this->maxFreq()) f2 = this->maxFreq();
|
---|
168 |
|
---|
169 | double f1Log = log10(f1);
|
---|
170 | double f2Log = log10(f2);
|
---|
171 | if(f1Log < -1.e99) f1Log = -1.e99;
|
---|
172 | if(f2Log > 1.e99) f2Log = 1.e99;
|
---|
173 | _rayIntLog = const_cast<RadSpectra *>(this);
|
---|
174 | TrpzInteg I(RadSpectra_logFluxFunction,f1Log,f2Log);
|
---|
175 | double value = (double)I * log(10.);
|
---|
176 | _rayIntLog = NULL;
|
---|
177 | return(value);
|
---|
178 | }
|
---|
179 |
|
---|
180 | /*!
|
---|
181 | same than logIntegratedFlux over the frequency range
|
---|
182 | of definition of the flux function
|
---|
183 | */
|
---|
184 | double
|
---|
185 | RadSpectra::logIntegratedFlux() const
|
---|
186 | {
|
---|
187 | return logIntegratedFlux(_numin,_numax);
|
---|
188 | }
|
---|
189 |
|
---|
190 | // the RadSpectra_filteredFlux function is used to call TrpzInteg double(double)
|
---|
191 | // (integration over a range of frequencies with a filter)
|
---|
192 | static SpectralResponse* _filter = NULL ;
|
---|
193 | static double RadSpectra_filteredFlux(double nu)
|
---|
194 | {
|
---|
195 | double flux = _raypourfinteg->flux(nu);
|
---|
196 | return(flux * _filter->transmission(nu));
|
---|
197 | }
|
---|
198 |
|
---|
199 | /*!
|
---|
200 | The filteredIntegratedFlux() function performs the integration
|
---|
201 | of the flux function in a frequency range <BR> defined by
|
---|
202 | f1 and f2 convolved by a SpectralResponse filter.
|
---|
203 | */
|
---|
204 | double
|
---|
205 | RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter, double f1, double f2) const
|
---|
206 | {
|
---|
207 | _raypourfinteg = const_cast<RadSpectra *>(this);
|
---|
208 | _filter = const_cast<SpectralResponse *>(&filter);
|
---|
209 | if(f1 < this->minFreq()) f1 = this->minFreq();
|
---|
210 | if(f2 > this->maxFreq()) f2 = this->maxFreq();
|
---|
211 |
|
---|
212 | TrpzInteg I(RadSpectra_filteredFlux,f1,f2);
|
---|
213 | double val = (double)I;
|
---|
214 | _raypourfinteg = NULL;
|
---|
215 | _filter = NULL;
|
---|
216 | return(val);
|
---|
217 | }
|
---|
218 |
|
---|
219 | /*!
|
---|
220 | Same than filteredIntegratedFlux() over the frequency range
|
---|
221 | defined as: <BR>
|
---|
222 | min_freq = MAX(minfreq_flux, minfreq_filter), <BR>
|
---|
223 | max_freq = MIN(maxfreq_flux, maxfreq_filter), <BR>
|
---|
224 | where:
|
---|
225 | <UL>
|
---|
226 | <LI> minfreq_flux is the minimum frequency of the flux definition
|
---|
227 | <LI> maxfreq_flux is the maximum frequency of the flux definition
|
---|
228 | <LI> minfreq_filter is the minimum frequency of the filter definition
|
---|
229 | <LI> maxfreq_filter is the maximum frequency of the filter definition
|
---|
230 | </UL>
|
---|
231 | */
|
---|
232 | double
|
---|
233 | RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter)
|
---|
234 | {
|
---|
235 | double minOfMin = filter.minFreq();
|
---|
236 | double maxOfMax = filter.maxFreq();
|
---|
237 | if(minOfMin < this->minFreq()) minOfMin = this->minFreq();
|
---|
238 | if(maxOfMax > this->maxFreq()) maxOfMax = this->maxFreq();
|
---|
239 | return(filteredIntegratedFlux(filter, minOfMin, maxOfMax ) );
|
---|
240 | }
|
---|
241 |
|
---|
242 |
|
---|
243 | // the RadSpectraVec_filteredFlux function is used to call TrpzInteg double(double)
|
---|
244 | // (integration over a range of frequencies with a filter)
|
---|
245 | static double RadSpectra_logFilteredFlux(double tau)
|
---|
246 | {
|
---|
247 | double nu = pow(10,tau);
|
---|
248 | double flux = _raypourfinteg->flux(nu)*nu;
|
---|
249 | double result = flux * _filter->transmission(nu);
|
---|
250 | return(result);
|
---|
251 | }
|
---|
252 |
|
---|
253 |
|
---|
254 | /*!
|
---|
255 | * The filteredIntegratedFlux() function performs the integration
|
---|
256 | * of the flux function in a frequency range <BR> defined by
|
---|
257 | * f1 and f2 convolved by a SpectralResponse filter (using the
|
---|
258 | * logarithm of the function).
|
---|
259 | */
|
---|
260 | double
|
---|
261 | RadSpectra::filteredLogIntFlux(SpectralResponse const& filter, double f1, double f2) const
|
---|
262 | {
|
---|
263 |
|
---|
264 | _raypourfinteg = NULL;
|
---|
265 | _filter = NULL;
|
---|
266 | if(f1 < this->minFreq()) f1 = this->minFreq();
|
---|
267 | if(f2 > this->maxFreq()) f2 = this->maxFreq();
|
---|
268 |
|
---|
269 | double f1Log = log10(f1);
|
---|
270 | double f2Log = log10(f2);
|
---|
271 | if(f1Log < -1.e99) f1Log = -1.e99;
|
---|
272 | if(f2Log > 1.e99) f2Log = 1.e99;
|
---|
273 | _raypourfinteg = const_cast<RadSpectra *>(this);
|
---|
274 | _filter = const_cast<SpectralResponse *>(&filter);
|
---|
275 | TrpzInteg I(RadSpectra_logFilteredFlux,f1Log,f2Log);
|
---|
276 | double val = (double)I;
|
---|
277 | _raypourfinteg = NULL;
|
---|
278 | _filter = NULL;
|
---|
279 | return(val* log(10.));
|
---|
280 | }
|
---|
281 |
|
---|
282 | double
|
---|
283 | RadSpectra::filteredLogIntFlux(SpectralResponse const& filter)
|
---|
284 | {
|
---|
285 | return(filteredLogIntFlux(filter, filter.minFreq(), filter.maxFreq() ) );
|
---|
286 | }
|
---|
287 |
|
---|
288 |
|
---|
289 |
|
---|
290 |
|
---|
291 | void
|
---|
292 | RadSpectra::Print(ostream& os) const
|
---|
293 | {
|
---|
294 | // os << "RadSpectra::Print (" << typeid(*this).name()
|
---|
295 | // << ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
|
---|
296 | os << "RadSpectra::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
|
---|
297 | os << "MeanFreq= " << meanFreq() << " Emission= " << flux(meanFreq()) << endl;
|
---|
298 | os << "PeakFreq= " << peakFreq() << " Emission= " << flux(peakFreq()) << endl;
|
---|
299 |
|
---|
300 | }
|
---|
301 |
|
---|
302 |
|
---|