function [WAVE] = AWACwvs2nc(baseFile,metadata)
%function WAVE = AWACwvs2nc(baseFile,metadata)
%
%   Companion function to AWAC2nc that will take processed waves data from
%   an AWAC and put calculated parameters into a netcdf File. Will also
%   create netCDF files of wave spectra, if these were output. At a
%   minimum, the .wap wave parameters file is necessary. Use the same
%   metafile for user entered metadata as with the currents processing in
%   AWAC2cdf.m - Outputs a structure of the wave spectral data
%
%   inputs - baseFile - base file name for AWAC output data
%
%            metadata - basic netCDF file metadata from the runmeta file
%
%   KJR, USGS, 2007.7.24

if nargin<1
    help(mfilename)
    return
end

wap = dir('*.wap');

if length(wap)==1
    wapfile = wap.name;
    data = load(wapfile);[N,M] = size(data);
elseif length(wap)<1
    error('There are no .wap wave parameters files in this directory')
elseif length(wap)>1
    wapfile = [baseFile '.wap'];
    data = load(wapfile);[N,M] = size(data);
end
% l = findstr('.',wapfile);awacfile = wapfile(1:l-1);
WAVE.metadata = metadata;
WAVE.instmeta = readWVSHeader(baseFile);
WAVE.filename = baseFile;

%get julian day and find data between deploy and recover times
JD = julian([data(:,3) data(:,1) data(:,2) data(:,4) data(:,5) data(:,6)]);
deploy = datenum2julian(datenum(metadata.Deployment_date));
recover = datenum2julian(datenum(metadata.Recovery_date));
ind = find(JD>=deploy & JD<=recover);
WAVE.data = data(ind,:);

%now move time to center of ensemble
instmeta.WaveSampleDuration = WAVE.instmeta.WaveNumberOfSamples/WAVE.instmeta.WaveSampleRate;
jd = JD(ind) + (instmeta.WaveSampleDuration/86400)/2;
WAVE.time = floor(jd); WAVE.time2 = (jd - floor(jd))*86400000;
WAVE.burst = ind;
%perform time correction if specified in metafile
if isfield(metadata,'TimeCorrection')
    newjd = jd + metadata.TimeCorrection/24;
    WAVE.time = floor(newjd); WAVE.time2 = (newjd - floor(newjd))*86400000;
else
end

%now look for the spectral parameters files
was = dir('*.was');if ~isempty(was),for i = 1:length(was),names{i} = was(i).name;end,end
wasfile = [baseFile '.was'];
if length(was)>1 || sum(strcmp(wasfile,names))>0
    disp('Reading in PUV directional wave parameters')
%     wasfile = was.name;
    WAS = load(wasfile);WDR = load([baseFile '.wdr']);
    WAVE.amp = WAS(ind+1,:);WAVE.dir = WDR(ind+1,:);
    WAVE.freq = WAS(1,:);m = length(WAVE.freq);
%     %convert results from m^2/Hz to mm^2/Hz
%     WAVE.amp = WAVE.amp*1000;
elseif isempty(was)
    disp('There are no .was Spectra files in this directory')
end

wds = dir('*.wds');if ~isempty(wds),for i = 1:length(wds),names{i} = wds(i).name;end,end
wdsfile = [baseFile '.wds'];
if length(wds)>1 || sum(strcmp(wdsfile,names))>0
    disp('Reading in Full Directional Wave spectra - this may take some time...')
%     wdsfile = wds.name;
    WDS = load(wdsfile);n = floor(length(WDS)/m);
    if n~=N,disp('Number of spectral bursts does not equal wave bursts in .wap file'),end
    dspec = reshape(WDS,n,m,90);WAVE.dspec = dspec(ind,:,:);
%     %convert results from m^2/Hz to mm^2/Hz
%     WAVE.dspec = WAVE.dspec*1000;
elseif isempty(wds)
    disp('There are no .wds Full Directional Spectra files in this directory')
end

disp('Writing processed wave parameters to netCDF...')
%write out the waves parameters and spectral output
ncwrite_AWACWvsp(WAVE);
clear WAS WDS
%now read in the raw data
wad = dir('*.wad');if ~isempty(wad),for i = 1:length(wad),names{i} = wad(i).name;end,end
wadfile = [baseFile '.wad'];
if length(wad)>1 || sum(strcmp(wadfile,names))>0
    disp('Reading in raw velocity and pressure burst timeseries - this will definitely take some time...')
    WAD = load(wadfile);[q,r] = size(WAD);
    samp = WAVE.instmeta.WaveNumberOfSamples;n = q/samp;
    if n~=N, 
        disp('Number of spectral bursts does not equal wave bursts in .wap file; cropping extra data off end of file')
        Q = N*samp;if Q>q,N= N-1;Q=N*samp;end       %the instrument may have cut out in the middle of a wave burst
        WAD = WAD(1:Q,:);
    end
    rawdata = reshape(WAD,N,samp,r);WAVE.rawdata = rawdata(ind,:,:);
elseif isempty(wad)
    disp('There are no .wad raw data files in this directory');
end

%now read in the raw surface track data
wst = dir('*.wst');if ~isempty(wst),for i = 1:length(wst),names{i} = wst(i).name;end,end
wstfile = [baseFile '.wst'];
if length(wst)>1 || sum(strcmp(wstfile,names))>0
    disp('Reading in raw acoustic surface track burst timeseries - this will definitely take some time...')
    wstfile = [baseFile '.wst'];
    WST = load(wstfile);[q,r] = size(WST);
    samp = WAVE.instmeta.WaveNumberOfSamples;n = q/samp;
    if n~=N, 
        disp('Number of spectral bursts does not equal wave bursts in .wap file; cropping extra data off end of file')
        Q = N*samp;if Q>q,N= N-1;Q=N*samp;end       %the instrument may have cut out in the middle of a wave burst
        WST = WST(1:Q,:);
    end
    ast = reshape(WST,N,samp,r);WAVE.ast = ast(ind,:,:);
elseif isempty(wst)
    disp('There are no .wst raw data files in this directory')
end

if isfield(WAVE,'rawdata') || isfield(WAVE,'ast')
    ncwrite_AWACWvsr(WAVE)
end