function MP=moon_phases(StartDate,EndDate)
% Return a list of moon phases in range of dates
% Package: celestial.SolarSys
% Description: Return a list of moon phases in range of dates.
% Input  : - Minimum date [D M Y] or JD.
%          - Maximum date [D M Y] or JD.
% Output : - List of Moon phases [Phase, JD].
%            Where Phase are: 0,1,2,3 for new moon, first quarter,
%            full moon, and last quarter, respectively.
%            JD is given in the TT time scale.
% Reference: Meeus 1991
% Tested : Matlab 7.3
%     By : Eran O. Ofek                    Jul 2007
%    URL : http://weizmann.ac.il/home/eofek/matlab/
% Example: MP=celestial.SolarSys.moon_phases([1 1 2015],[1 1 2016]);
% Reliable: 2
%--------------------------------------------------------------------------

RAD = 180./pi;

if (length(StartDate)==1),
   StartJD = StartDate;
else
   StartJD = celestial.time.julday(StartDate);
end
if (length(EndDate)==1),
   EndJD = EndDate;
else
   EndJD = celestial.time.julday(EndDate);
end

FTP = @(X) 2.*pi.*(X./(2.*pi)-floor(X./(2.*pi)));

%K   = (Year - 2000).* 12.3685;
StartK   = (StartJD - 2451545)./365.242198.*12.3685;
EndK     = (EndJD   - 2451545)./365.242198.*12.3685;

OrigK    = [floor(StartK - 3):1:ceil(EndK+3)].';
NK       = length(OrigK);

% Additional for all phases
C  = [0.000325, 0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060, 0.000056, 0.000047, 0.000042, 0.000040, 0.000037, 0.000035, 0.000023];

