clear
%close all

load_path
%addpath /automount1/mgoes/matlab/Hilbert_Huang
%addpath /automount1/mgoes/matlab/pcatool
%addpath /data/mgoes/matlab/mapcolor/colors/
%addpath /automount1/mgoes/matlab/m_map/
addpath /data/mgoes/SA_gyre/wavelet/Felipe
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 = {'ADT0','ADT1','ADT2','ADT3','ADT4'}
    varId = {'ADT3'}%,'ADT4'}
    varId = {'ADT1'}
   % filename =  '../wavelet/Felipe/SSA_4_filter.mat'
 %   filename =  '../wavelet/Felipe/SSA_4_filter_1.mat'  %Daily to 7 day
 %    varId = {'ADT2','ADT3'}%,'ADT3'} %ORIGINAL
 %   filename = '/phodnet/share/mgoes/SUDIP/ADT_0_5periods_good.mat'
    
     %  filename =  '../wavelet/Felipe/SSA_4_filter_HPS.mat'
 %   filename =  '../wavelet/Felipe/SSA_MESOSCALE_NEW.mat'

    % filename =  '../wavelet/Felipe/SSA_INTERANNUAL_ANOMALY.mat' %from 0.75-9
    filename =  '../wavelet/Felipe/SSA_INTERANNUAL_NEW.mat'        %DEFAULT
   varId = {'ADT0','ADT1','ADT2'} %ORIGINAL                        %DEFAUTL

 % filename = '/data/mgoes/SA_gyre/wavelet/Felipe/SSA_HP_AN_LP.mat'
   % varId = {'SSA_HP'}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    %  filename = '../wavelet/Felipe/ADT_5periods_new.mat';
      load(filename,'lon','lat');LAT=lat;LON=lon;
      load(filename,'time');month = julian2month(time);ydec = month2yfraction(month);
 %    varId = {'ADT5'};%'ADT2','ADT3'};%'ADT5'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    load(filename,'VecFilt')
func = sprintf('%s%s %s %s %s %s %s','load ',filename,varId{:,:})    
eval(func)
ADT = 0;
for dd = 1:length(varId)
  func=sprintf('%s%s%s','ADT = ADT + ',varId{dd},';');
  eval(func)
end

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 = 'basin';%'western';%'western'
%region = 'western'
%DEFINE NEW REGION
switch region
    case 'interior'
      vecx = 80:180; %interior
      vecy = 30:65;%28:38S
   %   vecy = 10:65;%28:45S
     
%      DT = 7;
    case 'eastern'
      vecx = 180:nx; %eastern  
      vecy = 30:65;%28:38S
%      DT=4;
    case 'western'
%        vecx = 20:80; %western
%       vecy = 50:90; % BC 20S-32S   %FICA LEGAL Western
%       % vecy = 30:65;%BC-Malvinas
%     vecx = 10:100;  vecy = 60:107
    
       vecx = find(LON>-60&LON<-25);
      % vecx = find(LON>-60&LON<5);
       vecy = find(LAT>-34.6&LAT<-20);
    % vecy = 70:100;  %17S-23S - more north centered at 23S
     %  vecx = 20:70;vecy = 50:100;  %Fernandes et al 2009
%      DT = 1;%2;%1;%Kelvin
      %DT = 4;%Rossby
     % DT = 25;
    case 'basin'
       vecx = 20:nx; %All basin
       vecy = 40:85;
       vecx = find(LON>-55&LON<20);%vecx=vecx(1:2:end);
   %    vecx = find(LON>-55&LON<15);%vecx=vecx(1:2:end);
       %vecy = find(LAT>-27&LAT<-21);%vecy=vecy(1:2:end);
       
       vecy = find(LAT>-34&LAT<-20);%vecy=vecy(1:2:end);   %FOR COEF
%       vecy = find(LAT>-34&LAT<-18);%vecy=vecy(1:2:end);
%        vecy = find(LAT>-45&LAT<-20);%vecy=vecy(1:2:end); 
      
 %      vecy = find(LAT>-36&LAT<-18);%vecy=vecy(1:2:end);
%            vecy = find(LAT>-36&LAT<-20);%vecy=vecy(1:2:end);   %FOR SLP
%            vecy = find(LAT>-36&LAT<-15);%vecy=vecy(1:2:end);

