/*
 * SUMMARY:      CalcSolar.c - inline solar calculations
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Mark Wigmosta
 * ORG:          Batelle Pacific Northwest Laboratories
 * E-MAIL:       ms_wigmosta@pnl.gov
 * ORIG-DATE:    Jul-96
 * Last Change: Wed Feb  3 23:11:44 1999 by Bart Nijssen <nijssen@u.washington.edu>
 * DESCRIPTION:  These functions make inline solar radiation calculations
 *               that take into slope and aspect of the pixel, but do not
 *               account for shadowing of neighbouring pixels.  These
 *               routines can be used instead of IPW, but it is up to the
 *               user to choose one method (IPW) or the other (inline
 *               calculations).  
 * DESCRIP-END.
 * FUNCTIONS:    SolarDay()
 *               SolarHour()
 *               SolarAngle()
 *               SolarConst()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: CalcSolar.c,v 1.2 2001/03/31 00:10:19 waichler Exp $";
#endif /* lint */

#include <math.h> 
/*#include "constants.h"*/
/*#include "settings.h" */
/*#include "Calendar.h"*/
#include "calcsolar.h"  
   
/*****************************************************************************
  Function name: SolarDay()

  Purpose:	 This subroutine calculates daily solar values 

  Required:
   int DayOfYear	  - day of year (January 1 is 1)
   float Longitude        - site longitude (rad)
   float Latitude         - site latitude (rad)
   float StandardMeridian - longitude of time zone of standard meridian (rad)

  Returns: void

  Modifies:
    float *NoonHour       - true solar noon (hr)
    float *Declination    - solar Declination (rad)
    float *HalfDayLength  - half-day length (hr)
    float *Sunrise        - time of Sunrise (hr)
    float *Sunset         - time of Sunset (hr)
    float *TimeAdjustment - required adjustment to local time (hr)
    float *SunEarthDist   - distance from sun to earth

  Comments     : EXECUTE AT START OF EACH DAY 
*****************************************************************************/
void SolarDay(int DayOfYear, float Longitude, float Latitude, 
              float StandardMeridian, float *NoonHour, float *Declination, 
              float *HalfDayLength, float *Sunrise, float *Sunset, 
              float *TimeAdjustment, float *SunEarthDist)
{
  float B;                      /* coefficient for equation of time */
  float EqnOfTime;              /* adjustment for equation of time (min) */
  float LongitudeAdjust;        /* adjustment for longitude (min) */
  float CosineHalfDayLength;    /* cosine of the half-day length */

 

  /* note need to check if day light savings time calculate adjustment for 
     true solar time longitude adjustment add 4 min per degree away from 
     StandardMeridian (4 min/degree * 180 degree/pi radian) */
  
  LongitudeAdjust = (MINPDEG * DEGPRAD)*(StandardMeridian - Longitude);
  
  /* equation of time */
  B = (2.0*PI*(DayOfYear - 81)) /  364.;
  EqnOfTime = 9.87*sin(2*B) - 7.53*cos(B) - 1.5*sin(B);
  
  /* adjustment factor to convert local time to solar time
     solar time = local time + TimeAdjustment  */
  *TimeAdjustment = -(LongitudeAdjust + EqnOfTime) / MINPHOUR;
  
  /* work in solar time  */
  *NoonHour = 12.0;
  
  /* solar Declinationation  */
  *Declination = .4098*sin(2*PI*(284 + DayOfYear) / DAYPYEAR);
  
  /* half-day length  */    
  CosineHalfDayLength = - tan(Latitude)*tan(*Declination);
  if (CosineHalfDayLength >= 1.0)
    *HalfDayLength = PI;
  else
    *HalfDayLength = acos(CosineHalfDayLength);
  
  /* convert HalfDayLength from radians to Hours 
     1 radian = (180 deg / PI) * (1 hr / 15 degrees rotation) */
  *HalfDayLength = *HalfDayLength / RADPHOUR;
  
  /* solar time of Sunrise and Sunset  */
  *Sunrise = *NoonHour - *HalfDayLength;
  *Sunset = *NoonHour + *HalfDayLength;
  
  /* calculate the sun-earth distance */
  *SunEarthDist = 1.0 + 0.033*cos( RADPDEG * (360. * DayOfYear  / 365 ) );
}

