clear all;

% load background field file
load background.dat;

% read data into individual arrays
xb_lat = background(:,1);
xb_lon = background(:,2);
xb_elev = background(:,3);
xb = background(:,4);

% set plotting range
xi = -117.0:0.05:-109.0;
yi = 39.0:0.05:44.0;
v = [-8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18];
w = [250 500 750 1000 1250 1500 1750 2000 2250 2500 2750 3000 3250 3500 3750];

[XI,YI] = meshgrid(xi,yi);

ZZI = griddata(xb_lon,xb_lat,xb,XI,YI);
contourf(XI,YI,ZZI,v);
colorbar;
hold;

ZI = griddata(xb_lon,xb_lat,xb_elev,XI,YI);
contour(XI,YI,ZI,w);

% load observation data
load merged.dat;
yo_lat = merged(:,1);
yo_lon = merged(:,2);
yo_elev = merged(:,3);
yo = merged(:,4);

% define observation error covariance matrix
sigo = 1;
for k=1:443
po(k,k)=sigo
end

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

