clear all; close all; clc;
%% Prepeparation
clr = flipud(brewermap(15,'Spectral'));
clri= interp1(1:1:15,clr,1:0.25:15,'linear');
clr1= brewermap(15,'YlOrRd');
clri1=interp1(1:1:15,clr1,1:0.25:15,'linear');

%% Set up
HOME = '/Users/aahmed78/Desktop/EEPS 1520/Paper 03/CODES/SWOT/';
cd(strcat(HOME))
fileName = 'SWOT_L2_LR_SSH_Expert_001_009_20160901T065134_20160901T074301_DG10_01.nc';
fileInfo = ncinfo(fileName)
ncdisp(fileName);
fileInfo.Variables.Name;      
ssa = 'simulated_true_ssha_karin';
SSA = ncread(fileName,ssa);

omega=2*pi/(24*60*60);          % Earth rotation angle velocity
g=9.807;                        % Gravitational acceleration

%% Lat Lon bounding Box
% Extracting the longitudes and latitudes
lon1 = ncread(fileName,'longitude');
lat1 = ncread(fileName,'latitude');
lat_idx=find(lat1>=25 & lat1<=50);
lon_idx=find(lon1>=298 & lon1<=325);

region=intersect(lat_idx,lon_idx);
[row2,col2]=ind2sub(size(SSA(:,:)),region);

row1 = min(row2);
row11 = max(row2);
col1 = min(col2);
col11 = max(col2);

lat=lat1(row1:row11,col1:col11);
lon=lon1(row1:row11,col1:col11);


SSH=SSA(row1:row11,col1:col11);
min_SSHA = min(min(SSH));
max_SSHA = max(max(SSH));
[x,y]=ll2psn(lat,lon);

%% Calculating geostrophic velocity
cd(strcat(HOME,'/m_map'))
u=nan(size(SSH));
v=nan(size(SSH));
for i=2:size(SSH,1)-1
    for j=2:size(SSH,2)-1
        dx=double(m_lldist([lon(i+1) lon(i-1)],[lat(j) lat(j)]))*1000;
        dy=double(m_lldist([lon(i) lon(i)],[lat(j+1) lat(j-1)]))*1000;
        f=2*omega*sin((lat(j)*(pi/180)));
        u(i,j)=-g/f*(SSH(i,j+1)-SSH(i,j-1))/dy; % u = -(g/f)*(dn/dy)
        v(i,j)= g/f*(SSH(i+1,j)-SSH(i-1,j))/dx; % v = (g/f)*(dn/dx)      
    end
end
geo_velocity = sqrt(u.^2 + v.^2);
cd('..')


%% Calculate vorticity
R=6400000;
u1=nan(size(u));
v1=nan(size(v));
 

for k=2:size(x,1)-1
    for l=2:size(y,2)-1
        dx1=(x(k+1,l)-x(k-1,l))*(R*cosd(y(k,l))*pi/180);
        dy1=(y(k,l+1)-y(k,l-1))*(R*pi/180);
        u1(k,l)=(u(k,l+1)-u(k,l-1))/dy1;
        v1(k,l)=(v(k+1,l)-v(k-1,l))/dx1;
    end
end
vorticity = v1-u1;                                 %vorticity = dv/dx-du/dy             
vorticity_normalized = vorticity./f;

%% Okubo-Weiss parameter (W = Sn^2 + Ss^2 - vorticity^2)
%Sn = du/dx-dv/dy
u_sn=nan(size(u));
v_sn=nan(size(v));

for i_sn=2:size(x,1)-1
    for j_sn=2:size(y,2)-1
        dx_sn=(x(i_sn+1,j_sn)-x(i_sn-1,j_sn))*(R*cosd(y(i_sn,j_sn))*pi/180);
        dy_sn=(y(i_sn,j_sn+1)-y(i_sn,j_sn-1))*(R*pi/180);
        u_sn(i_sn,j_sn)=(u(i_sn,j_sn+1)-u(i_sn,j_sn-1))/dx_sn;
        v_sn(i_sn,j_sn)=(v(i_sn+1,j_sn)-v(i_sn-1,j_sn))/dy_sn;
    end
end
Sn = u_sn - v_sn;          

% Ss = dv/dx+du/dy
u_ss=nan(size(u));
v_ss=nan(size(v));

