%ADCP_Whisk_Valley_Plotmaker
%m-file to create whisker plots of ADCP data from valley stations for the NY Bight
%Labels/Location of text in general was modified for appearance (Stations A-C, F)
%using Adobe Illustrator
%
%Soupy Alexander, 11/18/2001
%Requires NetCDF toolbox, "singleJD.m", "value2Index", "gregaxdNM.m", 
%"gregaxd.m", "gregaxd_fooled", and "NYCurr_Whisk.m"

y_range = [-80 80];
data_scale = 1;
plot_scale = 10;

%User input as to what depth to plot from each station
hit_depth = input('Input depth of interest, "B" for bottom ');

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

%Identify the files of interest
    ADCP_file_A = netcdf('5951adc-alp.nc','nowrite');
    ADCP_file_B = netcdf('5971adc-alp.nc','nowrite');
    ADCP_file_C = netcdf('5991adc-alp.nc','nowrite');
    ADCP_file_F = netcdf('6041adc-alp.nc','nowrite');

%Pull out data of interest
    A_t = ADCP_file_A{'time'}(:);
    A_t2 = ADCP_file_A{'time2'}(:);
    A_time = singleJD(A_t,A_t2);
    A_depth = ADCP_file_A{'depth'}(:);
    A_u_data = ADCP_file_A{'u_1205'}(:);
    A_v_data = ADCP_file_A{'v_1206'}(:);
    
    B_t = ADCP_file_B{'time'}(:);
    B_t2 = ADCP_file_B{'time2'}(:);
    B_time = singleJD(B_t,B_t2);
    B_depth = ADCP_file_B{'depth'}(:);
    B_u_data = ADCP_file_B{'u_1205'}(:);
    B_v_data = ADCP_file_B{'v_1206'}(:);
    
    C_t = ADCP_file_C{'time'}(:);
    C_t2 = ADCP_file_C{'time2'}(:);
    C_time = singleJD(C_t,C_t2);
    C_depth = ADCP_file_C{'depth'}(:);
    C_u_data = ADCP_file_C{'u_1205'}(:);
    C_v_data = ADCP_file_C{'v_1206'}(:);
    
    F_t = ADCP_file_F{'time'}(:);
    F_t2 = ADCP_file_F{'time2'}(:);
    F_time = singleJD(F_t,F_t2);
    F_depth = ADCP_file_F{'depth'}(:);
    F_u_data = ADCP_file_F{'u_1205'}(:);
    F_v_data = ADCP_file_F{'v_1206'}(:);
    
%Locate the index for each station of the depth closest to the hit depth
if hit_depth == 'B';
    indexA = 1;
    indexB = 1;
    indexC = 1;
    indexF = 1;
    dist_A = 0;
    dist_B = 0;
    dist_C = 0;
    dist_F = 0;
else    
    [indexA,dist_A] = value2Index(A_depth,hit_depth);
    [indexB,dist_B] = value2Index(B_depth,hit_depth);
    [indexC,dist_C] = value2Index(C_depth,hit_depth);
    [indexF,dist_F] = value2Index(F_depth,hit_depth);
end

if isnan(indexA) == 0;
u_A_depth = A_u_data(:,indexA);
v_A_depth = A_v_data(:,indexA);
end

u_B_depth = B_u_data(:,indexB);
v_B_depth = B_v_data(:,indexB);
u_C_depth = C_u_data(:,indexC);
v_C_depth = C_v_data(:,indexC);
u_F_depth = F_u_data(:,indexF);
v_F_depth = F_v_data(:,indexF);

%Convert the component velocities to current speed and direction
%        Cspeed_A = sqrt(A_u_data(:,indexA).^2 + ...
%            A_v_data(:,indexA).^2);
%        Cangle_A = atan2(A_v_data(:,indexA),...
%            A_u_data(:,indexA))*(180/pi);
%        Cspeed_B = sqrt(B_u_data(:,indexB).^2 + ...
%            B_v_data(:,indexB).^2);
%        Cangle_B = atan2(B_v_data(:,indexB),...
%            B_u_data(:,indexB))*(180/pi);
%        Cspeed_C = sqrt(C_u_data(:,indexC).^2 + ...
%            C_v_data(:,indexC).^2);
%        Cangle_C = atan2(C_v_data(:,indexC),...
%            C_u_data(:,indexC))*(180/pi);
%        Cspeed_F = sqrt(F_u_data(:,indexF).^2 + ...
%            F_v_data(:,indexF).^2);
%        Cangle_F = atan2(F_v_data(:,indexF),...
%            F_u_data(:,indexF))*(180/pi);

