function [ZEN,ALT,AZ] = solarzen_(DateTime,LAT,LON)

% SOLARZEN_ - Calculate Solar Zenith Angle
%
% Calculated the solar zenith angle for a given day and
% position. ALMANAC_.M is used to compute declination, GHA 
% distance and equation of time. These vaules are then input
% into ALTAZM_.M which calculates celestial altitude and azimuth 
% from declination, Greenwich Hour Angle, latitude and longitude. 
%
% Use As: ZEN = solarzen_(DateTime,LAT,LONG)
% Inputs: DateTime = GMT Date and Time (YYMMDD.HHMMSS)
%         LAT      = Latitude (DD.MMmm) where N+ and S-
%         LON      = Longitude (DD.MMmm) where W- and E+
% Output: ZEN      = Solar Zenith angle (degrees)
%         ALT      = Solar altitude     (degrees)
%         AZ       = Solar azimuth      (degrees)
% Notes:
% ALMANAC_.M
% Ref:  Van Flanderen and Pulkkinen (1979) Low Precision Formulae for Planetary 
%       Positions. Astrophysical Journal Supplement Series, 41:391-411. 
% ALTAZM_.M
% Algorithms taken from Wilson (1980). 

% 20 Dec 99; S. Flora - created from nderive_

ZEN = [];

if nargin == 0
  help solarzen_
  return
end

if nargin < 3
  disp('  SOLARZEN_ Error (3 Inputs Required)')
  return
end

LON = -LON;

[Year,Month,Day,Hour,Minute,Second,DOY,GMT] = timeupk_(DateTime);        % Normalize average for station
if Year >= 0 & Year <= 60                                                % Fix up the Year 
  Year = Year + 2000;
else
  Year = Year + 1900;
end

% ALMANAC_ and ALTAZM_ replaces SUNANG2_
[DECL,GHA,DIST,EQ_TIME] = almanac_(Year,Month,Day,GMT);
avg_latd = dms2deg_(LAT,0);                                    % Average and Convert to DD.dd
avg_lond = dms2deg_(LON,0);                                    % Average and Convert to DD.dd
[ALT,AZ] = altazm_(DECL,GHA,avg_latd,avg_lond);                % Calculate azimuth and altitude
ZEN = (90 - ALT);                                              % Zenith angle in degrees

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [DECL, GHA, DIST, EQ_TIME] = almanac_(YEAR,MONTH,DAY,HOUR) 
% ALMANAC_ - Moderate precision solar ephemeris computes declination, GHA 
%            distance and equation of time. 
% 
% Use As:   [DECL,GHA,DIST,EQ_TIME] = almanac_(y,mo,d,h) 
% 
% Input:    y         = year, as 1996 
%           mo        = month of the year (1-12) 
%           d         = GMT Day of the month (1-31)  
%           h         = GMT decimal hour of the day (0-24) 
% Output:   DECL      = Apparent Declination      (deg) 
%           GHA       = Greenwich Hour Angle      (deg) 
%           DIST      = Earth-Sun Distance        (AU) 
%           EQ_TIME   = Equation of Time          (deg) 
% 
% Example:  [a,b,c,d] = almanac_(1996,6,12,16.1) 
%                   a =  23.19; b = 61.52; c = 1.0156; d = -0.0287 
% 
% Note:     This algorithm is accurate to about 30 arc seconds.
 
%           The following internal values may be of interest  
%           CENT_1900 = Fractional century from 1.5 Jan 1900 
%           DAYS_2000 = Days from               1.5 Jan 2000   
 
%           GMST      = Greenwich Mean Siderial Time (hours) 
%           MEAN_RA   = Mean Right Ascension         (deg)                             
%           RA        = Apparent Right Ascension     (deg) 
%           SUN_LON   = Apparent Sun Longitude       (deg) 
%           SHA       = Siderial Hour Angle          (deg) 
%           GHA_ARIES = Greenwich Hour Angle Aries   (deg) 
% 
% See Also: SUNANG2_, DEG2DMS_, EPHEMS_, ALTAZM_, SUNALT_, SUNRISE_, SUNSET_, NOON_ 
% Ref:  Van Flanderen and Pulkkinen (1979) Low Precision Formulae for Planetary 
%       Positions. Astrophysical Journal Supplement Series, 41:391-411. 
 
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% Richard Reaves & William Broenkow 
% 21 March 1988 
% 28 May 96; W. Broenkow converted from NEWSUN.FOR           
% Adapted from PLOTephem.HPL and 
% 20 Jun 96; wwb changed degrad to deg2rad 
%  Internal parameter names (P10-P20) are identical to those in the  
%  HPL program, PLOTephem, from which this was extracted. 
 
