Feature tb from radiance iterative solution

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))

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.