Next: SPECULATION ON PILOT TRACE Up: Claerbout: Reflection tomography Previous: VALUE OF DATA EXTENSION

# CODING THE TOMOGRAPHY

I am using subroutine xtomo() .
```subroutine xtomo( adj,add,z0,dz,x0,dx,h0,dh,y0,dy,modl,nt,nz,nx,  data,   nh,ny)
real		                                  modl(nt,nz,nx), data(nt,nh,ny)
integer	                                               it,iz,ix,          ih,iy
real                      z0,dz,x0,dx,h0,dh,y0,dy
real		          z,    x,    h,    y,                zmax
do ix= 1, nx {  x= x0 + dx*(ix-1)
do ih= 1, nh {	h= h0 + dh*(ih-1)
do iz= 1, nz {	z= z0 + dz*(iz-1)

y = x +  h *  z/zmax    # maybe opposite sign on h or (1-z/zmax)

iy = 1.5 + (y-y0) / dy
if( 1 <= iy & iy <= ny)
do it= 1, nt
data(it,ih,iy) = data(it,ih,iy) + modl(it,iz,ix)
else
modl(it,iz,ix) = modl(it,iz,ix) + data(it,ih,iy)
}}}
return; end
```

The inversion problem in tomography was solved by using the conjugate-gradient subroutine cgstep() described in PVI invoked by subroutine tominv() . It allows for damping eps but I did not test damping.

```subroutine tominv( niter,modl,nt,nz,nx,  data,   nh,ny, z0,dz,x0,dx,h0,dh,y0,dy)
integer      iter, niter,     nt,nz,nx,          nh,ny
real	         modl(nt,nz,nx), data(nt,nh,ny), eps
real	                 z0,x0,          h0,y0
real	                 dz,dx,          dh,dy
integer                nmodl,           ndata,            ir
temporary real                             rr(nt*nh*ny + nt*nz*nx)
temporary real          dmodl(nt,nz,nx),   dr(nt*nh*ny + nt*nz*nx)
temporary real          smodl(nt,nz,nx),   sr(nt*nh*ny + nt*nz*nx)
nmodl= nt*nz*nx; ndata=nt*nh*ny;  ir= ndata + 1
call null(  modl, nmodl )
call copy(            ndata, data, rr( 1));	  eps= 0.  # Turn off.
call ident( 0,0, eps, nmodl, modl, rr(ir))
do iter= 0, niter {
call snap( 'rr.H',   nt*nh, ny, rr)			# snapshot
call xtomo( 1,0, z0,dz,x0,dx,h0,dh,y0,dy, dmodl, nt,nz,nx,rr( 1), nh,ny)
call ident( 1,1, eps, nmodl,              dmodl,          rr(ir)       )
call xtomo( 0,0, z0,dz,x0,dx,h0,dh,y0,dy, dmodl, nt,nz,nx,dr( 1), nh,ny)
call ident( 0,0, eps, nmodl,              dmodl,          dr(ir)       )
call cgstep( iter,       nmodl, modl, dmodl,smodl,
ndata+nmodl,   rr, dr,   sr    )
call snap( 'modl.H', nt*nz, nx, modl)			# snapshot
}
return; end
```

I omitted the listing of the main program because it merely serves to connect the listed subroutines to our file system.

Next: SPECULATION ON PILOT TRACE Up: Claerbout: Reflection tomography Previous: VALUE OF DATA EXTENSION
Stanford Exploration Project
11/16/1997