gmtmath

贡献者:

周茂

官方文档:

gmtmath

简介:

表数据的逆波兰表示法（RPN）计算

gmtmath 使用逆波兰表示法（Reverse Polish Notation）对一个或多个表数据文件或者常量进行加减乘除等操作。逆波兰表示法是一种后缀表示法，即将运算符写在操作数之后。该模块可以计算任意复杂的运算。模块默认逐元素进行计算操作。其中某些运算符仅需要一个参数。如果输入中不含表文件，则可以使用 -T 或 -N 选项。默认情况下，除了时间以外的列，所有的列都参与运算（见 -C ）。如果多个运算整体非常复杂或者频繁出现，可以编写为宏以便后续使用。

语法

gmt math [ -At_f(t)[+e][+r][+s|w] ] [ -Ccols ] [ -Eeigen ] [ -I ] [ -Nn_col[/t_col] ] [ -Q[c|i|p|n] ] [ -S[f|l] ] [ -T[min/max/inc[+b|i|l|n]|file|list] ] [ -V[level] ] [ -bbinary ] [ -dnodata[+ccol] ] [ -eregexp ] [ -fflags ] [ -ggaps ] [ -hheaders ] [ -iflags ] [ -oflags ] [ -qflags ] [ -sflags ] [ -wflags ] [ --PAR=value ] operand [ operand ] OPERATOR [ operand ] OPERATOR … = [ outfile ]

必选选项

operand

如果 operand 是一个文件，GMT 将试着以表文件的形式读取该文件。如果非文件，将被解析为 其他符号 。 若为 STDIN ，表示将读取标准输入并放入堆栈中。如果有必要， STDIN 参数可以多次出现。

= outfile

输出表数据文件名，不指定的情况下写到标准输出

可选选项

-At_f(t)[+e][+r][+s|w]

建立并求解线性方程。该选项需要 -Nn_col[/t_col] 选项，用来从给定的 t_f(t) 文件中进行初始化解方程所需的增广矩阵，其中 t_f(t) 文件只包含 t 和 f(t) 两列。在增广矩阵中，t 变量于编号为 t_col 指定的列， f(t) 位于编号为 n_col - 1 的列。

  • +e 使用 LSQFIT 或 SVDFIT 操作符时，使用该选项，可以求解拟合系数。结果包括 4 列：t，f(t)，t 位置的计算值，t 位置的残差 [默认只输出系数]

  • +r 选项只输出 f(t)

  • +s 指定 t_f(t) 的第三列为 1 倍中误差

  • +w 指定 t_f(t) 的第三列为权重

具体使用方法见倒数第二个示例。

-Ccols

选择要运算的列 cols 。多列以逗号分隔，类似 1,3-5,7 的表达形式也是允许的；此外，-Cx 等效于 -C0，-Cy 等效于 -C1。使用不带任何参数的 -C 选项可用来重置，即选择除时间列外的所有列（见 -N ）。-Ca 选择所有的列，包括时间列， -Cr 的作用是反向选择，即选择除当前已选定的列以外的其他列。

-Eeigen

设置 LSQFIT 和 SVDFIT 运算符的特征值，默认为 1e-7。比 1e-7 更小的特征值将被重置为 0 ，并在求解中不考虑。

-I

将输出数据从升序变为降序 [默认为升序]

-Nn_col[/t_col]

设置列数和时间所在的列，列从 0 开始计数，默认输入列最大编号 n_clo 为 2 ，时间列 t_col 位于第 0 列。

-Q[c|i|p|n]

实现不同长度单位之间的快速计算， -Ca -N1/0 -T0/0/1 的简写。在进行不同单位的长度计算时，在计算前都转换为 inch ，然后计算，并将结果转换为 PROJ_LENGTH_UNIT 设置的单位。使用 c ， i ， p 可设置对应的输出长度单位。 n 表示输出结果无量纲。

-S[f|l]

仅输出结果的第一行 [默认为所有行]。这在输出统计结果的时候可能会有用（比如 MODE 运算）。 l 和 f 分别表示第一行和最后一行。

-T[min/max/inc[+b|i|l|n]|file|list]

