Tronald / CoordinateSharp

A library designed to ease geographic coordinate format conversions, and determine sun/moon information in C#
Other
361 stars 59 forks source link

Calculate time change based on azimuth change from solar noon #222

Closed houstonhaynes closed 1 year ago

houstonhaynes commented 1 year ago

Is your feature request related to a problem? Please describe. I'm looking for a calculation that's similar to Get_Time_At_Solar_Altitude but instead using azimuth.

Describe the solution you'd like The idea is to abstract away the compensation necessary to adjust relative to the general "15 degrees per hour" case, per the user's date and position.

Describe alternatives you've considered I have considered looking at Get_Time_At_Solar_Altitude and creating a "Get_Time_At_Solar_Azimuth" myself - or perhaps building a functional equivalent in F#.

Additional context I thought it might be more useful for there to be a C#/CoordinateSharp "sibling function" to Get_Time_At_Solar_Altitude that everyone could use from the library.

Tronald commented 1 year ago

Hi @houstonhaynes

Thank you for the suggestion. We will look into this. I am not sure off the top of my head how accurate such a simple 15 degree per hour calculation will be (even though it sounds logical), but it may be within a acceptable tolerances if we can account for equation of time (EOT).

time = [ΞΈsunβˆ’ΞΈlong]Γ—.25hrs (15 minutes) +EoT+12hrs

Reference Meeus Astronomical Algorithms Ch 28 for EOT.

houstonhaynes commented 1 year ago

Thanks for your notes - I tried a few things around creating my own EoT function and it's just too off. I got some help on my helpers and noodled around with GitHub CoPilot and ChatGPT a bit -- but just couldn't "close the gap". For tomorrow at my location I should be seeing a pretty tight window (40 minutes) for +/-22 degrees around 'solarNoon'. A three-hour span that's off by a day is double trouble. 😞

#r "nuget: CoordinateSharp"

open System
open CoordinateSharp

let calculateSolarNoon latitude longitude (date: DateTime) =
    let celestialTimes = Celestial.CalculateCelestialTimes(latitude, longitude, date, 0.0)
    celestialTimes.SolarNoon.Value.ToLocalTime()

// Helper function to convert degrees to radians
let toRadians (degrees: float) = degrees * Math.PI / 180.0

// Helper for time zone offset
let localTimeZoneOffset = DateTimeOffset.Now.Offset
let timeZoneOffset = TimeSpan.FromTicks(localTimeZoneOffset.Ticks)

// Helper function to calculate the Equation of Time
let calculateEoT longitude (date: DateTime) =
    let julianDay = float date.DayOfYear
    let b = (julianDay - 81.0) * 360.0 / 364.0 |> toRadians
    let eot = 9.87 * sin(2.0 * b) - 7.53 * cos(b) - 1.5 * sin(b)
    let eotMinutes = eot * 4.0
    let eotTimeSpan = TimeSpan.FromMinutes(eotMinutes)
    let eotTimeSpanWithTimeZone = eotTimeSpan + timeZoneOffset

    eotTimeSpanWithTimeZone

// Helper function to calculate the time of a given azimuth
let calculateTimeFromAzimuth latitude longitude (date: DateTime) azimuthCompensation =
    let observer = Coordinate(latitude, longitude, date)
    let celestialInfo = observer.CelestialInfo
    let eot = calculateEoT longitude date
    let solarNoonTime = celestialInfo.SolarNoon.Value.ToLocalTime().TimeOfDay
    let timeDifferenceMinutes = (celestialInfo.SunAzimuth + azimuthCompensation) * 4.0
    let observerDateTime = date.Date + solarNoonTime + eot + TimeSpan.FromMinutes(timeDifferenceMinutes)

    observerDateTime

let date = DateTime.Today.AddDays(1.0).ToLocalTime() // Tomorrow's date
let latitude = 35.62740870 // Latitude of Asheville, NC
let longitude = -82.53818729 // Longitude of Asheville, NC
let solarNoon = calculateSolarNoon latitude longitude date

let minTime = calculateTimeFromAzimuth latitude longitude date (-22.0)
let maxTime = calculateTimeFromAzimuth latitude longitude date (22.0)
val solarNoon: DateTime = 7/19/2023 1:37:31 PM
val minTime: DateTime = 7/20/2023 11:19:01 AM
val maxTime: DateTime = 7/20/2023 2:15:01 PM
Tronald commented 1 year ago

