diff --git a/Code/Radiometry/otbImageToReflectanceImageFilter.h b/Code/Radiometry/otbImageToReflectanceImageFilter.h index fbf467ec45434c8f0a0ead28f2a1dccd8fa205e1..a3151a4ca7d3ade75f59ea7967f024b0af7d6d92 100644 --- a/Code/Radiometry/otbImageToReflectanceImageFilter.h +++ b/Code/Radiometry/otbImageToReflectanceImageFilter.h @@ -22,6 +22,7 @@ #ifndef __otbImageToReflectanceImageFilter_h #define __otbImageToReflectanceImageFilter_h +#include "otbVarSol.h" #include "otbUnaryImageFunctorWithVectorImageFilter.h" #include "otbImageToLuminanceImageFilter.h" #include "otbLuminanceToReflectanceImageFilter.h" @@ -315,12 +316,8 @@ protected: { if (m_Day * m_Month != 0 && m_Day < 32 && m_Month < 13) { - otb_6s_doublereal dsol = 0.; - otb_6s_integer day = static_cast<otb_6s_integer>(m_Day); - otb_6s_integer month = static_cast<otb_6s_integer>(m_Month); - //int cr(0); - otb_6s_varsol_(&day, &month, &dsol); - coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * static_cast<double>(dsol); + double dsol = VarSol::GetVarSol(m_Day, m_Month); + coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * dsol; } else { diff --git a/Code/Radiometry/otbLuminanceToReflectanceImageFilter.h b/Code/Radiometry/otbLuminanceToReflectanceImageFilter.h index 65dd878ff482c7a5b3fd964acf3f108fcea5ea36..2d29e7d3f2377c8cf083ddfbe5f85909689c1752 100644 --- a/Code/Radiometry/otbLuminanceToReflectanceImageFilter.h +++ b/Code/Radiometry/otbLuminanceToReflectanceImageFilter.h @@ -22,7 +22,7 @@ #ifndef __otbLuminanceToReflectanceImageFilter_h #define __otbLuminanceToReflectanceImageFilter_h -#include "otb_6S.h" +#include "otbVarSol.h" #include "otbUnaryImageFunctorWithVectorImageFilter.h" #include "otbVectorImage.h" #include "otbMath.h" @@ -297,12 +297,8 @@ protected: { if (m_Day * m_Month != 0 && m_Day < 32 && m_Month < 13) { - otb_6s_doublereal dsol = 0.; - otb_6s_integer day = static_cast<otb_6s_integer>(m_Day); - otb_6s_integer month = static_cast<otb_6s_integer>(m_Month); - //int cr(0); - otb_6s_varsol_(&day, &month, &dsol); - coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * static_cast<double>(dsol); + double dsol = VarSol::GetVarSol(m_Day, m_Month); + coefTemp = vcl_cos(m_ZenithalSolarAngle * CONST_PI_180) * dsol; } else { diff --git a/Code/Radiometry/otbVarSol.h b/Code/Radiometry/otbVarSol.h new file mode 100644 index 0000000000000000000000000000000000000000..6f772884c6147e9f44df7fe10e31cd374a2a000d --- /dev/null +++ b/Code/Radiometry/otbVarSol.h @@ -0,0 +1,72 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __otbVarSol_h +#define __otbVarSol_h + +#include "otbMath.h" + +namespace otb +{ + +/** \class VarSol + * \brief VarSol operations. + * + * Call VarSol main function + * + * [1] Vermote, E., Tanré, D., Deuze, J., Herman, M., Morcette, J., 1997. + * Second simulation of the satellite signal in the solar spectrum, 6S: An overview. + * IEEE Transactions on Geoscience and Remote Sensing 35 + * \ingroup Radiometry + * + */ + class ITK_EXPORT VarSol + { + public: + + /** Call the varSol function*/ + static double GetVarSol(const int day, const int month) + { + /* System generated locals */ + double d__1; + + /* Builtin functions */ + //double acos(double), cos(double); + + /* Local variables */ + int j; + double pi, om; + + /* calculation of the variability of the solar constant during the year. */ + /* day is the number of the day in the month */ + if (month <= 2) + j = (month - 1) * 31 + day; + else if (month > 8) + j = (month - 1) * 31 - (month - 2) / 2 - 2 + day; + else + j = (month - 1) * 31 - (month - 1) / 2 - 2 + day; + + pi = acos(0.) * 2.; + om = (double) (j - 4) * .9856 * pi / 180.; + /* Computing 2nd power */ + d__1 = 1. - cos(om) * .01673; + return 1. / (d__1 * d__1); + } + }; +} //end namespace otb +#endif +