function [ spacing, hcond, vcond, names ] = LIN_EM( thickness,conductivity )
%LIN_EM(thickness, conductivity) returns the apparent conductivity that a
%low-induction-number EM instrument would measure for a given conductivity
%structure. 'thickness' contains the thickness of all layers;
%'conductivity' contains the conductivities of all layers (including the
%half-space at the base, so it has one more element than 'thickness'

nlay=length(conductivity);
if (nlay-length(thickness)) ~= 1
    disp('There must be one more conductivity than thickness');
    return
end
z=[0,cumsum(thickness),99999999];
spacing=[1 1.21 3.66 10 20 40];
names=['EM38  ';'EMP400';'EM31  ';'EM34  ';'EM34  ';'EM34  '];
space=['  ';'  ';'  ';'  ';'  ';'  '];

hcond=zeros(size(spacing));
vcond=zeros(size(spacing));
for ilay=1:nlay
    hcond=hcond+(Rh(spacing,z(ilay))-Rh(spacing,z(ilay+1)))*conductivity(ilay);
    vcond=vcond+(Rv(spacing,z(ilay))-Rv(spacing,z(ilay+1)))*conductivity(ilay);
end

disp('Instr  Spcng Horizontal Vertical');
disp([names,space,num2str(spacing'),space,num2str(hcond'),space,num2str(vcond')]);

end

function R=Rv(r,z)
  a=z./r;
  R=1./sqrt(4*a.^2+1);
end

function R=Rh(r,z)
  a=z./r;
  R=sqrt(4*a.^2+1)-2*a;
end
