%
%  DipoleFieldGui graphs the magnetic anomaly of a dipole as
%  a contour map as well as a profile
%  dipole and field parameters are controlled via a GUI
%  default parameters are set and the default profile is South to North
%  The total field as well as X,Y and Z components can be displayed.
%  DipoleFieldGui calls DipoleFieldScript and DipoleFieldProfile

%***********************************************************************************
%
%                              Jim Merriam
%                       University of Saskatchewan
%                          jim.merriam@usask.ca
%                              (Nov, 2005)  
%                           (last - Jan, 2006
%***********************************************************************************
fig1=figure(1);clf,set(gcf,'position',[440 350 700 450]),
set(fig1,'name','                          INTERACTIVE MAGNETIC DIPOLE DEMONSTRATION')
clear i
haxmain=axes('position',[.05 .05 .92 .92],'visible','off');
swne=[-20 -20 20 20];
source=[0 0 1 4 0 15 100];

lat=52; 	% default latitude (Saskatoon)
field=[atan(2*tan(lat*pi/180))*180/pi 15];	% default field inclination

Component='T'; %default component is T Total Field
field(1)=atan(2*tan(lat*pi/180))*180/pi;


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
dipolelength=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

% uicontrol for depth
GuiColor=[.8 .8 .8];
GuiColor1=[.8 0 0];
GuiColor2=[0 .8 0];
hpopDepthValue=uicontrol('style','text','string',num2str(depth),'units','normalized', ...
    'position',[.93 .94 .039 .03],'backgroundcolor',GuiColor);
hpopDepthText=uicontrol('style','text','string','DEPTH TO SOUTH POLE','units','normalized', ...
    'position',[.725 .94 .205 .03],'backgroundcolor',GuiColor);
CALLBACKDepth=['depth=get(hpopDepth,''value'');set(hpopDepthValue,''string'',num2str(depth)); DipoleFieldScript'];
hpopDepth=uicontrol('style','slider','units','normalized','position',[.97 .88 .018 .115],'max',10,'min',0,'sliderstep',[.1 .2],'Callback','eval(CALLBACKDepth)');

yshift=.13;
hpopLengthValue=uicontrol('style','text','string',num2str(dipolelength),'units','normalized', ...
    'position',[.93 .94-yshift .039 .03],'backgroundcolor',GuiColor);
hpopLengthText=uicontrol('style','text','string','LENGTH OF DIPOLE','units','normalized', ...
    'position',[.73 .94-yshift .205 .03],'backgroundcolor',GuiColor);
CALLBACKLength=['dipolelength=get(hpopLength,''value'');set(hpopLengthValue,''string'',num2str(dipolelength)); DipoleFieldScript'];
hpopLength=uicontrol('style','slider','units','normalized','position',[.97 .88-yshift .018 .115],'max',10,'min',0,'sliderstep',[.1 .2],'Callback','eval(CALLBACKLength)');

hpopIncValue=uicontrol('style','text','string',num2str(incp),'units','normalized', ...
   'position',[.93 .94-yshift*2 .039 .03],'backgroundcolor',GuiColor);
hpopIncText=uicontrol('style','text','string','INCLINATION OF DIPOLE','units','normalized', ...
    'position',[.73 .94-yshift*2 .205 .03],'backgroundcolor',GuiColor);
CALLBACKInc=['incp=get(hpopInc,''value'');set(hpopIncValue,''string'',num2str(incp)); DipoleFieldScript'];
hpopInc=uicontrol('style','slider','units','normalized','position',[.97 .88-yshift*2 .018 .115],'max',90,'min',-90,'sliderstep',[5/180 .2],'Callback','eval(CALLBACKInc)');

hpopDecValue=uicontrol('style','text','string',num2str(decp),'units','normalized', ...
   'position',[.93 .94-yshift*3 .039 .03],'backgroundcolor',GuiColor);
hpopDecText=uicontrol('style','text','string','DECLINATION OF DIPOLE','units','normalized', ...
    'position',[.73 .94-yshift*3 .205 .03],'backgroundcolor',GuiColor);
CALLBACKDec=['decp=get(hpopDec,''value'');set(hpopDecValue,''string'',num2str(decp)); DipoleFieldScript'];
hpopDec=uicontrol('style','slider','units','normalized','position',[.97 .88-yshift*3 .018 .115],'max',360,'min',0,'sliderstep',[10/360 .1],'Callback','eval(CALLBACKDec)');

hpopLatValue=uicontrol('style','text','string',num2str(lat),'units','normalized', ...
   'position',[.93 .94-yshift*4 .039 .03],'backgroundcolor',GuiColor);
hpopLatText=uicontrol('style','text','string','LATITUDE','units','normalized', ...
    'position',[.73 .94-yshift*4 .205 .03],'backgroundcolor',GuiColor);

hpopIncfValue=uicontrol('style','text','string',num2str(floor(incf)),'units','normalized', ...
   'position',[.93 .94-yshift*4-.04 .039 .03],'backgroundcolor',GuiColor);
hpopIncfText=uicontrol('style','text','string','FIELD INCLINATION','units','normalized', ...
    'position',[.73 .94-yshift*4-.04 .205 .03],'backgroundcolor',GuiColor);
CALLBACKLat=['lat=get(hpopLat,''value'');set(hpopLatValue,''string'',num2str(lat));' ...
 'incf=atan(2*tan(lat*pi/180))*180/pi;set(hpopIncfValue,''string'',num2str(floor(incf))),DipoleFieldScript'];
hpopLat=uicontrol('style','slider','units','normalized','position',[.97 .88-yshift*4 .018 .115],'max',90,'min',-90,'sliderstep',[5/180 .1],'Callback','eval(CALLBACKLat)');

hpopDecfValue=uicontrol('style','text','string',num2str(decf),'units','normalized', ...
   'position',[.93 .94-yshift*5 .039 .03],'backgroundcolor',GuiColor);
hpopDecfText=uicontrol('style','text','string','FIELD DECLINATION','units','normalized', ...
    'position',[.73 .94-yshift*5 .205 .03],'backgroundcolor',GuiColor);
CALLBACKDecf=['decf=get(hpopDecf,''value'');set(hpopDecfValue,''string'',num2str(decf));field(2)=decf;DipoleFieldScript'];
hpopDecf=uicontrol('style','slider','units','normalized','position',[.97 .88-yshift*5 .018 .115],'max',360,'min',0,'sliderstep',[10/360 .1],'Callback','eval(CALLBACKDecf)');

hpopCompText2=uicontrol('style','text','string','TOTAL FIELD (T) EAST (X)','units','normalized', ...
    'position',[.73 .94-yshift*5.84 .205 .03],'backgroundcolor',GuiColor);
hpopCompText3=uicontrol('style','text','string','NORTH (Y) VERTICAL (Z)','units','normalized', ...
    'position',[.73 .94-yshift*6.04 .205 .03],'backgroundcolor',GuiColor);

CALLBACKComp=['Component=get(hpopComp,''value'');DipoleFieldScript'];
hpopComp=uicontrol('style','popup','STRING','T|X|Y|Z','units','normalized','position',[.94 .88-yshift*6 .05 .11],'Callback','eval(CALLBACKComp)');

hpopProfileText=uicontrol('style','text','string','EXTRACT A PROFILE?','units','normalized', ...
    'position',[.73 .94-7*yshift .24 .03],'backgroundcolor',GuiColor);
CALLBACKProfile=['DipoleFieldProfile'];
hpopProfile=uicontrol('style','radio','units','normalized','position',[.967 .94-yshift*7 .022 .03], ...
'Callback','eval(CALLBACKProfile)');



%  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):floor(nex)];
  y=[floor(swy):floor(ney)];
  [X,Y]=meshgrid(x,y); % set up X Y grids
  
xy=X+i*Y;
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 (-)
%
% 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

%F=fliplr(F);
axis('off')
 haxContour=axes('position',[.1 .45 .55 .5]); % axes for contour map
axis('square')
contour(X,Y,F,40);hcolorbar=colorbar('vert');
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',14);

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)'),title('DIPOLE MODEL')
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  


haxProfile=axes('position',[.1 .09 .48 .21]);

plot(Y(:,21),F(:,21),'.-')

StartTextHandle=text(-19.5,500,'START');          % print S(tart) of profile
EndTextHandle=text(16.8,500,'END');  % print E(nd) of profile  


xlabel('DISTANCE (grid units)'),ylabel('nT'),title('SOUTH TO NORTH 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);

axes(haxmain)