% PHDmagloop2   Plot the fields from a magnetic dipole in rock

% based on PHDmagloop1, this plots the fields in rock, using the
% approximate formulations that assume no boundary problem
% start at fig 4. Use ROCK1 fundtion instead of AIR

% later, plot field using the exact expression

%clear; clf; cla;

filetype = '-dtiff';

points=1001;

while 0,

% --------------------------- plot magnitude -----------------

fig=4; figure(fig);  

x=logspace(-1,2,points);
[Hax Hcp] = rock1(x);

h1 = loglog(x,abs(Hax),'k-','linewidth',3);
hold on;
h2 = loglog(x,abs(Hcp),'b-','linewidth',1.5);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Magnetic Field (Magnitude [ref axial @ T=1])','fontsize',12)
title('Variation of field with distance (rock)','fontsize',12)
legend('coaxial loops','coplanar loops')
axis([0.1 20 1E-10 1E4])
print(filetype, '', ['loop-' num2str(fig)] )

% --------------------------- plot phase -----------------

fig=5; figure(fig);  

x1=linspace(0.1,15,points/2);
[Hax Hcp] = rock1(x1);

h1 = plot(x1,(180/pi)*angle(Hax),'k.','markersize',12);
hold on;
h2 = plot(x1,(180/pi)*angle(Hcp),'b.','markersize',6);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Magnetic Field (Phase, degrees)','fontsize',12)
title('Variation of field with distance (rock)','fontsize',12)
legend('coaxial loops','coplanar loops')
print(filetype, '', ['loop-' num2str(fig)] )

% --------------------------- plot phase (zoom) -----------------

fig=6; figure(fig);  

x1=linspace(0.1,2,points);
[Hax Hcp] = rock1(x1);

h1 = plot(x1,(180/pi)*angle(Hax),'k-','linewidth',3);
hold on;
h2 = plot(x1,(180/pi)*angle(Hcp),'b-','linewidth',1.5);
%ax = axis;
%axis([0 2 ax(3:4)])

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Magnetic Field (Phase, degrees)','fontsize',12)
title('Variation of field with distance (rock)','fontsize',12)
legend('coaxial loops','coplanar loops')
print(filetype, '', ['loop-' num2str(fig)] )


% --------------------------- plot ratio of magnitudes

fig=7; figure(fig);  

x=logspace(-1,2,points);
[Hax Hcp] = rock1(x);

h1 = loglog(x,abs(Hax)./abs(Hcp),'k-','linewidth',3);
hold on
loglog([0.1 100],[1 1], 'b-','linewidth',1.5);
loglog([1 1],[0.001 10], 'b-','linewidth',1.5);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Magnetic Field (Magnitude)','fontsize',12)
title('Variation of field with distance: ratio (rock)','fontsize',12)
legend('axial/coplanar ratio')
axis([0.1 100 0.01 3])
print(filetype, '', ['loop-' num2str(fig)] )

% --------------------------- plot exact expression -----------------

fig=8; figure(fig);  

% use the exact expression, in subroutine PQintegral
% which takes arguments D,T,Z and MODE.
% For vertical (Q) field, mode=3
% Note that PQintegral uses Shope/Wait definition of T
% so need to multiply 'my' T by sqrt(2).

Z=0; % Z is normalised height of observation point
D=0; % D is normalised offset of observation point
T=0; % T is normalised depth of loop

points=51;
x=linspace(0.1,10,points);

[Hax Hcp] = rock2(x);
% this represents the coefficient in a similar way to PQ integral, 
% i.e. field is [Hax Hcp]/2.pi.x^3
% note factor of 2 required...
Hax = Hax / 2;

for i=1:points,
   T=x(i);
   
   % calling PQintegral with mode=3 returns Q field value
   y(i)=(PQintegral(D,T*sqrt(2),0,3));
   
end

h1 = semilogx(x,abs(Hax)./abs(y),'k-','linewidth',3);
hold on
%h1 = loglog(x,abs(Hax),'k-','linewidth',3);
%h1 = loglog(x,abs(Hexact),'k-','linewidth',3);
semilogx([0.1 100],[1 1], 'b-','linewidth',1.5);
semilogx([1 1],[1E-6 1E4], 'b-','linewidth',1.5);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Relative Magnetic Field (Magnitude)','fontsize',12)
title('whole-space/Half-space model: field v. distance','fontsize',12)
%legend('ratio', 'coaxial loops (half-space model)', ...
%    'coaxial loops (whole-space model)')
%axis([0.1 10 1E-6 1E4])
axis([0.1 10 0.7 1.1])

print(filetype, '', ['loop-' num2str(fig)] )

% --------------------------- plot magnitude without T^-3 ---------

fig=9; figure(fig);  

% repeat figure 4, but use ROCK2 which doesnt include the 1/R^3 term

x=logspace(-1,2,points);
[Hax Hcp] = rock2(x);

h1 = loglog(x,abs(Hax),'k-','linewidth',3);
hold on;
h2 = loglog(x,abs(Hcp),'b-','linewidth',1.5);

loglog([0.1 100], [0.71 0.71], 'b-','linewidth',1.5);
loglog([3.85 3.85], [0.001 10], 'b-','linewidth',1.5);

%peak = 1.6183; [Hax Hcp] = rock2(peak); abs(Hcp)
%plot(peak,abs(Hcp),'r.');

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Frequency','fontsize',12)
ylabel('Magnetic Field Coefficient (Magnitude)','fontsize',12)
title('Variation of field with distance (rock)','fontsize',12)
legend('coaxial loops','coplanar loops',3)
ax = axis;
axis([0.1 10 0.1 3])
print(filetype, '', ['loop-' num2str(fig)] )

    
% --------------------------- plot magnitude witht T^2 ---------

