% PHD_QvF2  plots diagrams of attenuation v frequency for different sigma

clear;
filetype = '-dtiff';

D=0; T=0; Z=0; mode=3;
mu = 4E-7 * pi;
k=1; % needs initialising here

% --- figure numbering carries on from QvF1

sigmarange = 10.^(-4 : 1/3 : -2);
fig = 3;
depth = 300;

figure(fig);  

% ------------------------- plot optimum frequency

for sigma = sigmarange,
   points=101;
   Amax = points; A2max = points;
   X = logspace(-1, 3, points); 
   for i=1:points,
      freq = X(i) * 1000;
      delta = sqrt(1/(pi*freq*mu*sigma));
      T = sqrt(2)* depth / delta;  % This is Wait's version of T
      Q(i) = abs(PQintegral(D,T,Z,mode));
      A(i) = Q(i) * freq * mu/(depth^3);
      A(i) = A(i) / (mu * 2 * pi);
      A(i) = 20 * log10(A(i));
      
      % A(i) is now attenuation in dB (arbitrary reference, used in early work)
                 
      % for the approximation, use different definition of T
      % and add the extra factor of 0.95
      
      T = T/sqrt(2) * 0.95;
      r = depth * sqrt((1+Z)^2 + D^2); 
      A2(i) = exp(-T) * sqrt(T^2 + (1+T)^2);
      A2(i) = A2(i) * freq * mu/(r^3);
      A2(i) = A2(i) / (mu * 2 * pi);
      A2(i) = 20 * log10(A2(i));
      
   end
   
   % -------------- now plot the data ----------------
   
   H1 = semilogx(X,A,'k-','linewidth',3);  
   hold on;
   [value index] = max(A);
   if index ~= points,
      H2 = plot(X(index),A(index),'kx','linewidth',2,'markersize',6);
   end

   H3 = semilogx(X,A2,'b-','linewidth',1.5);
   [value2 index2] = max(A2);
%  if index ~= points,
%     H4 = plot(X(index2),A2(index2),'bo','linewidth',2,'markersize',6);
%  end
    
   disp([depth X(index) A(index) X(index2) A2(index2)])
   optimum(k,:) = [depth X(index) A(index) X(index2) A2(index2)]; k=k+1;   
   
end

% [optimum(:,2)./optimum(:,4) optimum(:,3)-optimum(:,5)]
%
% plot(optimum(:,2),optimum(:,3),'k-','linewidth',2)
% plot(optimum(:,4),optimum(:,5),'b-','linewidth',1)

axis([0.1 1000 -180 0]);
set(gca, 'YTick', linspace(-180,0,11));
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
%parms = ['D=', num2str(D), ' depth=50:50:1000, Z=',num2str(Z), ', mode=',num2str(mode) ];
parms = ['D=', num2str(D), ' logsigma=-4:1/3:-2, Z=',num2str(Z), ', mode=',num2str(mode) ];
xlabel('Frequency (kHz)','fontsize',12)
ylabel('Induced voltage (dB)','fontsize',12)
title(['Optimum Frequency: ', parms],'fontsize',12)
legend([H1 H2 H3], ...
['exact solution for depth ', num2str(depth), ' m'],...
'optimum frequency', 'revised approximation',2)

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