% Plot all sorts of info on currents.
% Alex Wineteer
% Email: wineteer@jpl.nasa.gov, awinetee@calpoly.edu

%{
    This code will plot a number of plots. 
        Figure 1:   Geostrophic currents, MITgcm currents, and their
                    differences.

        Figure 2:   Calculated Ekman currents.

        Figure 3:   Geostrophic+ekman currents, MITgcm currents,
                    difference between them

        Figure 4:   Quiver plot of wind shear, calculated Ekman, difference
                    beetween MITgcm and calculated geostrophic currents.

    INPUTS: ii,  the averaging period desired (days).

    OUTPUTS: Plots as mentioned above. Also displays RMS differences.

%}

clc
close all
clear

siz = [480 795 1];
ii = 1.25;        % num days average desired

%% Calculated

% Original:
%Us = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/Eta/' num2str(ii) '_day_averages/surfGeoCur/Us_' num2str(ii) 'day_Eta_480x795x1.20110913T000000'],siz));
%Vs = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/Eta/' num2str(ii) '_day_averages/surfGeoCur/Vs_' num2str(ii) 'day_Eta_480x795x1.20110913T000000'],siz));

% Interp
Us = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/Eta/' num2str(ii) '_day_averages/surfGeoCur/mitGridU/UsMIT_Us_' num2str(ii) 'day_Eta_480x795x1.20110913T000000'],siz));
Vs = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/Eta/' num2str(ii) '_day_averages/surfGeoCur/mitGridV/VsMIT_Vs_' num2str(ii) 'day_Eta_480x795x1.20110913T000000'],siz));

%close all
figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
hh = tight_subplot(3,3);

