applan1 (tjd, ujd, l, n, topo, glon, glat, ht)

function [ra, dec, dis] = applan1 (tjd, ujd, l, n, topo, glon, glat, ht)

% this function computes the apparent geocentric or topocentric place
% of a planet or other solar system body. rectangular coordinates of
% solar system bodies are obtained from function solsys.

% jpl binary ephemeris

% input

%  tjd  = tdt julian date for apparent geocentric place
%  ujd  = ut1 julian date for apparent topocentric place
%  l    = body identification number for desired planet
%  n    = body identification number for the earth
%  topo = type of apparent place calculation
%       = 0 ==> geocentric
%       = 1 ==> topocentric
%  glon = geodetic longitude of observer (east +, degrees)
%  glat = geodetic latitude of observer (north +, degrees)
%  ht   = height of observer (meters)

% output

%  ra  = apparent geocentric or topocentric right ascension,
%        referred to true equator and equinox of date (hours)
%  dec = apparent geocentric or topocentric declination,
%        referred to true equator and equinox of date (degrees)
%  dis = true distance from earth to planet (astronomical units)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global tlasta peb veb pes ves

c = 173.1446333;

t0 = 2451545;

% compute t1, the tdb julian date corresponding to tjd

[x, secdif] = tdtimes (tjd);

t1 = tjd + secdif / 86400;

if (abs(tjd - tlasta) < 1e-6)
    
    for j = 1:1:3
        
        pb(j) = peb(j);
        vb(j) = veb(j);
        ps(j) = pes(j);
        vs(j) = ves(j);
        
    end
    
else
    
    % get position and velocity of the earth wrt barycenter
    % of solar system and wrt center of sun

    [peb, veb, ierr] = solsys (t1, n, 0);

    if (ierr ~= 0)
        ra = 0;
        dec = 0;
        dis = 0;
        tlasta = 0;
        return;
    end

    [pes, ves, ierr] = solsys (t1, n, 1);

    if (ierr ~= 0)
        ra = 0;
        dec = 0;
        dis = 0;
        tlasta = 0;
        return;
    end

    tlasta = 0;

    for j = 1:1:3
        pb(j) = peb(j);
        vb(j) = veb(j);
        ps(j) = pes(j);
        vs(j) = ves(j);
    end
end

lplan = l;

if (topo == 1)
    
    % get position and velocity of observer wrt center of earth

    st = gast2 (ujd, 0, 0);

    [x1, x2, eqeq, x3, x4] = etilt1 (t1);

    gast = st + eqeq / 3600;

    [pos1, vel1] = terra (glon, glat, ht, gast);

    pos2 = nutate1 (-t1, pos1);

    pog = precess (t1, pos2, t0);

    vel2 = nutate1 (-t1, vel1);

    vog = precess (t1, vel2, t0);

    % compute position and velocity of observer wrt barycenter of
    % solar system and wrt center of sun

    for j = 1:1:3
        pb(j) = peb(j) + pog(j);
        vb(j) = veb(j) + vog(j);
        ps(j) = pes(j) + pog(j);
        vs(j) = ves(j) + vog(j);
    end
    
end

% get position of planet wrt barycenter of solar system

[pos1, vel1, ierr] = solsys (t1, lplan, 0);

if (ierr ~= 0)
    
    ra = 0;
    dec = 0;
    dis = 0;
    tlasta = 0;
    return;
    
end

[pos2, tlight] = geocen (pos1, pb);

s = tlight * c;

t2 = t1 - tlight;

while(1)
    
    [pos1, vel1, ierr] = solsys (t2, lplan, 0);

    if (ierr ~= 0)
        ra = 0;
        dec = 0;
        dis = 0;
        tlasta = 0;
        return;
    end

    [pos2, tlight] = geocen (pos1, pb);

    t3 = t1 - tlight;

    if (abs(t3 - t2) < 1e-8)
        
        pos3 = sunfld (pos2, ps);

        pos4 = aberat (pos3, vb, tlight);

        pos5 = precess (t0, pos4, t1);

        pos6 = nutate1 (t1, pos5);

        [r, d] = angles (pos6);

        ra = r;

        dec = d;

        dis = s;

        break;
    end

    t2 = t3;

end