function [conductivity,synthetic]=LIN_EM_inv(v_spacing,v_cond,h_spacing,h_cond,thickness);
%LIN_EM(v_spacing,v_cond,h_spacing,h_cond,thickness) performs a flattened 
%least-squares linear inversion for a layered model, given a set of
%apparent conductivities for vertical and horizontal modes. For example,
%to invert for a two-layer model with an upper-layer thickness of 1 m using
%EM31 data of 120 mS/m in horizontal mode and 80 mS/m in vertical mode, the
%command would be
%    [conductivity,synthetic]=LIN_EM_inv([3.66],[80],[3.66],[120],[1]);
% (the square brackets aren't necessary, but indicate that these are
% row vectors). The result will be:
%    conductivity =
%      203.1585
%       62.8681
%    synthetic =
%       80.0497
%      119.9501
% If we also have a horizontal EM38 measurement of 150 mS/m andLIN want to
% invert for a series of 50 cm thick layers down to 3 m, we can do this:
%    [conductivity,synthetic]=LIN_EM_inv([3.66],[80],[3.66 1],[120 150],[0.5:0.5:3]);
% and we get:
%    conductivity =
%      172.8650
%      149.6994
%      103.1865
%       69.1387
%       48.5877
%       33.0424
%       16.3804
%    synthetic =
%       83.2058
%      116.6457
%      150.9645
% (note the imperfect fit to the data -- we can also adjust the damping and
% flattening weights in the function. For this run I raised the flattening
% weight to 0.1)

% Base model assumption: use average of all apparent conductivities
basemod=mean([v_cond,h_cond]);

z=[0,cumsum(thickness),99999999];
Nlay=(length(z)-1);
Nv=length(v_spacing);
Nh=length(h_spacing);

% Weights for flattening and damping
flat_weight=0.01;
damp_weight=0.0001;

% Build sensitivity matrix
for iv=1:Nv
    for ilay=1:Nlay
        sensitivity(iv,ilay)=Rv(v_spacing(iv),z(ilay))-Rv(v_spacing(iv),z(ilay+1));
    end
end
for ih=1:Nh
    for ilay=1:Nlay
        sensitivity(Nv+ih,ilay)=Rh(h_spacing(ih),z(ilay))-Rh(h_spacing(ih),z(ilay+1));
    end
end

% Build first-difference matrix (flattener)
flattener=zeros(Nlay-1,Nlay);
for irow=1:(Nlay-1)
    flattener(irow,irow)=-1;
    flattener(irow,irow+1)=1;
end

% Build identity matrix (damper)
damper=eye(Nlay);

% Assemble full matrix
B=[sensitivity; flat_weight*flattener; damp_weight*damper];

% Right-hand side for equation
rhs=[v_cond';h_cond';zeros(Nlay-1,1);ones(Nlay,1)*basemod];

% Invert
conductivity=B\rhs;

% Synthetic data
synthetic=sensitivity*conductivity;

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