% [CENT_1900,DAYS_2000,GMST] = EPHEMPAR(YEAR,MONTH,DAY,HOUR); 
% Do not use the needless function 
% This program agrees with Wilson (1980) and Duffet-Smith (1997) 
% W. Broenkow 21 Jul 1997 
 
JUL_HR    = juldate_(YEAR,MONTH,DAY) + HOUR/24.0;  % Julian date + hours 
DAYS_2000 = JUL_HR - 2451545.0; 
CENT_1900 = 1.0 + DAYS_2000/36525.0; 
GMST      = 6.6460656 + 2400.051262*CENT_1900 + 0.0000258*CENT_1900.^2; 
GMST      = rem(GMST,24); 
                         
P10 = 2*pi*(0.827362 + 0.03386319198*DAYS_2000); % originally FRAC360*() 
P11 = 2*pi*(0.347343 - 0.00014709391*DAYS_2000); % just convert to 2*pi radians 
P12 = 2*pi*(0.779072 + 0.00273790931*DAYS_2000); % and make angle positive 
P13 = 2*pi*(0.993126 + 0.00273777850*DAYS_2000); % Even without making angles 
P14 = 2*pi*(0.140023 + 0.00445036173*DAYS_2000); % the answer should be correct. 
P15 = 2*pi*(0.053856 + 0.00145561327*DAYS_2000); 
P16 = 2*pi*(0.056531 + 0.00023080893*DAYS_2000); 
 
P17 = (6910.0 - 17.0*CENT_1900).*sin(P13);       % all trig arguments are radians 
P17 =  P17  + 72.0.*sin(2.0*P13); 
P17 =  P17  -  7.0.*cos(P13 - P16); 
P17 =  P17  +  6.0.*sin(P10);          
P17 =  P17  +  5.0.*sin(4.0*P13 - 8.0*P15 + 3.0*P16); 
P17 =  P17  -  5.0.*cos(2.0*P13 - 2.0*P14);  
P17 =  P17  -  4.0.*sin(P13 - P14); 
P17 =  P17  +  4.0.*cos(4.0*P13 - 8.0*P15 + 3.0*P16); 
P17 =  P17  +  3.0.*sin(2.0*P13 - 2.0*P14);  
P17 =  P17  -  3.0.*sin(P16);              
P17 =  P17  -  3.0.*sin(2.0*P13 - 2.0*P16); 
                                            
P18 = (0.39785 - 0.00021*CENT_1900).*sin(P12); 
P18 =  P18 + (0.00003.*CENT_1900 - 0.01).*sin(P12-P13); 
P18 =  P18 +  0.00333.*sin(P12 + P13); 
P18 =  P18 +  0.00004.*sin(P12 + 2.0*P13); 
P18 =  P18 -  0.00004.*cos(P12); 
P18 =  P18 -  0.00004.*sin(P11 - P12); 
 
P19 = (-3349.0 + 8.0.*CENT_1900).*cos(P13); 
P19 = P19 - 14.0.*cos(2.0*P13); 
P19 = P19 -  3.0.*sin(P13 - P16); 
P19 = 1.0 + P19./1e5; 
 
P20 = (-0.04129 + 0.00005*CENT_1900).*sin(2*P12); 
P20 = P20 + (0.03211 - 0.00008*CENT_1900).*sin(P13); 
P20 = P20 + 0.00104.*sin(2*P12 - P13); 
P20 = P20 - 0.00035.*sin(2*P12 + P13) - 0.0001; 
P20 = P20 - 0.00008.*sin(P11); 
P20 = P20 + 0.00007.*sin(2*P13); 
P20 = P20 + 0.00003.*sin(P10); 
P20 = P20 - 0.00002.*cos(P13 - P16); 
P20 = P20 + 0.00002.*sin(4.0.*P13 - 8.0.*P15 + 3.0.*P16); 
P20 = P20 - 0.00002.*sin(P13 - P14); 
P20 = P20 - 0.00002.*cos(2.0.*P13 - 2.0.*P14); 
 
%  I went to some trouble to discover all of this... 
%  so keep the commented lines as a tutorial. 
 
