1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Misc statistical and special functions.
19 *
20 * Author: Walter H. F. Smith, P. Wessel, R. Scharroo
21 * Date: 1-JAN-2010
22 * Version: 5.x
23 *
24 * PUBLIC functions (59):
25 *
26 * gmt_sig_f : : Returns true if reduction in model misfit was significant
27 * gmt_bei: : Kelvin-Bessel function bei(x)
28 * gmt_ber: : Kelvin-Bessel function ber(x)
29 * gmt_kei: : Kelvin-Bessel function kei(x)
30 * gmt_ker: : Kelvin-Bessel function ker(x)
31 * gmt_plm: : Legendre polynomial of degree L order M
32 * gmt_plm_bar: : Normalized Legendre polynomial of degree L order M
33 * gmt_plm_bar_all :
34 * gmt_i0: : Modified Bessel function 1st kind order 0
35 * gmt_i1: : Modified Bessel function 1st kind order 1
36 * gmt_in: : Modified Bessel function 1st kind order N
37 * gmt_k0: : Modified Kelvin function 2nd kind order 0
38 * gmt_k1: : Modified Kelvin function 2nd kind order 1
39 * gmt_kn: : Modified Kelvin function 2nd kind order N
40 * gmt_dilog: : The dilog function
41 * gmt_erfinv: : The inverse error function
42 * gmt_rand: : Uniformly distributed random numbers 0 < x < 1
43 * gmt_nrand: : Normally distributed random numbers from N(0,1)
44 * gmt_lrand: : Laplace random number generator
45 * gmt_corrcoeff: : Correlation coefficient.
46 * gmt_psi: : Digamma (psi) function.
47 * gmt_PvQv: : Legendre functions Pv and Qv for imaginary v and real x (-1/+1).
48 * gmt_factorial: : Factorials.
49 * gmt_sinc :
50 * gmt_permutation :
51 * gmt_combination :
52 * gmt_f_pdf :
53 * gmt_f_cdf :
54 * gmt_t_pdf :
55 * gmt_t_cdf :
56 * gmt_weibull_pdf :
57 * gmt_weibull_cdf :
58 * gmt_weibull_crit :
59 * gmt_binom_pdf :
60 * gmt_binom_cdf :
61 * gmt_vonmises_pdf :
62 * gmt_zdist :
63 * gmt_zcrit :
64 * gmt_tcrit :
65 * gmt_chi2_pdf :
66 * gmt_chi2crit :
67 * gmt_Fcrit :
68 * gmt_chi2 :
69 * gmt_poissonpdf :
70 * gmt_poisson_cdf :
71 * gmt_mean_and_std :
72 * gmt_median :
73 * gmt_mean_weighted :
74 * gmt_quantile_weighted :
75 * gmt_median_weighted :
76 * gmt_mode_weighted :
77 * gmt_mode :
78 * gmt_mode_f :
79 * gmt_getmad :
80 * gmt_getmad_f :
81 * gmt_extreme :
82 * gmt_chebyshev :
83 * gmt_corrcoeff :
84 * gmt_corrcoeff_f :
85 * gmt_quantile :
86 * gmt_quantile_f :
87 */
88
89 #include "gmt_dev.h"
90 #include "gmt_internals.h"
91
92 enum gmtstat_enum_cplx {GMT_RE = 0, GMT_IM = 1}; /* Real and imaginary indices */
93
94 /* --------- Local functions to this file ------- */
95
96 #define ITMAX 100
97
gmtstat_ln_gamma(struct GMT_CTRL * GMT,double xx)98 GMT_LOCAL double gmtstat_ln_gamma (struct GMT_CTRL *GMT, double xx) {
99 /* Routine to compute natural log of Gamma(x)
100 by Lanczos approximation. Most accurate
101 for x > 1; fails for x <= 0. No error
102 checking is done here; it is assumed
103 that this is called by gmtstat_ln_gamma_r() */
104
105 static double cof[6] = {
106 76.18009173,
107 -86.50532033,
108 24.01409822,
109 -1.231739516,
110 0.120858003e-2,
111 -0.536382e-5
112 };
113
114 double x, tmp, ser;
115
116 int i;
117
118 x = xx - 1.0;
119 tmp = x + 5.5;
120 tmp = (x + 0.5) * d_log (GMT,tmp) - tmp;
121 ser = 1.0;
122 for (i = 0; i < 6; i++) {
123 x += 1.0;
124 ser += (cof[i]/x);
125 }
126 return (tmp + d_log (GMT, 2.50662827465*ser) );
127 }
128
gmtstat_ln_gamma_r(struct GMT_CTRL * GMT,double x,double * lngam)129 GMT_LOCAL int gmtstat_ln_gamma_r (struct GMT_CTRL *GMT, double x, double *lngam) {
130 /* Get natural logarithm of Gamma(x), x > 0.
131 To maintain full accuracy, this
132 routine uses Gamma(1 + x) / x when
133 x < 1. This routine in turn calls
134 gmtstat_ln_gamma(x), which computes the
135 actual function value. gmtstat_ln_gamma
136 assumes it is being called in a
137 smart way, and does not check the
138 range of x. */
139
140 if (x > 1.0) {
141 *lngam = gmtstat_ln_gamma (GMT, x);
142 return (0);
143 }
144 if (x > 0.0 && x < 1.0) {
145 *lngam = gmtstat_ln_gamma (GMT, 1.0 + x) - d_log (GMT,x);
146 return (0);
147 }
148 if (x == 1.0) {
149 *lngam = 0.0;
150 return (0);
151 }
152 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Ln Gamma: Bad x (x <= 0).\n");
153 return (GMT_NOTSET);
154 }
155
gmtstat_gamma_ser(struct GMT_CTRL * GMT,double * gamser,double a,double x,double * gln)156 GMT_LOCAL void gmtstat_gamma_ser (struct GMT_CTRL *GMT, double *gamser, double a, double x, double *gln) {
157 /* Returns the incomplete gamma function P(a,x) by series rep.
158 * Press et al, gser() */
159
160 int n;
161 double sum, del, ap;
162
163 if (gmtstat_ln_gamma_r (GMT, a, gln) == GMT_NOTSET) {
164 *gamser = GMT->session.d_NaN;
165 return;
166 }
167
168 if (x < 0.0) {
169 GMT_Report (GMT->parent, GMT_MSG_WARNING, "x < 0 in gmtstat_gamma_ser(x)\n");
170 *gamser = GMT->session.d_NaN;
171 return;
172 }
173 if (x == 0.0) {
174 *gamser = 0.0;
175 return;
176 }
177 ap = a;
178 del = sum = 1.0 / a;
179 for (n = 1; n <= ITMAX; n++) {
180 ap += 1.0;
181 del *= x / ap;
182 sum += del;
183 if (fabs (del) < fabs (sum) * DBL_EPSILON) {
184 *gamser = sum * exp (-x + a * log (x) - (*gln));
185 return;
186 }
187 }
188 GMT_Report (GMT->parent, GMT_MSG_WARNING, "a too large, ITMAX too small in gmtstat_gamma_ser(x)\n");
189 }
190
gmtstat_gamma_cf(struct GMT_CTRL * GMT,double * gammcf,double a,double x,double * gln)191 GMT_LOCAL void gmtstat_gamma_cf (struct GMT_CTRL *GMT, double *gammcf, double a, double x, double *gln) {
192 /* Returns the incomplete gamma function P(a,x) by continued fraction.
193 * Press et al, gcf() */
194 int n;
195 double gold = 0.0, g, fac = 1.0, b1 = 1.0;
196 double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
197
198 if (gmtstat_ln_gamma_r (GMT, a, gln) == GMT_NOTSET) {
199 *gln = GMT->session.d_NaN;
200 return;
201 }
202
203 a1 = x;
204 for (n = 1; n <= ITMAX; n++) {
205 an = (double) n;
206 ana = an - a;
207 a0 = (a1 + a0 * ana) * fac;
208 b0 = (b1 + b0 * ana) * fac;
209 anf = an * fac;
210 a1 = x * a0 + anf * a1;
211 b1 = x * b0 + anf * b1;
212 if (a1 != 0.0) {
213 fac = 1.0 / a1;
214 g = b1 * fac;
215 if (fabs ((g - gold) / g) < DBL_EPSILON) {
216 *gammcf = exp (-x + a * log (x) - (*gln)) * g;
217 return;
218 }
219 gold = g;
220 }
221 }
222 GMT_Report (GMT->parent, GMT_MSG_WARNING, "a too large, ITMAX too small in gmtstat_gamma_cf(x)\n");
223 }
224
gmtstat_gammq(struct GMT_CTRL * GMT,double a,double x)225 GMT_LOCAL double gmtstat_gammq (struct GMT_CTRL *GMT, double a, double x) {
226 /* Returns Q(a,x) = 1 - P(a,x) Inc. Gamma function */
227
228 double G = 0.0, gln;
229
230 if (x < 0.0 || a <= 0.0) {
231 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Invalid arguments to GMT_gammaq\n");
232 return (GMT->session.d_NaN);
233 }
234
235 if (x < (a + 1.0)) {
236 gmtstat_gamma_ser (GMT, &G, a, x, &gln);
237 return (1.0 - G);
238 }
239 gmtstat_gamma_cf (GMT, &G, a, x, &gln);
240 return (G);
241 }
242
243 #define BETA_EPS 3.0e-7
244
gmtstat_cf_beta(struct GMT_CTRL * GMT,double a,double b,double x)245 GMT_LOCAL double gmtstat_cf_beta (struct GMT_CTRL *GMT, double a, double b, double x) {
246 /* Continued fraction method called by gmtstat_inc_beta. */
247 double am = 1.0, bm = 1.0, az = 1.0;
248 double qab, qap, qam, bz, em, tem, d;
249 double ap, bp, app, bpp, aold;
250
251 int m = 0;
252
253 qab = a + b;
254 qap = a + 1.0;
255 qam = a - 1.0;
256 bz = 1.0 - qab * x / qap;
257
258 do {
259 m++;
260 em = (double)m;
261 tem = em + em;
262 d = em*(b-m)*x/((qam+tem)*(a+tem));
263 ap = az+d*am;
264 bp = bz+d*bm;
265 d = -(a+m)*(qab+em)*x/((a+tem)*(qap+tem));
266 app = ap+d*az;
267 bpp = bp+d*bz;
268 aold = az;
269 am = ap/bpp;
270 bm = bp/bpp;
271 az = app/bpp;
272 bz = 1.0;
273 } while (((fabs (az-aold) ) >= (BETA_EPS * fabs (az))) && (m < ITMAX));
274
275 if (m == ITMAX) GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_cf_beta: A or B too big, or ITMAX too small.\n");
276
277 return (az);
278 }
279
gmtstat_inc_beta(struct GMT_CTRL * GMT,double a,double b,double x,double * ibeta)280 GMT_LOCAL int gmtstat_inc_beta (struct GMT_CTRL *GMT, double a, double b, double x, double *ibeta) {
281 double bt, gama, gamb, gamab;
282
283 if (a <= 0.0) {
284 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_inc_beta: Bad a (a <= 0).\n");
285 return(GMT_NOTSET);
286 }
287 if (b <= 0.0) {
288 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_inc_beta: Bad b (b <= 0).\n");
289 return(GMT_NOTSET);
290 }
291 if (x > 0.0 && x < 1.0) {
292 gmtstat_ln_gamma_r(GMT, a, &gama);
293 gmtstat_ln_gamma_r(GMT, b, &gamb);
294 gmtstat_ln_gamma_r(GMT, (a+b), &gamab);
295 bt = exp(gamab - gama - gamb
296 + a * d_log (GMT, x) + b * d_log (GMT, 1.0 - x) );
297
298 /* Here there is disagreement on the range of x which
299 converges efficiently. Abramowitz and Stegun
300 say to use x < (a - 1) / (a + b - 2). Editions
301 of Numerical Recipes through mid 1987 say
302 x < ( (a + 1) / (a + b + 1), but the code has
303 x < ( (a + 1) / (a + b + 2). Editions printed
304 late 1987 and after say x < ( (a + 1) / (a + b + 2)
305 in text as well as code. What to do ? */
306
307 if (x < ( (a + 1) / (a + b + 2) ) )
308 *ibeta = bt * gmtstat_cf_beta (GMT, a, b, x) / a;
309 else
310 *ibeta = 1.0 - bt * gmtstat_cf_beta (GMT, b, a, (1.0 - x) ) / b;
311 return(0);
312 }
313 else if (x == 0.0) {
314 *ibeta = 0.0;
315 return (0);
316 }
317 else if (x == 1.0) {
318 *ibeta = 1.0;
319 return (0);
320 }
321 else if (x < 0.0) {
322 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_inc_beta: Bad x (x < 0).\n");
323 *ibeta = 0.0;
324 }
325 else if (x > 1.0) {
326 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_inc_beta: Bad x (x > 1).\n");
327 *ibeta = 1.0;
328 }
329 return (GMT_NOTSET);
330 }
331
gmtstat_f_q(struct GMT_CTRL * GMT,double chisq1,uint64_t nu1,double chisq2,uint64_t nu2,double * prob)332 GMT_LOCAL int gmtstat_f_q (struct GMT_CTRL *GMT, double chisq1, uint64_t nu1, double chisq2, uint64_t nu2, double *prob) {
333 /* Routine to compute Q(F, nu1, nu2) = 1 - P(F, nu1, nu2), where nu1
334 and nu2 are positive integers, chisq1 and chisq2 are random
335 variables having chi-square distributions with nu1 and nu2
336 degrees of freedom, respectively (chisq1 and chisq2 >= 0.0),
337 F = (chisq1/nu1)/(chisq2/nu2) has the F-distribution, and
338 P(F, nu1, nu2) is the cumulative F-distribution, that is,
339 the integral from 0 to (chisq1/nu1)/(chisq2/nu2) of the F-
340 distribution. Q = 1 - P is small when (chisq1/nu1)/(chisq2/nu2)
341 is large with respect to 1. That is, the value returned by
342 this routine is the likelihood that an F >= (chisq1/nu1)/
343 (chisq2/nu2) would occur by chance.
344
345 Follows Abramowitz and Stegun.
346 This is different from the method in Numerical Recipes, which
347 uses the incomplete beta function but makes no use of the fact
348 that nu1 and nu2 are known to be integers, and thus there is
349 a finite limit on the sum for their expression.
350
351 W H F Smith, August, 1999.
352
353 REVISED by W H F Smith, October 27, 2000 after GMT 3.3.6 release.
354 I found that the A&S methods overflowed for large nu1 and nu2, so
355 I decided to go back to the gmtstat_inc_beta way of doing things.
356
357 */
358
359 /* Check range of arguments: */
360
361 if (nu1 <= 0 || nu2 <= 0 || chisq1 < 0.0 || chisq2 < 0.0) {
362 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_f_q: Bad argument(s).\n");
363 return (GMT_NOTSET);
364 }
365
366 /* Extreme cases evaluate immediately: */
367
368 if (chisq1 == 0.0) {
369 *prob = 1.0;
370 return (0);
371 }
372 if (chisq2 == 0.0) {
373 *prob = 0.0;
374 return (0);
375 }
376
377 /* REVISION of Oct 27, 2000: This inc beta call here returns
378 the value. All subsequent code is not used. */
379
380 if (gmtstat_inc_beta (GMT, 0.5*nu2, 0.5*nu1, chisq2/(chisq2+chisq1), prob) ) {
381 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_f_q: Trouble in gmtstat_inc_beta call.\n");
382 return (GMT_NOTSET);
383 }
384 return (0);
385 }
386
gmtstat_f_test_new(struct GMT_CTRL * GMT,double chisq1,uint64_t nu1,double chisq2,uint64_t nu2,double * prob,int iside)387 GMT_LOCAL int gmtstat_f_test_new (struct GMT_CTRL *GMT, double chisq1, uint64_t nu1, double chisq2, uint64_t nu2, double *prob, int iside) {
388 /* Given chisq1 and chisq2, random variables distributed as chi-square
389 with nu1 and nu2 degrees of freedom, respectively, except that
390 chisq1 is scaled by var1, and chisq2 is scaled by var2, let
391 the null hypothesis, H0, be that var1 = var2. This routine
392 assigns prob, the probability that we can reject H0 in favor
393 of a new hypothesis, H1, according to iside:
394 iside=+1 means H1 is that var1 > var2
395 iside=-1 means H1 is that var1 < var2
396 iside=0 means H1 is that var1 != var2.
397 This routine differs from the old gmtstat_f_test() by adding the
398 argument iside and allowing one to choose the test. The old
399 routine in effect always set iside=0.
400 This routine also differs from gmtstat_f_test() in that the former
401 used the incomplete beta function and this one uses gmtstat_f_q().
402
403 Returns 0 on success, -1 on failure.
404
405 WHF Smith, 12 August 1999.
406 */
407
408 double q; /* The probability from gmtstat_f_q(), which is the prob
409 that H0 should be retained even though
410 chisq1/nu1 > chisq2/nu2. */
411
412 if (chisq1 <= 0.0 || chisq2 <= 0.0 || nu1 < 1 || nu2 < 1) {
413 *prob = GMT->session.d_NaN;
414 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_f_test_new: Bad argument(s).\n");
415 return (GMT_NOTSET);
416 }
417
418 gmtstat_f_q (GMT, chisq1, nu1, chisq2, nu2, &q);
419
420 if (iside > 0)
421 *prob = 1.0 - q;
422 else if (iside < 0)
423 *prob = q;
424 else if ((chisq1/nu1) <= (chisq2/nu2))
425 *prob = 2.0*q;
426 else
427 *prob = 2.0*(1.0 - q);
428
429 return (0);
430 }
431
gmtstat_factln(struct GMT_CTRL * GMT,int n)432 GMT_LOCAL double gmtstat_factln (struct GMT_CTRL *GMT, int n) {
433 /* returns log(n!) */
434 static double a[101]; /* Automatically initialized to zero */
435 if (n < 0) {
436 GMT_Report (GMT->parent, GMT_MSG_WARNING, "n < 0 in gmtstat_factln(n)\n");
437 return (GMT->session.d_NaN);
438 }
439 if (n <= 1) return 0.0;
440 if (n <= 100) return (a[n] ? a[n] : (a[n] = gmtstat_ln_gamma (GMT, n+1.0)));
441 else return gmtstat_ln_gamma (GMT, n+1.0);
442 }
443
gmtstat_beta(struct GMT_CTRL * GMT,double z,double w)444 GMT_LOCAL double gmtstat_beta (struct GMT_CTRL *GMT, double z, double w) {
445 double g1 = 0.0, g2 = 0.0, g3 = 0.0;
446 gmtstat_ln_gamma_r (GMT, z, &g1);
447 gmtstat_ln_gamma_r (GMT, w, &g2);
448 gmtstat_ln_gamma_r (GMT, z+w, &g3);
449 return exp (g1 + g2 - g3);
450 }
451
452 #if 0 /* Legacy code */
453 GMT_LOCAL int gmtstat_f_test (struct GMT_CTRL *GMT, double chisq1, uint64_t nu1, double chisq2, uint64_t nu2, double *prob) {
454 /* Routine to compute the probability that
455 two variances are the same.
456 chisq1 is distributed as chisq with
457 nu1 degrees of freedom; ditto for
458 chisq2 and nu2. If these are independent
459 and we form the ratio
460 F = max(chisq1,chisq2)/min(chisq1,chisq2)
461 then we can ask what is the probability
462 that an F greater than this would occur
463 by chance. It is this probability that
464 is returned in prob. When prob is small,
465 it is likely that the two chisq represent
466 two different populations; the confidence
467 that the two do not represent the same pop
468 is 1.0 - prob. This is a two-sided test.
469 This follows some ideas in Numerical Recipes, CRC Handbook,
470 and Abramowitz and Stegun. */
471
472 double f, df1, df2, p1, p2;
473
474 if (chisq1 <= 0.0) {
475 GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmtstat_f_test: Chi-Square One <= 0.0\n");
476 return(GMT_NOTSET);
477 }
478 if (chisq2 <= 0.0) {
479 GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmtstat_f_test: Chi-Square Two <= 0.0\n");
480 return(GMT_NOTSET);
481 }
482 if (chisq1 > chisq2) {
483 f = chisq1/chisq2;
484 df1 = (double)nu1;
485 df2 = (double)nu2;
486 }
487 else {
488 f = chisq2/chisq1;
489 df1 = (double)nu2;
490 df2 = (double)nu1;
491 }
492 if (gmtstat_inc_beta(GMT, 0.5*df2, 0.5*df1, df2/(df2+df1*f), &p1) ) {
493 GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmtstat_f_test: Trouble on 1st gmtstat_inc_beta call.\n");
494 return(GMT_NOTSET);
495 }
496 if (gmtstat_inc_beta(GMT, 0.5*df1, 0.5*df2, df1/(df1+df2/f), &p2) ) {
497 GMT_Report (GMT->parent, GMT_MSG_ERROR, "gmtstat_f_test: Trouble on 2nd gmtstat_inc_beta call.\n");
498 return(GMT_NOTSET);
499 }
500 *prob = p1 + (1.0 - p2);
501 return (0);
502 }
503 #endif
504
gmtstat_student_t_a(struct GMT_CTRL * GMT,double t,uint64_t n,double * prob)505 GMT_LOCAL int gmtstat_student_t_a (struct GMT_CTRL *GMT, double t, uint64_t n, double *prob) {
506 /* Probability integral called A(t,n) by Abramowitz &
507 Stegun for the student's t distribution with n degrees
508 of freedom. Uses expressions A&S 26.7.3 and 26.7.4
509
510 If X is distributed N(0,1) and V is distributed chi-
511 square with n degrees of freedom, then
512 tau = X / sqrt(V/n) is said to have Student's t-
513 distribution with n degrees of freedom. For example,
514 tau could be the sample mean divided by the sample
515 standard deviation, for a sample of N points; then
516 n = N - 1.
517
518 This function sets *prob = GMT->session.d_NaN and returns (-1)
519 if t < 0. Otherwise it sets *prob = the probability
520 fabs(tau) <= t and returns (0).
521
522 As n -> oo, we can replace this function with
523 erf (t / M_SQRT2). However, it isn't clear how large
524 n has to be to make this a good approximation. I
525 consulted six books; one of them suggested this
526 approximation for n >= 30, but all the others did not
527 say when to use this approximation (A&S, in particular,
528 does not say). I tried some numerical experiments
529 which suggested that the relative error in this
530 approximation would be < 0.01 for n > 30, all t, but
531 I also found that the expression here is stable to
532 large n and large t, so I decided to leave it as is.
533
534 W H F Smith, August 1999.
535 */
536
537 double theta, s, c, csq, term, sum;
538 int64_t k, kstop, kt, kb;
539
540 if (t < 0.0 || n == 0) {
541 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmtstat_student_t_a: Bad argument(s).\n");
542 *prob = GMT->session.d_NaN;
543 return (GMT_NOTSET);
544 }
545
546 if (t == 0.0) {
547 *prob = 0.0;
548 return (0);
549 }
550
551 theta = atan (t/sqrt ((double)n));
552
553 if (n == 1) {
554 *prob = 2.0 * theta / M_PI;
555 return (0);
556 }
557
558 sincos (theta, &s, &c);
559
560 csq = c * c;
561
562 kstop = n-2;
563 if (n%2 == 1) {
564 kt = 0;
565 kb = 1;
566 k = 1;
567 term = c;
568 }
569 else {
570 kt = -1;
571 kb = 0;
572 k = 0;
573 term = 1.0;
574 }
575 sum = term;
576 while (k < kstop) {
577 k += 2;
578 kt += 2;
579 kb += 2;
580 term *= (kt * csq)/kb;
581 sum += term;
582 }
583
584 sum *= s;
585
586 if (n%2 == 1)
587 *prob = 2.0 * (theta + sum) / M_PI;
588 else
589 *prob = sum;
590
591 /* Adjust in case of roundoff: */
592
593 if (*prob < 0.0) *prob = 0.0;
594 if (*prob > 1.0) *prob = 1.0;
595
596 return (0);
597 }
598
gmtstat_Cmul(double A[],double B[],double C[])599 GMT_LOCAL void gmtstat_Cmul (double A[], double B[], double C[]) {
600 /* Complex multiplication */
601 C[GMT_RE] = A[GMT_RE]*B[GMT_RE] - A[GMT_IM]*B[GMT_IM];
602 C[GMT_IM] = A[GMT_RE]*B[GMT_IM] + A[GMT_IM]*B[GMT_RE];
603 }
604
gmtstat_Cdiv(double A[],double B[],double C[])605 GMT_LOCAL void gmtstat_Cdiv (double A[], double B[], double C[]) {
606 /* Complex division */
607 double denom;
608 denom = B[GMT_RE]*B[GMT_RE] + B[GMT_IM]*B[GMT_IM];
609 C[GMT_RE] = (A[GMT_RE]*B[GMT_RE] + A[GMT_IM]*B[GMT_IM])/denom;
610 C[GMT_IM] = (A[GMT_IM]*B[GMT_RE] - A[GMT_RE]*B[GMT_IM])/denom;
611 }
612
613 #if 0 /* Unused */
614 GMT_LOCAL void gmtstat_Ccot (double Z[], double cotZ[]) {
615 /* Complex cot(z) */
616 double sx, cx, e, A[2], B[2];
617
618 sincos (2.0*Z[0], &sx, &cx);
619 e = exp (-2.0*Z[1]);
620 A[0] = -e * sx; A[1] = B[0] = e * cx;
621 A[1] += 1.0; B[0] -= 1.0; B[1] = -A[0];
622 gmtstat_Cdiv (A, B, cotZ);
623 }
624 #endif
625
gmtstat_Cabs(double A[])626 GMT_LOCAL double gmtstat_Cabs (double A[]) {
627 return (hypot (A[GMT_RE], A[GMT_IM]));
628 }
629
gmtstat_F_to_ch1_ch2(struct GMT_CTRL * GMT,double F,double nu1,double nu2,double * chisq1,double * chisq2)630 GMT_LOCAL void gmtstat_F_to_ch1_ch2 (struct GMT_CTRL *GMT, double F, double nu1, double nu2, double *chisq1, double *chisq2) {
631 /* Silly routine to break F up into parts needed for gmtstat_f_q */
632 gmt_M_unused(GMT);
633 *chisq2 = 1.0;
634 *chisq1 = F * nu1 / nu2;
635 }
636
gmtlib_compare_observation(const void * a,const void * b)637 int gmtlib_compare_observation (const void *a, const void *b) {
638 const struct GMT_OBSERVATION *obs_1 = a, *obs_2 = b;
639
640 /* Sorts observations into ascending order based on obs->value */
641 if (obs_1->value < obs_2->value)
642 return -1;
643 if (obs_1->value > obs_2->value)
644 return 1;
645 return 0;
646 }
647
gmt_sig_f(struct GMT_CTRL * GMT,double chi1,uint64_t n1,double chi2,uint64_t n2,double level,double * prob)648 int gmt_sig_f (struct GMT_CTRL *GMT, double chi1, uint64_t n1, double chi2, uint64_t n2, double level, double *prob) {
649 /* Returns true if chi1/n1 significantly less than chi2/n2
650 at the level level. Returns false if:
651 error occurs in gmtstat_f_test_new();
652 chi1/n1 not significantly < chi2/n2 at level.
653
654 Changed 12 August 1999 to use gmtstat_f_test_new() */
655
656 int trouble;
657
658 trouble = gmtstat_f_test_new (GMT, chi1, n1, chi2, n2, prob, -1);
659 if (trouble) return (0);
660 return ((*prob) >= level);
661 }
662
663 /*
664 * Kelvin-Bessel functions, ber, bei, ker, kei. ber(x) and bei(x) are even;
665 * ker(x) and kei(x) are defined only for x > 0 and x >=0, respectively.
666 * For x <= 8 we use polynomial approximations in Abramowitz & Stegun.
667 * For x > 8 we use asymptotic series of Russell, quoted in Watson (Theory
668 * of Bessel Functions).
669 */
670
gmt_ber(struct GMT_CTRL * GMT,double x)671 double gmt_ber (struct GMT_CTRL *GMT, double x) {
672 double t, rxsq, alpha, beta;
673 gmt_M_unused(GMT);
674
675 if (x == 0.0) return (1.0);
676
677 /* ber is an even function of x: */
678 x = fabs (x);
679
680 if (x <= 8.0) {
681 /* Telescoped power series from Abramowitz & Stegun */
682 t = x * 0.125;
683 t *= t;
684 t *= t; /* t = pow(x/8, 4) */
685 return (1.0 + t*(-64.0 + t*(113.77777774 + t*(-32.36345652 + t*(2.64191397 + t*(-0.08349609 + t*(0.00122552 - 0.00000901 * t)))))));
686 }
687 else {
688 /* Russell's asymptotic approximation, from Watson, p. 204 */
689
690 rxsq = 1.0 / (x * x);
691 t = x / M_SQRT2;
692
693 alpha = t;
694 beta = t - 0.125 * M_PI;
695 t *= 0.125 * rxsq;
696 alpha += t;
697 beta -= t;
698 beta -= 0.0625*rxsq;
699 t *= (25.0/48.0)*rxsq;
700 alpha -= t;
701 beta -= t;
702 alpha -= (13.0/128.0)*(rxsq*rxsq);
703
704 return (exp (alpha) * cos (beta) / sqrt (2.0 * M_PI * x) );
705 }
706 }
707
gmt_bei(struct GMT_CTRL * GMT,double x)708 double gmt_bei (struct GMT_CTRL *GMT, double x) {
709 double t, rxsq, alpha, beta;
710 gmt_M_unused(GMT);
711
712 if (x == 0.0) return (0.0);
713
714 /* bei is an even function of x: */
715 x = fabs(x);
716
717 if (x <= 8.0) {
718 /* Telescoped power series from Abramowitz & Stegun */
719 t = x * 0.125;
720 rxsq = t*t;
721 t = rxsq * rxsq; /* t = pow(x/8, 4) */
722 return (rxsq * (16.0 + t*(-113.77777774 + t*(72.81777742 + t*(-10.56765779 + t*(0.52185615 + t*(-0.01103667 + t*(0.00011346))))))));
723 }
724 else {
725 /* Russell's asymptotic approximation, from Watson, p. 204 */
726
727 rxsq = 1.0 / (x * x);
728 t = x / M_SQRT2;
729
730 alpha = t;
731 beta = t - 0.125 * M_PI;
732 t *= 0.125 * rxsq;
733 alpha += t;
734 beta -= t;
735 beta -= 0.0625*rxsq;
736 t *= (25.0/48.0)*rxsq;
737 alpha -= t;
738 beta -= t;
739 alpha -= (13.0/128.0)*(rxsq*rxsq);
740
741 return (exp (alpha) * sin (beta) / sqrt (2.0 * M_PI * x));
742 }
743 }
744
gmt_ker(struct GMT_CTRL * GMT,double x)745 double gmt_ker (struct GMT_CTRL *GMT, double x) {
746 double t, rxsq, alpha, beta;
747
748 if (x <= 0.0) {
749 GMT_Report (GMT->parent, GMT_MSG_WARNING, "x <= 0 in gmt_ker(x)\n");
750 return (GMT->session.d_NaN);
751 }
752
753 if (x <= 8.0) {
754 /* Telescoped power series from Abramowitz & Stegun */
755 t = 0.125 * x;
756 t *= t;
757 t *= t; /* t = pow(x/8, 4) */
758 return (-log (0.5 * x) * gmt_ber (GMT, x) + 0.25 * M_PI * gmt_bei (GMT, x) -M_EULER + \
759 t * (-59.05819744 + t * (171.36272133 + t * (-60.60977451 + t * (5.65539121 + t * (-0.199636347 + t * (0.00309699 + t * (-0.00002458 * t))))))));
760 }
761 else {
762 /* Russell's asymptotic approximation, from Watson, p. 204 */
763
764 rxsq = 1.0 / (x * x);
765 t = -x / M_SQRT2;
766
767 alpha = t;
768 beta = t - 0.125 * M_PI;
769 t *= 0.125 * rxsq;
770 alpha += t;
771 beta -= t;
772 beta -= 0.0625*rxsq;
773 t *= (25.0/48.0)*rxsq;
774 alpha -= t;
775 beta -= t;
776 alpha -= (13.0/128.0)*(rxsq*rxsq);
777
778 return (exp (alpha) * cos (beta) / sqrt (2.0 * x / M_PI));
779 }
780 }
781
gmt_kei(struct GMT_CTRL * GMT,double x)782 double gmt_kei (struct GMT_CTRL *GMT, double x) {
783 double t, rxsq, alpha, beta;
784
785 if (x <= 0.0) {
786 /* Zero is valid. If near enough to zero, return kei(0) */
787 if (x > -GMT_CONV8_LIMIT) return (-0.25 * M_PI);
788
789 GMT_Report (GMT->parent, GMT_MSG_WARNING, "x < 0 in gmt_kei(x)\n");
790 return (GMT->session.d_NaN);
791 }
792
793 if (x <= 8.0) {
794 /* Telescoped power series from Abramowitz & Stegun */
795 t = x * 0.125;
796 rxsq = t*t;
797 t = rxsq * rxsq; /* t = pow(x/8, 4) */
798 return (-log (0.5 * x) * gmt_bei (GMT, x) - 0.25 * M_PI * gmt_ber (GMT, x) +
799 rxsq * (6.76454936 + t * (-142.91827687 + t * (124.23569650 + t * (-21.30060904 + t * (1.17509064 + t * (-0.02695875 + t * (0.00029532 * t))))))));
800 }
801 else {
802 /* Russell's asymptotic approximation, from Watson, p. 204 */
803
804 rxsq = 1.0 / (x * x);
805 t = -x / M_SQRT2;
806
807 alpha = t;
808 beta = t - 0.125 * M_PI;
809 t *= 0.125 * rxsq;
810 alpha += t;
811 beta -= t;
812 beta -= 0.0625*rxsq;
813 t *= (25.0/48.0)*rxsq;
814 alpha -= t;
815 beta -= t;
816 alpha -= (13.0/128.0)*(rxsq*rxsq);
817
818 return (exp (alpha) * sin (beta) / sqrt (2.0 * x / M_PI));
819 }
820 }
821
gmt_i0(struct GMT_CTRL * GMT,double x)822 double gmt_i0 (struct GMT_CTRL *GMT, double x) {
823 /* Modified from code in Press et al. */
824 double y, res;
825 gmt_M_unused(GMT);
826
827 if (x < 0.0) x = -x;
828
829 if (x < 3.75) {
830 y = x * x / 14.0625;
831 res = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
832 }
833 else {
834 y = 3.75 / x;
835 res = (exp (x) / sqrt (x)) * (0.39894228 + y * (0.1328592e-1 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2 + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
836 }
837 return (res);
838 }
839
gmt_i1(struct GMT_CTRL * GMT,double x)840 double gmt_i1 (struct GMT_CTRL *GMT, double x) {
841 /* Modified Bessel function I1(x) */
842 double y, res;
843 gmt_M_unused(GMT);
844
845 if (x < 0.0) x = -x;
846 if (x < 3.75) {
847 y = pow (x / 3.75, 2.0);
848 res = x * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934 + y * (0.02658733 + y * (0.00301532 + y * 0.00032411))))));
849 }
850 else {
851 y = 3.75 / x;
852 res = (exp (x) / sqrt (x)) * (0.39894228 + y * (-0.03988024 + y * (-0.00362018 + y * (0.00163801+ y * (-0.01031555 + y * (0.02282967 + y * (-0.02895312 + y * (0.01787654 -y * 0.00420059))))))));
853 if (x < 0.0) res = - res;
854 }
855 return (res);
856 }
857
gmt_in(struct GMT_CTRL * GMT,unsigned int n,double x)858 double gmt_in (struct GMT_CTRL *GMT, unsigned int n, double x) {
859 /* Modified Bessel function In(x) */
860
861 unsigned int j, m, IACC = 40;
862 double res, tox, bip, bi, bim;
863 double BIGNO = 1.0e10, BIGNI = 1.0e-10;
864
865 if (n == 0) return (gmt_i0 (GMT, x));
866 if (n == 1) return (gmt_i1 (GMT, x));
867 if (x == 0.0) return (0.0);
868
869 tox = 2.0 / fabs (x);
870 bip = res = 0.0;
871 bi = 1.0;
872 m = 2 * (n + urint (sqrt ((double)(IACC * n))));
873 for (j = m; j >= 1; j--) {
874 bim = bip + ((double)j) * tox * bi;
875 bip = bi;
876 bi = bim;
877 if (fabs (bi) > BIGNO) {
878 res *= BIGNI;
879 bi *= BIGNI;
880 bip *= BIGNI;
881 }
882 if (j == n) res = bip;
883 }
884 res *= (gmt_i0 (GMT, x) / bi);
885 if (x < 0.0 && (n%2)) res = -res;
886
887 return (res);
888 }
889
gmt_k0(struct GMT_CTRL * GMT,double x)890 double gmt_k0 (struct GMT_CTRL *GMT, double x) {
891 /* Modified from code in Press et al. */
892 double y, z, res;
893 gmt_M_unused(GMT);
894
895 if (x < 0.0) x = -x;
896
897 if (x <= 2.0) {
898 y = 0.25 * x * x;
899 z = x * x / 14.0625;
900 res = (-log(0.5*x) * (1.0 + z * (3.5156229 + z * (3.0899424 + z * (1.2067492 + z * (0.2659732 + z * (0.360768e-1 + z * 0.45813e-2))))))) + (-M_EULER + y * (0.42278420 + y * (0.23069756 + y * (0.3488590e-1 + y * (0.262698e-2 + y * (0.10750e-3 + y * 0.74e-5))))));
901 }
902 else {
903 y = 2.0 / x;
904 res = (exp (-x) / sqrt (x)) * (1.25331414 + y * (-0.7832358e-1 + y * (0.2189568e-1 + y * (-0.1062446e-1 + y * (0.587872e-2 + y * (-0.251540e-2 + y * 0.53208e-3))))));
905 }
906 return (res);
907 }
908
gmt_k1(struct GMT_CTRL * GMT,double x)909 double gmt_k1 (struct GMT_CTRL *GMT, double x) {
910 /* Modified Bessel function K1(x) */
911
912 double y, res;
913
914 if (x < 0.0) x = -x;
915 if (x <= 2.0) {
916 y = x * x / 4.0;
917 res = (log (0.5 * x) * gmt_i1 (GMT, x)) + (1.0 / x) * (1.0 + y * (0.15443144 + y * (-0.67278579 + y * (-0.18156897 + y * (-0.01919402 + y * (-0.00110404 - y * 0.00004686))))));
918 }
919 else {
920 y = 2.0 / x;
921 res = (exp (-x) / sqrt (x)) * (1.25331414 + y * (0.23498619 + y * (-0.03655620 + y * (0.01504268 + y * (-0.00780353 + y * (0.00325614 - y * 0.00068245))))));
922 }
923 return (res);
924 }
925
gmt_kn(struct GMT_CTRL * GMT,unsigned int n,double x)926 double gmt_kn (struct GMT_CTRL *GMT, unsigned int n, double x) {
927 /* Modified Bessel function Kn(x) */
928
929 unsigned int j;
930 double bkm, bk, bkp, tox;
931
932 if (n == 0) return (gmt_k0 (GMT, x));
933 if (n == 1) return (gmt_k1 (GMT, x));
934
935 tox = 2.0 / x;
936 bkm = gmt_k0 (GMT, x);
937 bk = gmt_k1 (GMT, x);
938 for (j = 1; j <= (n-1); j++) {
939 bkp = bkm + j * tox * bk;
940 bkm = bk;
941 bk = bkp;
942 }
943
944 return (bk);
945 }
946
gmt_plm(struct GMT_CTRL * GMT,int l,int m,double x)947 double gmt_plm (struct GMT_CTRL *GMT, int l, int m, double x) {
948 /* Unnormalized associated Legendre polynomial of degree l and order m, including
949 * Condon-Shortley phase (-1)^m */
950 double fact, pll = 0, pmm, pmmp1, somx2;
951 int i, ll;
952
953 /* x is cosine of colatitude and must be -1 <= x <= +1 */
954 if (fabs(x) > 1.0) {
955 GMT_Report (GMT->parent, GMT_MSG_WARNING, "|x| > 1.0 in gmt_plm\n");
956 return (GMT->session.d_NaN);
957 }
958
959 if (m < 0 || m > l) {
960 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmt_plm requires 0 <= m <= l.\n");
961 return (GMT->session.d_NaN);
962 }
963
964 pmm = 1.0;
965 if (m > 0) {
966 somx2 = d_sqrt ((1.0 - x) * (1.0 + x));
967 fact = 1.0;
968 /* This loop used to go to i < m; corrected to i <= m by WHFS */
969 for (i = 1; i <= m; i++) {
970 pmm *= -fact * somx2;
971 fact += 2.0;
972 }
973 }
974 if (l == m) return (pmm);
975
976 pmmp1 = x * (2*m + 1) * pmm;
977 if (l == (m + 1)) return (pmmp1);
978
979 for (ll = (m+2); ll <= l; ll++) {
980 pll = (x * (2*ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
981 pmm = pmmp1;
982 pmmp1 = pll;
983 }
984 return (pll);
985 }
986
987 #ifndef M_SQRTPI
988 #define M_SQRTPI 1.77245385090551602729
989 #endif
990
gmt_plm_bar(struct GMT_CTRL * GMT,int l,int m,double x,bool ortho)991 double gmt_plm_bar (struct GMT_CTRL *GMT, int l, int m, double x, bool ortho) {
992 /* This function computes the normalized associated Legendre function of x for degree
993 * l and order m. x must be in the range [-1;1] and 0 <= |m| <= l.
994 * The routine is largely based on the second modified forward column method described in
995 * Holmes and Featherstone (2002). It is stable up to very high degree and order
996 * (at least 3000). This is achieved by individually computing the sectorials to P_m,m
997 * and then iterate up to P_l,m divided by P_m,m and the scale factor 1e280.
998 * Eventually, the result is multiplied again with these two terms.
999 *
1000 * When ortho=false, normalization is according to the geophysical convention.
1001 * - The Condon-Shortley phase (-1)^m is NOT included.
1002 * - The normalization factor is sqrt(k*(2l+1)*(l-m)!/(l+m)!) with k=2 except for m=0: k=1.
1003 * - The integral of Plm**2 over [-1;1] is 2*k.
1004 * - The integrals of (Plm*cos(m*lon))**2 and (Plm*sin(m*lon))**2 over a sphere are 4 pi.
1005 *
1006 * When ortho=true, the results are orthonormalized.
1007 * - The Condon-Shortley phase (-1)^m is NOT included.
1008 * - The normalization factor is sqrt((2l+1)*(l-m)!/(l+m)!/(4*pi)).
1009 * - The integral of Plm**2 over [-1;1] is 1/(2*pi).
1010 * - The integral of (Plm*exp(i*m*lon))**2 over a sphere is 1.
1011 *
1012 * When called with -m, the Condon-Shortley phase will be included.
1013 *
1014 * Note that the orthonormalized form produces an integral of 1 when using imaginary terms
1015 * for longitude. In geophysics, the imaginary components are split into cosine and sine terms
1016 * EACH of which have an average power of 1 over the sphere (i.e. an integral of 4 pi).
1017 *
1018 * This routine could be further expanded to produce an array of P_m,m through P_l,m. This
1019 * would be practical if Legendre functions of various degrees are to be computed at the same
1020 * latitude. Also, differentials could be added.
1021 *
1022 * Reference:
1023 * S. A. Holmes and W. E. Featherstone. A unified approach to the Clenshaw summation and the
1024 * recursive computation of very high degree and order normalised associated Legendre functions.
1025 * Journal of Geodesy, 76, 279-299, 2002. doi:10.1007/s00190-002-0216-2.
1026 */
1027 int i;
1028 bool csphase = false;
1029 double scalef = 1.0e280, u, r, pmm, pmm0, pmm1, pmm2;
1030
1031 /* x is cosine of colatitude (sine of latitude) and must be -1 <= x <= +1 */
1032
1033 if (fabs (x) > 1.0) {
1034 GMT_Report (GMT->parent, GMT_MSG_WARNING, "|x| > 1.0 in gmt_plm_bar\n");
1035 return (GMT->session.d_NaN);
1036 }
1037
1038 /* If m is negative, include Condon-Shortley phase */
1039
1040 if (m < 0) {
1041 csphase = true;
1042 m = -m;
1043 }
1044
1045 if (m > l) {
1046 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmt_plm_bar requires 0 <= m <= l.\n");
1047 return (GMT->session.d_NaN);
1048 }
1049
1050 /* u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1 */
1051
1052 u = d_sqrt ((1.0 - x)*(1.0 + x));
1053
1054 /* Initialize P_00 = 1 and compute up to P_mm using recurrence relation.
1055 The result is a normalisation factor: sqrt((2l+1)(l-m)!/(l+m)!) */
1056
1057 pmm = 1.0;
1058 for (i = 1; i <= m; i++) pmm = d_sqrt (1.0 + 0.5/i) * u * pmm;
1059
1060 /* If orthonormalization is requested: multiply by sqrt(1/4pi)
1061 In case of geophysical conversion : multiply by sqrt(2-delta_0m) */
1062
1063 if (ortho)
1064 pmm /= (M_SQRT2 * M_SQRTPI);
1065 else if (m != 0)
1066 pmm *= M_SQRT2;
1067
1068 /* If C-S phase is requested, apply it now */
1069
1070 if ((m & 1) && csphase) pmm = -pmm;
1071
1072 /* If l==m, we are done */
1073
1074 if (l == m) return (pmm);
1075
1076 /* In the next section all P_l,m are divided by (P_m,m * scalef).
1077 First compute P_m+1,m / P_m,m */
1078
1079 pmm0 = 1.0/scalef;
1080 pmm1 = pmm0 * x * d_sqrt ((double)(2*m + 3));
1081
1082 /* Use second modified column forward recurrence relation to compute P_m+2,m / P_m,m */
1083
1084 for (i = m+2; i <= l; i++) {
1085 r = (2*i+1.0) / (i+m) / (i-m);
1086 pmm2 = x * pmm1 * d_sqrt(r*(2*i-1)) - pmm0 * d_sqrt(r*(i-m-1)*(i+m-1)/(2*i-3));
1087 pmm0 = pmm1;
1088 pmm1 = pmm2;
1089 }
1090
1091 /* Return P_l,m */
1092
1093 pmm1 *= pmm;
1094 pmm1 *= scalef;
1095 return (pmm1);
1096 }
1097
gmt_plm_bar_all(struct GMT_CTRL * GMT,int lmax,double x,bool ortho,double * plm)1098 void gmt_plm_bar_all (struct GMT_CTRL *GMT, int lmax, double x, bool ortho, double *plm) {
1099 /* This function computes the normalized associated Legendre function of x for all degrees
1100 * l <= lmax and all orders m <= l. x must be in the range [-1;1] and 0 <= |m| <= l.
1101 * The routine is largely based on the second modified forward column method described in
1102 * Holmes and Featherstone (2002). It is stable up to very high degree and order
1103 * (at least 3000). This is achieved by individually computing the sectorials to P_m,m
1104 * and then iterate up to P_l,m divided by P_m,m and the scale factor 1e280.
1105 * Eventually, the result is multiplied again with these two terms.
1106 *
1107 * When ortho=false, normalization is according to the geophysical convention.
1108 * - The Condon-Shortley phase (-1)^m is NOT included.
1109 * - The normalization factor is sqrt(k*(2l+1)*(l-m)!/(l+m)!) with k=2 except for m=0: k=1.
1110 * - The integral of Plm**2 over [-1;1] is 2*k.
1111 * - The integrals of (Plm*cos(m*lon))**2 and (Plm*sin(m*lon))**2 over a sphere are 4 pi.
1112 *
1113 * When ortho=true, the results are orthonormalized.
1114 * - The Condon-Shortley phase (-1)^m is NOT included.
1115 * - The normalization factor is sqrt((2l+1)*(l-m)!/(l+m)!/(4*pi)).
1116 * - The integral of Plm**2 over [-1;1] is 1/(2*pi).
1117 * - The integral of (Plm*exp(i*m*lon))**2 over a sphere is 1.
1118 *
1119 * When called with -lmax, the Condon-Shortley phase will be included.
1120 *
1121 * Note that the orthonormalized form produces an integral of 1 when using imaginary terms
1122 * for longitude. In geophysics, the imaginary components are split into cosine and sine terms
1123 * EACH of which have an average power of 1 over the sphere (i.e. an integral of 4 pi).
1124 *
1125 * This routine produces an array of all values of P_l,m(x) in the array plm, in the order
1126 * P_0,0 P_1,0 P_1,1 P_2,0 P_2,1 P_2,2, etc. That means that plm[j] refers to:
1127 * tesserals P_l,m when j = l * (l+1) / 2 + m
1128 * zonals P_l,0 when j = l * (l+1) / 2
1129 * sectorials P_m,m when j = m * (m+3) / 2
1130 *
1131 * Reference:
1132 * S. A. Holmes and W. E. Featherstone. A unified approach to the Clenshaw summation and the
1133 * recursive computation of very high degree and order normalised associated Legendre functions.
1134 * Journal of Geodesy, 76, 279-299, 2002. doi:10.1007/s00190-002-0216-2.
1135 */
1136 int l, m, lm, mm;
1137 bool csphase = false;
1138 double scalef = 1.0e280, u, r, pmm, pmms, pmm0, pmm1, pmm2;
1139
1140 /* x is cosine of colatitude (sine of latitude) and must be -1 <= x <= +1 */
1141
1142 if (fabs (x) > 1.0) {
1143 GMT_Report (GMT->parent, GMT_MSG_WARNING, "|x| > 1.0 in gmt_plm_bar_all\n");
1144 return;
1145 }
1146
1147 /* If lmax is negative, include Condon-Shortley phase */
1148
1149 if (lmax < 0) {
1150 csphase = true;
1151 lmax = -lmax;
1152 }
1153
1154 /* u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1 */
1155
1156 u = d_sqrt ((1.0 - x)*(1.0 + x));
1157
1158 /* Initialize P_00 = 1 */
1159 mm = 0;
1160 plm[mm] = pmm = 1.0;
1161
1162 /* Loop over 0 <= m <= lmax */
1163
1164 for (m = 0; m <= lmax; m++, mm += m+1) {
1165
1166 /* Compute up P_mm using recurrence relation.
1167 The result is a normalisation factor: sqrt((2l+1)(l-m)!/(l+m)!) */
1168
1169 if (m != 0) pmm *= d_sqrt (1.0 + 0.5/m) * u;
1170
1171 /* If orthonormalization is requested: multiply by sqrt(1/4pi)
1172 In case of geophysical conversion : multiply by sqrt(2-delta_0m) */
1173
1174 if (ortho)
1175 plm[mm] = pmm / (M_SQRT2 * M_SQRTPI);
1176 else if (m != 0)
1177 plm[mm] = pmm * M_SQRT2;
1178
1179 /* If C-S phase is requested, apply it now */
1180
1181 if ((m & 1) && csphase) plm[mm] = -plm[mm];
1182
1183 /* In the next section all P_l,m are divided by (P_m,m * scalef).
1184 First compute P_m+1,m / P_m,m */
1185
1186 pmms = plm[mm] * scalef;
1187 pmm0 = 1.0 / scalef;
1188 pmm1 = pmm0 * x * d_sqrt ((double)(2*m + 3));
1189 lm = mm + m + 1;
1190 plm[lm] = pmm1 * pmms;
1191
1192 /* Use second modified column forward recurrence relation to compute P_m+2,m / P_m,m */
1193
1194 for (l = m+2; l <= lmax; l++) {
1195 r = (2*l+1.0) / (l+m) / (l-m);
1196 pmm2 = x * pmm1 * d_sqrt(r*(2*l-1)) - pmm0 * d_sqrt(r*(l-m-1)*(l+m-1)/(2*l-3));
1197 pmm0 = pmm1;
1198 pmm1 = pmm2;
1199 lm += l;
1200 plm[lm] = pmm1 * pmms;
1201 }
1202 }
1203 }
1204
1205 /* gmt_sinc (x) calculates the sinc function */
1206
gmt_sinc(struct GMT_CTRL * GMT,double x)1207 double gmt_sinc (struct GMT_CTRL *GMT, double x) {
1208 gmt_M_unused(GMT);
1209 if (x == 0.0) return (1.0);
1210 x *= M_PI;
1211 return (sin (x) / x);
1212 }
1213
1214 /* gmt_factorial (n) calculates the factorial n! */
1215
gmt_factorial(struct GMT_CTRL * GMT,int n)1216 double gmt_factorial (struct GMT_CTRL *GMT, int n) {
1217 static int ntop = 10; /* Initial portion filled in below */
1218 static double a[33] = {1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 3628800.0};
1219 int i;
1220
1221 if (n < 0) {
1222 GMT_Report (GMT->parent, GMT_MSG_WARNING, "n < 0 in gmt_factorial(n)\n");
1223 return (GMT->session.d_NaN);
1224 /* This could be set to return 0 without warning, to facilitate
1225 sums over binomial coefficients, if desired. -whfs */
1226 }
1227 if (n > 32) return (exp(gmtstat_ln_gamma (GMT, n+1.0)));
1228 while (ntop < n) {
1229 i = ntop++;
1230 a[ntop] = a[i] * ntop;
1231 }
1232 return (a[n]);
1233 }
1234
gmt_permutation(struct GMT_CTRL * GMT,int n,int r)1235 double gmt_permutation (struct GMT_CTRL *GMT, int n, int r) {
1236 /* Compute Permutations n_P_r */
1237 if (n < 0 || r < 0 || r > n) {
1238 GMT_Report (GMT->parent, GMT_MSG_WARNING, "n < 0 or r < 0 or r > n in gmt_permutation(n,r)\n");
1239 return (GMT->session.d_NaN);
1240 }
1241 return (floor (0.5 + exp (gmtstat_factln (GMT, n) - gmtstat_factln (GMT, n-r))));
1242 }
1243
gmt_combination(struct GMT_CTRL * GMT,int n,int r)1244 double gmt_combination (struct GMT_CTRL *GMT, int n, int r) {
1245 /* Compute Combinations n_C_r */
1246 if (n < 0 || r < 0 || r > n) {
1247 GMT_Report (GMT->parent, GMT_MSG_WARNING, "n < 0 or r < 0 or r > n in gmt_combination(n,r)\n");
1248 return (GMT->session.d_NaN);
1249 }
1250 return (floor (0.5 + exp (gmtstat_factln (GMT, n) - gmtstat_factln (GMT, r) - gmtstat_factln (GMT, n-r))));
1251 }
1252
gmt_dilog(struct GMT_CTRL * GMT,double x)1253 double gmt_dilog (struct GMT_CTRL *GMT, double x) {
1254 /* Compute dilog(x) (defined for x >= 0) by the method of Parker's
1255 Appendix A of his Geophysical Inverse Theory. The function
1256 is needed for x in the range 0 <= x <= 1 when solving the
1257 spherical spline interpolation in section 2.07 of Parker.
1258
1259 I tested this for x in range 0 to 1 by reproducing Figure
1260 2.07c of Parker's book. I also tested that the result was
1261 smooth for x out to 25. I also checked d[dilog(x)]/dx
1262 obtained numerically from this routine against an analytic
1263 expression, and I found the smallest angular separation
1264 between two points on a sphere such that their Gram
1265 matrix expression (see Parker) was different from 1; this
1266 is smallest epsilon such that dilog(epsilon) - dilog(0)
1267 != 0.0. This turned out to be very much smaller than I
1268 would have guessed from the apparent e-15 accuracy of
1269 Parker's polynomial. So I think this is very accurate.
1270
1271 If this is called with x < 0, it returns dilog(0) and
1272 prints a warning. It could also be set to return NaN.
1273 I did this because in applications it might happen that
1274 an x which should be zero is computed to be slightly off,
1275 in which case the desired result is for x = 0.
1276
1277 23 Oct 1996 WHFS. */
1278
1279 double pisqon6, y, ysq, z;
1280
1281 if (x < -GMT_CONV8_LIMIT) return (GMT->session.d_NaN); /* Tolerate minor slop before we are outside domain */
1282
1283 pisqon6 = M_PI * M_PI / 6.0;
1284 if (x <= 0.0) return (pisqon6); /* Meaning -GMT_CONV8_LIMIT < x <= 0 */
1285
1286 if (x < 0.5) {
1287 y = -log (1.0 - x);
1288 ysq = y * y;
1289 z = y * (1.0 + y * (-0.25 + y * (0.027777777777213 + ysq * (-2.7777776990e-04 + ysq * (4.724071696e-06 + ysq * (-9.1764954e-08 + 1.798670e-09 * ysq))))));
1290 return (pisqon6 - z + y * log (x));
1291 }
1292 if (x < 2.0) {
1293 y = -log (x);
1294 ysq = y * y;
1295 z = y * (1.0 + y * (-0.25 + y * (0.027777777777213 + ysq * (-2.7777776990e-04 + ysq * (4.724071696e-06 + ysq * (-9.1764954e-08 + 1.798670e-09 * ysq))))));
1296 return (z);
1297 }
1298 y = log (x);
1299 ysq = y * y;
1300 z = y * (1.0 + y * (-0.25 + y * (0.027777777777213 + ysq * (-2.7777776990e-04 + ysq * (4.724071696e-06 + ysq * (-9.1764954e-08 + 1.798670e-09 * ysq))))));
1301 return (-z - 0.5 * ysq);
1302 }
1303
1304 #ifndef M_2_SQRTPI
1305 #define M_2_SQRTPI 1.12837916709551257390
1306 #endif
1307
gmt_erfinv(struct GMT_CTRL * GMT,double y)1308 double gmt_erfinv (struct GMT_CTRL *GMT, double y) {
1309 double x = 0.0, fy, z;
1310
1311 /* Misc. efficients for expansion */
1312
1313 static double a[4] = {0.886226899, -1.645349621, 0.914624893, -0.140543331};
1314 static double b[4] = {-2.118377725, 1.442710462, -0.329097515, 0.012229801};
1315 static double c[4] = {-1.970840454, -1.624906493, 3.429567803, 1.641345311};
1316 static double d[2] = {3.543889200, 1.637067800};
1317
1318 fy = fabs (y);
1319
1320 if (fy > 1.0) return (GMT->session.d_NaN); /* Outside domain */
1321
1322 if (doubleAlmostEqual (fy, 1.0))
1323 return (copysign (DBL_MAX, y)); /* Close to +- Inf */
1324
1325 if (y > 0.7) { /* Near upper range */
1326 z = sqrt (-log (0.5 * (1.0 - y)));
1327 x = (((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
1328 }
1329 else if (y < -0.7) { /* Near lower range */
1330 z = sqrt (-log (0.5 * (1.0 + y)));
1331 x = -(((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0);
1332 }
1333 else { /* Central range */
1334 z = y * y;
1335 x = y * (((a[3] * z + a[2]) * z + a[1]) * z + a[0]) / ((((b[3] * z + b[2]) * z + b[1]) * z + b[0]) * z + 1.0);
1336 }
1337
1338 /* Apply two steps of Newton-Raphson correction to improve accuracy */
1339
1340 x -= (erf (x) - y) / (M_2_SQRTPI * exp (-x * x));
1341 x -= (erf (x) - y) / (M_2_SQRTPI * exp (-x * x));
1342
1343 return (x);
1344 }
1345
gmt_f_pdf(struct GMT_CTRL * GMT,double F,uint64_t nu1,uint64_t nu2)1346 double gmt_f_pdf (struct GMT_CTRL *GMT, double F, uint64_t nu1, uint64_t nu2) {
1347 /* Probability density distribution for F */
1348 double y;
1349
1350 y = (F < GMT_CONV8_LIMIT) ? 0.0 : sqrt (pow (nu1 * F, (double)nu1) * pow ((double)nu2, (double)nu2) / pow (nu1 * F + nu2, (double)(nu1+nu2))) / (F * gmtstat_beta (GMT, 0.5*nu1, 0.5*nu2));
1351 return (y);
1352 }
1353
gmt_f_cdf(struct GMT_CTRL * GMT,double F,uint64_t nu1,uint64_t nu2)1354 double gmt_f_cdf (struct GMT_CTRL *GMT, double F, uint64_t nu1, uint64_t nu2) {
1355 /* Cumulative probability density distribution for F */
1356 double y = 0.0;
1357
1358 if (F < GMT_CONV8_LIMIT) return 0.0;
1359 gmtstat_inc_beta (GMT, 0.5*nu1, 0.5*nu2, F*nu1/(F*nu1+nu2), &y);
1360
1361 return (y);
1362 }
1363
gmt_t_pdf(struct GMT_CTRL * GMT,double t,uint64_t nu)1364 double gmt_t_pdf (struct GMT_CTRL *GMT, double t, uint64_t nu) {
1365 /* Probability density distribution for Student t */
1366 double y, n = nu + 1.0, g1 = 0.0, g2 = 0.0;
1367
1368 gmtstat_ln_gamma_r (GMT, 0.5*n, &g1);
1369 gmtstat_ln_gamma_r (GMT, 0.5*nu, &g2);
1370 y = exp (g1 - g2) * pow (1.0 + t*t/nu, -0.5*n) / sqrt (M_PI * nu);
1371 return (y);
1372 }
1373
gmt_t_cdf(struct GMT_CTRL * GMT,double t,uint64_t nu)1374 double gmt_t_cdf (struct GMT_CTRL *GMT, double t, uint64_t nu) {
1375 double p;
1376 /* Cumulative Student t distribution */
1377 gmt_M_unused(GMT);
1378 gmtstat_student_t_a (GMT, fabs (t), nu, &p);
1379 p = 0.5 * p + 0.5;
1380 if (t < 0.0)
1381 p = 1.0 - p;
1382 return (p);
1383 }
1384
gmt_weibull_pdf(struct GMT_CTRL * GMT,double x,double scale,double shape)1385 double gmt_weibull_pdf (struct GMT_CTRL *GMT, double x, double scale, double shape) {
1386 double p, z;
1387 gmt_M_unused(GMT);
1388 /* Weibull distribution */
1389 if (x < 0.0) return 0.0;
1390 z = x / scale;
1391 p = (shape/scale) * pow (z, shape-1.0) * exp (-pow(z,shape));
1392 return (p);
1393 }
1394
gmt_weibull_cdf(struct GMT_CTRL * GMT,double x,double scale,double shape)1395 double gmt_weibull_cdf (struct GMT_CTRL *GMT, double x, double scale, double shape) {
1396 double p, z;
1397 gmt_M_unused(GMT);
1398 /* Cumulative Weibull distribution */
1399 if (x < 0.0) return 0.0;
1400 z = x / scale;
1401 p = 1.0 - exp (-pow (z, shape));
1402 return (p);
1403 }
1404
gmt_weibull_crit(struct GMT_CTRL * GMT,double p,double scale,double shape)1405 double gmt_weibull_crit (struct GMT_CTRL *GMT, double p, double scale, double shape) {
1406 double z;
1407 gmt_M_unused(GMT);
1408 /* Critical values for Weibull distribution */
1409 z = scale * pow (-log (1.0 - p), 1.0/shape);
1410 return (z);
1411 }
1412
gmt_binom_pdf(struct GMT_CTRL * GMT,uint64_t x,uint64_t n,double p)1413 double gmt_binom_pdf (struct GMT_CTRL *GMT, uint64_t x, uint64_t n, double p) {
1414 double c;
1415 /* Binomial distribution */
1416 c = gmt_combination (GMT, (int)n, (int)x) * pow (p, (int)x) * pow (1.0-p, (int)(n-x));
1417 return (c);
1418 }
1419
gmt_binom_cdf(struct GMT_CTRL * GMT,uint64_t x,uint64_t n,double p)1420 double gmt_binom_cdf (struct GMT_CTRL *GMT, uint64_t x, uint64_t n, double p) {
1421 /* Cumulative Binomial distribution */
1422 double c = 0.0;
1423 if (n > 12) /* Use Numerical Recipes fast way for larger n */
1424 gmtstat_inc_beta (GMT, (double)x, (double)(n-x+1), p, &c);
1425 else { /* Do the sum instead */
1426 uint64_t j;
1427 for (j = 0; j <= x; j++)
1428 c += gmt_binom_pdf (GMT, j, n, p);
1429 }
1430 return (c);
1431 }
1432
gmt_vonmises_pdf(struct GMT_CTRL * GMT,double x,double mu,double kappa)1433 double gmt_vonmises_pdf (struct GMT_CTRL *GMT, double x, double mu, double kappa) {
1434 /* Von Mises 1-D probability density function */
1435 double p;
1436 gmt_M_unused(GMT);
1437 /* Von Mises distribution, expects x and mu in degrees */
1438 p = exp (kappa * cosd (x - mu)) / (TWO_PI * gmt_i0 (GMT, kappa));
1439 return (p);
1440 }
1441
gmt_fisher_pdf(struct GMT_CTRL * GMT,double plon,double plat,double lon,double lat,double kappa)1442 double gmt_fisher_pdf (struct GMT_CTRL *GMT, double plon, double plat, double lon, double lat, double kappa) {
1443 /* Fisher 3-D probability density function centered on (plon, plat) evaluated at
1444 * another point (lon, lat), given kappa */
1445 double p, cos_psi;
1446 cos_psi = gmtlib_great_circle_dist_cos (GMT, plon, plat, lon, lat);
1447 p = kappa * exp (kappa * cos_psi) / (2.0 * TWO_PI * sinh (kappa));
1448 return (p);
1449 }
1450
gmt_zdist(struct GMT_CTRL * GMT,double x)1451 double gmt_zdist (struct GMT_CTRL *GMT, double x) {
1452 /* Cumulative Normal (z) distribution */
1453 gmt_M_unused(GMT);
1454 return (0.5 * (erf (x / M_SQRT2) + 1.0));
1455 }
1456
gmt_zcrit(struct GMT_CTRL * GMT,double alpha)1457 double gmt_zcrit (struct GMT_CTRL *GMT, double alpha) {
1458 double sign;
1459 /* Critical values for Normal (z) distribution */
1460
1461 /* Simple since z_a = M_SQRT2 * erf^-1 (1-a) */
1462
1463 if (alpha > 0.5) { /* right tail */
1464 alpha = (1.0 - alpha) * 2.0;
1465 sign = 1.0;
1466 }
1467 else { /* left tail */
1468 alpha *= 2.0;
1469 sign = -1.0;
1470 }
1471
1472 return (sign * M_SQRT2 * gmt_erfinv (GMT, 1.0 - alpha));
1473 }
1474
gmt_tcrit(struct GMT_CTRL * GMT,double alpha,double nu)1475 double gmt_tcrit (struct GMT_CTRL *GMT, double alpha, double nu) {
1476 /* Critical values for Student t-distribution */
1477
1478 int NU;
1479 bool done;
1480 double t_low, t_high, t_mid = 0.0, p_high, p_mid, p, sign;
1481
1482 if (alpha > 0.5) { /* right tail */
1483 p = 1 - (1 - alpha) * 2.0;
1484 sign = 1.0;
1485 }
1486 else {
1487 p = 1 - alpha * 2.0;
1488 sign = -1.0;
1489 }
1490 t_low = gmt_zcrit (GMT, alpha);
1491 if (gmt_M_is_dinf (t_low) || gmt_M_is_dnan (t_low)) return (t_low);
1492 t_high = 5.0;
1493 NU = irint(nu);
1494 gmtstat_student_t_a (GMT, t_high, NU, &p_high);
1495 while (p_high < p) { /* Must pick higher starting point */
1496 t_high *= 2.0;
1497 gmtstat_student_t_a (GMT, t_high, NU, &p_high);
1498 }
1499
1500 /* Now, (t_low, p_low) and (t_high, p_high) are bracketing the desired (t,p) */
1501
1502 done = false;
1503 while (!done) {
1504 t_mid = 0.5 * (t_low + t_high);
1505 gmtstat_student_t_a (GMT, t_mid, NU, &p_mid);
1506 if (doubleAlmostEqualZero (p_mid, p)) {
1507 done = true;
1508 }
1509 else if (p_mid > p) { /* Too high */
1510 t_high = t_mid;
1511 }
1512 else { /* p_mid < p */
1513 t_low = t_mid;
1514 }
1515 }
1516 return (sign * t_mid);
1517 }
1518
gmt_chi2_pdf(struct GMT_CTRL * GMT,double c,uint64_t nu)1519 double gmt_chi2_pdf (struct GMT_CTRL *GMT, double c, uint64_t nu) {
1520 /* Probability density distribution for chi-squared */
1521 double g = 0.0, y;
1522
1523 gmtstat_ln_gamma_r (GMT, 0.5*nu, &g);
1524 y = pow (c, 0.5*nu - 1.0) * exp (-0.5 * c - g) / pow (2.0, 0.5 * nu);
1525 return (y);
1526 }
1527
gmt_chi2crit(struct GMT_CTRL * GMT,double alpha,double nu)1528 double gmt_chi2crit (struct GMT_CTRL *GMT, double alpha, double nu) {
1529 /* Critical values for Chi^2-distribution */
1530
1531 bool done;
1532 double chi2_low, chi2_high, chi2_mid = 0.0, p_high, p_mid, p;
1533
1534 p = 1.0 - alpha;
1535 chi2_low = 0.0;
1536 chi2_high = 5.0;
1537 gmt_chi2 (GMT, chi2_high, nu, &p_high);
1538 while (p_high > p) { /* Must pick higher starting point */
1539 chi2_high *= 2.0;
1540 gmt_chi2 (GMT, chi2_high, nu, &p_high);
1541 }
1542
1543 /* Now, (chi2_low, p_low) and (chi2_high, p_high) are bracketing the desired (chi2,p) */
1544
1545 done = false;
1546 while (!done) {
1547 chi2_mid = 0.5 * (chi2_low + chi2_high);
1548 gmt_chi2 (GMT, chi2_mid, nu, &p_mid);
1549 if (doubleAlmostEqualZero (p_mid, p)) {
1550 done = true;
1551 }
1552 else if (p_mid < p) { /* Too high */
1553 chi2_high = chi2_mid;
1554 }
1555 else { /* p_mid > p */
1556 chi2_low = chi2_mid;
1557 }
1558 }
1559 return (chi2_mid);
1560 }
1561
gmt_Fcrit(struct GMT_CTRL * GMT,double alpha,double nu1,double nu2)1562 double gmt_Fcrit (struct GMT_CTRL *GMT, double alpha, double nu1, double nu2) {
1563 /* Critical values for F-distribution */
1564
1565 int NU1, NU2;
1566 bool done;
1567 double F_low, F_high, F_mid = 0.0, p_high, p_mid, p, chisq1, chisq2;
1568
1569 F_high = 5.0;
1570 p = 1.0 - alpha;
1571 F_low = 0.0;
1572 gmtstat_F_to_ch1_ch2 (GMT, F_high, nu1, nu2, &chisq1, &chisq2);
1573 NU1 = irint (nu1);
1574 NU2 = irint (nu2);
1575 gmtstat_f_q (GMT, chisq1, NU1, chisq2, NU2, &p_high);
1576 while (p_high > p) { /* Must pick higher starting point */
1577 F_high *= 2.0;
1578 gmtstat_F_to_ch1_ch2 (GMT, F_high, nu1, nu2, &chisq1, &chisq2);
1579 gmtstat_f_q (GMT, chisq1, NU1, chisq2, NU2, &p_high);
1580 }
1581
1582 /* Now, (F_low, p_low) and (F_high, p_high) are bracketing the desired (F,p) */
1583
1584 done = false;
1585 while (!done) {
1586 F_mid = 0.5 * (F_low + F_high);
1587 gmtstat_F_to_ch1_ch2 (GMT, F_mid, nu1, nu2, &chisq1, &chisq2);
1588 gmtstat_f_q (GMT, chisq1, NU1, chisq2, NU2, &p_mid);
1589 if (doubleAlmostEqualZero (p_mid, p)) {
1590 done = true;
1591 }
1592 else if (p_mid < p) { /* Too high */
1593 F_high = F_mid;
1594 }
1595 else { /* p_mid > p */
1596 F_low = F_mid;
1597 }
1598 }
1599 return (F_mid);
1600 }
1601
gmtstat_mix64(uint64_t a,uint64_t b,uint64_t c)1602 static inline uint64_t gmtstat_mix64 (uint64_t a, uint64_t b, uint64_t c) {
1603 /* Mix 3 64-bit values, from lookup8.c by Bob Jenkins
1604 * (http://burtleburtle.net/bob/index.html) */
1605 a -= b; a -= c; a ^= (c>>43);
1606 b -= c; b -= a; b ^= (a<<9);
1607 c -= a; c -= b; c ^= (b>>8);
1608 a -= b; a -= c; a ^= (c>>38);
1609 b -= c; b -= a; b ^= (a<<23);
1610 c -= a; c -= b; c ^= (b>>5);
1611 a -= b; a -= c; a ^= (c>>35);
1612 b -= c; b -= a; b ^= (a<<49);
1613 c -= a; c -= b; c ^= (b>>11);
1614 a -= b; a -= c; a ^= (c>>12);
1615 b -= c; b -= a; b ^= (a<<18);
1616 c -= a; c -= b; c ^= (b>>22);
1617 return c;
1618 }
1619
gmt_rand(struct GMT_CTRL * GMT)1620 double gmt_rand (struct GMT_CTRL *GMT) {
1621 /* Uniform random number generator. Will return values
1622 * x so that 0.0 < x < 1.0 occurs with equal probability. */
1623 static unsigned seed = 0;
1624 double random_val;
1625
1626 while (seed == 0) { /* repeat in case of unsigned overflow */
1627 /* Initialize random seed, idea from Jonathan
1628 * Wright (http://stackoverflow.com/q/322938) */
1629 seed = (unsigned) gmtstat_mix64 (clock(), time(NULL), getpid());
1630 srand (seed);
1631 }
1632
1633 /* generate random number */
1634 random_val = rand () / (double) RAND_MAX;
1635
1636 if (random_val == 0.0 || random_val >= 1.0)
1637 /* Ensure range (0.0,1.0) */
1638 return gmt_rand (GMT);
1639
1640 return random_val;
1641 }
1642
gmt_nrand(struct GMT_CTRL * GMT)1643 double gmt_nrand (struct GMT_CTRL *GMT) {
1644 /* Gaussian random number generator based on gasdev of
1645 * Press et al, Numerical Recipes, 2nd edition. Will
1646 * return values that have zero mean and unit variance.
1647 */
1648
1649 static int iset = 0;
1650 static double gset;
1651 double fac, r, v1, v2;
1652
1653 if (iset == 0) { /* We don't have an extra deviate handy, so */
1654 do {
1655 v1 = 2.0 * gmt_rand (GMT) - 1.0; /* Pick two uniform numbers in the -1/1/-1/1 square */
1656 v2 = 2.0 * gmt_rand (GMT) - 1.0;
1657 r = v1 * v1 + v2 * v2;
1658 } while (r >= 1.0 || r == 0.0); /* Keep trying until v1,v2 is inside unit circle */
1659
1660 fac = sqrt (-2.0 * log (r) / r);
1661
1662 /* Now make Box-Muller transformation to get two normal deviates. Return
1663 * one and save the other for the next time gmt_nrand is called */
1664
1665 gset = v1 * fac;
1666 iset = 1; /* Set flag for next time */
1667 return (v2 * fac);
1668 }
1669 else {
1670 iset = 0; /* Take old value, reset flag */
1671 return (gset);
1672 }
1673 }
1674
gmt_lrand(struct GMT_CTRL * GMT)1675 double gmt_lrand (struct GMT_CTRL *GMT) {
1676 /* Laplace random number generator. As nrand, it will
1677 * return values that have zero mean and unit variance.
1678 */
1679
1680 double rand_0_to_1;
1681
1682 rand_0_to_1 = gmt_rand (GMT); /* Gives uniformly distributed random values in 0-1 range */
1683 return (((rand_0_to_1 <= 0.5) ? log (2.0 * rand_0_to_1) : -log (2.0 * (1.0 - rand_0_to_1))) / M_SQRT2);
1684 }
1685
gmt_chi2(struct GMT_CTRL * GMT,double chi2,double nu,double * prob)1686 void gmt_chi2 (struct GMT_CTRL *GMT, double chi2, double nu, double *prob) {
1687 /* Evaluate probability that chi2 will exceed the
1688 * theoretical chi2 by chance. */
1689
1690 *prob = gmtstat_gammq (GMT, 0.5 * nu, 0.5 * chi2);
1691 }
1692
gmt_poissonpdf(struct GMT_CTRL * GMT,double k,double lambda)1693 double gmt_poissonpdf (struct GMT_CTRL *GMT, double k, double lambda) {
1694 /* Evaluate PDF for Poisson Distribution */
1695 return pow (lambda, k) * exp (-lambda) / gmt_factorial (GMT, irint (k));
1696 }
1697
gmt_poisson_cdf(struct GMT_CTRL * GMT,double k,double mu,double * prob)1698 void gmt_poisson_cdf (struct GMT_CTRL *GMT, double k, double mu, double *prob) {
1699 /* evaluate Cumulative Poisson Distribution */
1700
1701 *prob = (k == 0.0) ? exp (-mu) : gmtstat_gammq (GMT, k+1.0, mu);
1702 }
1703
gmt_mean_and_std(struct GMT_CTRL * GMT,double * x,uint64_t n,double * std)1704 double gmt_mean_and_std (struct GMT_CTRL *GMT, double *x, uint64_t n, double *std) {
1705 /* Return the standard deviation of the non-NaN values in x */
1706 uint64_t k, m = 0;
1707 double dx, mean = 0.0, sum2 = 0.0;
1708 for (k = 0; k < n; k++) { /* Use Welford (1962) algorithm to compute mean and corrected sum of squares */
1709 if (gmt_M_is_dnan (x[k])) continue;
1710 m++;
1711 dx = x[k] - mean;
1712 mean += dx / m;
1713 sum2 += dx * (x[k] - mean);
1714 }
1715 *std = (m > 1) ? sqrt (sum2 / (m-1.0)) : GMT->session.d_NaN;
1716 return ((m) ? mean : GMT->session.d_NaN);
1717 }
1718
gmt_std_weighted(struct GMT_CTRL * GMT,double * x,double * w,double wmean,uint64_t n)1719 double gmt_std_weighted (struct GMT_CTRL *GMT, double *x, double *w, double wmean, uint64_t n) {
1720 /* Return the weighted standard deviation of the non-NaN values in x.
1721 * If w == NULL then a regular std is computed. */
1722 uint64_t k, m = 0;
1723 double dx, sumw = 0.0, sum2 = 0.0, wk = 1.0;
1724 for (k = 0; k < n; k++) {
1725 if (gmt_M_is_dnan (x[k])) continue;
1726 if (w) {
1727 if (gmt_M_is_dnan (w[k])) continue;
1728 if (gmt_M_is_zero (w[k])) continue;
1729 wk = w[k];
1730 }
1731 dx = x[k] - wmean;
1732 sum2 += wk * dx * dx;
1733 sumw += wk;
1734 m++; /* Number of points with nonzero weights */
1735 }
1736 return (m > 1) ? sqrt (sum2 / (((m-1.0) * sumw) / m)) : GMT->session.d_NaN;
1737 }
1738
gmt_median(struct GMT_CTRL * GMT,double * x,uint64_t n,double xmin,double xmax,double m_initial,double * med)1739 int gmt_median (struct GMT_CTRL *GMT, double *x, uint64_t n, double xmin, double xmax, double m_initial, double *med) {
1740 double lower_bound, upper_bound, m_guess, t_0, t_1, t_middle;
1741 double lub, glb, xx, temp;
1742 uint64_t i;
1743 int64_t n_above, n_below, n_equal, n_lub, n_glb, one; /* These must be signed integers (PW: Why? -> (unsigned)x - (unsigned)y will never become negative) */
1744 int iteration = 0;
1745 bool finished = false;
1746
1747 if (n == 0) { /* No data, so no defined median */
1748 *med = GMT->session.d_NaN;
1749 return (1);
1750 }
1751 if (n == 1) {
1752 *med = x[0];
1753 return (1);
1754 }
1755 if (n == 2) {
1756 *med = 0.5 * (x[0] + x[1]);
1757 return (1);
1758 }
1759
1760 m_guess = m_initial;
1761 lower_bound = xmin;
1762 upper_bound = xmax;
1763 one = 1;
1764 t_0 = 0.0;
1765 t_1 = (double)(n - one);
1766 t_middle = 0.5 * t_1;
1767
1768 do {
1769
1770 n_above = n_below = n_equal = n_lub = n_glb = 0;
1771 lub = xmax;
1772 glb = xmin;
1773
1774 for (i = 0; i < n; i++) {
1775
1776 xx = x[i];
1777 if (xx == m_guess)
1778 n_equal++;
1779 else if (xx > m_guess) {
1780 n_above++;
1781 if (xx < lub) {
1782 lub = xx;
1783 n_lub = one;
1784 }
1785 else if (xx == lub)
1786 n_lub++;
1787 }
1788 else {
1789 n_below++;
1790 if (xx > glb) {
1791 glb = xx;
1792 n_glb = one;
1793 }
1794 else if (xx == glb)
1795 n_glb++;
1796 }
1797 }
1798
1799 iteration++;
1800
1801 /* Now test counts, watch multiple roots, think even/odd: */
1802
1803 if ((int64_abs (n_above - n_below)) <= n_equal) {
1804 *med = (n_equal) ? m_guess : 0.5 * (lub + glb);
1805 finished = true;
1806 }
1807 else if ((int64_abs ((n_above - n_lub) - (n_below + n_equal))) < n_lub) {
1808 *med = lub;
1809 finished = true;
1810 }
1811 else if ((int64_abs ((n_below - n_glb) - (n_above + n_equal))) < n_glb) {
1812 *med = glb;
1813 finished = true;
1814 }
1815 /* Those cases found the median; the next two will forecast a new guess: */
1816
1817 else if (n_above > (n_below + n_equal)) { /* Guess is too low */
1818 lower_bound = m_guess;
1819 t_0 = (double)(n_below + n_equal - one);
1820 temp = lower_bound + (upper_bound - lower_bound) * (t_middle - t_0) / (t_1 - t_0);
1821 m_guess = (temp > lub) ? temp : lub; /* Move guess at least to lub */
1822 }
1823 else if (n_below > (n_above + n_equal)) { /* Guess is too high */
1824 upper_bound = m_guess;
1825 t_1 = (double)(n_below + n_equal - one);
1826 temp = lower_bound + (upper_bound - lower_bound) * (t_middle - t_0) / (t_1 - t_0);
1827 m_guess = (temp < glb) ? temp : glb; /* Move guess at least to glb */
1828 }
1829 else { /* If we get here, I made a mistake! */
1830 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Internal goof in gmt_median; please report to developers!\n");
1831 *med = GMT->session.d_NaN;
1832 return GMT_NOTSET;
1833 }
1834
1835 } while (!finished);
1836
1837 /* That's all, folks! */
1838 return (iteration);
1839 }
1840
gmt_mean_weighted(struct GMT_CTRL * GMT,double * x,double * w,uint64_t n)1841 double gmt_mean_weighted (struct GMT_CTRL *GMT, double *x, double *w, uint64_t n) {
1842 /* Return the weighted mean of x given weights w */
1843 uint64_t k;
1844 double sum_xw = 0.0, sum_w = 0.0;
1845
1846 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined mean */
1847 for (k = 0; k < n; k++) {
1848 sum_w += w[k];
1849 sum_xw += x[k] * w[k];
1850 }
1851 if (sum_w == 0.0) return (GMT->session.d_NaN); /* No weights, so no defined mean */
1852 return (sum_xw / sum_w);
1853 }
1854
gmt_quantile_weighted(struct GMT_CTRL * GMT,struct GMT_OBSERVATION * data,uint64_t n,double quantile)1855 double gmt_quantile_weighted (struct GMT_CTRL *GMT, struct GMT_OBSERVATION *data, uint64_t n, double quantile) {
1856 uint64_t k;
1857 double weight_half = 0.0, weight_count;
1858
1859 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined mode */
1860
1861 /* First sort data on z */
1862
1863 qsort (data, n, sizeof (struct GMT_OBSERVATION), gmtlib_compare_observation);
1864
1865 /* Find weight sum, then get half-value */
1866
1867 for (k = 0; k < n; k++) weight_half += data[k].weight;
1868 weight_half *= quantile; /* Normally quantile = 0.5 hence the name "half" */
1869
1870 /* Determine the point where we hit the desired quantile */
1871
1872 k = 0; weight_count = data[k].weight;
1873 while (weight_count < weight_half) weight_count += data[++k].weight; /* Wind up until weight_count hits the mark */
1874
1875 return ((double)((weight_count == weight_half) ? 0.5 * (data[k].value + data[k+1].value) : data[k].value));
1876 }
1877
gmt_median_weighted(struct GMT_CTRL * GMT,struct GMT_OBSERVATION * data,uint64_t n)1878 double gmt_median_weighted (struct GMT_CTRL *GMT, struct GMT_OBSERVATION *data, uint64_t n) {
1879 /* Just hardwire quantile to 0.5 */
1880 return gmt_quantile_weighted (GMT, data, n, 0.5);
1881 }
1882
gmt_mode_weighted(struct GMT_CTRL * GMT,struct GMT_OBSERVATION * data,uint64_t n)1883 double gmt_mode_weighted (struct GMT_CTRL *GMT, struct GMT_OBSERVATION *data, uint64_t n) {
1884 /* Looks for the "shortest 50%". This means that when the cumulative weight
1885 (y) is plotted against the value (x) then the line between (xi,yi) and
1886 (xj,yj) should be the steepest for any combination where (yj-yi) is 50%
1887 of the total sum of weights */
1888
1889 double top, bottom, wsum, p, p_max, mode;
1890 uint64_t i, j;
1891
1892 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined mode */
1893
1894 /* First sort data on z */
1895 qsort (data, n, sizeof (struct GMT_OBSERVATION), gmtlib_compare_observation);
1896
1897 /* Compute the total weight */
1898 for (wsum = 0.0, i = 0; i < n; i++) wsum += data[i].weight;
1899
1900 /* Do some initializations */
1901 wsum = 0.5 * wsum; /* Sets the 50% range */
1902
1903 /* First check if any single point has 50% or more of the total weights; if so we are done */
1904 for (i = 0; i < n; i++) if (data[i].weight >= wsum) return data[i].value;
1905
1906 /* Some more initializations */
1907 top = p_max = 0.0;
1908 mode = 0.5 * (data[0].value + data[n-1].value);
1909
1910 for (i = j = 0; j < n; j++) {
1911 top += data[j].weight;
1912 if (top < wsum) continue;
1913 while (top > wsum && i < j) top -= data[i++].weight;
1914 bottom = data[j].value - data[i].value;
1915
1916 /* If all is comprised in one point or if a lot of values are the same,
1917 then we have a spike. Maybe another logic is needed to handle
1918 multiple spikes in the data */
1919 if (bottom == 0.0) return (data[i].value);
1920
1921 p = top / bottom;
1922 if (p > p_max) {
1923 p_max = p;
1924 mode = 0.5 * (data[i].value + data[j].value);
1925 }
1926 }
1927 return (mode);
1928 }
1929
gmt_mode(struct GMT_CTRL * GMT,double * x,uint64_t n,uint64_t j,bool sort,int mode_selection,unsigned int * n_multiples,double * mode_est)1930 int gmt_mode (struct GMT_CTRL *GMT, double *x, uint64_t n, uint64_t j, bool sort, int mode_selection, unsigned int *n_multiples, double *mode_est) {
1931 uint64_t i, istop;
1932 unsigned int multiplicity;
1933 double mid_point_sum = 0.0, length, short_length = DBL_MAX, this_mode;
1934
1935 if (n == 0) { /* No data, so no defined mode */
1936 *mode_est = GMT->session.d_NaN;
1937 return (0);
1938 }
1939 if (n == 1) {
1940 *mode_est = x[0];
1941 return (0);
1942 }
1943
1944 if (sort) gmt_sort_array (GMT, x, n, GMT_DOUBLE);
1945 while (n && gmt_M_is_dnan (x[n-1])) n--; /* Skip any NaNs at end of sorted array */
1946 istop = n - j;
1947 multiplicity = 0;
1948
1949 for (i = 0; i < istop; i++) {
1950 length = x[i + j] - x[i];
1951 if (length < 0.0) {
1952 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmt_mode: Array not sorted in non-decreasing order.\n");
1953 return (GMT_NOTSET);
1954 }
1955 else if (length == short_length) { /* Possibly multiple mode */
1956 switch (mode_selection) {
1957 case -1: /* Always pick lowest mode */
1958 this_mode = 0.5 * (x[i + j] + x[i]);
1959 if (this_mode < mid_point_sum) mid_point_sum = this_mode;
1960 break;
1961 case 0: /* Return average of all modes */
1962 multiplicity++;
1963 mid_point_sum += (0.5 * (x[i + j] + x[i]));
1964 break;
1965 case +1: /* Always pick highest mode */
1966 this_mode = 0.5 * (x[i + j] + x[i]);
1967 if (this_mode > mid_point_sum) mid_point_sum = this_mode;
1968 break;
1969 }
1970 }
1971 else if (length < short_length) { /* Update current best mode estimate */
1972 multiplicity = 1;
1973 mid_point_sum = (0.5 * (x[i + j] + x[i]));
1974 short_length = length;
1975 }
1976 }
1977
1978 if (multiplicity > 1) { /* Found more than 1 mode; return mean of them all */
1979 mid_point_sum /= multiplicity;
1980 (*n_multiples) += multiplicity;
1981 }
1982
1983 *mode_est = mid_point_sum;
1984 return (0);
1985 }
1986
gmt_mode_f(struct GMT_CTRL * GMT,gmt_grdfloat * x,uint64_t n,uint64_t j,bool sort,int mode_selection,unsigned int * n_multiples,double * mode_est)1987 int gmt_mode_f (struct GMT_CTRL *GMT, gmt_grdfloat *x, uint64_t n, uint64_t j, bool sort, int mode_selection, unsigned int *n_multiples, double *mode_est) {
1988 uint64_t i, istop;
1989 unsigned int multiplicity;
1990 double mid_point_sum = 0.0, length, short_length = FLT_MAX, this_mode;
1991
1992 if (n == 0) { /* No data, so no defined mode */
1993 *mode_est = GMT->session.d_NaN;
1994 return (0);
1995 }
1996 if (n == 1) {
1997 *mode_est = x[0];
1998 return (0);
1999 }
2000 if (sort) gmt_sort_array (GMT, x, n, GMT_FLOAT);
2001
2002 istop = n - j;
2003 multiplicity = 0;
2004
2005 for (i = 0; i < istop; i++) {
2006 length = x[i + j] - x[i];
2007 if (length < 0.0) {
2008 GMT_Report (GMT->parent, GMT_MSG_WARNING, "gmt_mode_f: Array not sorted in non-decreasing order.\n");
2009 return (GMT_NOTSET);
2010 }
2011 else if (length == short_length) { /* Possibly multiple mode */
2012 switch (mode_selection) {
2013 case -1: /* Always pick lowest mode */
2014 this_mode = 0.5 * (x[i + j] + x[i]);
2015 if (this_mode < mid_point_sum) mid_point_sum = this_mode;
2016 break;
2017 case 0: /* Return average of all modes */
2018 multiplicity++;
2019 mid_point_sum += (0.5 * (x[i + j] + x[i]));
2020 break;
2021 case +1: /* Always pick highest mode */
2022 this_mode = 0.5 * (x[i + j] + x[i]);
2023 if (this_mode > mid_point_sum) mid_point_sum = this_mode;
2024 break;
2025 }
2026 }
2027 else if (length < short_length) {
2028 multiplicity = 1;
2029 mid_point_sum = (0.5 * (x[i + j] + x[i]));
2030 short_length = length;
2031 }
2032 }
2033
2034 if (multiplicity > 1) { /* Found more than 1 mode; return mean of them all */
2035 mid_point_sum /= multiplicity;
2036 (*n_multiples) += multiplicity;
2037 }
2038
2039 *mode_est = mid_point_sum;
2040 return (0);
2041 }
2042
2043 /* Replacement slower functions until we figure out the problem with the algorithm */
2044
gmt_getmad(struct GMT_CTRL * GMT,double * x,uint64_t n,double location,double * scale)2045 void gmt_getmad (struct GMT_CTRL *GMT, double *x, uint64_t n, double location, double *scale) {
2046 uint64_t i;
2047 double med, *dev = NULL;
2048
2049 if (n == 0) { /* No data, so cannot define MAD */
2050 *scale = GMT->session.d_NaN;
2051 return;
2052 }
2053 else if (n == 1) { /* Single point has no deviation */
2054 *scale = 0.0;
2055 return;
2056 }
2057
2058 dev = gmt_M_memory (GMT, NULL, n, double);
2059 for (i = 0; i < n; i++) dev[i] = fabs (x[i] - location);
2060 gmt_sort_array (GMT, dev, n, GMT_DOUBLE);
2061 /* Eliminate any NaNs which would have congregated at the end of the array */
2062 for (i = n; i > 0 && gmt_M_is_dnan (dev[i-1]); i--);
2063 if (i)
2064 med = (i%2) ? dev[i/2] : 0.5 * (dev[(i-1)/2] + dev[i/2]);
2065 else
2066 med = GMT->session.d_NaN;
2067 gmt_M_free (GMT, dev);
2068 *scale = MAD_NORMALIZE * med;
2069 }
2070
gmt_getmad_f(struct GMT_CTRL * GMT,gmt_grdfloat * x,uint64_t n,double location,double * scale)2071 void gmt_getmad_f (struct GMT_CTRL *GMT, gmt_grdfloat *x, uint64_t n, double location, double *scale) {
2072 uint64_t i;
2073 gmt_grdfloat *dev = NULL;
2074 double med;
2075
2076 if (n == 0) { /* No data, so cannot define MAD */
2077 *scale = GMT->session.d_NaN;
2078 return;
2079 }
2080 else if (n == 1) { /* Single point has no deviation */
2081 *scale = 0.0;
2082 return;
2083 }
2084 dev = gmt_M_memory (GMT, NULL, n, double);
2085 for (i = 0; i < n; i++) dev[i] = (gmt_grdfloat) fabs (x[i] - location);
2086 gmt_sort_array (GMT, dev, n, GMT_FLOAT);
2087 for (i = n; i > 0 && gmt_M_is_fnan (dev[i-1]); i--);
2088 if (i)
2089 med = (i%2) ? dev[i/2] : 0.5 * (dev[(i-1)/2] + dev[i/2]);
2090 else
2091 med = GMT->session.d_NaN;
2092 gmt_M_free (GMT, dev);
2093 *scale = MAD_NORMALIZE * med;
2094 }
2095
gmt_extreme(struct GMT_CTRL * GMT,double x[],uint64_t n,double x_default,int kind,int way)2096 double gmt_extreme (struct GMT_CTRL *GMT, double x[], uint64_t n, double x_default, int kind, int way) {
2097 /* Returns the extreme value in the x array according to:
2098 * kind: -1 means only consider negative values, including zero.
2099 * kind: 0 means consider all values.
2100 * kind: +1 means only consider positive values, including zero.
2101 * way: -1 means look for minimum.
2102 * way: +1 means look for maximum.
2103 * If kind is non-zero we assign x_default is no values are found.
2104 * If kind == 0 we expect x_default to be set so that x[0] will reset x_select.
2105 */
2106
2107 uint64_t i, k;
2108 double x_select = GMT->session.d_NaN;
2109
2110 if (n == 0) return (x_select); /* No data, so no defined extreme value */
2111 for (i = k = 0; i < n; i++) {
2112 if (kind == -1 && x[i] > 0.0) continue;
2113 if (kind == +1 && x[i] < 0.0) continue;
2114 if (k == 0) x_select = x[i];
2115 if (way == -1 && x[i] < x_select) x_select = x[i];
2116 if (way == +1 && x[i] > x_select) x_select = x[i];
2117 k++;
2118 }
2119
2120 return ((k) ? x_select : x_default);
2121 }
2122
gmt_chebyshev(struct GMT_CTRL * GMT,double x,int n,double * t)2123 int gmt_chebyshev (struct GMT_CTRL *GMT, double x, int n, double *t) {
2124 /* Calculates the n'th Chebyshev polynomial at x */
2125
2126 double x2, a, b;
2127
2128 if (n < 0) gmt_M_err_pass (GMT, GMT_CHEBYSHEV_NEG_ORDER, "");
2129 if (fabs (x) > 1.0) gmt_M_err_pass (GMT, GMT_CHEBYSHEV_BAD_DOMAIN, "");
2130
2131 switch (n) { /* Testing the order of the polynomial */
2132 case 0:
2133 *t = 1.0;
2134 break;
2135 case 1:
2136 *t = x;
2137 break;
2138 case 2:
2139 *t = 2.0 * x * x - 1.0;
2140 break;
2141 case 3:
2142 *t = x * (4.0 * x * x - 3.0);
2143 break;
2144 case 4:
2145 x2 = x * x;
2146 *t = 8.0 * x2 * (x2 - 1.0) + 1.0;
2147 break;
2148 default: /* For higher degrees we do the recursion */
2149 gmt_chebyshev (GMT, x, n-1, &a);
2150 gmt_chebyshev (GMT, x, n-2, &b);
2151 *t = 2.0 * x * a - b;
2152 break;
2153 }
2154
2155 return (GMT_NOERROR);
2156 }
2157
gmt_corrcoeff(struct GMT_CTRL * GMT,double * x,double * y,uint64_t n,unsigned int mode)2158 double gmt_corrcoeff (struct GMT_CTRL *GMT, double *x, double *y, uint64_t n, unsigned int mode) {
2159 /* Returns plain correlation coefficient, r.
2160 * If mode = 1 we assume mean(x) = mean(y) = 0.
2161 */
2162
2163 uint64_t i, n_use;
2164 double xmean = 0.0, ymean = 0.0, dx, dy, vx, vy, vxy, r;
2165
2166 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined correlation */
2167 if (mode == 0) {
2168 for (i = n_use = 0; i < n; i++) {
2169 if (gmt_M_is_dnan (x[i]) || gmt_M_is_dnan (y[i])) continue;
2170 xmean += x[i];
2171 ymean += y[i];
2172 n_use++;
2173 }
2174 if (n_use == 0) return (GMT->session.d_NaN);
2175 xmean /= (double)n_use;
2176 ymean /= (double)n_use;
2177 }
2178
2179 vx = vy = vxy = 0.0;
2180 for (i = n_use = 0; i < n; i++) {
2181 if (gmt_M_is_dnan (x[i]) || gmt_M_is_dnan (y[i])) continue;
2182 dx = x[i] - xmean;
2183 dy = y[i] - ymean;
2184 vx += dx * dx;
2185 vy += dy * dy;
2186 vxy += dx * dy;
2187 }
2188
2189 r = vxy / sqrt (vx * vy);
2190 return (r);
2191 }
2192
gmt_corrcoeff_f(struct GMT_CTRL * GMT,gmt_grdfloat * x,gmt_grdfloat * y,uint64_t n,unsigned int mode)2193 double gmt_corrcoeff_f (struct GMT_CTRL *GMT, gmt_grdfloat *x, gmt_grdfloat *y, uint64_t n, unsigned int mode) {
2194 /* Returns plain correlation coefficient, r.
2195 * If mode = 1 we assume mean(x) = mean(y) = 0.
2196 */
2197
2198 uint64_t i, n_use;
2199 double xmean = 0.0, ymean = 0.0, dx, dy, vx, vy, vxy, r;
2200
2201 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined correlation */
2202 if (mode == 0) {
2203 for (i = n_use = 0; i < n; i++) {
2204 if (gmt_M_is_fnan (x[i]) || gmt_M_is_fnan (y[i])) continue;
2205 xmean += (double)x[i];
2206 ymean += (double)y[i];
2207 n_use++;
2208 }
2209 if (n_use == 0) return (GMT->session.d_NaN);
2210 xmean /= (double)n_use;
2211 ymean /= (double)n_use;
2212 }
2213
2214 vx = vy = vxy = 0.0;
2215 for (i = n_use = 0; i < n; i++) {
2216 if (gmt_M_is_fnan (x[i]) || gmt_M_is_fnan (y[i])) continue;
2217 dx = (double)x[i] - xmean;
2218 dy = (double)y[i] - ymean;
2219 vx += dx * dx;
2220 vy += dy * dy;
2221 vxy += dx * dy;
2222 }
2223
2224 r = vxy / sqrt (vx * vy);
2225 return (r);
2226 }
2227
gmt_quantile(struct GMT_CTRL * GMT,double * x,double q,uint64_t n)2228 double gmt_quantile (struct GMT_CTRL *GMT, double *x, double q, uint64_t n) {
2229 /* Returns the q'th (q in percent) quantile of x (assumed sorted).
2230 * q is expected to be 0 < q < 100 */
2231
2232 uint64_t i_f;
2233 double p, f, df;
2234
2235 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined quantile */
2236 if (q == 0.0) return (x[0]); /* 0% quantile == min(x) */
2237 while (n > 1 && gmt_M_is_dnan (x[n-1])) n--; /* Skip any NaNs at the end of x */
2238 if (q == 100.0) return (x[n-1]); /* 100% quantile == max(x) */
2239 f = (n - 1) * q / 100.0;
2240 i_f = (uint64_t)floor (f);
2241 if ((df = (f - (double)i_f)) > 0.0) /* Must interpolate between the two neighbors */
2242 p = x[i_f+1] * df + x[i_f] * (1.0 - df);
2243 else /* Exactly on a node */
2244 p = x[i_f];
2245
2246 return (p);
2247 }
2248
gmt_quantile_f(struct GMT_CTRL * GMT,gmt_grdfloat * x,double q,uint64_t n)2249 double gmt_quantile_f (struct GMT_CTRL *GMT, gmt_grdfloat *x, double q, uint64_t n) {
2250 /* Returns the q'th (q in percent) quantile of x (assumed sorted).
2251 * q is expected to be 0 < q < 100 */
2252
2253 uint64_t i_f;
2254 double p, f, df;
2255
2256 if (n == 0) return (GMT->session.d_NaN); /* No data, so no defined quantile */
2257 if (q == 0.0) return ((double)x[0]); /* 0% quantile == min(x) */
2258 while (n > 1 && gmt_M_is_fnan (x[n-1])) n--; /* Skip any NaNs at the end of x */
2259 if (q == 100.0) return ((double)x[n-1]); /* 100% quantile == max(x) */
2260 f = (n - 1) * q / 100.0;
2261 i_f = (uint64_t)floor (f);
2262 if ((df = (f - (double)i_f)) > 0.0) /* Must interpolate between the two neighbors */
2263 p = (double)(x[i_f+1] * df + x[i_f] * (1.0 - df));
2264 else /* Exactly on a node */
2265 p = (double)x[i_f];
2266
2267 return (p);
2268 }
2269
gmt_psi(struct GMT_CTRL * P,double zz[],double p[])2270 double gmt_psi (struct GMT_CTRL *P, double zz[], double p[]) {
2271 /* Psi Psi (or Digamma) function for complex arguments z.
2272 *
2273 * d
2274 * Psi(z) = --log(Gamma(z))
2275 * dz
2276 *
2277 * zz[GMT_RE] is real and zz[GMT_IM] is imaginary component; same for p on output (only if p != NULL).
2278 * We also return the real component as function result.
2279 */
2280 static double c[15] = { 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248,
2281 14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699e-4,
2282 0.46523628927048575665e-4, -0.98374475304879564677e-4, 0.15808870322491248884e-3,
2283 -0.21026444172410488319e-3, 0.21743961811521264320e-3, -0.16431810653676389022e-3,
2284 0.84418223983852743293e-4, -0.26190838401581408670e-4, 0.36899182659531622704e-5};
2285 double z[2], g[2], dx[2], dd[2], d[2], n[2], gg[2], f[2], x0, A[2], B[2], C[2], sx, cx, e;
2286 int k;
2287
2288 if (zz[GMT_IM] == 0.0 && rint(zz[GMT_RE]) == zz[GMT_RE] && zz[GMT_RE] <= 0.0) {
2289 if (p) { p[GMT_RE] = P->session.d_NaN; p[GMT_IM] = 0.0;}
2290 return (P->session.d_NaN); /* Singular points */
2291 }
2292
2293 z[GMT_RE] = zz[GMT_RE]; z[GMT_IM] = zz[GMT_IM];
2294 if ((x0 = z[GMT_RE]) < 0.5) { /* reflection point */
2295 z[GMT_RE] = 1.0 - z[GMT_RE];
2296 z[GMT_IM] = -z[GMT_IM];
2297 }
2298
2299 /* Lanczos approximation */
2300
2301 g[GMT_RE] = 607.0/128.0; g[GMT_IM] = 0.0; /* best results when 4<=g<=5 */
2302 n[GMT_RE] = d[GMT_RE] = n[GMT_IM] = d[GMT_IM] = 0.0;
2303 for (k = 14; k > 0; k--) {
2304 A[GMT_RE] = 1.0; A[GMT_IM] = 0.0;
2305 B[GMT_RE] = z[GMT_RE] + k - 1.0; B[GMT_IM] = z[GMT_IM];
2306 gmtstat_Cdiv (A, B, dx);
2307 dd[GMT_RE] = c[k] * dx[GMT_RE]; dd[GMT_IM] = c[k] * dx[GMT_IM];
2308 d[GMT_RE] += dd[GMT_RE]; d[GMT_IM] += dd[GMT_IM];
2309 gmtstat_Cmul (dd, dx, B);
2310 n[GMT_RE] -= B[GMT_RE]; n[GMT_IM] -= B[GMT_IM];
2311 }
2312 d[GMT_RE] += c[GMT_RE];
2313 gg[GMT_RE] = z[GMT_RE] + g[GMT_RE] - 0.5; gg[GMT_IM] = z[GMT_IM];
2314 gmtstat_Cdiv (n, d, A);
2315 gmtstat_Cdiv (g, gg, B);
2316 f[GMT_RE] = log (hypot(gg[GMT_RE], gg[GMT_IM])) + A[GMT_RE] - B[GMT_RE];
2317 f[GMT_IM] = atan2 (gg[GMT_IM], gg[GMT_RE]) + A[GMT_IM] - B[GMT_IM];
2318 if (x0 < 0.5) {
2319 C[GMT_RE] = M_PI * zz[GMT_RE]; C[GMT_IM] = M_PI * zz[GMT_IM];
2320 e = exp (-2*C[GMT_IM]); sx = sin (2*C[GMT_RE]); cx = cos (2*C[GMT_RE]);
2321 A[GMT_RE] = -e * sx; A[GMT_IM] = e * cx + 1.0;
2322 B[GMT_RE] = e * cx - 1.0; B[GMT_IM] = e * sx;
2323 gmtstat_Cdiv (A, B, C);
2324 f[GMT_RE] -= M_PI * C[GMT_RE]; f[GMT_IM] -= M_PI * C[GMT_IM];
2325 }
2326 if (p) {
2327 p[GMT_RE] = f[GMT_RE];
2328 p[GMT_IM] = f[GMT_IM];
2329 }
2330 return (f[GMT_RE]);
2331 }
2332
2333 #ifndef M_SQRT_PI
2334 #define M_SQRT_PI 1.772453850905516
2335 #endif
2336
2337 #define PV_RE 0
2338 #define PV_IM 1
2339 #define QV_RE 2
2340 #define QV_IM 3
2341
gmt_PvQv(struct GMT_CTRL * GMT,double x,double v_ri[],double pq[],unsigned int * iter)2342 void gmt_PvQv (struct GMT_CTRL *GMT, double x, double v_ri[], double pq[], unsigned int *iter) {
2343 /* Here, -1 <= x <= +1, v_ri is an imaginary number [r,i], and we return
2344 * the real amd imaginary parts of Pv(x) and Qv(x) in the pq array.
2345 * Based on recipe in An Atlas of Functions */
2346
2347 bool p_set, q_set;
2348 double M, L, K, Xn, x2, k, k1, ep, em, sx, cx, fact;
2349 double a[2], v[2], vp1[2], G[2], g[2], u[2], t[2], f[2];
2350 double R[2], r[2], z[2], s[2], c[2], w[2], tmp[2], X[2], A[2], B[2];
2351
2352 *iter = 0;
2353 pq[PV_RE] = pq[QV_RE] = pq[PV_IM] = pq[QV_IM] = 0.0; /* Initialize all answers to zero first */
2354 p_set = q_set = false;
2355 if (x == -1 && v_ri[GMT_IM] == 0.0) { /* Check special values for real nu when x = -1 */
2356 /* Special Pv(-1) values */
2357 if ((v_ri[GMT_RE] > 0.0 && fmod (v_ri[GMT_RE], 2.0) == 1.0) || (v_ri[GMT_RE] < 0.0 && fmod (v_ri[GMT_RE], 2.0) == 0.0)) { /* v = 1,3,5,.. or -2, -4, -6 */
2358 pq[PV_RE] = -1.0; p_set = true;
2359 }
2360 else if ((v_ri[GMT_RE] >= 0.0 && fmod (v_ri[GMT_RE], 2.0) == 0.0) || (v_ri[GMT_RE] < 0.0 && fmod (v_ri[GMT_RE]+1.0, 2.0) == 0.0)) { /* v = 0,2,4,6 or -1,-3,-5,-7 */
2361 pq[PV_RE] = 1.0; p_set = true;
2362 }
2363 else if (v_ri[GMT_RE] > 0.0 && ((fact = v_ri[GMT_RE]-2.0*floor(v_ri[GMT_RE]/2.0)) < 1.0 && fact > 0.0)) { /* 0 < v < 1, 2 < v < 3, 4 < v < 5, ... */
2364 pq[PV_RE] = GMT->session.d_NaN; p_set = true; /* -inf */
2365 }
2366 else if (v_ri[GMT_RE] < 1.0 && ((fact = v_ri[GMT_RE]-2.0*floor(v_ri[GMT_RE]/2.0)) < 1.0 && fact > 0.0)) { /* -2 < v < -1, -4 < v < -3, -6 < v < -5 */
2367 pq[PV_RE] = GMT->session.d_NaN; p_set = true; /* -inf */
2368 }
2369 else if (v_ri[GMT_RE] > 1.0 && (v_ri[GMT_RE]-2.0*floor(v_ri[GMT_RE]/2.0)) > 1.0) { /* 1 < v < 2, 3 < v < 4, 5 < v < 6, .. */
2370 pq[PV_RE] = GMT->session.d_NaN; p_set = true; /* +inf */
2371 }
2372 else if (v_ri[GMT_RE] < 2.0 && ( v_ri[GMT_RE]-2.0*floor(v_ri[GMT_RE]/2.0)) > 1.0) { /* -3 < v < -2, -5 < v < -4, -7 < v < -6, .. */
2373 pq[PV_RE] = GMT->session.d_NaN; p_set = true; /* +inf */
2374 }
2375 /* Special Qv(-1) values */
2376 if (v_ri[GMT_RE] > 0.0 && fmod (2.0 * v_ri[GMT_RE] - 1.0, 4.0) == 0.0) { /* v = 1/2, 5/2, 9/2, ... */
2377 pq[QV_RE] = -M_PI_2; q_set = true;
2378 }
2379 else if (v_ri[GMT_RE] > -1.0 && fmod (2.0 * v_ri[GMT_RE] + 1.0, 4.0) == 0.0) { /* v = -1/2, 3/2, 7/2, ... */
2380 pq[QV_RE] = M_PI_2; q_set = true;
2381 }
2382 else if (v_ri[GMT_RE] < -1.0 && fmod (2.0 * v_ri[GMT_RE] - 1.0, 4.0) == 0.0) { /* v = -3/2, -7/2, -11/2, ... */
2383 pq[QV_RE] = -M_PI_2; q_set = true;
2384 }
2385 else if (v_ri[GMT_RE] < -2.0 && fmod (2.0 * v_ri[GMT_RE] + 1.0, 4.0) == 0.0) { /* v = -5/2, -9/2, -13/2, ... */
2386 pq[QV_RE] = M_PI_2; q_set = true;
2387 }
2388 else { /* Either -inf or +inf */
2389 pq[QV_RE] = GMT->session.d_NaN; q_set = true;
2390 }
2391 }
2392 else if (x == +1 && v_ri[GMT_IM] == 0.0) { /* Check special values for real nu when x = +1 */
2393 pq[PV_RE] = 1.0; p_set = true;
2394 pq[QV_RE] = GMT->session.d_NaN; q_set = true;
2395 }
2396 if (p_set && q_set) return;
2397
2398 /* General case of |x| < 1 */
2399
2400 a[0] = a[1] = R[GMT_RE] = 1.0; R[GMT_IM] = 0.0;
2401 v[GMT_RE] = v_ri[GMT_RE]; v[GMT_IM] = v_ri[GMT_IM];
2402 gmtstat_Cmul (v, v, z);
2403 z[GMT_RE] = v[GMT_RE] - z[GMT_RE]; z[GMT_IM] = v[GMT_IM] - z[GMT_IM];
2404 K = 4.0 * sqrt (gmtstat_Cabs(z));
2405 vp1[GMT_RE] = v[GMT_RE] + 1.0; vp1[GMT_IM] = v[GMT_IM];
2406 if ((gmtstat_Cabs(vp1) + floor(vp1[GMT_RE])) == 0.0) {
2407 a[0] = GMT->session.d_NaN;
2408 a[1] = 0.0;
2409 v[GMT_RE] = -1 - v[GMT_RE];
2410 v[GMT_IM] = -v[GMT_IM];
2411 }
2412 z[GMT_RE] = 0.5 * M_PI * v[GMT_RE]; z[GMT_IM] = 0.5 * M_PI * v[GMT_IM];
2413 ep = exp (z[GMT_IM]); em = exp (-z[GMT_IM]);
2414 sincos (z[GMT_RE], &sx, &cx);
2415 s[GMT_RE] = 0.5 * sx * (em + ep);
2416 s[GMT_IM] = -0.5 * cx * (em - ep);
2417 c[GMT_RE] = 0.5 * cx * (em + ep);
2418 c[GMT_IM] = 0.5 * sx * (em - ep);
2419 tmp[GMT_RE] = 0.5 + v[GMT_RE]; tmp[GMT_IM] = v[GMT_IM];
2420 gmtstat_Cmul (tmp, tmp, w);
2421 z[GMT_IM] = v[GMT_IM];
2422 while (v[GMT_RE] <= 6.0) {
2423 v[GMT_RE] = v[GMT_RE] + 2.0;
2424 z[GMT_RE] = v[GMT_RE] - 1.0;
2425 gmtstat_Cdiv (z, v, tmp);
2426 gmtstat_Cmul (R,tmp,r);
2427 R[GMT_RE] = r[GMT_RE]; R[GMT_IM] = r[GMT_IM];
2428 }
2429 z[GMT_RE] = v[GMT_RE] + 1.0;
2430 tmp[GMT_RE] = 0.25; tmp[GMT_IM] = 0.0;
2431 gmtstat_Cdiv (tmp, z, X);
2432 tmp[GMT_RE] = 0.35 + 6.1 * X[GMT_RE]; tmp[GMT_IM] = 6.1*X[GMT_IM];
2433 gmtstat_Cmul (X, tmp, z);
2434 z[GMT_RE] = 1.0 - 3.0*z[GMT_RE]; z[GMT_IM] = -3.0*z[GMT_IM];
2435 gmtstat_Cmul (X, z, tmp);
2436 G[GMT_RE] = 1.0 + 5.0 * tmp[GMT_RE]; G[GMT_IM] = 5.0 * tmp[GMT_IM];
2437 z[GMT_RE] = 8.0 * X[GMT_RE]; z[GMT_IM] = 8.0 * X[GMT_IM];
2438 M = sqrt(hypot(z[GMT_RE], z[GMT_IM]));
2439 L = 0.5 * atan2 (z[GMT_IM], z[GMT_RE]);
2440 tmp[GMT_RE] = M * cos(L); tmp[GMT_IM] = M * sin(L);
2441 gmtstat_Cmul (G, X, z);
2442 z[GMT_RE] = 1.0 - 0.5*z[GMT_RE]; z[GMT_IM] = -0.5*z[GMT_IM];
2443 gmtstat_Cmul (X, z, r);
2444 r[GMT_RE] = 1.0 - r[GMT_RE]; r[GMT_IM] = -r[GMT_IM];
2445 gmtstat_Cmul (R, r, z);
2446 gmtstat_Cdiv (z, tmp, R);
2447 u[GMT_RE] = g[GMT_RE] = 2.0 * x; u[GMT_IM] = g[GMT_IM] = f[GMT_IM] = t[GMT_IM] = 0.0;
2448 f[GMT_RE] = t[GMT_RE] = 1.0;
2449 k = 0.5;
2450 x2 = x * x;
2451 Xn = 1.0 + (1e8/(1 - x2));
2452 k1 = k + 1.0;
2453 fact = x2 / (k1*k1 - 0.25);
2454 t[GMT_RE] = (k*k - w[GMT_RE]) * fact; t[GMT_IM] = -w[GMT_IM] * fact;
2455 k += 1.0;
2456 f[GMT_RE] += t[GMT_RE]; f[GMT_IM] += t[GMT_IM];
2457 k1 = k + 1.0;
2458 fact = u[GMT_RE] * x2 / (k1*k1 - 0.25);
2459 u[GMT_RE] = (k*k - w[GMT_RE]) * fact; u[GMT_IM] = -w[GMT_IM] * fact;
2460 k += 1;
2461 g[GMT_RE] += u[GMT_RE]; g[GMT_IM] += u[GMT_IM];
2462 tmp[GMT_RE] = Xn * t[GMT_RE]; tmp[GMT_IM] = Xn * t[GMT_IM];
2463 while (k < K || gmtstat_Cabs (tmp) > gmtstat_Cabs(f)) {
2464 (*iter)++;
2465 k1 = k + 1.0;
2466 tmp[GMT_RE] = k*k - w[GMT_RE]; tmp[GMT_IM] = -w[GMT_IM]; fact = x2 / (k1*k1 - 0.25);
2467 gmtstat_Cmul (t, tmp, A);
2468 t[GMT_RE] = A[GMT_RE] * fact; t[GMT_IM] = A[GMT_IM] * fact;
2469 k += 1.0;
2470 k1 = k + 1.0;
2471 f[GMT_RE] += t[GMT_RE]; f[GMT_IM] += t[GMT_IM];
2472 tmp[GMT_RE] = k*k - w[GMT_RE]; tmp[GMT_IM] = -w[GMT_IM]; fact = x2 / (k1*k1 - 0.25);
2473 gmtstat_Cmul (u, tmp, B);
2474 u[GMT_RE] = B[GMT_RE] * fact; u[GMT_IM] = B[GMT_IM] * fact;
2475 k += 1.0;
2476 g[GMT_RE] += u[GMT_RE]; g[GMT_IM] += u[GMT_IM];
2477 tmp[GMT_RE] = Xn * t[GMT_RE]; tmp[GMT_IM] = Xn * t[GMT_IM];
2478 }
2479 fact = x2 / (1.0 - x2);
2480 f[GMT_RE] += t[GMT_RE] * fact; f[GMT_IM] += t[GMT_IM] * fact;
2481 g[GMT_RE] += u[GMT_RE] * fact; g[GMT_IM] += u[GMT_IM] * fact;
2482 if (!p_set) {
2483 gmtstat_Cmul(s,R,z);
2484 gmtstat_Cdiv(c,R,tmp);
2485 gmtstat_Cmul (g, z, A); gmtstat_Cmul (f, tmp, B);
2486 pq[PV_RE] = (A[GMT_RE] + B[GMT_RE])/M_SQRT_PI;
2487 pq[PV_IM] = (A[GMT_IM] + B[GMT_IM])/M_SQRT_PI;
2488 }
2489 if (!q_set) {
2490 gmtstat_Cmul(c,R,z);
2491 gmtstat_Cdiv(s,R,tmp);
2492 gmtstat_Cmul (g, z, A); gmtstat_Cmul (f, tmp, B);
2493 pq[QV_RE] = a[0]*M_SQRT_PI*(A[GMT_RE] - B[GMT_RE])/2.0;
2494 pq[QV_IM] = a[1]*M_SQRT_PI*(A[GMT_IM] - B[GMT_IM])/2.0;
2495 }
2496 }
2497
gmt_grd_mean(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W)2498 double gmt_grd_mean (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W) {
2499 /* Compute the [weighted] mean of a grid. Handle geographic grids with spherical weights W [NULL for cartesian] */
2500 uint64_t node, n = 0;
2501 unsigned int row, col;
2502 double sum_zw = 0.0, sum_w = 0.0;
2503 if (W) { /* Weights provided */
2504 gmt_M_grd_loop (GMT, G, row, col, node) {
2505 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node])) continue;
2506 sum_zw += G->data[node] * W->data[node];
2507 sum_w += W->data[node];
2508 n++;
2509 }
2510 }
2511 else { /* Plain average */
2512 gmt_M_grd_loop (GMT, G, row, col, node) {
2513 if (gmt_M_is_fnan (G->data[node])) continue;
2514 sum_zw += G->data[node];
2515 n++;
2516 }
2517 sum_w = (double)n;
2518 }
2519 return (n == 0 || sum_w == 0.0) ? GMT->session.d_NaN : sum_zw / sum_w;
2520 }
2521
gmt_grd_std(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W)2522 double gmt_grd_std (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W) {
2523 /* Compute the [weighted] std of a grid. Handle geographic grids with spherical weights W [NULL for cartesian] */
2524 uint64_t node, n = 0;
2525 unsigned int row, col;
2526 double std, mean = 0.0, delta, sumw = 0.0;
2527 if (W) { /* Weights provided */
2528 double temp, R, M2 = 0.0;
2529 gmt_M_grd_loop (GMT, G, row, col, node) {
2530 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node])) continue;
2531 temp = W->data[node] + sumw;
2532 delta = G->data[node] - mean;
2533 R = delta * W->data[node] / temp;
2534 mean += R;
2535 M2 += sumw * delta * R;
2536 sumw = temp;
2537 n++;
2538 }
2539 std = (n <= 1 || sumw == 0.0) ? GMT->session.d_NaN : sqrt ((n * M2) / (sumw * (n - 1.0)));
2540 }
2541 else { /* Plain average */
2542 gmt_M_grd_loop (GMT, G, row, col, node) {
2543 if (gmt_M_is_fnan (G->data[node])) continue;
2544 n++;
2545 delta = G->data[node] - mean;
2546 mean += delta / n;
2547 sumw += delta * (G->data[node] - mean);
2548 }
2549 std = (n > 1) ? sqrt (sumw / (n - 1.0)) : GMT->session.d_NaN;
2550 }
2551 return (std);
2552 }
2553
gmt_grd_rms(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W)2554 double gmt_grd_rms (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W) {
2555 /* Compute the [weighted] rms of a grid. Handle geographic grids with spherical weights W [NULL for cartesian] */
2556 uint64_t node, n = 0;
2557 unsigned int row, col;
2558 double rms, sum_z2w = 0.0, sum_w = 0.0;
2559 if (W) { /* Weights provided */
2560 gmt_M_grd_loop (GMT, G, row, col, node) {
2561 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node])) continue;
2562 sum_z2w += W->data[node] * (G->data[node] * G->data[node]);
2563 sum_w += W->data[node];
2564 }
2565 }
2566 else { /* Plain average */
2567 gmt_M_grd_loop (GMT, G, row, col, node) {
2568 if (gmt_M_is_fnan (G->data[node])) continue;
2569 n++;
2570 sum_z2w += (G->data[node] * G->data[node]);
2571 }
2572 sum_w = (double)n;
2573 }
2574 rms = (sum_w > 0) ? sqrt (sum_z2w / sum_w) : GMT->session.d_NaN;
2575 return (rms);
2576 }
2577
gmt_grd_median(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W,bool overwrite)2578 double gmt_grd_median (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W, bool overwrite) {
2579 /* Compute the weighted median of a grid. Handle geographic grids with spherical weights W */
2580 /* Non-destructive: Original grid left as is unless overwrite = true */
2581 uint64_t node, n = 0;
2582 double wmed;
2583
2584 if (W) { /* Weights provided */
2585 unsigned int row, col;
2586 struct GMT_OBSERVATION *pair = gmt_M_memory (GMT, NULL, G->header->nm, struct GMT_OBSERVATION);
2587 /* 1. Create array of value,weight pairs, skipping NaNs */
2588 gmt_M_grd_loop (GMT, G, row, col, node) {
2589 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2590 continue;
2591 pair[n].value = G->data[node];
2592 pair[n++].weight = W->data[node];
2593 }
2594 /* 2. Find the weighted median */
2595 wmed = gmt_median_weighted (GMT, pair, n);
2596 gmt_M_free (GMT, pair);
2597 }
2598 else { /* Plain median */
2599 struct GMT_GRID *Z = (overwrite) ? G : gmt_duplicate_grid (GMT, G, GMT_DUPLICATE_DATA);
2600 gmt_grd_pad_off (GMT, Z); /* Undo pad if one existed so we can sort */
2601 gmt_sort_array (GMT, Z->data, Z->header->nm, GMT_FLOAT);
2602 for (n = Z->header->nm; n > 1 && gmt_M_is_fnan (Z->data[n-1]); n--);
2603 if (n)
2604 wmed = (n%2) ? Z->data[n/2] : 0.5f * (Z->data[(n-1)/2] + Z->data[n/2]);
2605 else
2606 wmed = GMT->session.f_NaN;
2607 if (!overwrite) gmt_free_grid (GMT, &Z, true);
2608 }
2609 return wmed;
2610 }
2611
gmt_grd_mad(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W,double * median,bool overwrite)2612 double gmt_grd_mad (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W, double *median, bool overwrite) {
2613 /* Compute the [weighted] MAD of a grid. Handle geographic grids with spherical weights W [NULL for cartesian].
2614 * Non-destructive: Original grid left as is unless overwrite = true.
2615 * If median == NULL then we must first compute the median, else we already have it */
2616 uint64_t node, n = 0;
2617 double wmed, wmad;
2618 if (W) { /* Weights provided */
2619 unsigned int row, col;
2620 struct GMT_OBSERVATION *pair = gmt_M_memory (GMT, NULL, G->header->nm, struct GMT_OBSERVATION);
2621 if (median) { /* Already have the median */
2622 wmed = *median;
2623 /* 3. Compute the absolute deviations from this median */
2624 gmt_M_grd_loop (GMT, G, row, col, node) {
2625 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2626 continue;
2627 pair[n].value = (gmt_grdfloat)fabs (G->data[node] - wmed);
2628 pair[n++].weight = W->data[node];
2629 }
2630 }
2631 else { /* Must first compute median */
2632 /* 1. Create array of value,weight pairs, skipping NaNs */
2633 gmt_M_grd_loop (GMT, G, row, col, node) {
2634 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2635 continue;
2636 pair[n].value = G->data[node];
2637 pair[n++].weight = W->data[node];
2638 }
2639 /* 2. Find the weighted median */
2640 wmed = gmt_median_weighted (GMT, pair, n);
2641 /* 3. Compute the absolute deviations from this median */
2642 for (node = 0; node < n; node++) pair[node].value = (gmt_grdfloat)fabs (pair[node].value - wmed);
2643 }
2644 /* 4. Find the weighted median absolute deviation */
2645 wmad = MAD_NORMALIZE * gmt_median_weighted (GMT, pair, n);
2646 gmt_M_free (GMT, pair);
2647 }
2648 else { /* Plain MAD */
2649 struct GMT_GRID *Z = (overwrite) ? G : gmt_duplicate_grid (GMT, G, GMT_DUPLICATE_DATA);
2650 gmt_grd_pad_off (GMT, Z); /* Undo pad if one existed so we can sort */
2651 if (median) { /* Already have the median */
2652 wmed = *median;
2653 n = Z->header->nm;
2654 }
2655 else { /* Must first compute the median */
2656 gmt_sort_array (GMT, Z->data, Z->header->nm, GMT_FLOAT);
2657 for (n = Z->header->nm; n > 1 && gmt_M_is_fnan (Z->data[n-1]); n--);
2658 if (n)
2659 wmed = (n%2) ? Z->data[n/2] : 0.5 * (Z->data[(n-1)/2] + Z->data[n/2]);
2660 else
2661 wmed = GMT->session.d_NaN;
2662 }
2663 gmt_getmad_f (GMT, Z->data, n, wmed, &wmad);
2664 if (!overwrite) gmt_free_grid (GMT, &Z, true);
2665 }
2666 return wmad;
2667 }
2668
gmt_grd_mode(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W,bool overwrite)2669 double gmt_grd_mode (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W, bool overwrite) {
2670 /* Compute the weighted mode (lms) of a grid. Handle geographic grids with spherical weights W */
2671 /* Non-destructive: Original grid left as is unless overwrite = true */
2672 uint64_t node, n = 0;
2673 double wmode;
2674
2675 if (W) { /* Weights provided */
2676 unsigned int row, col;
2677 struct GMT_OBSERVATION *pair = gmt_M_memory (GMT, NULL, G->header->nm, struct GMT_OBSERVATION);
2678 /* 1. Create array of value,weight pairs, skipping NaNs */
2679 gmt_M_grd_loop (GMT, G, row, col, node) {
2680 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2681 continue;
2682 pair[n].value = G->data[node];
2683 pair[n++].weight = W->data[node];
2684 }
2685 /* 2. Find the weighted lms mode */
2686 wmode = (gmt_grdfloat)gmt_mode_weighted (GMT, pair, n);
2687 gmt_M_free (GMT, pair);
2688 }
2689 else { /* Plain median */
2690 unsigned int gmt_mode_selection = 0, GMT_n_multiples = 0;
2691 struct GMT_GRID *Z = (overwrite) ? G : gmt_duplicate_grid (GMT, G, GMT_DUPLICATE_DATA);
2692 gmt_grd_pad_off (GMT, Z); /* Undo pad if one existed so we can sort */
2693 gmt_sort_array (GMT, Z->data, Z->header->nm, GMT_FLOAT);
2694 for (n = Z->header->nm; n > 1 && gmt_M_is_fnan (Z->data[n-1]); n--);
2695 if (n)
2696 gmt_mode_f (GMT, Z->data, n, n/2, 0, gmt_mode_selection, &GMT_n_multiples, &wmode);
2697 else
2698 wmode = GMT->session.d_NaN;
2699 if (!overwrite) gmt_free_grid (GMT, &Z, true);
2700 if (GMT_n_multiples > 0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "%d Multiple modes found in the grid\n", GMT_n_multiples);
2701 }
2702 return wmode;
2703 }
2704
gmt_grd_lmsscl(struct GMT_CTRL * GMT,struct GMT_GRID * G,struct GMT_GRID * W,double * mode,bool overwrite)2705 double gmt_grd_lmsscl (struct GMT_CTRL *GMT, struct GMT_GRID *G, struct GMT_GRID *W, double *mode, bool overwrite) {
2706 /* Compute the [weighted] LMSSCL of a grid. Handle geographic grids with spherical weights W [NULL for cartesian]
2707 * Non-destructive: Original grid left as is unless overwrite = true.
2708 * If mode == NULL then we must first compute the mode, else use this mode */
2709 uint64_t node, n = 0;
2710 double wmode, lmsscl;
2711 if (W) { /* Weights provided */
2712 unsigned int row, col;
2713 struct GMT_OBSERVATION *pair = gmt_M_memory (GMT, NULL, G->header->nm, struct GMT_OBSERVATION);
2714 if (mode) { /* Already got the mode */
2715 wmode = *mode;
2716 /* 3. Compute the absolute deviations from this mode */
2717 gmt_M_grd_loop (GMT, G, row, col, node) {
2718 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2719 continue;
2720 pair[n].value = (gmt_grdfloat)fabs (G->data[node] - wmode);
2721 pair[n++].weight = W->data[node];
2722 }
2723 }
2724 else { /* Must first find the mode */
2725 /* 1. Create array of value,weight pairs, skipping NaNs */
2726 gmt_M_grd_loop (GMT, G, row, col, node) {
2727 if (gmt_M_is_fnan (G->data[node]) || gmt_M_is_dnan (W->data[node]))
2728 continue;
2729 pair[n].value = G->data[node];
2730 pair[n++].weight = W->data[node];
2731 }
2732 /* 2. Find the weighted median */
2733 wmode = (gmt_grdfloat)gmt_mode_weighted (GMT, pair, n);
2734 /* 3. Compute the absolute deviations from this mode */
2735 for (node = 0; node < n; node++) pair[node].value = (gmt_grdfloat)fabs (pair[node].value - wmode);
2736 }
2737 /* 4. Find the weighted median absolute deviation and scale it */
2738 lmsscl = MAD_NORMALIZE * gmt_median_weighted (GMT, pair, n);
2739 gmt_M_free (GMT, pair);
2740 }
2741 else { /* Plain LMSSCL */
2742 unsigned int gmt_mode_selection = 0, GMT_n_multiples = 0;
2743 struct GMT_GRID *Z = (overwrite) ? G : gmt_duplicate_grid (GMT, G, GMT_DUPLICATE_DATA);
2744 gmt_grd_pad_off (GMT, Z); /* Undo pad if one existed so we can sort */
2745 if (mode) { /* Already got the mode */
2746 wmode = *mode;
2747 n = Z->header->nm;
2748 }
2749 else {/* First compute the mode */
2750 gmt_sort_array (GMT, Z->data, Z->header->nm, GMT_FLOAT);
2751 for (n = Z->header->nm; n > 1 && gmt_M_is_fnan (Z->data[n-1]); n--);
2752 if (n)
2753 gmt_mode_f (GMT, Z->data, n, n/2, 0, gmt_mode_selection, &GMT_n_multiples, &wmode);
2754 else
2755 wmode = GMT->session.d_NaN;
2756 }
2757 gmt_getmad_f (GMT, Z->data, n, wmode, &lmsscl);
2758 if (!overwrite) gmt_free_grid (GMT, &Z, true);
2759 if (GMT_n_multiples > 0) GMT_Report (GMT->parent, GMT_MSG_WARNING, "%d Multiple modes found in the grid\n", GMT_n_multiples);
2760 }
2761 return lmsscl;
2762 }
2763
gmtstat_get_geo_cellarea(struct GMT_CTRL * GMT,struct GMT_GRID * G)2764 GMT_LOCAL void gmtstat_get_geo_cellarea (struct GMT_CTRL *GMT, struct GMT_GRID *G) {
2765 /* Calculate geographic SPHERICAL area in km^2 and place in grid G.
2766 * Integrating the area between two parallels +/- yinc/2 to either side of latitude y on a sphere
2767 * and partition it amount all cells in longitude yields the exact area (angles in radians)
2768 * A(y) = 2 * R^2 * xinc * sin (0.5 * yinc) * cos(y)
2769 * except at the points at the pole (gridline registration) or near the pole (pixel reg).
2770 * There, the integration yields
2771 * A(pole) = R^2 * xinc * (1.0 - cos (f * yinc))
2772 * where f = 0.5 for gridline-registered and f = 1 for pixel-registered grids.
2773 * P.Wessel, July 2016.
2774 */
2775 uint64_t node;
2776 unsigned int row, col, j, first_row = 0, last_row = G->header->n_rows - 1, last_col = G->header->n_columns - 1, ltype;
2777 double lat, area, f, row_weight, col_weight = 1.0, R2 = pow (0.001 * GMT->current.proj.mean_radius, 2.0); /* squared mean radius in km */
2778 char *aux[6] = {"geodetic", "authalic", "conformal", "meridional", "geocentric", "parametric"};
2779 char *rad[5] = {"mean (R_1)", "authalic (R_2)", "volumetric (R_3)", "meridional", "quadratic"};
2780
2781 ltype = (GMT->current.setting.proj_aux_latitude == GMT_LATSWAP_NONE) ? 0 : 1+GMT->current.setting.proj_aux_latitude/2;
2782 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Compute spherical gridnode areas using %s radius [R = %.12g km] and %s latitudes\n",
2783 rad[GMT->current.setting.proj_mean_radius], GMT->current.proj.mean_radius, aux[ltype]);
2784 /* May need special treatment of pole points */
2785 f = (G->header->registration == GMT_GRID_NODE_REG) ? 0.5 : 1.0; /* Half pizza-slice for gridline regs with node at pole, full slice for grids */
2786 area = R2 * (G->header->inc[GMT_X] * D2R);
2787 if (doubleAlmostEqualZero (G->header->wesn[YHI], 90.0)) { /* North pole row */
2788 row_weight = 1.0 - cosd (f * G->header->inc[GMT_Y]);
2789 gmt_M_col_loop (GMT, G, first_row, col, node) {
2790 if (G->header->registration == GMT_GRID_NODE_REG) col_weight = (col == 0 || col == last_col) ? 0.5 : 1.0;
2791 G->data[node] = (gmt_grdfloat)(row_weight * col_weight * area);
2792 }
2793 first_row++;
2794 }
2795 if (doubleAlmostEqualZero (G->header->wesn[YLO], -90.0)) { /* South pole row */
2796 row_weight = 1.0 - cosd (f * G->header->inc[GMT_Y]);
2797 gmt_M_col_loop (GMT, G, last_row, col, node) {
2798 if (G->header->registration == GMT_GRID_NODE_REG) col_weight = (col == 0 || col == last_col) ? 0.5 : 1.0;
2799 G->data[node] = (gmt_grdfloat)(row_weight * col_weight * area);
2800 }
2801 last_row--;
2802 }
2803 /* Below we just use the standard graticule equation */
2804 area *= 2.0 * sind (0.5 * G->header->inc[GMT_Y]); /* Since no longer any special cases with poles below */
2805 for (row = first_row, j = first_row + G->header->pad[YHI]; row <= last_row; row++, j++) {
2806 lat = gmt_M_grd_row_to_y (GMT, row, G->header);
2807 lat = gmt_lat_swap (GMT, lat, GMT->current.setting.proj_aux_latitude); /* Convert to selected auxiliary latitude */
2808 row_weight = cosd (lat);
2809 gmt_M_col_loop (GMT, G, row, col, node) { /* Loop over cols; always save the next left before we update the array at that col */
2810 if (G->header->registration == GMT_GRID_NODE_REG) col_weight = (col == 0 || col == last_col) ? 0.5 : 1.0;
2811 G->data[node] = (gmt_grdfloat)(row_weight * col_weight * area);
2812 }
2813 }
2814 }
2815
gmtstat_get_cart_cellarea(struct GMT_CTRL * GMT,struct GMT_GRID * G)2816 GMT_LOCAL void gmtstat_get_cart_cellarea (struct GMT_CTRL *GMT, struct GMT_GRID *G) {
2817 /* Calculate Cartesian cell areas in user units */
2818 uint64_t node;
2819 unsigned int row, col, last_row = G->header->n_rows - 1, last_col = G->header->n_columns - 1;
2820 double row_weight = 1.0, col_weight = 1.0, area = G->header->inc[GMT_X] * G->header->inc[GMT_Y]; /* All whole cells have same area */
2821 gmt_M_unused(GMT);
2822 gmt_M_row_loop (GMT, G, row) { /* Loop over the rows */
2823 if (G->header->registration == GMT_GRID_NODE_REG) row_weight = (row == 0 || row == last_row) ? 0.5 : 1.0; /* half-cells in y */
2824 gmt_M_col_loop (GMT, G, row, col, node) { /* Now loop over the columns */
2825 if (G->header->registration == GMT_GRID_NODE_REG) col_weight = (col == 0 || col == last_col) ? 0.5 : 1.0; /* half-cells in x */
2826 G->data[node] = (gmt_grdfloat)(row_weight * col_weight * area);
2827 }
2828 }
2829 }
2830
gmt_get_cellarea(struct GMT_CTRL * GMT,struct GMT_GRID * G)2831 void gmt_get_cellarea (struct GMT_CTRL *GMT, struct GMT_GRID *G) {
2832 /* Calculate geographic spherical in km^2 or plain Cartesian area in user_unit^2 and place in grid G. */
2833 if (gmt_M_is_geographic (GMT, GMT_IN))
2834 gmtstat_get_geo_cellarea (GMT, G);
2835 else
2836 gmtstat_get_cart_cellarea (GMT, G);
2837 }
2838
gmt_von_mises_mu_and_kappa(struct GMT_CTRL * GMT,double * data,double * w,uint64_t n,double * kappa)2839 double gmt_von_mises_mu_and_kappa (struct GMT_CTRL *GMT, double *data, double *w, uint64_t n, double *kappa) {
2840 /* Return the mean and kappa for a von Mises fit to (possibly weighted) data.
2841 * It is assumed that data have been scaled to 0-360. Weights w is possibly NULL for no weights */
2842 uint64_t k;
2843 double mean = 0.0, x_r = 0.0, y_r = 0.0, s_w = 0.0, ww = 1.0;
2844 double lo, hi, midval, range2, delta_R, s, c, R_bar;
2845
2846 for (k = 0; k < n; k++) {
2847 if (gmt_M_is_dnan (data[k])) continue;
2848 if (w) ww = w[k]; /* Otherwise it is a constant 1 */
2849 sincosd (data[k], &s, &c);
2850 x_r += c * ww; y_r += s * ww;
2851 s_w += ww;
2852 }
2853 if (s_w > 0.0) { /* Can compute the statistics */
2854 x_r /= s_w; y_r /= s_w;
2855 mean = atan2d (y_r, x_r);
2856 }
2857 else { /* No data, basically */
2858 *kappa = GMT->session.d_NaN;
2859 return (GMT->session.d_NaN);
2860 }
2861 /* Here we have actual values */
2862 R_bar = hypot (x_r, y_r);
2863 if (R_bar >= 0.999) { /* Just return a big kappa for R almost = 1 */
2864 *kappa = 500.0; /* kappa = 500 gives R_bar = 0.998999916722 so close enough */
2865 return (mean);
2866 }
2867 /* Compute kappa by a dumb bisection search */
2868 lo = 0.0; hi = 500.0;
2869 while (fabs (hi - lo) > GMT_CONV8_LIMIT) {
2870 midval = 0.5 * (hi + lo);
2871 range2 = 0.5 * (hi - lo);
2872 delta_R = (gmt_i1 (GMT, midval) / gmt_i0 (GMT, midval)) - R_bar;
2873 if (delta_R > GMT_CONV8_LIMIT) /* Need a smaller R next time */
2874 hi -= range2;
2875 else if (delta_R < -GMT_CONV8_LIMIT) /* Need a larger R next time */
2876 lo += range2;
2877 else /* Got it, set lo = hi to exit loop */
2878 lo = hi;
2879 }
2880 *kappa = midval;
2881
2882 return (mean);
2883 }
2884
2885