function VLF_1

clear all;
close all;

%Top of conductor
y1 = 0.; z1 = 45;
%Bottom of conductor
y2 = 0.; z2 = 95;  w = 20.;
% Primary magnetic field
HYP = 1; quad_response=0.9

%Current flowing in conductor
i1 =300;

scale = 10.;  scale_quad = 50.; 
rlim = 100;

% Plotting paramters
y_surface = [-400,400]; z_surface = [-5.,-5.];
y_conductor = [y1-w,y1+w,y2+w,y2-w,y1-w]; 
z_conductor = [z1,z1,z2,z2,z1];

% Array where magnetic fields will be computed
y = [-300:20:300]; ny=length(y);
z = [-100:20:30]; nz=length(z);
limits =[-300,300,-100,100]
% z is positive in the ground

% Generate uniform primary field
for iy=1:ny;
   for iz=1:nz;
      hyp(iy,iz)=1; hzp(iy,iz)=0.;
   end
end

% Generate secondary field (assume in phase)
for iy=1:ny;
   for iz=1:nz;
      
      % Current at top of plate
      r1=sqrt((y(iy)-y1)^2+(z(iz)-z1)^2); hr1=i1/(2*pi*r1);   
      
      % Current at top of plate
      r2=sqrt((y(iy)-y2)^2+(z(iz)-z2)^2); hr2=i1/(2*pi*r2);          
      hys(iy,iz)= hr1*(z1-z(iz))/r1 -hr2*(z2-z(iz))/r2 ;      
      hzs(iy,iz)=-hr1*(y1-y(iy))/r1 +hr2*(y2-y(iy))/r2 ; 
      hms(iy,iz) = sqrt(hys(iy,iz)^2+hzs(iy,iz)^2);
      
      hys_quad(iy,iz)= quad_response*(hr1*(z1-z(iz))/r1 -hr2*(z2-z(iz))/r2) ;      
      hzs_quad(iy,iz)= quad_response*(-hr1*(y1-y(iy))/r1 +hr2*(y2-y(iy))/r2) ; 
      hm_quad(iy,iz) = sqrt(hys_quad(iy,iz)^2+hzs_quad(iy,iz)^2);
      
      hyt(iy,iz) = hys(iy,iz)+hyp(iy,iz);  hzt(iy,iz) = hzs(iy,iz)+hzp(iy,iz);
      hmt(iy,iz) = sqrt(hyt(iy,iz)^2+hzt(iy,iz)^2);
   end
end

%=============================================================
%Plot primary field
%=============================================================
figure(1);
subplot(2,1,1)
hold on
fill(y_conductor,-z_conductor,'r');
for iy=1:ny;
   for iz=1:nz;      
      y_plot = [y(iy),y(iy)+scale*hyp(iy,iz)];z_plot = -[z(iz),z(iz)+scale*hzp(iy,iz)];
      plot(y_plot,z_plot,'k'); plot(y_plot(1),z_plot(1),'k.')
   end
end
plot(y_surface,z_surface,'k'); hold off
xlabel('y(m)'); ylabel('z(m)'); title('VLF Primary magnetic field : In phase');
axis(limits)
pause
print (1,'-dpsc', 'VLF_1_inphase_primary.ps')

%=============================================================
% Plot secondary field
%=============================================================
figure(2);
subplot(2,1,1)
hold on
fill(y_conductor,-z_conductor,'r');
for iy=1:ny;
   for iz=1:nz;
      y_plot = [y(iy),y(iy)+scale*hys(iy,iz)/hms(iy,iz)]; z_plot = -[z(iz),z(iz)+scale*hzs(iy,iz)/hms(iy,iz)];
      plot(y_plot,z_plot,'k');      plot(y_plot(1),z_plot(1),'k.')
   end
end
plot(y_surface,z_surface,'k'); hold off
xlabel('y(m)')
ylabel('z(m)'); title('VLF Secondary magnetic field : In phase')
axis(limits)
pause
print (2,'-dpsc', 'VLF_1_inphase_secondary.ps')

%=============================================================
% Plot in-phase total field
%=============================================================
figure(3);
subplot(2,1,1)
hold on
fill(y_conductor,-z_conductor,'r');
for iy=1:ny;
   for iz=1:nz;
      %y_plot = [y(iy),y(iy)+scale*hyt(iy,iz)/hmt(iy,iz)];
      %z_plot = -[z(iz),z(iz)+scale*hzt(iy,iz)/hmt(iy,iz)];
      y_plot = [y(iy),y(iy)+scale*hyt(iy,iz)];
      z_plot = -[z(iz),z(iz)+scale*hzt(iy,iz)];
      plot(y_plot,z_plot,'k'); plot(y_plot(1),z_plot(1),'k.')
   end
end


plot(y_surface,z_surface,'k'); 
hold off

xlabel('y(m)'); ylabel('z(m)')
title('VLF Total magnetic field : In phase')
axis(limits)



print (3,'-dpsc', 'VLF_1_inphase_total.ps')

pause

%===============================================================
% Plot quad total field
%===============================================================
figure(4);
subplot(2,1,1)
hold on
fill(y_conductor,-z_conductor,'r');

for iy=1:ny;
   for iz=1:nz;
      %y_plot = [y(iy),y(iy)+scale*hys_quad(iy,iz)/hm_quad(iy,iz)];
      %z_plot = -[z(iz),z(iz)+scale*hzs_quad(iy,iz)/hm_quad(iy,iz)];
      y_plot = [y(iy),y(iy)+scale_quad*hys_quad(iy,iz)];
      z_plot = -[z(iz),z(iz)+scale_quad*hzs_quad(iy,iz)];
      plot(y_plot,z_plot,'r') ;  plot(y_plot(1),z_plot(1),'k.')
   end
end

plot(y_surface,z_surface,'k'); 
hold off

xlabel('y(m)'); ylabel('z(m)')
title('VLF Total magnetic field : Quadrature')
axis(limits)


print (4,'-dpsc', 'VLF_1_quad_total.ps')