next up previous [pdf]

Next: Compiling Up: Fomel: RSF API Previous: Compiling

Fortran-90 interface

The Fortran-90 clip function is listed below.

program Clipit
  use rsf

  implicit none
  type (file)                      :: in, out
  integer                          :: n1, n2, i1, i2
  real                             :: clip
  real, dimension (:), allocatable :: trace

  call sf_init()            ! initialize RSF
  in = rsf_input()
  out = rsf_output()

  if (sf_float /= gettype(in)) call sf_error("Need floats")

  call from_par(in,"n1",n1)
  n2 = filesize(in,1)

  call from_par("clip",clip) ! command-line parameter 

  allocate (trace (n1))

  do i2=1, n2                ! loop over traces
     call rsf_read(in,trace)
     
     where (trace >  clip) trace =  clip
     where (trace < -clip) trace = -clip

     call rsf_write(out,trace)
  end do
end program Clipit

Let us examine it in detail.

  use rsf
The program starts with importing the rsf module.

  call sf_init()            ! initialize RSF
A call to sf_init is needed to initialize the command-line interface.

  in = rsf_input()
  out = rsf_output()
The standard input and output files are initialized with rsf_input and rsf_output functions. Both functions accept optional arguments. For example, if the command line contains vel=velocity.rsf, then both rsf_input("velocity.rsf") and rsf_input("vel") are acceptable.

  if (sf_float /= gettype(in)) call sf_error("Need floats")

A call to from_par extracts the ``n1'' parameter from the input file. Conceptually, the RSF data model is a multidimensional hypercube. The n1 parameter refers to the fastest axis. If the input dataset is a collection of traces, n1 corresponds to the trace length. We could proceed in a similar fashion, extracting n2, n3, etc. If we are interested in the total number of traces, like in the clip example, a shortcut is to use the filesize function. Calling filesize(in) returns the total number of elements in the hypercube (the product of n1, n2, etc.), calling filesize(in,1) returns the number of traces (the product of n2, n3, etc.), calling filesize(in,2) returns the product of n3, n4, etc. By calling filesize, we avoid the need to extract additional parameters for the hypercube dimensions that we are not interested in.

  n2 = filesize(in,1)
The clip parameter is read from the command line, where it can be specified, for example, as clip=10. If we knew a good default value for clip, we could specify it with an optional argument, i.e. call from_par("clip",clip,default).

  allocate (trace (n1))

  do i2=1, n2                ! loop over traces
     call rsf_read(in,trace)
     
     where (trace >  clip) trace =  clip
     where (trace < -clip) trace = -clip

Finally, we do the actual work: loop over input traces, reading, clipping, and writing out each trace.



Subsections
next up previous [pdf]

Next: Compiling Up: Fomel: RSF API Previous: Compiling

2013-04-08