function fem3loop()
% program FEM3LOOP
%
% GPGN302 EM modeling lab
% EOSC350 EM lab (bugfixed/modified): GUI for inpout parameters
%
% Computes the quadrature and in-phase responses in a frequency-domain EM
% survey over a confined conductor using the three-loop model. The
% computation is meant for illustration only and it cannot be used 
% to represent a real-world problem
%
% Assumptions:
% (1) a horizontal coplanar loop-loop system with unit area for both
%     transmitter and receiver.
% (2) loop-loop system is aligned in east-west direction
% (3) a tilted loop representing a target with a user-input inductance,
%     resistance, and locarion of bruial

%
% Input parameters:
% L, R           !! inductance and resistance of the buried loop
% ec, nc, dc     !! easting, northing, and depth of the target loop
% inc, decl      !! inclination and declinartion of the normal vector of target
% S, h, f        !! separattion between Tx and Rx loops, height above the ground, and frequency
% emin, emax, de !! min and max position of obs grid in easting
% nmin, nmax, dn !! min and max position of obs grid in northing
%
%
clear all

prompt = {'Inductance (L) and resistance (R) of the target',...
          'Easting, northing, and depth of the target',...
          'Inclination and declination of the normal vector of the target',...
          'Tx-Rx loop separation, height above ground, and frequency',...
          'Min, max, and interval of observation grid in easting',...
          'Min, max, and interval of observation grid in northing'};
dlg_title = 'EM responses of 3-loop model: fem3loop';
num_lines = 1;
def = {'0.1    2000','0  0  2','0  90','4  0  10000','-10  10  0.25','-10  10  0.25'};
answer = inputdlg(prompt,dlg_title,num_lines,def);

if isempty(answer)
    return;
else
    temp = str2num(answer{1});
    L = temp(1); 
    R = temp(2);
    temp = str2num(answer{2});
    yc = temp(1);
    xc = temp(2);
    zc = temp(3);
    temp = str2num(answer{3});
    dincl = temp(1);
    ddecl = temp(2);
    temp = str2num(answer{4});
    S = temp(1);
    ht = temp(2);
    f = temp(3);
    temp = str2num(answer{5});
    ymin = temp(1);
    ymax = temp(2);
    dely = temp(3);
    temp = str2num(answer{6});
    xmin = temp(1);
    xmax = temp(2);
    delx = temp(3);
end

% generate the grid
xp=xmin:delx:xmax;
yp=ymin:dely:ymax;
[y,x]=meshgrid(yp,xp);
z=0.0*x-ht;

% set up the response arrays
real_response=0.0*x;
imag_response=0.0*x;

% frequency characteristics
alpha=2*pi*f*L/R;


f_factor=(alpha^2+i*alpha)/(1+alpha^2);

amin=0.01;
amax=100;
da=4/40;
alf=-2:da:2;
alf=10.^alf;

fre=alf.^2./(1+alf.^2);
fim=alf./(1+alf.^2);


% simulate anomalies
yt=y-S/2;
yr=y+S/2;

dm=-S/2;
dp= S/2;

M13=mind(0,dm,0,90,0, 0, dp, 0, 90,0);
M12=L*mind(x,yt,z,90.0,0.0,xc,yc,zc,dincl,ddecl);
M23=L*mind(xc,yc,zc,dincl,ddecl,x,yr,z,90.0,0.0);

c_response=-M12.*M23*f_factor/(M13*L);

%scaled to simulate a net volumetric effect
real_response=real(c_response)*1000;
imag_response=imag(c_response)*1000;

% plot output
figstr = [ 'fem3loop modeling parameters: L = ', num2str(L), ', R = ', num2str(R),...
           ', Target location = ', num2str([yc xc zc]),...
           ', Target inclination = ', num2str(dincl),...
           ', Target declination = ', num2str(ddecl),...
           ', Loop separation = ', num2str(S),...
           ', Height = ', num2str(ht),...
           ', Frequency = ', num2str(f)];
figure('numbertitle','off','name',figstr);
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure. 

subplot(2,2,1);
semilogx(alf,fre,'.-',alf,fim,'.--');%,alpha,real(f_factor), '+', alpha,imag(f_factor),'+')
hold on;
plot([alpha alpha],[0 1],'r-');
legend('Real','Quadrature');
xlabel('\alpha=\omega L /R');
ylabel('Frequency Response');
title('Plot 1: EM responses of loop model');

subplot(2,2,3);
pcolor(yp,xp,real_response);
shading interp;
axis equal; axis tight;
xlabel('Easting');
ylabel('Northing');
title('Plot 3: Real Component');
colorbar;

subplot(2,2,4);
pcolor(yp,xp,imag_response);
shading interp;
axis equal; axis tight;
xlabel('Easting');
ylabel('Northing');
title('Plot 4: Quadrature Component');
colorbar;

subplot(2,2,2);
[nx,ny]=size(x);
kx=ceil(nx/2);
plot(y(kx,:), real_response(kx,:),'.-',y(kx,:), imag_response(kx,:),'.--')
legend('Real','Quadrature');
xlabel('Easting');
ylabel('H_s/H_p');
title(['Plot 2: EW cross section along Northing = ',num2str(x(kx,1))]);

end


% end function fem3loop


function outa=mind(x,y,z,dincl,ddecl,x0,y0,z0, aincl,adecl)

cons=1.0;

di=pi*dincl/180.0;
dd=pi*ddecl/180.0;

cx=cos(di)*cos(dd);
cy=cos(di)*sin(dd);
cz=sin(di);


ai=pi*aincl/180.0;
ad=pi*adecl/180.0;

ax=cos(ai)*cos(ad);
ay=cos(ai)*sin(ad);
az=sin(ai);


% begin the calculation
a=x-x0;
b=y-y0;
h=z-z0;

rt=sqrt(a.^2+b.^2+h.^2).^5;

txy=3*a.*b./rt;
txz=3*a.*h./rt;
tyz=3*b.*h./rt;

txx=(2*a.^2-b.^2-h.^2)./rt;
tyy=(2*b.^2-a.^2-h.^2)./rt;
tzz=-(txx+tyy);

bx=cons*(txx*cx+txy*cy+txz*cz);
by=cons*(txy*cx+tyy*cy+tyz*cz);
bz=cons*(txz*cx+tyz*cy+tzz*cz);

outa=bx*ax+by*ay+bz*az;

% end of function mdipole
end