/******************************************************************************
 *
 * Purpose:  Implementation of ReflectanceCalculator class. Calculate
 *           reflectance values from radiance, for visual bands.
 * Author:   Bas Retsios, retsios@itc.nl
 *
 ******************************************************************************
 * Copyright (c) 2004, ITC
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************/

 #include "cpl_port.h"  // Must be first.

#include "reflectancecalculator.h"
#include <cmath>
#include <cstdlib>

CPL_CVSID("$Id: reflectancecalculator.cpp 3b0bbf7a8a012d69a783ee1f9cfeb5c52b370021 2017-06-27 20:57:02Z Even Rouault $")

using namespace std;

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

ReflectanceCalculator::ReflectanceCalculator(std::string sTimeStamp, double rRTOA)
: m_rRTOA(rRTOA)
{
  std::string sYear (sTimeStamp.substr(0, 4));
  std::string sMonth (sTimeStamp.substr(4, 2));
  std::string sDay (sTimeStamp.substr(6, 2));
  std::string sHours (sTimeStamp.substr(8, 2));
  std::string sMins (sTimeStamp.substr(10, 2));

  m_iYear = atoi(sYear.c_str());
  int iMonth = atoi(sMonth.c_str());
  m_iDay = atoi(sDay.c_str());
  for (int i = 1; i < iMonth; ++i)
      m_iDay += iDaysInMonth(i, m_iYear);
  int iHours = atoi(sHours.c_str());
  int iMins = atoi(sMins.c_str());

        m_rHours = iHours + iMins / 60.0;
}

ReflectanceCalculator::~ReflectanceCalculator() {}

double ReflectanceCalculator::rGetReflectance(double rRadiance, double rLat, double rLon) const
{
  double phi = rLat * M_PI / 180;
  //double lam = rLon * M_PI / 180;
  double rSunDist = rSunDistance();
  double ReflectanceNumerator = rRadiance*rSunDist*rSunDist;
  double zenithAngle = rZenithAngle(phi, rDeclination(), rHourAngle(rLon));
  double ReflectanceDenominator = m_rRTOA*cos(zenithAngle*M_PI/180);
  double Reflectance = ReflectanceNumerator / ReflectanceDenominator;
  return Reflectance;
}

double ReflectanceCalculator::rZenithAngle(double phi, double rDeclin, double l_rHourAngle)
{
  double rCosZen = (sin(phi) * sin(rDeclin) + cos(phi)
          * cos(rDeclin) * cos(l_rHourAngle));
  double zenithAngle = acos(rCosZen) * 180 / M_PI;
  return zenithAngle;
}

double ReflectanceCalculator::rDeclination() const
{
  double rJulianDay = m_iDay - 1;
  double yearFraction = (rJulianDay + m_rHours / 24) / iDaysInYear(m_iYear);
  double T = 2 * M_PI * yearFraction;

  double declin = 0.006918 - 0.399912 * cos(T) + 0.070257 * sin(T)
          - 0.006758 * cos(2 * T) + 0.000907 * sin(2 * T)
          - 0.002697 * cos(3 * T) + 0.00148 * sin(3 * T);
  return declin;
}

double ReflectanceCalculator::rHourAngle(double rLon) const
{
  // In: rLon (in degrees)
  // Out: hourAngle (in radians)
  double rJulianDay = m_iDay - 1;
  double yearFraction = (rJulianDay + m_rHours / 24) / iDaysInYear(m_iYear);
  double T = 2 * M_PI * yearFraction;

  double EOT2 = 229.18 * (0.000075 + 0.001868 * cos(T)- 0.032077 * sin(T));
  double EOT3 = 229.18 * (- 0.014615 * cos(2 * T) - 0.040849 * sin(2 * T));
  double EOT = EOT2 + EOT3;
  double TimeOffset = EOT + (4. * rLon);
  // True solar time in minutes:
  double TrueSolarTime = m_rHours * 60 + TimeOffset;
  // Solar hour angle in degrees and in radians:
  double HaDegr = (TrueSolarTime / 4. - 180.);
  double hourAngle = HaDegr * M_PI / 180;
  return hourAngle;
}

double ReflectanceCalculator::rSunDistance() const
{
  int iJulianDay = m_iDay - 1;
  double theta = 2*M_PI *(iJulianDay - 3) / 365.25;
        // rE0 is the inverse of the square of the sun-distance ratio
        double rE0 = 1.000110 + 0.034221*cos(theta)+0.00128*sin(theta) + 0.000719*cos(2*theta)+0.000077*sin(2*theta);
  // The calculated distance is expressed as a factor of the "average sun-distance" (on 1 Jan approx. 0.98, on 1 Jul approx. 1.01)
  return 1 / sqrt(rE0);
}

int ReflectanceCalculator::iDaysInYear(int iYear)
{
  bool fLeapYear = iDaysInMonth(2, iYear) == 29;

  if (fLeapYear)
      return 366;
  else
      return 365;
}

int ReflectanceCalculator::iDaysInMonth(int iMonth, int iYear)
{
  int iDays;

  if ((iMonth == 4) || (iMonth == 6) || (iMonth == 9) || (iMonth == 11))
    iDays = 30;
  else if (iMonth == 2)
  {
    iDays = 28;
    if (iYear % 100 == 0) // century year
    {
      if (iYear % 400 == 0) // century leap year
        ++iDays;
    }
    else
    {
      if (iYear % 4 == 0) // normal leap year
        ++iDays;
    }
  }
  else
    iDays = 31;

  return iDays;
}