function fem3loop()
% program FEM3LOOP
%
% GPGN302 EM modeling lab
%
% 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
close all

inpf=uigetfile('*.inp', 'Pick an input file');
fi=fopen(inpf,'rt');

%input target parameters
L=fscanf(fi, '%f',1);
R=fscanf(fi, '%f',1);
sdum=fgetl(fi);


yc=fscanf(fi, '%f',1);
xc=fscanf(fi, '%f',1);
zc=fscanf(fi, '%f',1);
sdum=fgetl(fi);

dincl=fscanf(fi,'%f',1);
ddecl=fscanf(fi,'%f',1);
sdum=fgetl(fi);

%input 
S=fscanf(fi, '%f',1);
ht=fscanf(fi, '%f',1);
f=fscanf(fi, '%f',1);
sdum=fgetl(fi);

%observation grid 
ymin=fscanf(fi,'%f',1);
ymax=fscanf(fi,'%f',1);
dely=fscanf(fi,'%f',1);
sdum=fgetl(fi);

xmin=fscanf(fi,'%f',1);
xmax=fscanf(fi,'%f',1);
delx=fscanf(fi,'%f',1);
sdum=fgetl(fi);

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

% 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);
figure(1)
semilogx(alf,fre,'-',alf,fim,'--',alpha,real(f_factor), '+', alpha,imag(f_factor),'+')
xlabel('\alpha=\omega L /R')
ylabel('Frequency Response')


% 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;

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

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

[nx,ny]=size(x);
kx=ceil(nx/2);

figure(3)
plot(y(kx,:), real_response(kx,:),'-',y(kx,:), imag_response(kx,:),'--')
legend('Real','Quadrature')
xlabel('Easting')
ylabel('B_s/B_p')

yl = get(get(3,'currentAxes'),'ylim');
if -1e-10 < yl(1)
    yl(1) = -1e-10;
end
if 1e-10 > yl(2)
    yl(2) = 1e-10;
end
ylim(yl);

fclose(fi);


% 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
