%callstart		% Cleans up all variables and figure windows
inst = ['   To make a hardcopy of any of the current figures type:'
        '        figure(#)     (# is the desired figure number)   '
        '        print         (send figure to default printer)   '
        '                                                         '
        '   To continue type: r-e-t-u-r-n                         '];

inst2= ['   At this point you may loop back to change mapping     '
        '   parameters without changing selected data points by   '
        '   entering a "1" below.                                 '
        '       If you enter a number other than 1, you will      '
        '       have the option of restarting "oa" to select      '
        '       different data points or exiting Matlab           '
        '           (note: you may also print hardcopies          '
        '            at this stage before exiting)                ']; 

x=-10:0.5:10;
y=x;
rtest=0:0.5:10;

% Get user name for plot identification
% hardwire height and width of analytical function for now
namstr=input('Enter your name for plot labels.............','s');
height=1.0;
width=20.0;

x0=0;
y0=0;
[X,Y]=meshgrid(x,y);
R=(X .^ 2 + Y .^ 2);
T = height .* exp(-R/width);
%T = height*(cos(2*pi*(X/width)) + cos(2*pi*(Y/width)));
[TM,TN]=size(T);
T_est = NaN .* zeros(TM,TN);
T_err = NaN .* zeros(TM,TN);

figure(1)
mesh(X,Y,T);
xlabel('X Location')
ylabel('Y Location')
title(['Analytical Function (' namstr ')'])

figure(2)
plot(X,Y,'g.')
set(gca,'aspectratio',[1,1])
xlabel('X Location')
ylabel('Y Location')
title(['Plan View with Obs Points (' namstr ')'])
hold on
cc=contour(X,Y,T);
clabel(cc)
hold off
%%plot(X,Y,'g+')

disp('  ')
disp('Left button for new observation, right button for last observation');
button=1;
xb=[];
yb=[];

while (button ~=3),

   [xx,yy,button]=ginput(1);
   
   xb=[xb;xx;];
   yb=[yb;yy];
   
   line('XData',xx,'YData',yy,'LineStyle','*','Color','r'); drawnow
   
end;

%loop here to select different mapping parameters using the
%same data points (restart to select new data points)
iloop=1;
while iloop == 1,

% Choose exponential or gaussian covariance function and
% Length scale parameter
Ctyp = input('Enter 1 (for Exp) or 2 (for Gaussian) Cov...');
L    = input('Enter covariance function e-folding scale...');
sl=int2str(L);
if Ctyp == 1
  disp(['Exponential Covariance Function Selected with L=' sl])
else
  disp(['Gaussian Covariance Function Selected with L=' sl])
end
amp = 1; % use unit amplitude for now (should be variance)

disp('Calculating Estimates........................')

% Calculate the distances between all observations.
obslen=length(xb);
s = zeros(obslen,obslen);
for j = 1:obslen
    s(j,:) = sqrt( (xb-xb(j)) .^2 + (yb-yb(j)) .^2 )'; 
end

% Calculate Cdd and Cddi and Cxx outside loop.

if Ctyp == 1
  Cdd = amp*exp(-s ./ L);
  Crr = amp*exp(-rtest ./ L);
else
  Cdd = amp*exp(-(s.^2) ./ L);
  Crr = amp*exp(-(rtest.^2) ./ L);
end
Cddi = inv(Cdd);
d = interp2(X,Y,T,xb,yb);

for xi = 1:1:41
    for yi = 1:1:41
    
% xs = sqrt( (xb - X(xi)) .^2 + (yb - Y(yi)) .^ 2 );  %%%ERROR LINE
        xs = sqrt( (xb - X(xi,yi)) .^2 + (yb - Y(xi,yi)) .^ 2 );
        xs = xs';

    if Ctyp == 1   
        Cxd = amp*exp(-xs ./ L);
    else
        Cxd = amp*exp(-(xs.^2) ./L);
    end
        A = Cxd * Cddi;
        T_est(xi,yi) = A * d;
		T_err(xi,yi) = 1 - A * Cxd';
        
%        if  rem(yi,10) == 0  &  rem(xi,10)  == 0
%           figure(3)
%           plot(Cxd)
%           disp([' ',int2str(X(xi,yi)),' ',int2str(Y(xi,yi))])

%           keyboard
%        end
        
    end
end

figure(3)
subplot(2,2,1)
cc=contour(X,Y,T);
clabel(cc)
title(['Plan View with Obs Points (' namstr ')'])
hold on
%%%%plot(X,Y,'g.')
xlabel('X Location')
ylabel('Y Location')
hold off
line('XData',xb,'YData',yb,'LineStyle','*','Color','r');
set(gca,'aspectratio',[1,1]);


subplot(2,2,2),plot(rtest,Crr)
xlabel('Separation Distance')
ylabel('Covariance Function')
if Ctyp ==1
  title(['Exp Cov Function, L=', sl ' (' namstr ')'])
else
  title(['Gauss Cov Function, L=', sl ' (' namstr ')'])
end

subplot(2,2,3)
mesh(X,Y,T_est);
xlabel('X Location')
ylabel('Y Location')
if Ctyp ==1
  title(['OA Map, Exp,L=' sl ' (' namstr ')'])
else
  title(['OA Map, Gauss,L=' sl ' (' namstr ')'])
end

subplot(2,2,4)
mesh(X,Y,T_err);
xlabel('X Location')
ylabel('Y Location')
if Ctyp ==1
  title(['OA Error, Exp,L=' sl ' (' namstr ')'])
else
  title(['OA Error, Gauss,L=' sl ' (' namstr ')'])
end

orient tall
% wysiwyg


% Compute rms difference between estimate and analytical fields
% as a measure of success
rmserr = sqrt(mean(mean((T-T_est).^2)));
trms = num2str(rmserr);
disp(['Rms difference between true (Fig1) & mapped (Fig3c) is: ' trms])

disp(' ')
disp(inst)
keyboard
disp(inst2)
iloop=input('Enter 1 (to select new mapping parameters)..');
end  %end of while loop

disp(' ')
disp('Enter "[xb yb]" to see listing of current obs points')
disp('Enter "oa" to restart and select new obs points')
disp('Enter e-x-i-t to end Matlab session')

