function Eq=equinox_solstice(Year)
% Approximate time of Equinox and Solstice
% Package: celestial.SolarSys
% Description: Calculate the approximate time of Equinox and Solstice
%              for a given list of years. Accurate to about 100s between
%              year -1000 to 3000.
% Input  : - Vector of years.
% Output : - Structure array with 4 fields, each one contains the JD
%            of the March Equinox, June Solstice, September Equinox
%            and December Solstice, respectively.
% Reference: Meeus 1991
% Tested : Matlab R2014a
%     By : Eran O. Ofek                    Jul 2014
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: Eq=celestial.SolarSys.equinox_solstice([1962;2000])
% Reliable: 2
%--------------------------------------------------------------------------
INV_RAD = pi./180;


F1 = Year>-1000 & Year<1000;
F2 = Year>=1000 & Year<3000;

Y   = floor(Year-2000)./1000;
D(1).JD = zeros(size(Year));
D(2).JD = zeros(size(Year));
D(3).JD = zeros(size(Year));
D(4).JD = zeros(size(Year));

%--- for -1000 to 1000 ---
% March equinox
D(1).JD(F1) = 1721139.29189 + 365242.13740.*Y(F1) + 0.06134.*Y(F1).^2 + ...
                     0.00111.*Y(F1).^3 - 0.00071.*Y(F1).^4;
                 
% June solstice
D(2).JD(F1) = 1721233.25401 + 365241.72562.*Y(F1) - 0.05323.*Y(F1).^2 + ...
                     0.00907.*Y(F1).^3 + 0.00025.*Y(F1).^4;

% September equinox
D(3).JD(F1) = 1721325.70455 + 365242.49558.*Y(F1) - 0.11677.*Y(F1).^2 - ...
                     0.00297.*Y(F1).^3 + 0.00074.*Y(F1).^4;
                 
% December solstice
D(4).JD(F1) = 1721414.39987 + 365242.88257.*Y(F1) - 0.00769.*Y(F1).^2 - ...
                       0.00933.*Y(F1).^3  - 0.00006.*Y(F1).^4;

%--- for 1000 to 3000 ---
% March equinox
D(1).JD(F2) = 2451623.80984 + 365242.37404.*Y(F2) + 0.05169.*Y(F2).^2 - ...
                     0.00411.*Y(F2).^3 - 0.00057.*Y(F2).^4;

% June solstice
D(2).JD(F2) = 2451716.56767 + 365241.62603.*Y(F2) + 0.00325.*Y(F2).^2 + ...
    0.00888.*Y(F2).^3 - 0.00030.*Y(F2).^4;
% September equinox
D(3).JD(F2) = 2451810.21715 + 365242.01767.*Y(F2) - 0.11575.*Y(F2).^2 + ...
    0.00337.*Y(F2).^3 + 0.00078.*Y(F2).^4;
% December solstice
D(4).JD(F2) = 2451900.05952 + 365242.74049.*Y(F2) - 0.06223.*Y(F2).^2 - ...
    0.00823.*Y(F2).^3 + 0.00032.*Y(F2).^4;

for I=1:1:4,
    D(I).T = (D(I).JD - 2451545.0)./36525;
    D(I).W = (35999.373.*D(I).T - 2.47).*INV_RAD;
    D(I).DLam = 1 + 0.0334.*cos(D(I).W) + 0.0007.*cos(2.*D(I).W);
                   
    D(I).S =485.*cosd(324.96 + 1934.136.*D(I).T) + ...
            203.*cosd(337.23 + 32964.467.*D(I).T) + ...
            199.*cosd(342.08 + 20.186.*D(I).T) + ...
            182.*cosd(27.85 + 445267.112.*D(I).T) + ...
            156.*cosd(73.14 + 45036.886.*D(I).T) + ...
            136.*cosd(171.52 + 22518.443.*D(I).T) + ...
             77.*cosd(222.54 + 65928.934.*D(I).T) + ...
             74.*cosd(296.72 + 3034.906.*D(I).T) + ...
             70.*cosd(243.58 + 9037.513.*D(I).T) + ...
             58.*cosd(119.81 + 33718.147.*D(I).T) + ...
             52.*cosd(297.17 + 150.678.*D(I).T) + ...
             50.*cosd(21.02 + 2281.226.*D(I).T) + ...
             45.*cosd(247.54 + 29929.562.*D(I).T) + ...
             44.*cosd(325.15 + 31555.956.*D(I).T) + ...
             29.*cosd(60.93 + 4443.417.*D(I).T) + ...
             18.*cosd(155.12 + 67555.328.*D(I).T) + ...
             17.*cosd(288.79 + 4562.452.*D(I).T) + ...
             16.*cosd(198.04 + 62894.029.*D(I).T) + ...
             14.*cosd(199.76 + 31436.921.*D(I).T) + ...
             12.*cosd(95.39 + 14577.848.*D(I).T) + ...
             12.*cosd(287.11 + 31931.756.*D(I).T) + ...
             12.*cosd(320.81 + 34777.259.*D(I).T) + ...
              9.*cosd(227.73 + 1222.114.*D(I).T) + ...
              8.*cosd(15.45 + 16859.074.*D(I).T);


    Eq(I).JD = D(I).JD + 0.00001.*D(I).S./D(I).DLam;
end

