%clear all
path_pg=pwd;
cd 'E:\Gatti2\golts\Survey\golts_2002'

path=uigetfile('*.mat','Quel fichier Thetis II tt bo ki marche?');
load(path)
[lignes colones]=size(ensemble);
deb=100;
bin=9;
fin=length(ensemble(:,13))-200;
long_nanise=ensemble(deb:fin,14);
lat_nanise=ensemble(deb:fin,13);
Vest_nanise=eastCurrentComponent(deb:fin,bin)/100;
Vnord_nanise=northCurrentComponent(deb:fin,bin)/100;
indices_propres=find([isnan(Vnord_nanise)]==0);%*************A rajouer sur le l'est:NON!!!*************

long_brut=long_nanise(indices_propres);
lat_brut=lat_nanise(indices_propres);
Vest_brut=Vest_nanise(indices_propres);
Vnord_brut=Vnord_nanise(indices_propres);
long_bolk=20;
long=long_brut;%blokaverage(long_brut,long_bolk);
lat=lat_brut;%blokaverage(lat_brut,long_bolk);
cd (path_pg)
Vest=blokaverage(Vest_brut,long_bolk);
Vnord=blokaverage(Vnord_brut,long_bolk);
%Vest=Vest-mean(Vest);Vnord=Vnord-mean(Vnord);
% deltaX=std(long);
% deltaY=std(lat);
% for i=1:fin-1
%     hx(i)=ensemble(i,14)-ensemble(i+1,14);
%     hy(i)=ensemble(i,13)-ensemble(i+1,13);
% end

% psi=hx.*Vest'+hy.*Vnord';%pb indice
% [PX,PY]=gradient(psi,hx,hy);
Ctyp = input('Enter 1 (for Exp) or 2 (for Gaussian) Cov...');
L    = input('Enter covariance function e-folding scale...');
namstr='alex';
%callstart		% Cleans up all variables and figure windows
inst = ['   To make a hardcopy of any of the current figures type:'
        '        figure(#)     (# is the desired figure number)   '
        '        print         (send figure to default printer)   '
        '                                                         '
        '   To continue type: r-e-t-u-r-n                         '];

inst2= ['   At this point you may loop back to change mapping     '
        '   parameters without changing selected data points by   '
        '   entering a "1" below.                                 '
        '       If you enter a number other than 1, you will      '
        '       have the option of restarting "oa" to select      '
        '       different data points or exiting Matlab           '
        '           (note: you may also print hardcopies          '
        '            at this stage before exiting)                ']; 

 
  
  
 
 x=linspace(4.6,6,40);
x=6378*cos(43*pi/180)*x*pi/180;X1=x(1);x=(x-x(1));
 y=linspace(42.4,43.3,40);
 y=6378*y*pi/180;Y1=y(1);y=(y-y(1));
 rtest=0:0.5:100;

% Get user name for plot identification
% hardwire height and width of analytical function for now



x0=0;
y0=0;
[X,Y]=meshgrid(x,y);
R=(X .^ 2 + Y .^ 2);
xI=4.8:0.01:5.9;
yI=42.1:0.01:43.2;
zI=0:-5:-200;
for direction=1:2
%T_est = NaN .* zeros(xI,yI);
%T_err = NaN .* zeros(xI,yI);
T_est = NaN .*X;
T_err = NaN .*X;

%%plot(X,Y,'g+')

xb=[];
yb=[];


   
   xb=6378*cos(43*pi/180)*long*pi/180;xb=xb-X1;
   yb=6378*lat*pi/180;yb=yb-Y1;
   
%  xb=[4.6 4.7 4.9 5 5.4 5.6]';
%  yb=[42.1 42.3 42.4 42.5 42.6 43]';
%  Vest=[1 0.2 0.3 0.6 0.4 0.8]';
%  Vnord=[0 0.4 0.2 0.5 0.6 0.2]';




%loop here to select different mapping parameters using the
%same data points (restart to select new data points)
iloop=1;
%while iloop == 1,

% Choose exponential or gaussian covariance function and
% Length scale parameter

sl=int2str(L);
if Ctyp == 1
  disp(['Exponential Covariance Function Selected with L=' sl])
else
  disp(['Gaussian Covariance Function Selected with L=' sl])
end
amp = 1; % use unit amplitude for now (should be variance)

disp('Calculating Estimates........................')