/*****************************************************************************
  Function name: SolarHour()  
  
  Purpose: This subroutine calculates position of the sun as a function of
           the time of day, the length of time the sun is above the horizon,
           and the maximum radiation.

  Required:
    float Latitude		- site laditude (rad)
    float LocalHour		- local time (hr)
    float Dt			- length of current timestep (hr)
    float NoonHour		- true solar noon (hr)
    float Declination		- solar Declination (rad)
    float Sunrise		- time of Sunrise (hr)
    float Sunset		- time of Sunset (hr)
    float TimeAdjustment	- required adjustment to convert local time
				  to solar time (hr)
    float SunEarthDist          - distance from Sun to Earth 

  Returns: void

  Modifies:
    float *SineSolarAltitude - sine of sun's SolarAltitude 
    int *DayLight	     - FALSE: measured solar radiation and the sun is
                                      below the horizon.  
			       TRUE: sun is above the horizon
    float *SolarTimeStep     - fraction of the timestep the sun is above the 
                               horizon 
    float *SunMax            - calculated solar radiation at the top of the 
                               atmosphere (W/m^2) 
                                                 
     float *SolarAzimuth    -  solar azimuth (rads eastward from north 
                                               N: 0 rads
                                               E: PI/2 rads
                                               S: PI rads
                                               W: 3PI/2 rads)

  Comments     : EXECUTE AT START OF EACH TIMESTEP 
*****************************************************************************/
void SolarHour(float Latitude, float LocalHour, float Dt, float NoonHour, 
               float Declination, float Sunrise, float Sunset, 
               float TimeAdjustment, float SunEarthDist, 
               float *SineSolarAltitude, int *DayLight, float *SolarTimeStep,
               float *SunMax, float *SolarAzimuth)
{
  float SolarAltitude;          /* SolarAltitude of sun from horizon (rads) */
  float SolarZenith;            /* sun zenith angle (rads) */
  float StartHour = 0;		/* currect Hour in solar time (hr) */
  float EndHour = 0;		/* mid-point of current solar Hour (hr) */
  float Hour;                   /* angle of current "halfhr" from solar noon
                                   (rads) */

  /* NOTE THAT HERE Dt IS IN HOURS NOT IN SECONDS */

  *SunMax = 0.0;
  
  *SolarTimeStep=1.0;
  
  /* all calculations in solar time correct local time */

  Hour = LocalHour + TimeAdjustment;
  
  *DayLight = FALSE;
  if ((Hour > Sunrise) && ((Hour-Dt) < Sunset))
    *DayLight = TRUE;
  
  if (*DayLight == TRUE) {
    /* sun is above the horizon */ 
    
    if( Dt > 0.0 ) {  
      /* compute average solar SolarAltitude over the timestep */
      StartHour = (Hour - Dt > Sunrise) ? Hour - Dt : Sunrise;
      EndHour = (Hour < Sunset) ? Hour : Sunset;
      
      /*  convert to radians  */
      StartHour = RADPHOUR * (StartHour - NoonHour);                   
      EndHour = RADPHOUR * (EndHour - NoonHour);   
      *SolarTimeStep = EndHour - StartHour;     
      
      /*  determine the average geometry of the sun angle  */
      *SineSolarAltitude = sin(Latitude)*sin(Declination)
        + (cos(Latitude)*cos(Declination)*(sin(EndHour) - sin(StartHour))
           / *SolarTimeStep);
    } 
    else {
      Hour = RADPHOUR * ( Hour - NoonHour );
      *SineSolarAltitude = sin(Latitude)*sin(Declination)
        + (cos(Latitude)*cos(Declination)*cos(Hour));
    }

    SolarAltitude = asin(*SineSolarAltitude);
    SolarZenith = PI/2-SolarAltitude;
    
    *SolarAzimuth = ((sin(Latitude)*(*SineSolarAltitude) - sin(Declination))
		    / (cos(Latitude)*sin(SolarZenith)));
    
    if (*SolarAzimuth > 1.) 
      *SolarAzimuth = 1.;
    
    *SolarAzimuth=acos(-(*SolarAzimuth));

    if (Dt > 0.0) {
      if (fabs(EndHour) > fabs(StartHour)) {
	*SolarAzimuth = 2*PI - (*SolarAzimuth);
      }
    }
    else {
      if (Hour > 0) {
	*SolarAzimuth = 2*PI - (*SolarAzimuth);
      }
    }
/*     printf("Solar Azimuth = %g\n", *SolarAzimuth*180./PI); */
    
    *SunMax = SOLARCON * SunEarthDist * *SineSolarAltitude;	
  }
}