I appreciate you attempting to go down this rabbit hole. I have noticed that LLMs still struggle with many astronomical algorithms for some reason. I was surprised by this, but I guess it make sense considering much of this knowledge remains in books/academic papers that are not published fully on the internet. Combine that with AI "hallucination" and a lack of internal testing of the algorithms it creates and you get terrible outputs more often than not. I am sure LLMs will get there eventually, but sadly we are still forced to work the trig.

Have you tried removing the EoT from your logic? If you do and the values are within -14 to +41 minutes of what they should be, you may be close and just need to adjust the EoT. The EoT in Meeus's book appears much more involved (just popped it open).

I plan to dive into this in my spare time, but please continue to provide insight you find it.

houstonhaynes commented 1 year ago

Great points. There were QUITE a few options when I handed the reins over to GitHub CoPilot. I looked at them and thought "well there's been some folks that have given this a stab and this thing just hasn't even bothered trying to consolidate them into one that works" 😸

I'll spend some time with the reference sometime - as this is something that I want to flex on being right "enough" to make an accurate solar tracker without direct sensing of the Sun's position in the sky. The design is such that the panel can only swing as far as +/- 22 degrees from "flat" (or as we refer to it "parking position"). Many studies of commercial solar tracking arrays is showed that eccentric repositioning of the panels to position perpendicular to the Sun at all hours of daylight doesn't make as much of a different in the early and late hours of the day. (as you might expect) The mechanical design I have prototyped "tips a hat" to that by looking to "spread the sweet spot" in the daylight hours that are "central" to the panel/array's position - to maximize the increase in collection for those hours where it more significantly "moves the needle" as it were.

So the basic concept for the mechanism would be for the panel to position at its east-facing height sometime around sunrise and hold that position until the Sun matched the azimuth that meets the -22 degree-from-solar-noon angle and track with it until it's to the +22 degree-from-solar-noon position. From there it would hold until again sometime around sunset and then it would "park" overnight.

The point is to maximize whatever advantage is there for this design. From that perspective I want to make this as accurate as I am able in order to illustrate that the methodology is "worth the trouble" 😁 I can use a site like sunfollower.net and "hand jam" a table of data for a week or so -- to set up the panel to test and check against a static/flat positioned solar cell. The important point at this early stage is to measure the actual difference in power collected by the "tracking" versus flat panel. But eventually this needs to be free-running, with the WildernessLabs Meadow "Feather" microcontroller doing all of the math. πŸ€– πŸ€“

That's a very long way of saying "I've got time - but I'd like to get it as right as I can" 🧠

houstonhaynes commented 1 year ago

Let me pivot for a moment and ask something that might be more orthogonal... would it be possible to involve solar altitude in the calculation for a given day? I'm not sure if that would avoid EoT completely, but it just occurred that if I can get the solar altitude for a given azimuth displacement, then that could give me the time. Now that I think of it, that would still involve EoT so never mind. 😺

houstonhaynes commented 1 year ago

This borders on truly awful 😈 but I decided to just "adjust" EoT on the fly and create a linear compensation for the next 90 days to give me a bit of "runway" until I can figure out a more elegant solution. And of course if I get to mid-October and just need to do it again. I'll - well - do it all over again. πŸ˜΅β€πŸ’«

#r "nuget: CoordinateSharp"

open System
open CoordinateSharp

let calculateSolarNoon latitude longitude (date: DateTime) =
    let celestialTimes = Celestial.CalculateCelestialTimes(latitude, longitude, date, 0.0)
    celestialTimes.SolarNoon.Value.ToLocalTime()

// Helper function to convert degrees to radians
let toRadians (degrees: float) = degrees * Math.PI / 180.0

// Helper for time zone offset
let localTimeZoneOffset = DateTimeOffset.Now.Offset

// helper to square a number
let sqr x = x ** 2.0