%        DT = 25;
%       vecy = find(LAT>-35&LAT<-18);%vecy=vecy(1:2:end);
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);
%ADT = bsxfun(@times,ADT,cos(lat(:,ones(1,nt))'*pi/180));  %area weight
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
disp('BEGIN LOOP...')

%COMPLEX EOFs
%M = map2mat(mask,ADT);N=6;METHOD = 2;
M = map2mat(mask,ADT);N=3;METHOD = 2;
M(isnan(M))=0;
[EOFs,PC,EXPVAR] = calCeof(M,N,METHOD);
%[PC2, r] = varimax( transpose(PC));
%PC = transpose(PC2);
%EOFs = transpose(EOFs)*r;

%[EOFs,PC,EXPVAR] = calEeof(M,N,METHOD,5,2);
EOFs = conj(EOFs);PC = conj(PC);

% oo = PC'*EOFs;
% reconstruction = real(oo);
% phase = imag(oo);
if 0
    load /data/mgoes/SEALEVEL/MONTHLY/mapa.mat
    %load /automount1/mgoes/SEALEVEL/MONTHLY/mapa.mat
   % ivv = 3;
 oo1 = transpose(PC([1 2 3],:))*conj(EOFs([1 2 3],:));
 rec1 = mat2map(mask,oo1);
 oo = transpose(PC([3],:))*conj(EOFs([3],:));
 rec = mat2map(mask,oo);
 [~,ilat] = min(abs(lat+34.5));
 [ilon] = find(lon>=-52.9 & lon<=-50.5);
 BC34 = real(rec(:,ilat,ilon(1))-rec(:,ilat,ilon(end)));
 BC341 = real(rec1(:,ilat,ilon(1))-rec1(:,ilat,ilon(end)));
 
figure(10)%,hold on
 plot(ydec,BC341)
 hold on
 load /data/mgoes/SA_gyre/AX18_shelf/Data/transport_from_altimetry_nomess.mat
 [bc] = lanczosfilter([tsp_bc2.transp_1000_minus]/1e6,1,1/680,501,'low');
 plot(year_total_mo,-1*(bc-nanmean(bc)));
 %stop
 if 1
 for ii=500:5:1200
     subplot(2,1,1)
 contourf(lon,lat,squeeze(real(rec(ii,:,:))),-10:.5:10)
 %contourf(lon,lat,squeeze((ADT(ii,:,:))),-5:.5:5)
 xlim([-55 20])
 colorbar;colormap(mapa)
 caxis([-2 2])
    subplot(2,1,2)
 contourf(lon,lat,squeeze(real(rec1(ii,:,:))),-10:.5:10)
 %contourf(lon,lat,squeeze((ADT(ii,:,:))),-5:.5:5)
 xlim([-55 20])
 colorbar
 title(month(ii,:))
 colormap(mapa)
  caxis([-3 3])
 drawnow
 end
 end %02
end
%stop
CHP = mat2map(mask,squeeze(EOFs));
theta = angle(CHP)*180/pi;
phi = angle(PC)*180/pi;
%CHP = ones(N,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
%SPATIAL PHASE
%theta = atan2d(imag(CHP),real(CHP));%*180/pi;

%theta = atan2d(real(CHP),imag(CHP));%this seems to be right
%r = abs(CHP);
%CHP2 = r.*exp(j*theta);

%TEMPORAL PHASE
%phi = atan2d(imag(PC),real(PC));%*180/pi;
%phi = atan(imag(PC)/real(PC))*180/pi;
%phi = atan2d(real(PC),imag(PC));%*180/pi; %THIS IS CORRECT
%stop
%%% CORRELATE TO SST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
correlate_SST
%%% CORRELATE TO SLP%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 stop
correlate_SLP
stop
correlate_indices
stop
 if 0
% addpath /data/mgoes/matlab/m_map
% oi=load('../wavelet/SLP/SLP_INTERANNUAL.mat');
% %oi=load('../wavelet/SLP/SLP_INTERANNUAL_ANOMALY.mat');
% load /data/mgoes/SA_gyre/wavelet/Felipe/mapa.mat
% p1 = linspace(0,1,length(mapa));
% mapa2 = interp1(p1,mapa,linspace(0,1,24));
% [X,Y]=meshgrid(oi.lon,oi.lat);
% 
% ygfrac=month2yfraction(julian2month(oi.time));
% 
% for nmod = 1:3
% [PC2,~,mo] = make_monthly_mean(real(PC(nmod,:)),month);
% CHP2 = squeeze(real(CHP(nmod,:,:)));
% ADT0=oi.ADT0(1:277,:,:);
% [nt0,ny0,nx0]=size(ADT0);
% C=zeros(ny0,nx0);P=nan*C;
% 
% % ADT0=ADT0(1:end-2,:,:);
% % PC2 = PC2(3:end);
% % ygfrac=month2yfraction(julian2month(oi.time(1:end-2)));
% 
% %DJF = find(str2num(mo(:,5:6))==12 | str2num(mo(:,5:6))<3);DJF=DJF(1:end-1);
% %MAM = find(str2num(mo(:,5:6))>=3 & str2num(mo(:,5:6))<6);
% for ii=1:ny0
% for jj=1:nx0
% [c,p] = corrcoef(PC2(:),squeeze(ADT0(:,ii,jj)));
% %[c,p] = corr2_AR(PC2',squeeze(ADT0(:,ii,jj)));
% c=c(2,1);p=p(2,1);
% C(ii,jj)=c;
% if p<0.05,P(ii,jj)=1;end
% end
% end
% 
% nfilt = 12;
% neff = nt/nfilt;
% t=abs(C).*sqrt( (neff-2)./(1-C.^2) );
% PVAL=2*tcdf( -t,neff-2 );
% P= double(PVAL<0.05); P(P==0)=nan;
% ui = find(~isnan(P));
% 
% figure;clf;
% set(gcf,'PaperPosition',[1 1 5 5]); 
% set(gcf,'DefaultLineLineWidth',1.5)
% wysiwyg;
% 
% %dois = subplot('Position',[0.1900 0.1100 0.6550 0.3412]);
% dois = subplot('Position',[0.1900 0.3100 0.6412 0.2 ]);%7550 ]);
% m_proj('miller','long',[-55 20],'lat',[-35 -20]);
% m_contourf(lon,lat,CHP2,'edgecolor','none')
% m_coast('patch',[.5 .5 .5])%,'linewidth',3);
% hold on
% m_contour(oi.lon,oi.lat,C,0:.1:1,'k')
% m_contour(oi.lon,oi.lat,C,-1:.1:0,'k--')
% m_contour(oi.lon,oi.lat,C,[0 0],'k','linewidth',3)
% m_text(-53,-23,'b.','backgroundcolor',[1 1 1],'fontweight','bold')
% m_grid('linest','none','tickdir','in','fontsize',12,'linewidth',2);
% colormap(mapa2);%brighten(0.5)
% caxis([-0.02 0.02])
% axis normal 
% %um = subplot('Position',[0.1900 0.5800 0.6550 0.3412]);%subplot(2,1,1)
% um = subplot('Position',[0.1900 0.5800 0.6412 0.3412]);%subplot(2,1,1)
% m_proj('miller','long',[-180 180],'lat',[-60 30]);
% m_contourf(oi.lon,oi.lat,C,-1:.1:1,'edgecolor','none')
% hold on
% m_scatter(X(ui),Y(ui),1,'k.')
% m_coast('color',[0 0 0],'linewidth',3);
% m_grid('linest','none','tickdir','out','fontsize',12,'linewidth',2)%,'box','fancy');
% colormap(mapa2);%brighten(0.5)
% caxis([-.601 .601])
% m_text(-175,18,'a.','backgroundcolor',[1 1 1],'fontweight','bold')
% col=colorbar('position',[.85 0.58 .02 .3412 ]);
% %col=colorbar('horizontal','position',[.2 .51 .6 .02]);
% set(col,'fontsize',12,'xtick',[-.5:.25:.5])
% axis normal
% %%m_line([28 28;15 15;15 28;28 15]',[-18 0;0 -18;0 0;-18 -18]','color',[0 0 0])
% %%slp_reg1 = nanmean(nanmean((ADT0(:,oi.lat>-18&oi.lat<0,oi.lon<28&oi.lon>15)),2),3);
% if nmod==1
% slp_reg = nanmean(nanmean((ADT0(:,oi.lat>-40&oi.lat<-22,oi.lon<-5&oi.lon>-25)),2),3);
% m_line([0 0;-25 -25;-25 0;0 -25]',[-40 -22;-22 -40;-22 -22;-40 -40]','color',[0 0 0])
% elseif nmod==2
% slp_reg = -nanmean(nanmean((ADT0(:,oi.lat>-35&oi.lat<-20,oi.lon<10&oi.lon>-10)),2),3);
% m_line([10 10;-10 -10;-10 10;10 -10]',[-35 -20;-20 -35;-20 -20;-35 -35]','color',[0 0 0])
% elseif nmod==3
%     llon = [-150 -180];
%     llat = [-20 0];   %original
%   %  llat = [-36 18];
%   %  llon = [-90 -150];
%   %  llat = [-60 -40];
% slp_reg = nanmean(nanmean((ADT0(:,oi.lat>llat(1)&oi.lat<llat(2),oi.lon<llon(1)&oi.lon>llon(2))),2),3);
% %slp_reg = -nanmean(nanmean((ADT0(:,oi.lat>-60&oi.lat<40,oi.lon<-90&oi.lon>-150)),2),3);
% %m_line([-160 -160;-180 -180;-180 -160;-160 -180]',[-20 0;0 -20;0 0;-20 -20]','color',[0 0 0])
% m_line([llon([1 1]);llon([2 2]);llon([2 1]);llon]',[llat;llat([2 1]);llat([2 2]);llat([1 1])]','color',[0 0 0])
% end
% R = corr(slp_reg,PC2');
% %slp_reg=slp_reg3;%-slp_reg1;
% um = subplot('Position',[0.19 0.0700 0.6412 0.18]);%figure
% plot(ygfrac,(PC2(:)-mean(PC2))/std(PC2))
% hold on
% plot(ygfrac,sign(R)*(slp_reg(:)-mean(slp_reg))/std(slp_reg))
% axis tight;ylim([-3.5 3.5])
% set(gca,'linewidth',1.5)
% text(1990,2.5,'c.','backgroundcolor',[1 1 1],'fontweight','bold')
% leg=legend(sprintf('%s%d','PC',nmod),'SLP','location','southwest');set(leg,'box','off','orientation','horizontal')%outside')
% %stop
% figname = sprintf('%s%d%s','CORR_SLP_MODE',nmod,'_b.jpg')
% %print('-dpng','-r400',figname)
% end %for nmod
% stop
end %if
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%HOVMOLLER
if 0
   [~,i22] = min(abs(lat+30));%22.5));%30));%22.5));
