%Curr_Scat_Plot
%
%Program to do a scatter plot of eastward component vs. northward component
%of current from the New York Bight BASS/MAVS data.  Superimposes the ellipse
% with the principal axes.  Needs to be given a specific
%station, and creates a four-paneled plot with hourly-averaged and low-passed
%data from 5 mbs and the bottom-most depth bin.
%
%Soupy Alexander, 1/14/02

%Select a station

station = input('Select a station. ');

%Pull out the current data for the proper station
if station == 'A';
    file_name_h = '5952v-a1h_d1.nc';
    file_name_lp = '5952v-alp_d1.nc';
    file_named2_h = '5952v-a1h_d2.nc';
    file_named2_lp = '5952v-alp_d2.nc';    
elseif station == 'B';
    file_name_h = '5972v-a1h_d1.nc';
    file_name_lp = '5972v-alp_d1.nc';
    file_named2_h = '5972v-a1h_d2.nc';
    file_named2_lp = '5972v-alp_d2.nc';   
elseif station == 'C';
    file_name_h = '5992v-a1h.nc';
    file_name_lp = '5992v-alp.nc';  
elseif station == 'D';
    file_name_h = '6012v-a1h_d1.nc';
    file_name_lp = '6012v-alp_d1.nc';
    file_named2_h = '6012v-a1h_d2.nc';
    file_named2_lp = '6012v-alp_d2.nc';  
elseif station == 'E';
    file_name_h = '6032v-a1h.nc';
    file_name_lp = '6032v-alp.nc';  
elseif station == 'F';
    error('Station F has no VMCM/BASS data!!');
end

file_h = netcdf(file_name_h,'nowrite');
file_lp = netcdf(file_name_lp,'nowrite');

%Take care of the hourly averaged data:  pull out current from proper bin
t_h = file_h{'time'}(:);
t2_h = file_h{'time2'}(:);
time_h_orig = singleJD(t_h,t2_h);

depth_h = file_h{'depth'};

u_h_orig = file_h{'u_1205'}(:);
v_h_orig = file_h{'v_1206'}(:);

[u_h,nh] = ridnan(u_h_orig);
[v_h,nh] = ridnan(v_h_orig);
time_h = time_h_orig(nh);

u_h = ridfill(u_h);
v_h = ridfill(v_h);


%Calculate the principal ellipses & mean currents.
theta=2*pi*(1:64)/64;
theta=[0 theta];

[majax_h,majaz_h,minax_h,minaz_h,ellip_h] = pcaben(u_h,v_h);
[x1_h,y1_h]=cmgspd2uv(majax_h,majaz_h);
[x2_h,y2_h]=cmgspd2uv(minax_h,minaz_h);
xx_h=(majax_h*cos(theta))';
yy_h=(minax_h*sin(theta))';
angle_h=-majaz_h*pi/180+pi/2;
xxx_h=xx_h(:).*cos(angle_h)-yy_h.*sin(angle_h);
yyy_h=xx_h(:).*sin(angle_h)+yy_h.*cos(angle_h);
meanu_h = nanmean(u_h);
meanv_h = nanmean(v_h);

[meansp_h,meandir_h] = UVtoSpDir(meanu_h,meanv_h);


%Now take care of the low-passed data
t_lp = file_lp{'time'}(:);
t2_lp = file_lp{'time2'}(:);
time_lp_orig = singleJD(t_lp,t2_lp);

u_lp_orig = file_lp{'u_1205'}(:);
v_lp_orig = file_lp{'v_1206'}(:);

[u_lp,nlp] = ridnan(u_lp_orig);
[v_lp,nlp] = ridnan(v_lp_orig);

[u_lp] = ridfill(u_lp);
[v_lp] = ridfill(v_lp);
time_lp = time_lp_orig(nlp);

ud2_lp_orig = file_lp{'u_1205'}(:);
vd2_lp_orig = file_lp{'v_1206'}(:);

[ud2_lp,nlp] = ridnan(ud2_lp_orig);
[vd2_lp,nlp] = ridnan(vd2_lp_orig);

ud2_lp = ridfill_nan(ud2_lp);
vd2_lp = ridfill_nan(vd2_lp);

%Calculate the principal ellipses and means.
theta=2*pi*(1:64)/64;
theta=[0 theta];

[majax_lp,majaz_lp,minax_lp,minaz_lp,ellip_lp] = pcaben(u_lp,v_lp);
[x1_lp,y1_lp]=cmgspd2uv(majax_lp,majaz_lp);
[x2_lp,y2_lp]=cmgspd2uv(minax_lp,minaz_lp);
xx_lp=(majax_lp*cos(theta))';
yy_lp=(minax_lp*sin(theta))';
angle_lp=-majaz_lp*pi/180+pi/2;
xxx_lp=xx_lp(:).*cos(angle_lp)-yy_lp.*sin(angle_lp);
yyy_lp=xx_lp(:).*sin(angle_lp)+yy_lp.*cos(angle_lp);

