Changeset 3347 in Sophya for trunk/Cosmo/SimLSS/planckspectra.cc
- Timestamp:
- Oct 11, 2007, 4:34:42 PM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cosmo/SimLSS/planckspectra.cc
r3325 r3347 105 105 //************************************************************************ 106 106 PlanckSpectra::PlanckSpectra(double T) 107 : T_(T) , approx_(0) , deriv_(0) , typvar_(0) , typspec_(0) , unitout_(0) 107 : T_(T) , spectraapprox_(PLANCK) , spectrafunc_(VALUE) , spectravar_(NU) 108 , spectraunit_(ANGSFLUX) , spectrapower_(POWER) 108 109 { 109 110 } 110 111 111 112 PlanckSpectra::PlanckSpectra(PlanckSpectra& s) 112 : T_(s.T_) , approx_(s.approx_) , deriv_(s.deriv_) , typvar_(s.typvar_) , typspec_(s.typspec_)113 , unitout_(s.unitout_)113 : T_(s.T_) , spectraapprox_(s.spectraapprox_) , spectrafunc_(s.spectrafunc_) , spectravar_(s.spectravar_) 114 , spectraunit_(s.spectraunit_), spectrapower_(s.spectrapower_) 114 115 { 115 116 } … … 120 121 121 122 //************************************************************************ 122 void PlanckSpectra::SetApprox(unsigned short approx) 123 // approx = 0 : Planck spectra 124 // = 1 : Rayleigh spectra 125 // = 2 : Wien spectra 126 { 127 if(approx>2) { 128 cout<<"PlanckSpectra::SetApprox : unknown approximation"<<endl; 129 throw ParmError("PlanckSpectra::SetApprox : unknown approximation"); 130 } 131 approx_ = approx; 132 } 133 134 void PlanckSpectra::SetDeriv(unsigned short deriv) 135 // deriv = 0 : spectra 136 // = 1 : derivative of spectra relative to temperature 137 { 138 if(deriv>1) { 139 cout<<"PlanckSpectra::SetDeriv : parameter out of range"<<endl; 140 throw ParmError("PlanckSpectra::SetDeriv : parameter out of range"); 141 } 142 deriv_ = deriv; 143 } 144 145 void PlanckSpectra::SetVar(unsigned short typvar) 146 // typvar = 0 : variable is frequency (Hz) 147 // = 1 : variable is wavelength (m) 148 { 149 if(typvar>1) { 150 cout<<"PlanckSpectra::SetVar : parameter out of range"<<endl; 151 throw ParmError("PlanckSpectra::SetVar : parameter out of range"); 152 } 153 typvar_ = typvar; 154 } 155 156 void PlanckSpectra::SetUnitOut(unsigned short unitout) 157 // unitout = 0 : in W/... 158 // = 1 : in photon/s/... 159 { 160 if(unitout>1) { 161 cout<<"PlanckSpectra::SetUnitOut : parameter out of range"<<endl; 162 throw ParmError("PlanckSpectra::SetUnitOut : parameter out of range"); 163 } 164 unitout_ = unitout; 165 } 166 167 void PlanckSpectra::SetTypSpectra(unsigned short typspec) 168 // typspec = 0 : radiance W/m^2/Sr/Hz or ... 169 // = 1 : emittance W/m^2/Hz or ... 170 // = 2 : density J/m^3/Hz or ... 171 // = 3 : just the function exp() or equivalent for Rayleigh 172 { 173 if(typspec>3) { 174 cout<<"PlanckSpectra::SetTypSpectra : parameter out of range"<<endl; 175 throw ParmError("PlanckSpectra::SetTypSpectra : parameter out of range"); 176 } 177 typspec_ = typspec; 123 void PlanckSpectra::SetSpectraApprox(SpectraApprox spectraapprox) 124 // spectraapprox = PLANCK : Planck spectra 125 // = RAYLEIGH : Rayleigh spectra 126 // = WIEN : Wien spectra 127 { 128 spectraapprox_ = spectraapprox; 129 } 130 131 void PlanckSpectra::SetSpectraFunc(SpectraFunc spectrafunc) 132 // spectrafunc = VALUE : spectra 133 // = DERIV : derivative of spectra relative to temperature 134 { 135 spectrafunc_ = spectrafunc; 136 } 137 138 void PlanckSpectra::SetSpectraVar(SpectraVar spectravar) 139 // spectravar = NU : variable is frequency (Hz) 140 // = LAMBDA : variable is wavelength (m) 141 { 142 spectravar_ = spectravar; 143 } 144 145 void PlanckSpectra::SetSpectraPower(SpectraPower spectrapower) 146 // spectrapower = POWER : in W/... 147 // = PHOTON : in photon/s/... 148 { 149 spectrapower_ = spectrapower; 150 } 151 152 void 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; 178 159 } 179 160 … … 190 171 double PlanckSpectra::PlanckExp(double fl) 191 172 { 192 double x = ( typvar_) ? HCsK_/(fl*T_) : HsK_*fl/T_;193 194 if( approx_ == 1) return 1./x; // Rayleigh173 double x = (spectravar_==LAMBDA) ? HCsK_/(fl*T_) : HsK_*fl/T_; 174 175 if(spectraapprox_ == RAYLEIGH) return 1./x; // Rayleigh 195 176 196 177 if(x>MAXARGEXP) return 0.; 197 178 198 if( approx_ == 2) return exp(-x); // Wien179 if(spectraapprox_ == WIEN) return exp(-x); // Wien 199 180 200 181 if(x<MINARGEXP) return 1./(x*(1.+x/2.)); … … 204 185 double PlanckSpectra::DPlanckExp_DT(double fl) 205 186 { 206 if(approx_ == 1) return (typvar_) ? fl/HCsK_: 1. / (HsK_*fl); // Rayleigh 207 208 double x = (typvar_) ? HCsK_/(fl*T_) : HsK_*fl/T_; 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_; 209 191 210 192 if(x>MAXARGEXP) return 0.; 211 193 212 if( approx_ == 2) return x/T_ * exp(-x); // Wien194 if(spectraapprox_ == WIEN) return x/T_ * exp(-x); // Wien 213 195 214 196 double v; … … 225 207 { 226 208 // --- type de spectre et d'approximation selon choix variable 227 double fex = ( deriv_) ? DPlanckExp_DT(x): PlanckExp(x);209 double fex = (spectrafunc_==DERIV) ? DPlanckExp_DT(x): PlanckExp(x); 228 210 229 211 // --- le coefficient multiplicateur devant l'expo 230 if(typspec_ == 3) return fex; // expo seule 231 232 // radiance en ph/... longeur d'onde frequency 233 double coeff = (typvar_) ? 2.*CinM_/(x*x*x*x): 2.*x*x/(CinM_*CinM_); 234 235 if(typspec_ == 1) coeff *= M_PI; // emittance en ph/... 236 else if(typspec_ == 2) coeff *= 4.*M_PI/CinM_; // densite en ph/... 237 238 // unite en W/... longeur d'onde frequency 239 if(unitout_ == 0) coeff *= (typvar_) ? h_Planck_Cst*CinM_/x: h_Planck_Cst*x; 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; 240 223 241 224 return coeff * fex; … … 274 257 { 275 258 double x=-1.; 276 if( typvar_ == 0) { // Frequence277 if( deriv_ == 0) { // Le spectre278 if( unitout_ == 0) // W/...279 {if( approx_==0) x=2.82143938074; else if(approx_==2) x=3.;}280 else if( unitout_ == 1) // Ph/s/...281 {if( approx_==0) x=1.59362426; else if(approx_==2) x=2.;}282 } else if( deriv_ == 1) { // La derivee du spectre283 if( unitout_ == 0) // W/...284 {if( approx_==0) x=3.830016095; else if(approx_==2) x=4.;}285 else if( unitout_ == 1) // Ph/s/...286 {if( approx_==0) x=2.575678888; else if(approx_==2) x=3.;}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.;} 287 270 } 288 271 } else { // Longeur d'onde 289 if( deriv_ == 0) { // Le spectre290 if( unitout_ == 0) // W/...291 {if( approx_==0) x=4.96511420; else if(approx_==2) x=5.;}292 else if( unitout_ == 1) // Ph/s/...293 {if( approx_==0) x=3.9206904099; else if(approx_==2) x=4.;}294 } else if( deriv_ == 1) { // La derivee du spectre295 if( unitout_ == 0) // W/...296 {if( approx_==0) x=5.969409190; else if(approx_==2) x=6.;}297 else if( unitout_ == 1) // Ph/s/...298 {if( approx_==0) x=4.928119345; else if(approx_==2) x=5.;}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.;} 299 282 } 300 283 } 301 284 302 285 if(x<0.) return x; // Rayleigh n'a pas de maximum 303 if( typvar_ == 0) return x / HsK_ *T_; // Frequence286 if(spectravar_ == NU) return x / HsK_ *T_; // Frequence 304 287 return HCsK_ / x /T_; // Longeur d'onde 305 288 } … … 308 291 // Retourne le maximum des spectres (a la precicion eps) 309 292 { 310 if( approx_ == 1) { // Pas de maximum pour les spectres de Raleigh293 if(spectraapprox_ == RAYLEIGH) { // Pas de maximum pour les spectres de Raleigh 311 294 cout<<"PlanckSpectra::FindMaximum : no maximum for Raleigh spectra"<<endl; 312 295 throw ParmError("PlanckSpectra::FindMaximum : no maximum for Raleigh spectra");
Note:
See TracChangeset
for help on using the changeset viewer.