%ADCP_Whisk_All_Plotmaker
%m-file to create whisker plots of ADCP data from all stations 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"
%
%Commented out parts relate to the original program, which ran using a plotting
%tool called sticksafe.m which I am no longer sure preserves angles/lengths.

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

%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_D = netcdf('6011adc-alp.nc','nowrite');
    ADCP_file_E = netcdf('6031adc-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'}(:);
    
    D_t = ADCP_file_D{'time'}(:);
    D_t2 = ADCP_file_D{'time2'}(:);
    D_time = singleJD(D_t,D_t2);
    D_depth = ADCP_file_D{'depth'}(:);
    D_u_data = ADCP_file_D{'u_1205'}(:);
    D_v_data = ADCP_file_D{'v_1206'}(:);
    
    E_t = ADCP_file_E{'time'}(:);
    E_t2 = ADCP_file_E{'time2'}(:);
    E_time = singleJD(E_t,E_t2);
    E_depth = ADCP_file_E{'depth'}(:);
    E_u_data = ADCP_file_E{'u_1205'}(:);
    E_v_data = ADCP_file_E{'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;
    indexD = 1;
    indexE = 1;
    indexF = 1;
    dist_A = 0;
    dist_B = 0;
    dist_C = 0;
    dist_D = 0;
    dist_E = 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);
    [indexD,dist_D] = value2Index(D_depth,hit_depth);
    [indexE,dist_E] = value2Index(E_depth,hit_depth);
    [indexF,dist_F] = value2Index(F_depth,hit_depth);
end

u_A_depth = A_u_data(:,indexA);
v_A_depth = A_v_data(:,indexA);
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_D_depth = D_u_data(:,indexD);
v_D_depth = D_v_data(:,indexD);
u_E_depth = E_u_data(:,indexE);
v_E_depth = E_v_data(:,indexE);
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_D = sqrt(D_u_data(:,indexD).^2 + ...
%            D_v_data(:,indexD).^2);
%        Cangle_D = atan2(D_v_data(:,indexD),...
%            D_u_data(:,indexD))*(180/pi);
%        Cspeed_E = sqrt(E_u_data(:,indexE).^2 + ...
%            E_v_data(:,indexE).^2);
%        Cangle_E = atan2(E_v_data(:,indexE),...
%            E_u_data(:,indexE))*(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        

badsA = find(u_A_depth > 1000 | u_A_depth < -1000 | v_A_depth > 1000 | v_A_depth < -1000);
u_A_depth(badsA) = NaN;
v_A_depth(badsA) = NaN;
badsB = find(u_B_depth > 1000 | u_B_depth < -1000 | v_B_depth > 1000 | v_B_depth < -1000);
u_B_depth(badsB) = NaN;
v_B_depth(badsB) = NaN;
badsC = find(u_C_depth > 1000 | u_C_depth < -1000 | v_C_depth > 1000 | v_C_depth < -1000);
u_C_depth(badsC) = NaN;
v_C_depth(badsC) = NaN;
badsD = find(u_D_depth > 1000 | u_D_depth < -1000 | v_D_depth > 1000 | v_D_depth < -1000);
u_D_depth(badsD) = NaN;
v_D_depth(badsD) = NaN;
badsE = find(u_E_depth > 1000 | u_E_depth < -1000 | v_E_depth > 1000 | v_E_depth < -1000);
u_E_depth(badsE) = NaN;
v_E_depth(badsE) = NaN;
badsF = find(u_F_depth > 1000 | u_F_depth < -1000 | v_F_depth > 1000 | v_F_depth < -1000);
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;
%badsD = find(Cspeed_D(:) > 80 | Cspeed_D(:) < -80);
%Cspeed_D(badsD) = NaN;
%Cangle_D(badsD) = NaN;
%badsE = find(Cspeed_E(:) > 80 | Cspeed_E(:) < -80);
%Cspeed_E(badsE) = NaN;
%Cangle_E(badsE) = NaN;
%badsF = find(Cspeed_F(:) > 80 | Cspeed_F(:) < -80);
%Cspeed_F(badsF) = NaN;
%Cangle_F(badsF) = NaN;

%Whisker plot the ADCP data from each station for the desired depths
subplot(6,1,1)
if dist_A > 5;
    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')
if hit_depth == 'B';
    title('Low-passed Current at All Stations from the Bottom');
else
    title(['Low-passed Current at All Stations - Depth ~ ' num2str(hit_depth) ' m'])
end

subplot(6,1,2)
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(6,1,3)
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(6,1,4)
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')    
    
subplot(6,1,5)
if dist_D > 5;
    axis([startjd endjd -80 80])
    if hit_depth > max(D_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station D');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station D 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(D_time));
%    sticksafe(D_time,dummy_y,Cspeed_D(:),Cangle_D(:))
%    axis([startjd endjd -80 80])

NYCurr_Whisk(u_D_depth(:),v_D_depth(:),D_time,data_scale,plot_scale,y_range);

text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station D, Depth = ' ...
            num2str(D_depth(indexD)) ' m'])
end
    ylabel('cm/s')    
%    set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);
    set(gca,'box','on')
    
subplot(6,1,6)
if dist_E > 5;
    axis([startjd endjd -80 80])
    if hit_depth > max(E_depth);
        text(julian(2000,1,1,00),0,'Depth is deeper than Station E');
    else
        text(julian(2000,1,1,00),0,['ADCP Data Unavailable at Station E for This Depth']);
    end
            gregaxd(A_time,5)
            set(gca,'ytick',[-80 -40 0 40 80],'yticklabel',[-80 -40 0 40 80]);

else
%    dummy_y = zeros(size(E_time));
%    sticksafe(E_time,dummy_y,Cspeed_E(:),Cangle_E(:))
%    axis([startjd endjd -80 80])

NYCurr_Whisk(u_E_depth(:),v_E_depth(:),E_time,data_scale,plot_scale,y_range);

text(julian(2000,1,1,00)*plot_scale,40*data_scale,['Station E, Depth = ' ...
            num2str(E_depth(indexE)) ' 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