%Mean_Flow_Ell_BASS
%m-file to create whisker plots of mean flow from BASS data; two plots will
%be created, one for valley stations and one for shelf stations.  
%Also includes low-passed wind stress.
%Labels/Location of text in general was modified for appearance
%using Adobe Illustrator
%
%Soupy Alexander, 03/05/2002
%Requires NetCDF toolbox, "singleJD.m", "value2Index", "gregaxdNM.m", 
%"gregaxd.m", "gregaxd_fooled", and "NYCurr_Whisk.m"

data_type = input('Hourly averaged (h) or low-passed data (l)? ');

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

data_scale_wind = 2;
plot_scale_wind = 1;

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

%Input the Ambrose data in the first plot of both shelf and valley
figure(1)
subplot(4,1,1)
title(['Low-passed Wind Stress and Current Mean Flow from Shelf Station BASS Data']);
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');

figure(2)
subplot(4,1,1)
title(['Low-passed Wind Stress and Current Mean Flow Valley Station BASS Data']);
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');

shelf_stations = ['A'; 'D'; 'E';];
valley_stations = ['A'; 'B'; 'C';];

for index_plot = 1:2;
    figure(index_plot)
    if index_plot == 1;
        interest = shelf_stations;
            yrange = [-20 20];
            data_scale = 1;
            plot_scale = 50;
    elseif index_plot == 2;
        interest = valley_stations;
            yrange = [-20 20];
            data_scale = 1;
            plot_scale = 50;
    end
for index_stat = 1:3;
    subplot(4,1,index_stat+1)
    station = interest(index_stat);

if station == 'A';
    BASS_file = netcdf(['5952v-' data_tag '_d1.nc']);
elseif station == 'B';
    BASS_file = netcdf(['5972v-' data_tag '_d1.nc']);
elseif station == 'C';
    BASS_file = netcdf(['5992v-' data_tag '.nc']);  
elseif station == 'D';
    BASS_file = netcdf(['6012v-' data_tag '_d1.nc']); 
elseif station == 'E';
    BASS_file = netcdf(['6032v-' data_tag '.nc']);  
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
BASS_t = BASS_file{'time'}(:);
BASS_t2 = BASS_file{'time2'}(:);
BASS_time = singleJD(BASS_t,BASS_t2);
BASS_depth = BASS_file{'depth'}(:);
BASS_u_data = BASS_file{'u_1205'}(:);
BASS_v_data = BASS_file{'v_1206'}(:);

%Eliminate unreasonable speeds
    u_depths = ridfill(BASS_u_data);
    v_depths = ridfill(BASS_v_data);

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

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

[u_scale,t_ind] = ridnan(u_all_scale);
[v_scale,t_ind] = ridnan(v_all_scale);

tt1 = BASS_t(t_ind);
tt2 = BASS_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(BASS_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')
    ylabel('cm/s')    
    set(gca,'box','on')  
end
orient landscape
end