fig=10; figure(fig);  

% repeat figure 4, but use ROCK3 which doesnt include the 1/R^3 term

x=logspace(-1,2,points);
[Hax Hcp] = rock3(x);

h1 = loglog(x,abs(Hax),'k-','linewidth',3);
hold on;
h2 = loglog(x,abs(Hcp),'b-','linewidth',1.5);

%peak = sqrt(8); [Hax Hcp] = rock3(peak); abs(Hcp)
%plot(peak,abs(Hax),'r.');

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Frequency','fontsize',12)
ylabel('Magnetic Field Coefficient (Magnitude)','fontsize',12)
title('Variation of field with distance (rock)','fontsize',12)
legend('coaxial loops','coplanar loops',3)
ax = axis;
axis([0.1 10 1E-1 20])
print(filetype, '', ['loop-' num2str(fig)] )

end; % matches while

% --------------------------- plot comparison of approxs to model

fig=11; figure(fig);  

x=logspace(-1,2,points);
[Hax Hcp] = rock1(x);

% compare the complete full-space conducting medium model
% with some approximations. Firstly, a basic near-field model
%scaled so that the value at x =approx= 0 is the same as
% that with the complete expression, which was normalised.
% (see ROCK1 function)

Htry0 = 1./x.^3;
Htry0 = Htry0 ./2 ; % for coplanar
Htry0 = Htry0 .* 2 .* exp(1) / (2 * sqrt(5));

Htry1 = exp(-x).* Htry0;
Htry2 = Htry1 .* sqrt(x.^2 +(1+x).^2);

% now plot the ratio of those approximations to the coaxial field 

h2 = loglog(x,Htry0./abs(Hcp),'k-','linewidth',1.5);
hold on
h2 = loglog(x,Htry2./abs(Hcp),'b-','linewidth',4);
h2 = loglog(x,Htry1./abs(Hcp),'k-','linewidth',2.5);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Error in Magnetic Field (Magnitude)','fontsize',12)
title('Error in coplanar approximations v. distance(rock)','fontsize',12)
legend('(a)' , '(b)' , '(c)',2)
axis([0.1 10 0.1 10])
print(filetype, '', ['loop-' num2str(fig)] )

% --------------------------- plot comparison of approxs to model
% as previous, but coaxial fields

fig=12; figure(fig);  

x=logspace(-1,2,points);
[Hax Hcp] = rock1(x);

% compare the complete full-space conducting medium model
% with some approximations. Firstly, a basic near-field model
%scaled so that the value at x =approx= 0 is the same as
% that with the complete expression, which was normalised.
% (see ROCK1 function)

Htry0 = 1./x.^3;
%Htry0 = Htry0 ./2 ; % for coplanar
Htry0 = Htry0 .* 2 .* exp(1) / (2 * sqrt(5));

Htry1 = exp(-x).* Htry0;
Htry2 = Htry1 .* sqrt(x.^2 +(1+x).^2);

% now plot the ratio of those approximations to the coaxial field 

h2 = loglog(x,Htry0./abs(Hax),'k-','linewidth',1.5);
hold on
h2 = loglog(x,Htry2./abs(Hax),'b-','linewidth',4);
h2 = loglog(x,Htry1./abs(Hax),'k-','linewidth',2.5);

set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('T: Normalised Distance','fontsize',12)
ylabel('Error in Magnetic Field (Magnitude)','fontsize',12)
title('Error in coaxial approximations v. distance(rock)','fontsize',12)
legend('(a)' , '(b)' , '(c)',2)
axis([0.1 10 0.1 10])
print(filetype, '', ['loop-' num2str(fig)] )

% ----------------------------- bandwidth calculations....
break

format compact

figure(9)

x=logspace(-1,2,points);
[Hax Hcp] = rock2(x);
Y = abs(Hax);  % or 'cp' as required
[v i] = max(Y);

ipeak = i; 
x_vpeak = x(ipeak)
vpeak = Y(i)
plot(x(ipeak),vpeak,'c.');
% that plots the peak of the curve, and prints its x&y value

% now look for the -3db point...
v = 2/sqrt(2); % or 1/sqrt(2) for 'cp'
clear f; f = find(Y>v);

imin = f(1); imax = f(length(f)); vmin = Y(imin); vmax = Y(imax);

% now we have min and max co-ords to plot...

x_vmin = x(imin)
y_vmin = vmin
x_vmax = x(imax)
y_vmax = vmax

plot(x(imin),vmin,'r.');
plot(x(imax),vmax,'y.');

% and the bandwidth is....
bandwidth = x(imax)-x(imin)

% ----------------------------- bandwidth calculations....

figure(10)

x=logspace(-1,2,points);
[Hax Hcp] = rock3(x);
Y = abs(Hax);  % or 'cp' as required
[v i] = max(Y);

ipeak = i; 
x_vpeak = x(ipeak)
vpeak = Y(i)
plot(x(ipeak),vpeak,'c.');
% that plots the peak of the curve, and prints its x&y value

% now look for the -3db point...
v = vpeak/sqrt(2);
clear f; f = find(Y>v);

imin = f(1); imax = f(length(f)); vmin = Y(imin); vmax = Y(imax);

% now we have min and max co-ords to plot...

x_vmin = x(imin)
y_vmin = vmin
x_vmax = x(imax)
y_vmax = vmax

plot(x(imin),vmin,'r.');
plot(x(imax),vmax,'y.');

% and the bandwidth is....
bandwidth = x(imax)-x(imin)


