DEFINITION MODULE SunPath;
(*******************************************************************
Module SunPath (Version 1.1)
Copyright (c) 1989-2006 by Thomas Nemecek and ETH Zurich.
Purpose Compute solar angles, daylength,
real solar time, etc.
Remarks Units:
angles: radians
time: fraction of day (24h = 1.0)
except where the variable name
"hours" is used
Initialization:
the module is initialised for julian day = 1;
position, zone time = Zürich, Switzerland;
winter time.
WARNING: Implementation is such that computations
are done with limited accuracy.
References:
Meteonorm - Theorie für den Solarplaner.
BA f. Energiewirtschaft, 3003 Bern; 1985.
France, J. & Thornley, J.H.M (1984):
Mathematical models in agriculture, Butterworths,
London.
Programming
o Design
Thomas Nemecek 24/10/1989
Andreas Fischlin 13/08/2007
o Implementation
Thomas Nemecek 25/10/1989
Andreas Fischlin 13/08/2007
ETH Zurich
Systems Ecology
CHN E 35.1
Universitaetstrasse 16
8092 Zurich
SWITZERLAND
URLs:
<mailto:RAMSES@env.ethz.ch>
<http://www.sysecol.ethz.ch>
<http://www.sysecol.ethz.ch/SimSoftware/RAMSES>
Last revision of definition: 13/08/2007 af
*******************************************************************)
(***************************************)
(*##### Angle transformations #####*)
(***************************************)
PROCEDURE Rad(deg: REAL): REAL;
PROCEDURE Deg(rad: REAL): REAL;
(**************************************)
(*##### Time transformations #####*)
(**************************************)
PROCEDURE DayFraction(hours, minutes, seconds: INTEGER): REAL;
PROCEDURE DayFractionToTime(dayFraction: REAL; VAR hours, minutes, seconds: INTEGER);
PROCEDURE FractionOfYear (dayOfYear: INTEGER; dayfraction: REAL): REAL;
(*
Returns fraction of year where dayOfYear = 1 ~ 1st January;
dayOfYear = 365 (not leap year) or 366 (leap year) ~ 31st
December and dayfraction is a number in [0.0..1.0].
*)
(***************************************)
(*##### Set location and time #####*)
(***************************************)
(*----------------------------*)
(*===== Set location =====*)
(*----------------------------*)
PROCEDURE SetGeogrLocation( latitude, longitude: REAL);
PROCEDURE GetGeogrLocation(VAR latitude, longitude: REAL);
(* latitude and longitude are to be defined as degrees and fractions
thereof.
Latitude: N ==> >0 S ==> <0
Longitude: W ==> >0 E ==> <0
*)
PROCEDURE SetTimeZone( hours: REAL);
PROCEDURE GetTimeZone(VAR hours: REAL);
(*
Difference between standard time (local time NOT adjusted
for daylight saving time or winter time) and UT or UTC.
Ex.: Switzerland (East of Greenwich): time zone = -1h.
*)
PROCEDURE SetTimeDelay( hours: REAL);
PROCEDURE GetTimeDelay(VAR hours: REAL);
(*
Time delay relative to time zone to correct for daylight
saving time. Typically summer time: +1.0; winter time:
0.0. The routines from this module know nothing about
daylight saving time rules. You need to call this
routine with hours = 0 at a switch from daylight saving
to standard or winter time or with typically hours = 1 at
a date when daylight saving time becomes active. Note,
this module knows nothing about daylight saving time
rules and will always add the specified time delay
'hours' in any computation involving local time. Ex.:
Switzerland: winter time (=standard time): hours = 0.0;
daylight saving time (from 3rd Sunday in March till last
Sunday in October): hours = 1.0. For other locations see
e.g. <http://en.wikipedia.org/wiki/Time_zone>.
*)
(*------------------------*)
(*===== Set time =====*)
(*------------------------*)
(*
Following routines are to be used to set the current day
within the year. You can freely use either SetJulianDay or
SetDayOfYear to achieve the same effect. You need to call
one of these set routines to use any solar position routine
for another day than 1st January (module initialization).
*)
PROCEDURE SetDayOfYear( dayOfYear: INTEGER);
PROCEDURE GetDayOfYear(VAR dayOfYear: INTEGER);
(* dayOfYear = 1 ~ 1st January;
dayOfYear = 365 (not leap year) or 366 (leap year) ~ 31st December
see also modules JulianDays (ScienceLib) and WriteDatTim (AuxLib)
*)
PROCEDURE SetJulianDay( julianDay: INTEGER);
PROCEDURE GetJulianDay(VAR julianDay: INTEGER);
(* same as SetDayOfYear, GetDayOfYear. Must not be confounded
with a true Julian Day (see also module JulianDays. Only kept
for upward compatibility reasons. *)
(********************************)
(*##### Sun's position #####*)
(********************************)
CONST
delta0 = 0.4092797096; (* declination of the Earth in radians *)
PROCEDURE SolarTime(localTime: REAL): REAL;
(*
Returns for the day as set by SetDayOfYear or
SetJulianDay the solar time given the localTime (=
standard time possibly adjusted for daylight saving time)
at a location with the time zone as defined by last call
to procedures and SetTimeDelay. The solar time is
defined to be 0 at noon (note, latitude as well as
longitude are ignored).
*)
PROCEDURE LocalTime(solarTime: REAL): REAL;
(*
Returns for the day as set by SetDayOfYear or
SetJulianDay the local time (= standard time possibly
adjusted for daylight saving time) given solarTime at a
location with the time zone as defined by last call to
procedures SetTimeZone and SetTimeDelay (note, latitude
as well as longitude are ignored). Standard time (=
winter time) is the local time without any adjustments
for daylight saving and defined by a multiple of half
hours offset from UT.
*)
PROCEDURE Noon(): REAL;
(*
Returns local time (= standard time) of true solar noon
for the day as set by SetDayOfYear or SetJulianDay at a
location with latitude and longitude as given in last
call to SetGeogrLocation and with the time zone as
defined by last call to procedures SetTimeZone and
SetTimeDelay.
*)
PROCEDURE SunHeightMax(): REAL; (* [1] Eq.13 *)
(* = SunHeight(Noon)); see procedure SunHeight *)
PROCEDURE SunHeight(localTime: REAL): REAL; (* [1] Eq.7 *)
(*
Returns the angle between the horizontal plane and the
sun (measured vertically from the point on the horizon
directly under the sun (cf. azimuth) up to the sun).
for the day as set by SetDayOfYear or SetJulianDay at a
location with latitude and longitude as given in last
call to SetGeogrLocation and with the time zone as
defined by last call to procedures SetTimeZone and
SetTimeDelay.
*)
PROCEDURE SunAzimuth(localTime: REAL): REAL; (* [1] Eq.8 *)
(*
Returns the south based azimuth, i.e. the angle between
the south direction and the sun for any given local time
(= standard time) for the day as set by SetDayOfYear or
SetJulianDay at a location with latitude and longitude as
given in last call to SetGeogrLocation and with the time
zone as defined by last call to procedures SetTimeZone
and SetTimeDelay. The returned angle lies between -Pi
and +Pi, solar time noon having an azimuth of 0. Often
azimuths are also given clockwise from true north to the
point on the horizon directly below the sun. To compute
the latter add Pi() (~180.0°) to the value returned by
SunAzimuth.
*)
(*
Zenith angle z
--------------
In following procedures the parameter z is the zenith angle
in radians between the sun and the zenith (The zenith being
the "top of the sky", i.e. the point on the imaginary
celestial sphere right above the head of the observer).
Depending on how exactly you define sunset and sunrise, the
following values are commonly used for z when computing day
length, sunset, sunrise:
- 90° (centre of the sun's disc is on the horizon, defining
astronomical sunset or sunrise);
- 90.8333° (upper margin of the sun's disc is on the horizon
given sun's radius or half of angular diameter corresponds
to 0.26667° and assuming a mean reffraction of 0.5667°, see
below z0);
- 96° (including civil twilight, which starts and ends when
the sun's centre is 6° below the horizon);
- 102° (nautical twilight, which starts and ends when the
sun's centre is more than 6° but less than 12° below the
horizon);
- 108° (astronomical twilight, which starts and ends when the
sun's centre is more than 12° but less than 18° below the
horizon)
*)
VAR z0: REAL; (* READ ONLY variable *)
(* = Rad(90.8333); see above "zenith angle z" *)
PROCEDURE DayLength(z: REAL): REAL; (* [2] Eq. 6.14 *)
(*
Returns the day length as a fraction of a day. Multiply
by 1440=360*4 (1 day ~ 360°, 1° ~ 4 minutes of elapsed
time) to obtain this value in minutes (of time).
<http://www.sunlit-design.com/infosearch/equationoftime.php>
*)
PROCEDURE SunRiseTime(z: REAL): REAL;
(*
Returns local time (= standard time) of sunrise for the
day as set by SetDayOfYear or SetJulianDay at a location
with latitude and longitude as given in last call to
SetGeogrLocation and with the time zone as defined by
last call to procedures SetTimeZone and SetTimeDelay.
*)
PROCEDURE SunSetTime(z: REAL): REAL;
(* As SunRiseTime but for sunset *)
(***************************)
(*##### Auxiliary #####*)
(***************************)
PROCEDURE EquationOfTime(): REAL;
(*
Difference between solar time and universal standard time
UT (or UTC) here returned as a fraction of a day
(multiply by 1440=360*4 (1 day ~ 360°, 1° ~ 4 minutes of
elapsed time) to obtain this value in minutes of time).
<http://www.sunlit-design.com/infosearch/equationoftime.php>
*)
PROCEDURE Delta(): REAL;
(*
Declination of the Sun in radian (angle between sun rays
and the plane of the earth equator). Delta is 0.0 at the
equinoxes and at a maximum at the solstices (+23°27' or
-23°27'). e.g. <http://en.wikipedia.org/wiki/Declination#Sun>
*)
END SunPath.
RAMSES
Modules:
A-Z
Function
Layer
QuickRefs:
DM
AuxLib
AuxLibE
SciLib
EasyMW
MW
ISIS
RMSLib