% DIpoleField Script makse calculations for DipoleFieldGui
% DipoleFieldScript is called by DipoleFieldGui


%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Dec, 1996)  
%                           (last - Jan, 2006
%***********************************************************************************
%

xsys=xsp+i*ysp; % coords of South pole
d1=depth; %       depth to South pole (+)
d2=depth+(dipolelength*sin(incp*pi/180));%    calc depth of North pole (-)
if (d2 < 0);disp(' '),disp('***************WARNING!!! ***************'),
    disp('NORTH POLE OF DIPOLE IS ABOVE GROUND'),
    disp('DECREASE LENGTH OF DIPOLE'),disp('          OR'),
    disp('CHANGE INCLINATION OF DIPOLE')
    disp(' '),disp('*****************************************'),
end
%
% calculate xy coords of north pole
%
xnyn=xsys+dipolelength*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/dipolelength; % 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/dipolelength; % converts  moment to monopole strength and T to nT
%   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); %TOTAL FIELD
Fx=F*cos(decf)*cos(incf);
Fy=F*sin(decf)*cos(incf);
Fz=F*sin(incf);
if Component==1;F=F;
elseif Component==2; F=Fxs+Fxn; % X component
elseif Component==3; F=Fys+Fyn; % Y component
elseif Component==4; F=Fzs+Fzn; % Z component
end
 axes(haxContour),cla
 contour(X,Y,F,40);hcolorbar=colorbar('vert');
htext=text(round(real(xsys)),round(imag(xsys))+1,'S','color','r');
set(htext,'fontname','times','fontsize',14);
htext=text(round(real(xnyn)),round(imag(xnyn))+1,'N','color','b');
set(htext,'fontname','times','fontsize',14);
hold on

harrow=plot([real(xnyn) real(xsys)],[imag(xnyn) imag(xsys)],'k-');
set(harrow,'LineWidth',2);
axes(haxContour); % return to contour axes
xlabel('E (grid units)'),ylabel('N (grid units)'),

if Component==1;title(' TOTAL FIELD DIPOLE MODEL');
elseif Component==2; title('X COMPONENT DIPOLE MODEL'); % X component
elseif Component==3; title('Y COMPONENT DIPOLE MODEL'); % Y component
elseif Component==4; title('Z COMPONENT DIPOLE MODEL'); % Z component
end
set(get(gca,'xlabel'),'fontname','times','fontsize',10);
set(get(gca,'ylabel'),'fontname','times','fontsize',10);
set(get(gca,'zlabel'),'fontname','times','fontsize',10);
set(get(gca,'title'),'fontname','times','fontsize',12);

ProfileHandle=plot([0 0],[-19 19],'-','LineWidth',2); % plot profile on grid
StartTextHandle1=text(1,-18,'START'); % print S(tart) of profile
EndTextHandle1=text(1,18,'END'); % print E(nd) of profile  


axes(haxmain)
%haxProfile=axes('position',[.1 .09 .5 .21]);
axes(haxProfile),cla

plot(Y(:,21),F(:,21),'.-')
xlabel('DISTANCE (grid units)'),ylabel('nT'),title('SOUTH TO NORTH PROFILE')
StartTextHandle=text(-19.5,max(F(:,21))/10,'START');          % print S(tart) of profile
EndTextHandle=text(16.8,max(F(:,21))/10,'END');  % print E(nd) of profile  

set(get(gca,'xlabel'),'fontname','times','fontsize',10);
set(get(gca,'ylabel'),'fontname','times','fontsize',10);
set(get(gca,'zlabel'),'fontname','times','fontsize',10);
set(get(gca,'title'),'fontname','times','fontsize',12);