for ii = 1:N
    oo = zeros(nt,size(EOFs,2));
    for jj=ii:-1:1%:ii
   %oo = oo + transpose(PC(jj,:))*conj(EOFs(jj,:)); %ADD MODES
   oo = transpose(PC(ii,:))*conj(EOFs(ii,:));       %INDIVIDUAL MODES
    end
   reconstruction = mat2map(mask,real(oo));
   figure(10+ii)
   contourf(lon,ydec,squeeze(reconstruction(:,i22,:)),-10:.5:10);
   caxis([-5 5])
end
end

%Phase velocity
if 0
[area,dx,dy]=sw_area(lon,lat);
clear u v
for ii=1:N
    %taux  = squeeze(theta(ii,:,:))*pi/180;
    taux = squeeze(theta(ii,:,:));taux(taux<0)=taux(taux<0)+360;
[dthetax,dthetay] = gradient(taux);
%[dthetax2,dthetay2] = gradient(taux2);
%dthetax = min(dthetax,dthetax2)
dt = 7*24*60*60;
dd = diff(phi(ii,:));
idd = find(abs(dd)>180);
phi2 = phi(ii,:);
for jj=1:length(idd)
    phi2(idd(jj)+1:end)=phi2(idd(jj)+1:end)+360;
end
dphidt = gradient(phi2*pi/180,dt);
for jj=1:nt
u(:,:,ii,jj) = -1*dphidt(jj)./(dthetax./dx);
v(:,:,ii,jj) = -1*dphidt(jj)./(dthetay./dy);
end
end
end %0

