testsunfun

  Calculation of local times of sunrise, solar noon, and sunset
  based on the calculation procedure by NOAA in the javascript in
  http : // www.srrb.noaa.gov/highlights/sunrise/sunrise.html and
  http : // www.srrb.noaa.gov/highlights/sunrise/azel.html

 Five functions are available for use from Excel worksheets :
  -sunrise (lat, lon, year, month, day, timezone, dlstime)
 -solarnoon (lat, lon, year, month, day, timezone, dlstime)
 -sunset (lat, lon, year, month, day, timezone, dlstime)
 -solarazimuth (lat, lon, year, month, day, hour, minute, second, timezone, dlstime)
 -solarelevation (lat, lon, year, month, day, hour, minute, second, timezone, dlstime)

 The sign convention for inputs to the functions named sunrise, solarnoon, sunset, solarazimuth, and solarelevationis :
  -positive latitude decimal degrees for northern hemisphere
 -negative longitude degrees for western hemisphere
 -negative time zone hours for western hemisphere

 The other functions in the VBA module use the original
  NOAA sign convention of positive longitude in the western hemisphere. The calculations in the NOAA Sunrise/Sunset and Solar Position
  Calculators are based on equations from Astronomical Algorithms, by Jean Meeus.NOAA also included atmospheric refraction effects.
  The sunrise and sunset results were reported by NOAA
  to be accurate to within + /-1 minute for locations between + /-72 \[Degree]
  latitude, and within ten minutes outside of those latitudes. This translation was tested for selected locations
  and found to provide results within + /-1 minute of the
  original Javascript code. This translation does not include calculation of prior or next
  susets for locations above the Arctic Circle and below the
  Antarctic Circle, when a sunrise or sunset does not occur. Translated from NOAAs Javascript to Excel VBA by : Greg Pelletier
  Department of Ecology
  P.O.Box 47600
  Olympia, WA 98504 - 7600
  email : gpel461@ecy.wa.gov


Name : calcJD
'*Type : Function
'*Purpose : Julian day from calendar day
    '*Arguments : '*year : 4 digit year
'*month : January = 1
'*day : 1 - 31
'*Return value : '*The Julian day corresponding to the date
    '*Note : '*Number is returned for start of day.Fractional days should be
'*added later.' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/

BeginPackage["Sunpath`"];
dawn::usage =
  "dawn[lat_,lon_,year_,month_,day_,timezone_,dlstime_,solardepression_] \
gives the time of the dawn for various solar depressions (6 deg is civil \
dawn, 12 deg is nautical dawn, 18 deg is astronomical dawn). The number \
returned is the time in a single number ranging from 0 to 24; hours, minutes \
and seconds can be derived from this. Lat and long coordinates are in \
degrees, with directions east and north counted positive. Timezone is the \
offset with respect to GMT. dlstime indicates whether daylight savings time \
is in effect (0=false, 1=true).";
sunrise::usage =
  "sunrise[lat_,lon_,year_,month_,day_,timezone_,dlstime_] gives the time of \
the sunrise taking into account atmospheric refraction. The number returned \
is the time in days; hours, minutes and seconds can be derived from this. Lat \
and long coordinates are in degrees, with directions east and north counted \
positive. Timezone is the offset with respect to GMT. dlstime indicates \
whether daylight savings time is in effect (0=false, 1=true).";
solarnoon::usage =
  "solarnoon[lat_,lon_,year_,month_,day_,timezone_,dlstime_] gives the time \
at which the sun reaches its highest elevation. The number returned is the \
time in days; hours, minutes and seconds can be derived from this. Lat and \
long coordinates are in degrees, with directions east and north counted \
positive. Timezone is the offset with respect to GMT. dlstime indicates \
whether daylight savings time is in effect (0=false, 1=true).";
sunset::usage =
  "sunset[lat_,lon_,year_,month_,day_,timezone_,dlstime_] gives the time of \
the sunset taking into account atmospheric refraction. The number returned is \
the time in days; hours, minutes and seconds can be derived from this. Lat \
and long coordinates are in degrees, with directions east and north counted \
positive. Timezone is the offset with respect to GMT. dlstime indicates \
whether daylight savings time is in effect (0=false, 1=true).";
dusk::usage =
  "dusk[lat_,lon_,year_,month_,day_,timezone_,dlstime_,solardepression_] \
gives the time of the sunset taking into account atmospheric refraction.The \
number returned is the time in a single number ranging from 0 to 24; hours, \
minutes and seconds can be derived from this. Lat and long coordinates are in \
degrees, with directions east and north counted positive. Timezone is the \
offset with respect to GMT. dlstime indicates whether daylight savings time \
is in effect (0=false, 1=true).";
solarposition::usage =
  "solarposition[lat_,lon_,year_,month_,day_,hours_,minutes_,seconds_,\
timezone_,dlstime_] gives the time of the sunset taking into account \
atmospheric refraction. Lat and long coordinates are in degrees, with \
directions east and north counted positive. Timezone is the offset with \
respect to GMT. dlstime indicates whether daylight savings time is in effect \
(0=false, 1=true). The function returns the sumposition as a list consisting \
of the azimuth and elevation is degrees";
Begin["`Private`"]