/*****************************************************************************
  Function name: SolarAngle()  
  
  Purpose      : this subroutine uses solar radiation measured on a horizontal
                 surface to calculate direct, diffuse, and incoming reflected 
                 solar radiation for a sloping surface.

  Required     :
    float Latitude	        - site laditude (rad)
    float Albedo                - surface Albedo
    float Declination	        - solar Declination (rad)
    float CellAspect	        - cell aspect (rads eastward from north)
    float CellSlope	        - ground surface slope (rads)
    float SineSolarAltitude	- sine of sun's SolarAltitudede 
    float SunMax                - calculated solar radiation at the top of 
                                  the atmosphere (W/m^2) 
    float SolarTimeStep         - fraction of the timestep the sun is above 
                                  the horizon 
    int DayLight	        - FALSE: measured solar radiation and the sun
                                  is below the horizon.  
                                  TRUE: sun is above the horizon
     float SolarAzimuth    -  solar azimuth (rads eastward from north 
                                             N: 0 rads
                                             E: PI/2 rads
                                             S: PI rads
                                             W: 3PI/2 rads)
    float Dt			- length of current timestep (hr)

  Returns      : void

  Modifies     :
    float *Direct  - direct beam solar radiation (W/m^2)
    float *Diffuse - diffuse solar radiation (W/m^2)

  Comments     : EXECUTE EACH TIMESTEP FOR EACH GRID CELL 
                 Source:  Gates, D.M., "Biophysical ecology",
                 Springer-Verlag,  New York, etc.,  1980.
                 Especially Chapter 6, p. 136 ->


*****************************************************************************/
void SolarAngle(float Latitude, float Albedo, float Declination, 
                float CellAspect, float CellSlope, float SunMax, 
                float SineSolarAltitude, float SolarTimeStep, int DayLight,
                float SolarAzimuth, float Dt, float *Direct, float *Diffuse) 
{
  float SolarAltitude;          /* SolarAltitude of sun from horizon (rads) */
  float CosineIncidenceAngle;   /* cosine of the incidence angle between
                                   solar rays and the normal to the surface */ 
  float DiffuseSkyView;         /* sky view factor for diffuse radation 
                                   (0.0 - 1.0)  */
  float ReflectedSkyView;       /* view factor for incoming refected solar 
                                   radiation (0.0 - 1.0) */
  float Reflect;                 /* incoming reflected solar radiation (W/m^2)
                                 */ 

  /* NOTE THAT HERE Dt IS IN HOURS, NOT IN SECONDS */

  if (DayLight == TRUE) {
    
    /* calculate shortwave components horizontal surface */
    
    if (CellSlope == 0.0) {
      *Direct = SunMax;
      *Diffuse = SunMax;
      Reflect = 0.0;
    }
    else {
      /*  sloping surface  */
      DiffuseSkyView=(PI - CellSlope) / PI;
      ReflectedSkyView = CellSlope / PI;
      SolarAltitude = asin(SineSolarAltitude);
      
      CosineIncidenceAngle = cos(SolarAltitude) * sin(CellSlope) *
        cos(SolarAzimuth-CellAspect) + cos(CellSlope) * sin(SolarAltitude);
      
      *Direct = SunMax * CosineIncidenceAngle / SineSolarAltitude;
      if (CosineIncidenceAngle <= 0.0) 
        *Direct = 0.0;
      
      *Diffuse = SunMax * DiffuseSkyView;
      Reflect = Albedo * SunMax * ReflectedSkyView;
    } /* end sloping surface*/
  } /* end daylight is true*/
  else {   
    /* sun is below the horizon */
    
    /* all measured solar radiation is diffuse, thus total rad on CellSlope 
       is (measured)*DiffuseSkyView */
    *Direct = 0.0;
    Reflect = 0.0;
    *Diffuse = (PI - CellSlope) / PI;
  }
  
  /* Average over timestep */    
  if (Dt > 0.0) {
    
    /* The assumption is made that the incoming reflected radiation is diffuse 
       radiation for the "receiving" gridcell */
    
    *Diffuse += Reflect;
    *Direct *= SolarTimeStep / (Dt * RADPHOUR);
    *Diffuse *= SolarTimeStep / (Dt * RADPHOUR);
  }  
}