不加任何选项和文件，表明没有输入文件。其他详见 生成一维数组

-V[level] (more …)

设置 verbose 等级 [w]

-bi[ncols][type][w][+l|b] (more …)

设置二进制输入数据的格式

-bo[ncols][type][w][+l|b] (more …)

设置二进制输出的数据格式

-d[i|o]nodata (more …)

将输入数据中等于 nodata 的记录替换为 NaN，或将输出数据中值为 NaN 的记录替换为 nodata

-e[~]“pattern” | -e[~]/regexp/[i] (more …)

筛选或剔除匹配指定模式的数据记录

-f[i|o]colinfo (more …)

指定输入或输出列的数据类型

-g[a]x|y|d|X|Y|D|[col]zgap[+n|p] (more …)

确定数据或线段的间断

-h[i|o][n][+c][+d][+msegheader][+rremark][+ttitle] (more …)

跳过或生成指定数目的头段记录

-icols[+l][+sscale][+ooffset][,…][,t[word]] (more …)

设置输入数据列及简单变换（0表示第一列，t 表示文本列）

-ocols[,…][,t[word]] (more …)

设置输出数据列（0表示第一列，t 表示文本列）

-q[i|o][~]rows[+ccol][+a|f|s] (more …)

筛选输入或输出的行或数据范围

-s[cols][+a|+r] (more …)

设置 NaN 记录的处理方式

-wy|a|w|d|h|m|s|cperiod[/phase][+ccol] (more …)

将输入坐标转换为循环坐标

-^ 或 -

显示简短的帮助信息，包括模块简介和基本语法信息（Windows下只能使用 -）

-+ 或 +

显示帮助信息，包括模块简介、基本语法以及模块特有选项的说明

-? 或无参数

显示完整的帮助信息，包括模块简介、基本语法以及所有选项的说明

--PAR=value

临时修改GMT参数的值，可重复多次使用。参数列表见 配置参数

生成一维数组

下面将展示如何使用 gmtmath 生成一维数组（其中大部分操作也可通过 linux 中的 seq 命令方便地实现）

不使用任何子选项

以 0.1 为步长，生成 3.1 到 4.2 等等间隔序列

gmt math -o0 -T3.1/4.2/0.1 T =
3.1
3.2
...
4.2

+a 选项

该选项不对生成的数列进行运算操作，而是将该数列以列的形式追加到输出表数据，类似 linux 中的 paste 命令。

+b 选项

以 3 和 20 分别为起点和终点，创建一个 2 的整数幂的列表，将生成的整数幂以 1 为步长取值 -T3/20/1+b

gmt math -o0 -T3/20/1+b T =
4
8
16

+l 选项

以 7 和 135 分别为起点和终点，先创建一个 10 的整数幂的列表，当 inc 为 1 时，输出该列表；当 inc 为 2 时，同时输出该列表中每个元素的 1，2，5 倍值（不包括超出原始数据范围的部分）； 当 inc 为 3 时，同时输出该列表中每个元素的 1，2，…，9 倍值（同样不包括超出原始数据范围的部分）。 -T7/135/2+l 将会生成如下数列

gmt math -o0 -T7/135/2+l T =
10
20
50
100

inc 为负整数时，将实现如下效果

gmt math -o0 -T1e-4/1e4/-2+l T =
0.0001
0.01
1
100
10000

+i 选项

该选项将以 1 作为默认步长，第三个数字做为每步长中的数字个数 length。假设需要在 1 分钟内每 24 秒生成一帧

gmt math -o0 -T0/60/24+i T =
0
0.0416666666667
0.0833333333333
0.125
0.166666666667
...
60

+n 选项

生成固定长度的数列。以 3.44 和 7.82 为起点和终点，生成长度为 5 的等间隔序列

gmt math -o0 -T3.44/7.82/5+n T =
3.44
4.535
5.63
6.725
7.82

-T 后直接加文件或逗号分隔的数列

-T 后可以直接加文件，文件中即为要生成的列表；另外，可以直接使用逗号分隔，将列表附加到 -T 选项后，如下为 Fibonacci 数列前 6 项

