function [L,B,R]=ple_mercury(Date)
% Low accuracy ephemerides for Mercury
% Package: celestial.SolarSys
% Description: Low accuracy ephemerides for Mercury. Calculate
%              Mercury heliocentric longitude, latitude and radius
%              vector referred to mean ecliptic and equinox of date.
%              Accuarcy: better than 1' in long/lat, ~0.001 au in dist.
% Input  : - matrix of dates, [D M Y frac_day] per line,
%            or JD per line. In TT time scale.
% Output : - Longitude in radians.
%          - Latitude in radians.
%          - Radius vector in au.
% Reference: VSOP87
% See also: ple_planet.m
% Tested : Matlab 5.3
%     By : Eran O. Ofek                    Oct 2001
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: [L,B,R]=celestial.SolarSys.ple_mercury(2451545+(0:1:10)')
% Reliable: 2
%------------------------------------------------------------------------------


RAD = 180./pi;
FunTPI = @(X) (X./(2.*pi) - floor(X./(2.*pi))).*2.*pi;

SizeDate = size(Date);
N        = SizeDate(1);
ColN     = SizeDate(2);

if (ColN==4),
   JD = celestial.time.julday(Date);
elseif (ColN==1),
   JD = Date;
else
   error('Illigal number of columns in date matrix');
end

Tau = (JD(:) - 2451545.0)./365250.0;

SumL0 = 440250710 ...
   + 40989415.*cos(1.48302034 + 26087.90314157.*Tau) ...
   + 5046294.*cos(4.4778549 + 52175.8062831.*Tau) ...
   + 855347.*cos(1.165203 + 78263.709425.*Tau) ...
   + 165590.*cos(4.119692 + 104351.612566.*Tau) ...
   + 34562.*cos(0.77931 + 130439.51571.*Tau) ...
   + 7583.*cos(3.7135 + 156527.4188.*Tau) ...
   + 3560.*cos(1.5120 + 1109.3786.*Tau) ...
   + 1803.*cos(4.1033 + 5661.3320.*Tau) ...
   + 1726.*cos(0.3583 + 182615.3220.*Tau) ...
   + 1590.*cos(2.9951 + 25028.5212.*Tau) ...
   + 1365.*cos(4.5992 + 27197.2817.*Tau) ...
   + 1017.*cos(0.8803 + 31749.2352.*Tau);

SumL1 = 2608814706223 ...
   + 1126008.*cos(6.2170397 + 26087.9031416.*Tau) ...
   + 303471.*cos(3.055655 + 52175.806283.*Tau) ...
   + 80538.*cos(6.10455 + 78263.70942.*Tau) ...
   + 21245.*cos(2.83532 + 104351.61257.*Tau) ...
   + 5592.*cos(5.8268 + 130439.5157.*Tau) ...
   + 1472.*cos(2.5185 + 156527.4188.*Tau);

SumL2 = 53050 ...
   + 16904.*cos(4.69072 + 26087.90314.*Tau) ...
   + 7397.*cos(1.3474 +52175.8063.*Tau) ...
   + 3018.*cos(4.4564 + 78263.7094.*Tau) ...
   + 1107.*cos(1.2623 + 104351.6126.*Tau);

SumL3 = 188.*cos(0.035 + 52175.806.*Tau);

SumL4 = -114;

SumL5 = -1;


L = SumL0 + SumL1.*Tau + SumL2.*Tau.^2 ...
          + SumL3.*Tau.^3 + SumL4.*Tau.^4 + SumL5.*Tau.^5;
L = L.*1e-8;

L = FunTPI(L);

SumB0 = 11737529.*cos(1.98357499 + 26087.90314157.*Tau) ...
   + 2388077.*cos(5.0373896 + 52175.8062831.*Tau) ...
   - 1222840 ...
   + 543252.*cos(1.796444 + 78263.709425.*Tau) ...
   + 129779.*cos(4.832325 + 104351.612566.*Tau) ...
   + 31867.*cos(1.58088 + 130439.51571.*Tau) ...
   + 7963.*cos(4.6097 + 156527.4188.*Tau) ...
   + 2014.*cos(1.3532 + 182615.3220.*Tau);

SumB1 = 429151.*cos(3.501698 + 26087.903142.*Tau) ...
   - 146234 ...
   + 22675.*cos(0.01515 + 52175.80628.*Tau) ...
   + 10895.*cos(0.48540 + 78263.70942.*Tau) ...
   + 6353.*cos(3.4294 + 104351.6126.*Tau) ...
   + 2496.*cos(0.1605 + 130439.5157.*Tau);

SumB2 = 11831.*cos(4.79066 + 26087.90314.*Tau) ...
   + 1914 ...
   + 1045.*cos(1.2122 + 52175.8063.*Tau);

SumB3 = 235.*cos(0.354 + 26087.903.*Tau) + 161;

SumB4 = 4.*cos(1.75 + 26087.90.*Tau) - 1;

B = SumB0 + SumB1.*Tau + SumB2.*Tau.^2 ...
          + SumB3.*Tau.^3 + SumB4.*Tau.^4;
B = B.*1e-8;

%B = FunTPI(B);

SumR0 = 39528272 ...
   + 7834132.*cos(6.1923372 + 26087.9031416.*Tau) ...
   + 795526.*cos(2.959897 + 52175.806283.*Tau) ...
   + 121282.*cos(6.010642 + 78263.709425.*Tau) ...
   + 21922.*cos(2.77820 + 104351.61257.*Tau);

SumR1 = 217348.*cos(4.656172 + 26087.903142.*Tau) ...
   + 44142.*cos(1.42386 + 52175.806283.*Tau) ...
   + 10094.*cos(4.47466 + 78263.70942.*Tau);

SumR2 = 3118.*cos(3.0823 + 26087.9031.*Tau);

SumR3 = 33.*cos(1.68 + 26087.90.*Tau);

R = SumR0 + SumR1.*Tau + SumR2.*Tau.^2 ...
          + SumR3.*Tau.^3;
R = R.*1e-8;