REQ_TIME   = asin(P20./sqrt(P19 - P18.*P18)); % Equation of Time (radians) 
EQ_TIME    = rad2deg_(REQ_TIME);              % Equation of Time (degrees) 
DIST       = 1.00021.*sqrt(P19);              % Earth-Sun Distance (A.U.) 
RDECL      = asin(P18./sqrt(P19));            % Declination (radians) 
DECL       = rad2deg_(RDECL);                 % Declination (degrees) 
RMEAN_RA   = P12;                             % Mean Right Ascension (radians) 
MEAN_RA    = rad2deg_(P12);                   % Mean Right Ascension (degrees) 
RRA        = RMEAN_RA + REQ_TIME;             % Right Ascension (radians) 
RA         = rad2deg_(RRA);                   % Right Ascension (degrees) 
GHA_ARIES  = 15*(GMST + HOUR);                % Greenwich Hour Angle Aries (degrees) 
GHA        = GHA_ARIES - RA;                  % Greenwich Hour Angle Sun (degrees) 
GHA        = rem(GHA,360); 
SHA        = 360 + RA;                        % Siderial  Hour Angle Sun (degrees) 
SUN_LONG   = MEAN_RA + P17./3600;             % Sun Longitude 
 
if EQ_TIME > 180 
  EQ_TIME =  EQ_TIME - 360; 
end 
 
% The following four quantities are returned: 
% print *, 'Equation of Time  ', eq_time 
% print *, 'Declination       ', decl 
% print *, 'Earth-Sun Dist    ', dist 
% print *, 'GHA               ', gha 
 
% The following tutorial quantities are not returned in this subroutine: 
% print *, 'Mean Right Ascen  ', mean_ra 
% print *, 'Right Ascension   ', dmod(ra,circle) 
% print *, 'GHA Aries         ', dmod(gha_aries,circle) 
% print *, 'Sun Longitude     ', dmod(sun_long,circle) 
% print *, 'Siderial Hr Angle ', dmod(sha,circle)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [altitude,azimuth] = altazm_(decl,gha,lat,long) 
% ALTAZM_ - Celestial altitude and azimuth from declination,  
%           Greenwich Hour Angle, latitude and longitude. 
% 
% Use As:  [altitude,azimuth] = altazm_(decl,gha,lat,long) 
% 
% Input:    decl     = Declination          (degrees) 
%           gha      = Greenwich Hour Angle (degrees) 
%           lat      = Latitude             (decimal degrees +N/-S) 
%           long     = Longitude            (decimal degrees +W/-E) 
% Output:   altitude = Solar altitude       (degrees) 
%           azimuth  = Solar azimuth        (degrees) 
% 
% Example:  [alt,azm] = altazm_(-0.42,118.13,36,121) 
%                       alt = 53.48; azm = 175.17 
% 
% See Also: ALMANAC_, SUNANG2_, EPHEMS_ 
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% 22 Jul 1997; W. Broenkow 
 
% Much effort was required to adapt this from sunang2_.m 
% which has a quadrant ambiguity. The algorithm here is 
% taken from Wilson (1980). 
 
Nargs = nargin; 
if Nargs == 0 
  help altazm 
  return 
end 
 
latR  = deg2rad_(lat);                 % latitude    in radians 
lonR  = deg2rad_(long);                % longitude   in radians 
decR  = deg2rad_(decl);                % declination in radians 
GHAR  = deg2rad_(gha);                 % GHA         in radians 
LHAR  = GHAR - lonR;                   % LHA         in radians 
 
% solar zenith angle   theta-0 
zenR  = acos(sin(latR).*sin(decR) + cos(latR).*cos(decR).*cos(LHAR)); 
 
% solar azimuth angle  Phi-0 
% sunang2_.m algorithm 
% azmR  = asin(sin(abs(LHAR)).*cos(decR)./sin(zenR)); 
 
% wilson_ (1980) Algorithm 
azmR  = acos((-sin(latR).*cos(LHAR).*cos(decR) + sin(decR).*cos(latR))./sin(zenR)); 
 
sign1 = sign(-cos(decR).*sin(LHAR)./sin(zenR)); 
sign2 = sign(sin(azmR)); 
indx = find(sign1 ~= sign2); 
 
if exist('indx') 
  azmR(indx) = 2*pi - azmR(indx); 
end 
 
