Here is the source code for POTS. You can read the website here

https://sites.google.com/site/positionofthesun/

Here is an example of how you use it

```
using ModPOTS
# enter the time in UTC
# enter the location in longtitude and latitude
yeardata = (2013, 6, 9)
hourdata = (23,24,0.0)
longtitudedeg = 145.15805555
latitudedeg = -37.862638888
CalcPositionOfSun(yeardata, hourdata, longtitudedeg, latitudedeg) |> println
```

Here is the source code in a module. You dont have to use the module structure, you can simply put the source code directly into your program.

```
module ModPOTS
# exporting all the functions
export displaydatetimetype, deg2hrs, rad2hrs, hrs2rad
export CalcJulianDayNumber, CalcJulianDay00, CalcJulianDayWithHMS
export CalcJulianDayWithHMS_Since_2451545, ObliquityOfTheEcliptic
export CalcGMTSiderealTime, CalcLocalSiderealTime
export CalcHourAngle, JulianDayNumberToGregorianDate
export DecimalToHMS, DecimalToDMS, HMSToDecimal, DMSToDecimal
export JulianDayNumberToGregorianDateHour
export UTC2LocalTime, LocalTime2UTC
export EclipticalCoor2EquatorialCoor
export EquatorialCoor2EclipticalCoor
export Sun_Mean_Ecliptic_Long_Since_2451545
export Sun_Perigee_Ecliptic_Long_Since_2451545
export HADec2AzAlt, CalcSunMeanAnomaly, CalcVee
export CalcPositionOfSun, MinAngleDegDiff
# Type aliases
const datearray = Tuple{Int64,Int64,Int64}
const HMStype = Tuple{Int64,Int64,Float64}
const DMStype = Tuple{Int64,Int64,Float64}
const datetimetype = Tuple{datearray,HMStype}
export datearray,HMStype,DMStype,datetimetype
# Constants
const unitdegree = Base.pi / 180.0
export unitdegree
# Function displaydatetimetype
function displaydatetimetype(x::datetimetype)
println("[",x[1][1],",",x[1][2],",",x[1][3],"] (",x[2][1],",",x[2][2],",",x[2][3],")" )
end
# Function sign
function sign(x::Float64)::Float64
if x < 0.0
return -1.0
elseif x > 0.0
return 1.0
elseif x == 0.0
return 0.0
else
return (1.0/0.0)
end
end
# Function deg2rad
@inline function deg2rad(x::Float64)::Float64
return (x * Base.pi / 180.0)
end
# Function rad2deg
@inline function rad2deg(x::Float64)::Float64
return (x * 180.0 / Base.pi)
end
# Function deg2hrs
@inline function deg2hrs(x::Float64)::Float64
return (x / 15.0)
end
# Function rad2hrs
@inline function rad2hrs(x::Float64)::Float64
return (x * 12.0 / Base.pi)
end
# Function hrs2rad
@inline function hrs2rad(x::Float64)::Float64
return (x * Base.pi / 12.0)
end
# Function CalcJulianDayNumber
# Calculate Julian Day Number when given date in UTC at midday
function CalcJulianDayNumber(date::datearray)
year = date[1]
month = date[2]
day = date[3]
a = (14 - month) Ă· 12
y = year + 4800 - a
m = month + 12 * a - 3
return ( day + ((153 * m + 2) Ă· 5) + 365 * y + (y Ă· 4) - (y Ă· 100) + (y Ă· 400) - 32045 )
end
# Function CalcJulianDay00
# Calculate Julian Day Number when given date in UTC at midnight
function CalcJulianDay00(date::datearray)
return Float64(CalcJulianDayNumber(date)) - 0.5
end
# Function CalcJulianDayWithHMS
# Calculate Julian Day Number when given date and time in UTC
function CalcJulianDayWithHMS(date::datearray,time::HMStype)
return CalcJulianDay00(date) + Float64(time[1])/24.0 + Float64(time[2])/(24.0 * 60.0) + time[3]/(24.0 * 60.0 * 60.0)
end
# Function CalcJulianDayWithHMS_Since_2451545
# Calc JuliaDay since 2451545
function CalcJulianDayWithHMS_Since_2451545(date::datearray,time::HMStype)
return CalcJulianDayWithHMS(date,time) - 2451545.0
end
# Function ObliquityOfTheEcliptic
# Where is n = JulianDay - 2451545 (1 st Jan 2000)
function ObliquityOfTheEcliptic(n::Float64)
return 23.43927944444444444444 - 3.565569769564225416381E-7 * n - 3.812460864886209963630E-17 * n^2.0 + 1.142074657903306174053E-20 * n^3.0 - 8.9899874245583241720389E-29 * n^4.0 - 1.854539576326295103555837E-34 * n^5.0
end
# Function CalcGMTSiderealTime
# Calc GMT Sidereal Time when given date and time in UTC
function CalcGMTSiderealTime(date::datearray, time::HMStype)
JD00 = CalcJulianDay00(date)
UT = Float64(time[1]) + Float64(time[2])/60.0 + time[3]/(60.0 * 60.0)
return mod((6.656306 + 0.0657098242 * (JD00 - 2445700.5 ) + 1.0027379093 * UT) , 24.0)
end
# Function CalcLocalSiderealTime
# Calc Local Sidereal Time in hours given GMTSidereal Time and longitude (in radians)
function CalcLocalSiderealTime(GMTSiderealTime::Float64, LongRadian::Float64)
return mod(( GMTSiderealTime + LongRadian/(15.0 * unitdegree) ) , 24.0)
end
# Function CalcHourAngle
# Calc Hour Angle (in hours) when given LocalSiderealTime (in hours) and RightAscension
function CalcHourAngle(LocalSiderealTime::Float64,RAinHours::Float64)
return mod((LocalSiderealTime - RAinHours) , 24.0)
end
# Function JulianDayNumberToGregorianDate
# Calc Gregorian {yyyy, mm, dd} when given JulianDay Number
function JulianDayNumberToGregorianDate(juliandaynum::Int64)
jnum = juliandaynum + 32044
gnum = jnum Ă· 146097
dgnum = mod(jnum , 146097)
cnum = (((dgnum Ă· 36524) + 1) * 3) Ă· 4
dcnum = dgnum - (cnum * 36524)
bnum = dcnum Ă· 1461
dbnum = mod(dcnum , 1461)
anum = ((dbnum Ă· 365) + 1 ) * 3
anum = anum Ă· 4
danum = dbnum - (anum * 365)
ynum = gnum * 400 + cnum * 100 + bnum * 4 + anum
mnum = ((danum*5 + 308) Ă· 153) - 2
dd = danum - ( ((mnum + 4) * 153) Ă· 5 ) + 123
yyyy = ynum - 4800 + ( (mnum + 2) Ă· 12 )
mm = ( mod((mnum + 2) , 12) ) + 1
return SVector{3,Int64}(yyyy, mm, dd)
end
# DecimalToHMS
# Convert a decimal number to a HMS tuple
function DecimalToHMS(x::Float64)
signnum = ModPOTS.sign(x)
absx = abs(x)
temp = floor(absx)
h = Int64(signnum * temp)
y = (absx - temp) * 60.0
temp = floor(y)
return (h,Int64(signnum * temp),signnum * (y - temp) * 60.0)
end
# DecimalToDMS
# Convert a decimal number to a DMS tuple
DecimalToDMS = DecimalToHMS
# Function HMSToDecimal
function HMSToDecimal(x::HMStype)
return Float64(x[1]) + (Float64(x[2])/60.0) + (x[3]/(60.0*60.0))
end
# Function DMSToDecimal
function DMSToDecimal(x::DMStype)
return Float64(x[1]) + (Float64(x[2])/60.0) + (x[3]/(60.0*60.0))
end
# JulianDayNumberToGregorianDateHour
#
function JulianDayNumberToGregorianDateHour(juliandaynum::Float64)
juliandaynumint = Int64( floor(juliandaynum) )
juliandaynumfrac = juliandaynum - Float64( juliandaynumint )
if juliandaynumfrac >= 0.5
juliandaynumint = juliandaynumint + 1
juliandaynumfrac = juliandaynumfrac - 0.5
HMStemp = DecimalToHMS(juliandaynumfrac*24.0)
else
HMStemp = DecimalToHMS((juliandaynumfrac+0.5)*24.0)
end
datetemp = JulianDayNumberToGregorianDate(juliandaynumint)
return (datetemp,HMStemp)
end
# Function UTC2LocalTime
# Convert UTC date and time into local date and time when given the TimeZone offset and DaylightSaving (1.0 is on , 0.0 is off)
function UTC2LocalTime(date::datearray, time::HMStype, TZoffset::Float64, DaylightSaving::Float64)
JD = CalcJulianDayWithHMS(date, time) + (TZoffset + DaylightSaving)/24.0
Wholeday = floor(JD)
Fractionday = JD - Wholeday
if Fractionday >= 0.5
Wholeday = Wholeday + 1
Fractionday = Fractionday - 0.5
else
Fractionday = Fractionday + 0.5
end
resultdate = JulianDayNumberToGregorianDate(Int64(Wholeday))
resulttime = DecimalToHMS(Fractionday * 24.0)
return (resultdate,resulttime)
end
# Function LocalTime2UTC
# Convert local date and time into UTC date and time when given the TimeZone offset and DaylightSaving (1 is on , 0 is off)
function LocalTime2UTC(date::datearray,time::HMStype, TZoffset::Float64, DaylightSaving::Float64)
JD = CalcJulianDayWithHMS(date,time) - (TZoffset + DaylightSaving)/24.0
Wholeday = floor(JD)
Fractionday = JD - Wholeday
if Fractionday >= 0.5
Wholeday = Wholeday + 1
Fractionday = Fractionday - 0.5
else
Fractionday = Fractionday + 0.5
end
resultdate = JulianDayNumberToGregorianDate(Int64(Wholeday))
resulttime = DecimalToHMS(Fractionday * 24.0)
return (resultdate,resulttime)
end
# Function EclipticalCoor2EquatorialCoor
#
function EclipticalCoor2EquatorialCoor(lambda::Float64, beta::Float64, epsilon::Float64)
alpha = atan(sin(lambda) * cos(epsilon) - tan(beta) * sin(epsilon),cos(lambda))
delta = asin(sin(beta) * cos(epsilon) + cos(beta) * sin(epsilon) * sin(lambda));
return [alpha,delta]
end
function EquatorialCoor2EclipticalCoor(alpha::Float64, delta::Float64, epsilon::Float64)
lambda = atan(sin(alpha) * cos(epsilon) + tan(delta) * sin(epsilon), cos(alpha) )
beta = asin( sin(delta) * cos(epsilon) - cos(delta) * sin(epsilon) * sin(alpha) )
return SVector{2,Float64}(lambda,beta)
end
# Function Sun_Mean_Ecliptic_Long_Since_2451545
# Calculate The Mean Longitude of the Sun (in degrees) given number of days since midday 1st Jan 2000
function Sun_Mean_Ecliptic_Long_Since_2451545(n::Float64)
return mod(280.465900300 + 0.985647351813826146 * n + 2.26748764711145967E-13 * n^2.0 , 360.0)
end
# Function Sun_Perigee_Ecliptic_Long_Since_2451545
# Calculate The Perigee Longitude of the Sun (in degrees) given number of days since midday 1 st Jan 2000
function Sun_Perigee_Ecliptic_Long_Since_2451545(n::Float64)
return mod((282.940472178 + 0.0000470932390417522245 * n + 3.39394552688870243E-13 * n^2.0) , 360.0)
end
# Function HADec2AzAlt
# Input : Hour Angle, Declination and Latitude. All three in radians
# Output: Azimuth and Altitude in radians
function HADec2AzAlt(HArad::Float64, delta::Float64, lat::Float64)
Azimuth = mod( 2*Base.pi + Base.pi + atan(sin(HArad),cos(HArad)*sin(lat) - tan(delta)*cos(lat) ) , 2 * Base.pi)
Altitude = asin( sin(lat)*sin(delta) + cos(lat)*cos(delta)*cos(HArad) )
return (Azimuth,Altitude)
end
# Convert Equatorial Coordinates to Horizontal Coordinates
#(* Input: alpha RA, delta DEC, lat, long, {utcdate}, {utctime}
# In radians are alpha RA, delta DEC, lat and long.
# {utcdate} and {utctime} in UTC
# Output Azimuth and Altitude in radians.
function EquatorialCoor2HorizontalCoor(alpha::Float64, delta::Float64,
lat::Float64, long::Float64, utcdate::datearray, utctime::HMStype)
alphahours = alpha / (15.0 * unitdegree)
GMTSideRealTime = CalcGMTSiderealTime(utcdate, utctime)
LocalSideRealTime = CalcLocalSiderealTime(GMTSideRealTime, long)
hourangle = CalcHourAngle(LocalSideRealTime, alphahours)
HArad = 2.0 * Base.pi * hourangle/24.0
Azimuth_Altitude = HADec2AzAlt(HArad, delta, lat)
return Azimuth_Altitude
end
# Convert Equatorial Coordinates to Horizontal Coordinates
# Input: alpha RA, delta DEC, lat, long, {date}, {time}
# In radians are alpha RA, delta DEC, lat and long.
# {date} and {time} in Local time
# Output Azimuth and Altitude in radians.
function EquatorialCoor2HorizontalCoorĚ†localtime(alpha, delta, lat,
long, Localdate::datearray, Localtime::HMStype, TZoffset, DaylightSaving)
(utcdate, utctime) = LocalTime2UTC(Localdate, Localtime, TZoffset, DaylightSaving)
return EquatorialCoor2HorizontalCoor(alpha, delta, lat, long,
utcdate, utctime)
end
function CalcSunMeanAnomaly(yeardata::datearray, hourdata::HMStype)
JD = CalcJulianDayWithHMS(yeardata, hourdata)
n = CalcJulianDayWithHMS_Since_2451545(yeardata, hourdata)
SunMeanLong = deg2rad(Sun_Mean_Ecliptic_Long_Since_2451545(n))
perigeelongitude = deg2rad(Sun_Perigee_Ecliptic_Long_Since_2451545(n))
SunMeanAnomaly = mod( SunMeanLong - perigeelongitude , 2 * Base.pi )
return (SunMeanAnomaly, perigeelongitude, JD, n)
end
# ecc = 0.01671123
# SunMeanAnomalyDeg is 155.31965461588
# SunMeanAnomaly = 2.710839366107514
# solve x - ecc * sin(x) == SunMeanAnomaly
# Solution x is 2.71771269701826570
# Function CalcVee
# Input: SunMeanAnomaly in unit radians
# ecc (eccentricity of the orbit)
# Output: v in unit radians
function CalcVee(SunMeanAnomaly, ecc)
e = 0.0
loop = true
while loop
old_e = e
e = SunMeanAnomaly + ecc * sin(e)
if abs(e - old_e) < 1.0e-14
loop = false
end
end
v = 2 * atan( sqrt((1 + ecc)/(1 - ecc)) * tan(e/2) ) # in radians
v = mod(v, 2 * Base.pi)
return v
end
function CalcPositionOfSun(yeardata, hourdata, longtitudedeg, latitudedeg)
longtitude = deg2rad(longtitudedeg)
latitude = deg2rad(latitudedeg)
(SunMeanAnomaly, perigeelongitude, JD, n) = CalcSunMeanAnomaly(yeardata, hourdata)
ecc = 0.01671123
v = CalcVee(SunMeanAnomaly, ecc)
lambda = mod(v + perigeelongitude, 2.0 * Base.pi)
epsilon = deg2rad(ObliquityOfTheEcliptic(n))
beta = 0.0
(alpha, delta) = EclipticalCoor2EquatorialCoor(lambda, beta, epsilon)
alphahours = rad2hrs(alpha)
theta_gmt = CalcGMTSiderealTime(yeardata, hourdata)
theta_local = CalcLocalSiderealTime(theta_gmt, longtitude)
hourangle = CalcHourAngle(theta_local, alphahours)
HArad = hrs2rad(hourangle)
(Azimuth, Altitude) = HADec2AzAlt(HArad, delta, latitude)
Azimuthdeg = rad2deg(Azimuth)
Altitudedeg = rad2deg(Altitude)
return (Azimuthdeg,Altitudedeg)
end
function MinAngleDegDiff(x,y)
(xval,yval) = (x,y)
while yval<xval
yval += 360.0
end
diff = mod(yval - xval,360)
if diff > 180.0
diff = 360.0 - diff
end
return diff
end
end
```