/*****************************************************************************
  Function name: SolarConst()  
  
  Purpose      : this subroutine reads in site data, calculates constants,
                 and converts data from degrees to radians.

  Required     :
    float lat_deg           - latitude of site in degrees
    float lat_min           - latitude of site in minutes
    float lng_deg           - longitude of site in degrees
    float lng_min           - longitude of site in minutes
    float *StandardMeridian - standard meridian for local time zone in 
                              degrees
  
  Returns      : void

  Modifies     :
    float *Latitude         - latitude of site in radians
    float *Longitude        - longitude of site in radians
    float *StandardMeridian - standard meridian for local time zone (rads)

  Comments     : EXECUTE ONCE BEFORE TIME LOOP
*****************************************************************************/
void SolarConst(float lat_deg, float lat_min, float lng_deg, 
                float lng_min, float *StandardMeridian, float *Latitude, 
                float *Longitude)
{
  /*  convert to radians  */      
  *Latitude=(lat_deg+(lat_min/60.))*PI/180.;
  *Longitude=(lng_deg+(lng_min/60.))*PI/180.;
  *StandardMeridian=*StandardMeridian*PI/180.;
}

/*****************************************************************************
  Function name: CalcSolarDay()

  Purpose      : Initialize the Earth-Sun geometry variables at the beginning
                 of each day
*****************************************************************************/
void CalcSolarDay(int DayOfYear, SOLARGEOMETRY *SolarGeo)
{
  const char *Routine = "InitNewDay";

  SolarDay(DayOfYear, SolarGeo->Longitude, SolarGeo->Latitude,
           SolarGeo->StandardMeridian, &(SolarGeo->NoonHour), 
           &(SolarGeo->Declination), &(SolarGeo->HalfDayLength),
           &(SolarGeo->Sunrise), &(SolarGeo->Sunset),
           &(SolarGeo->TimeAdjustment), &(SolarGeo->SunEarthDistance));
}


/*****************************************************************************
  Function name: CalcSolarHour()

  Purpose      : Calculate variables related to the position of the sun above the
                 horizon
*****************************************************************************/
void CalcSolarHour(float StartHour, float Dt, SOLARGEOMETRY *SolarGeo)
{

  SolarHour(SolarGeo->Latitude, StartHour + Dt, Dt, 
            SolarGeo->NoonHour, SolarGeo->Declination, SolarGeo->Sunrise, 
            SolarGeo->Sunset, SolarGeo->TimeAdjustment, 
            SolarGeo->SunEarthDistance, &(SolarGeo->SineSolarAltitude), 
            &(SolarGeo->DayLight), &(SolarGeo->SolarTimeStep), 
            &(SolarGeo->SunMax), &(SolarGeo->SolarAzimuth));
}

