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

data_scale = 1;
plot_scale = 3;
range_y = [-40 40];

%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_D = netcdf('6011adc-alp.nc','nowrite');
    ADCP_file_E = netcdf('6031adc-alp.nc','nowrite');

%Pull out data of interest
    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'}(:);

%Locate the index for each station of the depth closest to the hit depth
if hit_depth == 'B';
    indexD = 1;
    indexE = 1;
    dist_D = 0;
    dist_E = 0;
else    
    [indexD,dist_D] = value2Index(D_depth,hit_depth);
    [indexE,dist_E] = value2Index(E_depth,hit_depth);
end

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


%Convert the component velocities to current speed and direction
%        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);

%Eliminate bad data        

badsD = find(u_D_depth > 80 | u_D_depth < -80 | v_D_depth > 80 | v_D_depth < -80);
u_D_depth(badsD) = NaN;
v_D_depth(badsD) = NaN;
badsE = find(u_E_depth > 80 | u_E_depth < -80 | v_E_depth > 80 | v_E_depth < -80);
u_E_depth(badsE) = NaN;
v_E_depth(badsE) = NaN;

%Cspeed_D(badsD) = NaN;
%Cangle_D(badsD) = NaN;
%badsE = find(Cspeed_E(:) > 80 | Cspeed_E(:) < -80);
%Cspeed_E(badsE) = NaN;
%Cangle_E(badsE) = 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 = 1;

subplot(3,1,1)
if hit_depth == 'B';
    title('Low-passed Wind Stress and Current at Shelf Stations Bottom Depth');
else
    title(['Low-passed Wind Stress and Current at Shelf 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(3,1,2)
if dist_D > 5;
    axis([startjd endjd -40 40])
    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
            gregaxd(D_time,5)
            set(gca,'ytick',[-40 -20 0 20 40],'yticklabel',[-40 -20 0 20 40]);

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,range_y);

text(julian(2000,1,1,00)*plot_scale,20*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(3,1,3)
if dist_E > 5;
    axis([startjd endjd -40 40])
    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(D_time,5)
            set(gca,'ytick',[-40 -20 0 20 40],'yticklabel',[-40 -20 0 20 40]);

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,range_y);

text(julian(2000,1,1,00)*plot_scale,20*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