function jd= julian(mon,day,year,hhmm,method70)
%	Calculate Julian day.
%
%	function jd= julian(mon,day,year[,hhmm,[method70]])
%	  hhmm= time {optional}
%	    if the hours are given, then 1,1,85,0000 returns 0, while
%	    if no hours are given, then  1,1,85      returns 1.  Get it?
%	  method70= calculate jd as # days since 1/1/70. (y or n) 
%	            default= n {optional}
%

%	Latest revision by P. Dragos 7-Jun-91.
%
%	P. Dragos, Battelle Ocean Sciences
%	Duxbury, MA 02332
%	(617) 934-0571
%
%	Copyright (c) 1991 by Paul Dragos.
%	This program is copyrighted.  However, you are encouraged to
%	copy and share it, so long as no charge is imposed.

%				---------------------------------------
%				data
yr= [70 : 1 : 99];
dyr= [	 0  365  730 1096 1461 1826 2191 2557 2922 3287 ...
     3652 4018 4383 4748 5113 5479 5844 6209 6574 6940 ...
     7305 7670 8035 8401 8766 9131 9496 9862 10227 10592];
d= [0 31 28 31 30 31 30 31 31 30 31 30 31];
dmon= cumsum(d)';
m30=[4 6 9 11];
%				---------------------------------------
%				condition the input if necessary
sz= size(mon);
if (sz(2) > sz(1)), mon=mon'; day=day'; year=year'; sz=size(mon); end;
if ((length(mon) ~= length(day)) ...
  | (length(mon) ~= length(year)) ...
  | (length(mon) ~= length(hhmm))),
  error(' Vectors not the same length.'); end;
if year(1) > 99, year= year-1900; end;
%				---------------------------------------
%				check for impossible dates
if any(year > 99 | year < 56), error(' Years must be 1970 -1999!'); end;
er=0;
if any(mon > 12 | mon < 1), er=1; end;
if any(day > 31 | day < 1), er=1; end;
j=find(mon == 4 | mon == 6 | mon == 9 | mon == 11);
if any(day(j) > 30), er=1; end;
j= find (year./4 ~= fix(year./4) & mon == 2);
if any(day(j) > 28), er=1; end;
j= find (year./4 == fix(year./4) & mon == 2);
if any(day(j) > 29), er=1; end;
if (er == 1), error(' Impossible Date!'); end;
%				---------------------------------------
%				defaults
if nargin == 4
  if isstr(hhmm)
    method70= hhmm;
    hhmm= ones(sz(1),sz(2)).*2400;
  else
    method70= 'n';
  end;
elseif nargin == 3
   hhmm= ones(sz(1),sz(2)).*2400;
  method70= 'n';
end;  
hour= fix(hhmm/100);
min= hhmm-hour*100;
hr= -1+(hour+min/60)/24;
%				---------------------------------------
%				calculate day of year, leap day, etc
leap= zeros(sz(1),sz(2));
i= find ((year./4 == fix(year./4)) & (mon > 2));
if (length(i) > 0); leap(i)= 1; end;

dyr70= zeros(sz(1),sz(2));
if (method70 == 'Y'|method70 =='y'),   dyr70= dyr(year-yr(1)+1)'; end;

jd= dmon(mon) + day + leap + hr + dyr70;

return;
