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;

% hardwire height and width of analytical function for now
height=1.0;
width=20.0;

x0=0;
y0=0;
[X,Y] = meshgrid(x,y);
R=(X .^ 2 + Y .^ 2);
% pick an analytical function:
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 ')
pause(2)
%%
figure(2)
plot(X,Y,'g.')
set(gca,'aspectratio',[1,1])
xlabel('X Location')
ylabel('Y Location')
title('Plan View with Observations')
hold on
cc=contour(X,Y,T);
clabel(cc)
hold off

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

%% The button loop
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 to pick different mapping techniques.
%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 or de-correlation 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. Please wait.')
    
    %% Calculate the distances between all observations.
    obslen = length(xb); %number of buttonpresses
    s = zeros(obslen,obslen);
    for j = 1:obslen
        s(j,:) = sqrt( (xb-xb(j)) .^2 + (yb-yb(j)) .^2 )';
    end
    
    %%
    % Calculate Cdd and Cxx outside loop.
    
    if Ctyp == 1 % Exponential
        Cdd = amp*exp(-s ./ L);
        Crr = amp*exp(-rtest ./ L);
    else % Gauss
        Cdd = amp*exp(-(s.^2) ./ L); % distance
        Crr = amp*exp(-(rtest.^2) ./ L); % r
    end
    %Cddinv = inv(Cdd); instead use left devide
    % interpolates the analytical function onto 'observation'
    d = interp2(X,Y,T,xb,yb);
    
    %%
    for xi = 1:1:length(x)
        for yi = 1:1:length(y)
            
            % xs = sqrt( (xb - X(xi)) .^2 + (yb - Y(yi)) .^ 2 );  %%%ERROR LINE
            % Mixed covariance
            xs = sqrt( (xb - X(xi,yi)) .^2 + (yb - Y(xi,yi)) .^ 2 );
            xs = xs';
            
            if Ctyp == 1  % Exponential
                Cxd = amp*exp(-xs ./ L);
            else % Gauss
                Cxd = amp*exp(-(xs.^2) ./L);
            end
            
            T_est(xi,yi) = (Cxd \ Cdd) * d;
            T_err(xi,yi) = 1 - (Cxd \ Cdd )* Cxd';
            
        end
    end
    %%
    figure(3)
    subplot(2,2,1)
    cc=contour(X,Y,T);
    clabel(cc)
    title('Plan View with Obs.')
    hold on
    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 ''])
    else
        title(['Gauss Cov Function, L=', sl ''])
    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 ''])
    else
        title(['OA Map, Gauss,L=' sl ''])
    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 ''])
    else
        title(['OA Error, Gauss,L=' sl ''])
    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')