close all
clear all


addpath('../matlab_toolbox');
osave = true;

indir='../exp_nemed/';
outdir = '../doc/ocean_modelling/fig/';
load ../../data/coastlines-split-4326/east_med.mat


robs=load([ indir 'obs.dat']);
rsim=load([ indir 'rfloat_nodiv.dat']);
%rsim=load([ indir 'rfloat_div2.dat']);
uvb=load([ indir 'uv_back_aviso.dat']);

uvr=load([indir 'uv_tot_nodiv.dat']);
%uvr=load([indir 'uv_tot_div2.dat']);
meshg=load([indir 'meshgrid_aviso.dat']);
mask = load([indir 'mask_aviso.dat']);

iech=1;
jech=1;
%LIM=[34.5 36 32.5 35];
LON = [33.5 35.6];
LAT = [32.5 35];
LIM= [LON LAT];

%nlon=87;
%nlat=58;
lon=unique(meshg(:,1));
lat=unique(meshg(:,2));
nlon=length(lon);
nlat=length(lat);

[Lon,Lat]=meshgrid(lon',lat);
Mask = reshape(mask,[nlon,nlat]);
Mask = Mask';

nfloat=length(unique(robs(:,2)));

%% To change unit
%change grid from unitless to m/s

R_earth=6371229; %in meters
delta_x=zeros(nlat,nlon);
delta_y=zeros(nlat,nlon);

for ilon=1:nlon-1;
    for jlat=1:nlat-1;
        
    delta_x(jlat,ilon)=R_earth*(2*pi/360)*(Lon(jlat,ilon+1)-Lon(jlat,ilon))...
        .*cos(2*pi*Lat(jlat,ilon)/360); 
    delta_y(jlat,ilon)=R_earth*(2*pi/360)*(Lat(jlat+1,ilon)-Lat(jlat,ilon));  
        
end 
end



delta_x(1:nlat,nlon)=delta_x(1:nlat,nlon);
delta_y(1:nlat,nlon)=delta_y(1:nlat,nlon);


delta_x(nlat,1:nlon)=delta_x(nlat-1,1:nlon);
delta_y(nlat,1:nlon)=delta_y(nlat-1,1:nlon);


%% Reshape wind data
%it, ilon,ilat,u,v
ntime=length(unique(uvb(:,1)));
Ub=reshape(uvb(:,end-1),nlat,nlon,ntime);
Vb=reshape(uvb(:,end),nlat,nlon,ntime);

Ubms = Ub .* repmat(delta_x,1,1,ntime);
Vbms = Vb .* repmat(delta_x,1,1,ntime);

%it, ilon,ilat,u,v
ntime=length(unique(uvr(:,1)));
Ur=reshape(uvr(:,end-1),nlat,nlon,ntime);
Vr=reshape(uvr(:,end),nlat,nlon,ntime);
Ur(repmat(Mask,1,1,ntime)==0)=0;
Vr(repmat(Mask,1,1,ntime)==0)=0;

%Ur(Ub<1e-15)=0;
%Vr(Vb<1e-15)=0;
Urms = Ur .* repmat(delta_x,1,1,ntime);
Vrms = Vr .* repmat(delta_x,1,1,ntime);

%For the quiver legend
[~,jleg] = min(abs(lon-35.3));
[~,ileg] = min(abs(lat-32.8));
Urms(ileg,jleg,:) = 1;
%Ubms(ileg-1,jleg,:) = 1;

%% Plots 
figure(1)
clf
Hr=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Urms(1:jech:end,1:iech:end,1:end-1),3), ...
          mean(Vrms(1:jech:end,1:iech:end,1:end-1),3),'r');
hold on
Hb=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ubms(1:jech:end,1:iech:end,1:end-1),3),mean(Vbms(1:jech:end,1:iech:end,1:end-1),3),'b');




for j=1:nfloat     
  
  rfloat=rsim(rsim(:,2)==j-1&rsim(:,end)>0,[end-1 end]);
  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
  %Hsim=plot(rfloat_lon,rfloat_lat,'.-r');
    rfloat=robs(robs(:,2)==j-1&robs(:,end)>0,[end-1 end]);
  rfloat_lon=interp1(1:nlon,lon,rfloat(:,2)+1);
  rfloat_lat=interp1(1:nlat,lat',rfloat(:,1)+1);
  Hobs=plot(rfloat_lon,rfloat_lat,'m','Color',[0.3 0.3 0.3]);
  
end






%Ajout des trajectoires réelles
axis(LIM);
mapshow(S,'Color',[0 0 0]);
xlabel('Degrees East')
ylabel('Degrees North')

legend([Hb Hr Hobs],'Averaged over 9 days : background','Corrected','Drifter trajectory','Location','SouthWest');
text(35.33,32.75,'1 m/s')
%print -dpng ../exp_forw_dan/simu_drifter_dan.png

if osave
    print('-dpng','-r300',[outdir 'Eddy_velocity.png']);
    savefig([outdir 'Eddy_velocity.fig']);
end


%% eddy
cax = [-10 10].*1e-10;
cmap = scol_dif_lin(100,cax);

figure(2);
clf
[sX,sY] = meshgrid(34.3, 33:34.5);
%streamline(Lon,Lat,mean(Ut(:,:,end-1),3),mean(Vt(:,:,end-1),3),sX,sY)
[Uxb,Uyb]=gradient(mean(Ub(:,:,1:end-1),3),1,1);
[Vxb,Vyb]=gradient(mean(Vb(:,:,1:end-1),3),1,1);

%Uxt = Uxt .* (pi/180) .* R .* cos(Lat * pi/180);
%Vxt = Vxt .* (pi/180) .* R .* cos(Lat * pi/180);
%Uyt = Uyt .* (pi/180) .* R ;
%Vyt = Vyt .* (pi/180) .* R ;


owb = okubo_weiss(Uxb,Uyb,Vxb,Vyb);

imagesc(lon',lat,owb(1:iech:end,1:iech:end));
caxis(cax);
axis xy
shading interp
hold on
    mapshow(S,'Color',[0 0 0]);
    colorbar
axis(LIM);
colormap(cmap);
[Uxr,Uyr]=gradient(mean(Ur(:,:,1:end-1),3),1,1);
[Vxr,Vyr]=gradient(mean(Vr(:,:,1:end-1),3),1,1);

%Uxr = Uxr .* (pi/180) .* R .* cos(Lat * pi/180);
%Vxr = Vxr .* (pi/180) .* R .* cos(Lat * pi/180);
%Uyr = Uyr .* (pi/180) .* R ;
%Vyr = Vyr .* (pi/180) .* R ;

owr = okubo_weiss(Uxr,Uyr,Vxr,Vyr);
title('Okubo-Weiss parameter of the Background');
xlabel('Degrees East')
ylabel('Degrees North')

if osave
    print('-dpng','-r300',[outdir 'okubo_weiss_aviso.png']);
    savefig([outdir 'okubo_weiss_aviso.fig']);
end

figure(3)
clf
imagesc(lon',lat,owr(1:iech:end,1:iech:end));
caxis(cax);
shading interp
axis xy
hold on
    mapshow(S,'Color',[0 0 0])
    colorbar
axis(LIM);
colormap(cmap);
title('Okubo-Weiss parameter of the corrected field');
xlabel('Degrees East')
ylabel('Degrees North')

if osave
    print('-dpng','-r300',[outdir 'okubo_weiss_analyse.png']);
    savefig([outdir 'okubo_weiss_analyse.fig']);
end