%Curr_Sp_Dir_Plotmaker
%
%Program which creates a plot of the percentage of all of the current data
%from a given site at a given depth which fall into various bins of current
%speed and direction.
%Requires singleJD.m, varify.m, binsort.m
%Modified from a program (sc_bins) received from Ben Guitterez
%
%Soupy Alexander
%12/5/2001

ncclear
figure

%Inputs to determine what to plot
station = input('Enter the station of interest: ');
hit_depth = input('Enter the depth (mbs) of interest (B = Bottom): ');
data_type = input('Enter sampling interval (h = hour averaged, l = low-passed): ');

%Differentation of hour averaged vs. low-pass filtered
if data_type == 'h';
    dataTag = '-a1h';
elseif data_type == 'l';
    dataTag = '-alp';
else
    error('Non-existent data type.  Be sure to use lower case letter and single quotes')
end

%Selection of data files
if station == 'A' & hit_depth ~= 'B';
    fileName = ['5951adc' dataTag '.nc'];
elseif station == 'A' & hit_depth == 'B';
    fileName = ['5952v' dataTag '_d1.nc'];
elseif station == 'B' & hit_depth ~= 'B';
    fileName = ['5971adc' dataTag '.nc'];
elseif station == 'B' & hit_depth == 'B';
    fileName = ['5972v' dataTag '_d1.nc'];
elseif station == 'C' & hit_depth ~= 'B';
    fileName = ['5991adc' dataTag '.nc'];
elseif station == 'C' & hit_depth == 'B';
    fileName = ['5992v' dataTag '.nc'];
elseif station == 'D' & hit_depth ~= 'B';
    fileName = ['6011adc' dataTag '.nc'];
elseif station == 'D' & hit_depth == 'B';
    fileName = ['6012v' dataTag '_d1.nc'];
elseif station == 'E' & hit_depth ~= 'B';
    fileName = ['6031adc' dataTag '.nc'];
elseif station == 'E' & hit_depth == 'B';
    fileName = ['6032v' dataTag '.nc'];
elseif station == 'F';
    fileName = ['6041adc' dataTag '.nc'];
else
    error('Non-existent station.  Be sure to use capital letter and single quotes')
end    

%Extract useful information
file = netCDF(fileName,'nowrite');
time = file{'time'};
time2 = file{'time2'};
jdTime = singleJD(time,time2);
u_data = file{'u_1205'}(:);
v_data = file{'v_1206'}(:);
depth_vector = file{'depth'}(:);

%Determine the index of the depth closest to the desired depth
if hit_depth ~= 'B';
    [depth_index,dist_to_depth] = value2Index(depth_vector,hit_depth);
        if dist_to_depth > 5;
            error('No depth data available within 5 meters of selected depth')
        end
elseif hit_depth == 'B';
    depth_index = 1;
end    
        
u_data_hit = u_data(:,depth_index);
v_data_hit = v_data(:,depth_index);
act_depth = depth_vector(depth_index);

u_data_hit = ridfill(u_data_hit);
v_data_hit = ridfill(v_data_hit);

%Calculate the current speed and direction
cspd = sqrt(u_data_hit.^2 + v_data_hit.^2);
cdir = atan2(v_data_hit,u_data_hit)*(180/pi);

%Convert to proper compass angles
for index = 1:length(cdir);
    if cdir(index) <= 90 ,
         cdir(index) = abs(cdir(index) - 90);
   	elseif  cdir(index) > 90,
         cdir(index) = 360 - (cdir(index) - 90);
    end
end

%Set up bins
xbin = [0 22.5 45 67.5 90 112.5 135 157.5 180 202.5 225 247.5 270 292.5 315 337.5 360 382.5];
ybin = [0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60];

d_xbin = (mean(diff(xbin))/2);
d_ybin = (mean(diff(ybin))/2);

hspotx = zeros(length(ybin)-1,length(xbin)-1);
hspoty = zeros(length(ybin)-1,length(xbin)-1);

for j = 1:length(ybin)-1,
   for i = 1:length(xbin)-1,
        hspotx(j,i) = xbin(i);
		hspoty(j,i) = ybin(j);
    end
end

z = cspd;
z = varify(z);

%Define bins, set plotting parameters, plot results
[a,h] = binsort(cdir,cspd,z,xbin,ybin);

[m,n] = size(hspotx);

a2 = a(2:m+1,2:n+1);
h2 = h(2:m+1,2:n+1);


hmax = max(max(h));
hmin = min(min(h));

caxis = ([0.01 hmax]);
axis([0 360 0 50])
cm = colormap;
cm(1,:) = [1 1 1];
colormap(cm)

axis([0 360 0 60])

H = surface(hspotx,hspoty,h2)

shading faceted
ax1 = gca;
set(ax1,'ydir','reverse')
set(ax1,'Box','on')
set(ax1,'XTick',xbin,'FontSize',[8])
set(ax1,'YTick',ybin,'FontSize',[8])


hpos = find(h2 > 0.03);
hfin = h2(hpos);
hxpos = hspotx(hpos) + (d_xbin);
hypos = hspoty(hpos) + (d_ybin);

for i = 1:length(hpos),
   hname = hfin(i);
   hval = sprintf('%0.2g',hname);
   
   if hname > 1.3 & hname < 3.5,
      text(hxpos(i),hypos(i),hfin(i)+1,hval,'color',[0 0 0],'HorizontalAlignment','center',...
      'VerticalAlignment','middle','FontSize',[7],'FontWeight','bold')
   else
	   text(hxpos(i),hypos(i),hfin(i)+1,hval,'color',[1 1 1],'HorizontalAlignment','center',...
      'VerticalAlignment','middle','FontSize',[7],'FontWeight','bold')
   end
end

title(['TITLE:  NOTE STATION F, BOTTOM IS ADCP! Depth = ' num2str(depth_vector(depth_index))])

cb = colorbar

xlabel('Current Direction (compass degrees)')
ylabel('Current Speed (cm/s)')
axes(cb)
ylabel('Percent of data')