% 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

%grid 1km
%rush valley, oquirrhs, salt lake valley

vs = input('Enter contribution ratio of vert distance to horiz distance>'); % vertical scaling factor
i = 1:101;
rg = i;
zg = (i >= 0 & i <= 30).*(200/225 * ( i - 16 ).^2 + 1300) ...
    +(i >=31 & i <= 50).*(2000 - 500/100 * ( i - 41 ).^2) ...
    +(i >=51 & i <=101).*(250/625 * ( i - 76 ).^2 + 1200);

xt = (i >= 0 & i <= 30).*(3/225 * ( i - 16 ).^2 + 3 ) ...
    +(i >=31 & i <= 50).*(4/100 * ( i - 41 ).^2 + 2 ) ...
    +(i >=51 & i <=101).*(1/625 * ( i - 76 ).^2 + 5.);

xb = 2 +4/800 * (2000-zg);

% arbitrarily defined observations
yo = [   4   3   3.5    5  6.1  1.5    5  5.5  5.8  6.1];
ro = [   5  15    20   25   31   40   75   85   95  100];
zo = [1400 1310 1325 1380 1490 1990 1210 1222 1334 1440];
oscale = input('Enter the observation error scale>');
yo = xt(ro) + (yo - xt(ro))*oscale;

d = yo-xb(ro);

% polynomial solution to observations
p = input('Enter the polynomial order>');
poly = polyfit(ro,yo,p);
xp = polyval(poly,rg);

% polynomial solution to observation corrections
polyj = polyfit(ro,d,p);
xpc = xb + polyval(polyj,rg);

% barnes solution
radb = input('Enter the radius of influence for Barnes>');
xba = xb;
for n = 1:10
    for i = 1:101        
        distx = rg(i) - ro;
        distz = zg(i) - zo;
        wba =exp(-(distx/radb).^2-(distz/radb*vs).^2);
        for k = 1:10
            xba(i) = xba(i) + wba(k) * (yo(k) - xba(i));
        end
    end 
end

sigb = var(xb - xt);        %background error variance
sigo = var(yo - xt(ro));    %observation error
sigma = sigo / sigb;

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

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

% oi solution
a = po + h * pb * h';
ai = inv(a);
pbht = pb * h';
w = pbht * ai;
xoi = xb' + w * d' ;

% bratseth solution
radbr = input('Enter the Bratseth radius of influence>');
for k=1:10
    r = ro(k) - ro;
    rz = zo(k) - zo;
    rho(:,k) = exp(-(r/radbr).^2-(rz/radbr*vs).^2)';
end

m = sigma + sum(rho);
for k=1:10
    wo(:,k) = (rho(:,k) + sigma*(k==1) )/m(k);
end

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

dbr = d;

for n = 1:10
    db = d - dbr*wo' + dbr;
    dbr = db;
end

xbr = xb;
xbr = xbr + db*wg;

% bratseth matrix solution
mm = diag(m);
mi = inv(mm);

dbrm = d';
for j=1:10
    dbrm = (eye(10,10) - a * mi)*dbrm + dbrm;
end
xbrm = xb' + pbht * mi * dbrm;

%3DVar
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,zg)           % elevations
% plot(ro,zo,'r*')      
plot(rg,xt,'b-x')       %Truth
hold on
plot(rg,xb,'g--o')      %Background
plot(ro,yo,'r*')        %Obs
%plot(rg,xp,'y^')        %Poly O
%plot(rg,xpc,'b--')      %Poly OC
plot(rg,xba,'md-')      %Barnes 
plot(rg,xoi,'cd')       %OI
plot(rg,xbr,'ks-')      %Bratseth
plot(rg,xbrm,'rs-')     %Bratseth Matrix
plot(rg,x3dvcg,'gd')    %3D Var solution

legend('Truth','Background','Obs','Barnes','OI','Bratseth','Bratseth Matrix', '3DVar', 0)
% legend('Truth','Background','Obs','Poly O', 'Poly OC','Barnes','OI','Bratseth','Bratseth Matrix', 0)
hold off
fprintf(1,'Analysis quality: %f\t%f\t%f\t%f\t%f\n', ...
    1/var(xt - xba  ), ...
    1/var(xt - xoi' ), ...
    1/var(xt - xbr  ), ...
    1/var(xt - xbrm'), ...
    1/var(xt - x3dvcg'));