gmt math -o0 -T0,1,1,2,3,5 T =
0
1
1
2
3
5

注 ： 如果数列只包含一个值，必须在其后加逗号以表明仍是一个数列

+u 选项

如果文件或者逗号分隔的数列中有重复数字或未排序，可以使用 +u 选项去重并排序。

+t 选项

生成绝对时间序列。在 inc 后分别添加 y ，o ，d ，h ， m 和 s 表示时间步长的单位为年，月，日，时，分，秒。在其后附加 +t 选项，可以进一步强调生成时间序列，也可以不加

gmt math -o0 -T2020-03-01T/2020-03-07T/1d T =
2020-03-01T00:00:00
2020-03-02T00:00:00
2020-03-03T00:00:00
2020-03-04T00:00:00
2020-03-05T00:00:00
2020-03-06T00:00:00
2020-03-07T00:00:00

生成距离序列

如果输入文件中包含两列以上的数据，可以使用前两列计算距离，并生成等距序列。在 inc 后分别添加 d，m，s，e，f，M，n 和 u 表示距离步长的单位为度，分，秒，米，英尺，公里，英里，海里，英尺。如果为笛卡尔坐标，使用 c 作为距离步长单位。

+e 选项

如果只给定 inc 而从数据中获取最大值和最小值，则 (max - min)/inc 可能不是整数，GMT 讲会自动对 inc 进行一定的调整。如果不想调整 inc ，则可以使用 +e 选项，GMT 会固定最小值，适当调整最大值。

运算符

Operator 为运算符名称； args 分别表示输入参数和输出参数个数； Returns 中的 A B C … 等字符表示输入参数。所有的三角函数运算符默认输入为弧度，若运算符最后一个字母为 D 表示使用角度计算。

Operator

args

Returns

ABS

1 1

abs (A)

ACOS

1 1

acos (A)

ACOSD

1 1

acosd (A)

ACOSH

1 1

acosh (A)

ACSC

1 1

acsc (A)

ACSCD

1 1

acscd (A)

ACOT

1 1

acot (A)

ACOTD

1 1

acotd (A)

ADD

2 1

A + B

AND

2 1

B if A == NaN, else A

ASEC

1 1

asec (A)

ASECD

1 1

asecd (A)

ASIN

1 1

asin (A)

ASIND

1 1

asind (A)

ASINH

1 1

asinh (A)

ATAN

1 1

atan (A)

ATAND

1 1

atand (A)

ATAN2

2 1

atan2 (A, B)

ATAN2D

2 1

atan2d (A, B)

ATANH

1 1

atanh (A)

BCDF

3 1

Binomial cumulative distribution function for p = A, n = B, and x = C

BPDF

3 1

Binomial probability density function for p = A, n = B, and x = C

BEI

1 1

Kelvin function bei (A)

BER

1 1

Kelvin function ber (A)

BITAND

2 1

A & B (bitwise AND operator)

BITLEFT

2 1

A << B (bitwise left-shift operator)

BITNOT

1 1

~A (bitwise NOT operator, i.e., return two’s complement)

BITOR

2 1

A | B (bitwise OR operator)

BITRIGHT

2 1

A >> B (bitwise right-shift operator)

BITTEST

2 1

1 if bit B of A is set, else 0 (bitwise TEST operator)

BITXOR

2 1

A ^ B (bitwise XOR operator)

CEIL

1 1

ceil (A) (smallest integer >= A)

CHICRIT

2 1

Chi-squared distribution critical value for alpha = A and nu = B

CHICDF

2 1

Chi-squared cumulative distribution function for chi2 = A and nu = B

CHIPDF

2 1

Chi-squared probability density function for chi2 = A and nu = B

COL

1 1

Places column A on the stack

COMB

2 1

Combinations n_C_r, with n = A and r = B

CORRCOEFF

2 1

Correlation coefficient r(A, B)

COS

1 1

cos (A) (A in radians)

COSD

1 1

cos (A) (A in degrees)

COSH

1 1

cosh (A)

COT

1 1

cot (A) (A in radians)