meanu_lp = nanmean(u_lp);
meanv_lp = nanmean(v_lp);

[meansp_lp,meandir_lp] = UVtoSpDir(meanu_lp,meanv_lp);


if station == 'A' | station == 'B' | station == 'D';
    filed2_h = netcdf(file_named2_h,'nowrite');
    filed2_lp = netcdf(file_named2_lp,'nowrite');
    
%Take care of the hourly averaged data:  pull out current from proper bin
    td2_h = filed2_h{'time'}(:);
    t2d2_h = filed2_h{'time2'}(:);
    timed2_h_orig = singleJD(td2_h,t2d2_h);
    
    depth_d2 = filed2_h{'depth'}(:);
    
    ud2_h_orig = filed2_h{'u_1205'}(:);
    vd2_h_orig = filed2_h{'v_1206'}(:);

    [ud2_h,nh_d] = ridnan(ud2_h_orig);
    [vd2_h,nh_d] = ridnan(vd2_h_orig);
    time_h_d2 = timed2_h_orig(nh_d);

    ud2_h = ridfill(ud2_h);
    vd2_h = ridfill(vd2_h);
   
    
    [majaxd2_h,majazd2_h,minaxd2_h,minazd2_h,ellipd2_h] = pcaben(ud2_h,vd2_h);
    [x1d2_h,y1d2_h]=cmgspd2uv(majaxd2_h,majazd2_h);
    [x2d2_h,y2d2_h]=cmgspd2uv(minaxd2_h,minazd2_h);
    xxd2_h=(majaxd2_h*cos(theta))';
    yyd2_h=(minaxd2_h*sin(theta))';
    angled2_h=-majazd2_h*pi/180+pi/2;
    xxxd2_h=xxd2_h(:).*cos(angled2_h)-yyd2_h.*sin(angled2_h);
    yyyd2_h=xxd2_h(:).*sin(angled2_h)+yyd2_h.*cos(angled2_h);
    
    meanud2_h = nanmean(ud2_h);
    meanvd2_h = nanmean(vd2_h);
    
    [meanspd2_h,meandird2_h] = UVtoSpDir(meanud2_h,meanvd2_h);

    
%Now take care of the low-passed data
    td2_lp = filed2_lp{'time'}(:);
    t2d2_lp = filed2_lp{'time2'}(:);
    timed2_lp_orig = singleJD(td2_lp,t2d2_lp);

    ud2_lp_orig = filed2_lp{'u_1205'}(:);
    vd2_lp_orig = filed2_lp{'v_1206'}(:);

    [ud2_lp,nlp] = ridnan(ud2_lp_orig);
    [vd2_lp,nlp] = ridnan(vd2_lp_orig);

    [ud2_lp] = ridfill(ud2_lp);
    [vd2_lp] = ridfill(vd2_lp);
    timed2_lp = timed2_lp_orig(nlp);

[majaxd2_lp,majazd2_lp,minaxd2_lp,minazd2_lp,ellipd2_lp] = pcaben(ud2_lp,vd2_lp);
[x1d2_lp,y1d2_lp]=cmgspd2uv(majaxd2_lp,majazd2_lp);
[x2d2_lp,y2d2_lp]=cmgspd2uv(minaxd2_lp,minazd2_lp);
xxd2_lp=(majaxd2_lp*cos(theta))';
yyd2_lp=(minaxd2_lp*sin(theta))';
angled2_lp=-majazd2_lp*pi/180+pi/2;
xxxd2_lp=xxd2_lp(:).*cos(angled2_lp)-yyd2_lp.*sin(angled2_lp);
yyyd2_lp=xxd2_lp(:).*sin(angled2_lp)+yyd2_lp.*cos(angled2_lp);

meanud2_lp = nanmean(ud2_lp);
meanvd2_lp = nanmean(vd2_lp);

[meanspd2_lp,meandird2_lp] = UVtoSpDir(meanud2_lp,meanvd2_lp);
 
end

%Now plot the data
%Hourly Averaged, Depth 1
subplot(2,2,1)
plot(u_h,v_h,'.')
hold on
plot([0 meanu_h],[0 meanv_h],'r','linewidth',2);
plot(xxx_h+meanu_h,yyy_h+meanv_h,'r','linewidth',2);
axis equal
axis tight
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
text(-45,-23,['Mean U = ' num2str(meanu_h,'%0.2f') ...
        '; Mean V = ' num2str(meanv_h,'%0.2f')]);
text(-45,-30,['Vector speed = ' num2str(meansp_h,'%0.1f') ... 
        '; Direction = ' num2str(angleXtoCompass(meandir_h),'%0.0f') '^o']);
text(-45,-38,['Major Axis = ' num2str(majax_h,'%0.1f') ...
        '; Minor Axis = ' num2str(minax_h,'%0.1f')]);
