function G = vsp(nz,zg);
% Makes the tomographic  VSP Matrix  G needed for your problem
% nz: number of layers
% zg: position of geophones 
% msacchi@ualberta.ca 

ng = length(zg);
zmax = zg(end);
dz = zmax/nz;


G = zeros(ng,nz);

for ig = 1:ng
   nl = floor(zg(ig)/dz);
   dif = zg(ig)-nl*dz;
     for iz = 1:nl
      G(ig,iz) = dz;
     end
  if dif>0; G(ig,nl+1)=dif;end
  end