zenith   =  rad2deg_(zenR); 
altitude =  90 - zenith; 
azimuth  =  rad2deg_(azmR); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function JD = juldate_(YEAR,MONTH,DAY) 
% JULDATE_ - Converts calendar date (Year, Month, Day) to Julian date 
% 
% Use As:   juldate2(y,mo,d) 
% 
% Input:    y  = integer year as in 1979 
%           mo = integer month 1..12 
%           d  = integer day of month 1..31 
% Output:   days from julian date 0: sometime in 4712 BC  
%       
% Example:  juldate_(1979,6,14) -> 2444038.500 
% 
% Note:     Notice that the julian date starts at noon_, not midnight. 
%           This function returns the date starting at midnight, hence 
%           the fractional date in this example.  
% 
% See Also: CALDATE_, DEG2DMS_, DMS2DEG_, JULDAY_, JULSEC_, NOW_ 
%           JULDATE_.MEX
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% Adapted from NUMERICAL RECIPIES but vectorized                    
% 14 Nov 97; WWB start from scratch
 
if nargin == 2 & size(YEAR) ~= size(MONTH) 
  disp('JULDATE Error: input vectors must be same size') 
  return 
end 
if nargin == 3 & size(YEAR) ~= size(DAY) 
  disp('JULDATE Error: input vectors must be same size') 
  return 
end 
[a,b] = size(YEAR); 
if a > b 
  YEAR = YEAR'; 
 if nargin > 1 
  MONTH = MONTH'; 
 end 
 if nargin > 2 
  DAY = DAY'; 
 end 
end 

IGREG = 588829; % date when gregorian calendar was adopted

ix = find(YEAR < 0);
YEAR(ix) = YEAR(ix) + 1;

JY = YEAR  - 1;
JM = MONTH + 13;

ix = find(MONTH > 2);
JY(ix) = YEAR(ix);
JM(ix) = (MONTH(ix) + 1);

JD = floor(365.25*JY) + floor(30.6001*JM) + DAY + 1720995 - 0.5;  % Numerical Recipies does not use the 0.5

JJ = (DAY + 31*(MONTH + 12*YEAR));

ix = find(JJ >= IGREG);
JA(ix) = floor(0.01*JY(ix));
JD(ix) = JD(ix) + 2 - JA(ix) + floor(0.25*JA(ix));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yyy = dms2deg_(in_angle,option,minutes,seconds) 
% DMS2DEG_ - Convert degrees, minutes, seconds (DD.MMSSss) to  
%            decimal equivalent DD.ddd 
% 
% Use As:   out_angle = dms2deg_(in_angle,option,minutes,seconds) 
% 
% Input:    in_angle  = angle in the form DD.MMSS or HH.MMSS 
%           option    = 0 for inputs of  DD.MMm 
%           option    = 1 for inputs of  DD.MMSSs 
%           option    = 2 for inputs of  DD, MM.m 
%           option    = 3 for inputs of  DD, MM,  SS.s 
%           minutes   for option 2 is decimal minutes 
%                         option 3 is whole minutes 
%           seconds   for option 3 is decimal seconds  
% Output:   out_angle = Angle or Time as decimal fraction 
%   
% Example:  dms2deg_( 36.59599,0)           ->  36.99331  DD.ddddd 
%           dms2deg_(-36.59599,1)           -> -36.99997  DD.ddddd 
%           dms2deg_(-36,      2, 59.599)   -> -36.99 
%           dms2deg_(-36,      3, 59, 59.9) -> -36.999972  
% 
% Note:     When using option = 2 or option = 3; the absolute value of the 
%           minutes and seconds are used. 
%           A warning message is issued if minutes or seconds > 60. 
% 
% See Also: DEG2DMS_
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% adapted from DTO_DEG.FOR 
% Notice the addition of a small angle to foil truncation errors 
% 22 Jan 96; W. Broenkow 
% 25 Mar 96; W. Broenkow vectorized 
%  7 Apr 96; W. Broenkow added DD MM SS.s modes 
% 23 May 96; W. Broenkow debugged negative inputs 
% 17 Jun 96; W. Broenkow changed form to option, min to minutes
% 11 Oct 2000; S. Flora - fixed so option 2 and 3 work correctly when in_angle = 0

if nargin == 0
  help dms2deg_
  return
end

if nargin == 1 
  option = 0;                   % default DD.MMmm 
end 
 
if (option == 0)                % from DD.MMmmm to DD.dddd 
  P1 = in_angle + sign(in_angle)*2e-14; 
  P2 = fix(P1);               % degrees (hours) 
  P3 = P1 - P2;               % .minutes 
  P4 = 100.*P3;               % minutes             
  if any(P4 > 60.0); 
    disp('DMS2DEG_ BAD PARAMETER Minutes > 60......................................................') 
  end 
  yyy = P2 + P4/60.0;        
  return 