axes(hh(1))
quikpcolor(Vs')
title(['Vs ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

axes(hh(2))
quikpcolor(Us')
title(['Us ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.3 .6])

calcMag = sqrt(Vs.^2 + Us.^2);
axes(hh(3))
quikpcolor(calcMag')
title(['Calc Mag ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .7])

%% MIT 

Vmit =(readbin(['~dmenemen/llc_4320/regions/Osmosis3/V/' num2str(ii) '_day_averages/' num2str(ii) 'day_V_480x795x1.20110913T000000'],siz));
Umit = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/U/' num2str(ii) '_day_averages/' num2str(ii) 'day_U_480x795x1.20110913T000000'],siz));

axes(hh(4))
quikpcolor(Vmit')
title(['V MIT ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

axes(hh(5))
quikpcolor(Umit')
title(['U MIT ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.3 .6])

mitMag = sqrt(Vmit.^2 + Umit.^2);

axes(hh(6))
quikpcolor(mitMag')
title(['MIT Mag ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .7])

diffU = Umit - Us;
diffV = Vmit - Vs;
diffMag = mitMag - calcMag;

axes(hh(8))
quikpcolor(diffU')
title(['Diff U (MIT - Calc) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh(7))
quikpcolor(diffV')
title(['Diff V (MIT - Calc) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh(9))
quikpcolor(diffMag')
title(['Diff Mag (MIT - Calc) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

img = getframe(gcf);
imwrite(img.cdata, ['figures/geoMITDiff_' num2str(ii) 'day.png']);

%% Ekman

Vek =(readbin(['~dmenemen/llc_4320/regions/Osmosis3/oceTAUY/' num2str(ii) '_day_averages/ekmanV/V_ekman' num2str(ii) 'day_oceTAUY_480x795.20110913T000000'],siz));
Uek = (readbin(['~dmenemen/llc_4320/regions/Osmosis3/oceTAUX/' num2str(ii) '_day_averages/ekmanU/U_ekman' num2str(ii) 'day_oceTAUX_480x795.20110913T000000'],siz));

figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
hh2 = tight_subplot(2,3);

axes(hh2(2));
quikpcolor(Uek')
title(['U Ekman ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh2(1));
quikpcolor(Vek')
title(['V Ekman ' num2str(ii) ' day averages'])
colorbar('southoutside')
caxis([-.1 .1])

magEk = sqrt(Vek.^2 + Uek.^2);
axes(hh2(3));
quikpcolor(magEk')
title(['Ekman Mag ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .2])

axes(hh2(5));
quikpcolor(diffU')
title(['U Diff ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh2(4));
quikpcolor(diffV')
title(['V Diff ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh2(6));
quikpcolor(diffMag')
title(['Diff in Mag  ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .2])

img = getframe(gcf);
imwrite(img.cdata, ['figures/ekmanDiff_' num2str(ii) 'day.png']);
%% Total Currents Estimated

Vtot = Vek + Vs;
Utot = Uek + Us;
magTot  = sqrt(Vtot.^2 + Utot.^2);

diffU2 = Umit - Utot;
diffV2 = Vmit - Vtot;
diffMag2 = mitMag - magTot;


fig = figure('color','white','units','normalized','outerposition',[0 0 1 1],'name','1');
hh3 = tight_subplot(3,3);

axes(hh3(1));
quikpcolor(Vtot')
title(['V Calc Total ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

axes(hh3(2));
quikpcolor(Utot')
title(['U Calc Total ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

magEk = sqrt(Vek.^2 + Uek.^2);
axes(hh3(3));
quikpcolor(magTot')
title(['Calc Total Mag ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .7])

axes(hh3(4))
quikpcolor(Vmit')
title(['V MIT' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

axes(hh3(5))
quikpcolor(Umit')
title(['U MIT ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.4 .4])

mitMag = sqrt(Vmit.^2 + Umit.^2);
axes(hh3(6))
quikpcolor(mitMag')
title(['MIT Mag ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([0 .7])

axes(hh3(8))
quikpcolor(diffU2')
title(['Diff U (MIT - CalcTot) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh3(7))
quikpcolor(diffV2')
title(['Diff V (MIT - CalcTot) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

axes(hh3(9))
quikpcolor(diffMag2')
title(['Diff Mag (MIT - CalcTot) ' num2str(ii) ' day averages']);
colorbar('southoutside')
caxis([-.1 .1])

img = getframe(gcf);
imwrite(img.cdata, ['figures/calcTotalMITDiff_' num2str(ii) 'day.png']);

%% Statistics

% rms
rmsV = sqrt(mean(diffV(~isnan(diffV)).^2));
rmsV2 = sqrt(mean(diffV2(~isnan(diffV2)).^2));
rmsU = sqrt(mean(diffU(~isnan(diffU)).^2));
rmsU2 = sqrt(mean(diffU2(~isnan(diffU2)).^2));
rmsMag = sqrt(mean(diffMag(~isnan(diffMag)).^2));
rmsMag2 = sqrt(mean(diffMag2(~isnan(diffMag2)).^2));

format short
disp('---------RMS Differences----------')
disp(['RMS diff V No Ekman     = ' num2str(rmsV) ' m/s']);
disp(['RMS diff V With Ekman   = ' num2str(rmsV2) ' m/s']);
disp(['RMS diff U No Ekman     = ' num2str(rmsU) ' m/s']);
disp(['RMS diff U With Ekman   = ' num2str(rmsU2) ' m/s']);
disp(['RMS diff Mag No Ekman   = ' num2str(rmsMag) ' m/s']);
disp(['RMS diff Mag With Ekman = ' num2str(rmsMag2) ' m/s']);

% angles
Ty =(readbin(['~dmenemen/llc_4320/regions/Osmosis3/oceTAUY/' num2str(ii) ...
    '_day_averages/' num2str(ii) 'day_oceTAUY_480x795.20110913T000000'],siz));
Tx =(readbin(['~dmenemen/llc_4320/regions/Osmosis3/oceTAUX/' num2str(ii) ...
    '_day_averages/' num2str(ii) 'day_oceTAUX_480x795.20110913T000000'],siz));

UVek = cat(3,Uek,Vek);
T = cat(3,Tx,Ty);
normUV = zeros(siz(1),siz(2));
normT = normUV;
avgN = ii;
for ii = 1 : siz(1)
    for jj = 1 : siz(2)
        
        normUV(ii,jj) = sqrt(UVek(ii,jj,1).^2 + UVek(ii,jj,2).^2);
        normT(ii,jj) = sqrt(T(ii,jj,1).^2 + T(ii,jj,2).^2);
    end
end

ang = real(acosd(dot(UVek,T,3)./(normUV.*normT)));
lats = (readbin('~dmenemen/llc_4320/regions/Osmosis3/grid/YC_480x795',siz));
lons = (readbin('~dmenemen/llc_4320/regions/Osmosis3/grid/XC_480x795',siz));

figure('color','white')
set(gcf,'renderer','zbuffer')
pp = pcolor(lons',lats',ang');
set(pp,'edgecolor','none');
shading interp
colorbar
caxis([20,60])
colormap cool

xx = 1:25:siz(1);
yy = 1:25:siz(2);

hold on
quiver(lons(xx,yy)',lats(xx,yy)',T(xx,yy,1)',T(xx,yy,2)','r')
quiver(lons(xx,yy)',lats(xx,yy)',UVek(xx,yy,1)',UVek(xx,yy,2)','k')
quiver(lons(xx,yy)',lats(xx,yy)',diffU(xx,yy)',diffV(xx,yy)','g')
legend('Calc Angle Mag','Wind Shear','Ekman Calculated','MIT-Calc Diff')

img = getframe(gcf);
imwrite(img.cdata, ['figures/quiverSurfaceCurrents_' num2str(ii) 'day.png']);