Sun's trajectory around earth

Hi All,

My father (82) has decided to put a photovoltaic panels to his flat roof. And since he is a die hard engineer and wants to optimizing everything, he decided that he will allow to change angle of panels, but only along one axis, since two axis system is too complicated. The axis of tilting the panels would east-west and according to his computations, he has figured out that at morning, the panels should be horizontal, then they should turn 45 to face south and then in the evening, they should be horizontal.

I would like to verify his model, hence I would like to ask, if someone has a model of sunâ€™s trajectory around earth, which would allow me to compute angle from north-south axis and also the elevation angle from the horizon.

Tomas

I have Stephan Buchert / Sun position Â· GitLab which my conversion of

Roberto Grena (2012), Five new algorithms for the computation of sun position
from 2010 to 2110, Solar Energy, 86(5):1323â€“1337, doi:10.1016/j.solener.2012.01.024.

from C++ to Julia. It is not a package, download `SunPosition.jl` and make it available in your installation as a module.

3 Likes

Thanks of you. This is what I have needed. Just add and run.

it really isnâ€™t worth the optimization

The cost of the components to move the panels outweighs any gains, especially as they have to be strong enough to resist wind

A fixed position installation can have far more robust fittings at a fraction of the cost

Iâ€™m pretty sure you can do it with routines in `AstroLib.jl`. So an example of use at Drawing the analemma with Julia to draw the analemma in a given place.

What about the fun of it?

7 Likes

I have a hard time believing that an east-west axis is better than a north-south axis, and that 45 degrees is optimal for any place besides 45 degrees latitude. For example, at the equator, youâ€™d want the axis to be aligned north-south and parallel to the ground, so that the panels faced due east at dawn and due west at sunset. And at the north pole, youâ€™d want it orthogonal to the ground with the panels rotating once per day keeping them facing the sun as it circles the horizon. (This would be exactly optimal only at the equinoxes, but with a fixed axis itâ€™ll give you the best mean flux over the year).

In general, for a fixed-axis panel, youâ€™ll want the axis aligned north-south at an angle to the ground equal to the latitude, so that the normal of the panel traces out the equator on the sphere of fixed stars. Thatâ€™s using the equator as an mean path of the sun in the sky over the course of a year.

1 Like

itâ€™s all good fun until your panels get blown down

1 Like

GMTâ€™s module `solar` gives info like this

``````using GMT

julia> solar(sun=(pos=(-8,37), date=("2021-5-23")))
9-element Vector{String}:
"Sun current position:"
"\tLongitude = -172.814944"
"\tLatitude  = 20.579669"
"\tAzimuth   = 352.0617"
"\tElevation = -32.0227"
"\n\tSunrise   = 05:18"
"\tSunset    = 19:39"
"\tNoon      = 12:29"
"\tDuration  = 14:21"
``````
2 Likes

Not all locations are stormyâ€¦

A 1m^2 panel experiencing a 10mph headwind is the same as a 6kg / 13lb weight being attached

Iâ€™ve always been told that earth is on a trajectory around the sunâ€¦!?

4 Likes

Look at the sky and see yourselfâ€¦

Seriously, the JPL ephemeris which the Astrolib package apparently is wrapping is based on numerical integration of the â€śtrajectoriesâ€ť of the solar-system n-bodies. This is heavy computations and the highest available precision, allowing to navigate the different spacecraft to sling-shot around planets in order to reach other planets, comets, astroids, etc. Written in C++?

For solar energy the JPL ephemeris is of course an overkill, so they use SPA (Sun Position Algorithm) or with further simplifications by Grena. This just models the Sun trajectory over the 2-D sky of the Earth using suitable functions. Much less precise, but sufficient and good especially if it is needed for billions of points in time and geographic locations, or in embedded controllers of solar panelsâ€¦

1 Like

Thatâ€™s what I said GMTâ€™s `solar` see refs at the end of solar â€” GMT 6.4.0 documentation

Can Makie do shadows? It would be really cool if someone with mad Makie skills could create a â€śsundialâ€ť example, where the sun position functions control the light position in Makie, and show the shadowâ€™s cast by a cylinder (or better yet, model Stonehenge :-)). (I could be the first person to optimize my Burning Man Geodesic shade structure with Juliaâ€¦ ).

I also remember reading about a pure Julia ray-tracer, a sundial could be a nice example for that as wellâ€¦

Thanks a lot for wonderful answers, I new the community will let me down.

My intention is to revisit calculations my father did to confirm it. Frankly, i wish he succeeds, but what is more important to me is that he has the project. He is 82, 11 years on dialysis. This project brings fresh blood to his veins.

9 Likes

I have wrote a julia program called POTS.jl which stands for
â€śPosition of the sunâ€ť. I even have a website to explain the calculations. Let me find the details.

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

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 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 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

return (x * Base.pi / 180.0)
end

return (x * 180.0 / Base.pi)
end

# Function deg2hrs
@inline function deg2hrs(x::Float64)::Float64
return (x / 15.0)
end

return (x * 12.0 / Base.pi)
end

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
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)
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

# Input :  Hour Angle, Declination and Latitude. All three in radians
# Output:  Azimuth and Altitude in radians
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
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)
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)
(SunMeanAnomaly, perigeelongitude, JD, n) = CalcSunMeanAnomaly(yeardata, hourdata)
ecc = 0.01671123
v = CalcVee(SunMeanAnomaly, ecc)
lambda = mod(v + perigeelongitude, 2.0 * Base.pi)
beta = 0.0
(alpha, delta) = EclipticalCoor2EquatorialCoor(lambda, beta, epsilon)
theta_gmt = CalcGMTSiderealTime(yeardata, hourdata)
theta_local = CalcLocalSiderealTime(theta_gmt, longtitude)
hourangle = CalcHourAngle(theta_local, alphahours)
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
``````