function [X,Y,F]=DipoleGrid(source,field,Xin,Yin)
%
%   [X,Y,F]=DipoleGrid(source,field,swne)    grid is  1m by 1m
%                                         computed from swne corners
%                    or
%
%   [X,Y,Z]=DipoleGrid(source,field,Xin,Yin) grid is input as (meshgrid) Xin Yin
%
%
% swne(1) x coordinate sw corner of contour map
% swne(2) y coordinate sw corner 
% swne(3) x coordinate ne corner
% swne(4) y coordinate ne corner
% source(1) x cooordinate of south pole of source
% source(2) y coordinate of south pole
% source(3) depth to south pole
% source(4) length of dipole 
% source(5) inclination of polarization
% source(6) declination of polarization
% source(7) dipole moment
% field(1) inclination of field
% field(2) declination of field

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Dec, 1996)  
%                           (last - Nov, 2000
%***********************************************************************************
%
%  DipoleGrid graphs the total field anomaly from a dipole whose South pole is
%  at depth 'depth', and x,y=xsp,ysp. The length of the dipole is 'length'. 
%  The inclination of the Earth's field is 'incf' (deg)
%  and the declination of the Earth's field is  'decf'.
%  The inclination of polarization is 'incp', and the declination of 
%  polarization is 'decp'.
%  North is the top of the figure.
%  moment is the magnetic dipole moment in amp m^2 
%
%
clf
xsp=source(1);     % x cooordinate of south pole of source
ysp=source(2);     % y coordinate of south pole
depth=source(3);   % depth to south pole
length=source(4);  % length of dipole
incp=source(5);    % inclination of polarization
decp=source(6);    % declination of polarization
moment=source(7);  % dipole moment
incf=field(1);     % inclination of field
decf=field(2);     % declination of field

  if nargin==3     % if only three arguments then grid is computed from corners
  swne=Xin;
  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=[floor(swx):.1:floor(nex)];
  y=[floor(swy):.1:floor(ney)];
  [X,Y]=meshgrid(x,y); % set up X Y grids
  elseif nargin==4
  X=Xin;
  Y=Yin;
  end
  
xy=X+i*Y;
xsys=xsp+i*ysp; % coords of South pole
d1=depth; %       depth to South pole (+)
d2=depth+(length*sin(incp*pi/180));%    calc depth of North pole (-)
%
% calculate xy coords of north pole
%
xnyn=xsys+length*cos(incp*pi/180)*(sin(decp*pi/180)+i*cos(-decp*pi/180));
%
r1sq=(real(xy-xsys).^2+imag(xy-xsys).^2+d1^2);%  calc dist^2 from obs pt to xsys  
r2sq=(real(xy-xnyn).^2+imag(xy-xnyn).^2+d2^2);% calc dist^2 from obs pt to xnyn 
Fms=moment*ones(size(xy))./(r1sq);%    calc field strength due to South pole.
Fms=Fms*100/length; % converts moment to monopole strength and T to nT
Fmn=-moment*ones(size(xy))./(r2sq);%    calc field strength due to North pole
Fmn=Fmn*100/length; % converts  moment to monopole strength and T to nT
%Fxy=f1.*(xy./sqrt(r1sq))+f2.*(xy+length*cos(incp*pi/180))./sqrt(r2sq);
%Fms=f1./(r1sq)
%Fmn=f2./(r2sq);%  calc the mag of total field vect
%
%   calculate the total field strength in the direction of the inducing field
%
% x y z coords of survey point
%
vpx=real(xy);
vpy=imag(xy);
vpz=zeros(size(xy));
%
% x y z coords of south pole
%
vsx=real(xsys);
vsy=imag(xsys);
vsz=d1;
%
% direction cosines of vector from south pole to field point
%
vsy=(vsy-vpy)./sqrt(r1sq);
vsz=(vsz)./sqrt(r1sq);
vsx=(vsx-vpx)./sqrt(r1sq);
%
%
% x y z coords of north pole
%
vnx=real(xnyn);
vny= imag(xnyn);
vnz=(d2);
%
% direction cosine of vector from north pole to field point
%
vnx=(vnx-vpx)./sqrt(r2sq);
vny=(vny-vpy)./sqrt(r2sq);
vnz=(vnz)./sqrt(r2sq);
%
%  calculate the unit vector in the main field direction
%
vf=[cos(incf*pi/180)*sin(decf*pi/180) cos(incf*pi/180)*cos(decf*pi/180) sin(incf*pi/180)];
%
% project total field vector onto main field direction
%
Fxn=Fmn*vf(1).*vnx;
Fyn=Fmn*vf(2).*vny;
Fzn=Fmn*vf(3).*vnz;
Fxs=Fms*vf(1).*vsx;
Fys=Fms*vf(2).*vsy;
Fzs=Fms*vf(3).*vsz;
F=(Fxn+Fyn+Fzn+Fxs+Fys+Fzs); 
%F=fliplr(F);
hy=gca;
axis('off')
  axes('position',[.4 .1 .55 .8]); % axes for contour map
axis('square')
hx=gca;
contour(X,Y,F,40);
h=text(round(real(xsys)),round(imag(xsys))+1,'S');
set(h,'fontname','times','fontsize',16);
h=text(round(real(xnyn)),round(imag(xnyn))+1,'N');
set(h,'fontname','times','fontsize',16);
hold on
hl=plot([real(xnyn) real(xsys)],[imag(xnyn) imag(xsys)],'w-');
set(hl,'LineWidth',2);
%
% this section creates the table of model parameters
%
  axes(hy); % return to original axes for model table
  axis([0 10 0 10]);
  hold on
  h0=text(.9,9,'MODEL');,set(h0,'fontname','times','fontsize',14);
  h1=text(.5,8,'POLARIZATION');
  set(h1,'fontname','times','fontsize',12);
  h1=text(.2,7.5,'INCLINATION');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,7.5,num2str(floor(incp)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,7,'DECLINATION');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,7,num2str(floor(decp)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,6.5,'MAXIMUM FIELD');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,6.5,num2str(floor(max(max(F)))));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,6,'MINIMUM FIELD');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,6,num2str(floor(min(min(F)))));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,5.5,'MAXIMUM DEPTH');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,5.5,num2str(floor(d2)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,5,'MINIMUM DEPTH');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,5,num2str(floor(d1)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,4.5,'LENGTH');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,4.5,num2str(floor(length)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(1,3,'FIELD');
  set(h1,'fontname','times','fontsize',12);
  h1=text(.2,2.5,'INCLINATION');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,2.5,num2str(floor(incf)));
  set(h1,'fontname','times','fontsize',9);
  h1=text(.2,2,'DECLINATION');
  set(h1,'fontname','times','fontsize',9);
  h1=text(2,2,num2str(floor(decf)));
  set(h1,'fontname','times','fontsize',9);
  vector(9.5,8.5,1,(90-decf)*pi/180,0,0,0,'g-') % vector for dec
  boxtop=8.5; % top of box around table
  boxbot=1.5; % bottom of box around table
  boxleft=.05;% left side of box
  boxright=2.5;% right side of box
  hb=plot([boxleft boxright boxright boxleft boxleft],[boxbot boxbot boxtop ...
  boxtop boxbot],'w-');
  set(hb,'LineWidth',1.5)
  axis('off')
  hold on
%
%
axes(hx); % return to contour axes
xlabel('X'),ylabel('Y'),title('DIPOLE MODEL')
pretty