// Calculate the Equation of Time (EOT) for a given date
let calculateEoT (date: DateTime) =
    let n = float date.DayOfYear
    let b = (n - 1.0) * 2.0 * Math.PI / 365.0
    let g = 357.0 + b
    let e = 0.0167
    let sinG = Math.Sin g
    let cosG = Math.Cos g
    let c = sinG * sinG + cosG * cosG * e * e
    let d = Math.Sqrt c
    let eot = 720.0 * (d / Math.PI)

    // Helper to hack the difference in days from July 19, 2023
    let dayDifference (date: DateTime) : int =
        let staticDate = DateTime(2023, 7, 19)
        let daysPassed = (date - staticDate).Days
        daysPassed

    // Helper to hack EoT to be more accurate for 90 days after July 19, 2023
    let interpolateEmpiricalCorrection (dayOfYear: int) =
       let days = [| 0, -10.0; 15, 40.0; 30, 100.0; 45, 130.0; 60, 145.0; 75, 100.0; 90, 57.0; |]

        let rec findDataPoints (days: (int * float)[]) prevDay currentDay =
            match days with
            | [||] | [| (_, _)|] -> (prevDay, currentDay)
            | _ ->
                match days.[1] with
                | (day, _) when day < dayOfYear ->
                    findDataPoints (Array.tail days) days.[0] days.[1]
                | _ -> (prevDay, currentDay)

        let (prevDay, currentDay) = findDataPoints days days.[0] days.[1]

        let prevCorrection = snd prevDay
        let currentCorrection = snd currentDay

        let prevDayOfYear = fst prevDay
        let currentDayOfYear = fst currentDay

        let weight = float (dayOfYear - prevDayOfYear) / float (currentDayOfYear - prevDayOfYear)

        prevCorrection + weight * (currentCorrection - prevCorrection)

    // Empirical correction factor from July 2023 to October 2023
    let empiricalCorrection = interpolateEmpiricalCorrection (dayDifference date) 

    TimeSpan.FromMinutes (eot + empiricalCorrection)

// Calculate the time from azimuth with Equation of Time (EOT) compensation
let calculateTimeFromAzimuth latitude longitude (date: DateTime) azimuthCompensation =
    let timeZoneOffset = localTimeZoneOffset
    let observer = Coordinate(latitude, longitude, date)
    let celestialInfo = observer.CelestialInfo
    let eot = calculateEoT date
    let solarNoonTime = celestialInfo.SolarNoon.Value.ToLocalTime().TimeOfDay
    let azimuthOffset = azimuthCompensation / 15.0 // Convert azimuth compensation to hours

    // Adjust azimuth offset based on latitude to refine the result
    let azimuthMultiplier =
        if latitude > 0.0 then 1.15
        else if latitude < 0.0 then 0.85
        else 1.0

    let adjustedAzimuthOffset = azimuthOffset * azimuthMultiplier

    let observerDateTime = date.Date + solarNoonTime + eot + TimeSpan.FromHours(adjustedAzimuthOffset) + timeZoneOffset

    // Adjust the observerDateTime if it falls on a different day
    if observerDateTime.Date <> date.Date then
        observerDateTime - TimeSpan.FromDays(1.0)
    else
        observerDateTime

// Example usage
let latitude = 35.5951 // Latitude of Asheville, NC
let longitude = -82.5515 // Longitude of Asheville, NC
let date = DateTime.Today.AddDays(0.0)

let maxTime = calculateTimeFromAzimuth latitude longitude date (22.0)
let solarNoon = calculateSolarNoon latitude longitude date

let minTime =
    let positiveTimeDifference = maxTime - solarNoon
    let minDateTime = solarNoon - positiveTimeDifference
    minDateTime
Tronald commented 1 year ago

Yikes, just dove in and this is a bit more complex than I had imaged. The below EOT compares within one minute (at first glance) to https://planetcalc.com/9198/

You could translate to to F# as needed if you wish to continue investigating. It returns EOT in total minutes.

using CoordinateSharp;

 public static double Equation_of_Time(DateTime d)
 {
         var jde = JulianConversions.GetJulian(d);          
         var r = (jde - 2451545.0) / 365250;

         var l0 = 280.4664567 + 360007.6982779 * r + .03032028 * Math.Pow(r, 2) + Math.Pow(r, 3) / 49931 - Math.Pow(r, 4) / 15300 - Math.Pow(r, 5) / 2000000;
         l0 = l0.NormalizeDegrees360();

         var sc = Celestial.Get_Solar_Coordinate(d);
         var rightAscension = sc.RightAscension;

         var nutationInLong = -17.20 * Math.Sin(sc.Longitude) - 1.32 * Math.Sin(2 * r) - .23 * Math.Sin(2 * r) + .21 * Math.Sin(2 * sc.Longitude); //22.1
         var nutDeg = nutationInLong / 3600;// convert arcseconds to degrees.

         var ob = sc.ObliquityOfEcliptic;

         var E = l0 - .0057183 - rightAscension + nutDeg * Math.Cos(ob);

         return E * (Math.PI / 180) / (Math.PI / 720);
 }

