function [L,B,R]=ple_neptune(Date)
% Low accuracy planetray ephemeris for Neptune
% Package: celestial.SolarSys
% Description: Low accuracy planetray ephemeris for Neptune. Calculate
%              Neptune 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_neptune(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 = 531188633 ...
   + 1798476.*cos(2.9010127 + 38.1330356.*Tau) ...
   + 1019728.*cos(0.4858092 + 1.4844727.*Tau) ...
   + 124532.*cos(4.830081 + 36.648563.*Tau) ...
   + 42064.*cos(5.41055 + 2.96895.*Tau) ...
   + 37715.*cos(6.09222 + 35.16409.*Tau) ...
   + 33785.*cos(1.24489 + 76.26607.*Tau) ...
   + 16483.*cos(0.00008 + 491.55793.*Tau) ...
   + 9199.*cos(4.9375 + 39.6175.*Tau) ...
   + 8994.*cos(0.2746 + 175.1661.*Tau) ...
   + 4216.*cos(1.9871 + 73.2971.*Tau) ...
   + 3365.*cos(1.0359 + 33.6796.*Tau) ...
   + 2285.*cos(4.2061 + 4.4534.*Tau) ...
   + 1434.*cos(2.7834 + 74.7816.*Tau);

SumL1 = 3837687717 ...
   + 16604.*cos(4.86319 + 1.48447.*Tau) ...
   + 15807.*cos(2.27923 + 38.13304.*Tau) ...
   + 3335.*cos(3.6820 + 76.2661.*Tau) ...
   + 1306.*cos(3.6732 + 2.9689.*Tau) ...
   + 605.*cos(1.505 + 35.164.*Tau);

SumL2 = 53893 ...
   + 296.*cos(1.855 + 1.484.*Tau) ...
   + 281.*cos(1.191 + 38.133.*Tau) ...
   + 270.*cos(5.721 + 76.266.*Tau);

SumL3 = 31 ...
   + 15.*cos(1.35 + 76.27.*Tau);

SumL4 = -114;


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

L = FunTPI(L);

SumB0 = 3088623.*cos(1.4410437 + 38.1330356.*Tau) ...
   + 27780.*cos(5.91272 + 76.26607.*Tau) ...
   + 27624 ...
   + 15448.*cos(3.50877 + 39.61751.*Tau) ...
   + 15355.*cos(2.52124 + 36.64856.*Tau) ...
   + 2000.*cos(1.5100 + 74.7816.*Tau) ...
   + 1968.*cos(4.3778 + 1.4845.*Tau) ...
   + 1015.*cos(3.2156 + 35.1641.*Tau) ...
   + 606.*cos(2.802 + 73.297.*Tau);

SumB1 = 227279.*cos(3.807931 + 38.133036.*Tau) ...
   + 1803.*cos(1.9758 + 76.2661.*Tau) ...
   - 1433 ...
   + 1386.*cos(4.8256 + 36.6486.*Tau) ...
   + 1073.*cos(6.0805 + 39.6175.*Tau);

SumB2 = 9691.*cos(5.5712 + 38.1330.*Tau) ...
   + 79.*cos(3.63 + 76.27.*Tau) ...
   + 72.*cos(0.45 + 36.65.*Tau);

SumB3 = 273.*cos(1.017 + 38.133.*Tau);

SumB4 = 6.*cos(2.67 + 38.13.*Tau);

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

%B = FunTPI(B);

SumR0 = 3007013206 ...
   + 27062259.*cos(1.32999459 + 38.13303564.*Tau) ...
   + 1691764.*cos(3.2518614 + 36.6485629.*Tau) ...
   + 807831.*cos(5.185928 + 1.484473.*Tau) ...
   + 537761.*cos(4.521139 + 35.164060.*Tau) ...
   + 495726.*cos(1.571057 + 491.557929.*Tau) ...
   + 274572.*cos(1.845523 + 175.166060.*Tau) ...
   + 135134.*cos(3.372206 + 39.617508.*Tau) ...
   + 121802.*cos(5.797544 + 76.266071.*Tau) ...
   + 100895.*cos(0.377027 + 73.297126.*Tau) ...
   + 69792.*cos(3.79617 + 2.96895.*Tau) ...
   + 46688.*cos(5.74938 + 33.67962.*Tau) ...
   + 24594.*cos(0.50802 + 109.94569.*Tau) ...
   + 16939.*cos(1.59422 + 71.81265.*Tau) ...
   + 14230.*cos(1.07786 + 74.78160.*Tau) ...
   + 12012.*cos(1.92062 + 1021.24889.*Tau) ...
   + 8395.*cos(0.6782 + 146.5943.*Tau) ...
   + 7572.*cos(1.0715 + 388.4652.*Tau) ...
   + 5721.*cos(2.5906 + 4.4534.*Tau);

SumR1 = 236339.*cos(0.704980 + 38.133036.*Tau) ...
   + 13220.*cos(3.32015 + 1.48447.*Tau) ...
   + 8622.*cos(6.2163 + 35.1641.*Tau) ...
   + 2702.*cos(1.8814 + 39.6175.*Tau) ...
   + 2155.*cos(2.0943 + 2.9689.*Tau) ...
   + 2153.*cos(5.1687 + 76.2661.*Tau) ...
   + 1603 ...
   + 1464.*cos(1.1842 + 33.6796.*Tau) ...
   + 1136.*cos(3.9189 + 36.6486.*Tau);

SumR2 = 4247.*cos(5.8991 + 38.1330.*Tau) ...
   + 218.*cos(0.346 + 1.484.*Tau);

SumR3 = 166.*cos(4.552 + 38.133.*Tau);

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

