Skip to content

Instantly share code, notes, and snippets.

@andrewxhill
Last active September 14, 2016 12:07
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save andrewxhill/b2b7b7da7dff71e17b7c to your computer and use it in GitHub Desktop.
Test of SunCalc in SQL
-- Test of porting SunCalc to SQL
-- Original SunCalc https://github.com/mourner/suncalc
-- Original SunCalc License https://github.com/mourner/suncalc/blob/master/LICENSE
CREATE TYPE suncalc_position AS (
azimuth decimal,
altitude decimal
);
CREATE TYPE suncalc_coords AS (
declination decimal,
rightAscension decimal
);
CREATE TYPE suncalc_positions AS (
name text,
time timestamptz
);
CREATE OR REPLACE FUNCTION SunCalc_ToJulian(date TIMESTAMPTZ) RETURNS numeric as $$
DECLARE
dayMs NUMERIC = 60 * 60 * 24;
J1970 NUMERIC = 2440588;
J2000 NUMERIC = 2451545;
BEGIN
RETURN (EXTRACT(EPOCH FROM date)) / dayMs - 0.5 + J1970;
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_FromJulian(j NUMERIC) RETURNS timestamptz as $$
DECLARE
dayMs NUMERIC = 60 * 60 * 24;
J1970 NUMERIC = 2440588;
J2000 NUMERIC = 2451545;
BEGIN
RETURN TIMESTAMP WITH TIME ZONE 'epoch' + (j + 0.5 - J1970) * dayMs * INTERVAL '1 second';
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_ToDays(date TIMESTAMPTZ) RETURNS NUMERIC as $$
DECLARE
dayMs NUMERIC = 60 * 60 * 24;
J1970 NUMERIC = 2440588;
J2000 NUMERIC = 2451545;
BEGIN
RETURN SunCalc_ToJulian(date) - J2000;
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_RightAscension(l NUMERIC, b NUMERIC) RETURNS NUMERIC as $$
DECLARE
rad NUMERIC = pi() / 180.0;
e NUMERIC = rad * 23.4397;
BEGIN
RETURN atan2(sin(l) * cos(e) - tan(b) * sin(e), cos(l));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_Declination(l NUMERIC, b NUMERIC) RETURNS NUMERIC as $$
DECLARE
rad NUMERIC = pi() / 180.0;
e NUMERIC = rad * 23.4397;
BEGIN
RETURN asin(sin(b) * cos(e) + cos(b) * sin(e) * sin(l));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_Azimuth(H NUMERIC, phi NUMERIC, deci NUMERIC) RETURNS NUMERIC as $$
BEGIN
RETURN atan2(sin(H), cos(H) * sin(phi) - tan(deci) * cos(phi));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_Altitude(H NUMERIC, phi NUMERIC, deci NUMERIC) RETURNS NUMERIC as $$
BEGIN
RETURN asin(sin(phi) * sin(deci) + cos(phi) * cos(deci) * cos(H));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_SiderealTime(d NUMERIC, lw NUMERIC) RETURNS NUMERIC as $$
DECLARE
rad NUMERIC = pi() / 180.0;
BEGIN
RETURN rad * (280.16 + 360.9856235 * d) - lw;
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_SolarMeanAnomaly(d NUMERIC) RETURNS NUMERIC as $$
DECLARE
rad NUMERIC = pi() / 180.0;
BEGIN
RETURN rad * (357.5291 + 0.98560028 * d);
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_EclipticLongitude(M NUMERIC) RETURNS NUMERIC as $$
DECLARE
rad NUMERIC = pi() / 180.0;
c NUMERIC;
p NUMERIC;
BEGIN
c = rad * (1.9148 * sin(M) + 0.02 * sin(2 * M) + 0.0003 * sin(3 * M));
p = rad * 102.9372;
RETURN M + c + p + pi();
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_SunCoords(d NUMERIC) RETURNS suncalc_coords as $$
DECLARE
m NUMERIC;
l NUMERIC;
BEGIN
m = SunCalc_SolarMeanAnomaly(d);
l = SunCalc_EclipticLongitude(M);
-- Returns (dec, ra)
RETURN (SunCalc_Declination(L, 0), SunCalc_RightAscension(L, 0));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_julianCycle(d NUMERIC, lw NUMERIC) RETURNS NUMERIC as $$
DECLARE
jO NUMERIC = 0.0009;
BEGIN
RETURN round(d - jO - lw / (2 * pi()));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_approxTransit(Ht NUMERIC, lw NUMERIC, n NUMERIC) RETURNS NUMERIC as $$
DECLARE
jO NUMERIC = 0.0009;
BEGIN
RETURN jO + (Ht + lw) / (2 * pi()) + n;
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_solarTransitJ(ds NUMERIC, M NUMERIC, L NUMERIC) RETURNS NUMERIC as $$
DECLARE
jO NUMERIC = 0.0009;
J2000 NUMERIC = 2451545;
BEGIN
RETURN J2000 + ds + 0.0053 * sin(M) - 0.0069 * sin(2 * L);
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_hourAngle(h NUMERIC, phi NUMERIC, d NUMERIC) RETURNS NUMERIC as $$
BEGIN
RETURN acos((sin(h) - sin(phi) * sin(d)) / (cos(phi) * cos(d)));
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_getSetJ(h NUMERIC, lw NUMERIC, phi NUMERIC, deci NUMERIC, n NUMERIC, M NUMERIC, L NUMERIC) RETURNS NUMERIC as $$
DECLARE
w NUMERIC;
a NUMERIC;
BEGIN
w = SunCalc_hourAngle(h, phi, deci);
a = SunCalc_approxTransit(w, lw, n);
RETURN SunCalc_solarTransitJ(a, M, L);
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_GetPosition(date TIMESTAMPTZ, coord GEOMETRY) RETURNS suncalc_position as $$
DECLARE
rad NUMERIC = pi() / 180.0;
lw NUMERIC;
phi NUMERIC;
d NUMERIC;
c suncalc_coords;
h NUMERIC;
azimuth NUMERIC;
altitude NUMERIC;
BEGIN
lw = rad * -1 * ST_X(coord);
phi = rad * ST_Y(coord);
d = SunCalc_ToDays(date);
c = SunCalc_SunCoords(d);
h = SunCalc_SiderealTime(d, lw) - c.rightAscension;
-- returns (azimuth, altitude)
RETURN (SunCalc_Azimuth(h, phi, c.declination), SunCalc_Altitude(h, phi, c.declination) );
END;
$$ language plpgsql IMMUTABLE;
CREATE OR REPLACE FUNCTION SunCalc_GetTimes(date TIMESTAMPTZ, coord GEOMETRY) RETURNS SETOF suncalc_positions as $$
DECLARE
rad NUMERIC = pi() / 180.0;
lw NUMERIC;
phi NUMERIC;
d NUMERIC;
n NUMERIC;
ds NUMERIC;
M NUMERIC;
L NUMERIC;
deci NUMERIC;
Jnoon NUMERIC;
i NUMERIC;
len NUMERIC;
tt NUMERIC;
Jset NUMERIC;
Jrise NUMERIC;
degs NUMERIC[] = ARRAY[-0.833, -0.3, -6, -12, -18, 6];
rises TEXT[] = ARRAY['sunrise','sunriseEnd','dawn','nauticalDawn','nightEnd','goldenHourEnd'];
sets TEXT[] = ARRAY['sunset','sunsetStart','dusk','nauticalDusk','night','goldenHour'];
BEGIN
lw = rad * -1 * ST_X(coord);
phi = rad * ST_Y(coord);
d = SunCalc_ToDays(date);
n = SunCalc_JulianCycle(d, lw);
ds = SunCalc_ApproxTransit(0, lw, n);
M = SunCalc_SolarMeanAnomaly(ds);
L = SunCalc_EclipticLongitude(M);
deci = SunCalc_Declination(L, 0);
Jnoon = SunCalc_SolarTransitJ(ds, M, L);
RETURN NEXT ('solarnoon'::text, SunCalc_fromJulian(Jnoon));
RETURN NEXT ('nadir'::text, SunCalc_fromJulian(Jnoon - 0.5));
FOR i IN 1 .. array_upper(degs, 1)
LOOP
Jset = SunCalc_getSetJ(degs[i] * rad, lw, phi, deci, n, M, L);
Jrise = Jnoon - (Jset - Jnoon);
RETURN NEXT (rises[i], SunCalc_fromJulian(Jrise));
RETURN NEXT (sets[i], SunCalc_fromJulian(Jset));
END LOOP;
RETURN;
END;
$$ language plpgsql IMMUTABLE;
@wojnicki
Copy link

Does not work for some timestamps, eg.:
select * from SunCalc_GetTimes('2015-06-01 13:00'::timestamp with time zone,ST_PointFromText('POINT(19.938333 50.061389)', 4326));

ERROR: input is out of range
CONTEXT: PL/pgSQL function suncalc_hourangle(numeric,numeric,numeric) line 3 at RETURN
PL/pgSQL function suncalc_getsetj(numeric,numeric,numeric,numeric,numeric,numeric,numeric) line 6 at assignment
PL/pgSQL function suncalc_gettimes(timestamp with time zone,geometry) line 39 at assignment

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment