function gr_mercury_precession(varargin) %#codegen

if nargin > 0
  dt = varargin{1};
else
  dt = 0.5;   % Time step (for standard iterations)
end
n = 100;  % Trail history (visual effect only)

coder.extrinsic('plot3', 'findobj', 'set', 'drawnow', 'sprintf', 'fprintf');

%   Copyright 2010-2013 The MathWorks, Inc.

c = 299792458;          % Speed of light
G = 6.6742e-11;         % Newton's gravitational constant

mr = 6.982e10;          % Mercury max radius
mv = 38.86e3/mr;        % Mercury max velocity (in units of c)

init; % Initialize the figure windows

first = true;

mass = 1.9891e30;
geometric_mass = G*mass/c/c;

mercury = zeros(1,8);
mercury(1) = 0;             % Internal clock
mercury(2) = mr/c;          % Maximum radius in units of c
mercury(3) = pi/2;          % Ecliptic latitude
mercury(4) = 0;             % Ecliptic longitude
mercury(5) = 1/sqrt(1-mv*mv); % Speed of clock
mercury(6) = 0;             % Delta radius = radial speed
mercury(7) = 0;             % Delta theta = ecliptic latitude speed
mercury(8) = mv;            % Angular velocity

% To plot a trail of dots for the visual effect
dots_x = zeros(n,1);
dots_y = zeros(n,1);
dots_z = zeros(n,1);

internal_cnt = int32(65536*2);
index = int32(0);

% History matrix. Each column captures the state of Mercury for a single time step.
hist = zeros(numel(mercury),internal_cnt);
t = 0;
lastProgress = 0;
while true
    % Initialize history matrix
    for i = 1:internal_cnt
        mercury = gr_trajectory(geometric_mass/c,mercury,dt);
        hist(:,i) = mercury(:);
    end
    % Get the current angle/location of the planet
    phi = mercury(4);
    progress = round(100 * phi / (2*pi));
    if progress ~= lastProgress && (mod(progress,5) == 0 || abs(lastProgress-progress) >= 5)
        fprintf(1, 'Progress: %d%%\n', progress);
        lastProgress = progress;
    end
    
    % Have we tried to complete a revolution?
    if phi > 2*pi
        % Reduce time step
        dt = dt / 2;
        % If time step is too big then reset the state of the planet
        % using its history. Otherwise we're trying to find the exact point
        % where the planet is doing a complete revolution.
        if dt > 0.1
            mercury(:) = hist(:,1);
            continue;
        end
    end

    % Have surpassed a complete revolution?
    if phi > 2*pi
        % Find the maximum distance/radius away from sun in the history
        % matrix
        [~,k] = max(hist(2,:));
        % Get the angle/location of the planet at this distance
        phi = hist(4,k);
        % Normalize angle w.r.t. revolutions
        while phi > 2*pi
            phi = phi - 2*pi;
        end
        % Update time to this position
        t = t + dt*k;
        years = int32(floor(t/3600/24/365));
        days = t/3600/24 - double(years*365);
        prec = (100/(t/3600/24/365))*(phi/2/pi)*360*3600;
        fprintf('precession: %3.5f" (%d years %3.5f days) => %3.3f"/century\n', ...
                (phi/2/pi)*360*3600, years, days, prec);
        return;
    end
    
    % Update universal clock
    t = t + dt*double(internal_cnt);
    
    % Update trail of dots
    if index == numel(dots_x)
        for i = 1:index-1
            dots_x(i) = dots_x(i+1);
            dots_y(i) = dots_y(i+1);
            dots_z(i) = dots_z(i+1);
        end
    else
        index = index + 1;
    end

    radius = mercury(2);
    theta = mercury(3);
    phi = mercury(4);
    
    % Add current position of planet to the trail of dots
    dots_x(index) = radius*sin(theta)*cos(phi);
    dots_y(index) = radius*sin(theta)*sin(phi);
    dots_z(index) = radius*cos(theta);

    % First iteration? Initialize the plot.
    if first
        first = false;
        plot3(0,0,0,'y*','Tag','mercury_plot_obj1');
        plot3(dots_x,dots_y,dots_z,'w-','Tag','mercury_plot_obj2');
    end
    
    plotObj2 = findobj('Tag','mercury_plot_obj2');

    % Update plot data
    set(plotObj2,'XData', dots_x);
    set(plotObj2,'YData', dots_y);
    set(plotObj2,'ZData', dots_z);

    % Draw display string with current information
    fig = findobj('Tag','mercury_plot');
    years = int32(floor(t/3600/24/365));            
    days = t/3600/24 - double(years*365);
    days0 = int32(days);
    coder.varsize('str');
    str = ' ';
    str = repmat(str,1,256);
    str = sprintf('mercury_plot: %3d years %3d days', years, days0);
    set(fig, 'Name', str); 

    drawnow;

    % Avoid non termination. It shouldn't take this long for planet Mercury
    % to complete a revolution.
    if years >= 100
        break
    end
end

% Sub-function to initialize the figure windows.
function init
coder.extrinsic('findobj', 'figure', 'clf', 'set', 'hold', 'plot3', 'mat2cell');

fig = findobj('Tag','mercury_plot');
if feval('isempty',fig)
    fig = figure;
end
clf(fig);
set(fig, 'Name', 'mercury_plot');
set(fig, 'DoubleBuffer', 'on', 'Tag', 'mercury_plot');
set(fig, 'Renderer', 'painters');
set(fig, 'Color', 'black');
set(fig, 'InvertHardcopy', 'off');
pos = get(fig, 'Position');
% No cell array support in MATLAB Coder, so we use these tricks to avoid
% it.
s1 = feval('struct','type','()','subs',mat2cell(mat2cell(1,1),1));
s2 = feval('struct','type','()','subs',mat2cell(mat2cell(2,1),1));
set(fig, 'Position', feval('horzcat',feval('subsref',pos,s1),feval('subsref',pos,s2),400,400));
hold('on');
fig_axes = get(fig, 'CurrentAxes');
init_axes(fig_axes);

% Sub-function to initialize the coordinate axes of the figure window.
function init_axes(a)
coder.extrinsic('set');
LightMinute = 60;
set(a, 'CameraTarget', [0,0,0] );
set(a, 'CameraPosition', [100*LightMinute,0*LightMinute,100*LightMinute]);
set(a, 'CameraViewAngle', 2 );
set(a, 'CameraUpVector', [0,1,0]);
set(a, 'Visible', 'off' ); 
set(a, 'XLim', [-100*LightMinute,100*LightMinute]);
set(a, 'YLim', [-100*LightMinute,100*LightMinute]);
set(a, 'ZLim', [-100*LightMinute,100*LightMinute]);
set(a, 'PlotBoxAspectRatio', [2 1 1]);