Unfortunately the formula I suggested for acquiring the Time of Day based on the EoT doesn't work. This should be doable considering this is how sundials work and are adjusted for more accuracy. Will look into the second portion later.

houstonhaynes commented 1 year ago

This is great - thanks! I was just sitting with some of the explanations in "Astronomical Algorithms" and looking for a visual hook to see where the represented formulas might map out to some of the samples that GitHub Copilot showed me. There's a funny little quip in the text that reads

The equation of time is always smaller than 20 minutes in absolute value. If |E| appears to be too large, add 24 hours to or subtract it from your result. (p. 184, right before the example)

That's not confusing AT ALL. 😏

I'll be glad to take a look at your code first - it already looks like it has the right stuff. πŸ‘

houstonhaynes commented 1 year ago

I couldn't access NormalizeDegrees360 or ObliquityOfEcliptic. Are those private helpers in your library that aren't exposed to the API?

I wrote a normalizer and let CoPilot help me with ObliquityOfEcliptic. Other than returning calculateEoT as a TimeSpan it's pretty much one-for-one. That said - the results I'm getting for today don't line up (and I'm calculating my -22 degree-from-solar-noon value directly not simply reflecting it from the maxTime value)

[fixed below]

houstonhaynes commented 1 year ago

Output of F# Interactive

val latitude: float = 35.5951
val longitude: float = -82.5515
val date: DateTime = 7/19/2023 12:00:00 AM
val maxTime: DateTime = 7/19/2023 11:12:29 AM
val solarNoon: DateTime = 7/19/2023 1:37:34 PM
val minTime: DateTime = 7/19/2023 7:50:05 AM

I suspect the OoE helper might have led me astray but wanted to at least throw this over the fence to you before I started another dive.

Tronald commented 1 year ago

Oh shoot. Sorry ObliquityOfEcliptic isn't accessible yet my apologies (forgot I just added the property to expose it in my working update). .NormalizeDegrees360 requires reference to CoordinateSharp.Formatters.

I think your values are off however so you may want to check them. CoordinateSharp should be producing within a minute of NOAA (calculation methods are different so there is slight variation). https://gml.noaa.gov/grad/solcalc/

image

Astronomical Algorithms actually has a Sundial section (chapter 58). I am going to read it over later and see if there is a simplistic way to calculate this (as you are essentially creating a sundial), but it looks involved and complex and different areas of the world appear alter the calculation. It may be easier to create a simple generic iteration function to essentially iterate the azimuth/time associations for the given date/time between rise and set. This wouldn't be as efficient, but I don't suspect it needs to be as a function like this wouldn't be called in bulk so the loss in efficiency would be negligable.

I will look into this a bit later.

houstonhaynes commented 1 year ago

That's great. And for sure - given the relatively short span of time between -22 and +22 degrees around solar noon I'm not planning on getting to picky about plotting the change in angle of the solar panel with too fine of a grain. This is both to not "over optimize" but also simply not activating the motor that controls the position so often that the drain ends up soaking up an excess of 'value' gained through the additional power collected.

In effect - my thought was to simply take the two min/max values and plot three intermediate points on each side of solar noon. So an "expensive" call before sunrise to determine those times (plus the time of sunrise and sunset, naturally) falls well within the "no big deal" category.

If anything - it's those periods - from sunrise to solar-noon-minus-22-degrees and its counterpart later in the day which really "move the needle" in terms of the additional power that can be collected for a given day. All of the fancy pivoting in the "middle" of each day is designed to not squander the opportunity, as it were. It's a mild irony that most inverters actually "clip" at those highest levels of collection so there's a figurative "lid" on those middle hours in either case. Being slightly off in that window likely won't make that much difference in actual power collected overall. 🌞

houstonhaynes commented 1 year ago

Is there something off about returning the final calculation with E as a TimeSpan (in minutes?) Should I be using some other factor? Do you take that double and convert the result when it's implemented?

In your example you return a double - but I'm using it here [eot] where the date is adding eot and the adjustedAzimuthOffset plus timeZoneOffset. Any thoughts on whether that might be where I'm going wrong?