COTD

1 1

cot (A) (A in degrees)

CSC

1 1

csc (A) (A in radians)

CSCD

1 1

csc (A) (A in degrees)

DDT

1 1

d(A)/dt Central 1st derivative

D2DT2

1 1

d^2(A)/dt^2 2nd derivative

D2R

1 1

Converts degrees to radians

DEG2KM

1 1

Converts spherical degrees to kilometers

DENAN

2 1

Replace NaNs in A with values from B

DILOG

1 1

dilog (A)

DIFF

1 1

Forward difference between adjacent elements of A (A[1]-A[0], A[2]-A[1], …, NaN)

DIV

2 1

A / B

DUP

1 2

Places duplicate of A on the stack

ECDF

2 1

Exponential cumulative distribution function for x = A and lambda = B

ECRIT

2 1

Exponential distribution critical value for alpha = A and lambda = B

EPDF

2 1

Exponential probability density function for x = A and lambda = B

ERF

1 1

Error function erf (A)

ERFC

1 1

Complementary Error function erfc (A)

ERFINV

1 1

Inverse error function of A

EQ

2 1

1 if A == B, else 0

EXCH

2 2

Exchanges A and B on the stack

EXP

1 1

exp (A)

FACT

1 1

A! (A factorial)

FCDF

3 1

F cumulative distribution function for F = A, nu1 = B, and nu2 = C

FCRIT

3 1

F distribution critical value for alpha = A, nu1 = B, and nu2 = C

FLIPUD

1 1

Reverse order of each column

FLOOR

1 1

floor (A) (greatest integer <= A)

FMOD

2 1

A % B (remainder after truncated division)

FPDF

3 1

F probability density function for F = A, nu1 = B, and nu2 = C

GE

2 1

1 if A >= B, else 0

GT

2 1

1 if A > B, else 0

HSV2LAB

3 3

Convert h,s,v triplets to l,a,b triplets, with h = A (0-360), s = B and v = C (0-1)

HSV2RGB

3 3

Convert h,s,v triplets to r,g,b triplets, with h = A (0-360), s = B and v = C (0-1)

HSV2XYZ

3 3

Convert h,s,v triplets to x,t,z triplets, with h = A (0-360), s = B and v = C (0-1)

HYPOT

2 1

hypot (A, B) = sqrt (A*A + B*B)

I0

1 1

Modified Bessel function of A (1st kind, order 0)

I1

1 1

Modified Bessel function of A (1st kind, order 1)

IFELSE

3 1

B if A != 0, else C

IN

2 1

Modified Bessel function of A (1st kind, order B)

INRANGE

3 1

1 if B <= A <= C, else 0

INT

1 1

Numerically integrate A

INV

1 1

1 / A

ISFINITE

1 1

1 if A is finite, else 0

ISNAN

1 1

1 if A == NaN, else 0

J0

1 1

Bessel function of A (1st kind, order 0)

J1

1 1

Bessel function of A (1st kind, order 1)

JN

2 1

Bessel function of A (1st kind, order B)

K0

1 1

Modified Kelvin function of A (2nd kind, order 0)

K1

1 1

Modified Bessel function of A (2nd kind, order 1)

KM2DEG

1 1

Converts kilometers to spherical degrees

KN

2 1

Modified Bessel function of A (2nd kind, order B)

KEI

1 1

Kelvin function kei (A)

KER

1 1

Kelvin function ker (A)

KURT

1 1

Kurtosis of A

LAB2HSV

3 3

Convert l,a,b triplets to h,s,v triplets

LAB2RGB

3 3

Convert l,a,b triplets to r,g,b triplets

LAB2XYZ

3 3

Convert l,a,b triplets to x,y,z triplets

LCDF

1 1

Laplace cumulative distribution function for z = A

LCRIT

1 1

Laplace distribution critical value for alpha = A

LE

2 1

1 if A <= B, else 0

LMSSCL

1 1

LMS (Least Median of Squares) scale estimate (LMS STD) of A

LMSSCLW

2 1

Weighted LMS scale estimate (LMS STD) of A for weights in B

LOG

1 1

