% ttide_argn_comb.m  An mfile to run T_TIDE on South Carolina Argonaut 
%                    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.
%                       Also found and fixed a few typos that would have
%                       produced weird results.
% C. Sullivan 01/31/05, This version use velocity data from surface,
%                       bottom, and user-selected mid-depths.  The pressure
%                       used is the bottom pressure (P_4).  The harmonic
%                       analysis is run using the entire timeseries
%                       thus combining the data from deployments 1 and 2.


more off

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

% Define Argonaut data filenames.  
dataFiles={'7263arc-a.nc','7463arc-a.nc'};

% The following loop loads the entire timeseries of 
% pressure and velocities at site 6 and runs harmonic
% analysis on the data.
for f = 2:2:length(dataFiles)
    
    nc1=netcdf(fullfile(dataDirectory,dataFiles{f-1}), 'nowrite');
    nc2=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    
    site=nc1.DESCRIPT(6);
    latitude=nc1.latitude(:);
    binDepths1=nc1{'depth'}(:);
    binDepths2=nc2{'depth'}(:);
    
    %time in julian days
    jd1=nc1{'time'}(:) + (nc1{'time2'}(:)/3600/1000/24);
    jd2=nc2{'time'}(:) + (nc2{'time2'}(:)/3600/1000/24);
    rawTime=[jd1;jd2];
    
    %time intervals for deployments 1 and 2 should
    %be 15 minutes each
    interval1=mean(diff(jd1));
    interval2=mean(diff(jd2));
    if abs(interval1-interval2)>15/60/24
        disp('Time intervals don''t match!')
    end
    interval=interval1;
    
    %will interpolate data onto an array that is time
    %in julian days.
    interpTime=rawTime(1):interval:rawTime(end);
    
    %%%%%%run t_tide on pressure data (bottom pressure)
    %bottom pressure converted to height using H=P/(rho*g)
    %where rho=1025 km/m^3 and g=9.8 m/s^2
    P_41=nc1{'P_4'}(:);
    P_41(P_41>=1e35)=nan; 
    Height1=P_41/(1025*9.8);
    
    P_42=nc2{'P_4'}(:);
    P_42(P_42>=1e35)=nan; 
    Height2=P_42/(1025*9.8);
    
    rawHeight=[Height1;Height2];
    
    %interpolate data  
    [interpHeight]=smart_interp(rawTime,rawHeight,interpTime,3); 
    
    %%%%%%run t_tide on pressure data (bottom pressure)
    if ~isempty(interpHeight)
        startInd = find(~isnan(interpHeight),1,'first');
        startTime = julian2datenum(interpTime(startInd));
        interval = mean(diff(julian2datenum(interpTime)))*24;
        ttideOutFile = fullfile(plotDirectory,['ttide_pres_site',site,'.txt']);
        [nameu,fu,tidecon,xout]=t_tide(interpHeight,'start time',startTime,...
                                                    'interval',interval,...
                                                    'latitude',latitude,...
                                                    'output',ttideOutFile);
    end

    %%%%%%run t_tide on velocity data at the bottom
    ubot1 = nc1{'u_1205'}(:,1);
    vbot1 = nc1{'v_1206'}(:,1);
    bads=find(ubot1+vbot1 >= 1e35);
    ubot1(bads) = nan;
    vbot1(bads) = nan;
    
    ubot2 = nc2{'u_1205'}(:,1);
    vbot2 = nc2{'v_1206'}(:,1);
    bads=find(ubot2+vbot2 >= 1e35);
    ubot2(bads) = nan;
    vbot2(bads) = nan;
    
    raw_uvbot=[ubot1 + i*vbot1; ubot2 + i*vbot2];
    
    %interpolate data
    [interp_uvbot]=smart_interp(rawTime,raw_uvbot,interpTime,3); 

    startInd = find(~isnan(interp_uvbot),1,'first');
    startTime = julian2datenum(interpTime(startInd));
    interval = mean(diff(julian2datenum(interpTime)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_bottom.txt']);
    [nameu,fu,tidecon,xout]=t_tide(interp_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
    questBins1 = nc1{'u_1205'}.bins_questioned(:);
    if isempty(questBins1)
        lastGoodBin1 = size(nc1{'u_1205'}(:),2);
    else
        lastGoodBin1 = questBins1(1)-1;
    end
    usfc1 = nc1{'u_1205'}(:,lastGoodBin1);
    vsfc1 = nc1{'v_1206'}(:,lastGoodBin1);
    bads=find(usfc1+vsfc1 >= 1e35);
    usfc1(bads) = nan;
    vsfc1(bads) = nan;
    
    questBins2 = nc2{'u_1205'}.bins_questioned(:);
    if isempty(questBins2)
        lastGoodBin2 = size(nc2{'u_1205'}(:),2);
    else
        lastGoodBin2 = questBins2(1)-1;
    end
    usfc2 = nc2{'u_1205'}(:,lastGoodBin2);
    vsfc2 = nc2{'v_1206'}(:,lastGoodBin2);
    bads=find(usfc2+vsfc2 >= 1e35);
    usfc2(bads) = nan;
    vsfc2(bads) = nan;
    
    raw_uvsfc=[usfc1 + i*vsfc1; usfc2 + i*vsfc2];

    %interpolate data
    [interp_uvsfc]=smart_interp(rawTime,raw_uvsfc,interpTime,3); 

    startInd = find(~isnan(interp_uvsfc),1,'first');
    startTime = julian2datenum(interpTime(startInd));
    interval = mean(diff(julian2datenum(interpTime)))*24;
    ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_surface.txt']);
    [nameu,fu,tidecon,xout]=t_tide(interp_uvsfc,'start time',startTime,...
                                                'interval',interval,...
                                                'latitude',latitude,...
                                                'output',ttideOutFile);
                                   
    %%%%%%run t_tide on velocity data at 4 mbs and 6 mbs
    targetDepths=[4,6]; %mbs
    for b=1:size(targetDepths,2)
        [index_num,dist_to_hit] = value2Index(binDepths1,targetDepths(b),0.5);
        u1 = nc1{'u_1205'}(:,index_num);
        v1 = nc1{'v_1206'}(:,index_num);
        bads=find(u1+v1 >= 1e35);
        u1(bads) = nan;
        v1(bads) = nan;
        
        [index_num,dist_to_hit] = value2Index(binDepths2,targetDepths(b),0.5);
        u2 = nc2{'u_1205'}(:,index_num);
        v2 = nc2{'v_1206'}(:,index_num);
        bads=find(u2+v2 >= 1e35);
        u2(bads) = nan;
        v2(bads) = nan;
        
        raw_uv=[u1 + i*v1; u2 + i*v2];

        %interpolate data
        [interp_uv]=smart_interp(rawTime,raw_uv,interpTime,3); 
        
        startInd = find(~isnan(interp_uv),1,'first');
        startTime = julian2datenum(interpTime(startInd));
        interval = mean(diff(julian2datenum(interpTime)))*24;
        ttideOutFile = fullfile(plotDirectory,['ttide_curr_site',site,'_',num2str(targetDepths(b)),'mbs.txt']);
        [nameu,fu,tidecon,xout]=t_tide(interp_uv,'start time',startTime,...
                                                 'interval',interval,...
                                                 'latitude',latitude,...
                                                 'output',ttideOutFile);
    end

    nc1=close(nc1);
    nc2=close(nc2);
   
end
   
    
    
    
    