[postgis-users] Solar azimuth calculator
Alexander Gataric
gataric at usa.net
Sat Apr 17 11:57:17 PDT 2021
Thanks!
Get BlueMail for Android
On Apr 17, 2021, 9:12 AM, at 9:12 AM, Bruce Rindahl <bruce.rindahl at gmail.com> wrote:
>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
>
>All rights reserved.
>
>
>
>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/
>>
>>
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>postgis-users mailing list
>postgis-users at lists.osgeo.org
>https://lists.osgeo.org/mailman/listinfo/postgis-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210417/d5300654/attachment.html>
More information about the postgis-users
mailing list