clear all;

% read entire file
load xbfile_full.dat;

% read data into individual arrays
xb_lat_full = xbfile_full(:,1);
xb_lon_full = xbfile_full(:,2);
xb_elev_full = xbfile_full(:,3);
xb_full = xbfile_full(:,4);
% set plotting range
xi_full = -117.0:0.05:-109.0;
yi_full = 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_full,YI_full] = meshgrid(xi_full,yi_full);

% plot terrain over entire domain
ZI_full = griddata(xb_lon_full,xb_lat_full,xb_elev_full,XI_full,YI_full);
contour(XI_full,YI_full,ZI_full,w);
pause;

% plot background temperature over entire domain
ZZI_full = griddata(xb_lon_full,xb_lat_full,xb_full,XI_full,YI_full);
contourf(XI_full,YI_full,ZZI_full,v);
colorbar;
pause;

%read truncated file
load xbfile.dat;

% read data into individual arrays
xb_lat = xbfile(:,1);
xb_lon = xbfile(:,2);
xb_elev = xbfile(:,3);
xb = xbfile(:,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;
pause;

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

% load observation data
[yo_lat,yo_lon,yo_elev,yo_elev2,yo,yo_nam,yo_nam2] = textread('merged.dat','%f %f %f %f %f %s %s');

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

% define background error covariance matrix

sigb = ?????;
                                                                                    
radius = 6370;
rad = input('Enter the background error decorrelation scale>');
rad2 = rad * rad;
pb=zeros(4284,4284);
for j=1:4284
for i=1:4284
dx = radius * cos(pi * xb_lat(i)/180.)*pi*(xb_lon(i)-xb_lon(j))/180.; 
dy = radius * pi * (xb_lat(i) - xb_lat(j))/180.;
r2 = dx*dx + dy*dy;
if(r2<5*rad2)
pb(???) =  ????
end
end
end

% define forward transform matrix

h = zeros(443,4284);
for  k=1:443
min = 6370;
loc=1;
for i=1:4284
dx = radius * cos(pi * yo_lat(k)/180.)*pi*(yo_lon(k)-xb_lon(i))/180. ;
dy = radius * pi * (yo_lat(k) - xb_lat(i))/180.;
dist(i) = sqrt(dx*dx + dy*dy);
if(dist(i)<min)
min=dist(i);
loc=i;
end
end
h(k,loc) = 1;
end

% define the observation innovation matrix
d = ????

% 3Dvar conjugate gradient solution
                                                                                
x3dv = xb + pb * cgs(???,?????,.001,200);

hold off;
ZZV = griddata(xb_lon,xb_lat,x3dv,XI,YI);
contourf(XI,YI,ZZV,v);
colorbar;
