APPENDIX

Several Mathematica packages and run streams are listed here as samples for seismic processing and display.

The following routine generates the Hilbert transform of a trace.

BeginPackage["hilbert`"]

Hilbert::usage =
             "Hilbert- create the Hilbert transform of a trace."
Begin["`Private`"]
Hilbert::Badarg = " Bad input to Hilbert"

Hilbert[x_] := (
 Block[
        {i,half,junk123},
        half= Floor[ N[ Dimensions[x][[1]]/2 +1 ] ];
        junk123 =RotateRight[ Fourier[N[x]],half-1 ];
        Re[InverseFourier[
                         RotateLeft[
                          Table[
                                +Sign[ i-half ] I junk123[[i]],
                                {i,1,Dimensions[x][[1]]}]
                            ,half-1]
                         ]
          ]
      ]
                )
End[]
EndPackage[]

The run stream used to generate the Hilbert transforms in the paper follows. Note that the YPrint.m package produces Postscript files on disk for use in this paper. The GraphicsArray function is used to get around what appears to be a bug in the Sun system's display of these Postscript files. The DEC and IBM do not appear to be troubled by the output not processed with GraphicsArray.

ia = <<WSstack.m;
<<YPrint.m
<<~/amath/Hilbert.m
<<~/amath/Wiggle.m
<<Graphics`Graphics3D`
<<Graphics`PlotField`
a = ia;
ah = Map[Hilbert,a];
env = N[Sqrt[ah^2+a^2]];
p1 = GraphicsArray[{ListDensityPlot[a,Mesh->False]}]
YPrint[p1]
!mv out.ps origden.ps
p2 = GraphicsArray[{ListDensityPlot[ah,Mesh->False]}]
YPrint[p2]
!mv out.ps hilbertden.ps
p3 = GraphicsArray[{ListDensityPlot[env,Mesh->False]}]
YPrint[p3]
!mv out.ps envden.ps

The wiggle trace plot and the vector plot were generated with the following run.

<<YPrint.m
<<~/amath/Hilbert.m
<<~/amath/Wiggle.m
<<Graphics`PlotField`

ia = <<WSstack.m;
a = Table[ia[[i,j]],{i,1,60,2},{j,31,90,2}];
ah = Map[Hilbert,a];
env = N[Sqrt[ah^2+a^2]];

p5 = Wiggle[ia,2 MRMS[ia]]
YPrint[p5]
!mv out.ps wiggle.ps

t = Table[{a[[i,j]],ah[[i,j]]},{i,30},{j,30}];
p4 = ListPlotVectorField[t,ScaleFactor->1.0]
YPrint[p4]
!mv out.ps vector.ps

The Fourier transform of the data in Figure 6 was generated with the following run stream.

<<YPrint.m

a = <<WSstack.m;
bb = Fourier[a];
b=N[ Table[ bb[[i,j]],{i,1,Dimensions[bb][[1]]},{j,1,Dimensions[bb][[2]]/2} ] ];
p1 = GraphicsArray[{ListDensityPlot[Re[N[b]],Mesh->False]}]
YPrint[p1]
!mv out.ps rs.ps
p2 = GraphicsArray[{ListDensityPlot[Im[N[b]],Mesh->False]}]
YPrint[p2]
!mv out.ps imag.ps
p3 = GraphicsArray[{ListDensityPlot[Abs[N[b]],Mesh->False]}]
YPrint[p3]
!mv out.ps abs.ps
p4 = GraphicsArray[{ListPlot3D[Abs[N[b]],Mesh->False,PlotRange->All]}]
YPrint[p4]
!mv out.ps FT3D.ps

A small package that allows a plot to be output to disk in PostScript format is the YPrint.m routine listed below. It outputs the plot as a file out.ps that can be renamed from within Mathematica as seen in the previous example.

System`YPrint::usage =
             "YPrint[-graphics-] sends graphics to out.ps file."
Begin["System`Private`"]
YPrint[x_] := ( Display["!psfix > out.ps", x]; x )
End[]

Finally a FORTRAN program to convert SEP files to Mathematica format is given here. Perhaps a better method of converting the numbers would be to built the output in terms of powers of ten, to avoid the many format statements.

# Convert SEP format to Mathematica format
#
# Main.x head=/dev/null < input.H | cat outfile
#
#%
integer nt, nx
from history:    integer n1 : nt
from history:    integer n2 : nx
allocate:       integer gather(nt,nx)

subroutine doit( nt,nx, gather )
integer it,in,nt, nx
real gather(nt,nx),rnum

#        ----   call hclose()
100     format(a8)

        call sreed( 'in', gather, nt*nx*4)
        write(*,200) '{'
200     format(a)
        do in=1,nx
        {
          write(*,200) '{'
          do it=1,nt
          {
#          write(*,*) gather(it,in)
#          print *, gather(it,in)
           rnum = gather(it,in)
#           write(*,300) gather(it,in)
           if (abs(rnum) .lt. 1)
            write(*,298) rnum
           else if (abs(rnum) .lt. 1000)
           write(*,299) rnum
           else if (abs(rnum) .lt. 1000000)
            write(*,300)  rnum
           else if (abs(rnum) .lt. 1.E9)
                write(*,301)  rnum
           else if (abs(rnum) .lt. 1.E12 )
                write(*,302)  rnum
           else if (abs(rnum) .lt. 1.E16 )
                write(*,303)  rnum
           else if (abs(rnum) .lt. 1.E20 )
                write(*,304)  rnum
           else if (abs(rnum) .lt. 1.E25 )
                write(*,305)  rnum
           else if (abs(rnum) .lt. 1.E30 )
                write(*,306)  rnum
                else {
                       write(*,306)  sign(1.E30,rnum)
                     }
298        format(f14.10)
299        format(f14.7)
300        format(f14.4)
301        format(f17.4)
302        format(f20.4)
303        format(f24.0)
304        format(f28.0)
305        format(f32.0)
306        format(f34.0)

           if (it .ne. nt) write(*,200) ','
          }
          write(*,200) '}'
           if (in .ne. nx) write(*,200) ','
        }
        write(*,200) '}'

#       call srite( 'out', gather, nt*nx*4)
return; end