Sun's trajectory around earth

General Usage
1

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.

Thanks in advance,
Tomas

2

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
3

Please read: https://www.researchgate.net/profile/Chee-Woon-Wong/publication/221905606_General_Formula_for_On-Axis_Sun-Tracking_System/links/09e41510713662528d000000/General-Formula-for-On-Axis-Sun-Tracking-System.pdf

4

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

5

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

6

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.

7

What about the fun of it?

7 Likes
8

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
9

it’s all good fun until your panels get blown down :slight_smile:

1 Like
10

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
11

Not all locations are stormy…

12

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

13

I’ve always been told that earth is on a trajectory around the sun…!?

4 Likes
14

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
15

That’s what I said GMT’s solar see refs at the end of solar — GMT 6.4.0 documentation

16

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…

17

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
18

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.

19

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