function [day,month,year]=dayjulian(JD)
%   Program for finding the Gregorian date for a given julian day
%   The algorithm is valid for any Gregorian date producing a Julian date
%   greater than zero.
%   (From Fliegel H.F. and van Flandern T.C. 1968. Comm. ACM, 11, p.657)
%
% years B.C. are decreased by 1, as the year 0 does not exist
% in this calculation year 0 is effectivly the same as year -1

L = fix( JD + 68569 );
N = fix( 4*L/146097 );
L = L - fix( (146097*N + 3)/4 );
I = fix( 4000 * (L + 1)/1461001 );
L = L - fix( 1461*I/4 ) + 31 ;
J = fix( 80*L/2447 ) ;
day = L - fix( 2447*J/80 );
L = fix( J/11 );
month = fix( J + 2 - 12*L);
year = fix(100*(N - 49) + I + L);

I=year<=0; 
if sum(I)>0; year(I)=year(I)-1; end
if nargout<=1;
  ii=size(day);
  if ii(2)==1;
     day=[day,month,year];
  else
     day=[day;month;year]';
  end    
end