log (A) (natural log)

LOG10

1 1

log10 (A) (base 10)

LOG1P

1 1

log (1+A) (accurate for small A)

LOG2

1 1

log2 (A) (base 2)

LOWER

1 1

The lowest (minimum) value of A

LPDF

1 1

Laplace probability density function for z = A

LRAND

2 1

Laplace random noise with mean A and std. deviation B

LSQFIT

1 0

Let current table be [A | b] return least squares solution x = A \ b

LT

2 1

1 if A < B, else 0

MAD

1 1

Median Absolute Deviation (L1 STD) of A

MADW

2 1

Weighted Median Absolute Deviation (L1 STD) of A for weights in B

MAX

2 1

Maximum of A and B

MEAN

1 1

Mean value of A

MEANW

2 1

Weighted mean value of A for weights in B

MEDIAN

1 1

Median value of A

MEDIANW

2 1

Weighted median value of A for weights in B

MIN

2 1

Minimum of A and B

MOD

2 1

A mod B (remainder after floored division)

MODE

1 1

Mode value (Least Median of Squares) of A

MODEW

2 1

Weighted mode value (Least Median of Squares) of A for weights in B

MUL

2 1

A * B

NAN

2 1

NaN if A == B, else A

NEG

1 1

-A

NEQ

2 1

1 if A != B, else 0

NORM

1 1

Normalize (A) so max(A)-min(A) = 1

NOT

1 1

NaN if A == NaN, 1 if A == 0, else 0

NRAND

2 1

Normal, random values with mean A and std. deviation B

OR

2 1

NaN if B == NaN, else A

PCDF

2 1

Poisson cumulative distribution function for x = A and lambda = B

PERM

2 1

Permutations n_P_r, with n = A and r = B

PPDF

2 1

Poisson distribution P(x,lambda), with x = A and lambda = B

PLM

3 1

Associated Legendre polynomial P(A) degree B order C

PLMg

3 1

Normalized associated Legendre polynomial P(A) degree B order C (geophysical convention)

POP

1 0

Delete top element from the stack

POW

2 1

A ^ B

PQUANT

2 1

The B’th quantile (0-100%) of A

PQUANTW

3 1

The C’th weighted quantile (0-100%) of A for weights in B

PSI

1 1

Psi (or Digamma) of A

PV

3 1

Legendre function Pv(A) of degree v = real(B) + imag(C)

QV

3 1

Legendre function Qv(A) of degree v = real(B) + imag(C)

R2

2 1

R2 = A^2 + B^2

R2D

1 1

Convert radians to degrees

RAND

2 1

Uniform random values between A and B

RCDF

1 1

Rayleigh cumulative distribution function for z = A

RCRIT

1 1

Rayleigh distribution critical value for alpha = A

RGB2HSV

3 3

Convert r,g,b triplets to h,s,v triplets, with r = A, g = B, and b = C (in 0-255 range)

RGB2LAB

3 3

Convert r,g,b triplets to l,a,b triplets, with r = A, g = B, and b = C (in 0-255 range)

RGB2XYZ

3 3

Convert r,g,b triplets to x,y,x triplets, with r = A, g = B, and b = C (in 0-255 range)

RINT

1 1

rint (A) (round to integral value nearest to A)

RMS

1 1

Root-mean-square of A

RMSW

1 1

Weighted root-mean-square of A for weights in B

RPDF

1 1

Rayleigh probability density function for z = A

ROLL

2 0

Cyclicly shifts the top A stack items by an amount B

ROTT

2 1

Rotate A by the (constant) shift B in the t-direction

SEC

1 1

sec (A) (A in radians)

SECD

1 1

sec (A) (A in degrees)

SIGN

1 1

sign (+1 or -1) of A

SIN

1 1

sin (A) (A in radians)

SINC

1 1

sinc (A) (sin (pi*A)/(pi*A))

SIND

1 1

sin (A) (A in degrees)

SINH

1 1

sinh (A)

SKEW

1 1

Skewness of A

SQR

1 1

A^2

SQRT

1 1

sqrt (A)

STD

1 1

Standard deviation of A