MP = zeros(0,2);
% for each Moon Phase
for PhaseN=0:1:3,

   Phase = PhaseN./4;
   K = OrigK + Phase;

   T   = K./1236.85;

   JDE = 2451550.09765 + 29.530588853.*K + ...
                          0.0001337  .*T.^2 - ...
                          0.000000150.*T.^3 + ...
                          0.00000000073.*T.^4;

   M  =   2.5534 + 29.10535669.*K  - 0.0000218.*T.^2 - 0.00000011.*T.^3;
   Mt = 201.5643 + 385.81693528.*K + 0.0107438.*T.^2 + 0.00001239.*T.^3 - 0.000000058.*T.^4;
   F  = 160.7108 + 390.67050274.*K - 0.0016341.*T.^2 - 0.00000227.*T.^3 + 0.000000011.*T.^4;
   Om = 124.7746 - 1.56375580.*K   + 0.0020691.*T.^2 + 0.00000216.*T.^3;
   M  = M./RAD;
   Mt = Mt./RAD;
   F  = F./RAD;
   Om = Om./RAD;
   M  = FTP(M);
   Mt = FTP(Mt);
   F  = FTP(F);
   Om = FTP(Om);
   E  = 1 - 0.002516.*T - 0.0000074.*T.^2;

   A(:,1) = (299.77 + 0.107408.*K - 0.009173.*T.^2)./RAD;
   A(:,2) = (251.88 + 0.016321.*K)./RAD;
   A(:,3) = (251.83 + 26.651886.*K)./RAD;
   A(:,4) = (349.42 + 36.412478.*K)./RAD;
   A(:,5) = ( 84.66 + 18.206239.*K)./RAD;
   A(:,6) = (141.74 + 53.303771.*K)./RAD;
   A(:,7) = (207.14 + 2.453732.*K)./RAD;
   A(:,8) = (154.84 + 7.306860.*K)./RAD;
   A(:,9) = ( 34.52 + 27.261239.*K)./RAD;
   A(:,10)= (207.19 + 0.121824.*K)./RAD;
   A(:,11)= (291.34 + 1.844379.*K)./RAD;
   A(:,12)= (161.72 + 24.198154.*K)./RAD;
   A(:,13)= (239.56 + 25.513099.*K)./RAD;
   A(:,14)= (331.55 + 3.592518.*K)./RAD;

   PAC = zeros(NK,1);
   for IK=1:1:NK,
      PAC(IK) = sum(C.*sin(A(IK,:)));
   end

   switch PhaseN
    case 0
       % New Moon
       PC = -0.40720.*   sin(Mt)...
            +0.17241.*E.*sin(M)...
            +0.01608.*   sin(2.*Mt)...
            +0.01039.*   sin(2.*F)...
            +0.00739.*E.*sin(Mt-M)...
            -0.00514.*E.*sin(Mt+M)...
            +0.00208.*E.*E.*sin(2.*M)...
            -0.00111.*sin(Mt-2.*F)...
            -0.00057.*sin(Mt+2.*F)...
            +0.00056.*E.*sin(2.*Mt+M)...
            -0.00042.*sin(3.*Mt)...
            +0.00042.*E.*sin(M+2.*F)...
            +0.00038.*E.*sin(M-2.*F)...
            -0.00024.*E.*sin(2.*Mt-M)...
            -0.00017.*sin(Om)...
            -0.00007.*sin(Mt+2.*M)...
            +0.00004.*sin(2.*(Mt-F))...
            +0.00004.*sin(3.*M)...
            +0.00003.*sin(Mt+M-2.*F)...
            +0.00003.*sin(2.*(Mt+F))...
            -0.00003.*sin(Mt+M+2.*F)...
            +0.00003.*sin(Mt-M+2.*F)...
            -0.00002.*sin(Mt-M-2.*F)...
            -0.00002.*sin(3.*Mt+M)...
            +0.00002.*sin(4.*Mt);

       W = 0;
    case 2
       % Full Moon
       PC = -0.40614.*sin(Mt)...
            +0.17302.*E.*sin(M)...
            +0.01614.*sin(2.*Mt)...
            +0.01043.*sin(2.*F)...
            +0.00734.*E.*sin(Mt-M)...
            -0.00515.*E.*sin(Mt+M)...
            +0.00209.*E.*E.*sin(2.*M)...
            -0.00111.*sin(Mt-2.*F)...
            -0.00057.*sin(Mt+2.*F)...
            +0.00056.*E.*sin(2.*Mt+M)...
            -0.00042.*sin(3.*Mt)...
            +0.00042.*E.*sin(M+2.*F)...
            +0.00038.*E.*sin(M-2.*F)...
            -0.00024.*E.*sin(2.*Mt-M)...
            -0.00017.*sin(Om)...
            -0.00007.*sin(Mt+2.*M)...
            +0.00004.*sin(2.*(Mt-F))...
            +0.00004.*sin(3.*M)...
            +0.00003.*sin(Mt+M-2.*F)...
            +0.00003.*sin(2.*(Mt+F))...
            -0.00003.*sin(Mt+M+2.*F)...
            +0.00003.*sin(Mt-M+2.*F)...
            -0.00002.*sin(Mt-M-2.*F)...
            -0.00002.*sin(3.*Mt+M)...
            +0.00002.*sin(4.*Mt);

       W = 0;
   
    case {1,3}
       % First and last Quarters
       PC = -0.62801.*sin(Mt)...
            +0.17172.*E.*sin(M)...
            -0.01183.*E.*sin(Mt+M)...
            +0.00862.*sin(2.*Mt)...
            +0.00804.*sin(2.*F)...
            +0.00454.*E.*sin(Mt-M)...
            +0.00204.*E.*E.*sin(2.*M)...
            -0.00180.*sin(Mt-2.*F)...
            -0.00070.*sin(Mt+2.*F)...
            -0.00040.*sin(3.*Mt)...
            -0.00034.*E.*sin(2.*Mt-M)...
            +0.00032.*E.*sin(M+2.*F)...
            +0.00032.*E.*sin(M-2.*F)...
            -0.00028.*E.*E.*sin(Mt+2.*M)...
            +0.00027.*E.*sin(2.*Mt+M)...
            -0.00017.*sin(Om)...
            -0.00005.*sin(Mt-M-2.*F)...
            +0.00004.*sin(2.*(Mt+F))...
            -0.00004.*sin(Mt+M+2.*F)...
            +0.00004.*sin(Mt-2.*M)...
            +0.00003.*sin(Mt+M-2.*F)...
            +0.00003.*sin(3.*M)...
            +0.00002.*sin(2.*(Mt-F))...
            +0.00002.*sin(Mt-M+2.*F)...
            -0.00002.*sin(3.*Mt+M);

       % For the quarter phases only:
       W = 0.00306 - 0.00038.*E.*cos(M) + 0.00026.*cos(Mt) - 0.00002.*cos(Mt-M) + 0.00002.*cos(Mt+M) + 0.00002.*cos(2.*F);
       switch PhaseN
        case 1
           W = +W;
        case 3
           W = -W;
        otherwise
           error('Unknown PhaseN Option');
       end
    otherwise
       error('Unknown PhaseN Option');
   end

   MP = [MP; [PhaseN.*ones(NK,1), JDE + PC + W + PAC]];

end

FlagI  = MP(:,2)>=StartJD & MP(:,2)<=EndJD;
MP = MP(FlagI,:);
MP = sortrows(MP,2);
