getf("maketimes.sci");
getf("makejacobian.sci");



svalsused = input("number of singular values used")
if svalsused == [] then
  svalsused = 5;
  printf("using default svalsused")
end



// parameters of synthetic noise
exactvar = 1;
corlen = 0;
dist = 'U';

ncols = 200;   // size of model used in reconstruction.  can be changed
               // arbitrarily.  The "true" model is picewise continuous
               // with a low-velocity zone.  The true data were simulated 
	       // using 1000 layers, as if piecewise continuous.
	          
nrows = 100;   // number of data
ztop = 0;      // physical dimensions of the model
zbottom = 40;

deltaz = (zbottom - ztop)/ncols;
firstrec = ztop + deltaz/2;
lastrec = zbottom ;
deltarec = (lastrec - firstrec)/nrows;
rec = firstrec:deltarec:lastrec-deltarec;


// create the velocity model.  The exact model has a linear increase
// in velocity and a low velocity zone in the middle (irrespective
// of parameterization).

vmax = 4;
vmin = .5;
dv = (vmax - vmin)/ncols;
exactvelmodel=vmin:dv:vmax-dv;
mid = floor(ncols/2);
deltam = floor(.1 * ncols);
exactvelmodel(mid-deltam:mid+deltam) = 1;
exactslowmodel = 1 ./ exactvelmodel;


// this is just for zero-offset, straight ray tomography
jacobian =  makejacobian(nrows,ncols,deltaz,rec);  


// starting model
backgroundvelmodel(1:ncols) = 2;
backgroundslowmodel = 1 ./ backgroundvelmodel;

backgroundtt = jacobian * backgroundslowmodel;

// creat syntheic travel times.  junk is the sample variance of the 
// synthetic noise
[times,junk] = maketimes(jacobian,exactslowmodel,1000,100,dist,0,exactvar,corlen);     
deltatt = times - backgroundtt;

//////////////////////////end of setup//////////////////////////////////////


// do an svd of the jacobian
[U,s,V] = svd(jacobian);

// construct truncates svd approximation
lambda = diag(s);
tsvd = zeros(1,ncols);
for i = 1:svalsused
  factor = U(:,i)' * deltatt / lambda(i);
  tsvd = tsvd + (V(:,i) * factor)';
end

// update the background model, vel or slowness
tsvdslowupdate = backgroundslowmodel + tsvd';
tsvdvelupdate = 1 ./ tsvdslowupdate;
tsvdresidual = times - jacobian * tsvdslowupdate;
meanres = mean(tsvdresidual);

tsvdsigma = sqrt(norm(tsvdresidual-meanres));

subplot(2,1,1)
plot2d(exactslowmodel)
plot2d(tsvdslowupdate)

subplot(2,1,2)
plot2d(times)
plot2d(jacobian * tsvdslowupdate)

subplot(2,1,1)
xstring(40,1.6,"exact vs computed model")
subplot(2,1,2)
xstring(5,24,"data vs computed response")


xbasimp(0,'vsp.ps',0
// this iwll make a file called vsp.ps.0
// make viewable eps file by doing
// BEpsf vsp.ps.0