STDW

2 1

Weighted standard deviation of A for weights in B

STEP

1 1

Heaviside step function H(A)

STEPT

1 1

Heaviside step function H(t-A)

SUB

2 1

A - B

SUM

1 1

Cumulative sum of A

TAN

1 1

tan (A) (A in radians)

TAND

1 1

tan (A) (A in degrees)

TANH

1 1

tanh (A)

TAPER

1 1

Unit weights cosine-tapered to zero within A of end margins

TN

2 1

Chebyshev polynomial Tn(-1<A<+1) of degree B

TCRIT

2 1

Student’s t distribution critical value for alpha = A and nu = B

TPDF

2 1

Student’s t probability density function for t = A, and nu = B

TCDF

2 1

Student’s t cumulative distribution function for t = A, and nu = B

UPPER

1 1

The highest (maximum) value of A

VAR

1 1

Variance of A

VARW

2 1

Weighted variance of A for weights in B

VPDF

3 1

Von Mises density distribution V(x,mu,kappa), with angles = A, mu = B, and kappa = C

WCDF

3 1

Weibull cumulative distribution function for x = A, scale = B, and shape = C

WCRIT

3 1

Weibull distribution critical value for alpha = A, scale = B, and shape = C

WPDF

3 1

Weibull density distribution P(x,scale,shape), with x = A, scale = B, and shape = C

XOR

2 1

B if A == NaN, else A

XYZ2HSV

3 3

Convert x,y,z triplets to h,s,v triplets

XYZ2LAB

3 3

Convert x,y,z triplets to l,a,b triplets

XYZ2RGB

3 3

Convert x,y,z triplets to r,g,b triplets

Y0

1 1

Bessel function of A (2nd kind, order 0)

Y1

1 1

Bessel function of A (2nd kind, order 1)

YN

2 1

Bessel function of A (2nd kind, order B)

ZCDF

1 1

Normal cumulative distribution function for z = A

ZPDF

1 1

Normal probability density function for z = A

ZCRIT

1 1

Normal distribution critical value for alpha = A

ROOTS

2 1

Treats col A as f(t) = 0 and returns its roots

其他符号

以下符号具有特殊的含义：

PI

3.1415926…

E

2.7182818…

EULER

0.5772156…

PHI

1.6180339… (golden ratio)

EPS_F

1.192092896e-07 (sgl. prec. eps)

EPS_D

2.2204460492503131e-16 (dbl. prec. eps)

TMIN

Minimum t value

TMAX

Maximum t value

TRANGE

Range of t values

TINC

t increment

N

The number of records

T

Table with t-coordinates

TNORM

Table with normalized t-coordinates

TROW

Table with row numbers 1, 2, …, N-1

上述符号均可以作为变量使用，当其为多个数时，逐元素操作。以 生成一维数组 中第一个实例中的 T 符号为例，其含义为将 -T 选项生成的序列转换为表，以便后续使用 = 输出。

运算符注意事项

  1. PLM 和 PLMg 运算符用来计算 L 阶 M 次 x 的缔合勒让德函数；各参数的范围应该满足 -1 <= x <= +1 、0 <= M <= L。 PLM 运算符没有经过标准化，并且乘以 phase (-1)^M。 PLMg 使用大地测量/地球物理常见的标准化。使用 -M 参数可以附加输出球谐系数 C 和 S。 PLM 在较高的阶次就会出现溢出，具体和纬度相关， PLMg 则可以保证在 3000 阶以下都不会溢出

  2. 不同参数文件名相同时，应给定相对路径或绝对路径以区分不同文件

  3. 该模块计算过程保存在堆栈中，栈中保存结果最大为 100，即不能叠加太多的操作以防溢出

  4. 所有需要半径的运算符为保证其为正值，均自动取绝对值后计算

  5. DDT 和 DDT2 函数仅适用于等间隔数据

  6. 所有导数都基于 central finite differences 和 natural boundary conditions

  7. ROOTS 必须是栈上最后一个操作，后面只能加 =

  8. 位运算符 BITAND，BITLEFT，BITNOT，BITOR，BITRIGHT， BITTEST 和 BITXOR 会将表数据中的双精度数转换为无符号的 64 位整数，然后按位运算。因此，可储存在双精度数中的最大整数为 2^53 ，更大的数都会被截断。如果被比较的数中包括 NaN ，则最终结果也为 NaN

  9. TAPER will interpret its argument to be a width in the same units as the time-axis, but if no time is provided (i.e., plain data tables) then the width is taken to be given in number of rows.

  10. 颜色转换函数，例如 RGB2HSV 等，不仅包括 rgb 到 hsv 等三元数的转换，还包括 lab 到 sRGB 等四元数的转换。这些函数在是否使用 -Q 选项下的表现不同。使用 -Q 时，输入应该为三元数，并将三个输出结果放在栈中。由于只有栈顶的元素被打印，必须使用 POP 或 ROLL 等操作符来获取其他感兴趣的项。如果没有 -Q ，GMT 将三元素视为整体，转换后，同样以整体的方式返回到堆栈中。

  11. VPDF 运算符的输入数据单位为角度