calcJD[year0_, month0_, day0_] :=
  Module[{year = year0, month = month0, day = day0, A, B},
   If [month <= 2, year = year - 1; month = month + 12];
   A = Floor[year/100];
   B = 2 - A + Floor[A/4];
   Floor[365.25*(year + 4716)] + Floor[30.6001*(month + 1)] + day + B - 1524.5
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***
/'*Name : calcTimeJulianCent
'*Type : Function
'*Purpose : convert Julian Day to centuries since J2000 .0.'*Arguments : '*jd : the Julian Day to convert
        '*Return value : '*the T value corresponding to the Julian Day
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim t As Double

calcTimeJulianCent[JD_] := (JD - 2451545)/36525;

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*
Name : calcJDFromJulianCent
'*Type : Function
'*Purpose : convert centuries since J2000 .0 to Julian Day.'*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*the Julian Day corresponding to the t value
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim JD As Double

calcJDFromJulianCent[t_] := t*36525 + 2451545;

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/
'*Name : calGeomMeanLongSun
'*Type : Function
'*Purpose : calculate the Geometric Mean Longitude of the Sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*the Geometric Mean Longitude of the Sun in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim l0 As Double

calcGeomMeanLongSun[t_] :=
  Mod[280.46646 + t*(36000.76983 + 0.0003032*t), 360];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*
Name : calGeomAnomalySun
'*Type : Function
'*Purpose : calculate the Geometric Mean Anomaly of the Sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*the Geometric Mean Anomaly of the Sun in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim m As Double

calcGeomMeanAnomalySun[t_] := 357.52911 + t*(35999.05029 - 0.0001537*t);

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'
*Name : calcEccentricityEarthOrbit
'*Type : Function
'*Purpose : calculate the eccentricity of earth' s orbit
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*the unitless eccentricity
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim e As Double

calcEccentricityEarthOrbit[t_] :=
  0.016708634 - t*(0.000042037 + 0.0000001267*t);

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'
*Name : calcSunEqOfCenter
'*Type : Function
'*Purpose : calculate the equation of center for the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim m As Double, mrad As Double, sinm As Double, sin2m As Double, sin3m As Double
Dim c As Double

calcSunEqOfCenter[t_] :=
  Module[{m},
   m = calcGeomMeanAnomalySun[t];
   Sin[m \[Degree]]*(1.914602 - t*(0.004817 + 0.000014*t)) +
    Sin[2 m \[Degree]]*(0.019993 - 0.000101*t) + Sin[3 m \[Degree]]*0.000289
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunTrueLong
'*Type : Function
'*Purpose : calculate the true longitude of the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun' s true longitude in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim l0 As Double, c As Double, O As Double

calcSunTrueLong[t_] := calcGeomMeanLongSun[t] + calcSunEqOfCenter[t];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunTrueAnomaly (not used by sunrise, solarnoon, sunset)
    '*Type : Function
'*Purpose : calculate the true anamoly of the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun' s true anamoly in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim m As Double, c As Double, v As Double

calcSunTrueAnomaly[t_] := calcGeomMeanAnomalySun[t] + calcSunEqOfCenter[t];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunRadVector (not used by sunrise, solarnoon, sunset)
    '*Type : Function
'*Purpose : calculate the distance to the sun in AU
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun radius vector in AUs
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim v As Double, e As Double, R As Double

calcSunRadVector[
   t_] := (1.000001018*(1 - calcEccentricityEarthOrbit[t]^2))/(1 +
     calcEccentricityEarthOrbit[t]*Cos[calcSunTrueAnomaly[t] \[Degree]]);

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunApparentLong (not used by sunrise, solarnoon, sunset)
    '*Type : Function
'*Purpose : calculate the apparent longitude of the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun' s apparent longitude in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim O As Double, omega As Double, lambda As Double

calcSunApparentLong[t_] :=
  calcSunTrueLong[t] - 0.00569 - 0.00478*Sin[(125.04 - 1934.136*t) \[Degree]];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcMeanObliquityOfEcliptic
'*Type : Function
'*Purpose : calculate the mean obliquity of the ecliptic
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*mean obliquity in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim seconds As Double, e0 As Double

calcMeanObliquityOfEcliptic[t_] :=
  Module[{seconds},
   seconds = 21.448 - t*(46.815 + t*(0.00059 - t*(0.001813)));
   23 + (26 + (seconds/60))/60
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcObliquityCorrection
'*Type : Function
'*Purpose : calculate the corrected obliquity of the ecliptic
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*corrected obliquity in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim e0 As Double, omega As Double, e As Double

calcObliquityCorrection[t_] :=
  Module[{e0, omega},
   e0 = calcMeanObliquityOfEcliptic[t];
   omega = 125.04 - 1934.136*t;
   e0 + 0.00256*Cos[omega \[Degree]]
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunRtAscension (not used by sunrise, solarnoon, sunset)
    '*Type : Function
'*Purpose : calculate the right ascension of the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun' s right ascension in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim e As Double, lambda As Double, tananum As Double, tanadenom As Double
Dim alpha As Double

calcSunRtAscension[t_] :=
  Module[{e, lambda, tananum, tanadenom},
   e = calcObliquityCorrection[t];
   lambda = calcSunApparentLong[t];
   tananum = (Cos[e \[Degree]]*Sin[lambda \[Degree]]);
   tanadenom = (Cos[lambda \[Degree]]);
   ArcTan[tanadenom, tananum]/\[Degree]
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunDeclination
'*Type : Function
'*Purpose : calculate the declination of the sun
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*sun' s declination in degrees
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim e As Double, lambda As Double, sint As Double, theta As Double

calcSunDeclination[t_] :=
  Module[{e, lambda},
   e = calcObliquityCorrection[t];
   lambda = calcSunApparentLong[t];
   ArcSin[Sin[e \[Degree]]*Sin[lambda \[Degree]]]/\[Degree]
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcEquationOfTime
'*Type : Function
'*Purpose : calculate the difference between true solar time and mean
    '*solar time
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*Return value : '*equation of time in minutes of time
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim epsilon As Double, l0 As Double, e As Double, m As Double
Dim y As Double, sin2l0 As Double, sinm As Double
Dim cos2l0 As Double, sin4l0 As Double, sin2m As Double, Etime As Double

calcEquationOfTime[t_] :=
  Module[{epsilon, l0, e, m, y, sinm, sin2l0, cos2l0, sin4l0, sin2m, Etime},
   epsilon = calcObliquityCorrection[t];
   l0 = calcGeomMeanLongSun[t];
   e = calcEccentricityEarthOrbit[t];
   m = calcGeomMeanAnomalySun[t];
   y = Tan[epsilon \[Degree]/2]^2;
   sin2l0 = Sin[2*l0 \[Degree]];
   sinm = Sin[m \[Degree]];
   cos2l0 = Cos[2*l0 \[Degree]];
   sin4l0 = Sin[4*l0 \[Degree]];
   sin2m = Sin[2*m \[Degree]];
   Etime =
    y*sin2l0 - 2*e*sinm + 4*e*y*sinm*cos2l0 - 0.5*y*y*sin4l0 - 1.25*e*e*sin2m;
   Etime/\[Degree]*4
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcHourAngleDawn
'*Type : Function
'*Purpose : calculate the hour angle of the sun at dawn for the
    '*latitude
    '*for user selected solar depression below horizon
    '*Arguments : '*lat : latitude of observer in degrees
        '*solarDec : declination angle of sun in degrees
          '*solardepression : angle of the sun below the horizion in degrees
            '*Return value : '*hour angle of dawn in radians
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

calcHourAngleDawn[lat_, solarDec_, solardepression_] :=
  ArcCos[Cos[(90 + solardepression) \[Degree]] /(Cos[lat \[Degree]]*
       Cos[solarDec \[Degree]]) - Tan[lat \[Degree]]*Tan[solarDec \[Degree]]];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcHourAngleSunrise
'*Type : Function
'*Purpose : calculate the hour angle of the sun at sunrise for the
    '*latitude
    '*Arguments : '*lat : latitude of observer in degrees
        '*solarDec : declination angle of sun in degrees
          '*Return value : '*hour angle of sunrise in radians
'*'*Note : For sunrise and sunset calculations, we assume 0.833 \[Degree] of atmospheric refraction
'*For details about refraction see http : // www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
'*' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

calcHourAngleSunrise[lat_, solarDec_] :=
  ArcCos[Cos[90.833 \[Degree]]/(Cos[lat \[Degree]]*Cos[solarDec \[Degree]]) -
    Tan[lat \[Degree]]*Tan[solarDec \[Degree]]];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcHourAngleSunset
'*Type : Function
'*Purpose : calculate the hour angle of the sun at sunset for the
    '*latitude
    '*Arguments : '*lat : latitude of observer in degrees
        '*solarDec : declination angle of sun in degrees
          '*Return value : '*hour angle of sunset in radians
'*'*Note : For sunrise and sunset calculations, we assume 0.833 \[Degree] of atmospheric refraction
'*For details about refraction see http : // www.srrb.noaa.gov/highlights/sunrise/calcdetails.html
'*' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

calcHourAngleSunset[lat_,
   solarDec_] := -ArcCos[
    Cos[90.833 \[Degree]]/(Cos[lat \[Degree]]*Cos[solarDec \[Degree]]) -
     Tan[lat \[Degree]]*Tan[solarDec \[Degree]]];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcHourAngleDusk
'*Type : Function
'*Purpose : calculate the hour angle of the sun at dusk for the
    '*latitude
    '*for user selected solar depression below horizon
    '*Arguments : '*lat : latitude of observer in degrees
        '*solarDec : declination angle of sun in degrees
          '*solardepression : angle of sun below horizon in degrees
            '*Return value : '*hour angle of dusk in radians
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double


calcHourAngleDusk[lat_, solarDec_,
   solardepression_] := -ArcCos[
    Cos[(90 + solardepression) \[Degree]]/(Cos[lat \[Degree]]*
        Cos[solarDec \[Degree]]) -
     Tan[lat \[Degree]]*Tan[solarDec \[Degree]]];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcDawnUTC
'*Type : Function
'*Purpose : calculate the Universal Coordinated Time (UTC) of dawn
    '*for the given day at the given location on earth
    '*for user selected solar depression below horizon
    '*Arguments : '*JD : julian day
        '*latitude : latitude of observer in degrees
          '*longitude : longitude of observer in degrees
            '*solardepression : angle of sun below the horizon in degrees
              '*Return value : '*time in minutes from zero Z
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double

calcDawnUTC[JD_, latitude_, longitude_, solardepression_] :=
  Module[{t, eqtime, solarDec, hourangle, delta, timeDiff, timeUTC, newt},
   t = calcTimeJulianCent[JD];

   (* First pass to approximate sunrise *)
   eqtime = calcEquationOfTime[t];
   solarDec = calcSunDeclination[t];
   hourangle = calcHourAngleSunrise[latitude, solarDec];

   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   (* in minutes of time *)
   timeUTC = 720 + timeDiff - eqtime;
   (* in minutes *)

   (*Second pass includes fractional jday in gamma calc *)

   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + timeUTC/1440];
   eqtime = calcEquationOfTime[newt];
   solarDec = calcSunDeclination[newt];
   hourangle = calcHourAngleDawn[latitude, solarDec, solardepression];
   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime(* in minutes *)
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunriseUTC
'*Type : Function
'*Purpose : calculate the Universal Coordinated Time (UTC) of sunrise
    '*for the given day at the given location on earth
    '*Arguments : '*JD : julian day
        '*latitude : latitude of observer in degrees
          '*longitude : longitude of observer in degrees
            '*Return value : '*time in minutes from zero Z
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double

calcSunriseUTC[JD_, latitude_, longitude_] :=
  Module[{t, eqtime, solarDec, hourangle, delta, timeDiff, timeUTC, newt},
   t = calcTimeJulianCent[JD];
   (* First pass to approximate sunrise *)
   eqtime = calcEquationOfTime[t];
   solarDec = calcSunDeclination[t];
   hourangle = calcHourAngleSunrise[latitude, solarDec];
   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   (* in minutes of time *)
   timeUTC = 720 + timeDiff - eqtime;
   (* in minutes *)

   (* Second pass includes fractional jday in gamma calc *)

   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + timeUTC/1440];
   eqtime = calcEquationOfTime[newt];
   solarDec = calcSunDeclination[newt];
   hourangle = calcHourAngleSunrise[latitude, solarDec];
   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime
   (* in minutes *)
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSolNoonUTC
'*Type : Function
'*Purpose : calculate the Universal Coordinated Time (UTC) of solar
    '*noon for the given day at the given location on earth
    '*Arguments : '*t : number of Julian centuries since J2000 .0
        '*longitude : longitude of observer in degrees
          '*Return value : '*time in minutes from zero Z
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim newt As Double, eqtime As Double, solarNoonDec As Double, solNoonUTC As Double

calcSolNoonUTC[t_, longitude_] :=
  Module[{newt, eqtime, solarNoonDec},
   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + 0.5 + longitude/360];
   eqtime = calcEquationOfTime[newt];
   (*solarNoonDec=calcSunDeclination[newt];*)
   720 + (longitude*4) - eqtime
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcSunsetUTC
'*Type : Function
'*Purpose : calculate the Universal Coordinated Time (UTC) of sunset
    '*for the given day at the given location on earth
    '*Arguments : '*JD : julian day
        '*latitude : latitude of observer in degrees
          '*longitude : longitude of observer in degrees
            '*Return value : '*time in minutes from zero Z
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double

calcSunsetUTC[JD_, latitude_, longitude_] :=
  Module[{t, eqtime, solarDec, hourangle, delta, timeDiff, timeUTC, newt},
   t = calcTimeJulianCent[JD];

   (*First calculates sunrise and approx length of day*)

   eqtime = calcEquationOfTime[t];
   solarDec = calcSunDeclination[t];
   hourangle = calcHourAngleSunset[latitude, solarDec];

   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime;

   (* first pass used to include fractional day in gamma calc *)

   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + timeUTC/1440];
   eqtime = calcEquationOfTime[newt];
   solarDec = calcSunDeclination[newt];
   hourangle = calcHourAngleSunset[latitude, solarDec];

   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime
   (* in minutes *)
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : calcDuskUTC
'*Type : Function
'*Purpose : calculate the Universal Coordinated Time (UTC) of dusk
'*for the given day at the given location on earth
'*for user selected solar depression below horizon
'*Arguments : '*JD : julian day
'*latitude : latitude of observer in degrees
'*longitude : longitude of observer in degrees
'*solardepression : angle of sun below horizon
'*Return value : '*time in minutes from zero Z
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
Dim delta As Double, timeDiff As Double, timeUTC As Double
Dim newt As Double

calcDuskUTC[JD_, latitude_, longitude_, solardepression_] :=
  Module[{t, eqtime, solarDec, hourangle, delta, timeDiff, timeUTC, newt},
   t = calcTimeJulianCent[JD];

   (* First calculates sunrise and approx length of day *)

   eqtime = calcEquationOfTime[t];
   solarDec = calcSunDeclination[t];
   hourangle = calcHourAngleSunset[latitude, solarDec];

   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime;

   (*first pass used to include fractional day in gamma calc *)

   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + timeUTC/1440];
   eqtime = calcEquationOfTime[newt];
   solarDec = calcSunDeclination[newt];
   hourangle = calcHourAngleDusk[latitude, solarDec, solardepression];

   delta = longitude - hourangle/\[Degree];
   timeDiff = 4*delta;
   timeUTC = 720 + timeDiff - eqtime
   (*in minutes*)
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : dawn
'*Type : Main Function called by spreadsheet
    '*Purpose : calculate time of dawn for the entered date
'*and location.'*For latitudes greater than 72 degrees N and S, calculations are
'*accurate to within 10 minutes.For latitudes less than + /-72 \[Degree]
'*accuracy is approximately one minute.'*Arguments : ' latitude = latitude (decimal degrees)
' longitude = longitude (decimal degrees)
   ' NOTE : longitude is negative for western hemisphere for input cells
     ' in the spreadsheet for calls to the functions named
' sunrise, solarnoon, and sunset.Those functions convert the
  ' longitude to positive for the western hemisphere for calls to
  ' other functions using the original sign convention
' from the NOAA javascript code.' year = year
' month = month
' day = day
' timezone = time zone hours relative to GMT/UTC (hours)
' dlstime = daylight savings time (0 = no, 1 = yes) (hours)
' solardepression = angle of sun below horizon in degrees
   '*Return value : '*dawn time in local time (days)
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double, JD As Double
Dim riseTimeGMT As Double, riseTimeLST As Double

dawn[lat_, lon_, year_, month_, day_, timezone_, dlstime_, solardepression_] :=


  Module[{longitude, latitude, JD, riseTimeGMT, riseTimeLST},
   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   JD = calcJD[year, month, day];

   (* Calculate sunrise for this date *)
   riseTimeGMT = calcDawnUTC[JD, latitude, longitude, solardepression];

   (* adjust for time zone and daylight savings time in minutes *)
   riseTimeLST = riseTimeGMT + (60*timezone) + (dlstime*60);

   (* convert to days *)
   riseTimeLST/1440
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : sunrise
'*Type : Main Function called by spreadsheet
    '*Purpose : calculate time of sunrise for the entered date
'*and location.'*For latitudes greater than 72 degrees N and S, calculations are
'*accurate to within 10 minutes.For latitudes less than + /-72 \[Degree]
'*accuracy is approximately one minute.'*Arguments : ' latitude = latitude (decimal degrees)
' longitude = longitude (decimal degrees)
   ' NOTE : longitude is negative for western hemisphere for input cells
     ' in the spreadsheet for calls to the functions named
' sunrise, solarnoon, and sunset.Those functions convert the
  ' longitude to positive for the western hemisphere for calls to
  ' other functions using the original sign convention
' from the NOAA javascript code.' year = year
' month = month
' day = day
' timezone = time zone hours relative to GMT/UTC (hours)
' dlstime = daylight savings time (0 = no, 1 = yes) (hours)
   '*Return value : '*sunrise time in local time (days)
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double, JD As Double
Dim riseTimeGMT As Double, riseTimeLST As Double

sunrise[lat_, lon_, year_, month_, day_, timezone_, dlstime_] :=
  Module[{longitude, latitude, JD, riseTimeGMT, riseTimeLST},
   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   JD = calcJD[year, month, day];

   (* Calculate sunrise for this date *)
   riseTimeGMT = calcSunriseUTC[JD, latitude, longitude];

   (* adjust for time zone and daylight savings time in minutes *)
   riseTimeLST = riseTimeGMT + (60*timezone) + (dlstime*60);

   (* convert to days *)
   riseTimeLST/1440
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : solarnoon
'*Type : Main Function called by spreadsheet
    '*Purpose : calculate the Universal Coordinated Time (UTC) of solar
      '*noon for the given day at the given location on earth
      '*Arguments : ' year
        ' month
        ' day
        '*longitude : longitude of observer in degrees
          ' NOTE : longitude is negative for western hemisphere for input cells
            ' in the spreadsheet for calls to the functions named
' sunrise, solarnoon, and sunset.Those functions convert the
  ' longitude to positive for the western hemisphere for calls to
  ' other functions using the original sign convention
  ' from the NOAA javascript code.'*Return value : '*time of solar noon in local time days
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double, JD As Double
Dim t As Double, newt As Double, eqtime As Double
Dim solarNoonDec As Double, solNoonUTC As Double

solarnoon[lat_, lon_, year_, month_, day_, timezone_, dlstime_] :=
  Module[{longitude, latitude, JD, t, newt, eqtime, solarNoonDec, solNoonUTC,
    solarnoon},
   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   JD = calcJD[year, month, day];

   t = calcTimeJulianCent[JD];

   newt = calcTimeJulianCent[calcJDFromJulianCent[t] + 0.5 + longitude/360];

   eqtime = calcEquationOfTime[newt];
   (* solarNoonDec=calcSunDeclination[newt]; *)
   solNoonUTC = 720 + (longitude*4) - eqtime;

   (* adjust for time zone and daylight savings time in minutes *)
   solarnoon = solNoonUTC + (60*timezone) + (dlstime*60);

   (* convert to days *)
   solarnoon/1440
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : sunset
'*Type : Main Function called by spreadsheet
    '*Purpose : calculate time of sunrise and sunset for the entered date
'*and location.'*For latitudes greater than 72 degrees N and S, calculations are
'*accurate to within 10 minutes.For latitudes less than + /-72 \[Degree]
'*accuracy is approximately one minute.'*Arguments : ' latitude = latitude (decimal degrees)
' longitude = longitude (decimal degrees)
   ' NOTE : longitude is negative for western hemisphere for input cells
     ' in the spreadsheet for calls to the functions named
' sunrise, solarnoon, and sunset.Those functions convert the
  ' longitude to positive for the western hemisphere for calls to
  ' other functions using the original sign convention
' from the NOAA javascript code.' year = year
' month = month
' day = day
' timezone = time zone hours relative to GMT/UTC (hours)
' dlstime = daylight savings time (0 = no, 1 = yes) (hours)
   '*Return value : '*sunset time in local time (days)
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double, JD As Double
Dim setTimeGMT As Double, setTimeLST As Double

sunset[lat_, lon_, year_, month_, day_, timezone_, dlstime_] :=
  Module[{longitude, latitude, JD, setTimeGMT, setTimeLST},
   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   JD = calcJD[year, month, day];

   (* Calculate sunset for this date *)
   setTimeGMT = calcSunsetUTC[JD, latitude, longitude];

   (* adjust for time zone and daylight savings time in minutes *)
   setTimeLST = setTimeGMT + (60*timezone) + (dlstime*60);

   (* convert to days *)
   setTimeLST/1440
   ];


' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : dusk
'*Type : Main Function called by spreadsheet
'*Purpose : calculate time of sunrise and sunset for the entered date
'*and location.'*For latitudes greater than 72 degrees N and S, calculations are
'*accurate to within 10 minutes.For latitudes less than + /-72 \[Degree]
'*accuracy is approximately one minute.'*Arguments : ' latitude = latitude (decimal degrees)
' longitude = longitude (decimal degrees)
' NOTE : longitude is negative for western hemisphere for input cells
' in the spreadsheet for calls to the functions named
' sunrise, solarnoon, and sunset.Those functions convert the
' longitude to positive for the western hemisphere for calls to
' other functions using the original sign convention
' from the NOAA javascript code.' year = year
' month = month
' day = day
' timezone = time zone hours relative to GMT/UTC (hours)
' dlstime = daylight savings time (0 = no, 1 = yes) (hours)
' solardepression = angle of sun below horizon in degrees
'*Return value : '*dusk time in local time (days)
' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double, JD As Double
Dim setTimeGMT As Double, setTimeLST As Double


dusk[lat_, lon_, year_, month_, day_, timezone_, dlstime_, solardepression_] :=


  Module[{longitude, latitude, JD, setTimeGMT, setTimeLST},

   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   JD = calcJD[year, month, day];

   (* Calculate sunset for this date *)
   setTimeGMT = calcDuskUTC[JD, latitude, longitude, solardepression];

   (* adjust for time zone and daylight savings time in minutes *)
   setTimeLST = setTimeGMT + (60*timezone) + (dlstime*60);

   (* convert to days *)
   setTimeLST/1440
   ];

' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/'*Name : solarazimuth
'*Type : Main Function
'*Purpose : calculate solar azimuth (deg from north) for the entered
'*date, time and location.Returns - 999999 if darker than twilight
'*'*Arguments : '*latitude, longitude, year, month, day, hour, minute, second, '*timezone, daylightsavingstime
'*Return value : '*solar azimuth in degrees from north
'*'*Note : solarelevation and solarazimuth functions are identical
'*and could converted to a VBA subroutine that would return
'*both values.'*' ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ***/Dim longitude As Double, latitude As Double
Dim zone As Double, daySavings As Double
Dim hh As Double, mm As Double, ss As Double, timenow As Double
Dim JD As Double, t As Double, R As Double
Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double
Dim solarDec As Double, earthRadVec As Double, solarTimeFix As Double
Dim trueSolarTime As Double, hourangle As Double, harad As Double
Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double
Dim azimuth As Double, exoatmElevation As Double
Dim step1 As Double, step2 As Double, step3 As Double
Dim refractionCorrection As Double, te As Double, solarzen As Double

solarposition[lat_, lon_, year_, month_, day_, hours_, minutes_, seconds_,
   timezone_, dlstime_] :=
  Module[{longitude, latitude, zone, daySavings, hh, mm, ss, timenow, JD, t,
    R, alpha, theta, Etime, eqtime, solarDec, earthRadVec, solarTimeFix,
    trueSolarTime, hourangle, harad, csz, zenith, azDenom, azRad, azimuth,
    exoatmElevation, refractionCorrection, te, step1, step2, step3, solarzen,
    solarazimuth, solarelevation},
   (* change sign convention for longitude from negative to positive in \
western hemisphere *)
   longitude = -lon;
   latitude = lat;
   If [latitude > 89.8, latitude = 89.8];
   If [latitude < -89.8, latitude = -89.8];

   (* change time zone to ppositive hours in western hemisphere *)
   zone = -timezone;
   daySavings = dlstime*60;
   hh = hours - (daySavings/60);
   mm = minutes;
   ss = seconds;

   (* timenow is GMT time for calculation in hours since 0Z *)
   timenow = hh + mm/60 + ss/3600 + zone;

   JD = calcJD[year, month, day];
   t = calcTimeJulianCent[JD + timenow/24];
   R = calcSunRadVector[t];
   alpha = calcSunRtAscension[t];
   theta = calcSunDeclination[t];
   Etime = calcEquationOfTime[t];

   eqtime = Etime;
   solarDec = theta; (* in degrees *)
   earthRadVec = R;

   solarTimeFix = eqtime - 4*longitude + 60*zone;
   trueSolarTime = hh*60 + mm + ss/60 + solarTimeFix;
   (* in minutes *)

   While [trueSolarTime > 1440, trueSolarTime = trueSolarTime - 1440];

   hourangle = trueSolarTime/4 - 180;
   (* Thanks to Louis Schwarzmayr for the next line: *)
   If[ (hourangle < -180) , hourangle = hourangle + 360];

   harad = hourangle \[Degree];

   csz = Sin[latitude \[Degree]]*Sin[solarDec \[Degree]] +
     Cos[latitude \[Degree]]*Cos[solarDec \[Degree]]*Cos[harad];

   If[csz > 1 ,
    csz = 1,
    If [csz < -1 , csz = -1]
    ];

   zenith = ArcCos[csz]/\[Degree];

   azDenom = Cos[latitude \[Degree]]*Sin[zenith \[Degree]];

   If [Abs[azDenom] > 0.001,
    azRad = (Sin[latitude \[Degree]]*Cos[zenith \[Degree]] -
        Sin[solarDec \[Degree]])/azDenom;
    If[Abs[azRad] > 1,
     If [azRad < 0, azRad = -1, azRad = 1]
     ];
    azimuth = 180 - ArcCos[azRad]/\[Degree];
    If[hourangle > 0, azimuth = -azimuth],
    (*else*)
    If[latitude > 0, azimuth = 180, azimuth = 0]
    ];
   If [azimuth < 0, azimuth = azimuth + 360];

   exoatmElevation = 90 - zenith;

   (* beginning of simplified expression *)
   If [exoatmElevation > 85,
    refractionCorrection = 0,
    te = Tan[exoatmElevation \[Degree]];
    If[exoatmElevation > 5,
     refractionCorrection = 58.1/te - 0.07/te^3 + 0.000086/te^5,
     (*else*)
     If [exoatmElevation > -0.575,
      step1 = (-12.79 + exoatmElevation*0.711);
      step2 = (103.4 + exoatmElevation*(step1));
      step3 = (-518.2 + exoatmElevation*(step2));
      refractionCorrection = 1735 + exoatmElevation*(step3),
      (*else*)
      refractionCorrection = -20.774/te
      ]
     ];
    refractionCorrection = refractionCorrection/3600;
    ];
   (* end of simplified expression *)

   solarzen = zenith - refractionCorrection;
   solarazimuth = azimuth;
   solarelevation = 90 - solarzen;
   Return[{solarazimuth, solarelevation}];
   ];

End[ ]
EndPackage[ ]