for i_ss=2:size(x,1)-1
    for j_ss=2:size(y,2)-1
        dx_ss=(x(i_ss+1,j_ss)-x(i_ss-1,j_ss))*(R*cosd(y(i_ss,j_ss))*pi/180);
        dy_ss=(y(i_ss,j_ss+1)-y(i_ss,j_ss-1))*(R*pi/180);
        u_ss(i_ss,j_ss)=(u(i_ss,j_ss+1)-u(i_ss,j_ss-1))/dy_ss;
        v_ss(i_ss,j_ss)=(v(i_ss+1,j_ss)-v(i_ss-1,j_ss))/dx_ss;
    end
end
Ss = v_ss + u_ss;

W = Sn.^2 + Ss.^2 - vorticity.^2;
W01 = std(W,0,2,'omitnan');
W01_domain = W01';
W0 = 0.2*std(W01_domain,0,2,'omitnan')
%% Plotting the SSH(m) and Geostrophic velocity(m/s)
close all

figure(1),
ax(1) = subplot(1,1,1);
[a1, hContour]=contourf(lon,lat,SSH,25,'edgecolor','none');
title('SSH (m)','fontweight','bold')
xlabel('Longitude (°W)','fontweight','bold'); ylabel('Latitude (°N)','fontweight','bold')
colormap(gca,parula); colorbar; shading interp
caxis([min(SSH(:)), max(SSH(:))]);
setup_plots(ax(1))

figure(2),
ax(2) = subplot(1,1,1);
[a2, hContour]=contourf(lon,lat,geo_velocity./15,25,'edgecolor','none');
title('V_g(m/s)','fontweight','bold')
xlabel('Longitude (°E)','fontweight','bold'); ylabel('Latitude (°N)','fontweight','bold')
colormap(gca,hot); colorbar; shading interp
hold all;
setup_plots(ax(2))
%%
figure(3), 
close all
ax(3) = subplot(1,1,1);
[a3, hContour]=contourf(lon,lat,vorticity_normalized,25,'edgecolor','none');
title('Vorticity Normalized by the Coriolis Factor, f','fontweight','bold')
xlabel('Longitude (°E)','fontweight','bold'); ylabel('Latitude (°N)','fontweight','bold')
colormap(gca,hot); colorbar; shading interp
hold all; 
caxis([-1.67e-4 -1.65e-3])
setup_plots(ax(3))
%%
figure(4),
ax(4) = subplot(1,1,1);
[a4, hContour]=contourf(lon,lat,W,25,'edgecolor','none');
title('Okubo-Weiss field','fontweight','bold')
xlabel('Longitude (°E)','fontweight','bold'); ylabel('Latitude (°N)','fontweight','bold')
colormap(gca,hot); colorbar; shading interp
hold all;
setup_plots(ax(4))



%% Projection
cd(strcat(HOME,'m_map'))
lon = lon-180;
figure(2),
bndry_lon=[128.8 128.8 128.3 128 126.8 126.6 128.8];
    bndry_lat=[49      50.33  50.33  50   49.5   49     49];
    
    clf;
    m_proj('lambert','long',[-134 -108],'lat',[20 50],'rectbox','on');

    m_gshhs_i('color','k');              % Coastline...
    m_gshhs_i('speckle','color','k');    % with speckle added

    m_line(bndry_lon,bndry_lat,'linewi',2,'color','b');     % Area outline ...
    m_hatch(bndry_lon,bndry_lat,'single',30,5,'color','b'); % ...with hatching added.

    m_grid('linewi',2.5,'linest','none','tickdir','out','fontsize',18);
    %title({'Location of California current system and SWOT swath'});
    m_text(-128,24,{'Pacific','Ocean'},'fontsize',18);
    
    m_northarrow(-113,42,5,'type',3,'linewi',2);    
    clc;


%% Oblique Mercator Projection with quiver and contour data
figure(3),
clf
%[lon,lat]=meshgrid([-136:2:-114],[36:2:54]);

m_proj('oblique','lat',[50 0],'lon',[-144 -108],'aspect',.8);
m_coast('patch',[.7 .7 .7]);
m_grid('tickdir','out','yaxislocation','right',...
            'xaxislocation','top','xlabeldir','end','ticklen',.01);

%% Bathymetry
figure(4),
clf
 
m_proj('lambert','lat',[20 50],'long',[-134 -108]);

set(gcf,'color','w')   % Set background colour before m_image call

caxis([-6000 0]);
colormap(flipud([flipud(m_colmap('blues',10));m_colmap('water',118)]));
m_etopo2('shadedrelief','gradient',3);
 
m_gshhs_i('patch',[.8 .8 .8]);
 
m_grid('box','fancy');

ax=m_contfbar(.97,[.5 .9],[-6000 0],[-6000:100:000],'edgecolor','none','endpiece','no');
xlabel(ax,'meters','color','k');

