let observerDateTime = date + solarNoonTime + eot + TimeSpan.FromHours(adjustedAzimuthOffset) + timeZoneOffset
houstonhaynes commented 1 year ago

Here's the updated script with use of CoordinateSharp.Formatters and a few other tweaks.

#r "nuget: CoordinateSharp"

open System
open CoordinateSharp
open CoordinateSharp.Formatters

// Helper function to convert degrees to radians
let toRadians (degrees: float) = degrees * Math.PI / 180.0
// Helper for time zone offset
let localTimeZoneOffset = DateTimeOffset.Now.Offset
// helper to square a number
let sqr x = x ** 2.0

// Calculate the obliquity of the ecliptic for a given date
let calculateObliquityOfEcliptic (d: DateTime) =
    let jde = JulianConversions.GetJulian d
    let t = (jde - 2451545.0) / 365250.0
    let epsilonDeg = 23.43929111 - 0.013004167 * t - 1.63889E-7 * (t ** 2.0) + 5.036111E-7 * (t ** 3.0)
    epsilonDeg

// Calculate the Equation of Time (EOT) for a given date
let calculateEoT (d: DateTime) =
    let jde = JulianConversions.GetJulian d
    let r = (jde - 2451545.0) / 365250.0
    let l0 = 280.4664567 + 360007.6982779 * r + 0.03032028 * (r ** 2.0) + (r ** 3.0) / 49931.0 - (r ** 4.0) / 15300.0 - (r ** 5.0) / 2000000.0
    let l0Normalized = l0.NormalizeDegrees360()
    let sc = Celestial.Get_Solar_Coordinate d
    let rightAscension = sc.RightAscension
    let nutationInLong = -17.20 * sin sc.Longitude - 1.32 * sin (2.0 * r) - 0.23 * sin (2.0 * r) + 0.21 * sin (2.0 * sc.Longitude) // 22.1
    let nutDeg = nutationInLong / 3600.0 // convert arcseconds to degrees.
    let ob = calculateObliquityOfEcliptic d // Calculate the obliquity of the ecliptic
    let E = l0Normalized - 0.0057183 - rightAscension + nutDeg * cos (ob * Math.PI / 180.0)
    let minutes = E * (Math.PI / 180.0) / (Math.PI / 720.0)
    TimeSpan.FromMinutes(minutes)

// Calculate the time from azimuth with Equation of Time (EOT) compensation
let calculateTimeFromAzimuth latitude longitude (date: DateTime) azimuthCompensation =
    let timeZoneOffset = localTimeZoneOffset
    let observer = Coordinate(latitude, longitude, date)
    let celestialInfo = observer.CelestialInfo
    let eot = calculateEoT date
    let solarNoonTime = celestialInfo.SolarNoon.Value.ToLocalTime().TimeOfDay
    let azimuthOffset = azimuthCompensation / 15.0 // Convert azimuth compensation to hours

    // Adjust azimuth offset based on latitude to refine the result
    let azimuthMultiplier =
        if latitude > 0.0 then 1.15
        else if latitude < 0.0 then 0.85
        else 1.0

    let adjustedAzimuthOffset = azimuthOffset * azimuthMultiplier

    let observerDateTime = date + solarNoonTime + eot + TimeSpan.FromHours(adjustedAzimuthOffset) + timeZoneOffset

    observerDateTime

let calculateSolarNoon latitude longitude (date: DateTime) =
    let celestialTimes = Celestial.CalculateCelestialTimes(latitude, longitude, date, 0.0)
    celestialTimes.SolarNoon.Value.ToLocalTime()

// Example usage
let latitude = 35.5951 // Latitude of Asheville, NC
let longitude = -82.5515 // Longitude of Asheville, NC
let date = DateTime.Today.AddDays(0.0)
let maxTime = calculateTimeFromAzimuth latitude longitude date (22.0) // last value is degrees from solar noon
let solarNoon = calculateSolarNoon latitude longitude date
let minTime = calculateTimeFromAzimuth latitude longitude date (-22.0) // last value is degrees from solar noon
Tronald commented 1 year ago

Is there something off about returning the final calculation with E as a TimeSpan (in minutes?) Should I be using some other factor? Do you take that double and convert the result when it's implemented?