存储、调用和清空

用户可以将中间计算结果储存到一个变量中，并在后续计算中调用该变量。这在需要对某部分进行多次重复计算时可以提高效率和可读性。保存结果需要使用特殊的运算符 STO@label ， 其中 label 是变量的名称。调用该变量时，使用 [RCL]@label ， RCL 是可选的。使用后要清除该变量，可以使用 CLR@label ， STO 和 CLR 均不影响计算中的堆栈。

宏

用户可以将特定的运算符组合保存为宏文件 gmtmath.macros 。文件中可以包含任意数量的宏， # 开头的行为注释。宏的格式为

name = arg1 arg2 ... arg2[ : comment]

其中，name 是宏名，当此运算符出现在命令中时，则将其简单替换为参数列表。宏不可以互相调用。下面将给出一个宏的例子：宏 DEPTH 将在时间列以 Myr 为单位给出海底构造的年龄，并计算预测的 half-space 海深

DEPTH = SQRT 350 MUL 2500 ADD NEG : usage: DEPTH to return half-space seafloor depths

由于在宏中可能使用地理或时间常数，因此可使用 : 后加一个空格的形式作为注释的开端。下面是一个 GPSWEEK 宏，用于确定给定的时间处于的 GPS 周数（相对于某个参考时刻）:

GPSWEEK = 1980-01-06T00:00:00 SUB 86400 DIV 7 DIV FLOOR : usage: GPS week without rollover

参与运算列的选择

-Ccols 设置后，任何操作，包括从文件中加载数据都会受到影响。因此，当 -Ccols 放在数据加载前时，只有设定的列会更新，其他的列会置为 0，若不想其他列置为 0，则需把文件放在 -Ccols 前。

绝对时间列

如果输入数据有多列并且包含绝对时间列（通过 -N [0] 设置 id)，则时间输出为相对时间。在标量计算模式（-Q）中，同样会将绝对时间转换为相对时间。如果 -C 选定的列中包含绝对时间，同样会转换为相对时间。用户可以设置合适的 -f 或 -fo 选项来避免这种转换。例如：如果需要计算时间差，则相对时间更加合适；如果要计算绝对时间，则需要使用 -fo 选项设置输出时间格式为绝对时间。

对单位进行缩放

使用 -Q 计算时，请注意 GMT 内部将 c、i 或 p 等长度单位转换为 inch 计算，并对结果进行缩放到设置的单位，但这仅在线性运算能得到正确的结果。例如 gmt math -Qp 1c 0.5i ADD =，其中 1c = 28.364p ，0.5i = 36p，且 ADD 为线性运算，因此会得到正确的结果 64.346p。考虑另外一个例子：gmt math -Qc 1c 1c MUL = ，结果却不是 1，而是 0.393。原因是在非线性运算时， math 模块无法跟踪计算栈上单位，因此在内部均假定以英寸为单位，在最后将结果缩放为 cm。这里在 math 内部计算的结果是 1 英寸的平方，math 将其转换为 0.3937 厘米* 英寸 作为单位，因此结果是错误的。