elseif (option == 1)            % from DD.MMSSs to DD.dddd 
  P1 = in_angle + sign(in_angle)*2e-14; 
  P2 = fix(P1);               % degrees (hours) replace aint w/ fix 
  P3 = (P1 - P2);             % .minutesseconds 
  P4 = 100*P3;                % minutes.seconds 
  P5 = fix(P4);               % minutes 
  P6 = P4 - P5;               % .seconds 
  P7 = 100*P6;                % seconds  
elseif (option == 2) 
  if nargin < 3 
    disp('DMS2DEG_ requires input argument for minutes') 
    disp('i.e. dms2deg_(D,2,M)') 
    return 
  end 
  P2 = abs(in_angle); 
  P5 = abs(minutes); 
  P7 = 0; 
elseif (option == 3) 
  if nargin < 4 
    disp('DMS2DEG_ requires input arguments for minutes and seconds') 
    disp('i.e.  dms2deg_(D,3,M,S)') 
    return 
  end 
  P2 = abs(in_angle); 
  P5 = abs(minutes);
  P7 = abs(seconds); 
end 
if any((abs(P5) >= 60.0) | (abs(P7) > 60.0)) 
  disp('DMS2DEG_  BAD PARAMETER minutes > 59 or seconds > 60') 
end 

yyy = P2 + P5./60.0 + P7./3600.0;
 
if in_angle < 0
  yyy = abs(yyy).*(-1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = deg2rad_(x) 
% DEG2RAD_ - Converts angle in degrees to radians  
%  
% Use As:   deg2rad_(angle) 
% 
% Input:    angle in degrees 
% Output:   angle in radians 
% 
% Example:  deg2rad_(270) -> 4.7124 
% 
% See Also: RAD2DEG_  
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% 28 Jan 1996; W. Broenkow 
% 20 Jun 96; changed name wwb 
 
y = x*pi/180; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = rad2deg_(x) 
% rad2deg_ - Converts angle in radians to degrees 
%  
% Use As:   rad2deg_(angle) 
% 
% Input:    angle in radians 
% Output:   angle in degrees 
% 
% Example:  rad2deg_(3*pi/2) -> 270 
% 
% See Also: DEG2RAD_
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% 28 Jan 1996; W. Broenkow 
% 20 Jun 96; changed name wwb 
 
y = x*180/pi;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Year,Month,Day,Hour,Minute,Second,DOY,HR] = timeupk_(T1,option) 
 
% timeupk_ - Unpack time from the YYMMDD.HHMMSS or DD.MMYY packed form into Y,M,D,H,M,S 
% 
% Use As:  [Year,Month,Day,Hour,Minute,Second,DOY,HR] = timeupk_(T1,option) 
%  
% For Option = 0: 
% Inputs:   T1     = time as YYMMDD.HHMMSS 
% Output:   Year   = 4-digit year 
%           Month  = month (1-12) 
%           Day    = day of month (1-31) 
%           Hour   = hour of day (0-23) 
%           Minute = minute (0-59)  
%           Second = second (0-59) 
%           DOY    = julian day of year (1-365) 
%           HR     = decimal hour of day (1-24) 
%           for option = 1: 
%           Year   = 4-digit year 
%           Month  = month (1-12) 
%           Day    = day of month (1-31) 
% For Option = 1: 
% Input:    T1     = time as DD.MMYY 
% Output:   Year   = 4-digit year 
%           Month  = month (1-12) 
%           Day    = day of month (1-31) 
%
% Example: [Year,Month,Day,Hour,Minute,Second,DOY,HR]   = timeupk_(891017.170400,0) ->
%           89   10    17  17   04     00     290 17.07
%          [Year,Month,Day,Hour,Minute,Second,DOY,HR] = timeupk_(17.1089,1) ->
%           1989 10    17  NaN  NaN    NaN    290 NaN  
% 
% See Also: CALDATE_, DEG2DMS_, DMS2DEG_, JULDAY_, TIMEPK_ 
% NOTE: This function is not vectorized. Use it in a loop. 
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% 06 Jun 96; W.Broenkow 
% 09 Sep 96; S. Flora - changed to use a string to unpack the time 
%                       this eliminates the rounding errors. So that 
%                       941114.220000 does not return Hour = 21, Min = 99 
%                       sec = 99.99997206032276 
% 21 Dec 96; W. Broenkow added unpacking of DD.MMYY 
 
if nargin < 2 
 option = 0; 
end 
 
