% ttide_adcp_indiv.m  An mfile to run T_TIDE on South Carolina ADCP 
%                     velocity and pressure data.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% C. Sullivan 05/15/06, Water flow data in bins likely exposed to air has
%                       been set to the NetCDF FillValue_.  Editing code to
%                       consider the surface bin the upward-most good bin.
% C. Sullivan 01/11/05, Use data from surface, bottom, and user-selected
%                       mid-depths.  Run the analysis on each 3-month
%                       timeseries using that timeseries' start time and
%                       time step.

more off

% Define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\PLOTS\TTIDE';

% Define adcp data filenames.  
dataFiles={'7201adc-a.nc','7401adc-a.nc',...
           '7221adc-a.nc','7421adc-a.nc',...
           '7241adc-a.nc','7441adc-a.nc',...
           '7331adc-a.nc','7531adc-a.nc',...
           '7321adc-a.nc','7521wh-a.nc'};

for f = 1:length(dataFiles)
    
    nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    
    %get the site, the latitude, and the depth of
    %each velocity bin
    site=nc.DESCRIPT(6);
    if strcmp(site,'7');
        site='7B';
    elseif strcmp(site,'') & f==8
        site='4';
    end
    latitude=nc.latitude(:);
    binDepths=nc{'depth'}(:);
    
    %get time in julian days
    jd=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
     
    %get bottom pressure and convert it to height
    %with H=P/(rho*g) where rho=1025 km/m^3 and
    %g=9.8 m/s^2
    P_4=nc{'P_4'}(:);
    P_4(P_4>=1e35)=nan; %nan bad points
    Height=P_4/(1025*9.8);
    
    %sites 3 and 7 have bad pressure during deployment 1.
    %we will not perform harmonic analysis on this data.
    if f==5  %mooring 724
        Height=[];
    elseif f==9  %mooring 732
        Height=[];
    end

    %%%%%%run t_tide on pressure data (bottom pressure)
    if ~isempty(Height)
        startInd = find(~isnan(Height),1,'first');
        startTime = julian2datenum(jd(startInd));
        interval = mean(diff(julian2datenum(jd)))*24;
        ttideOutFile = fullfile(plotDirectory,['ttide_pres_',dataFiles{f}(1:3),'.txt']);
        [nameu,fu,tidecon,xout]=t_tide(Height,'start time',startTime,...
                                              'interval',interval,...
                                              'latitude',latitude,...
                                              'output',ttideOutFile);
    end

    %%%%%%run t_tide on velocity data at the bottom
    ubot = [nc{'u_1205'}(:,1)];
    vbot = [nc{'v_1206'}(:,1)];
    bads=find(ubot+vbot >= 1e35);
    ubot(bads) = nan;
    vbot(bads) = nan;

    uvbot = ubot + i*vbot;

    startInd = find(~isnan(uvbot),1,'first');
    startTime = julian2datenum(jd(startInd));
    interval = mean(diff(julian2datenum(jd)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_',dataFiles{f}(1:3),'_bottom.txt']);
    [nameu,fu,tidecon,xout]=t_tide(uvbot,'start time',startTime,...
                                         'interval',interval,...
                                         'latitude',latitude,...
                                         'output',ttideOutFile);

    %%%%%%run t_tide on velocity data at the surface
    %%%%%%C. Sullivan, 05/15/06: the surface is the upward-most good bin
    questBins = nc{'u_1205'}.bins_questioned(:);
    if isempty(questBins)
        lastGoodBin = size(nc{'u_1205'}(:),2);
    else
        lastGoodBin = questBins(1)-1;
    end
    depthLastGoodBin = binDepths(lastGoodBin)
    usfc = [nc{'u_1205'}(:,lastGoodBin)];
    vsfc = [nc{'v_1206'}(:,lastGoodBin)];
    bads=find(usfc+vsfc >= 1e35);
    usfc(bads) = nan;
    vsfc(bads) = nan;

    uvsfc = usfc + i*vsfc;

    startInd = find(~isnan(uvsfc),1,'first');
    startTime = julian2datenum(jd(startInd));
    interval = mean(diff(julian2datenum(jd)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_',dataFiles{f}(1:3),'_surface.txt']);
    [nameu,fu,tidecon,xout]=t_tide(uvsfc,'start time',startTime,...
                                         'interval',interval,...
                                         'latitude',latitude,...
                                         'output',ttideOutFile);
                                   
    %%%%%%run t_tide on velocity data at user-specified mid-depths
    targetDepths=[4,6]; %mbs
    for b=1:size(targetDepths,2)
        [index_num,dist_to_hit] = value2Index(binDepths,targetDepths(b),0.5);

        u = [nc{'u_1205'}(:,index_num)];
        v = [nc{'v_1206'}(:,index_num)];
        bads=find(u+v >= 1e35);
        u(bads) = nan;
        v(bads) = nan;

        uv = u + i*v;
        
        startInd = find(~isnan(uv),1,'first');
        startTime = julian2datenum(jd(startInd));
        interval = mean(diff(julian2datenum(jd)))*24;
        ttideOutFile = fullfile(plotDirectory,['ttide_curr_',dataFiles{f}(1:3),'_',num2str(targetDepths(b)),'mbs.txt']);
        [nameu,fu,tidecon,xout]=t_tide(uv,'start time',startTime,...
                                           'interval',interval,...
                                           'latitude',latitude,...
                                           'output',ttideOutFile);
    end

    nc=close(nc);
    
   
end