% PHD_QvF3  plots diagrams of attenuation v frequency for different depths
% based on QvF1 but with additions for atmospheric noise
clear
filetype = '-dtiff';

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

% --- edit program for...

resolution = 401;
% 21 points is OK for draft, but 401 points is needed to get accurate 
% data on optimum frequency. With 401 points the program line...
% [optimum(:,2)./optimum(:,4) optimum(:,3)-optimum(:,5)]
% produces an output of, typically,
%    1.09647819614319   1.41897172987125
% showing that the exact method produces an optimum frequency 9.6% 
% higher, with an attenuation of 1.4db less than the approximate method.

%sigma = 0.01; fig = 1;
sigma = 0.001; fig = 2;

fig=3; % for noise version

figure(3); clf;
figure(4); clf;

points=resolution;
Amax = points; A2max = points;
X = logspace(-1, 3, points); 


   % --------------------- add in the effect of atmospheric noise
   
   % we just need to divide A by Fk.
   % We note that Fk is proportional to freq^2 . sqrt(Fa)
   % and that Fa takes the extremes, as written in the thesis, of...
   % Max: 172dB at 10kHz, falling at 36dB/decade
   % Min: Max of: 148dB at 10kHz, falling at 74dB/decade 
   %         and: 123dB at 10kHz, falling at 28dB/decade 
   % We wont plot the 'no noise' curve as it will be very similar to
   % the max curve (40dB/decade corresponding to 'no noise effect with freq'
   % and we'll scale all the values so that the range of SNR is similar
   % to the Y-axis of the original no noise plot, so...
   
   rel_freq = X/10;
   
   FKmax = -172 +       172 + (40 - 36)*log10(rel_freq);
   FKmin = -172 + max ( 148 + (40 - 72)*log10(rel_freq) ,...
                        123 + (40 - 28)*log10(rel_freq) );
   
   % now add the modification for <10kHz
   
   for i=1:points,
       FKmax(i) = -172 + 172 + (40 - 23)*log10(rel_freq(i));
       FKmin(i) = -172 + 148 + (40 - 31)*log10(rel_freq(i));
       if rel_freq(i) > 1
           break
       end
   end    


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

for depth= [linspace(50,500,10) linspace(600,1000,5)], % 50:50:1000,
%for depth= 50:50:100,
 
   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
      
      % skip the exact calculation as it takes too long
      
      while 0,
      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)
      end; % matches while
      
      % for the approximation, use different definition of T
      
      T = T/sqrt(2);
      r = depth * sqrt((1+Z)^2 + D^2); 
      A(i) = exp(-T) * sqrt(T^2 + (1+T)^2);
      A(i) = A(i) * freq * mu/(r^3);
      A(i) = A(i) / (mu * 2 * pi);
      A(i) = 20 * log10(A(i));
      
      % to get a better indication of SNR, add the missing terms,
      % 2pi.mu.moment/sqrt(4kTB), using Momeent of 30 and B = 1000
      
      A(i) = A(i) + 95;
      
   end
   
   
   % -------------- now plot the data ----------------
 fig=3; figure(fig);  

   H1 = semilogx(X,A,'k-','linewidth',1.2);  
   hold on;

   H5 = semilogx(X,A - FKmax,'b-','linewidth',3);  
   [value index] = max(A - FKmax);
   if index ~= points,
      H2 = plot(X(index),A(index) - FKmax(index),'ko','linewidth',2,'markersize',6);
   end

   [value2 index2] = max(A - FKmin);

   disp([depth X(index) A(index) X(index2) A(index2)])
   optimum(k,:) = [depth X(index) A(index)- FKmax(index) X(index2) A(index2)- FKmin(index2)]; k=k+1;   
   
   % -------------- now re-plot the data on a new fig ----------------

  fig=4; figure(fig);  
   
%   H1 = semilogx(X,A,'k-','linewidth',1.2);  
%   hold on;

   H7 = semilogx(X,A - FKmax,'b-','linewidth',1.2);  
   hold on
   [value index] = max(A - FKmax);
   if index ~= points,
      H2 = plot(X(index),A(index) - FKmax(index),'ko','linewidth',1.2,'markersize',4);
   end

   H8 = semilogx(X,A - FKmin,'b-','linewidth',3);  
   [value2 index2] = max(A - FKmin);
   if index ~= points,
      H2 = plot(X(index2),A(index2) - FKmin(index2),'kx','linewidth',2.5,'markersize',6);
   end

end

   % -------------- tidy up the plots ----------------

fig=3; figure(fig);  
   
plot(optimum(:,2),optimum(:,3),'k-','linewidth',1.5); % MAX NOISE
%plot(optimum(:,4),optimum(:,5),'k-','linewidth',1.5)

axis([0.1 1000 -60 120]);
set(gca, 'YTick', linspace(-60, 120,10));
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
parms = ['D=', num2str(D), ' depth=50:50:500 & 600:100:1000, Z=',num2str(Z), ', mode=',num2str(mode) ];
xlabel('Frequency (kHz)','fontsize',12)
ylabel('Signal/Noise ratio (dB)','fontsize',12)
title(['Opt Freq: ', parms],'fontsize',12)
legend([H1 H5], ...
['no noise, ', num2str(sigma), ' S/m'],...
'maximum noise',2)

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


fig=4; figure(fig);  
   
plot(optimum(:,2),optimum(:,3),'k-','linewidth',1.5)
plot(optimum(:,4),optimum(:,5),'k-','linewidth',1.5)

axis([0.1 1000 -60 120]);
set(gca, 'YTick', linspace(-60, 120,10));
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
parms = ['D=', num2str(D), ' depth=50:50:500 & 600:100:1000, Z=',num2str(Z), ', mode=',num2str(mode) ];
xlabel('Frequency (kHz)','fontsize',12)
ylabel('Signal/Noise ratio (dB)','fontsize',12)
title(['Opt Freq: ', parms],'fontsize',12)
legend([H7 H8], ...
'maximum noise', 'minimum noise',2)

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


% --------------- plot FK and FA --------------------------

fig=5; figure(fig); 

H1 = semilogx(X,FKmax+38,'b-','linewidth',3);  
% i.e. scale to match table in thesis
hold on
H2 = semilogx(X,FKmin+38,'b-','linewidth',1.5);  
 
axis([0.1 1000 -20 60]);
set(gca, 'YTick', linspace(-20,60,9));
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('Frequency (kHz)','fontsize',12)
ylabel('Fk Noise Term (dB)','fontsize',12)
title('Atmospheric & Man-made noise','fontsize',12)
legend([H1 H2 ], ...
    'maximum Fk', 'minimum Fk',3)

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


% ---------------------now convert to Fa and re-plot

fig=6; figure(fig);  
   
H1 = semilogx(X,FKmax+38 + 294 - 160 - 40*log10(rel_freq),'k-','linewidth',3);  
hold on
H2 = semilogx(X,FKmin+38 + 294 - 160 - 40*log10(rel_freq),'k-','linewidth',1.5);  

axis([0.1 1000 60 220]);
set(gca, 'YTick', linspace(60,220,9));
set(gca,'LineWidth',1.5,'FontSize',12,'TickDir','out')
xlabel('Frequency (kHz)','fontsize',12)
ylabel('Fa Noise Term (dB)','fontsize',12)
title('Atmospheric & Man-made noise','fontsize',12)
legend([H1 H2], ...
    'maximum Fa', 'minimum Fa',3)


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