function [gtime]=gregorian(jd)
% GREGORIAN: Converts Julian day numbers to Gregorian calendar.
%
% USAGE: [gtime]=gregorian(jd)
%
% DESCRIPTION: Converts decimal Julian days to Gregorian dates using the
% astronomical convension, but with time zero starting
% at midnight instead of noon. In this convention,
% Julian day 2440000 begins at 0000 hours, May 23, 1968.
% The Julian day does not have to be an integer, and with
% Matlab's double precision, the accuracy of decimal days
% is about 0.1 milliseconds.
%
%
% INPUT: jd = decimal Julian days
%
% OUTPUT: gtime = six column Gregorian time matrix, where each row is
% [yyyy mo da hr mi sec].
% yyyy = year (e.g., 1979)
% mo = month (1-12)
% da = day (1-31)
% hr = hour (0-23)
% mi = minute (0-59)
% sec = decimal seconds
% example: [1990 12 12 0 0 0] is midnight on Dec 12, 1990.
%
% AUTHOR: Rich Signell (rsignell@usgs.gov)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is self described.
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty,
% expressed or implied, is made by the USGS or the United States
% Government as to the accuracy and functioning of the program
% and related program material nor shall the fact of distribution
% constitute any such warranty, and no responsibility is assumed
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Add 0.2 milliseconds before Gregorian calculation to prevent
% roundoff error resulting from math operations on time
% from occasionally representing midnight as
% (for example) [1990 11 30 23 59 59.99...] instead of [1990 12 1 0 0 0]);
% If adding a 0.2 ms to time (each time you go back and forth between
% Julian and Gregorian) bothers you more than the inconvenient representation
% of Gregorian time at midnight you can comment this line out...
jd=jd+2.e-9;
% if you want Julian Days to start at noon...
% h=rem(jd,1)*24+12;
% i=(h >= 24);
% jd(i)=jd(i)+1;
% h(i)=h(i)-24;
secs=rem(jd,1)*24*3600;
j = floor(jd) - 1721119;
in = 4*j -1;
y = floor(in/146097);
j = in - 146097*y;
in = floor(j/4);
in = 4*in +3;
j = floor(in/1461);
d = floor(((in - 1461*j) +4)/4);
in = 5*d -3;
m = floor(in/153);
d = floor(((in - 153*m) +5)/5);
y = y*100 +j;
mo=m-9;
yr=y+1;
i=(m<10);
mo(i)=m(i)+3;
yr(i)=y(i);
[hour,min,sec]=s2hms(secs);
gtime=[yr(:) mo(:) d(:) hour(:) min(:) sec(:)];