%CHP = mat2map(mask,EOFs);
NLAG = 5        %6;%4;%DT = 2;%7;Dt = 7;
if 0
PHI = imag(PC(:,1:NLAG));
else
%FIND PHI'S for each 45degree (veci):
PHI = zeros(N,NLAG);
for ii = 1:N
    dd = diff(phi(ii,:));
    %idd = find(abs(dd)>180); %dd<-180
    idd = find(abs(dd)>200); %dd<-180
 %  veci = -1*[180 90  0 -90 -180];
    veci = -1*[180 135 90 45 0];
    odd = find(diff(idd)>5);%odd=odd(1)
    vec = interp1(phi(ii,idd(odd(1))+1:idd(odd(2))),idd(odd(1))+1:idd(odd(2)),veci,'nearest','extrap');
    if sum(abs(phi(ii,vec)-veci))>100
     vec = interp1(phi(ii,idd(odd(2))+1:idd(odd(3))),idd(odd(2))+1:idd(odd(3)),veci,'nearest','extrap');
    end
   PHI(ii,1:NLAG) = PC(ii,vec(1:NLAG));
end
end
%
%PHI = (PC(:,1:NLAG));
dispCEOF_aviso(lon,lat,ydec,CHP,EXPVAR,PHI,1:N,PC)  %PC timeseries to implement still
%lon,lat,CHP,PC,EXPVAR,Dt,NLAG,1:N)
%dispEEOF_aviso_figures(lon,lat,ydec,CHP,PC,EXPVAR,Dt*DT,NLAG,1:N)
stop

if 0
   % for ii=1:N
    save2eps(sprintf('%s%d%s%d%s%d%s','maps_interannual_',N,'modes_20_34.eps'),1,300)
 %    save2eps(sprintf('%s%d%s%d%s%d%s','CEOF',ii,'_LAG',NLAG,region,wave,'.eps'),ii,300)
 %   end
end

wave_global(real(PC(1,:)),time,7,1);
wave_global(real(PC(2,:)),time,7,0);
wave_global(real(PC(3,:)),time,7,0);
%wave_global(real(PC(4,:)),time,7,0);
h=get(gca);
legend(h.Children([6 4 2]),'PC1','PC2','PC3','location','northeast')
%print -djpeg -r300 wavelet_spectrum_interannual.jpg
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

