% j horel 04/09
% examples of different types of analysis solutions
% for temperature in the early morning extending from
% the Rush Valley to the Salt Lake Valley

% arbitrarily defined input matrices

% rg- distance in km from left edge of gridpoints
%     grid simply 101 km from left to right
% zg- terrain elevation (in m) of grid
% xt- truth at grid points: assumed arbitrary analytic solution 
% xto- truth at observation points
% xb- background. here assumed to be very poor 
% ro- distance in km from left edge of observation points
% zo- terrain elevation (in m) of observations

% other matrices

%p - polynomial order
%radb - radius of influence of barnes analysis

% d- observation innovation: difference between observations
%    and background at observation points
% rho- background error correlation matrix
% m- diagonal matrix for bratseth solution
% wba-  weights for barnes solution
% w- optimal weights for oi solution
% wo- weights at obs locations for bratseth
% wg- weights at grid locations for bratseth

% error matrices

% sigb- background error variance
% sigo- observation error variance
% sigma - ration of sigb to sigo
% h- linear observation operator
% po- observation error covariance matrix
% pb- background error covariance matrix
%     assumed to decorrelate as a function of gridpoint separation

% solutions

% xp- polynomial fit using observations
% xpc- polynomial fit using observation innovations
% xba- barnes analysis
% xoi- optimum interpolation solution
% xbr- bratseth iterative solution
% xbrm- bratseth matrix solution
% x3dvcg- 3dvar counter gradient solution

%grid 1km
%rush valley
for i = 1:31
rg(i) = i;
zg(i) = 200/225 * ( i - 16 )^2 + 1300;
xt(i) = 3/225 * ( i - 16 )^2 + 3;
xb(i) = 2 +4/800 * (2000-zg(i)) ;
end

%oquirrhs
for i = 31:51
rg(i) = i;
zg(i) =  2000 - 500/100 * ( i - 41 )^2;
xt(i) = 4/100 * ( i - 41 )^2 +  2;
xb(i) = 2 +4/800 * (2000-zg(i)) ;
end

%salt lake valley
for i = 51:101
rg(i) = i;
zg(i) = 250/625 * ( i - 76 )^2 + 1200;
xt(i) =  1/625 * ( i - 76 )^2 + 5.;
xb(i) = 2 +4/800 * (2000-zg(i)) ;
end


% arbitrarily defined observations
yo(1) = 4;
ro(1) = 5;
zo(1) = 1400.;
yo(2) = 3;
ro(2) = 15;
zo(2) = 1310.;
yo(3) = 3.5;
ro(3) = 20;
zo(3) = 1325.;
yo(4) = 5;
ro(4) = 25;
zo(4) = 1380.;
yo(5) = 6.1;
ro(5) = 31;
zo(5) = 1490;
yo(6) = 1.5;
ro(6) = 40;
zo(6) = 1990;
yo(7) = 5;
ro(7) = 75;
zo(7) = 1210;
yo(8) = 5.5;
ro(8) = 85;
zo(8) = 1222.;
yo(9) = 5.8;
ro(9) = 95 ;
zo(9) = 1334;
yo(10) = 6.1;
ro(10) = 100;
zo(10) = 1440;

% plot terrain
plot(rg,zg)
hold on
plot(ro,zo,'r*')
legend('Grid','Obs', 0)
pause
hold off
plot(rg,xt,'b-x')
hold on
plot(rg,xb,'g--o')
plot(ro,yo,'r*')
legend('Truth','Background','Obs', 0)

for k=1:10
d(k) = yo(k)-xb(ro(k));
end

% polynomial solution to observations
p = input('Enter the polynomial order>');
poly = polyfit(ro,yo,p);
xp = polyval(poly,rg);
plot(rg,xp,'y^')
qual_xp = 1/var(xp - xt)

% polynomial solution to observation corrections
polyj = polyfit(ro,d,p);
xpc = xb + polyval(polyj,rg);
qual_xpc = 1/var(xpc - xt)
plot(rg,xpc,'y--')
legend('Truth','Background','Obs','Poly O', 'Poly OC', 0)

pause
% barnes solution
radb = input('Enter the radius of influence for Barnes>');
% set initial values
for i = 1:101
xba(i) = xb(i);
end
for n = 1:10
for i = 1:101
for k = 1:10
distx = (rg(i) - ro(k))^2;
wba(k,i)=exp(-distx/radb^2);
xba(i) = xba(i) + wba(k,i) * ( yo(k) - xba(i));
end
end  
end
plot(rg,xba,'md-')
qual_xba = 1/var(xba - xt)
legend('Truth','Background','Obs','Poly O', 'Poly OC','Barnes', 0)
pause

%background error variance
sigb = var(xb - xt);
qual_b = 1/sigb

%observation error
for k = 1:10
xto(k) = xt(ro(k));
end

sigo = var(yo - xto);
qual_o = 1/sigo

sigma = sigo / sigb;


h = zeros(10,101);
for  k=1:10
h(k,ro(k)) = 1;
po(k,k)=sigo;
end

rad = input('Enter the background error decorrelation scale>');
for j=1:101
for i=1:101
r = i-j;
pb(i,j) =  sigb* exp (-r^2/rad^2);
end
end

% oi solution
a = po + h * pb * h';
ai = inv(a);
pbht = pb * h';
w = pbht * ai;
xoi = xb' + w * d' ;
plot(rg,xoi,'cd')
qual_xoi = 1/var(xoi - xt')
legend('Truth','Background','Obs','Poly O', 'Poly OC','Barnes','OI', 0)
pause

% bratseth solution
radbr = input('Enter the Bratseth radius of influence>');
% set initial values
for k=1:10
for l=1:10
r = ro(k) - ro(l);
rho(l,k) =  exp (-r^2/radbr^2);
end
end
m = sigma + sum(rho);
for k=1:10
for l=1:10
if k==l
wo(l,k) = (rho(l,k) + sigma )/m(k);
else
wo(l,k) = rho(l,k) /m(k);
end
end
end

for i=1:101
for k=1:10
r = rg(i) - ro(k);
wg(k,i) = exp (-r^2/rad^2)/m(k);
end
end

for k = 1:10
dbr(k) = d(k);
end

% iterate n times
for n = 1:10

for k = 1:10
db(k) = d(k);
for l = 1:10
db(k) = db(k) - wo(k,l)*dbr(l);
end
db(k) = db(k) + dbr(k);
dbr(k)=db(k);
end  
end

for i = 1:101
xbr(i) = xb(i);
for k = 1:10
xbr(i) = xbr(i) + wg(k,i) * db(k);
end
end

plot(rg,xbr,'ks-')
qual_xbr = 1/var(xbr - xt)

% bratseth matrix solution

mm = diag(m);
mi = inv(mm);

for i=1:101

xbrm = xb' + pbht * mi * d';
dbrm = (eye(10,10) - a * mi)*d' + d';
end
plot(rg,xbrm,'rs-')
qual_xbrm = 1/var(xbrm - xt')
legend('Truth','Background','Obs','Poly O', 'Poly OC','Barnes','OI','Bratseth','Bratseth Matrix', 0)

% 3Dvar conjugate gradient solution

pbhpo = pb' * h' * inv(po) * d';
pbt = pb' + pb' * h' * inv(po) * h * pb;
x3dvcg = xb' + pb * cgs(pbt,pbhpo,1e-6,100); 
plot(rg,x3dvcg,'gd')
qual_x3dvcg = 1/var(x3dvcg - xt')
legend('Truth','Background','Obs','Poly O', 'Poly OC','Barnes','OI','Bratseth','Bratseth Matrix','3DVAR',0)
pause