% Calculate the distances between all observations.
obslen=length(xb);
s = zeros(obslen,obslen);
for j = 1:obslen
    s(j,:) = sqrt( (xb-xb(j)) .^2 + (yb-yb(j)) .^2 )'; 
end

% Calculate Cdd and Cddi and Cxx outside loop.

if Ctyp == 1
  Cdd = amp*exp(-s ./ L);
  Crr = amp*exp(-rtest ./ L);
else
  Cdd = amp*exp(-(s.^2) ./ L);
  Crr = amp*exp(-(rtest.^2) ./ L);
end
Cddi = inv(Cdd);
if direction==1
d=Vest;
else
   d=Vnord; 
end
for xi = 1:1:length(x)
    for yi = 1:1:length(y)
    
% xs = sqrt( (xb - X(xi)) .^2 + (yb - Y(yi)) .^ 2 );  %%%ERROR LINE
        xs = sqrt( (xb - X(xi,yi)) .^2 + (yb - Y(xi,yi)) .^ 2 );
        xs = xs';

    if Ctyp == 1   
        Cxd = amp*exp(-xs ./ L);
    else
        Cxd = amp*exp(-(xs.^2) ./L); 
    end
        A = Cxd * Cddi;
        T_est(xi,yi) = A * d;
		T_err(xi,yi) = 1 - A * Cxd';
        
%        if  rem(yi,10) == 0  &  rem(xi,10)  == 0
%           figure(3)
%           plot(Cxd)
%           disp([' ',int2str(X(xi,yi)),' ',int2str(Y(xi,yi))])

%           keyboard
%        end
        
    end
end

figure(3)
% subplot(2,2,1)
% %cc=contour(X,Y,T);
% %clabel(cc)
% %title(['Plan View with Obs Points (' namstr ')'])
% hold on
% %%%%plot(X,Y,'g.')
% xlabel('X Location')
% ylabel('Y Location')
% hold off
% line('XData',xb,'YData',yb,'LineStyle','*','Color','r');
% set(gca,'aspectratio',[1,1]);


subplot(2,2,2),plot(rtest,Crr)
xlabel('Separation Distance')
ylabel('Covariance Function') 
if Ctyp ==1
  title(['Exp Cov Function, L=', sl ' (' namstr ')'])
else
  title(['Gauss Cov Function, L=', sl ' (' namstr ')'])
end

subplot(2,2,3)
mesh(X,Y,T_est);
xlabel('X Location')
ylabel('Y Location')
if Ctyp ==1
  title(['OA Map, Exp,L=' sl ' (' namstr ')'])
else
  title(['OA Map, Gauss,L=' sl ' (' namstr ')'])
end

subplot(2,2,4)
mesh(X,Y,T_err);
xlabel('X Location')
ylabel('Y Location')
if Ctyp ==1
  title(['OA Error, Exp,L=' sl ' (' namstr ')'])
else
  title(['OA Error, Gauss,L=' sl ' (' namstr ')'])
end

orient tall
% wysiwyg
figure 
mesh(X,Y,T_est)

if direction==1
u=T_est;
erru=T_err;
else
 v=T_est; 
 errv=T_err;
end

% Compute rms difference between estimate and analytical fields
% as a measure of success
end
xb=6378*cos(43*pi/180)*long*pi/180;xb=xb-X1;
   yb=6378*lat*pi/180;yb=yb-Y1;
X=(X+X1)/(6378*cos(43*pi/180)*pi/180);
Y=(Y+Y1)/(6378*pi/180);
xb=(xb+X1)/(6378*cos(43*pi/180)*pi/180);
yb=(yb+Y1)/(6378*pi/180);

figure; 
golts
hold on
err=abs(erru)+abs(errv);
%pcolor(X,Y,err),colorbar

quiver(X,Y,u,v,2) 
quiver(xb,yb,Vest,Vnord,'r')
[k,kk]=flowfun(u,v);
contour(X,Y,kk)
text(5,43.25,path(length(path)-17:length(path)-4))
figure 
pcolor(X,Y,divergence(X,Y,u,v)),colorbar
hold on
quiver(xb,yb,Vest,Vnord,'r')
figure
pcolor(X,Y,curl(u,v)),colorbar
hold on
quiver(xb,yb,Vest,Vnord,'r')
figure 
[c,h] =contour(X,Y,err),colorbar;
clabel(c,h), colorbar
%quiver(long_nanise,lat_nanise,Vest_nanise,Vnord_nanise,'r')

%rmserr = sqrt(mean(mean((T-T_est).^2)));
%trms = num2str(rmserr);
