function [r,r0t,r0pt,dr00pp] = kergauss(mxddifxd,mxgdifxd,theta) 
% KERGAUSS, covariance KERnel, GAUSSian
%     subroutine to be called by SPAtial smoothing version of SiZer
%
%     Creates a series of images,
%     (as postscript files, or mpeg movies)
%     with family of smooths at top, colored SiZer map at bottom
%
%   Can use first 3, 4, 5, 6, 7, 8 or 9 arguments.
% Inputs:
%    mxddifxd  - matrix of pairwise differences, xdata - xdata
%          Computed in SpaSiZer1.m as:
%              mxddifxd = vec2mat(xdat,ndat) - vec2mat(xdat',ndat) ;
%    mxgdifxd  - matrix of pairwise differences, xgrid - xdata
%          Computed in SpaSiZer1.m as:
%              mxgdifxd = vec2mat(xgrid,ndat) - vec2mat(xdat',nxgrid) ;
% Outputs:
%           r  - matrix sig11 / rho 
%         r0t  - matrix sig21 / rho = transpose of R0 matrix
%                                   = transpose of matrix sig12 / rho
%     (Optional:)
%        r0pt  - matrix sig31 / rho = transpose of R0p matrix
%                                   = transpose of matrix sig13 / rho
%      dr00pp  - diagonal entries of matrix sig33 / rho =
%                                   = matrix R00pp




r = exp(-mxddifxd.^2 / (2*theta^2)) ;
          %  matrix sig11 / rho 


r0t = exp(-mxgdifxd.^2 / (2*theta^2)) ;
          %  matrix sig21 / rho = transpose of R0 matrix
          %                     = transpose of matrix sig12 / rho


if nargout > 2 ;

  r0pt = r0t .* (-mxgdifxd / theta^2) ;  
          %  matrix sig31 / rho = transpose of R0p matrix
          %                     = transpose of matrix sig13 / rho

  if nargout > 3 ;

    nxgrid = size(mxgdifxd,1) ;

    dr00pp = (1 / theta^2) * ones(nxgrid,1) ;
          %ds33  diagonal entries of sig33 matrix (suitably rescaled)
          %  (this matrix is actually exactly r00pp)

  end ;

end ;


