%% Marine Navigation Using Planetary Ephemerides
% This example shows how to use the planetary ephemerides and a Earth
% Centered Inertial to Earth Centered Earth Fixed (ECI to ECEF)
% transformation to perform celestial navigation of a marine vessel.
%
% This example uses the Mapping Toolbox(TM). You must also download data
% for the example using the <matlab:aeroDataPackage aeroDataPackage>
% command.
%
% <<astCelestialNavigationExampleMap.jpg>>
%
% This example uses the route followed by the 1947 expedition across the
% Pacific Ocean Kon-Tiki. The expedition, led by Thor Heyerdahl, aimed to
% prove the theory that the Polynesian islands were populated by people
% from South America in pre-Columbian times. The expedition took 101 days
% and sailed from the port of Callao, Peru to the Raroia atoll, French
% Polynesia. 
%
% <<astCelestialNavigationExampleVessel.jpg>>
%
% *Notes:* This example loosely recreates the expedition route. It takes
% some liberties to show the planetary ephemerides and ECI to ECEF
% transformation simply.

% Copyright 2012-2019 The MathWorks, Inc.

%% Load Vessel Track
% 
% Load the astKonTikiData.mat file. It contains the ship trajectory,
% velocity and course for this example. This file stores, the latitude and
% longitude for the different track points in the variables "lat" and
% "long", respectively. The variables contain enough data for one track
% point per day from the port of Callao to the Raroia atoll. Additionally,
% this file also stores values for each day's vessel velocity in knots per
% day "V" and course in deg "T".

load astKonTikiData

%% Create an Observation Structure
%
% The nautical reduction process is a series of steps that a navigator
% follows to determine the latitude and longitude of his vessel. It is
% based on the theory described in The American Practical Navigator[1], the
% Nautical Almanac[2], and the Explanatory Supplement to the Astronomical
% Almanac[3]. The process uses observational data obtained from a sextant,
% clock, compass, and navigational charts. It returns the intercept (p) and
% true azimuth (Z) for each of the observed objects. This example uses an
% observation structure array, obs, to contain the observational data. The
% fields for the structure array are:
%
% * h:              Height of eye level of the observer, in m.
% * IC:             Sextant's index correction, in deg.
% * P:              Local ambient pressure, in mb.
% * T:              Local temperature, in C.
% * year:           Local year at the time of the observation.
% * month:          Local month at the time of the observation.
% * day:            Local day at the time of the observation.
% * hour:           Local hour at the time of the observation.
% * UTC:            Coordinated Universal Time for the observation,
%                   represented as a six element vector with year, month,
%                   day, hour, minutes, and seconds.
% * Hs:             Sextant altitude of the celestial object above the
%                   horizon, in deg.
% * object:         Celestial object used for the measurement (i.e.,
%                   Jupiter, Neptune, Saturn, etc). 
% * latitude:       Estimated latitude of the vessel at the time of the 
%                   observation, in deg.
% * longitude:      Estimated longitude of the vessel at the time of the 
%                   observation, in deg.
% * declination:    Declination of the celestial object, in deg.
% * altitude:       Distance from the surface of the earth to the celestial
%                   object, in km.
% * GHA:            Greenwich Hour Angle, which is the angle in degrees of
%                   the celestial object relative to the Greenwich
%                   meridian.
%
% For simplicity, assume that all measurements are taken at the same
% location in the boat, with the same sextant, at the same ambient
% temperature and pressure: 

obs.h = 4; 
obs.IC = 0; 
obs.P = 982;
obs.T = 15;

%%
% The expedition departed on April 28th 1947. As a result, initialize the
% structure for the observation for this date:

obs.year = 1947;
obs.month = 4; 
obs.day = 28;

%% Initialize Dead Reckoning Process for Navigation
% To start the dead reckoning process, define the initial conditions for
% the position of the vessel. Store the latitude for a fixed solution for
% the latitude and longitude in the latFix  and longFix variables,
% respectively. In this example, for the first fix location, use the
% latitude and longitude of Callao, Peru:

longFix = zeros(size(long));
latFix = zeros(size(lat));
longFix(1) = long(1);
latFix(1) = lat(1);

%% Daily Dead Reckoning
% For this example, assume that a fix is obtained daily using the
% observational data. Therefore, this example uses a "for loop" for each
% observation. The variable "m" acts as a counter representing every
% elapsed day since the departure from the port:

for m = 1:size(lat,1)-1
    %%
    % Increment the day and make day adjustments for the months of June and
    % April, both of which have only 30 days:
    
    obs.day = obs.day + 1;
    [obs.month,obs.day] = astHelperDayCheck(obs.year,obs.month,obs.day);
    
    %% Actual Latitude and Longitude
    % Extract the vessel actual position for each day from the track points
    % loaded earlier. The example uses this value to calculate the local
    % time zone and the position of the selected planets in the sky:
    
    longActual = long(m+1);
    latActual = lat(m+1);
        
    %% Planet Selection
    % Select the planets for observation if they are visible to the
    % vessel for the given latitude and longitude. The following code uses
    % precalculated data:
    
    if longActual>-90
        obs.object = {'Saturn';'Neptune'};
    elseif longActual<=-90 && longActual>-95
        obs.object = {'Saturn';'Neptune';'Jupiter'}; 
    elseif longActual<=-95
        obs.object = {'Neptune';'Jupiter'};
    end
    
    %% UTC Time Calculation
    % Adjust local time to UTC depending on the assumed longitude. For this
    % example, assume that all observations are taken at the same time
    % every day at 8 p.m. local time.
    
    obs.hour = 20;
    
    %%
    % For the dead reckoning process, update the observation structure with
    % the estimate of the current location. In this case, the location is
    % estimated using the previous fix, the vessel's velocity V, and
    % course T.

    obs.longitude = longFix(m)+(1/60)*V(m)*sind(T(m))/cosd(latFix(m));
    obs.latitude = latFix(m)+(1/60)*V(m)*cosd(T(m));

    
    %%
    % Adjust the local time to UTC using the helper function
    % astHelperLongitudeHour. This function adjusts the UTC observation
    % time depending on the estimated longitude of the vessel.
    
    obs.UTC = astHelperLongitudeHour(obs);
    
    %% Sextant Altitude Calculation
    % For each of the planets, the astHelperNauticalCalculation helper
    % function calculates the sextant measurement that the crew in the
    % Kon-Tiki would have measured. This function models the actual
    % behavior of the planets, while compensating for the local conditions.
    % This function uses the planetary ephemeris and the ECI to ECEF
    % transformation matrix. The analysis doesn't include planetary
    % aberration, gravitational light deflection, and aberration of light
    % phenomena.
    
    obs.Hs = astHelperNauticalCalculation(obs,latActual,longActual);
          
    %% Calculate position
    % The following calculations replace the use of the nautical almanac.
    % They include the use of the planetary ephemerides and the ECI to ECEF
    % transformation matrix.
    %
    % Initialize declination, Greenwich Hour Angle (GHA), and altitude for
    % the observed object:
    
    obs.declination = zeros(size(obs.Hs));
    obs.GHA = zeros(size(obs.Hs));
    obs.altitude = zeros(size(obs.Hs));
    
    %%
    % Calculate the modified Julian date for the measurement time:
    
    mjd = mjuliandate(obs.UTC);
    
    %%
    % Calculate the difference between UT1 and UTC:
    
    dUT1 = deltaUT1(mjd,'Action','None');
    
    %%
    % Calculate the ECI to ECEF transformation matrix using the values of
    % TAI-UTC (dAT) from the U.S. Naval Observatory:
    
    dAT = 1.4228180;
    TM = dcmeci2ecef('IAU-76/FK5',obs.UTC,dAT,dUT1);
    
    %%
    % Calculate Julian date for Terrestrial Time to approximate the
    % Barycentric Dynamical Time:
    jdTT = juliandate(obs.UTC)+(dAT+32.184)/86400;
    
    %%
    % Calculate declination, Greenwich Hour Angle, and altitude for each of
    % the celestial objects:
    
    for k =1:length(obs.object)
        %%
        % Calculate the ECI position for every planet:
        
        posECI = planetEphemeris(jdTT,'Earth',obs.object{k},'405','km');
        
        %%
        % Calculate the ECEF position of every planet:
        
        posECEF = TM*posECI';

        %%
        % Calculate the Greenwich Hour Angle (GHA) and declination using
        % the ECEF position:
        
        obs.GHA(k) = -atan2d(posECEF(2),posECEF(1));
        obs.declination(k) = atan2d(posECEF(3),sqrt(posECEF(1)^2 + ...
            posECEF(2)^2));
        %%
        % Calculate the distance from the surface of the Earth to the
        % center of the planet using the ECEF to LLA transformation
        % function:
        
        posLLA = ecef2lla(1000*posECEF');
        obs.altitude(k) = posLLA(3);
 
    end
    %% Sight Reduction for Planetary Objects
    % Reduce sight for each of the planets specified in the observation
    % structure array:
        
    [p,Z] = astHelperNauticalReduction(obs);
    
    %%
    % Calculate the increments to the latitude and longitude from the last
    % fix for the current fix using the following equations. These
    % equations are based on the Nautical Almanac.
       
    Ap = sum(cosd(Z).^2);
    Bp = sum(cosd(Z).*sind(Z));
    Cp = sum(sind(Z).^2);
    Dp = sum(p.*cosd(Z));
    Ep = sum(p.*sind(Z));
    Gp = Ap*Cp-Bp^2;
    
    %%
    % Calculate the increments for latitude and longitude according to the
    % reduction:
    
    deltaLongFix = (Ap*Ep-Bp*Dp)/(Gp*cosd(latFix(m)));
    deltaLatFix = (Cp*Dp-Bp*Ep)/Gp;
    
    %%
    % After the increments for the latitude and longitude are calculated,
    % add them to the estimated location, obtaining the fix for the time of
    % observation:
      
    longFix(m+1) = obs.longitude + deltaLongFix;
    latFix(m+1) = obs.latitude + deltaLatFix;

end

%% Navigation Solution Visualization
% The following figure shows the actual track and sight reduction
% solutions:

astHelperVisualization(long,lat,longFix,latFix,'Plot')

%%
% You can use the Mapping Toolbox to obtain a more detailed graph
% depicting the sight reduction solutions with the American continent and
% French Polynesia.

astHelperVisualization(long,lat,longFix,latFix,'Map')

%%
% The relative error in longitude and latitude is accumulated as the vessel
% sails from Callao to Raroia. This error is due to small errors in the
% measurement of the sextant altitude. For June 9th, the reduction method
% calculates a true azimuth (Z) for Neptune and Jupiter. The true azimuths
% for Neptune and Jupiter are close to 180 deg apart. This difference
% causes a small peak in the relative error. This error, however, is still
% within the reduction method error boundaries.

%% References
% [1] Bowditch, N. The American Practical Navigator. National Geospatial
% Intelligence Agency, 2012.
%%
% [2] United Kingdom Hydrographic Office. Nautical Almanac 2012 Commercial
% Edition. Paradise Publications, Inc. 2011.
%%
% [3] Urban, Sean E. and P. Kenneth Seidelmann. Explanatory Supplement to
% the Astronomical Almanac. 3rd Ed., University Science Books, 2013.
%%
% [4] United States Naval Observatory.
%
displayEndOfDemoMessage(mfilename)
