From d21af1b445edf644f001819d57265507e26bd226 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Tue, 6 Oct 2020 08:29:14 +0200 Subject: [PATCH 1/2] Add function to help Newton-Raphson convergence Signed-off-by: Adam.Dybbroe --- pyspectral/blackbody.py | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/pyspectral/blackbody.py b/pyspectral/blackbody.py index efbf0b10..aa473c58 100644 --- a/pyspectral/blackbody.py +++ b/pyspectral/blackbody.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2013-2019 Adam.Dybbroe +# Copyright (c) 2013-2020 Adam.Dybbroe # Author(s): @@ -140,7 +140,7 @@ def planck(wave, temperature, wavelength=True): Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1 Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1: - 1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1 + 1.0 W/m^2 sr^-1 (m^-1)^-1 = 1.0e5 mW/m^2 sr^-1 (cm^-1)^-1 """ units = ['wavelengths', 'wavenumbers'] @@ -230,7 +230,7 @@ def blackbody_wn(wavenumber, temp): Unit = W/m^2 sr^-1 (m^-1)^-1 = W/m sr^-1 Converting from SI units to mW/m^2 sr^-1 (cm^-1)^-1: - 1.0 W/m^2 sr^-1 (m^-1)^-1 = 0.1 mW/m^2 sr^-1 (cm^-1)^-1 + 1.0 W/m^2 sr^-1 (m^-1)^-1 = 1.0e5 mW/m^2 sr^-1 (cm^-1)^-1 """ return planck(wavenumber, temp, wavelength=False) @@ -249,3 +249,20 @@ def blackbody(wavel, temp): """ return planck(wavel, temp, wavelength=True) + + +def ratio_planck_vs_planck_first_derivative(wavel, temp): + """Derive the ratio of the Planck radiation over the first derivative (in T) of the Planck radiation. + + SI units. + + wavel = Wavelength or a sequence of wavelengths (m) + temp = Temperature (scalar) or a sequence of temperatures (K) + + Output: Temperature in Kelvin + Unit = K + + """ + + const = H_PLANCK * C_SPEED / (K_BOLTZMANN * wavel) + return temp**2 / const * (1. - 1./np.exp(const/temp)) From dc2af95dae06d7e1ce14681cf66b045079682577 Mon Sep 17 00:00:00 2001 From: "Adam.Dybbroe" Date: Tue, 6 Oct 2020 08:29:53 +0200 Subject: [PATCH 2/2] Add new function for radiance2tb conversion Signed-off-by: Adam.Dybbroe --- pyspectral/radiance_tb_conversion.py | 46 ++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/pyspectral/radiance_tb_conversion.py b/pyspectral/radiance_tb_conversion.py index 1be5daad..9c7b2505 100644 --- a/pyspectral/radiance_tb_conversion.py +++ b/pyspectral/radiance_tb_conversion.py @@ -31,6 +31,7 @@ import numpy as np from pyspectral.blackbody import (H_PLANCK, K_BOLTZMANN, C_SPEED) from pyspectral.blackbody import blackbody, blackbody_wn +from pyspectral.blackbody import ratio_planck_vs_planck_first_derivative from pyspectral.utils import WAVE_NUMBER from pyspectral.utils import WAVE_LENGTH from pyspectral.utils import BANDNAMES @@ -257,6 +258,51 @@ def radiance2tb(self, rad): """ return radiance2tb(rad, self.rsr[self.bandname][self.detector]['central_wavelength'] * 1e-6) + def band_radiance2tb(self, rad): + """Get the Tb given the observed (band integrated) radiance and the spectral response function. + + Using an iterative approach: + + 1) set uncertainty parameter delta-L + + 2) set T_j = T_{first guess} + + 3) Calculate B(T_j) + + 4) if (B(T_j) - L_{measure}) > delta-L then adjust T_j and go back to 3) + + 5) T_j matches the measurement within the defined uncertainty + + rad: + Radiance in SI units + """ + if self.wavespace == WAVE_NUMBER: + raise NotImplementedError( + 'Getting the band Tb from the band radiance is not supported in wavenumber space yet!') + + delta_rad = 1000.0 + first_guess_temp = self.radiance2tb(rad) + temp_i = first_guess_temp + rad_dev = 10000.0 + + import ipdb + while(rad_dev > delta_rad): + + ratio_ = ratio_planck_vs_planck_first_derivative( + self.wavelength_or_wavenumber, temp_i) * self.response + ratio_ = integrate.trapz(ratio_, self.wavelength_or_wavenumber) / self.rsr_integral + + ipdb.set_trace() + + temp_i = temp_i - ratio_ + + #ratio_ = ratio_planck_vs_planck_first_derivative(self.wavelength_or_wavenumber, temp_i) + + temp_i = temp_i - ratio_ + rad_dev = abs(rad - self.tb2radiance(temp_i)['radiance']) + + return temp_i + def radiance2tb(rad, wavelength): """Get the Tb from the radiance using the Planck function.