Code covered by the BSD License  

Highlights from
Lunar Azimuth and Altitude Estimation Algorithm

image thumbnail

Lunar Azimuth and Altitude Estimation Algorithm

by

Darin Koblick (view profile)

 

15 Feb 2009 (Updated )

Predict the Lunar Azimuth and Altitude within +/- .2 deg of any lat and lon for a given UTC

LunarAzEl(UTC,Lat,Lon,Alt)
function [Az h] = LunarAzEl(UTC,Lat,Lon,Alt)

% Programed by Darin C. Koblick 2/14/2009
%
% Updated on 03/04/2009 to clean up code and add quadrant check to Azimuth
% Thank you Doug W. for your help with the test code to find the quadrant check 
% error.
%
% Updated on 04/13/2009 to add Lunar perturbation offsets (this will
% increase the accuracy of calculations)
%
% Updated on 08/17/2010 to make use of the site altitude (this will affect
% the elevation angle)

% External Function Call Sequence:
% [Az El] = LunarAzEl('1991/05/19 13:00:00',50,10,0)


% Function Description:
% LunarAzEl will ingest a Universal Time, and specific site location on earth
% it will then output the lunar Azimuth and Elevation angles relative to that
% site.

% External Source References:

% Basics of Positional Astronomy and Ephemerides
% http://jgiesen.de/elevazmoon/basics/index.htm

% Computing planetary positions - a tutorial with worked examples
% http://stjarnhimlen.se/comp/tutorial.html

%Input Description:
% UTC (Coordinated Universal Time YYYY/MM/DD hh:mm:ss)
% Lat (Site Latitude in degrees -90:90 -> S(-) N(+))
% Lon (Site Longitude in degrees -180:180 W(-) E(+))
% Altitude of the site above sea level (km)

%Output Description:
%Az (Azimuth location of the moon in degrees)
%El (Elevation location of the moon in degrees)

%Verified output by comparison with the following source data:
%http://aa.usno.navy.mil/data/docs/AltAz.php

% Code Sequence
%--------------------------------------------------------------------------
%Do initial Longitude Latitude check

while Lon > 180
   Lon = Lon - 360;
end
while Lon < -180
   Lon = Lon + 360; 
end
while Lat > 90
    Lat = Lat - 360;
end
while Lat < -90
   Lat = Lat + 360; 
end


%Declare Earth Equatorial Radius Measurements in km
EarthRadEq = 6378.1370;

%Convert Universal Time to Ephemeris Time
jd = juliandate(UTC,'yyyy/mm/dd HH:MM:SS');

%Find the Day Number
d = jd - 2451543.5;

