NAME

     grdmath - Reverse Polish Notation calculator for grd files


SYNOPSIS

     grdmath  [  -Ixinc[m|c][/yinc[m|c]]  -Rwest/east/south/north
     -V]  operand [ operand ] OPERATOR [ operand ] OPERATOR ... =
     outgrdfile


DESCRIPTION

     grdmath will perform operations like add,  subtract,  multi-
     ply,  and divide on one or more grd files or constants using
     Reverse Polish Notation (RPN) syntax (e.g.,  Hewlett-Packard
     calculator-style).   Arbitrarily complicated expressions may
     therefore be evaluated; the final result is  written  to  an
     output  grd  file. When two grd files are on the stack, each
     element in file A is modified by the  corresponding  element
     in file B.  However, some operators only require one operand
     (see below).  If no grdfiles are used in the expression then
     options -R, -I must be set (and optionally -F).

     operand
          If operand can be opened as a file it will be read as a
          grd file.  If not a file, it is interpreted as a numer-
          ical constant or a special symbol (see below).

     outgrdfile is a 2-
          D grd file that will hold the final result.

     OPERATORS
          Choose among the following operators:
          Operator       n_args    Returns

          ABS       1    abs (A).
          ACOS      1    acos (A).
          ACOSH          1    acosh (A).
          ADD(+)         2    A + B.
          AND       2    NaN if A and B == NaN, B if  A  ==  NaN,
          else A.
          ASIN      1    asin (A).
          ASINH          1    asinh (A).
          ATAN      1    atan (A).
          ATAN2          2    atan2 (A, B).
          ATANH     1    atanh (A).  ATANH     1    atanh (A).
          BEI       1    bei (A).
          BER       1    ber (A).
          CDIST          2    Cartesian  distance  between   grid
          nodes and stack x,y.
          CEIL      1    ceil (A) (smallest integer >= A).
          COS       1    cos (A) (A in radians).
          COSD      1    cos (A) (A in degrees).
          COSH      1    cosh (A).
          CURV      1    Curvature of A (Laplacian).
          D2DX2          1    d^2(A)/dx^2 2nd derivative.
          D2DY2          1    d^2(A)/dy^2 2nd derivative.
          D2R       1    Converts Degrees to Radians.
          DDX       1    d(A)/dx 1st derivative.
          DDY       1    d(A)/dy 1st derivative.
          DILOG          1    Dilog (A).
          DIV(/)         2    A / B.
          DUP       1    Places duplicate of A on the stack.
          ERF       1    Error function of A.
          ERFC      1    Complimentory Error function of A.
          ERFINV         1    Inverse error function of A.
          EXCH      2    Exchanges A and B on the stack.
          EXP       1    exp (A).
          FLOOR          1    floor (A) (greatest integer <= A).
          FMOD      2    A % B (remainder).
          GDIST          2    Great distance (in degrees) between
          grid nodes and stack lon,lat.
          HYPOT          2    hypot (A, B).
          I0        1    Modified Bessel function of A (1st kind,
          order 0).
          I1        1    Modified Bessel function of A (1st kind,
          order 1).
          IN        2    Modified Bessel function of A (1st kind,
          order B).
          INV       1    1 / A.
          J0        1    Bessel function of A  (1st  kind,  order
          0).
          J1        1    Bessel function of A  (1st  kind,  order
          1).
          JN        2    Bessel function of A  (1st  kind,  order
          B).
          K0        1    Modified Kelvin function of A (2nd kind,
          order 0).
          K1        1    Modified Bessel function of A (2nd kind,
          order 1).
          KN        2    Modified Bessel function of A (2nd kind,
          order B).
          KEI       1    kei (A).
          KER       1    ker (A).
          LOG       1    log (A) (natural log).
          LOG10          1    log10 (A).
          LOG1P          1    log (1+A) (accurate for small A).
          MAX       2    Maximum of A and B.
          MEAN      1    Mean value of A.
          MED       1    Median value of A.
          MIN       2    Minimum of A and B.
          MUL(x)         2    A * B.
          NEG       1    -A.
          OR        2    NaN if A or B == NaN, else A.
          PLM       3    Associated   Legendre   polynomial   P(-
          1<A<+1) degree B order C.
          POP       1    Delete top element from the stack.
          POW(^)         2    A ^ B.
          R2        2    R2 = A^2 + B^2.
          R2D       1    Convert Radians to Degrees.
          RINT      1    rint (A) (nearest integer).
          SIGN      1    sign (+1 or -1) of A.
          SIN       1    sin (A) (A in radians).
          SIND      1    sin (A) (A in degrees).
          SINH      1    sinh (A).
          SQRT      1    sqrt (A).
          STD       1    Standard deviation of A.
          STEPX          1    Heaviside step function in x:  H(x-
          A).
          STEPY          1    Heaviside step function in y:  H(y-
          A).
          SUB(-)         2    A - B.
          TAN       1    tan (A) (A in radians).
          TAND      1    tan (A) (A in degrees).
          TANH      1    tanh (A).
          Y0        1    Bessel function of A  (2nd  kind,  order
          0).
          Y1        1    Bessel function of A  (2nd  kind,  order
          1).
          YLM       2    Re and Im normalized  surface  harmonics
          (degree A, order B).
          YN        2    Bessel function of A  (2nd  kind,  order
          B).

     SYMBOLS
          The following symbols have special meaning:

          PI   3.1415926...
          E    2.7182818...
          X    Grid with x-coordinates
          Y    Grid with y-coordinates


