function xy = readCDIPxy(fname)
% function xy = readCDIPxy(fname)
%
% T. Hilmer, derived from J. Aucan's "cdipxyz2mw"
%
% hdr:  header string
% time: time of original samples. corrected, see below
% x:    North (cm?)
% y:    East (cm?)  different from text file; y is West
% z:    vertical displacement (cm?)
% Hz:   sampling frequency, read from header
% declination:  magnetic declination, read from header
% percnt_good:  percent of vectors that are good. read from header
% H:    individual wave heights from zero up-crossing method
% P:    period (s) of each aforementioned wave height
% T:    centered time of each aforementioned wave height
%
% !!! WARNING !!!
% this fcn has a forced correction to the time.  The data stored in the
% text file has 1 sec resolution, whereas the sampling rate for the Waimea
% buoy is 1.28 Hz.  This is evident in the repeated times for sample
% points within the text file.  Sampling rate is read from header and
% correction applied.
%
% corrects for magnetic declination, read from header. might want to check
% that I'm doing this right. for a 10 E declination, the corrected vectors
% are ccw of the original
%
% changes centimeters to meters

% fname = '/media/MAIN/D/data/waimea_buoy/xy10601200112160007';

try
    drr = dir(fname);
    if isempty(drr), error([fname ': not found']), end
    [fid,stat] = fopen(fname);
    if fid == -1, error([fname ': ' stat]), end

    xy.hdr = '';
    for n = 1:20 % 20 header lines
        xy.hdr = strvcat(xy.hdr,fgetl(fid));
    end

    data = textscan(fid,'%f');
    data = data{1};

    if rem(length(data),4)
        warning('found incorrect number of data points')
    end

    data = reshape(data,4,length(data)/4)';

    xy.time = datenum(num2str(data(:,1)),'yyyymmddHHMMSS');
    xy.x = data(:,2);
    xy.y = -data(:,3);% !!!! EAST POSITIVE !!!
    xy.z = data(:,4);

    fcn = @(x)strfind(x,':')+2;

    xy.name = xy.hdr(1,fcn(xy.hdr(1,:)):end);
    xy.station = sscanf(xy.hdr(2,:),'%*s %f');
    tmp = sscanf(xy.hdr(3:4,:)','%*s %*s %f %f %*s %s');
    xy.lat = tmp(1) + tmp(2)/60;
    xy.lon = tmp(4) + tmp(5)/60;
    if strcmp('S',char(tmp(3))), xy.lat = -xy.lat;
    elseif ~strcmp('N',char(tmp(3))), error('N or S in latitude??'), end
    if strcmp('W',char(tmp(6))), xy.lon = -xy.lon;
    elseif ~strcmp('E',char(tmp(6))), error('E or W in longitude??'), end
    xy.water_depth = sscanf(xy.hdr(5,:),'%*s %*s %f',1);
    xy.declination = sscanf(xy.hdr(6,:),'%*s %*s %*s %f');
    tmp = sscanf(xy.hdr(6,:),'%*s %*s %*s %*s %s');
    if strcmp(tmp,'E'), elseif strcmp(tmp,'W')% The declination is positive when the magnetic north is east of true north
        xy.declination = - xy.declination;
    else warning('could not determine magnetic declination from header'), end
    %     xy.data_type = xy.hdr(7,fcn(xy.hdr(7,:)):end);
    %     xy.gauge_type = xy.hdr(8,fcn(xy.hdr(8,:)):end);
    xy.rate = sscanf(xy.hdr(9,:),'%*s %*s %f');

    tmp = sscanf(xy.hdr(14,:),'%*s %*s %s');  xy.start = datenum(tmp,'yyyymmddHHMMSS');
    tmp = sscanf(xy.hdr(15,:),'%*s %*s %s');  xy.stop = datenum(tmp,'yyyymmddHHMMSS');
    tmp = sscanf(xy.hdr(16,:),'%*s %*s %s');
    tmp = str2double({tmp([1 2]) tmp([4 5]) tmp([7 8])});
    xy.sample_time = tmp(1)*3600 + tmp(2)*60 + tmp(3);
    xy.nvectors = sscanf(xy.hdr(17,:),'%*s %*s %*s %*s %f');
    xy.percnt_good = sscanf(xy.hdr(18,:),'%*s %*s %f');


    %----------------------
    % Modifications: up to this point, everything is as CDIP made it
    %----------------------

    %----------------------
    % Centimeters to meters
    %----------------------
    xy.x = xy.x/100; xy.y = xy.y/100; xy.z = xy.z/100;

    %----------------------
    % Forced correction to time
    %----------------------
    if ~isnumeric(xy.rate), warning('could not determine sampling rate from header'), end
    t = xy.time(1) + 1/86400/xy.rate*(0:length(xy.time)-1)';
    %     figure, plot(xy.time,'or'), hold on, plot(t,'.k')
    %     figure, plot(diff(xy.time)*86400,'or'), hold on, plot(diff(t)*86400,'.k')
    xy.time = t;

    %---------------------------------
    % Correct for magnetic declination
    %---------------------------------
    mag = sqrt(xy.x.^2 + xy.y.^2);% magnitude
    ang = atan2(xy.y,xy.x);% angle
    ang = ang + xy.declination * pi/180;% correct for declination
    xy.x = mag .* cos(ang);
    xy.y = mag .* sin(ang);

    %---------------------------------
    % Calcuate some basic wave statistics using the zero-crossing method
    %---------------------------------
    [xy.H,xy.P,xy.T] = waveStats(xy.z,xy.time);

    if fclose(fid), error(['error closing: ' fname]), end

catch
    try fclose(fid); end
end

function [H,P,T] = waveStats(z,t)
%
% H:    wave heights, from zero up-crossing method
% P:    period of each wave
% T:    time (centered) of each wave

N = length(z);

zc = find(diff(z <= 0) == -1);% using +1 here would find down-crossings

if zc(1)~=1,   zc = [1; zc]; end
if zc(end)~=N, zc = [zc; N]; end

%     plot(z(1:100),'.-k')
%     hold on
%     plot(zc(1:10),z(zc(1:10)),'or')

Nw = length(zc)-1;% # of waves

H = nan(Nw,1); P = H; T = H;

for n = 1:Nw
    ind = zc(n):zc(n+1);
    H(n) = max(z(ind)) - min(z(ind));
    P(n) = 86400* (t(ind(end)) - t(ind(1)));
    T(n) = mean(t(ind([1 end])));% don't average across all times. would be incorrect for sporadic sampling.
end

%  function [mwh,swh,tswh,H,tH] = cdipxyz2mw(filelist);
%     name=textread(filelist,'%s');
%
%     n_name = length(name);
%
%     H = NaN*ones(800,n_name); %wave height from upcrossings
%     tH = H;                     %time of wave heights (centered)
%     mwh = NaN*ones(n_name,1);     %max wave height
%     swh = mwh;                  %significant wave height
%     tswh = mwh;                 %time of significant wave height
%
%
%     for j = 1:n_name
%         [tz,x,y,zz]=textread(char(name(j)),'%n %n %n %n', 'headerlines', 21);
%         clear x y
%         %get time in Julian days
%         yr = floor(tz/1e10);
%         mm = floor((tz-yr*1e10)/1e8);
%         dd = floor((tz-yr*1e10-mm*1e8)/1e6);
%         hh = floor((tz-yr*1e10-mm*1e8-dd*1e6)/1e4);
%         mi = floor((tz-yr*1e10-mm*1e8-dd*1e6-hh*1e4)/1e2);
%         ss = floor((tz-yr*1e10-mm*1e8-dd*1e6-hh*1e4-mi*1e2)/1);
%         date= datenum(yr,mm,dd,hh,mi,ss);
%
%         %compute zero upcrossing wave heights
%         m = length(zz);
%         k = find(zz(1:m-1) < 0 & zz(2:m) >=0);
%         nk = length(k)-1;
%         for ik = 1:nk
%             H(ik,j) = range(zz(k(ik):k(ik+1)));
%             T(ik,j) = k(ik+1)-k(ik);
%             tH(ik,j) = mean(date(k(ik):k(ik+1)));
%         end
%         mwh(j) = nanmax(H(:,j));
%         Hs = sort(H(1:nk,j));
%         swh(j) = mean(Hs(floor(2*nk/3):nk));
%         tswh(j) = mean(date);
%     end