search

kylebarron / suncalc-py

A Python port of suncalc.js for calculating sun position and sunlight phases
MIT License
60 stars 8 forks source link

Email bug report #6

Open kylebarron opened 3 years ago

kylebarron commented 3 years ago

From email:

I've been playing around with your suncalc module and think I have encountered a systematic error, here are the details and I've attached my program.

  1. Use get_times to get the 'solar_noon' data for the June solstice

  2. Feed that time into get_position to get Azimuth & Altitude data

The Azimuth at such a time should either be 0 or PI but that was not what I was seeing

Doing the same for Dec solstice didn't produce such an error [at least for the value of long & lat I was using to understand the issue - 0,0]

Eventually, I hypothesised that there might be a screwup in interpreting the time [by me or the function] and tried altering the hour of the solar noon as given and adjusting by +1 hour seemed to fix the problem. I ported the fix in to my original intention for the program [check difference between altitude angles at solar noon on solstices is a constant where ever you are] and could demonstrate that this is an issue at all the locations in the set of co-ords I was generating.

In the program you'll find fix_offset which is currently set to 1 and thus correctly generates correct values. Change this to 0 and the error manifests itself.

Python Version: 3.8.5

Suncalc: 0.1.2 installed via pip3 on Mar 17th

Code:

#!/usr/bin/env python3
from suncalc import get_position, get_times
from datetime import datetime
from collections import namedtuple
import random, math

Sun_dates = namedtuple('Sun_dates','march_equinox june_solstice sept_equinox dec_solstice')

sun_dates_2021 = Sun_dates(datetime(2021,3,20),datetime(2021,6,21),datetime(2021,9,22),datetime(2021,12,22))

fix_offset = 1 # added to solar noon time to correct systematic error in summer soltice positon?

for lat in range(5): #range(20,70,10):
    lon = random.randint(0,100)
    lat = random.randint(0,60)
    print(lon,lat)

    june_solar_noon = get_times(sun_dates_2021.june_solstice,lon,lat)['solar_noon']
    print(june_solar_noon)

    dec_solar_noon = get_times(sun_dates_2021.dec_solstice,lon,lat)['solar_noon']
    print(dec_solar_noon)

    june_pos = get_position(june_solar_noon.replace(hour=june_solar_noon.hour + fix_offset),lng=lon,lat=lat)
    print(june_pos,math.degrees(june_pos['azimuth']),math.degrees(june_pos['altitude']))

    dec_pos = get_position(dec_solar_noon,lng=lon,lat=lat)
    print(dec_pos,math.degrees(dec_pos['azimuth']),math.degrees(dec_pos['altitude']))

    if abs(june_pos['azimuth']) > 1:
        if abs(dec_pos['azimuth']) > 1:
            max_elev = max(june_pos['altitude'],dec_pos['altitude'])
            min_elev = min(june_pos['altitude'],dec_pos['altitude'])
            print(math.degrees(max_elev - min_elev))
        else: 
            print(math.degrees(math.pi - june_pos['altitude'] - dec_pos['altitude']))
    elif abs(dec_pos['azimuth']) > 1:
        print(math.degrees(math.pi - june_pos['altitude'] - dec_pos['altitude']))
    else:
        max_elev = max(june_pos['altitude'],dec_pos['altitude'])
        min_elev = min(june_pos['altitude'],dec_pos['altitude'])
        print(math.degrees(max_elev - min_elev))
    print('===')