%Mean_Flow_Ellipse_Plot
%m-file to create whisker plots of mean flow from ADCP data at each station 
%for the NY Bight.  Also includes low-passed wind stress.
%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 ');
data_type = input('Hourly-averaged (h) or low-passed (l): ');

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

if data_type == 'h';
    data_tag = 'a1h.nc';
elseif data_type == 'l';
    data_tag = 'alp.nc';
end

%Pull the ADCP data
if station == 'A';
    ADCP_file = netcdf(['5951adc-' data_tag],'nowrite');
elseif station == 'B';
    ADCP_file = netcdf(['5971adc-' data_tag],'nowrite');
elseif station == 'C';
    ADCP_file = netcdf(['5991adc-' data_tag],'nowrite');
elseif station == 'D';
    ADCP_file = netcdf(['6011adc-' data_tag],'nowrite');
elseif station == 'E';
    ADCP_file = netcdf(['6031adc-' data_tag],'nowrite');
elseif station == 'F';
    ADCP_file = netcdf(['6041adc-' data_tag],'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

%Select out the u and v's of the depths of interest
for index = 1:length(hit_depths);
    if isnan(index_hit(index))==1;
        u_depths(:,index) = 0;
        v_depths(:,index) = 0;
    else
    u_depths(:,index) = ADCP_u_data(:,index_hit(index));
    v_depths(:,index) = ADCP_v_data(:,index_hit(index));
end
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);
    u_depths(:,index) = ridfill(u_depths(:,index));
    v_depths(:,index) = ridfill(v_depths(:,index));
end

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

theta=2*pi*(1:64)/64;
theta=[0 theta];

u_all_scale = data_scale*u_depths;
v_all_scale = data_scale*v_depths;


%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 later (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 Mean Flow 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_wind,0,'Estimated 10 m Wind Stress, Ambrose Light');


for index_plot = 2:num_plots;
    subplot(num_plots,1,index_plot)
        if isnan(index_hit(index_plot-1))==1;
axis([0*plot_scale 6*plot_scale round(min(yrange)*data_scale) round(max(yrange)*data_scale)])
axis equal
box('on')
set(gca,'ytick',floor(data_scale.*[min(yrange):10:max(yrange)]));
ylabel('cm/s')
set(l,'color','k')
set(gca,'xtick',plot_scale.*[1:5]);
set(gca,'xticklabel',['Dec';'Jan';'Feb';'Mar';'Apr';]);
            if hit_depths(index_plot-1) > max(ADCP_depth);
    text(1*plot_scale,0,['Depth ' hit_depths(index_plot-1) ' m is deeper than this station']);
            else
    text(1*plot_scale,0,'Data unavailable from this depth');
            end
    gregaxd(ADCP_time,5)        
        else
u_hit = u_all_scale(:,index_plot-1);
v_hit = v_all_scale(:,index_plot-1);
[u_scale,t_ind] = ridnan(u_hit);
[v_scale,t_ind] = ridnan(v_hit);

tt1 = ADCP_t(t_ind);
tt2 = ADCP_t2(t_ind);

dnum = ep_datenum([tt1 tt2]);
[Y,M,DD,H,MI,S] = datevec(dnum);

month_list = [12 1 2 3 4];

for index = 1:5;
    this_month = find(M == month_list(index));
diffs = mean(diff(ADCP_time));    
        if (diffs > 0.1 & length(this_month) >= 60) | ...
                (diffs < 0.1 & length(this_month) >= 360);
    mean_u(index) = nanmean(u_scale(this_month));
    mean_v(index) = nanmean(v_scale(this_month));
    [majax(index),majaz(index),minax(index),minaz(index),ellip(index)] = ...
        pcaben(u_scale(this_month),v_scale(this_month));
    [x1(index),y1(index)]=cmgspd2uv(majax(index),majaz(index));
    [x2(index),y2(index)]=cmgspd2uv(minax(index),minaz(index));
    xx(:,index)=(majax(index)*cos(theta))';
    yy(:,index)=(minax(index)*sin(theta))';
    angle(index)=-majaz(index)*pi/180+pi/2;
    xxx(:,index)=xx(:,index).*cos(angle(index))-yy(:,index).*sin(angle(index));
    yyy(:,index)=xx(:,index).*sin(angle(index))+yy(:,index).*cos(angle(index));
        else
    month_list(index) = NaN;        
        end
end    

hold on
keep_index = [1:5];
bads = find(isnan(month_list)==1);
month_list(bads) = [];
keep_index(bads) = [];
clist = ['r';'g';'b';'c';'m'];

for index = keep_index;

plot_loc = plot_scale*index;  
    
%This bit could add the axis lines to the ellipse circles    
    
%hh1=plot([mean_u(index)+index+x1(index) mean_u(index)+index-x1(index)],...
%    [mean_v(index)+y1(index) mean_v(index)-y1(index)],'-.r',...
%    [mean_u(index)+index+x2(index) mean_u(index)+index-x2(index)],...
%    [mean_v(index)+y2(index) mean_v(index)-y2(index)],'-.r','linewidth',0.5);

hh1(3)=plot(mean_u(index)+plot_loc+xxx(:,index),mean_v(index)+yyy(:,index),clist(index),'linewidth',1);
hh1(4)=plot([plot_loc; mean_u(index)+plot_loc],[0; mean_v(index)],clist(index),'linewidth',2);
end

axis equal
axis([0*plot_scale 6*plot_scale round(min(yrange)*data_scale) round(max(yrange)*data_scale)])
set(gca,'ytick',floor(data_scale.*[min(yrange):10:max(yrange)]));
ylabel('cm/s')
l = line([0 6*plot_scale],[0 0]);
set(l,'color','k')
set(gca,'xtick',plot_scale.*[1:5]);
set(gca,'xticklabel',['Dec';'Jan';'Feb';'Mar';'Apr';]);
box('on')

    if index_plot == num_plots;
        text(0.25*plot_scale,-25,['Bottom, Depth = '...
                num2str(ADCP_depth(index_hit(index_plot-1))) ' m']);
    else
    text(0.25*plot_scale,-25,['Depth = ' ...
            num2str(ADCP_depth(index_hit(index_plot-1))) ' m'])
    end
        end
    ylabel('cm/s')    
    set(gca,'box','on')  
end

orient landscape