In your example you return a double - but I'm using it here [eot] where the date is adding eot and the adjustedAzimuthOffset plus timeZoneOffset. Any thoughts on whether that might be where I'm going wrong?

let observerDateTime = date + solarNoonTime + eot + TimeSpan.FromHours(adjustedAzimuthOffset) + timeZoneOffset

Hmmm, I am not really sure of the top of my head with this one sorry. I plan to patch up a bug reported by another user, and with that update you'll have access to the ObliquityOfEcliptic within the solar coordinates (hopefully in a couple days) btw.

I am hoping I can dive into this one with a bit more focus soon. It's surprisingly a much larger puzzle than I expected.

houstonhaynes commented 1 year ago

Thanks so much for batting this around with me. It's pretty wonky, a much deeper dive than I was originally prepared for... but one that I might end up blogging about or making a YouTube video for my "Glad Scientist" channel. πŸ˜„

Tronald commented 1 year ago

@houstonhaynes Are you referring to the last number in the 28.2 formula for L0 (2,000,000)? If so it is six zeros. The "space" in the text sometimes creates the illusion of seven if you don't look carefully and your eyes are shot from coding all day. It has got me before too lol.

houstonhaynes commented 1 year ago

Yup that's exactly what happened - after I took a screen shot I looked at it again and realized my misread.

I've "switched tacks" to focus on the accelerometer on the Meadow - to give X/Y/Z position of the board. But I keep going back to the book and reading bits until my brain goes on "tilt" 🚨 🀭

Tronald commented 1 year ago

2.19.1.1 was released and exposes the Obliquity of Ecliptic value within the SolarCoordinate class. I will keep this issue open until we can find a solution to the original problem. Unfortunately, we may not be able to get to it for a few weeks.

houstonhaynes commented 1 year ago

Thanks! My current hack will give me a few months of runway so long as I say at the same coordinates 😁 🌐 I'm looking forward to taking your update for a spin! I sorted out the accelerometer's Radian output conversion to degrees and I feel almost guilty for how reasonably easy it was - comparatively speaking. 😊

Tronald commented 1 year ago

Here is a "brute" program that can calculate Time at Azimuth using CoordinateSharp. This method allows you to operate with a UTC offset as well.

This method doesn't have terrible overhead, but it could present latency if doing bulk calculations (which I can't imagine why anyone would for this problem). This algorithm hasn't been tested heavy, nor do I know how it will behave at poles, but it should work for your case. There should also be validations built in to ensure azimuth degree limits aren't surpassed.

This should get us going in the right direction though.

Time at Azimuth CoordinateSharp

houstonhaynes commented 1 year ago

Great! I'll give it a look tonight, and probably a bit more early next week. I'm about to go do a few days of (relatively) basic camping in a national forest in prep for a full week during the height of the Perseid meteor event. πŸ˜† My partner and I can go to a nearby town with 5G for a few hours each weekday. That is what makes the camping "relatively" basic. πŸ˜…

houstonhaynes commented 1 year ago

DUDE that is SPOT ON! πŸ”­ πŸ‘ :dart: I just took it for a test drive in the fiddle you linked and it's really tight! Without an IDE to point "back" into your library to see how you pulled it off I'm a bit mystified, but I'm looking forward to giving it a deeper look this evening. Thank you!!

Tronald commented 1 year ago

Enjoy your camping trip! Glad the logic seems to work.

houstonhaynes commented 1 year ago

Thanks! I plan to F#-ify this code in order to put it into a "module" (or even a purpose-build library) that I can plug into the WildernessLabs Meadow board and/or the "backplane" I'm running in OpenFaaS. I want the device to be as "standalone" as possible but since I have some designs that do or don't have a GPS/GLONASS onboard I might just do that work by calling into an endpoint.

And now that I think about it - since the mobile device would likely be providing the lat/lon to the Meadow in cases where there's no onboard location capability I might also shift it to the mobile to do that work and simply hand it off (right now working on BLE services and could add that in). Anyway - this is a long LONG way of saying thank you again - this is really great work. πŸ†

houstonhaynes commented 1 year ago

FYI if you're interested heres' my F# refactor for the sample you provided. This is an .fsx file that can be run in F# interactive. Pretty boss! Thanks again!

#r "nuget: CoordinateSharp"

open System

open CoordinateSharp

