
% Computes induced magnetic field variation across a 
% slab
% inc = 90 (won't change .... need to generalize)

clear all
close all

% Set initial inclination (degrees)
inc =0;
rad = 180/pi;
% inc = 0 is a vertical H field

% Define slab
z1 = 2;   % Upper surface
z2 = 4;   % Lower surface
xl =-5;   % Left edge
xr = 5;    % Right edge
B=50000;  % Background field
k=0.015;

%Define vectors for geometry plot
ground_x = [-20,20]
ground_z = [0,0]
h_bg =[B,B]

% Start loop to vary inclination and depth
irun =1
while irun ~= -1   
   
  % Define vector to draw slab
  slab_x =[xl,xr,xr,xl,xl];
  slab_z =[-z1,-z1,-z2,-z2,-z1];
  h_earth_z =  -B*cos(inc);
  h_earth_x =  -B*sin(inc);
     
  % Earth's field at surface
  hplotx2 = [10+1.5*sin(inc),10-1.5*sin(inc)];
  hplotz2 = [2+1.5*cos(inc),2-1.5*cos(inc)];

  
  % Define points where H is computed on z=0 surface
  ix=0;
  for xdum = -20:0.1:20.0
    ix=ix+1;
    x(ix) =xdum;
  end
  nx=ix;
  
  % Compute Z at each point
  for ix =1:nx     
     theta1_right = atan((xr - x(ix))/z1);
     theta1_left =  atan((xl - x(ix))/z1);
     theta1 = theta1_right-theta1_left;
     theta2_right = atan((xr - x(ix))/z2);
     theta2_left =  atan((xl - x(ix))/z2);
     theta2 = theta2_right-theta2_left;
     Z(ix)= B+2*B*k*(theta1-theta2);  
   end
  
  % Plot magnetic field response
  subplot (2,1,1)
  plot(x,Z)
  hold on
  plot(ground_x,h_bg,':')
  hold off
  xlabel ('x(m)')
  ylabel ('Z (nT)')
  title ('Vertical magnetic field over a slab :  k=0.015')
  axis([-20,20,48000,52000])
  
  %Plot model geometry and inclination
  subplot (2,1,2)
  plot(ground_x,ground_z)
  axis([-20,20, -10,5])
  hold on
    plot(slab_x,slab_z)
    plot(hplotx2,hplotz2)  
    plot(hplotx2(2),hplotz2(2),'.')
    xlabel ('x(m)')
    ylabel ('z(m)')
  hold off
  
  
% Change these parameters if desired 
 choice=menu('Options', 'Change i','Change z1', 'Change z2','Change left','Change right','Quit');
  if choice ==1  
   %inc*rad
   %inc=input('Enter inc:')
   %inc = inc/rad
  end
  
  if choice ==2
     z1
     z1=input('Enter z1:')
  end
  
  if choice ==3
     z2
     z2=input('Enter z2:')
  end
  
  if choice ==4
     xl
     xl=input('Enter xl:')
  end

  if choice ==5
     xr
     xr=input('Enter xr:')
  end

  if choice > 5
   irun=-1;
   print -dpsc  mag_slab.ps
  end
end