[postgis-users] Solar azimuth calculator

Sat Apr 17 07:11:15 PDT 2021

Here it is in psql.  Doesn't need PostGIS.   Fun little exercise.

CREATE OR REPLACE FUNCTION public.solarazelev(

    obs_time timestamp with time zone,

    latitude numeric,

    longitude numeric,

    alt numeric)

    RETURNS TABLE(azimuth numeric, elevation numeric)

    LANGUAGE 'plpgsql'

AS $BODY$

DECLARE copyright varchar :=

' Programed by Darin C.Koblick 2 / 17 / 2009

              Darin C.Koblick 4 / 16 / 2013 Vectorized for Speed

                                         Allow for MATLAB Datevec input in

                                         addition to a UTC string.

                                         Cleaned up comments and code to

                                         avoid warnings in MATLAB editor.

                Kevin Godden 9/1/2020      Ported from Matlab to C++, tried
to change as little as possible.

                                         this is a non-vectorised port.

                Bruce Rindahl 4/17/2021        Ported to PostgreSQL

--------------------------------------------------------------------------

Copyright(c) 2010, Darin Koblick

Redistribution and use in source and binary forms, with or without

modification, are permitted provided that the following conditions are met :

*Redistributions of source code must retain the above copyright notice, this

list of conditions and the following disclaimer.

* Redistributions in binary form must reproduce the above copyright notice,

this list of conditions and the following disclaimer in the documentation

and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"

AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE

DISCLAIMED.IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE

FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

DAMAGES(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR

SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY,

OR TORT(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE

OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.';

    DECLARE M_PI numeric := pi();

    DECLARE jd numeric := to_char(obs_time,'J')::numeric;

    DECLARE d numeric := jd - 2451543.5;

    -- Keplerian Elements for the Sun(geocentric)

    DECLARE w numeric := 282.9404 + 4.70935e-5*d; -- (longitude of
perihelion degrees)

    -- a = 1.000000; % (mean distance, a.u.)

    DECLARE e numeric := 0.016709 - 1.151e-9*d; -- (eccentricity)

    DECLARE M numeric := mod(356.0470 + 0.9856002585*d, 360.0); -- (mean
anomaly degrees)

    DECLARE L numeric := w + M; -- (Sun's mean longitude degrees)

    DECLARE oblecl numeric := 23.4393 - 3.563e-7*d; -- (Sun's obliquity of
the ecliptic)

    -- auxiliary angle

    DECLARE  EE numeric := M + (180 / M_PI)*e*sin(M*(M_PI / 180))*(1 +
e*cos(M*(M_PI / 180)));

    -- rectangular coordinates in the plane of the ecliptic(x axis toward
perhilion)

    DECLARE x numeric := cos(EE*(M_PI / 180)) - e;

    DECLARE y numeric := sin(EE*(M_PI / 180))*sqrt(1 - pow(e, 2));

    -- find the distance and true anomaly

    DECLARE r numeric := sqrt(pow(x,2) + pow(y,2));

    -- Is this y,x or x,y? opposite from Excel  - use C++ version

    DECLARE v numeric := atan2(y, x)*(180 / M_PI);

    -- find the longitude of the sun

    DECLARE lon numeric := v + w;

    -- compute the ecliptic rectangular coordinates

    DECLARE xeclip numeric := r*cos(lon*(M_PI / 180));

    DECLARE yeclip numeric := r*sin(lon*(M_PI / 180));

    DECLARE zeclip numeric := 0.0;

    --rotate these coordinates to equitorial rectangular coordinates

    DECLARE xequat numeric := xeclip;

    DECLARE yequat numeric := yeclip*cos(oblecl*(M_PI / 180)) + zeclip *
sin(oblecl*(M_PI / 180));

    DECLARE zequat numeric := yeclip*sin(23.4406*(M_PI / 180)) + zeclip *
cos(oblecl*(M_PI / 180));

    -- convert equatorial rectangular coordinates to RA and Decl:

    DECLARE RR numeric := sqrt(pow(xequat, 2) + pow(yequat, 2) +
pow(zequat, 2)) - (Alt / 149598000); -- roll up the altitude correction

    DECLARE RA numeric := atan2(yequat, xequat)*(180 / M_PI);

    DECLARE delta numeric := asin(zequat / RR)*(180 / M_PI);

    DECLARE UTH numeric := date_part('hour',obs_time)

        + date_part('minute',obs_time)/60

        + date_part('second',obs_time) /3600

        - extract( timezone from now())/3600;

    -- Calculate local siderial time

    DECLARE GMST0 numeric := mod(L + 180, 360.0) / 15;

    DECLARE SIDTIME numeric := GMST0 + UTH + Longitude / 15;

    -- Replace RA with hour angle HA

    DECLARE HA numeric := (SIDTIME*15 - RA);

    -- convert to rectangular coordinate system

    DECLARE XX numeric := cos(HA*(M_PI / 180))*cos(delta*(M_PI / 180));

    DECLARE YY numeric := sin(HA*(M_PI / 180))*cos(delta*(M_PI / 180));

    DECLARE z numeric := sin(delta*(M_PI / 180));

    -- rotate this along an axis going east - west.

    DECLARE xhor numeric := XX*cos((90 - Latitude)*(M_PI / 180)) -
z*sin((90 - Latitude)*(M_PI / 180));

    DECLARE yhor numeric := YY;

    DECLARE zhor numeric := XX*sin((90 - Latitude)*(M_PI / 180)) +
z*cos((90 - Latitude)*(M_PI / 180));

    -- Find the h and AZ

    DECLARE Az numeric := round((atan2(yhor, xhor)*(180 / M_PI) +
180)::numeric,3);

    DECLARE El numeric := round((asin(zhor)*(180 / M_PI))::numeric,3);

BEGIN

RETURN QUERY select Az,El;

END

$BODY$;

On Thu, Apr 15, 2021, 1:46 PM Bruce Rindahl <bruce.rindahl at gmail.com> wrote:

> Here is the code for it in c:
>
>
> https://www.ridgesolutions.ie/index.php/2020/01/14/c-code-to-estimate-solar-azimuth-and-elevation-given-gps-position-and-time/
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210417/0e7e1b7f/attachment.html>