clear; close all
load Feb98_SST

utl_mfig;clf
utl_contourfill(lon,lat,sst,50);  colormap(utl_pmap(7)); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on; 
ax=caxis;
title 'Pacific SST Feb. 1998'


in=find(~isnan(sst));
lon1=lon(in);
lat1=lat(in);
y=sst(in);


% fit a plane
%==========================================================
%   y = Ex
%   J = (y-Ex)' * (y-Ex)
%==========================================================
utl_mfig;
E = zeros(length(in),3);

E(:,1)=lon1;
E(:,2)=lat1;
E(:,3)=1;


CME = E'*E;
[I,J]=size(CME);
cond = [1:I]*0 + 1.0e-6;
lambda=eig(CME);
CME = CME + diag(cond);
conditional_number = min(lambda) 
xhat = CME \ E'*y;

yhat=sst*nan;
yhat(in) = E*xhat;
clf
subplot(2,1,1)
utl_contourfill(lon,lat,yhat,50); colormap(utl_pmap(7)); 
caxis(ax); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on; 
title 'Mean SST = a + b*lon + c*lat     [ J = (y-Ex)'' (y-Ex) ]  '

subplot(2,1,2)
utl_contourfill(lon,lat,sst-yhat,50); colormap(utl_pmap(7)); 
caxis('auto'); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on; 
title ' y - yhat '




%==========================================================
%   y = Ex
%   J = (y-Ex)' *W* (y-Ex)   California
%==========================================================
utl_mfig
in2 = find ( lon > -250 & lon < -225 & lat <45 );
mask=sst; mask(~isnan(mask))=1;
area=sst*nan; area(in2)=1; area=area.*mask;

% South China Sea
%in2 = find ( lon1 > -250 & lon1 < -225 & lat1 <45 & lat1 > 25);

% California
in2 = find ( lon1 > -133 & lon1 < -90 & lat1 <50 );

%in2 = find (( lon1 > -250 & lon1 < -225 & lat1 <45 & lat1 > 25) ...
%             | ( lon1 > -133 & lon1 < -90 & lat1 <50 ) );


E = zeros(length(in),3);

E(:,1)=lon1;
E(:,2)=lat1;
E(:,3)=1;

W =zeros(length(in),1) + 0.00001; W(in2) = 1;
W = diag(W);

 
CME = E'*W*E;
[I,J]=size(CME);
cond = [1:I]*0 + 1.0e-6;
lambda=eig(CME);
CME = CME + diag(cond);
conditional_number = min(lambda) 
xhat = CME \ E'*W*y;

yhat=sst*nan;
%xhat = [ 36 0 0]';
yhat(in) = E*xhat;
clf
subplot(2,1,1)
utl_contourfill(lon,lat,yhat,50); colormap(utl_pmap(7)); 
caxis(ax); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on;
title 'Mean SST = a + b*lon + c*lat    [ J = (y-Ex)'' *W* (y-Ex) ]  W = California '

subplot(2,1,2)
utl_contourfill(lon,lat,sst-yhat,50); colormap(utl_pmap(7)); 
caxis('auto'); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on;
title ' y - yhat'


%==========================================================
%   y = Ex
%   J = (y-Ex)' *W* (y-Ex)   S. China
%==========================================================
utl_mfig
in2 = find ( lon > -250 & lon < -225 & lat <45 );
mask=sst; mask(~isnan(mask))=1;
area=sst*nan; area(in2)=1; area=area.*mask;

% South China Sea
in2 = find ( lon1 > -250 & lon1 < -225 & lat1 <45 & lat1 > 25);


E = zeros(length(in),3);

E(:,1)=lon1;
E(:,2)=lat1;
E(:,3)=1;

W =zeros(length(in),1) + 0.00001; W(in2) = 1;
W = diag(W);

 
CME = E'*W*E;
[I,J]=size(CME);
cond = [1:I]*0 + 1.0e-6;
lambda=eig(CME);
CME = CME + diag(cond);
conditional_number = min(lambda) 
xhat = CME \ E'*W*y;

yhat=sst*nan;
%xhat = [ 36 0 0]';
yhat(in) = E*xhat;
clf
subplot(2,1,1)
utl_contourfill(lon,lat,yhat,50); colormap(utl_pmap(7)); 
caxis(ax); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on;
title 'Mean SST = a + b*lon + c*lat    [ J = (y-Ex)'' *W* (y-Ex) ]  W = S. China '

subplot(2,1,2)
utl_contourfill(lon,lat,sst-yhat,50); colormap(utl_pmap(7)); 
caxis('auto'); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on;
title ' y - yhat'




%==========================================================
%   y = Ex
%   J = (y-Ex)' *W* (y-Ex) + x'W2x
%==========================================================
utl_mfig;
in2 = find ( lon > -250 & lon < -225 & lat <45 );
mask=sst; mask(~isnan(mask))=1;
area=sst*nan; area(in2)=1; area=area.*mask;

% South China Sea
in2 = find ( lon1 > -250 & lon1 < -225 & lat1 <45 & lat1 > 25);


%xhat_bad=[         0.463666967682766
%         -1.24393472594894
%          162.749383952116];
	    
%W2 = [ 1 .2 .2]'./xhat_bad;

%W2 = [72.992700729927
%         1.52695067949305
%        0.0250634732459955];
	  
%W2 = [ 0 31826950000000000000 0.0000]';
%W2 = [ 0 0 31826950000000000000 ]';

W2=[1000 1000 1000]';	    
W2=diag(W2);

E = zeros(length(in),3);

E(:,1)=lon1;
E(:,2)=lat1;
E(:,3)=1;

W =zeros(length(in),1) + 0.00001; W(in2) = 1;
W = diag(W);

 
CME = E'*W*E + W2;
[I,J]=size(CME);
cond = [1:I]*0 + 1.0e-6;
lambda=eig(CME);
CME = CME + diag(cond);
conditional_number = min(lambda) 
xhat = CME \ E'*W*y;

yhat=sst*nan;
%xhat = [ 36 0 0]';
yhat(in) = E*xhat;

clf
subplot(2,1,1)
utl_contourfill(lon,lat,yhat,50); colormap(utl_pmap(7)); 
caxis(ax); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on; 
title 'Mean SST = a + b*lon + c*lat    [ J = (y-Ex)'' *W* (y-Ex) + x'' W2 x]  W = S. China '

subplot(2,1,2)
utl_contourfill(lon,lat,sst-yhat,50); colormap(utl_pmap(7)); 
caxis('auto'); colorbar('h')
axis('equal'); set(gca,'ylim',[20 60],'xlim',[-250 -90]);
hold on; 
title 'Mean SST = a + b*lon + c*lat    [ J = (y-Ex)'' *W* (y-Ex) + x'' W2 x]  W = S. China '