source: altifloat/matlab_toolbox/script_perso.m @ 160

close all
clear all
indir='../exp_hr/';
%indir='../exp_twin_aviso/';
%indir='../../code_leila/altifloat/obs_float/';
%indir='../obs_float/';
%ut=load([ indir 'uzero.dat']);
%vt=load([ indir 'vzero.dat']);
robs=load([ indir 'obs_jum_cut_18.dat']);
rsim=load([ indir 'rfloat_total_jum1.dat']);
uvb=load([ indir 'uv_back_aviso_0901_0902_sp_t2.dat']);
uvt=load([ indir 'uv_truthdan_0901_0903_lin_t2.dat']);
uvr=load([indir 'uv_tot_jum1.dat']);
uvgeo=load([indir 'uv_geos_jum1.dat']);
meshg=load([indir 'meshgrid_aviso_t2.dat']);
Div_ana_file=[indir 'Div_ana.dat'];
Div_fg_file=[indir 'Div_fg.dat'];

%drifterdir='../altifloat/';
%fname={'a300234011043340.mat','a300234011120770.mat','a300234011249440.mat'};
%X1=ut(:,end-3);
%Y1=ut(:,end-2);
iech=1;
jech=1;
LIM=[34.5 36 32.5 35];
%Format : t/s, ilon, ilat, val

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

[Lon,Lat]=meshgrid(lon',lat);

%X1=X1+2;
%Y1=Y1+2

%Ut=reshape(ut(:,end),nlat,nlon);
%Vt=reshape(vt(:,end),nlat,nlon);



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



%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);

ntime=length(unique(uvt(:,1)));
Ut=reshape(uvt(:,end-1),nlat,nlon,ntime);
Vt=reshape(uvt(:,end),nlat,nlon,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);

if exist('uvgeo')
    ntime=length(unique(uvgeo(:,1)));
Ugeo=reshape(uvgeo(:,end-1),nlat,nlon,ntime);
Vgeo=reshape(uvgeo(:,end),nlat,nlon,ntime);
end


Hr=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ur(1:jech:end,1:iech:end,1:end-1),3), ...
          mean(Vr(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(Ub(1:jech:end,1:iech:end,1:end-1),3),mean(Vb(1:jech:end,1:iech:end,1:end-1),3),'b');

%Ht=quiver(Ut,Vt,'g');


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');

  %   plot(rfloat(:,2)+1,rfloat(:,1)+1,'.-m');

  % Hreal=plot(data.lon(1:4),data.lat(1:4),'.-k')
  
end


Ht=quiver(Lon(1:jech:end,1:iech:end),Lat(1:jech:end,1:iech:end),mean(Ut(1:jech:end,1:iech:end,1:end-1),3),mean(Vt(1:jech:end,1:iech:end,1:end-1),3),'m');
axis(LIM);



%Ajout des trajectoires réelles



%legend([Hb Hr],'background','retrieved','Location','SouthEast');
legend([Hb Hr Ht],'background.','retrieved','truth','Location','SouthEast');

%print -dpng ../exp_forw_dan/simu_drifter_dan.png


if (false)
delta_b=sqrt(mean((Ut(:,:,1:end-1)-Ub(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vb(:,:,1:end-1)),3).^2);
delta_r=sqrt(mean((Ut(:,:,1:end-1)-Ur(:,:,1:end-1)),3).^2+mean((Vt(:,:,1:end-1)-Vr(:,:,1:end-1)),3).^2);
delta_inc=sqrt(mean((Ub(:,:,1:end-1)-Ur(:,:,1:end-1)),3).^2+mean((Vb(:,:,1:end-1)-Vr(:,:,1:end-1)),3).^2);

cscale=[0 5e-4];

figure(2)
subplot(2,2,1)
imagesc(lon,lat,delta_b,cscale);
axis xy;axis(LIM);
hold on
for j=1:3
   data= load([drifterdir fname{j}]);
  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);
plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
end
title('Backround error');


subplot(2,2,2)
imagesc(lon,lat,delta_r,cscale);
axis xy;axis(LIM);
hold on
for j=1:3
   data= load([drifterdir fname{j}]);
  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);
plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
end

title('Analysis error');



subplot(2,2,4)
imagesc(lon,lat,delta_inc);
axis xy;axis(LIM);
colorbar
hold on
for j=1:3
   data= load([drifterdir fname{j}]);
  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);
plot([data.lon(1);rfloat_lon],[data.lat(1);rfloat_lat],'o-m');
end
title('Increment');


%print -f1 -dpng ../exp_forw_dan/twin_exp_uv_field.png
%print -f2 -dpng ../exp_forw_dan/twin_exp_errors.png
end


%% Erreur le long du trajet
if(false)
Id=0;
for id=Id
  figure(5+id)
  I=find(robs(:,2)==id-1);
  Xq=robs(I,end)+1;
  Yq=robs(I,end-1)+1;
  Tq=robs(I,1)+1;
  
  err_b=interp3(1:nlon,1:nlat,1:ntime,delta_b,Xq,Yq,Tq);
  err_r=interp3(1:nlon,1:nlat,1:ntime,delta_r,Xq,Yq,Tq);
  incr=interp3(1:nlon,1:nlat,1:ntime,delta_inc,Xq,Yq,Tq);
end %id
end %if false

%% Divergence
Div_ana=load(Div_ana_file);
Div_fg=load(Div_fg_file);

Div2_ana=reshape(Div_ana(:,end),nlat,nlon,ntime-1);
Div2_fg=reshape(Div_fg(:,end),nlat,nlon,ntime-1);

figure(3)
subplot(3,1,1)
imagesc(lon,lat,mean(Div2_fg,3),[-2 2]*1e-5);
axis xy;
axis(LIM);
title('Mean divergence Background');
colorbar


figure(3)
subplot(3,1,3)
imagesc(lon,lat,mean(Div2_ana,3),[-2 2]*1e-5);
axis xy;
axis(LIM);
title('Mean divergence Analyse');
colorbar

subplot(3,1,2);
imagesc(lon,lat,mean(Div2_ana,3)-mean(Div2_fg,3),[-2 2]*1e-5);
axis xy;
axis(LIM);
title('Mean difference');
 colorbar