% The trouble with the string method is that it is not vectorized 
 
if option == 0  
  if isnan(T1) 
    Year  = NaN; 
    Month = NaN; 
    Day   = NaN; 
    Hour  = NaN; 
    Minute = NaN; 
    Second = NaN; 
    DOY    = NaN; 
    Hour   = NaN; 
    return 
  end 
  str    = sprintf('%013.6f',T1); 
  Year   = str2num(str(1:2)); 
  Month  = str2num(str(3:4)); 
  Day    = str2num(str(5:6)); 
  Hour   = str2num(str(8:9)); 
  Minute = str2num(str(10:11)); 
  Second = str2num(str(12:13)); 
 
  DOY = julday_(Year,Month,Day,0); 
  HR  = Hour + Minute/60 + Second/3600; 
else 
  if isnan(T1) 
    Year  = NaN; 
    Month = NaN; 
    Day   = NaN; 
    Hour  = NaN; 
    Minute = NaN; 
    Second = NaN; 
    DOY    = NaN; 
    Hour   = NaN; 
    return 
  end 
  str    = sprintf('%07.4f',T1); 
  Year   = str2num(str(6:7)); 
  Month  = str2num(str(4:5)); 
  Day    = str2num(str(1:2)); 
  Hour   = NaN; 
  Minute = NaN; 
  Second = NaN; 
  if Year < 30 
    Year = 2000 + Year; 
  else 
    Year = 1900 + Year; 
  end 
  DOY = julday_(Year,Month,Day,0); 
  HR  = NaN; 
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [jd,dow] = julday_(y,mo,d,option) 
% JULDAY_ - Determine the julian day of the year (1..365) or julian date 
% 
% Use As:   [jday,dow] = julday_(y,mo,d,option) 
% 
% Input:    y  = year         30..99 is used as 1930..  
%                             00..29 is used as 2000.. 
%                             1582 any 4-digit year         
%           mo = month        1..12 
%           d  = day of month 1..31 
%           option = 0 or 1  day of year or Julian Date option 
% Output:   jday = days from 0 Jan current year if option = 0 Julian day of year 
%      or   jday = days from 4713 BC            if option = 1 Julian Date 
%           dow  = day of week 1 = Sunday  2 = Monday .. 7 = Saturday 
% 
% Example:  julday_(79,06,14,0)    -> 165 
%           julday_(79,06,14,1)    -> 24440385 
%           julday_(1582,10,15,1)  -> 2299170.5 
%           julday_(1968,5,23.5,1) -> 2440000.0 
%           julday(1582,10,16) - julday_(1582,10,15) -> -9 !!  
%              This is because the Julian calendar 
%              started at noon_ on 15 October 1582 
%              Julian Date 0 is 4712 BC 
% 
% See Also: CALDATE_, DEG2DMS_, DMS2DEG_, HELPDATE, JULDATE_, JULSEC_, NOW_ 
 
% Copyright (c) 1996 by Moss Landing Marine Laboratories 
% Uses function juldate_.m 
% 15 Apr 95; W. Broenkow 
% 29 Jan 96; replaced juldate_.mex by juldate_.m 
% 25 Mar 96; WWB vectorized added day of week, debugged 
% 21 Jun 96; WWB changed help  
 
if ~exist('option'), 
option = 0; 
end; 
if (nargin < 3),  
  jd = NaN;  
  disp('ERROR:  JULDAY_ Must input 3 arguments:  day, month, year') 
else 
  % Parse shorthand years 
  yr = y;                              % do all years first 
  ifind1 = find(y < 30*ones(size(y)) & y >= zeros(size(y)));   % then correct shorthand years  
  if exist('ifind1') 
    yr(ifind1) = 2000 + y(ifind1);   % up to 2030 
  end 
  ifind2 = find(y >= 30*ones(size(y)) & y < 100*ones(size(y)));
 
  if exist('ifind2') 
    yr(ifind2) = 1900 + y(ifind2);           % from 1931..1999 
  end 
end; 
% look for Julian Date option 
if (option > 0), 
  jd  = juldate_(yr,mo,d); 
  dow = 1+rem(floor(jd)+2,7); 
else 
  jd1  = juldate_(yr,ones(size(yr)),zeros(size(yr)));  % julian date for first of given year 
  jd2  = juldate_(yr,mo,d);                            % julian date 
  jd   = jd2 - jd1;              % number of days from Jan 0 = day of year 
  dow  = 1+ rem(floor(jd2+2),7); % day of the week 1=Sunday 
end 