%Eliminate bad data        
if isnan(indexA) == 0;
badsA = find(u_A_depth > 80 | u_A_depth < -80 | v_A_depth > 80 | v_A_depth < -80);
u_A_depth(badsA) = NaN;
v_A_depth(badsA) = NaN;
end
badsB = find(u_B_depth > 80 | u_B_depth < -80 | v_B_depth > 80 | v_B_depth < -80);
u_B_depth(badsB) = NaN;
v_B_depth(badsB) = NaN;
badsC = find(u_C_depth > 80 | u_C_depth < -80 | v_C_depth > 80 | v_C_depth < -80);
u_C_depth(badsC) = NaN;
v_C_depth(badsC) = NaN;
badsF = find(u_F_depth > 80 | u_F_depth < -80 | v_F_depth > 80 | v_F_depth < -80);
u_F_depth(badsF) = NaN;
v_F_depth(badsF) = NaN;

%badsA = find(Cspeed_A(:) > 80 | Cspeed_A(:) < -80);
%Cspeed_A(badsA) = NaN;
%Cangle_A(badsA) = NaN;
%badsB = find(Cspeed_B(:) > 80 | Cspeed_B(:) < -80);
%Cspeed_B(badsB) = NaN;
%Cangle_B(badsB) = NaN;
%badsC = find(Cspeed_C(:) > 80 | Cspeed_C(:) < -80);
%Cspeed_C(badsC) = NaN;
%Cangle_C(badsC) = NaN;
%badsF = find(Cspeed_F(:) > 80 | Cspeed_F(:) < -80);
%Cspeed_F(badsF) = NaN;
%Cangle_F(badsF) = NaN;

%Wind Stress

%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);

    data_scale_wind = 3;
    plot_scale_wind = 2;
    
subplot(5,1,1)
if hit_depth == 'B';
    title('Low-passed Wind Stress and Current at Valley Stations Bottom Depth');
else
    title(['Low-passed Wind Stress and Current at Valley Stations Depth ~ ' num2str(hit_depth) ' m'])
end
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_wind,0,'Estimated 10 m Wind Stress, Ambrose Light');


%Whisker plot the ADCP data from each station for the desired depths
subplot(5,1,2)
if isnan(dist_A) == 1;
    axis([startjd endjd -80 80])
    if hit_depth > max(A_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station A');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station A for This Depth']); 
    end
            gregaxdNM(A_time,5)
            set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
else
%    dummy_y = zeros(size(A_time));
%    sticksafe(A_time,dummy_y,Cspeed_A(:),Cangle_A(:))
%    axis([startjd endjd -80 80])
NYCurr_Whisk(u_A_depth(:),v_A_depth(:),A_time,data_scale,plot_scale,y_range);

text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station A, Depth = ' ...
            num2str(A_depth(indexA)) ' m'])
end
    ylabel('cm/s')    

 %   set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')


subplot(5,1,3)
if dist_B > 5;
    axis([startjd endjd -80 80])
    if hit_depth > max(B_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station B');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station B for This Depth']);
    end
            gregaxdNM(A_time,5)
            set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);

else
%    dummy_y = zeros(size(B_time));
%    sticksafe(B_time,dummy_y,Cspeed_B(:),Cangle_B(:))
%    axis([startjd endjd -80 80])

NYCurr_Whisk(u_B_depth(:),v_B_depth(:),B_time,data_scale,plot_scale,y_range);

    text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station B, Depth = ' ...
            num2str(B_depth(indexB)) ' m'])
end
    ylabel('cm/s')    
%    set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')

subplot(5,1,4)
if dist_C > 5;
    axis([startjd endjd -80 80])
    if hit_depth > max(C_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station C');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station C for This Depth']);
    end
            gregaxdNM(A_time,5)
            set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);

else
%    dummy_y = zeros(size(C_time));
%    sticksafe(C_time,dummy_y,Cspeed_C(:),Cangle_C(:))
%    axis([startjd endjd -80 80])

NYCurr_Whisk(u_C_depth(:),v_C_depth(:),C_time,data_scale,plot_scale,y_range);

text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station C, Depth = ' ...
            num2str(C_depth(indexC)) ' m'])
end
    ylabel('cm/s')    
%    set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')

subplot(5,1,5)
if dist_F > 5;
    axis([startjd endjd -80 80])
    if hit_depth > max(F_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station F');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station F for This Depth']);
    end
            gregaxdNM(A_time,5)
            set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);

else
%    dummy_y = zeros(size(F_time));
%    sticksafe(F_time,dummy_y,Cspeed_F(:),Cangle_F(:))
%    axis([startjd endjd -80 80])

NYCurr_Whisk(u_F_depth(:),v_F_depth(:),F_time,data_scale,plot_scale,y_range);

text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station F, Depth = ' ...
            num2str(F_depth(indexF)) ' m'])
end
    ylabel('cm/s')    
%    set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')    
     
orient landscape