let find_Closest_TimeUnit addTimeFn rangeStart rangeEnd (c : Coordinate) azimuth adjustNegatively =
    let el = EagerLoad(EagerLoadType.Celestial)
    el.Extensions <- EagerLoad_Extensions(EagerLoad_ExtensionsType.Solar_Cycle)
    let mutable closestTime = 0
    let mutable minDiff = 999.0

    for direction in [1; -1] do
        if direction = 1 || adjustNegatively then
            for x in rangeStart .. rangeEnd do
                let newTime = addTimeFn (float (x * direction))
                let nc = Coordinate(c.Latitude.ToDouble(), c.Longitude.ToDouble(), newTime, el)
                nc.Offset <- c.Offset
                let currentAzimuth = nc.CelestialInfo.SunAzimuth
                let diff = Math.Abs(currentAzimuth - azimuth)
                if diff < minDiff then
                    minDiff <- diff
                    closestTime <- x * direction

    addTimeFn (float closestTime)

let get_Time_At_Solar_Azimuth azimuth (c : Coordinate) =
    let hour = find_Closest_TimeUnit (c.GeoDate.AddHours) 0 23 c azimuth false
    let minute = find_Closest_TimeUnit (hour.AddMinutes) 0 59 c azimuth true
    let second = find_Closest_TimeUnit (minute.AddSeconds) 0 59 c azimuth true
    second

let main () =
    let lateAzimuth = 202.00
    let earlyAzimuth = 158.00
    let lat = 35.0
    let lon = -82.5
    let offset = -5.0
    let d = DateTime(2024, 1, 1)
    let c = Coordinate(lat, lon, d)
    c.Offset <- offset
    let timeAtEarlyAzimuth = get_Time_At_Solar_Azimuth earlyAzimuth c
    let timeAtLateAzimuth = get_Time_At_Solar_Azimuth lateAzimuth c
    printfn "Time at Sunrise: %A" c.CelestialInfo.SunRise
    printfn "Time at early azimuth: %A" timeAtEarlyAzimuth
    printfn "Solar noon: %A" c.CelestialInfo.SolarNoon
    printfn "Time at late azimuth: %A" timeAtLateAzimuth
    printfn "Time at Sunset: %A" c.CelestialInfo.SunSet
    if c.CelestialInfo.SolarNoon.HasValue then
        let ts = timeAtLateAzimuth - c.CelestialInfo.SolarNoon.Value
        printfn "Time range from Solar noon: %d:%d" ts.Hours (abs ts.Minutes)
    let c = Coordinate(lat, lon, timeAtEarlyAzimuth)
    c.Offset <- offset
    printfn "Early azimuth at time verification: %A" c.CelestialInfo.SunAzimuth

main ()
Tronald commented 1 year ago

Thanks @houstonhaynes

We actually have the logic built and ready to go for this feature, and are mainly awaiting implementing unit tests. We have a large project ongoing at the moment however, which has unfortunately feature delayed updates for CoordinateSharp for a bit. We are still rushing any emergency fixes that are identified still however. Hoping to have something out soon though.

houstonhaynes commented 1 year ago

That's great. I really enjoyed it and took it as an exercise to understand the API - and am really impressed with what's under the hood.

FWIW I gave you a shout-out in the "Amplifying F#" channel as well as follow up posts after the post record was available for "skunkworks" presentation. 😊

Tronald commented 11 months ago

Officially added with 2.20.1.1 pushed today. Feature will return a nullable DateTime that must be checked. This is because time at azimuth is not accurate if the sun is down. Also, if an azimuth doesn't exist during the day specified, a time should not be returned.

You may pass a specified azimuth delta however to set an allowable error tolerance. So lets say the minimum azimuth for the day is 23, and you set a delta of 1.5 degrees, an azimuth of 21.5 would still return a time. The greater the error delta, the greater the chance for inaccurate information however. The default error delta is 1 degree, but this may be reduced with the constructor.

//Create a coordinate and specify a date.
Coordinate c = new Coordinate(49, -122, new DateTime(2023, 9, 30));

//Set local UTC offset as desired.
c.Offset = -7;

//Set current sun azimuth in degrees E of N
double az = 120;

//Determine time of day. Default azimuth accuracy error delta is 1 degree by default, 
//but it is set at .5 for this example.
DateTime? t = Celestial.Get_Time_At_Solar_Azimuth(az, c, .5);

Console.WriteLine($"{t}"); //9/30/2023 9:21:44 AM