source: Sophya/trunk/Cosmo/SimLSS/planckspectra.cc@ 3617

Last change on this file since 3617 was 3347, checked in by cmv, 18 years ago
  • changements des noms des methodes.
  • definition des options par enum. cmv 11/10/2007
File size: 11.3 KB
Line 
1#include "machdefs.h"
2#include <iostream>
3#include <stdlib.h>
4#include <stdio.h>
5#include <string.h>
6#include <math.h>
7
8#include "pexceptions.h"
9
10#include "constcosmo.h"
11#include "planckspectra.h"
12
13namespace SOPHYA {
14
15/*
16----------------------------------------------------------------
17- UNITES:
18 Temperature en Kelvin "K"
19 Frequence en GigaHertz "Hz"
20 Longueur d'onde en Metre "m"
21----------------------------------------------------------------
22- Frequence / Longeur d'onde
23 l = c/f ; f = c/l
24 dl = c/f^2 df = l^2/c df
25 df = c/l^2 dl = f^2/c dl
26 P[f] df = P[l] dl = P[f] c/l^2 dl = P[l] c/f^2 df
27----------------------------------------------------------------
28- Lois de Planck Wien et Rayleigh
29 Planck H*Nu/(K*T) = H*C/(Lambda*K*T) quelconque
30 Wien H*Nu/(K*T) = H*C/(Lambda*K*T) >> 1
31 Rayleigh H*Nu/(K*T) = H*C/(Lambda*K*T) << 1
32----------------------------------------------------------------
33- Les termes en facteur pour les lois de Planck Wien et Rayleigh
34 (la partie contenant l'exponentielle)
35 ainsi que pour les derivees des ces termes par rapport a la temperature
36 Planck: 1/[exp(H*f/(K*t))-1] = 1/[exp(HC/(l*K*t))-1]
37 der= H*f/(K*t) /t * exp(HC/(l*K*t))/[exp(HC/(l*K*t))-1]^2
38 Wien: 1/exp(H*f/(K*t))
39 der= H*f/(K*t) /t /exp(H*f/(K*t))
40 Rayleigh: 1/(H*f/(K*t))
41 der= K/(H*f)
42 Pour x->0: compute exp(x)=1+x+x^2/2 for 1/(exp(x)-1)
43 soit 1/(exp(x)-1) -> 1/(x+x^2/2)
44 car exp(x)-1 vaut zero plus vite que son DL a cause de la soustraction
45----------------------------------------------------------------
46 - RADIANCE (Brillance) en fonction de nu ou lambda
47 nombre de watts (ou photons par seconde)
48 radies par unite de surface par steradian et par hertz (ou metre)
49 2*f^2/c^2 * 1/(exp()-1) ph/s/m^2/Sr/Hz
50 2*h*f^3/c^2 * 1/(exp()-1) W/m^2/Sr/Hz
51 2*c/l^4 * 1/(exp()-1) ph/s/m^2/Sr/m
52 2*h*c^2/l^5 * 1/(exp()-1) W/m^2/Sr/m
53----------------------------------------------------------------
54 - EMITTANCE en fonction de nu ou lambda
55 nombre de watts (ou photons par seconde)
56 radies dans un hemisphere (2*PI sr)
57 par unite de surface et par hertz (ou metre)
58 ..quand on integre sur 2Pi sr il faut
59 ..multiplier par cos(theta) -> facteur Pi uniquement
60 EMITTANCE = Pi * RADIANCE
61 2*Pi*f^2/c^2 * 1/(exp()-1) ph/s/m^2/Hz
62 2*Pi*h*f^3/c^2 * 1/(exp()-1) W/m^2/Hz
63 2*Pi*c/l^4 * 1/(exp()-1) ph/s/m^2/m
64 2*Pi*h*c^2/l^5 * 1/(exp()-1) W/m^2/m
65----------------------------------------------------------------
66 - DENSITE en fonction de nu ou lambda
67 nombre de joules (ou photons) par m^3 et par hertz (ou metre)
68 DENSITE = 4*Pi/c * RADIANCE
69 8*Pi*f^2/c^3 * 1/(exp()-1) ph/m^3/Hz
70 8*Pi*h*f^3/c^3 * 1/(exp()-1) J/m^3/Hz
71 8*Pi/l^4 * 1/(exp()-1) ph/m^3/m
72 8*Pi*h*c/l^5 * 1/(exp()-1) J/m^3/m
73----------------------------------------------------------------
74 - DENSITE d'energie totale (loi de Stefan-Boltzman)
75 4*sigma*T^4/C J/m^3
76 Sigma = 2*PI^5*K^4/(15.*C^2*H^3) = PI^2*K^4/(60.*Hbar^3*C2)
77----------------------------------------------------------------
78*/
79
80/*
81POUR TESTER LES LIMITES DE L'EXPONENTIELLE:
821-/ copier/coller les lignes suivantes dans toto.icc
83double x=10.;
84for(int i=0;i<35;i++) {
85 double ex = exp(x);
86 double exm1a = x*(1.+x/2.*(1.+x/3.)),;
87 printf("x=%g exp(x)=%e exp(x)-1=%e dl=%e (d=%e)\n",x,ex,ex-1.,exm1a,ex-1.-exm1a);
88 x /= 10.;
89}
902-/ executer sous spiapp: > c++execfrf toto.icc
91*/
92
93
94#define MAXARGEXP 100. /* maxi = 709.7827128933840 */
95#define MINARGEXP 1.e-10 /* Replace exp(x) by DL for 1/(exp(x)-1) when x->0 */
96
97//************************************************************************
98// vitesse de la lumiere en m/s
99double CinM_ = SpeedOfLight_Cst*1e+3;
100// h/k pour une frequence en Hz
101double HsK_ = h_Planck_Cst/k_Boltzman_Cst;
102// hc/k pour une longeur d'onde en metre
103double HCsK_ = h_Planck_Cst*CinM_/k_Boltzman_Cst;
104
105//************************************************************************
106PlanckSpectra::PlanckSpectra(double T)
107 : T_(T) , spectraapprox_(PLANCK) , spectrafunc_(VALUE) , spectravar_(NU)
108 , spectraunit_(ANGSFLUX) , spectrapower_(POWER)
109{
110}
111
112PlanckSpectra::PlanckSpectra(PlanckSpectra& s)
113 : T_(s.T_) , spectraapprox_(s.spectraapprox_) , spectrafunc_(s.spectrafunc_) , spectravar_(s.spectravar_)
114 , spectraunit_(s.spectraunit_), spectrapower_(s.spectrapower_)
115{
116}
117
118PlanckSpectra::~PlanckSpectra(void)
119{
120}
121
122//************************************************************************
123void PlanckSpectra::SetSpectraApprox(SpectraApprox spectraapprox)
124// spectraapprox = PLANCK : Planck spectra
125// = RAYLEIGH : Rayleigh spectra
126// = WIEN : Wien spectra
127{
128 spectraapprox_ = spectraapprox;
129}
130
131void PlanckSpectra::SetSpectraFunc(SpectraFunc spectrafunc)
132// spectrafunc = VALUE : spectra
133// = DERIV : derivative of spectra relative to temperature
134{
135 spectrafunc_ = spectrafunc;
136}
137
138void PlanckSpectra::SetSpectraVar(SpectraVar spectravar)
139// spectravar = NU : variable is frequency (Hz)
140// = LAMBDA : variable is wavelength (m)
141{
142 spectravar_ = spectravar;
143}
144
145void PlanckSpectra::SetSpectraPower(SpectraPower spectrapower)
146// spectrapower = POWER : in W/...
147// = PHOTON : in photon/s/...
148{
149 spectrapower_ = spectrapower;
150}
151
152void PlanckSpectra::SetSpectraUnit(SpectraUnit spectraunit)
153// spectraunit = ANGSFLUX : radiance W/m^2/Sr/Hz or ...
154// = SFLUX : emittance W/m^2/Hz or ...
155// = DENSENERG : density J/m^3/Hz or ...
156// = EXPON : just the function exp() or equivalent for Rayleigh
157{
158 spectraunit_ = spectraunit;
159}
160
161
162//************************************************************************
163double PlanckSpectra::F2L(double nu_or_lambda)
164// Convert frequency (GHz) to/from Wavelength (m)
165// Lambda(m) * Nu(GHz) = C(km/s) * 10^-6
166{
167 return CinM_ / nu_or_lambda;
168}
169
170//************************************************************************
171double PlanckSpectra::PlanckExp(double fl)
172{
173 double x = (spectravar_==LAMBDA) ? HCsK_/(fl*T_) : HsK_*fl/T_;
174
175 if(spectraapprox_ == RAYLEIGH) return 1./x; // Rayleigh
176
177 if(x>MAXARGEXP) return 0.;
178
179 if(spectraapprox_ == WIEN) return exp(-x); // Wien
180
181 if(x<MINARGEXP) return 1./(x*(1.+x/2.));
182 return 1./(exp(x)-1.);
183}
184
185double PlanckSpectra::DPlanckExp_DT(double fl)
186{
187 if(spectraapprox_ == RAYLEIGH)
188 return (spectravar_==LAMBDA) ? fl/HCsK_: 1. / (HsK_*fl); // Rayleigh
189
190 double x = (spectravar_==LAMBDA) ? HCsK_/(fl*T_) : HsK_*fl/T_;
191
192 if(x>MAXARGEXP) return 0.;
193
194 if(spectraapprox_ == WIEN) return x/T_ * exp(-x); // Wien
195
196 double v;
197 if(x<MINARGEXP) { // 1/(exp(x)-1) -derive-> -exp(x)/(exp(x)-1)^2
198 v = x*(1.+x/2.);
199 return x/T_ * (1.+v)/(v*v);
200 }
201 v = exp(x);
202 return x/T_ * v/((v-1.)*(v-1.));
203}
204
205//------------------------------------------------------------
206double PlanckSpectra::operator() (double x)
207{
208 // --- type de spectre et d'approximation selon choix variable
209 double fex = (spectrafunc_==DERIV) ? DPlanckExp_DT(x): PlanckExp(x);
210
211 // --- le coefficient multiplicateur devant l'expo
212 if(spectraunit_ == EXPON) return fex; // expo seule
213
214 // radiance en ph/... longeur d'onde frequency
215 double coeff = (spectravar_==LAMBDA) ? 2.*CinM_/(x*x*x*x): 2.*x*x/(CinM_*CinM_);
216
217 if(spectraunit_ == SFLUX) coeff *= M_PI; // emittance en ph/...
218 else if(spectraunit_ == DENSENERG) coeff *= 4.*M_PI/CinM_; // densite en ph/...
219
220 // unite en W/... longeur d'onde frequency
221 if(spectrapower_ == POWER)
222 coeff *= (spectravar_==LAMBDA) ? h_Planck_Cst*CinM_/x: h_Planck_Cst*x;
223
224 return coeff * fex;
225}
226
227//------------------------------------------------------------
228double PlanckSpectra::PlanckEnergie(void)
229// densite d'energie: energie par m^3 en "J / m^3"
230// cste de Stefan-Boltzmann en W/m^2/K^4
231{
232 return 4.*StefBoltz_Cst*T_*T_*T_*T_/(SpeedOfLight_Cst*1000.);
233}
234
235//------------------------------------------------------------
236double PlanckSpectra::PlanckPhoton(void)
237// densite de photons: nombre de photons par m^3 en "ph / m^3"
238{
239 double x = k_Boltzman_Cst*T_/(h_Planck_Cst*SpeedOfLight_Cst*1000.);
240 return 16.*M_PI*Zeta_3_Cst* x*x*x;
241}
242
243double PlanckSpectra::WienLaw(void)
244// Retourne la frequence (ou longeur d'onde) du maximum
245// du spectre de radiance (ou emittance ou densite).
246// La valeur du maximum depend de l'unite du spectre (W/... ou Ph/s/...)
247/*
248--- Les fonctions:
249En frequence: n = 3 en W/... , n = 2 en Ph/...
250En longueur d'onde: n = 5 en W/... , n = 4 en Ph/...
251f[x_] := x^n*Exp[-x]; Solve[f'[x]==0,x] // Wien
252f[x_] := x^n/(Exp[x]-1); NSolve[f'[x]==0,x,17] // Planck
253La derivee par rapport a T des fonctions en frequence:
254df/dT[x_] := x^(n+1)*Exp[-x]; Solve[f'[x]==0,x] // Wien
255df/dT[x_] := x^(n+1)*Exp[x]/(Exp[x]-1)^2; NSolve[f'[x]==0,x,17] // Planck
256*/
257{
258 double x=-1.;
259 if(spectravar_ == NU) { // Frequence
260 if(spectrafunc_ == VALUE) { // Le spectre
261 if(spectrapower_ == POWER) // W/...
262 {if(spectraapprox_==PLANCK) x=2.82143938074; else if(spectraapprox_==WIEN) x=3.;}
263 else if(spectrapower_ == PHOTON) // Ph/s/...
264 {if(spectraapprox_==PLANCK) x=1.59362426; else if(spectraapprox_==WIEN) x=2.;}
265 } else if(spectrafunc_ == DERIV) { // La derivee du spectre
266 if(spectrapower_ == POWER) // W/...
267 {if(spectraapprox_==PLANCK) x=3.830016095; else if(spectraapprox_==WIEN) x=4.;}
268 else if(spectrapower_ == PHOTON) // Ph/s/...
269 {if(spectraapprox_==PLANCK) x=2.575678888; else if(spectraapprox_==WIEN) x=3.;}
270 }
271 } else { // Longeur d'onde
272 if(spectrafunc_ == VALUE) { // Le spectre
273 if(spectrapower_ == POWER) // W/...
274 {if(spectraapprox_==PLANCK) x=4.96511420; else if(spectraapprox_==WIEN) x=5.;}
275 else if(spectrapower_ == PHOTON) // Ph/s/...
276 {if(spectraapprox_==PLANCK) x=3.9206904099; else if(spectraapprox_==WIEN) x=4.;}
277 } else if(spectrafunc_ == DERIV) { // La derivee du spectre
278 if(spectrapower_ == POWER) // W/...
279 {if(spectraapprox_==PLANCK) x=5.969409190; else if(spectraapprox_==WIEN) x=6.;}
280 else if(spectrapower_ == PHOTON) // Ph/s/...
281 {if(spectraapprox_==PLANCK) x=4.928119345; else if(spectraapprox_==WIEN) x=5.;}
282 }
283 }
284
285 if(x<0.) return x; // Rayleigh n'a pas de maximum
286 if(spectravar_ == NU) return x / HsK_ *T_; // Frequence
287 return HCsK_ / x /T_; // Longeur d'onde
288}
289
290double PlanckSpectra::FindMaximum(double eps)
291// Retourne le maximum des spectres (a la precicion eps)
292{
293 if(spectraapprox_ == RAYLEIGH) { // Pas de maximum pour les spectres de Raleigh
294 cout<<"PlanckSpectra::FindMaximum : no maximum for Raleigh spectra"<<endl;
295 throw ParmError("PlanckSpectra::FindMaximum : no maximum for Raleigh spectra");
296 }
297 if(eps<=0.) eps = 0.001;
298
299 double xmax = WienLaw(), wmax=(*this)(xmax);
300 double xlim[2]={xmax,xmax};
301
302 // Encadrement du maximum par xlim[0] et xlim[1]
303 while( (*this)(xlim[1])>=0.01*wmax ) xlim[1] *= 2.;
304 while( (*this)(xlim[0])>=0.01*wmax ) xlim[0] /= 2.;
305 //cout<<"FindMaximum: wmax="<<wmax<<" pour x=["<<xlim[0]<<","<<xlim[1]<<"]"<<endl;
306
307 while( (xlim[1]-xlim[0])/xmax > eps ) {
308 wmax=-1.;
309 double dx = (xlim[1]-xlim[0])/10.;
310 for(double x=xlim[0]; x<=xlim[1]; x+=dx) {
311 double f=(*this)(x);
312 if(f>wmax) {wmax=f; xmax=x;}
313 }
314 //cout<<"...wmax="<<wmax<<" pour xmax="<<xmax<<" dx="<<dx<<endl;
315 xlim[0] = xmax-dx*1.1; xlim[1] = xmax+dx*1.1;
316 }
317
318
319 //cout<<"wmax="<<wmax<<" pour x="<<xmax<<" dx="<<xlim[1]-xlim[0]<<endl;
320 return xmax;
321}
322
323} // Fin namespace SOPHYA
Note: See TracBrowser for help on using the repository browser.