1GMTMATH(1) Generic Mapping Tools GMTMATH(1)
2
3
4
6 gmtmath - Reverse Polish Notation calculator for data tables
7
9 gmtmath [ -At_f(t).d ] [ -Ccols ] [ -Fcols ] [ -H[i][nrec] ] [ -I ] [
10 -Nn_col/t_col ] [ -Q ] [ -S[f|l] ] [ -Tt_min/t_max/t_inc[+]|tfile ] [
11 -V ] [ -b[i|o][s|S|d|D[ncol]|c[var1/...]] ] [ -f[i|o]colinfo ] [
12 -m[i|o][flag] ] operand [ operand ] OPERATOR [ operand ] OPERATOR ... =
13 [ outfile ]
14
16 gmtmath will perform operations like add, subtract, multiply, and
17 divide on one or more table data files or constants using Reverse Pol‐
18 ish Notation (RPN) syntax (e.g., Hewlett-Packard calculator-style).
19 Arbitrarily complicated expressions may therefore be evaluated; the
20 final result is written to an output file [or standard output]. When
21 two data tables are on the stack, each element in file A is modified by
22 the corresponding element in file B. However, some operators only
23 require one operand (see below). If no data tables are used in the
24 expression then options -T, -N can be set (and optionally -b to indi‐
25 cate the data domain). If STDIN is given, <stdin> will be read and
26 placed on the stack as if a file with that content had been given on
27 the command line. By default, all columns except the "time" column are
28 operated on, but this can be changed (see -C).
29
30 operand
31 If operand can be opened as a file it will be read as an ASCII
32 (or binary, see -bi) table data file. If not a file, it is
33 interpreted as a numerical constant or a special symbol (see
34 below). The special argument STDIN means that stdin will be
35 read and placed on the stack; STDIN can appear more than once if
36 necessary.
37
38 outfile
39 The name of a table data file that will hold the final result.
40 If not given then the output is sent to stdout.
41
42 OPERATORS
43 Choose among the following 131 operators. "args" are the number
44 of input and output arguments.
45
46 Operator args Returns
47
48 ABS 1 1 abs (A).
49 ACOS 1 1 acos (A).
50 ACOSH 1 1 acosh (A).
51 ACOT 1 1 acot (A).
52 ACSC 1 1 acsc (A).
53 ADD 2 1 A + B.
54 AND 2 1 NaN if A and B == NaN, B if A == NaN, else A.
55 ASEC 1 1 asec (A).
56 ASIN 1 1 asin (A).
57 ASINH 1 1 asinh (A).
58 ATAN 1 1 atan (A).
59 ATAN2 2 1 atan2 (A, B).
60 ATANH 1 1 atanh (A).
61 BEI 1 1 bei (A).
62 BER 1 1 ber (A).
63 CEIL 1 1 ceil (A) (smallest integer >= A).
64 CHICRIT 2 1 Critical value for chi-squared-distribution, with
65 alpha = A and n = B.
66 CHIDIST 2 1 chi-squared-distribution P(chi2,n), with chi2 = A
67 and n = B.
68 COL 1 1 Places column A on the stack.
69 CORRCOEFF 2 1 Correlation coefficient r(A, B).
70 COS 1 1 cos (A) (A in radians).
71 COSD 1 1 cos (A) (A in degrees).
72 COSH 1 1 cosh (A).
73 COT 1 1 cot (A) (A in radians).
74 COTD 1 1 cot (A) (A in degrees).
75 CPOISS 2 1 Cumulative Poisson distribution F(x,lambda), with
76 x = A and lambda = B.
77 CSC 1 1 csc (A) (A in radians).
78 CSCD 1 1 csc (A) (A in degrees).
79 D2DT2 1 1 d^2(A)/dt^2 2nd derivative.
80 D2R 1 1 Converts Degrees to Radians.
81 DDT 1 1 d(A)/dt Central 1st derivative.
82 DILOG 1 1 dilog (A).
83 DIV 2 1 A / B.
84 DUP 1 2 Places duplicate of A on the stack.
85 EQ 2 1 1 if A == B, else 0.
86 ERF 1 1 Error function erf (A).
87 ERFC 1 1 Complementary Error function erfc (A).
88 ERFINV 1 1 Inverse error function of A.
89 EXCH 2 2 Exchanges A and B on the stack.
90 EXP 1 1 exp (A).
91 FACT 1 1 A! (A factorial).
92 FCRIT 3 1 Critical value for F-distribution, with alpha =
93 A, n1 = B, and n2 = C.
94 FDIST 3 1 F-distribution Q(F,n1,n2), with F = A, n1 = B,
95 and n2 = C.
96 FLIPUD 1 1 Reverse order of each column.
97 FLOOR 1 1 floor (A) (greatest integer <= A).
98 FMOD 2 1 A % B (remainder after truncated division).
99 GE 2 1 1 if A >= B, else 0.
100 GT 2 1 1 if A > B, else 0.
101 HYPOT 2 1 hypot (A, B) = sqrt (A*A + B*B).
102 I0 1 1 Modified Bessel function of A (1st kind, order
103 0).
104 I1 1 1 Modified Bessel function of A (1st kind, order
105 1).
106 IN 2 1 Modified Bessel function of A (1st kind, order
107 B).
108 INRANGE 3 1 1 if B <= A <= C, else 0.
109 INT 1 1 Numerically integrate A.
110 INV 1 1 1 / A.
111 ISNAN 1 1 1 if A == NaN, else 0.
112 J0 1 1 Bessel function of A (1st kind, order 0).
113 J1 1 1 Bessel function of A (1st kind, order 1).
114 JN 2 1 Bessel function of A (1st kind, order B).
115 K0 1 1 Modified Kelvin function of A (2nd kind, order
116 0).
117 K1 1 1 Modified Bessel function of A (2nd kind, order
118 1).
119 KEI 1 1 kei (A).
120 KER 1 1 ker (A).
121 KN 2 1 Modified Bessel function of A (2nd kind, order
122 B).
123 KURT 1 1 Kurtosis of A.
124 LE 2 1 1 if A <= B, else 0.
125 LMSSCL 1 1 LMS scale estimate (LMS STD) of A.
126 LOG 1 1 log (A) (natural log).
127 LOG10 1 1 log10 (A) (base 10).
128 LOG1P 1 1 log (1+A) (accurate for small A).
129 LOG2 1 1 log2 (A) (base 2).
130 LOWER 1 1 The lowest (minimum) value of A.
131 LRAND 2 1 Laplace random noise with mean A and std. devia‐
132 tion B.
133 LSQFIT 1 0 Let current table be [A | b]; return least
134 squares solution x = A \ b.
135 LT 2 1 1 if A < B, else 0.
136 MAD 1 1 Median Absolute Deviation (L1 STD) of A.
137 MAX 2 1 Maximum of A and B.
138 MEAN 1 1 Mean value of A.
139 MED 1 1 Median value of A.
140 MIN 2 1 Minimum of A and B.
141 MOD 2 1 A mod B (remainder after floored division).
142 MODE 1 1 Mode value (Least Median of Squares) of A.
143 MUL 2 1 A * B.
144 NAN 2 1 NaN if A == B, else A.
145 NEG 1 1 -A.
146 NEQ 2 1 1 if A != B, else 0.
147 NOT 1 1 NaN if A == NaN, 1 if A == 0, else 0.
148 NRAND 2 1 Normal, random values with mean A and std. devia‐
149 tion B.
150 OR 2 1 NaN if A or B == NaN, else A.
151 PLM 3 1 Associated Legendre polynomial P(A) degree B
152 order C.
153 PLMg 3 1 Normalized associated Legendre polynomial P(A)
154 degree B order C (geophysical convention).
155 POP 1 0 Delete top element from the stack.
156 POW 2 1 A ^ B.
157 PQUANT 2 1 The B'th Quantile (0-100%) of A.
158 PSI 1 1 Psi (or Digamma) of A.
159 PV 3 1 Legendre function Pv(A) of degree v = real(B) +
160 imag(C).
161 QV 3 1 Legendre function Qv(A) of degree v = real(B) +
162 imag(C).
163 R2 2 1 R2 = A^2 + B^2.
164 R2D 1 1 Convert Radians to Degrees.
165 RAND 2 1 Uniform random values between A and B.
166 RINT 1 1 rint (A) (nearest integer).
167 ROOTS 2 1 Treats col A as f(t) = 0 and returns its roots.
168 ROTT 2 1 Rotate A by the (constant) shift B in the t-
169 direction.
170 SEC 1 1 sec (A) (A in radians).
171 SECD 1 1 sec (A) (A in degrees).
172 SIGN 1 1 sign (+1 or -1) of A.
173 SIN 1 1 sin (A) (A in radians).
174 SINC 1 1 sinc (A) (sin (pi*A)/(pi*A)).
175 SIND 1 1 sin (A) (A in degrees).
176 SINH 1 1 sinh (A).
177 SKEW 1 1 Skewness of A.
178 SQR 1 1 A^2.
179 SQRT 1 1 sqrt (A).
180 STD 1 1 Standard deviation of A.
181 STEP 1 1 Heaviside step function H(A).
182 STEPT 1 1 Heaviside step function H(t-A).
183 SUB 2 1 A - B.
184 SUM 1 1 Cumulative sum of A.
185 TAN 1 1 tan (A) (A in radians).
186 TAND 1 1 tan (A) (A in degrees).
187 TANH 1 1 tanh (A).
188 TCRIT 2 1 Critical value for Student's t-distribution, with
189 alpha = A and n = B.
190 TDIST 2 1 Student's t-distribution A(t,n), with t = A, and
191 n = B.
192 TN 2 1 Chebyshev polynomial Tn(-1<A<+1) of degree B.
193 UPPER 1 1 The highest (maximum) value of A.
194 XOR 2 1 B if A == NaN, else A.
195 Y0 1 1 Bessel function of A (2nd kind, order 0).
196 Y1 1 1 Bessel function of A (2nd kind, order 1).
197 YN 2 1 Bessel function of A (2nd kind, order B).
198 ZCRIT 1 1 Critical value for the normal-distribution, with
199 alpha = A.
200 ZDIST 1 1 Cumulative normal-distribution C(x), with x = A.
201
202 SYMBOLS
203 The following symbols have special meaning:
204
205 PI 3.1415926...
206 E 2.7182818...
207 EULER 0.5772156...
208 TMIN Minimum t value
209 TMAX Maximum t value
210 TINC t increment
211 N The number of records
212 T Table with t-coordinates
213
215 -A Requires -N and will partially initialize a table with values
216 from the given file containing t and f(t) only. The t is placed
217 in column t_col while f(t) goes into column n_col - 1 (see -N).
218
219 -C Select the columns that will be operated on until next occur‐
220 rence of -C. List columns separated by commas; ranges like
221 1,3-5,7 are allowed. -C (no arguments) resets the default
222 action of using all columns except time column (see -N). -Ca
223 selects all columns, including time column, while -Cr reverses
224 (toggles) the current choices.
225
226 -F Give a comma-separated list of desired columns or ranges that
227 should be part of the output (0 is first column) [Default out‐
228 puts all columns].
229
230 -H Input file(s) has header record(s). If used, the default number
231 of header records is N_HEADER_RECS. Use -Hi if only input data
232 should have header records [Default will write out header
233 records if the input data have them]. Blank lines and lines
234 starting with # are always skipped.
235
236 -I Reverses the output row sequence from ascending time to descend‐
237 ing [ascending].
238
239 -N Select the number of columns and the column number that contains
240 the "time" variable. Columns are numbered starting at 0 [2/0].
241
242 -Q Quick mode for scalar calculation. Shorthand for -Ca -N1/0
243 -T0/0/1.
244
245 -S Only report the first or last row of the results [Default is all
246 rows]. This is useful if you have computed a statistic (say the
247 MODE) and only want to report a single number instead of numer‐
248 ous records with identical values. Append l to get the last row
249 and f to get the first row only [Default].
250
251 -T Required when no input files are given. Sets the t-coordinates
252 of the first and last point and the equidistant sampling inter‐
253 val for the "time" column (see -N). Append + if you are speci‐
254 fying the number of equidistant points instead. If there is no
255 time column (only data columns), give -T with no arguments; this
256 also implies -Ca. Alternatively, give the name of a file whose
257 first column contains the desired t-coordinates which may be
258 irregular.
259
260 -V Selects verbose mode, which will send progress reports to stderr
261 [Default runs "silently"].
262
263 -bi Selects binary input. Append s for single precision [Default is
264 d (double)]. Uppercase S or D will force byte-swapping.
265 Optionally, append ncol, the number of columns in your binary
266 input file if it exceeds the columns needed by the program. Or
267 append c if the input file is netCDF. Optionally, append
268 var1/var2/... to specify the variables to be read.
269
270 -bo Selects binary output. Append s for single precision [Default
271 is d (double)]. Uppercase S or D will force byte-swapping.
272 Optionally, append ncol, the number of desired columns in your
273 binary output file. [Default is same as input, but see -F]
274
275 -m Multiple segment file(s). Segments are separated by a special
276 record. For ASCII files the first character must be flag
277 [Default is '>']. For binary files all fields must be NaN and
278 -b must set the number of output columns explicitly. By default
279 the -m setting applies to both input and output. Use -mi and
280 -mo to give separate settings to input and output.
281
283 The ASCII output formats of numerical data are controlled by parameters
284 in your .gmtdefaults4 file. Longitude and latitude are formatted
285 according to OUTPUT_DEGREE_FORMAT, whereas other values are formatted
286 according to D_FORMAT. Be aware that the format in effect can lead to
287 loss of precision in the output, which can lead to various problems
288 downstream. If you find the output is not written with enough preci‐
289 sion, consider switching to binary output (-bo if available) or specify
290 more decimals using the D_FORMAT setting.
291
293 (1) The operators PLM and PLMg calculate the associated Legendre poly‐
294 nomial of degree L and order M in x which must satisfy -1 <= x <= +1
295 and 0 <= M <= L. x, L, and M are the three arguments preceding the
296 operator. PLM is not normalized and includes the Condon-Shortley phase
297 (-1)^M. PLMg is normalized in the way that is most commonly used in
298 geophysics. The C-S phase can be added by using -M as argument. PLM
299 will overflow at higher degrees, whereas PLMg is stable until ultra
300 high degrees (at least 3000).
301
302 (2) Files that have the same names as some operators, e.g., ADD, SIGN,
303 =, etc. should be identified by prepending the current directory (i.e.,
304 ./LOG).
305
306 (3) The stack depth limit is hard-wired to 100.
307
308 (4) All functions expecting a positive radius (e.g., LOG, KEI, etc.)
309 are passed the absolute value of their argument.
310
311 (5) The DDT and D2DT2 functions only work on regularly spaced data.
312
313 (6) All derivatives are based on central finite differences, with natu‐
314 ral boundary conditions.
315
316 (7) ROOTS must be the last operator on the stack, only followed by =.
317
319 To take the square root of the content of the second data column being
320 piped through gmtmath by process1 and pipe it through a 3rd process,
321 use
322
323 process1 | gmtmath STDIN SQRT = | process3
324
325 To take log10 of the average of 2 data files, use
326
327 gmtmath file1.d file2.d ADD 0.5 MUL LOG10 = file3.d
328
329 Given the file samples.d, which holds seafloor ages in m.y. and
330 seafloor depth in m, use the relation depth(in m) = 2500 + 350 * sqrt
331 (age) to print the depth anomalies:
332
333 gmtmath samples.d T SQRT 350 MUL 2500 ADD SUB = | lpr
334
335 To take the average of columns 1 and 4-6 in the three data sets
336 sizes.1, sizes.2, and sizes.3, use
337
338 gmtmath -C1,4-6 sizes.1 sizes.2 ADD sizes.3 ADD 3 DIV = ave.d
339
340 To take the 1-column data set ages.d and calculate the modal value and
341 assign it to a variable, try
342
343 set mode_age = `gmtmath -S -T ages.d MODE =`
344
345 To evaluate the dilog(x) function for coordinates given in the file
346 t.d:
347
348 gmtmath -Tt.d T DILOG = dilog.d
349
350 To use gmtmath as a RPN Hewlett-Packard calculator on scalars (i.e., no
351 input files) and calculate arbitrary expressions, use the -Q option.
352 As an example, we will calculate the value of Kei (((1 + 1.75)/2.2) +
353 cos (60)) and store the result in the shell variable z:
354
355 set z = `gmtmath -Q 1 1.75 ADD 2.2 DIV 60 COSD ADD KEI =`
356
357 To use gmtmath as a general least squares equation solver, imagine that
358 the current table is the augmented matrix [ A | b ] and you want the
359 least squares solution x to the matrix equation A * x = b. The opera‐
360 tor LSQFIT does this; it is your job to populate the matrix correctly
361 first. The -A option will facilitate this. Suppose you have a 2-col‐
362 umn file ty.d with t and b(t) and you would like to fit a the model
363 y(t) = a + b*t + c*H(t-t0), where H is the Heaviside step function for
364 a given t0 = 1.55. Then, you need a 4-column augmented table loaded
365 with t in column 0 and your observed y(t) in column 3. The calculation
366 becomes
367
368 gmtmath -N4/1 -Aty.d -C0 1 ADD -C2 1.55 STEPT ADD -Ca LSQFIT = solu‐
369 tion.d
370
371 Note we use the -C option to select which columns we are working on,
372 then make active all the columns we need (here all of them, with -Ca)
373 before calling LSQFIT. The second and fourth columns (col numbers 1
374 and 3) are preloaded with t and y(t), respectively, the other columns
375 are zero. If you already have a precalculated table with the augmented
376 matrix [ A | b ] in a file (say lsqsys.d), the least squares solution
377 is simply
378
379 gmtmath -T lsqsys.d LSQFIT = solution.d
380
382 Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical Func‐
383 tions, Applied Mathematics Series, vol. 55, Dover, New York.
384 Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the
385 Clenshaw summation and the recursive computation of very high degree
386 and order normalised associated Legendre functions. Journal of
387 Geodesy, 76, 279-299.
388 Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
389 1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.
390 Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere
391 Publishing Corp.
392
394 GMT(1), grdmath(1)
395
396
397
398GMT 4.5.6 10 Mar 2011 GMTMATH(1)