clear
close all

load_path
addpath /automount1/mgoes/matlab/Hilbert_Huang
addpath /automount1/mgoes/matlab/pcatool
addpath /automount1/mgoes/matlab/mapcolor/colors/
addpath /automount1/mgoes/matlab/m_map/
global lat lon

display('LOAD data...')
nc=netcdf('../../AVISO/adt_week_SA.cdf');

LON = nc{'NBLONGITUDESN179_61'}(:);

LAT = nc{'NBLATITUDES307_413'}(:);

time = nc{'TIME'}(:);
jul_yr = time/24 + julian(1950,1,1,0);
month = gregorian2month(gregorian(jul_yr));
ydec = month2yfraction(month);
days = (time-time(1))/24;
Ts = 7;%days(2)

wave =1;
if wave
    filename = '../wavelet/ADT_wavelet.mat';
    filename = '../wavelet/ADT_0_5periods.mat';
    varId = 'ADT1'%'ADT0'
%     load('../wavelet/ADT_wavelet.mat',varId);
%     eval(['ADT = ',varId,';']);clear(varId)
%       varId = 'ADT1'%'ADT0'
%load('../wavelet/ADT_wavelet.mat','ADT3')
%ADT = ADT3;
%DT= 25;
%    load(filename,'ADT3','ADT4');
%    ADT = ADT3 +ADT4; clear ADT3 ADT4
    
    load(filename,'ADT4'); ADT = ADT4;
   % load(filename,'ADT1','ADT2');
   % ADT = ADT1 + ADT2; clear ADT1 ADT2
    %load(filename,'ADT0')
    %ADT = ADT0 +ADT1 +ADT2;clear ADT0 ADT1 ADT2
  %  DT = 2;%2; %ADT0
    %DT = 6; %ADT1
    DT = 13; %ADT2
else
ADT = nc{'ADT'}(:);
MISS = nc{'ADT'}.('_FillValue')(:);
clear nc 

ADT(ADT==MISS)=nan;
a = apply_filter(ADT,Ts,365);
ADT = ADT - a; clear a
adt_clim = make_clim(ADT,month);
ADT = ADT - adt_clim;
DT = 7;%4; 
end
  %  ADT = bsxfun(@rdivide,ADT,std(ADT));

[nt,ny,nx] = size(ADT);
display('SELECT REGION...')
region = 'western';%'western'

%DEFINE NEW REGION
switch region
    case 'interior'
      vecx = 80:180; %interior
      vecy = 30:65;%28:38S
      DT = 7;
    case 'eastern'
      vecx = 180:nx; %eastern  
      vecy = 30:65;%28:38S
      DT=4;
    case 'western'
       vecx = 20:80; %western
        vecx = 20:70
      vecy = 50:90; % BC 20S-32S   %FICA LEGAL Western
      % vecy = 30:65;%BC-Malvinas
     %vecy = 70:100;  %17S-23S - more north centered at 23S
    %  vecx = 10:100;  vecy = 20:107

    %   vecx = 20:70;vecy = 50:100;  %Fernandes et al 2009
      DT = 1;%2;%1;%Kelvin
      %DT = 4;%Rossby
     % DT = 52;%25;
    case 'basin'
       vecx = 20:nx; %All basin
       vecy = 40:85;
        DT = 25;
end

        
% %vecx = 10:70;
% vecx = 20:80; %western
% %vecx = 20:180;
% vecx = 20:nx; %All basin
% %vecx = 180:nx; %eastern
% vecx = 80:180; %interior
% 
% vecy = 30:65;%BC-Malvinas
% %vecy = 40:85;
% %vecy = 50:90; % BC 20S-32S   %FICA LEGAL Western

vect = 1:length(time);



%ADT = ADT(1:100,:,:);
ADT = ADT(vect,vecy,vecx);
lon = LON(vecx);
lat = LAT(vecy);

nt = size(ADT,1);
nl = numel(ADT(1,:,:));
vec = find(~isnan(ADT(1,:)));

mask = ~isnan(squeeze(ADT(1,:,:)));
mask2 = isnan(squeeze(ADT));
%mask2 = squeeze(prod(double(mask)));
[nlat,nlon] = size(mask);
% for ii=1:nt
%    mask3(ii,:,:) = mask2; 
% end
    
%mask(vecy,vecx)=0;



%keyboard
display('BEGIN LOOP...')

%EXTENDED EOFs
M = map2mat(mask,ADT);N=6;METHOD = 2;NLAG = 4;%DT = 2;%7;
M(isnan(M))=0;
[EOFs,PC,EXPVAR] = calEeof(M,N,METHOD,NLAG,DT);
%[EOFs,PC,EXPVAR] = calCeof(M,N,METHOD);
CHP = ones(N,NLAG,nlat,nlon);
for ii = 1:N
CHP(ii,:,:,:) = mat2map(mask,squeeze(EOFs(ii,:,:)));
    %    for jj = 1:NLAG
%        mask2(ii,jj,:,:) = squeeze(mask2(ii,jj,:,:)).*squeeze(mask(1,:,:));
%    end
end
ydec = ydec(1:nt-(NLAG*DT)+1);
%CHP = mat2map(mask,EOFs);
Dt = 7;
%dispEEOF_aviso(lon,lat,CHP,PC,EXPVAR,Dt,NLAG,1:N)
dispEEOF_aviso_figures(lon,lat,ydec,CHP,PC,EXPVAR,Dt*DT,NLAG,1:N)
if 0
    for ii=1:N
     save2eps(sprintf('%s%d%s%d%s%d%s','EEOF',ii,'_LAG',NLAG,region,wave,'.eps'),ii,300)
    end
end
return

IMF = cell(nl,6);
%V   = nan*ones(1,nl);
B   = nan*ones(nl,6);

for ii=vec %1:nl
ii
   x = ADT(:,ii);
   imf = emd(x);
   ni = length(imf);
   b=[];d=[];
   for k = 1:ni
     b(k) = sum(imf{k}.*imf{k});
     th   = angle(hilbert(imf{k}));
     d{k} = diff(th)/Ts/(2*pi);
   end
 %   D(1:ni,ii) = d;
%[~,v] = sort(-b);
%b     = 1-b/max(b);
      B(ii,1:length(b)) = b;
   IMF(ii,1:length(imf)) = imf;

end