OPTIONS

     -I   x_inc [and  optionally  y_inc]  is  the  grid  spacing.
          Append m to indicate minutes or c to indicate seconds.

     -R   west, east, south, and  north  specify  the  Region  of
          interest.  To specify boundaries in degrees and minutes
          [and seconds], use the dd:mm[:ss] format.  Append r  if
          lower  left  and  upper right map coordinates are given
          instead of wesn.

     -F   Select pixel registration.  [Default is grid  registra-
          tion].

     -V   Selects verbose mode, which will send progress  reports
          to stderr [Default runs "silently"].



BEWARE

     The operator GDIST calculates  spherical  distances  bewteen
     the  (lon, lat) point on the stack and all node positions in
     the grid.  The grid domain and  the  (lon,  lat)  point  are
     expected  to be in degrees.  The operator YLM calculates the
     fully normalized spherical harmonics for degree L and  order
     M  for  all positions in the grid, which is assumed to be in
     degrees.  YLM returns two grids, the Real (cosine) and  Ima-
     ginary  (sine)  component of the complex spherical harmonic.
     Use the POP operator (and EXCH) to get rid of one  of  them.
     The  operator PLM calculates the associated Legendre polyno-
     mial of degree L and order M, and its argument is the cosine
     of  the  colatitude which must satisfy -1 <= x <= +1. Unlike
     YLM, PLM is not normalized.
     All the derivatives are based on central finite differences,
     with natural boundary conditions.


EXAMPLES

     To take log10 of the average of 2 files, use
          grdmath  file1.grd  file2.grd  ADD  0.5  MUL  LOG10   =
     file3.grd

     Given the file ages.grd, which holds seafloor ages in  m.y.,
     use  the  relation  depth(in m) = 2500 + 350 * sqrt (age) to
     estimate normal seafloor depths:
          grdmath ages.grd SQRT 350 MUL 2500 ADD = depths.grd

     To find the angle a (in degrees) of  the  largest  principal
     stress  from  the  stress  tensor  given  by the three files
     s_xx.grd s_yy.grd, and s_xy.grd from the relation tan  (2*a)
     = 2 * s_xy / (s_xx - s_yy), try
          grdmath 2 s_xy.grd MUL s_xx.grd s_yy.grd SUB DIV  ATAN2
     2 DIV = direction.grd

     To calculate the  fully  normalized  spherical  harmonic  of
     degree 8 and order 4 on a 1 by 1 degree world map, using the
     real amplitude 0.4 and the imaginary amplitude 1.1, try
          grdmath -R0/360/-90/90 -I1 8 4 YML 1.1 MUL EXCH 0.4 MUL
     ADD = harm.grd


BUGS

     Files that has the same name as some operators,  e.g.,  ADD,
     SIGN,  =, etc. cannot be read and must not be present in the
     current directory.  Piping of files are  not  allowed.   The
     stack limit is hard-wired to 50.  Bessel and error functions
     may not be available  on  all  systems.   The  Kelvin-Bessel
     functions  (bei,  ber, kei, ker) are based on the polynomial
     approximations by Abramowitz and Stegun for  r  <=  8.   All
     functions expecting a positive radius (e.g., log, kei, etc.)
     are passed the absolute value of their argument.



REFERENCES

     Abramowitz,  M.,  and  I.  A.  Stegun,  1964,  Handbook   of
     Mathematical Functions, Applied Mathematics Series, vol. 55,
     Dover, New York.
     Press, W. H.,  S. A. Teukolsky,  W.  T.  Vetterling,  B.  P.
     Flannery,  1992,  Numerical  Recipes, 2nd edition, Cambridge
     Univ., New York.


SEE ALSO

     gmt(l),  gmtmath(l),  grd2xyz(l),  grdedit(l),   grdinfo(l),
     xyz2grd(l)