%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Program to perform a simple Cressman analysis
%
%  A.S. Lawless   14/01/02
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% User options
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_pts=100;   
  % Number of grid points
obs_position = [3,15,18,21,43,44,45,74,75,76,77,78,81,82,95]; 
  % Grid points at which there are observations
obs_sign=[1]; 
  % Position in observation vector of obs which have wrong sign
obs_big=[5]; 
  % Position in observation vector of obs which are too large
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of user options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Close any open figure windows
h=gcf;
for i=1:h
  close(i)
end
%
xb=zeros(n_pts,1);
xa=zeros(n_pts,1);
xt=zeros(n_pts,1);
%
%%%%%%%%%%%%
% Generate truth and obs
delta_x=2*pi/(n_pts-1);
for i=1:n_pts
  xt(i)=sin((i-1)*delta_x);
  xb(i)=0.5*xt(i);
end
;
n_obs=length(obs_position);
obs=zeros(n_obs,1);
for k=1:n_obs
  obs(k) = xt(obs_position(k));
end
%
% Make some obs incorrect
n_wrong=length(obs_sign);
for k=1:n_wrong
  obs(obs_sign(k))=-obs(obs_sign(k));
end
n_big=length(obs_big);
for k=1:n_big
  obs(obs_big(k))=obs(obs_big(k))*1.5;
end
%
%%%%%%%%%%%%
% Plot truth and obs
h=plot(obs_position,obs,'o','LineWidth',2);
set(gca,'YLim',[-1.1 1.1])
set(gca,'XLim',[1 n_pts])
hold on
plot(xt,'k','LineWidth',1.5);
hold on
plot(xb,'b','LineWidth',1.5);
hold on
xlabel('Distance')
ylabel('Temperature')
title('Plot of temperature using Cressman analysis')
legend('Observations','Truth','Background')
%
%%%%%%%%%%%%
% Decide on weighting function
l_w=input('Which weighting function? (1=Square, 2=Exponential) ');
if (l_w~=1 & l_w~=2)
  error('Incorrect weighting function requested')
end
%
% Set up weighting function
R=input('Radius of influence? ');
Rsq=R*R;
w=zeros(n_pts,n_obs);
if l_w == 1  % Square weighting function
  for k=1:n_obs
    for i=1:n_pts
      d=abs(obs_position(k)-i);
      dsq=d*d;
      w(i,k) = (Rsq-dsq)/(Rsq+dsq);
      if (w(i,k)<0)
        w(i,k) = 0.;
      end
    end
  end
else % Exponential weighting function
  for k=1:n_obs
    for i=1:n_pts
      d=abs(obs_position(k)-i);
      dsq=d*d;
      w(i,k) = exp(-dsq/(2*Rsq));
    end
  end
end
%
%%%%%%%%%%%%
% Perform analysis
for i=1:n_pts
  delta_x = 0.;
  scalar = 0.;
  for k=1:n_obs
    i_ob = obs_position(k);
    delta_x = delta_x + w(i,k) * (obs(k)-xb(i_ob));
    scalar = scalar + w(i,k);
  end
  if scalar==0
    delta_x = 0.;
  else
    delta_x = delta_x / scalar;
  end
  xa(i) = xb(i) + delta_x;
end
%%%%%%%%%%%%
% Plot analysis
plot(xa,'r','LineWidth',1.5);
hold off
%
legend('Observations','Truth','Background','Analysis')