#!/usr/bin/python """Solar position calculations, adapted from Suncalc, a JS lib by Vladimir Agafonkin (https://github.com/mourner/suncalc)""" from math import sin, cos, tan, asin, atan2, pi from datetime import datetime rad = pi / 180.0 epochStart = datetime(1970, 1, 1) J1970 = 2440588 J2000 = 2451545 dayMs = 24 * 60 * 60 * 1000 e = rad * 23.4397 # obliquity of the Earth def toMillis(date): return (date - epochStart).total_seconds() * 1000 def toJulian(date): return toMillis(date) / dayMs - 0.5 + J1970 def toDays(date): return toJulian(date) - J2000 # function fromJulian(j) { return new Date((j + 0.5 - J1970) * dayMs); } def rightAscension(l, b): return atan2(sin(l) * cos(e) - tan(b) * sin(e), cos(l)) def declination(l, b): return asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l)) def azimuth(H, phi, dec): return atan2(sin(H), cos(H) * sin(phi) - tan(dec) * cos(phi)) def altitude(H, phi, dec): return asin(sin(phi) * sin(dec) + cos(phi) * cos(dec) * cos(H)) def siderealTime(d, lw): return rad * (280.16 + 360.9856235 * d) - lw def solarMeanAnomaly(d): return rad * (357.5291 + 0.98560028 * d) def eclipticLongitude(M): C = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M)) # equation of center P = rad * 102.9372 # perihelion of the Earth return M + C + P + pi def sunCoords(d): M = solarMeanAnomaly(d) L = eclipticLongitude(M) return { 'dec': declination(L, 0), 'ra': rightAscension(L, 0) } def solarPosition(date, lat, lng): lw = rad * -lng phi = rad * lat d = toDays(date) c = sunCoords(d) H = siderealTime(d, lw) - c['ra'] return { 'azimuth': azimuth(H, phi, c['dec']), 'altitude': altitude(H, phi, c['dec']) } def getFinalSolarPosition(latitude, longitude, currentTime=datetime.utcnow()): position = solarPosition(currentTime, latitude, longitude) azimuth = position['azimuth'] * 180.0 / pi altitude = position['altitude'] * 180.0 / pi pos = { 'azimuth': round(azimuth + 180.0, 1), 'elevation': round(altitude, 1) } return pos def main(): latitude = -37.94514 longitude = 145.07791 now = datetime.utcnow() position = solarPosition(now, latitude, longitude) azimuth = position['azimuth'] * 180.0 / pi altitude = position['altitude'] * 180.0 / pi print(180 + azimuth) print(altitude) if __name__ == '__main__': main()