function [Xout,Yout,Zout]=PointMass(depth,gmax,xo,yo,Xin,Yin,deltax,deltay)
%
%
% creates gridded values of the gravity anomaly due to a point mass at depth
%  (in grid units) and position xo yo. The peak value is gmax.
%
% Xout -
% Yout - 
% Zout - output grid of values 
% Xin - grid
% Yin - grid
% swne - four element vector with xy coordinates of south west and north east corners
%        of grid to be created.
% xo - x position of point mass
% yo - y position of point mass
% depth - depth of point mass
% gmax - maximum value of gravity anomaly (mass is adjusted to force this)
%
%    [Xout,Yout,Zout]=PointMass(depth,gmax,xo,yo,Xin,Yin)
%
%                       or
%
%    [Xout,Yout,Zout]=PointMass(depth,gmax,xo,yo,swne,deltax,deltay)

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (June, 1995)  
%                              (last Nov 1998)
%***********************************************************************************

  if nargin==7     % if only five arguments then grid is computed from corners
  swne=Xin;
  deltay=deltax;
  deltax=Yin;
  swx=swne(1);       % x coordinate sw corner of contour map  
  swy=swne(2);       % y coordinate sw corner
  nex=swne(3);       % x coordinate ne corner
  ney=swne(4);       % y coordinate ne corner
  x=[swx:deltax:nex];
  y=[swy:deltay:ney];
  [Xout,Yout]=meshgrid(x,y); % set up X Y grids
  elseif nargin==6
  Xout=Xin;
  Yout=Yin;
  end

Zout=gmax*depth^3*ones(size(Xout))./((Xout-xo).^2+(Yout-yo).^2+depth^2).^(3/2);