%Keplerian Elements of the Moon
%This will also account for the Sun's perturbation
 N = 125.1228-0.0529538083.*d; %    (Long asc. node deg)
 i = 5.1454; %                      (Inclination deg)
 w = 318.0634 + 0.1643573223.*d; %  (Arg. of perigee deg)
 a =  60.2666;%                     (Mean distance (Earth's Equitorial Radii)
 e = 0.054900;%                     (Eccentricity)
 M = mod(115.3654+13.0649929509.*d,360);%    (Mean anomaly deg)
  
 LMoon =  mod(N + w + M,360);                 %(Moon's mean longitude deg)
 FMoon =  mod(LMoon - N,360);                 %(Moon's argument of latitude)

 %Keplerian Elements of the Sun
 wSun = mod(282.9404 + 4.70935E-5.*d,360);    % (longitude of perihelion)
 MSun = mod(356.0470 + 0.9856002585.*d,360);  % (Sun mean anomaly)
 LSun = mod(wSun + MSun,360);                 % (Sun's mean longitude)
     
 DMoon =  LMoon - LSun;                     % (Moon's mean elongation)  


 %Calculate Lunar perturbations in Longitude
 LunarPLon = [ -1.274.*sin((M - 2.*DMoon).*(pi/180)); ...
     .658.*sin(2.*DMoon.*(pi/180)); ...
     -0.186.*sin(MSun.*(pi/180)); ...
     -0.059.*sin((2.*M-2.*DMoon).*(pi/180)); ...
     -0.057.*sin((M-2.*DMoon + MSun).*(pi/180)); ...
     .053.*sin((M+2.*DMoon).*(pi/180)); ...
     .046.*sin((2.*DMoon-MSun).*(pi/180)); ...
     .041.*sin((M-MSun).*(pi/180)); ...
    -0.035.*sin(DMoon.*(pi/180)); ...           
    -0.031.*sin((M+MSun).*(pi/180)); ...
    -0.015.*sin((2.*FMoon-2.*DMoon).*(pi/180)); ...
    .011.*sin((M-4.*DMoon).*(pi/180))];
 
 %Calculate Lunar perturbations in Latitude 
 LunarPLat = [ -0.173.*sin((FMoon-2.*DMoon).*(pi/180)); ...
    -0.055.*sin((M-FMoon-2.*DMoon).*(pi/180)); ...
    -0.046.*sin((M+FMoon-2.*DMoon).*(pi/180)); ...
    +0.033.*sin((FMoon+2.*DMoon).*(pi/180)); ...
    +0.017.*sin((2.*M+FMoon).*(pi/180))];

%Calculate perturbations in Distance
 LunarPDist = [ -0.58*cos((M-2.*DMoon).*(pi/180)); ...
    -0.46.*cos(2.*DMoon.*(pi/180))];

% Compute E, the eccentric anomaly

%E0 is the eccentric anomaly approximation estimate 
%(this will initially have a relativly high error)
E0 = M+(180./pi).*e.*sin(M.*(pi/180)).*(1+e.*cos(M.*(pi/180)));

%Compute E1 and set it to E0 until the E1 == E0
E1 = E0-(E0-(180/pi).*e.*sin(E0.*(pi/180))-M)./(1-e*cos(E0.*(pi/180)));
while E1-E0 > .000005
    E0 = E1;
    E1 = E0-(E0-(180/pi).*e.*sin(E0.*(pi/180))-M)./(1-e*cos(E0.*(pi/180)));    
end
E = E1;

%Compute rectangular coordinates (x,y) in the plane of the lunar orbit
x = a.*(cos(E.*(pi/180))-e);
y = a.*sqrt(1-e.*e).*sin(E.*(pi/180));

%convert this to distance and true anomaly
r = sqrt(x.*x + y.*y);
v = atan2(y.*(pi/180),x.*(pi/180)).*(180/pi);

%Compute moon's position in ecliptic coordinates
xeclip = r.*(cos(N.*(pi/180)).*cos((v+w).*(pi/180))-sin(N.*(pi/180)).*sin((v+w).*(pi/180)).*cos(i.*(pi/180)));
yeclip = r.*(sin(N.*(pi/180)).*cos((v+w).*(pi/180))+cos(N.*(pi/180))*sin(((v+w).*(pi/180)))*cos(i.*(pi/180)));
zeclip = r.*sin((v+w).*(pi/180)).*sin(i.*(pi/180));

%Add the calculated lunar perturbation terms to increase model fidelity
[eLon eLat eDist] = cart2sph(xeclip,yeclip,zeclip);
[xeclip yeclip zeclip] = sph2cart(eLon + sum(LunarPLon).*(pi/180), ...
                                  eLat + sum(LunarPLat).*(pi/180), ...
                                  eDist + sum(LunarPDist));
clear eLon eLat eDist;
                              
%convert the latitude and longitude to right ascension RA and declination
%delta
T = (jd-2451545.0)/36525.0;

%Generate a rotation matrix for ecliptic to equitorial
%RotM=rotm_coo('E',jd);
%See rotm_coo.m for obl and rotational matrix transformation
Obl = 23.439291 - 0.0130042.*T - 0.00000016.*T.*T + 0.000000504.*T.*T.*T;
Obl = Obl.*(pi/180);
RotM = [1 0 0; 0 cos(Obl) sin(Obl); 0 -sin(Obl) cos(Obl)]';

%Apply the rotational matrix to the ecliptic rectangular coordinates
%Also, convert units to km instead of earth equatorial radii
sol = RotM*[xeclip yeclip zeclip]'.*EarthRadEq;

%Find the equatorial rectangular coordinates of the location specified
[xel yel zel] = sph2cart(Lon.*(pi/180),Lat.*(pi/180),Alt+EarthRadEq);

%Find the equatorial rectangular coordinates of the location @ sea level
[xsl ysl zsl] = sph2cart(Lon.*(pi/180),Lat.*(pi/180),EarthRadEq);

%Find the Angle Between sea level coordinate vector and the moon vector
theta1 = 180 - acosd(dot([xsl ysl zsl],[sol(1)-xsl sol(2)-ysl sol(3)-zsl]) ...
        ./(sqrt(xsl.^2 + ysl.^2 + zsl.^2) ...
         .*sqrt((sol(1)-xsl).^2 + (sol(2)-ysl).^2 + (sol(3)-zsl).^2)));

%Find the Angle Between the same coordinates but at the specified elevation
theta2 = 180 - acosd(dot([xel yel zel],[sol(1)-xel sol(2)-yel sol(3)-zel]) ...
    ./(sqrt(xel.^2 + yel.^2 + zel.^2) ...
         .*sqrt((sol(1)-xel).^2 + (sol(2)-yel).^2 + (sol(3)-zel).^2)));
     
%Find the Difference Between the two angles (+|-) is important
thetaDiff = theta2 - theta1;

% equatorial to horizon coordinate transformation
 [RA,delta] = cart2sph(sol(1),sol(2),sol(3));
 delta = delta.*(180/pi);
 RA = RA.*(180/pi);
 
%Following the RA DEC to Az Alt conversion sequence explained here:
%http://www.stargazing.net/kepler/altaz.html

%Find the J2000 value
J2000 = jd - 2451545.0;
hourvec = datevec(UTC,'yyyy/mm/dd HH:MM:SS');
UTH = hourvec(4) + hourvec(5)/60 + hourvec(6)/3600;

%Calculate local siderial time
LST = mod(100.46+0.985647.*J2000+Lon+15*UTH,360);

%Replace RA with hour angle HA
HA = LST-RA;

%Find the h and AZ at the current LST
h = asin(sin(delta.*(pi/180)).*sin(Lat.*(pi/180)) + cos(delta.*(pi/180)).*cos(Lat.*(pi/180)).*cos(HA.*(pi/180))).*(180/pi);
Az = acos((sin(delta.*(pi/180)) - sin(h.*(pi/180)).*sin(Lat.*(pi/180)))./(cos(h.*(pi/180)).*cos(Lat.*(pi/180)))).*(180/pi);

%Add in the angle offset due to the specified site elevation
h = h + thetaDiff;

if sin(HA.*(pi/180)) >= 0
   Az = 360-Az; 
end

%Apply Paralax Correction if we are still on earth
if Alt < 100
    horParal = 8.794/(r*6379.14/149.59787e6);
    p = asin(cos(h.*(pi/180))*sin((horParal/3600).*(pi/180))).*(180/pi);
    h = h-p;
end

function jd = juliandate(varargin)
% This sub function is provided in case juliandate does not come with your 
% distribution of Matlab

[year month day hour min sec] = datevec(datenum(varargin{:}));

for k = length(month):-1:1
    if ( month(k) <= 2 ) % january & february
        year(k)  = year(k) - 1.0;
        month(k) = month(k) + 12.0;
    end
end

jd = floor( 365.25*(year + 4716.0)) + floor( 30.6001*( month + 1.0)) + 2.0 - ...
    floor( year/100.0 ) + floor( floor( year/100.0 )/4.0 ) + day - 1524.5 + ...
    (hour + min/60 + sec/3600)/24;












Contact us