text(-45,-45,['Prin. Orient. = ' num2str(majaz_h,'%0.0f') '^o']);
l = line([-50 50], [0 0]);
l2 = line([0 0], [-50 50]);
set(l,'color','k', 'linestyle','-.');
set(l2,'color','k', 'linestyle','-.');
set(gca,'xtick',labels);
set(gca,'ytick',labels);
set(gca,'xticklabel','');
set(gca,'yticklabel','');

%Low-Passed, Depth 1
subplot(2,2,2)
plot(u_lp,v_lp,'.')
hold on
plot([0 meanu_lp],[0 meanv_lp],'r','linewidth',2);
plot(xxx_lp+meanu_lp,yyy_lp+meanv_lp,'r','linewidth',2);
axis equal
axis tight
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
text(-45,-23,['Mean U = ' num2str(meanu_lp,'%0.2f') ...
        '; Mean V = ' num2str(meanv_lp,'%0.2f')]);
text(-45,-30,['Vector speed = ' num2str(meansp_lp,'%0.1f') ... 
        '; Direction = ' num2str(angleXtoCompass(meandir_lp),'%0.0f') '^o']);
text(-45,-38,['Major Axis = ' num2str(majax_lp,'%0.1f') ...
        '; Minor Axis = ' num2str(minax_lp,'%0.1f')]);
text(-45,-45,['Prin. Orient. = ' num2str(majaz_lp,'%0.0f') '^o']);
l = line([-50 50], [0 0]);
l2 = line([0 0], [-50 50]);
set(l,'color','k', 'linestyle','-.');
set(l2,'color','k', 'linestyle','-.');
set(gca,'xtick',labels);
set(gca,'ytick',labels);
set(gca,'xticklabel','');
set(gca,'yticklabel','');


if station == 'A' | station == 'B' | station == 'D';
%Hourly Averaged, Depth 2
subplot(2,2,3)
plot(ud2_h,vd2_h,'.')
hold on
plot([0 meanud2_h],[0 meanvd2_h],'r','linewidth',2);
plot(xxxd2_h+meanud2_h,yyyd2_h+meanvd2_h,'r','linewidth',2);
axis equal
axis tight
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
text(-45,-23,['Mean U = ' num2str(meanud2_h,'%0.2f') ...
        '; Mean V = ' num2str(meanvd2_h,'%0.2f')]);
text(-45,-30,['Vector speed = ' num2str(meanspd2_h,'%0.1f') ... 
        '; Direction = ' num2str(angleXtoCompass(meandird2_h),'%0.0f') '^o']);
text(-45,-38,['Major Axis = ' num2str(majaxd2_h,'%0.1f') ...
        '; Minor Axis = ' num2str(minaxd2_h,'%0.1f')]);
text(-45,-45,['Prin. Orient. = ' num2str(majazd2_h,'%0.0f') '^o']);
l = line([-50 50], [0 0]);
l2 = line([0 0], [-50 50]);
set(l,'color','k', 'linestyle','-.');
set(l2,'color','k', 'linestyle','-.');
set(gca,'xtick',labels);
set(gca,'ytick',labels);
set(gca,'xticklabel','');
set(gca,'yticklabel','');


%Low-Passed, Depth 2
subplot(2,2,4)
plot(ud2_lp,vd2_lp,'.')
hold on
plot([0 meanud2_lp],[0 meanvd2_lp],'r','linewidth',2);
plot(xxxd2_lp+meanud2_lp,yyyd2_lp+meanvd2_lp,'r','linewidth',2);
axis equal
axis tight
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
text(-45,-23,['Mean U = ' num2str(meanud2_lp,'%0.2f') ...
        '; Mean V = ' num2str(meanvd2_lp,'%0.2f')]);
text(-45,-30,['Vector speed = ' num2str(meanspd2_lp,'%0.1f') ... 
        '; Direction = ' num2str(angleXtoCompass(meandird2_lp),'%0.0f') '^o']);
text(-45,-38,['Major Axis = ' num2str(majaxd2_lp,'%0.1f') ...
        '; Minor Axis = ' num2str(minaxd2_lp,'%0.1f')]);
text(-45,-45,['Prin. Orient. = ' num2str(majazd2_lp,'%0.0f') '^o']);
l = line([-50 50], [0 0]);
l2 = line([0 0], [-50 50]);
set(l,'color','k', 'linestyle','-.');
set(l2,'color','k', 'linestyle','-.');
set(gca,'xtick',labels);
set(gca,'ytick',labels);
set(gca,'xticklabel','');
set(gca,'yticklabel','');

else
%Hourly Averaged, Depth 2
subplot(2,2,3)
axis equal
axis tight
box on
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
set(gca,'xtick',labels);
set(gca,'ytick',labels);


%Low-Passed, Depth 2
subplot(2,2,4)
box on
axis equal
axis tight
axis([-50 50 -50 50])
labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
set(gca,'xtick',labels);
set(gca,'ytick',labels);

end