%ADCP_Whisk_Plotmaker
%m-file to create whisker plots of ADCP data from each station for the NY Bight
%Labels/Location of text in general was modified for appearance
%using Adobe Illustrator
%
%Soupy Alexander, 11/18/2001
%Requires NetCDF toolbox, "singleJD.m", "value2Index", "gregaxdNM.m", 
%"gregaxd.m", "gregaxd_fooled", and "NYCurr_Whisk.m"



%Determine the ADCP file corresponding to the station of interest
station = input('Enter the station of interest, use single quotes ');

%Determine the depths to plot for valley stations (A,B,C,F) vs. shelf stations
if station == 'B' | station == 'C' | station == 'F';
    hit_depths = [5 15 30 45];
    y_range = [-80 80];
    data_scale = 1;
    plot_scale = 12;
    data_scale_wind = 3;
    plot_scale_wind = 2;
else
    hit_depths = [5 15];
    y_range = [-40 40];
    data_scale = 1;
    plot_scale = 5;
    data_scale_wind = 2;
    plot_scale_wind = 1;
end    

%Pull the ADCP data
if station == 'A';
    ADCP_file = netcdf('5951adc-alp.nc','nowrite');
elseif station == 'B';
    ADCP_file = netcdf('5971adc-alp.nc','nowrite');
elseif station == 'C';
    ADCP_file = netcdf('5991adc-alp.nc','nowrite');
elseif station == 'D';
    ADCP_file = netcdf('6011adc-alp.nc','nowrite');
elseif station == 'E';
    ADCP_file = netcdf('6031adc-alp.nc','nowrite');
elseif station == 'F';
    ADCP_file = netcdf('6041adc-alp.nc','nowrite');
end

%Set the time limits for the plot
startjd = julian(1999,12,01,00);
endjd = julian(2000,4,20,00);

%Pull out data of interest
ADCP_t = ADCP_file{'time'}(:);
ADCP_t2 = ADCP_file{'time2'}(:);
ADCP_time = singleJD(ADCP_t,ADCP_t2);
ADCP_depth = ADCP_file{'depth'}(:);
ADCP_u_data = ADCP_file{'u_1205'}(:);
ADCP_v_data = ADCP_file{'v_1206'}(:);

%Add the bottom depth to the depths of interest
hit_depths(end+1) = max(ADCP_depth);

%Find the indices of the depths of interest
for index = 1:length(hit_depths);
    [index_hit(index),dist_to_hit(index)] = ...
        value2Index(ADCP_depth,hit_depths(index));
end

index_hit(find(isnan(index_hit)==1))=1;
dist_to_hit(find(isnan(index_hit)==1)) = 100;

%Select out the u and v's of the depths of interest
for index = 1:length(hit_depths);
    u_depths(:,index) = ADCP_u_data(:,index_hit(index));
    v_depths(:,index) = ADCP_v_data(:,index_hit(index));
end

%%Convert the u and v components to speed and direction
%for index = 1:length(hit_depths);
%        Cspeed(:,index) = sqrt(ADCP_u_data(:,index_hit(index)).^2 + ...
%            ADCP_v_data(:,index_hit(index)).^2);
%        Cangle(:,index) = atan2(ADCP_v_data(:,index_hit(index)),...
%            ADCP_u_data(:,index_hit(index)))*(180/pi);
%end

%Eliminate unreasonable speeds

for index = 1:length(hit_depths);
    bads = find(u_depths(:,index) > 80 | u_depths(:,index) < -80 | ...
        v_depths(:,index) > 80 | v_depths(:,index) < -80);
    u_depths(bads,index) = NaN;
    v_depths(bads,index) = NaN;
end

%for index = 1:length(hit_depths);    
%    bads = find(Cspeed(:,index) > 80 | Cspeed(:,index) < -80);
%    Cspeed(bads,index) = NaN;
%    Cangle(bads,index) = NaN;
%end

%Wind Stress from Ambrose Lighthouse
ambrose = netcdf('alsn6-a.nc','nowrite');
wu_ambrose = ambrose{'WU_422'}(:);
wv_ambrose = ambrose{'WV_423'}(:);
time_ambrose = ambrose{'time'}(:);
time2_ambrose = ambrose{'time2'}(:);

jdtime_ambrose = singleJD(time_ambrose,time2_ambrose);

%Create the plots
%Determine the number of subplots (one each for depth + 1 wind stress)
num_depths = length(hit_depths);
num_plots = num_depths + 1;

subplot(num_plots,1,1)
title(['Low-passed Wind Stress and Current from Station ' station]);
whiskWindNY(jdtime_ambrose,wu_ambrose,wv_ambrose,data_scale_wind,plot_scale_wind);
ylabel('dynes/cm^2')
text(julian(2000,1,1,00)*plot_scale,0,'Estimated 10 m Wind Stress, Ambrose Light');

for index = 2:num_plots;
    subplot(num_plots,1,index)
        if dist_to_hit(index-1) > 5;
    axis([startjd endjd -80 80])
            if hit_depths(index-1) > max(ADCP_depth);
    text(julian(2000,1,1,00),0,['Depth ' hit_depths(index-1) ' m is deeper than this station']);
            else
    text(julian(2000,1,1,00),0,'Data unavailable from this depth');
            end
    gregaxd(ADCP_time,5)        
        else
%%I was originally using sticksafe to plot velocity data; I am no longer
%convinced at this program preserves angles and magnitudes so I have
%replaced it with a site/data specific whisker plotting program.
%    dummy_y = zeros(size(ADCP_time));
%    sticksafe(ADCP_time,dummy_y,Cspeed(:,index-1),...
%        Cangle(:,index-1));
%    axis([startjd endjd -80 80])
NYCurr_Whisk(u_depths(:,index-1),v_depths(:,index-1),ADCP_time,data_scale,plot_scale,y_range);

    if index == num_plots;
        text(julian(2000,1,1,00)*plot_scale,40,['Bottom, Depth = '...
                num2str(ADCP_depth(index_hit(index-1))) ' m']);
    else
    text(julian(2000,1,1,00)*plot_scale,40,['Depth = ' ...
            num2str(ADCP_depth(index_hit(index-1))) ' m'])
    end
        end
    ylabel('cm/s')    
%    gregaxd(ADCP_time,5)
%    set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')  
end

orient landscape