%ADCP_Whisk_Rot_Val_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, vectors were rotated to along and across axis
%
%Soupy Alexander, 11/18/2001
%Requires NetCDF toolbox, "singleJD.m", "value2Index", "gregaxdNM.m", 
%"gregaxd.m", "gregaxd_fooled", and "NYCurr_Whisk.m"
%
%Modified 04/09/02 to rotate along valley axis

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_depth1 = A_u_data(:,indexA);
v_A_depth1 = A_v_data(:,indexA);
end

u_B_depth1 = B_u_data(:,indexB);
v_B_depth1 = B_v_data(:,indexB);
u_C_depth1 = C_u_data(:,indexC);
v_C_depth1 = C_v_data(:,indexC);
u_F_depth1 = F_u_data(:,indexF);
v_F_depth1 = F_v_data(:,indexF);

%Eliminate bad data & rotate along axis
%Angles of rotation chosen by Brad Butman (looked at current and bathy data)
if isnan(indexA) == 0;
badsA = find(u_A_depth1 > 80 | u_A_depth1 < -80 | v_A_depth1 > 80 | v_A_depth1 < -80);
u_A_depth1(badsA) = NaN;
v_A_depth1(badsA) = NaN;
[v_A_depth,u_A_depth] = transformNY(u_A_depth1,v_A_depth1,0);
end
badsB = find(u_B_depth1 > 80 | u_B_depth1 < -80 | v_B_depth1 > 80 | v_B_depth1 < -80);
u_B_depth1(badsB) = NaN;
v_B_depth1(badsB) = NaN;
[v_B_depth,u_B_depth] = transformNY(u_B_depth1,v_B_depth1,15);
badsC = find(u_C_depth1 > 80 | u_C_depth1 < -80 | v_C_depth1 > 80 | v_C_depth1 < -80);
u_C_depth1(badsC) = NaN;
v_C_depth1(badsC) = NaN;
[v_C_depth,u_C_depth] = transformNY(u_C_depth1,v_C_depth1,65);
badsF = find(u_F_depth1 > 80 | u_F_depth1 < -80 | v_F_depth1 > 80 | v_F_depth1 < -80);
u_F_depth1(badsF) = NaN;
v_F_depth1(badsF) = NaN;
[v_F_depth,u_F_depth] = transformNY(u_F_depth1,v_F_depth1,80);



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


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


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


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