示例

将两种不同单位的长度相减，并赋值给 length 变量

length=`gmt math -Q 15c 2i SUB =`

计算不同单位的长度之间的比，输出无量纲

ratio=`gmt math -Qn 15c 2i DIV =`

对 process1 产生的数据进行开方，然后使用管道传递给 process3

process1 | gmt math STDIN SQRT = | process3

对两个数据文件进行平均，并对结果取 log10

gmt math file1.txt file2.txt ADD 0.5 MUL LOG10 = file3.txt

对包含海底地形年代（以 m.y. 为单位）和深度（以 m 为单位）的数据文件 samples.txt ，使用公式 depth(单位为 m) = 2500 + 350 * sqrt (age) 计算正常深度，并最终获得深度异常

gmt math samples.txt T SQRT 350 MUL 2500 ADD SUB =

取三个文件对应列（第 1 列，第 4 到 6列）的平均值

gmt math -C1,4-6 sizes.1 sizes.2 ADD sizes.3 ADD 3 DIV = ave.txt

从含有一列数据的 ages.txt 文件中计算模值，并赋值给变量

mode_age=`gmt math -S -T ages.txt MODE =`

计算 t.txt 中坐标的 dilog(x) 函数值

gmt math -Tt.txt T DILOG = dilog.txt

下面展示一个使用存储变量的例子。将 (2*pi*T/360) 作为变量储存，对该变量分别乘以 2 和 3 ，并对上述结果求余弦，然后加和

gmt math -T0/360/1 2 PI MUL 360 DIV T MUL STO@kT COS @kT 2 MUL COS ADD @kT 3 MUL COS ADD = harmonics.txt

使用 gmtmath 实现标量计算（不含输入文件）可以使用 -Q 选项，计算 kei(((1 + 1.75)/2.2) + cos (60)) 并将结果赋值给 shell 中的变量 z

z=`gmt math -Q 1 1.75 ADD 2.2 DIV 60 COSD ADD KEI =`

将黄色的 RGB 值转换为 HSV 并保存为 hue

set hue=`gmt math -Q 255 255 0 RGB2HSV POP POP =`

下面将展示使用 gmtmath 求解方程组。假设当前存在一个增广矩阵 [ A ｜b ]，A 为系数矩阵，b 为等号右边的常数列，求解的方程组为 A * x = b 。LSQFIT 操作符可实现方程的求解，但前提是使用 -A 选项正确放置参数。假定表数据文件 ty.txt 包含 t 和 y(t) 两列，拟合模型为 y(t) = a + b*t + c*H(t-t0) ，其中 H 为给定 t0 = 1.55 的 Heaviside 阶跃函数。这时，需要四列的增广矩阵，第二列为 t ，第四列为 y(t)（或者说 b） ，则计算命令为

gmt math -N4/1 -Aty.txt -C0 1 ADD -C2 1.55 STEPT ADD -Ca LSQFIT = solution.txt

上面的例子中使用了 -C 选项来激活操作的列，并在 LSQFIT 求解前激活了所有列（ -Ca ）; 第二列和第四列（列编号为 1 和 3）被设置为 t 和 y(t)。如果用户自己已经准备了增广矩阵文件 lsqsys.txt ， 求解将更加简单

gmt math -T lsqsys.txt LSQFIT = solution.txt

下面展示 -C 选项的位置对于数据加载的影响

echo 1 2 3 4 | gmt math STDIN -C3 1 ADD =
1    2    3    5

echo 1 2 3 4 | gmt math -C3 STDIN 1 ADD =
0    0    0    5

参考文献

Abramowitz, M., and I. A. Stegun, 1964, Handbook of Mathematical Functions, Applied Mathematics Series, vol. 55, Dover, New York.

Holmes, S. A., and W. E. Featherstone, 2002, A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalized associated Legendre functions. Journal of Geodesy, 76, 279-299.

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, 1992, Numerical Recipes, 2nd edition, Cambridge Univ., New York.

Spanier, J., and K. B. Oldman, 1987, An Atlas of Functions, Hemisphere Publishing Corp.

相关模块

grdmath