Logo Search packages:      
Sourcecode: saga version File versions  Download package

SADO_SolarRadiation.cpp


///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                    Module Library:                    //
//                     ta_lighting                       //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                SADO_SolarRadiation.cpp                //
//                                                       //
//                 Copyright (C) 2006 by                 //
//                      Olaf Conrad                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'. SAGA is free software; you   //
// can redistribute it and/or modify it under the terms  //
// of the GNU General Public License as published by the //
// Free Software Foundation; version 2 of the License.   //
//                                                       //
// SAGA is distributed in the hope that it will be       //
// useful, but WITHOUT ANY WARRANTY; without even the    //
// implied warranty of MERCHANTABILITY or FITNESS FOR A  //
// PARTICULAR PURPOSE. See the GNU General Public        //
// License for more details.                             //
//                                                       //
// You should have received a copy of the GNU General    //
// Public License along with this program; if not,       //
// write to the Free Software Foundation, Inc.,          //
// 59 Temple Place - Suite 330, Boston, MA 02111-1307,   //
// USA.                                                  //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     oconrad@saga-gis.org                   //
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "SADO_SolarRadiation.h"


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSADO_SolarRadiation::CSADO_SolarRadiation(void)
{
      CSG_Parameter     *pNode;

      //-----------------------------------------------------
      Set_Name          (_TL("Insolation"));

      Set_Author        (SG_T("(c) 2006 by O.Conrad"));

      Set_Description   (_TW(
            "Calculation of incoming solar radiation (insolation). "
            "Based on the SADO (System for the Analysis of Discrete Surfaces) routines "
            "developed by Boehner & Trachinow. "
            "\n\n"
            "References:\n"
      ));

      //-----------------------------------------------------
      Parameters.Add_Grid(
            NULL  , "GRD_DEM"             , _TL("Elevation"),
            _TL(""),
            PARAMETER_INPUT
      );

      Parameters.Add_Grid(
            NULL  , "GRD_VAPOUR"          , _TL("Water Vapour Pressure [mbar]"),
            _TL(""),
            PARAMETER_INPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "GRD_LAT"             , _TL("Latitude [degree]"),
            _TL(""),
            PARAMETER_INPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "GRD_LON"             , _TL("Longitude [degree]"),
            _TL(""),
            PARAMETER_INPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "GRD_DIRECT"          , _TL("Direct Insolation"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "GRD_DIFFUS"          , _TL("Diffuse Insolation"),
            _TL(""),
            PARAMETER_OUTPUT_OPTIONAL
      );

      Parameters.Add_Grid(
            NULL  , "GRD_TOTAL"           , _TL("Total Insolation"),
            _TL("Total insolation, the sum of direct and diffuse incoming solar radiation."),
            PARAMETER_OUTPUT_OPTIONAL
      );

      //-----------------------------------------------------
      Parameters.Add_Value(
            NULL  , "BHORIZON"            , _TL("Horizontal Projection"),
            _TL(""),
            PARAMETER_TYPE_Bool                 , false
      );

      Parameters.Add_Value(
            NULL  , "SOLARCONST"          , _TL("Solar Constant [J]"),
            _TL(""),
            PARAMETER_TYPE_Double         , 8.164, 0.0, true
      );

      Parameters.Add_Value(
            NULL  , "ATMOSPHERE"          , _TL("Height of Atmosphere [m]"),
            _TL(""),
            PARAMETER_TYPE_Double         , 12000.0, 0.0, true
      );

      Parameters.Add_Value(
            NULL  , "VAPOUR"              , _TL("Water Vapour Pressure [mbar]"),
            _TL("This value is used as constant if no vapour pressure grid is given."),
            PARAMETER_TYPE_Double         , 10, 0.0, true
      );

      //-----------------------------------------------------
      pNode = Parameters.Add_Node(
            NULL  , "NODE_TIME"           , _TL("Time"),
            _TL("")
      );

      Parameters.Add_Choice(
            pNode , "PERIOD"              , _TL("Time Period"),
            _TL(""),

            CSG_String::Format(SG_T("%s|%s|%s|%s|"),
                  _TL("moment"),
                  _TL("day"),
                  _TL("range of days"),
                  _TL("same moment for a range of days")
            ), 1
      );

      Parameters.Add_Value(
            pNode , "DHOUR"               , _TL("Daily Time Resolution [h]"),
            _TL("Time resolution for a day's calculation."),
            PARAMETER_TYPE_Double         , 1.0, 0.0, true, 24.0, true
      );

      Parameters.Add_Value(
            pNode , "DDAYS"               , _TL("Range of Days Time Resolution [days]"),
            _TL("Time resolution for a range of days calculation."),
            PARAMETER_TYPE_Int                  , 1, 1, true
      );

      Parameters.Add_Value(
            pNode , "DAY_A"               , _TL("Day of Year"),
            _TL(""),
            PARAMETER_TYPE_Int                  , 81
      );

      Parameters.Add_Value(
            pNode , "DAY_B"               , _TL("Day of Year (Range End)"),
            _TL(""),
            PARAMETER_TYPE_Int                  , 87
      );

      Parameters.Add_Value(
            pNode , "MOMENT"              , _TL("Moment [h]"),
            _TL(""),
            PARAMETER_TYPE_Double         , 12.0, 0.0, true, 24.0, true
      );

      //-----------------------------------------------------
      pNode = Parameters.Add_Node(
            NULL  , "NODE_LOCATION" , _TL("Location"),
            _TL("")
      );

      Parameters.Add_Value(
            pNode , "LATITUDE"            , _TL("Latitude"),
            _TL(""),
            PARAMETER_TYPE_Degree         , 53.0, -90.0, true, 90.0, true
      );

      Parameters.Add_Value(
            pNode , "BENDING"             , _TL("Planetery Bending"),
            _TL(""),
            PARAMETER_TYPE_Bool                 , false
      );

      Parameters.Add_Value(
            pNode , "RADIUS"              , _TL("Planetary Radius"),
            _TL(""),
            PARAMETER_TYPE_Double         , 6366737.96, 0.0, true
      );

      Parameters.Add_Choice(
            pNode , "LAT_OFFSET"          , _TL("Latitude relates to grid's..."),
            _TL(""),

            CSG_String::Format(SG_T("%s|%s|%s|%s|"),
                  _TL("bottom"),
                  _TL("center"),
                  _TL("top"),
                  _TL("user defined reference")
            ), 3
      );

      Parameters.Add_Choice(
            pNode , "LON_OFFSET"          , _TL("Local time relates to grid's..."),
            _TL(""),

            CSG_String::Format(SG_T("%s|%s|%s|%s|"),
                  _TL("left"),
                  _TL("center"),
                  _TL("right"),
                  _TL("user defined reference")
            ), 1
      );

      Parameters.Add_Value(
            pNode , "LAT_REF_USER"  , _TL("Latitude (user defined reference)"),
            _TL(""),
            PARAMETER_TYPE_Double         , 0.0
      );

      Parameters.Add_Value(
            pNode , "LON_REF_USER"  , _TL("Local Time (user defined reference)"),
            _TL(""),
            PARAMETER_TYPE_Double         , 0.0
      );

      //-----------------------------------------------------
      pNode = Parameters.Add_Node(
            NULL  , "NODE_UPDATEVIEW"     , _TL("Update View"),
            _TL("")
      );

      Parameters.Add_Value(
            pNode , "UPD_DIRECT"          , _TL("Direct Insolation"),
            _TL(""),
            PARAMETER_TYPE_Bool                 , false
      );

      Parameters.Add_Value(
            pNode , "UPD_DIFFUS"          , _TL("Diffuse Insolation"),
            _TL(""),
            PARAMETER_TYPE_Bool                 , false
      );

      Parameters.Add_Value(
            pNode , "UPD_TOTAL"           , _TL("Total Insolation"),
            _TL(""),
            PARAMETER_TYPE_Bool                 , false
      );
}

//---------------------------------------------------------
CSADO_SolarRadiation::~CSADO_SolarRadiation(void)
{}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSADO_SolarRadiation::On_Execute(void)
{
      //-----------------------------------------------------
      m_pDEM                  = Parameters("GRD_DEM")       ->asGrid();
      m_pVP             = Parameters("GRD_VAPOUR")    ->asGrid();

      m_pSumDirect      = Parameters("GRD_DIRECT")    ->asGrid();
      m_bUpdateDirect   = Parameters("UPD_DIRECT")    ->asBool() && m_pSumDirect;
      m_pSumDiffus      = Parameters("GRD_DIFFUS")    ->asGrid();
      m_bUpdateDiffus   = Parameters("UPD_DIFFUS")    ->asBool() && m_pSumDiffus;
      m_pSumTotal       = Parameters("GRD_TOTAL")     ->asGrid();
      m_bUpdateTotal    = Parameters("UPD_TOTAL")     ->asBool() && m_pSumTotal;

      m_Solar_Const     = Parameters("SOLARCONST")    ->asDouble();
      m_Atmosphere      = Parameters("ATMOSPHERE")    ->asDouble();
      m_VP              = Parameters("VAPOUR")        ->asDouble();
      m_bHorizon        = Parameters("BHORIZON")      ->asBool();

      m_Latitude        = Parameters("LATITUDE")      ->asDouble() * M_DEG_TO_RAD;
      m_bBending        = Parameters("BENDING")       ->asBool();

      m_Hour                  = Parameters("MOMENT")        ->asDouble();
      m_dHour                 = Parameters("DHOUR")         ->asDouble();
      m_dDays                 = Parameters("DDAYS")         ->asInt();
      m_Day_A                 = Parameters("DAY_A")         ->asInt();
      m_Day_B                 = Parameters("DAY_B")         ->asInt();

      switch( Parameters("PERIOD")->asInt() )
      {
      case 0:     // moment
            m_bMoment         = true;
            m_bUpdateDirect   = m_bUpdateDiffus = m_bUpdateTotal  = false;
            break;

      case 1:     // day
            m_bMoment         = false;
            m_Day_B                 = m_Day_A;
            m_Hour                  = 0.0;
            break;

      case 2:     // range of days
            m_bMoment         = false;
            m_Hour                  = 0.0;
            break;

      case 3:     // same moment for a range of days
            m_bMoment         = false;
            m_dHour                 = 24.0;
            break;
      }

      //-----------------------------------------------------
      if( !m_pSumDirect && !m_pSumDiffus && !m_pSumTotal )
      {
            Message_Dlg(_TL("No output grid has been choosen."), Get_Name());
      }
      else if( m_pSumDirect == m_pDEM || m_pSumDiffus == m_pDEM || m_pSumTotal == m_pDEM )
      {
            Message_Dlg(_TL("Output must not overwrite elevation grid."), Get_Name());
      }
      else
      {
            return( Get_Insolation() );
      }

      //-----------------------------------------------------
      return( false );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSADO_SolarRadiation::Initialise(void)
{
      int         x, y;

      Process_Set_Text(_TL("initialising..."));

      //-----------------------------------------------------
      CSG_Colors  c(100, SG_COLORS_YELLOW_RED, true);

      c.Set_Ramp(SG_GET_RGB(  0,   0,  64), SG_GET_RGB(255, 159,   0),  0, 50);
      c.Set_Ramp(SG_GET_RGB(255, 159,   0), SG_GET_RGB(255, 255, 255), 50, 99);

      if( m_pSumDirect )
      {
            m_pSumDirect->Assign(0.0);
            m_pSumDirect->Set_Unit(_TL("Joule"));
            DataObject_Set_Colors(m_pSumDirect, c);

            if( m_bUpdateDirect )
            {
                  m_TmpDirect.Create(*Get_System(), SG_DATATYPE_Float);
                  DataObject_Update(m_pSumDirect, true);
            }
      }

      if( m_pSumDiffus )
      {
            m_pSumDiffus->Assign(0.0);
            m_pSumDiffus->Set_Unit(_TL("Joule"));
            DataObject_Set_Colors(m_pSumDiffus, c);

            if( m_bUpdateDiffus )
            {
                  m_TmpDiffus.Create(*Get_System(), SG_DATATYPE_Float);
                  DataObject_Update(m_pSumDiffus, true);
            }
      }

      if( m_pSumTotal )
      {
            m_pSumTotal ->Assign(0.0);
            m_pSumTotal ->Set_Unit(_TL("Joule"));
            DataObject_Set_Colors(m_pSumTotal , c);

            if( m_bUpdateTotal )
            {
                  m_TmpTotal.Create(*Get_System(), SG_DATATYPE_Float);
                  DataObject_Update(m_pSumTotal , true);
            }
      }

      //-----------------------------------------------------
      Process_Set_Text(_TL("initialising gradient..."));

      m_Shade .Create(*Get_System(), SG_DATATYPE_Byte);
      m_Slope .Create(*Get_System(), SG_DATATYPE_Float);
      m_Aspect.Create(*Get_System(), SG_DATATYPE_Float);

      for(y=0; y<Get_NY() && Set_Progress(y); y++)
      {
            for(x=0; x<Get_NX(); x++)
            {
                  double      s, a;

                  if( m_pDEM->Get_Gradient(x, y, s, a) )
                  {
                        m_Slope .Set_Value(x, y, s);
                        m_Aspect.Set_Value(x, y, a);
                  }
                  else
                  {
                        m_Slope .Set_NoData(x, y);
                        m_Aspect.Set_NoData(x, y);
                  }
            }
      }

      //-----------------------------------------------------
      if( m_bBending )
      {
            Process_Set_Text(_TL("initialising planetary bending..."));

            CSG_Grid    *pLat = Parameters("GRD_LAT")->asGrid(),
                        *pLon = Parameters("GRD_LON")->asGrid();

            m_Lat       .Create(*Get_System(), SG_DATATYPE_Float);
            m_Lon       .Create(*Get_System(), SG_DATATYPE_Float);
            m_Decline   .Create(*Get_System(), SG_DATATYPE_Float);
            m_Azimuth   .Create(*Get_System(), SG_DATATYPE_Float);

            if( pLat || pLon )
            {
                  if( pLat )
                  {
                        m_Lat = *pLat;
                        m_Lat *= M_DEG_TO_RAD;
                  }

                  if( pLon )
                  {
                        m_Lon = *pLon;
                        m_Lon *= M_DEG_TO_RAD;
                  }
            }

            //-------------------------------------------------
            else
            {
                  double      d, dx, dy, dxA, dyA;

                  d     = M_DEG_TO_RAD / (Parameters("RADIUS")->asDouble() * M_PI / 180.0);

                  switch( Parameters("LON_OFFSET")->asInt() )
                  {
                  case 0:     dxA   = Get_System()->Get_Extent().Get_XMin();  break;      // left
                  case 1:     dxA   = Get_System()->Get_Extent().Get_XCenter();     break;      // center
                  case 2:     dxA   = Get_System()->Get_Extent().Get_XMax();  break;      // right
                  case 3:     dxA   = Parameters("LON_REF_USER")->asDouble(); break;      // user defined coordinate
                  }

                  switch( Parameters("LAT_OFFSET")->asInt() )
                  {
                  case 0:     dyA   = Get_System()->Get_Extent().Get_YMin();  break;      // bottom
                  case 1:     dyA   = Get_System()->Get_Extent().Get_YCenter();     break;      // center
                  case 2:     dyA   = Get_System()->Get_Extent().Get_YMax();  break;      // top
                  case 3:     dyA   = Parameters("LAT_REF_USER")->asDouble(); break;      // user defined coordinate
                  }

                  dxA    = d * (Get_XMin() - dxA);
                  dyA    = d * (Get_YMin() - dyA) + m_Latitude;
                  d     *= Get_Cellsize();

                  for(y=0, dy=dyA; y<Get_NY() && Set_Progress(y); y++, dy+=d)
                  {
                        for(x=0, dx=dxA; x<Get_NX(); x++, dx+=d)
                        {
                              m_Lat.Set_Value(x, y, dy);
                              m_Lon.Set_Value(x, y, dx);
                        }
                  }
            }
      }

      //-----------------------------------------------------
      return( true );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Finalise(double SumFactor)
{
      //-----------------------------------------------------
      if( m_pSumDirect )
      {
            if( m_bUpdateDirect )   m_pSumDirect->Assign(&m_TmpDirect);
            if( SumFactor != 1.0 )  *m_pSumDirect     *= SumFactor;
      }

      if( m_pSumDiffus )
      {
            if( m_bUpdateDiffus )   m_pSumDiffus->Assign(&m_TmpDiffus);
            if( SumFactor != 1.0 )  *m_pSumDiffus     *= SumFactor;
      }

      if( m_pSumTotal )
      {
            if( m_bUpdateTotal )    m_pSumTotal->Assign(&m_TmpTotal);
            if( SumFactor != 1.0 )  *m_pSumTotal      *= SumFactor;
      }

      //-----------------------------------------------------
      m_Shade           .Destroy();
      m_Slope           .Destroy();
      m_Aspect    .Destroy();

      m_Lat       .Destroy();
      m_Lon       .Destroy();
      m_Decline   .Destroy();
      m_Azimuth   .Destroy();

      m_TmpDirect .Destroy();
      m_TmpDiffus .Destroy();
      m_TmpTotal  .Destroy();

      //-----------------------------------------------------
      return( true );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Insolation(void)
{
      //-----------------------------------------------------
      if( Initialise() )
      {
            if( m_bMoment )
            {
                  Get_Insolation(m_Day_A, m_Hour);

                  Finalise();
            }

            //-------------------------------------------------
            else
            {
                  for(int Day=m_Day_A; Day<=m_Day_B && Process_Get_Okay(false); Day+=m_dDays)
                  {
                        for(double Hour=m_Hour; Hour<24.0 && Process_Get_Okay(false); Hour+=m_dHour)
                        {
                              Process_Set_Text(CSG_String::Format(SG_T("%s: %d(%d-%d), %s: %f"), _TL("Day"), Day, m_Day_A, m_Day_B, _TL("Hour"), Hour));

                              if( m_bUpdateDirect )   m_pSumDirect->Assign(0.0);
                              if( m_bUpdateDiffus )   m_pSumDiffus->Assign(0.0);
                              if( m_bUpdateTotal )    m_pSumTotal ->Assign(0.0);

                              if( Get_Insolation(Day, Hour) )
                              {
                                    if( m_bUpdateDirect )
                                    {
                                          m_TmpDirect += *m_pSumDirect;
                                          DataObject_Update(m_pSumDirect);
                                    }

                                    if( m_bUpdateDiffus )
                                    {
                                          m_TmpDiffus += *m_pSumDiffus;
                                          DataObject_Update(m_pSumDiffus);
                                    }

                                    if( m_bUpdateTotal )
                                    {
                                          m_TmpTotal  += *m_pSumTotal;
                                          DataObject_Update(m_pSumTotal);
                                    }
                              }
                        }
                  }

                  Finalise(m_dHour / (24.0 * (1 + m_Day_B - m_Day_A))); // *m_pSumDirect  *= m_dHour / D->size();
            }
      }

      //-----------------------------------------------------
      return( true );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Insolation(int Day, double Hour)
{
      double      Azimuth, Decline;

      if( m_bBending )
      {
            bool  bGo   = false;

            for(int y=0; y<Get_NY() && Set_Progress(y); y++)
            {
                  for(int x=0; x<Get_NX(); x++)
                  {
                        Get_Solar_Position(Day, Hour, m_Lat.asDouble(x, y), m_Lon.asDouble(x, y), Decline, Azimuth);

                        m_Azimuth.Set_Value(x, y, Azimuth);
                        m_Decline.Set_Value(x, y, Decline);

                        if( Decline > 0.0 )
                        {
                              bGo   = true;
                        }
                  }
            }

            if( bGo )
            {
                  return( Set_Insolation(Decline, Azimuth) );
            }
      }
      else if( Get_Solar_Position(Day, Hour, m_Latitude, 0.0, Decline, Azimuth) )
      {
            return( Set_Insolation(Decline, Azimuth) );
      }

      return( false );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Vapour_Exponent(double VapourPressure)
{
      return( 0.916  - (VapourPressure > 0.0 ? 0.05125 * sqrt(VapourPressure) : 0.0) );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Vapour_A(double VapourPressure)
{
      return( 0.4158 + (VapourPressure > 0.0 ? 0.0399  * sqrt(VapourPressure) : 0.0) );
}

//---------------------------------------------------------
bool CSADO_SolarRadiation::Set_Insolation(double Decline, double Azimuth)
{
      int         x, y;
      double      a, e, Direct, Diffus;

      Get_Shade(Decline, Azimuth);

      e     = Get_Vapour_Exponent   (m_VP);
      a     = Get_Vapour_A                (m_VP);

      for(y=0; y<Get_NY() && Set_Progress(y); y++)
      {
            for(x=0; x<Get_NX(); x++)
            {
                  if( m_pDEM->is_NoData(x, y) )
                  {
                        if( m_pSumDirect )      m_pSumDirect->Set_NoData(x, y);
                        if( m_pSumDiffus )      m_pSumDiffus->Set_NoData(x, y);
                        if( m_pSumTotal ) m_pSumTotal ->Set_NoData(x, y);
                  }
                  else
                  {
                        if( m_pVP != NULL )
                        {
                              e     = Get_Vapour_Exponent   (m_pVP->asDouble(x, y));
                              a     = Get_Vapour_A                (m_pVP->asDouble(x, y));
                        }

                        if( m_bBending )
                        {
                              Azimuth     = m_Azimuth.asDouble(x, y);
                              Decline     = m_Decline.asDouble(x, y);
                        }

                        if( Decline > 0.0 )
                        {
                              Direct      = Get_Solar_Direct(x, y, Decline, Azimuth, e);
                              Diffus      = Get_Solar_Diffus(x, y, Decline, a      , e);

                              if( m_pSumDirect )      m_pSumDirect->Add_Value(x, y, Direct);
                              if( m_pSumDiffus )      m_pSumDiffus->Add_Value(x, y, Diffus);
                              if( m_pSumTotal  )      m_pSumTotal ->Add_Value(x, y, Direct + Diffus);
                        }
                  }
            }
      }

      return( true );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Reduction(double Elevation, double Decline, double Reduction)
{
      return( Decline > 0.0 ? (m_Solar_Const * pow(Reduction, (1.0 - Elevation / m_Atmosphere) / sin(Decline))) : 0.0 );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Direct(int x, int y, double Decline, double Azimuth, double Exponent)
{
      if( m_Shade.asInt(x, y) == 0 )
      {
            double      Angle, Slope;

            Slope = m_bHorizon ? 0.0 : m_Slope.asDouble(x, y);

            Angle = cos(Slope) * cos(Decline - M_PI_090)
                        + sin(Slope) * sin(M_PI_090 - Decline)
                        * cos(Azimuth - m_Aspect.asDouble(x, y));

            if( Angle > 0.0 )
            {
                  return( Angle * Get_Solar_Reduction(m_pDEM->asDouble(x, y), Decline, Exponent) );
            }
      }

      return( 0.0 );
}

//---------------------------------------------------------
inline double CSADO_SolarRadiation::Get_Solar_Diffus(int x, int y, double Decline, double A, double Exponent)
{
      double      Dh, Ih, d, Slope, srad;

      Ih          = Get_Solar_Reduction(m_pDEM->asDouble(x, y), Decline, Exponent) * sin(Decline);
//    Decline     *= M_RAD_TO_DEG;  // Fehler in SADO ????!!!!!
      d           = m_Atmosphere / (m_Atmosphere - m_pDEM->asDouble(x, y))
                  * (0.0001165 * Decline*Decline - 0.0152 * Decline + A);
      d           = 1.0 - d;
      Dh          = Ih / d - Ih;

      Slope = m_bHorizon ? 0.0 : m_Slope.asDouble(x, y);

      srad  = (Dh + Dh * cos(Slope)) / 2.0;

      return( srad < 0.0 || srad > m_Solar_Const ? m_Solar_Const : srad );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Shade(double Decline, double Azimuth)
{
      //-----------------------------------------------------
      m_Shade.Assign(0.0);

      if( !m_bBending )
      {
            int         i, x, y;
            double      dx, dy, dz;

            Get_Shade_Params(Decline, Azimuth, dx, dy, dz);

            for(i=0; i<Get_NCells() && Set_Progress_NCells(i); i++)
            {
                  if( m_pDEM->Get_Sorted(i, x, y) && !Get_Shade_Complete(x, y) )
                  {
                        Set_Shade(x, y, dx, dy, dz);
                  }
            }
      }

      //-----------------------------------------------------
      else
      {
            int         i, x, y, iLock;

            for(i=0, iLock=1; i<Get_NCells() && Set_Progress_NCells(i); i++, iLock++)
            {
                  if( m_pDEM->Get_Sorted(i, x, y) && !Get_Shade_Complete(x, y) )
                  {
                        if( iLock >= 255 )
                              iLock = 1;

                        if( iLock == 1 )
                              Lock_Create();

                        Set_Shade_Bended(x, y, iLock);
                  }
            }
      }

      return( true );
}

//---------------------------------------------------------
void CSADO_SolarRadiation::Set_Shade(int x, int y, double dx, double dy, double dz)
{
      for(double ix=x+0.5, iy=y+0.5, iz=m_pDEM->asDouble(x, y); ; )
      {
            x     = (int)(ix  += dx);
            y     = (int)(iy  += dy);
                              iz    -= dz;

            if( !is_InGrid(x, y) || m_pDEM->asDouble(x, y) > iz )
            {
                  return;
            }

            m_Shade.Set_Value(x, y, 1);
      }
}

//---------------------------------------------------------
void CSADO_SolarRadiation::Set_Shade_Bended(int x, int y, char iLock)
{
      double      dx, dy, dz;

      Get_Shade_Params(m_Decline.asDouble(x, y), m_Azimuth.asDouble(x, y), dx, dy, dz);

      for(double ix=x+0.5, iy=y+0.5, iz=m_pDEM->asDouble(x, y); ; )
      {
            x     = (int)(ix  += dx);
            y     = (int)(iy  += dy);
                              iz    -= dz;

            if( !is_InGrid(x, y) || m_pDEM->asDouble(x, y) > iz || Lock_Get(x, y) == iLock )
            {
                  return;
            }

            m_Shade.Set_Value(x, y, 1);

            //---------------------------------------------
            Lock_Set(x, y, iLock);

            Get_Shade_Params(m_Decline.asDouble(x, y), m_Azimuth.asDouble(x, y), dx, dy, dz);
      }
}

//---------------------------------------------------------
inline bool CSADO_SolarRadiation::Get_Shade_Complete(int x, int y)
{
      if( m_Shade.asInt(x, y) == 1 )
      {
            for(int iy=y-1; iy<=y+1; iy++)
            {
                  for(int ix=x-1; ix<x+1; ix++)
                  {
                        if( is_InGrid(ix, iy) && m_Shade.asInt(ix, iy) == 0 )
                        {
                              return( false );
                        }
                  }
            }

            return( true );
      }

      return( false );
}

//---------------------------------------------------------
inline void CSADO_SolarRadiation::Get_Shade_Params(double Decline, double Azimuth, double &dx, double &dy, double &dz)
{
      dz    = Azimuth + M_PI_180;
      dx    = sin(dz);
      dy    = cos(dz);

      if( fabs(dx) > fabs(dy) )
      {
            dy    /= fabs(dx);
            dx    = dx < 0 ? -1 : 1;
      }
      else if( fabs(dy) > fabs(dx) )
      {
            dx    /= fabs(dy);
            dy    = dy < 0 ? -1 : 1;
      }
      else
      {
            dx    = dx < 0 ? -1 : 1;
            dy    = dy < 0 ? -1 : 1;
      }

      dz    = tan(Decline) * sqrt(dx*dx + dy*dy) * Get_Cellsize();
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
inline bool CSADO_SolarRadiation::Get_Solar_Position(int Day, double Hour, double LAT, double LON, double &Declination, double &Azimuth)
{
      static const int  Day2Month[13]     = {   0, 31, 49, 80, 109, 140, 170, 201, 232, 272, 303, 333, 366 };

      static const double     ECLIPTIC_OBL      = M_DEG_TO_RAD * 23.43999;    // obliquity of ecliptic

      int         i;

      double      JD, T, M, L, X, Y, Z, R, UTime,
                  DEC, RA, theta, tau,
                  Month, Year       = 2000;

      //-----------------------------------------------------

      for(Month=1, i=0; i<=12; i++)
      {
            if( Day < Day2Month[i] )
            {
                  Month = i;
                  Day         -= Day2Month[i - 1];
                  break;
            }
      }

      if( Month < 1 || Month > 12 ) // nur Tag (1 - 365) zaehlt...
      {
            Month       = 1;
      }

      if( Month <= 2 )
      {
            Month       += 12;
            Year        -= 1;
      }


      //-----------------------------------------------------

//    UTime       = Hour - LON * 12.0 / M_PI;
      UTime       = Hour;


      //-----------------------------------------------------
      // 1. Julian Date...

      JD                = (int)(365.25 * Year) + (int)(30.6001 * (Month + 1)) - 15 + 1720996.5 + Day + UTime / 24.0;
      T                 = (JD - 2451545.0 ) / 36525.0;      // number of Julian centuries since 2000/01/01 at 12 UT (JD = 2451545.0)


      //-----------------------------------------------------
      // 2. Solar Coordinates (according to: Jean Meeus: Astronomical Algorithms), accuracy of 0.01 degree

      M                 = M_DEG_TO_RAD * (357.52910 + 35999.05030 * T - 0.0001559 * T * T - 0.00000048 * T * T * T);    // mean anomaly
      L                 = M_DEG_TO_RAD * (      (280.46645 + 36000.76983 * T + 0.0003032 * T * T)                                   // mean longitude
                                          +     (     (1.914600 - 0.004817 * T - 0.000014  * T * T) * sin(M)
                                                +     (0.019993 - 0.000101 * T) * sin(2 * M) + 0.000290 * sin(3 * M)          // true longitude
                                                )
                                          );


      //-----------------------------------------------------
      // 3. convert ecliptic longitude to right ascension RA and declination delta

      X                 = cos(L);
      Y                 = cos(ECLIPTIC_OBL) * sin(L);
      Z                 = sin(ECLIPTIC_OBL) * sin(L);
      R                 = sqrt(1.0 - Z*Z); 

      DEC               = atan2(Z, R);
      RA                = 2.0 * atan2(Y, (X + R));


      //-----------------------------------------------------
      // 4. compute sidereal time (degree) at Greenwich local sidereal time at longitude (Degree)

      theta       = LON + M_DEG_TO_RAD * (280.46061837 + 360.98564736629 * (JD - 2451545.0) + T*T * (0.000387933 - T / 38710000.0));


      //-----------------------------------------------------
      // 5. compute local hour angle (degree)

      tau               = theta - RA;


      //-----------------------------------------------------
      // 6. convert (tau, delta) to horizon coordinates (h, az) of the observer

      Declination = asin ( sin(LAT) * sin(DEC) + cos(LAT) * cos(DEC) * cos(tau) );

      Azimuth           = atan2( -sin(tau) * cos(DEC), cos(LAT) * sin(DEC) - sin(LAT) * cos(DEC) * cos(tau) );
//    Azimuth           = atan2( -sin(Tau), cos(LAT) * tan(DEC) - sin(LAT) * cos(Tau) );  // previous formula gives same result but is better because of division by zero effects...

      return( Declination > 0.0 );
}


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

/*/---------------------------------------------------------
#define CHECK_SHADOW(X, Y)    if( (d = SG_Get_Distance(x, y, X, Y) * Get_Cellsize()) > 0.0 )\
                                          {\
                                                if( (t = (m_pDEM->asDouble(X, Y) - z) / d) > tMax )\
                                                {\
                                                      tMax  = t;\
                                                }\
                                                \
                                                if( (d = zMax / d) < tanDec || d < tMax )\
                                                {\
                                                      return( false );\
                                                }\
                                                else if( tMax > tanDec )\
                                                {\
                                                      return( true );\
                                                }\
                                          }

//---------------------------------------------------------
bool CSADO_SolarRadiation::Get_Solar_Shadow(int x, int y, double Decline, double Azimuth)
{
      static double     Zenith      = 89.9999 * M_DEG_TO_RAD;

      //-----------------------------------------------------
      if( Decline >= Zenith || Decline > m_Slope.Get_ZMax() )
      {
            return( false );
      }

      if( Decline <= 0.0 )
      {
            return( true );
      }

      //-----------------------------------------------------
      bool  bGo;
      int         dx, nx, ix, iy, iy_Last, jy;
      double      dy, d, t, tMax, z, zMax, tanDec;

      //-----------------------------------------------------
      Azimuth     = fmod(Azimuth, M_PI_360);    if( Azimuth < 0.0 )     Azimuth     += M_PI_360;

      if( Azimuth < M_PI_180 )
      {
            dx    = 1;
            nx    = Get_NX() - 1;
      }
      else
      {
            dx    = -1;
            nx    = 0;
      }

      Azimuth     = M_PI_090 - Azimuth;

      if(          Azimuth    ==  M_PI_090 )
                   Azimuth    =   89.9999 * M_DEG_TO_RAD;
      else if( Azimuth  == -M_PI_090 )
                   Azimuth    = - 89.9999 * M_DEG_TO_RAD;
      else if( Azimuth  == -M_PI_180 )
                   Azimuth    = -179.9999 * M_DEG_TO_RAD;
      else if( Azimuth  == -M_PI_270)
                   Azimuth    = -269.9999 * M_DEG_TO_RAD;

      dy          = tan(Azimuth);

      //-----------------------------------------------------
      z           = m_pDEM->asDouble(x, y);
      zMax  = m_pDEM->Get_ZMax();
      tanDec      = tan(Decline);
      tMax  = 0.0;
  
      //-----------------------------------------------------
      for(ix=x, iy=y, bGo=true; ix!=nx && bGo; ix+=dx)
      {
            iy_Last     = iy;
            iy          = y + (int)floor(dy * (ix - x + (dx > 0 ? 0.5 : -0.5)));

            if( iy < 0 )
            {
                  bGo   = false;
                  iy    = 0;
            }
            else if( iy >= Get_NY() )
            {
                  bGo   = false;
                  iy    = Get_NY() - 1;
            }

            //-------------------------------------------------
            if( iy_Last < iy )
            {
                  for(jy=iy_Last+1; jy<=iy; jy++)
                  {
                        CHECK_SHADOW(ix, jy);
                  }
            }
            else if( iy_Last > iy )
            {
                  for(jy=iy_Last; jy>iy; jy--)
                  {
                        CHECK_SHADOW(ix, jy);
                  }
            }
            else if( iy_Last == iy ) 
            {
                  CHECK_SHADOW(ix, iy);
            }
      }

      return( false );
}/**/


///////////////////////////////////////////////////////////
//                                                                                   //
//                                                                                   //
//                                                                                   //
///////////////////////////////////////////////////////////

//---------------------------------------------------------

Generated by  Doxygen 1.6.0   Back to index