-
Notifications
You must be signed in to change notification settings - Fork 375
/
Copy pathsurface_old.c
1922 lines (1678 loc) · 79.6 KB
/
surface_old.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*--------------------------------------------------------------------
*
* Copyright (c) 1991-2025 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
* See LICENSE.TXT file for copying and redistribution conditions.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; version 3 or any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* Contact info: www.generic-mapping-tools.org
*--------------------------------------------------------------------*/
/*
* surface.c: a gridding program.
* reads xyz triples and fits a surface to the data.
* surface satisfies (1 - T) D4 z - T D2 z = 0,
* where D4 is the biharmonic operator,
* D2 is the Laplacian,
* and T is a "tension factor" between 0 and 1.
* End member T = 0 is the classical minimum curvature
* surface. T = 1 gives a harmonic surface. Use T = 0.25
* or so for potential data; something more for topography.
*
* Program includes overrelaxation for fast convergence and
* automatic optimal grid factorization.
*
* See Smith & Wessel (Geophysics, 3, 293-305, 1990) for details.
*
* Authors: Walter H. F. Smith and Paul Wessel
* Date: 1-JAN-2010
* Version: 6 API
* Tech Note: 1. surface uses old grid-structure (transpose of current GMT grid)
* and must therefore waste time transposing grid before writing.
* However, this is a tiny fraction of total run time.
*
* Note on KEYS: DD(= means -D takes an optional input Dataset as argument which may be followed by optional modifiers.
*/
#include "gmt_dev.h"
#define THIS_MODULE_CLASSIC_NAME "surface_old"
#define THIS_MODULE_MODERN_NAME "surface_old"
#define THIS_MODULE_LIB "core"
#define THIS_MODULE_PURPOSE "Grid table data using adjustable tension continuous curvature splines"
#define THIS_MODULE_KEYS "<D{,DD(=,LG(=1,GG}"
#define THIS_MODULE_NEEDS "R"
#define THIS_MODULE_OPTIONS "-:RVabdefhirs" GMT_OPT("FH")
struct SURFACE_CTRL {
struct SURFACE_OLD_A { /* -A<aspect_ratio>|m */
bool active;
unsigned int mode; /* 1 if we want cos(mid_latitude) */
double value;
} A;
struct SURFACE_OLD_C { /* -C<converge_limit> */
bool active;
unsigned int mode; /* 1 if given as fraction */
double value;
} C;
struct SURFACE_OLD_D { /* -D<line.xyz> */
bool active;
char *file; /* Name of file with breaklines */
} D;
struct SURFACE_OLD_G { /* -G<file> */
bool active;
char *file;
} G;
struct SURFACE_OLD_L { /* -Ll|u<limit> */
bool active;
char *low, *high;
double min, max;
unsigned int lmode, hmode;
} L;
struct SURFACE_OLD_N { /* -N<max_iterations> */
bool active;
unsigned int value;
} N;
struct SURFACE_OLD_Q { /* -Q */
bool active;
} Q;
struct SURFACE_OLD_S { /* -S<radius>[m|s] */
bool active;
double radius;
char unit;
} S;
struct SURFACE_OLD_T { /* -T<tension>[i][b] */
bool active;
double b_tension, i_tension;
} T;
struct SURFACE_OLD_Z { /* -Z<over_relaxation_parameter> */
bool active;
double value;
} Z;
};
#define SURFACE_OUTSIDE LONG_MAX /* Index number indicating data is outside usable area */
#define SURFACE_CONV_LIMIT 0.0001 /* Default is 100 ppm of data range as convergence criterion */
#define SURFACE_MAX_ITERATIONS 500 /* Default iterations at final grid size */
#define SURFACE_OVERRELAXATION 1.4 /* Default over-relaxation value */
#define SURFACE_CLOSENESS_FACTOR 0.05 /* A node is considered known if the nearest data is within 0.05 of a gridspacing of the node */
#define SURFACE_IS_UNCONSTRAINED 0 /* Various constants used to flag nearest node */
#define SURFACE_DATA_IS_IN_QUAD1 1
#define SURFACE_DATA_IS_IN_QUAD2 2
#define SURFACE_DATA_IS_IN_QUAD3 3
#define SURFACE_DATA_IS_IN_QUAD4 4
#define SURFACE_IS_CONSTRAINED 5
static unsigned int p[5][4] = { /* Indices into C->offset for each of the 4 quadrants, i.e. C->offset[kase][p[quadrant][k]]] */
{ 0, 0, 0, 0}, /* This row is not used */
{10, 9, 5, 1}, /* Indices for 1st quadrant */
{ 8, 9, 6, 3}, /* Indices for 2nd quadrant */
{ 1, 2, 6, 10}, /* Indices for 3rd quadrant */
{ 3, 2, 5, 8} /* Indices for 4th quadrant */
};
struct SURFACE_DATA { /* Data point and index to node it currently constrains */
gmt_grdfloat x;
gmt_grdfloat y;
gmt_grdfloat z;
uint64_t index;
#ifdef DEBUG /* For debugging purposes only - it is the original input data point number before sorting */
uint64_t number;
#endif
};
struct SURFACE_BRIGGS { /* Coefficients in Taylor series for Laplacian(z) a la I. C. Briggs (1974) */
double b[6];
};
struct SURFACE_GLOBAL { /* Things needed inside compare function must be global for now */
int block_ny; /* Number of nodes in y-dir for a given grid factor */
double grid_xinc, grid_yinc; /* size of each grid cell for a given grid factor */
double x_min, y_min; /* Lower left corner of grid */
} GMT_Surface_Global;
struct SURFACE_INFO { /* Control structure for surface setup and execution */
unsigned char *iu; /* Pointer to grid info array */
char mode_type[2]; /* D means include data points when iterating
* I means just interpolate from larger grid */
char format[GMT_BUFSIZ];
char *low_file, *high_file; /* Pointers to grids with low and high limits, if selected */
int grid, old_grid; /* Node spacings */
unsigned int n_fact; /* Number of factors in common (n_rows-1, n_columns-1) */
unsigned int factors[32]; /* Array of common factors */
unsigned int set_low; /* 0 unconstrained,1 = by min data value, 2 = by user value */
unsigned int set_high; /* 0 unconstrained,1 = by max data value, 2 = by user value */
size_t n_alloc;
uint64_t npoints; /* Number of data points */
uint64_t ij_sw_corner, ij_se_corner,ij_nw_corner, ij_ne_corner;
uint64_t n_empty; /* No of unconstrained nodes at initialization */
int n_columns; /* Number of nodes in x-dir. */
int n_rows; /* Number of nodes in y-dir. (Final grid) */
uint64_t nxny; /* Total number of grid nodes without boundaries */
int mx;
int my;
uint64_t mxmy; /* Total number of grid nodes with boundaries */
int block_nx; /* Number of nodes in x-dir for a given grid factor */
int block_ny; /* Number of nodes in y-dir for a given grid factor */
unsigned int max_iterations; /* Max iter per call to iterate */
unsigned int converge_limit_mode; /* 1 if -C set fractional convergence limit */
uint64_t total_iterations;
bool periodic; /* true if geographic grid and west-east == 360 */
int grid_east;
int offset[25][12]; /* Indices of 12 nearby points in 25 cases of edge conditions */
bool constrained; /* true if set_low or set_high is true */
double low_limit, high_limit; /* Constrains on range of solution */
double grid_xinc, grid_yinc; /* size of each grid cell for a given grid factor */
double r_grid_xinc, r_grid_yinc; /* Reciprocals */
double converge_limit; /* Convergence limit */
double radius; /* Search radius for initializing grid */
double tension; /* Tension parameter on the surface */
double boundary_tension;
double interior_tension;
double a0_const_1, a0_const_2; /* Constants for off grid point equation */
double e_2, e_m2, one_plus_e2;
double eps_p2, eps_m2, two_plus_ep2, two_plus_em2;
double x_edge_const, y_edge_const;
double l_epsilon;
double z_mean;
double z_scale; /* Root mean square range of z after removing planar trend */
double r_z_scale; /* reciprocal of z_scale */
double plane_c0, plane_c1, plane_c2; /* Coefficients of best fitting plane to data */
double coeff[2][12]; /* Coefficients for 12 nearby points, constrained and unconstrained */
double relax_old, relax_new; /* Coefficients for relaxation factor to speed up convergence */
double wesn_orig[4]; /* Original -R domain as we might have shifted it due to -r */
struct SURFACE_DATA *data;
struct SURFACE_BRIGGS *briggs;
struct GMT_GRID *Grid; /* The final grid */
struct GMT_GRID *Low, *High; /* arrays for minmax values, if set */
};
GMT_LOCAL void surface_set_coefficients (struct SURFACE_INFO *C) {
/* These are the coefficients in the finite-difference expressions */
double e_4, loose, a0;
loose = 1.0 - C->interior_tension;
C->e_2 = C->l_epsilon * C->l_epsilon;
e_4 = C->e_2 * C->e_2;
C->eps_p2 = C->e_2;
C->eps_m2 = 1.0/C->e_2;
C->one_plus_e2 = 1.0 + C->e_2;
C->two_plus_ep2 = 2.0 + 2.0*C->eps_p2;
C->two_plus_em2 = 2.0 + 2.0*C->eps_m2;
C->x_edge_const = 4 * C->one_plus_e2 - 2 * (C->interior_tension / loose);
C->e_m2 = 1.0 / C->e_2;
C->y_edge_const = 4 * (1.0 + C->e_m2) - 2 * (C->interior_tension * C->e_m2 / loose);
a0 = 1.0 / ( (6 * e_4 * loose + 10 * C->e_2 * loose + 8 * loose - 2 * C->one_plus_e2) + 4*C->interior_tension*C->one_plus_e2);
C->a0_const_1 = 2 * loose * (1.0 + e_4);
C->a0_const_2 = 2.0 - C->interior_tension + 2 * loose * C->e_2;
C->coeff[1][4] = C->coeff[1][7] = -loose;
C->coeff[1][0] = C->coeff[1][11] = -loose * e_4;
C->coeff[0][4] = C->coeff[0][7] = -loose * a0;
C->coeff[0][0] = C->coeff[0][11] = -loose * e_4 * a0;
C->coeff[1][5] = C->coeff[1][6] = 2 * loose * C->one_plus_e2;
C->coeff[0][5] = C->coeff[0][6] = (2 * C->coeff[1][5] + C->interior_tension) * a0;
C->coeff[1][2] = C->coeff[1][9] = C->coeff[1][5] * C->e_2;
C->coeff[0][2] = C->coeff[0][9] = C->coeff[0][5] * C->e_2;
C->coeff[1][1] = C->coeff[1][3] = C->coeff[1][8] = C->coeff[1][10] = -2 * loose * C->e_2;
C->coeff[0][1] = C->coeff[0][3] = C->coeff[0][8] = C->coeff[0][10] = C->coeff[1][1] * a0;
C->e_2 *= 2; /* We will need these in boundary conditions */
C->e_m2 *= 2;
C->ij_sw_corner = 2 * C->my + 2; /* Corners of array of actual data */
C->ij_se_corner = C->ij_sw_corner + (C->n_columns - 1) * C->my;
C->ij_nw_corner = C->ij_sw_corner + C->n_rows - 1;
C->ij_ne_corner = C->ij_se_corner + C->n_rows - 1;
}
GMT_LOCAL void surface_set_offset (struct SURFACE_INFO *C) {
/* Because of the multigrid approach the distance from the central node to
* its neighbors needed in the finite difference expressions varies. E.g.,
* when C->grid is 8 then the next neighbor to the right is 8 columns over.
* But at the boundaries the spacing is always 1. Thus, as the current node
* moves over the interior of the grid the distances change as we get close
* to any of the 4 boundaries. The offset array is used to determine what
* the offset in rows and columns are relative to the current point, given
* what "kase" we are examining. kase is a combination of x and y information
* related to how close we are to the left/right or top/bottom boundary.
*/
int add_w[5], add_e[5], add_s[5], add_n[5], add_w2[5], add_e2[5], add_s2[5], add_n2[5];
unsigned int i, j, kase;
add_w[0] = -C->my; add_w[1] = add_w[2] = add_w[3] = add_w[4] = -C->grid_east;
add_w2[0] = -2 * C->my; add_w2[1] = -C->my - C->grid_east; add_w2[2] = add_w2[3] = add_w2[4] = -2 * C->grid_east;
add_e[4] = C->my; add_e[0] = add_e[1] = add_e[2] = add_e[3] = C->grid_east;
add_e2[4] = 2 * C->my; add_e2[3] = C->my + C->grid_east; add_e2[2] = add_e2[1] = add_e2[0] = 2 * C->grid_east;
add_n[4] = 1; add_n[3] = add_n[2] = add_n[1] = add_n[0] = C->grid;
add_n2[4] = 2; add_n2[3] = C->grid + 1; add_n2[2] = add_n2[1] = add_n2[0] = 2 * C->grid;
add_s[0] = -1; add_s[1] = add_s[2] = add_s[3] = add_s[4] = -C->grid;
add_s2[0] = -2; add_s2[1] = -C->grid - 1; add_s2[2] = add_s2[3] = add_s2[4] = -2 * C->grid;
for (i = 0, kase = 0; i < 5; i++) {
for (j = 0; j < 5; j++, kase++) {
C->offset[kase][0] = add_n2[j];
C->offset[kase][1] = add_n[j] + add_w[i];
C->offset[kase][2] = add_n[j];
C->offset[kase][3] = add_n[j] + add_e[i];
C->offset[kase][4] = add_w2[i];
C->offset[kase][5] = add_w[i];
C->offset[kase][6] = add_e[i];
C->offset[kase][7] = add_e2[i];
C->offset[kase][8] = add_s[j] + add_w[i];
C->offset[kase][9] = add_s[j];
C->offset[kase][10] = add_s[j] + add_e[i];
C->offset[kase][11] = add_s2[j];
}
}
}
GMT_LOCAL void fill_in_forecast (struct SURFACE_INFO *C) {
/* Fills in bilinear estimates into new node locations
after grid is divided. These new nodes are marked as
unconstrained while the coarser data are considered
constraints in the next iteration.
*/
uint64_t index_0, index_1, index_2, index_3, index_new;
int ii, jj, i, j;
unsigned char *iu = C->iu;
double delta_x, delta_y, a0, a1, a2, a3, old_size, a0_plus_a1_dx, a2_plus_a3_dx;
gmt_grdfloat *u = C->Grid->data;
old_size = 1.0 / (double)C->old_grid;
/* First do from southwest corner */
for (i = 0; i < (C->n_columns-1); i += C->old_grid) {
for (j = 0; j < (C->n_rows-1); j += C->old_grid) {
/* Get indices of bilinear square */
index_0 = C->ij_sw_corner + i * C->my + j;
index_1 = index_0 + C->old_grid * C->my;
index_2 = index_1 + C->old_grid;
index_3 = index_0 + C->old_grid;
/* Get coefficients */
a0 = u[index_0];
a1 = u[index_1] - a0;
a2 = u[index_3] - a0;
a3 = u[index_2] - a0 - a1 - a2;
/* Find all possible new fill ins */
for (ii = i; ii < (i + C->old_grid); ii += C->grid) {
delta_x = (ii - i) * old_size;
a0_plus_a1_dx = a0 + a1 * delta_x;
a2_plus_a3_dx = a2 + a3 * delta_x;
for (jj = j; jj < (j + C->old_grid); jj += C->grid) {
index_new = C->ij_sw_corner + ii * C->my + jj;
if (index_new == index_0) continue;
delta_y = (jj - j) * old_size;
u[index_new] = (gmt_grdfloat)(a0_plus_a1_dx + delta_y * a2_plus_a3_dx);
iu[index_new] = SURFACE_IS_UNCONSTRAINED;
}
}
iu[index_0] = SURFACE_IS_CONSTRAINED;
}
}
/* Now do linear guess along east edge */
for (j = 0; j < (C->n_rows-1); j += C->old_grid) {
index_0 = C->ij_se_corner + j;
index_3 = index_0 + C->old_grid;
for (jj = j; jj < j + C->old_grid; jj += C->grid) {
index_new = C->ij_se_corner + jj;
delta_y = (jj - j) * old_size;
u[index_new] = u[index_0] + (gmt_grdfloat)(delta_y * (u[index_3] - u[index_0]));
iu[index_new] = SURFACE_IS_UNCONSTRAINED;
}
iu[index_0] = SURFACE_IS_CONSTRAINED;
}
/* Now do linear guess along north edge */
for (i = 0; i < (C->n_columns-1); i += C->old_grid) {
index_0 = C->ij_nw_corner + i * C->my;
index_1 = index_0 + C->old_grid * C->my;
for (ii = i; ii < i + C->old_grid; ii += C->grid) {
index_new = C->ij_nw_corner + ii * C->my;
delta_x = (ii - i) * old_size;
u[index_new] = u[index_0] + (gmt_grdfloat)(delta_x * (u[index_1] - u[index_0]));
iu[index_new] = SURFACE_IS_UNCONSTRAINED;
}
iu[index_0] = SURFACE_IS_CONSTRAINED;
}
/* Now set northeast corner to fixed and we're done */
iu[C->ij_ne_corner] = SURFACE_IS_CONSTRAINED;
}
GMT_LOCAL int surface_compare_points (const void *point_1v, const void *point_2v) {
/* Routine for qsort to sort data structure for fast access to data by node location.
Sorts on index first, then on radius to node corresponding to index, so that index
goes from low to high, and so does radius.
*/
uint64_t block_i, block_j, index_1, index_2;
double x0, y0, dist_1, dist_2;
const struct SURFACE_DATA *point_1 = point_1v, *point_2 = point_2v;
index_1 = point_1->index;
index_2 = point_2->index;
if (index_1 < index_2) return (-1);
if (index_1 > index_2) return (1);
if (index_1 == SURFACE_OUTSIDE) return (0);
/* Points are in same grid cell, find the one who is nearest to grid point */
block_i = point_1->index/GMT_Surface_Global.block_ny;
block_j = point_1->index%GMT_Surface_Global.block_ny;
x0 = GMT_Surface_Global.x_min + block_i * GMT_Surface_Global.grid_xinc;
y0 = GMT_Surface_Global.y_min + block_j * GMT_Surface_Global.grid_yinc;
dist_1 = (point_1->x - x0) * (point_1->x - x0) + (point_1->y - y0) * (point_1->y - y0);
dist_2 = (point_2->x - x0) * (point_2->x - x0) + (point_2->y - y0) * (point_2->y - y0);
if (dist_1 < dist_2) return (-1);
if (dist_1 > dist_2) return (1);
return (0);
}
GMT_LOCAL void surface_smart_divide (struct SURFACE_INFO *C) {
/* Divide grid by its largest prime factor */
C->grid /= C->factors[C->n_fact - 1];
C->n_fact--;
}
GMT_LOCAL void surface_set_index (struct SURFACE_INFO *C) {
/* recomputes data[k].index for new value of grid,
sorts data on index and radii, and throws away
data which are now outside the usable limits. */
int i, j;
uint64_t k, k_skipped = 0;
struct GMT_GRID_HEADER *h = C->Grid->header;
for (k = 0; k < C->npoints; k++) {
i = irint (floor(((C->data[k].x-h->wesn[XLO])*C->r_grid_xinc) + 0.5));
j = irint (floor(((C->data[k].y-h->wesn[YLO])*C->r_grid_yinc) + 0.5));
if (i < 0 || i >= C->block_nx || j < 0 || j >= C->block_ny) {
C->data[k].index = SURFACE_OUTSIDE;
k_skipped++;
}
else
C->data[k].index = i * C->block_ny + j;
}
qsort (C->data, C->npoints, sizeof (struct SURFACE_DATA), surface_compare_points);
C->npoints -= k_skipped;
}
GMT_LOCAL void find_nearest_point (struct SURFACE_INFO *C) {
/* Determines the nearest data point epr bin and sets the
* Briggs parameters or, if really close, sets the node value */
uint64_t ij_v2, k, last_index, iu_index, briggs_index;
int i, j, block_i, block_j;
double x0, y0, dx, dy, dxpdy, xys, xy1, btemp, *b = NULL;
gmt_grdfloat z_at_node, *u = C->Grid->data;
unsigned char *iu = C->iu;
struct GMT_GRID_HEADER *h = C->Grid->header;
last_index = UINTMAX_MAX;
for (i = 0; i < C->n_columns; i += C->grid) /* Reset grid info */
for (j = 0; j < C->n_rows; j += C->grid)
iu[C->ij_sw_corner + i*C->my + j] = SURFACE_IS_UNCONSTRAINED;
briggs_index = 0;
for (k = 0; k < C->npoints; k++) { /* Find constraining value */
if (C->data[k].index != last_index) {
block_i = (int)C->data[k].index/C->block_ny;
block_j = (int)C->data[k].index%C->block_ny;
last_index = C->data[k].index;
iu_index = C->ij_sw_corner + (block_i * C->my + block_j) * C->grid;
x0 = h->wesn[XLO] + block_i*C->grid_xinc;
y0 = h->wesn[YLO] + block_j*C->grid_yinc;
dx = (C->data[k].x - x0)*C->r_grid_xinc;
dy = (C->data[k].y - y0)*C->r_grid_yinc;
if (fabs(dx) < SURFACE_CLOSENESS_FACTOR && fabs(dy) < SURFACE_CLOSENESS_FACTOR) { /* Close enough to assign value to node */
iu[iu_index] = SURFACE_IS_CONSTRAINED;
/* v3.3.4: NEW CODE
* Since point is basically moved from (dx, dy) to (0,0) we must adjust for
* the small change in the planar trend between the two locations, and then
* possibly clip the range if constraining surfaces were given. Note that
* dx, dy is in -1/1 range normalized by (grid * x|y_inc) so to recover the
* dx,dy in final grid fractions we must scale by grid */
z_at_node = C->data[k].z + (gmt_grdfloat) (C->r_z_scale * C->grid * (C->plane_c1 * dx + C->plane_c2 * dy));
if (C->constrained) { /* Must use ij_v2 since constrained grids are in standard scanline format */
ij_v2 = gmt_M_ijp (C->Grid->header, C->n_rows - block_j * C->grid - 1, block_i * C->grid);
if (C->set_low && !gmt_M_is_fnan (C->Low->data[ij_v2]) && z_at_node < C->Low->data[ij_v2])
z_at_node = C->Low->data[ij_v2];
else if (C->set_high && !gmt_M_is_fnan (C->High->data[ij_v2]) && z_at_node > C->High->data[ij_v2])
z_at_node = C->High->data[ij_v2];
}
u[iu_index] = z_at_node;
}
else { /* We have a nearby data point in one of the quadrants */
if (dx >= 0.0) {
if (dy >= 0.0)
iu[iu_index] = SURFACE_DATA_IS_IN_QUAD1;
else {
iu[iu_index] = SURFACE_DATA_IS_IN_QUAD4;
//dy = -dy; gmt_M_double_swap (dx, dy);
}
}
else {
//dx = -dx;
if (dy >= 0.0) {
iu[iu_index] = SURFACE_DATA_IS_IN_QUAD2;
//gmt_M_double_swap (dx, dy);
}
else {
iu[iu_index] = SURFACE_DATA_IS_IN_QUAD3;
//dy = -dy;
}
}
/* Evaluate the Briggs coefficients */
dx = fabs (dx);
dy = fabs (dy);
dxpdy = dx + dy;
xys = 1.0 + dxpdy;
btemp = 2 * C->one_plus_e2 / ( dxpdy * xys );
b = C->briggs[briggs_index].b; /* Shorthand to this Briggs-array */
b[0] = 1.0 - 0.5 * (dx + (dx * dx)) * btemp;
b[3] = 0.5 * (C->e_2 - (dy + (dy * dy)) * btemp);
xy1 = 1.0 / xys;
b[1] = (C->e_2 * xys - 4 * dy) * xy1;
b[2] = 2 * (dy - dx + 1.0) * xy1;
b[4] = b[0] + b[1] + b[2] + b[3] + btemp;
b[5] = btemp * C->data[k].z;
b[4] = 1.0 / (C->a0_const_1 + C->a0_const_2 * b[4]); /* Do this calculation here instead of inside iterate loop */
briggs_index++;
}
}
}
}
GMT_LOCAL void surface_set_grid_parameters (struct SURFACE_INFO *C) {
/* Updates the grid space parameters given the new C->grid setting */
GMT_Surface_Global.block_ny = C->block_ny = (C->n_rows - 1) / C->grid + 1;
C->block_nx = (C->n_columns - 1) / C->grid + 1;
GMT_Surface_Global.grid_xinc = C->grid_xinc = C->grid * C->Grid->header->inc[GMT_X];
GMT_Surface_Global.grid_yinc = C->grid_yinc = C->grid * C->Grid->header->inc[GMT_Y];
C->grid_east = C->grid * C->my;
C->r_grid_xinc = 1.0 / C->grid_xinc;
C->r_grid_yinc = 1.0 / C->grid_yinc;
}
GMT_LOCAL void surface_initialize_grid (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
/*
* For the initial gridsize, compute weighted averages of data inside the search radius
* and assign the values to u[i,j] where i,j are multiples of gridsize.
*/
uint64_t index_1, index_2, k, k_index;
int irad, jrad, i, j, imin, imax, jmin, jmax, ki, kj;
double r, rfact, sum_w, sum_zw, weight, x0, y0;
gmt_grdfloat *u = C->Grid->data;
struct GMT_GRID_HEADER *h = C->Grid->header;
irad = irint (ceil(C->radius/C->grid_xinc));
jrad = irint (ceil(C->radius/C->grid_yinc));
rfact = -4.5/(C->radius*C->radius);
for (i = 0; i < C->block_nx; i ++ ) {
x0 = h->wesn[XLO] + i*C->grid_xinc;
for (j = 0; j < C->block_ny; j ++ ) {
y0 = h->wesn[YLO] + j*C->grid_yinc;
imin = i - irad;
if (imin < 0) imin = 0;
imax = i + irad;
if (imax >= C->block_nx) imax = C->block_nx - 1;
jmin = j - jrad;
if (jmin < 0) jmin = 0;
jmax = j + jrad;
if (jmax >= C->block_ny) jmax = C->block_ny - 1;
index_1 = (uint64_t)imin*C->block_ny + jmin;
index_2 = (uint64_t)imax*C->block_ny + jmax + 1;
sum_w = sum_zw = 0.0;
k = 0;
while (k < C->npoints && C->data[k].index < index_1) k++;
for (ki = imin; k < C->npoints && ki <= imax && C->data[k].index < index_2; ki++) {
for (kj = jmin; k < C->npoints && kj <= jmax && C->data[k].index < index_2; kj++) {
k_index = (uint64_t)ki*C->block_ny + kj;
while (k < C->npoints && C->data[k].index < k_index) k++;
while (k < C->npoints && C->data[k].index == k_index) {
r = (C->data[k].x-x0)*(C->data[k].x-x0) + (C->data[k].y-y0)*(C->data[k].y-y0);
weight = exp (rfact*r);
sum_w += weight;
sum_zw += weight*C->data[k].z;
k++;
}
}
}
if (sum_w == 0.0) {
sprintf (C->format, "No data inside search radius at: %s %s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
GMT_Report (GMT->parent, GMT_MSG_ERROR, C->format, x0, y0);
u[C->ij_sw_corner + (i * C->my + j) * C->grid] = (gmt_grdfloat)C->z_mean;
}
else {
u[C->ij_sw_corner + (i*C->my+j)*C->grid] = (gmt_grdfloat)(sum_zw/sum_w);
}
}
}
}
GMT_LOCAL int surface_read_data (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, struct GMT_OPTION *options) {
/* Procdss input data into data structure */
int i, j, error;
uint64_t k, kmax = 0, kmin = 0, n_dup = 0;
double *in, half_dx, zmin = DBL_MAX, zmax = -DBL_MAX, wesn_lim[4];
struct GMT_GRID_HEADER *h = C->Grid->header;
struct GMT_RECORD *In = NULL;
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Processing input table data\n");
C->data = gmt_M_memory (GMT, NULL, C->n_alloc, struct SURFACE_DATA);
/* Read in xyz data and computes index no and store it in a structure */
if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
return (error);
}
if (GMT_Init_IO (GMT->parent, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
return (GMT->parent->error);
}
k = 0;
C->z_mean = 0.0;
/* Initially allow points to be within 1 grid spacing of the grid */
wesn_lim[XLO] = h->wesn[XLO] - C->grid_xinc; wesn_lim[XHI] = h->wesn[XHI] + C->grid_xinc;
wesn_lim[YLO] = h->wesn[YLO] - C->grid_yinc; wesn_lim[YHI] = h->wesn[YHI] + C->grid_yinc;
half_dx = 0.5 * C->grid_xinc;
if (GMT_Begin_IO (GMT->parent, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
return (GMT->parent->error);
}
do { /* Keep returning records until we reach EOF */
if ((In = GMT_Get_Record (GMT->parent, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */
return (GMT_RUNTIME_ERROR);
else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
break;
continue; /* Go back and read the next record */
}
/* Data record to process */
in = In->data; /* Only need to process numerical part here */
if (gmt_M_is_dnan (in[GMT_Z])) continue;
if (gmt_M_y_is_outside (GMT, in[GMT_Y], wesn_lim[YLO], wesn_lim[YHI])) continue; /* Outside y-range */
if (gmt_x_is_outside (GMT, &in[GMT_X], wesn_lim[XLO], wesn_lim[XHI])) continue; /* Outside x-range (or longitude) */
if (C->periodic && (h->wesn[XHI]-in[GMT_X] < half_dx)) { /* Push all values to the western nodes */
in[GMT_X] -= 360.0; /* Make this point constraining the western node value and then duplicate later */
i = 0;
}
else
i = irint (floor(((in[GMT_X]-h->wesn[XLO])*C->r_grid_xinc) + 0.5));
if (i < 0 || i >= C->block_nx) continue;
j = irint (floor(((in[GMT_Y]-h->wesn[YLO])*C->r_grid_yinc) + 0.5));
if (j < 0 || j >= C->block_ny) continue;
C->data[k].index = i * C->block_ny + j;
#ifdef DEBUG
C->data[k].number = k + 1;
#endif
C->data[k].x = (gmt_grdfloat)in[GMT_X];
C->data[k].y = (gmt_grdfloat)in[GMT_Y];
C->data[k].z = (gmt_grdfloat)in[GMT_Z];
if (zmin > in[GMT_Z]) zmin = in[GMT_Z], kmin = k;
if (zmax < in[GMT_Z]) zmax = in[GMT_Z], kmax = k;
k++;
C->z_mean += in[GMT_Z];
if (k == C->n_alloc) {
C->n_alloc <<= 1;
C->data = gmt_M_memory (GMT, C->data, C->n_alloc, struct SURFACE_DATA);
}
if (C->periodic && i == 0) { /* Replicate information to eastern boundary */
i = C->block_nx - 1;
C->data[k].index = i * C->block_ny + j;
#ifdef DEBUG
C->data[k].number = k + 1;
#endif
C->data[k].x = (gmt_grdfloat)(in[GMT_X] + 360.0);
C->data[k].y = (gmt_grdfloat)in[GMT_Y];
C->data[k].z = (gmt_grdfloat)in[GMT_Z];
if (zmin > in[GMT_Z]) zmin = in[GMT_Z], kmin = k;
if (zmax < in[GMT_Z]) zmax = in[GMT_Z], kmax = k;
k++;
C->z_mean += in[GMT_Z];
if (k == C->n_alloc) {
C->n_alloc <<= 1;
C->data = gmt_M_memory (GMT, C->data, C->n_alloc, struct SURFACE_DATA);
}
n_dup++;
}
} while (true);
if (GMT_End_IO (GMT->parent, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
return (GMT->parent->error);
}
C->npoints = k;
if (C->npoints == 0) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, " No datapoints inside region, aborts\n");
return (GMT_RUNTIME_ERROR);
}
C->z_mean /= k;
if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
sprintf (C->format, "%s %s %s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Minimum value of your dataset x,y,z at: ");
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, C->format, (double)C->data[kmin].x, (double)C->data[kmin].y, (double)C->data[kmin].z);
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Maximum value of your dataset x,y,z at: ");
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, C->format, (double)C->data[kmax].x, (double)C->data[kmax].y, (double)C->data[kmax].z);
if (C->periodic && n_dup) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Number of input values shared between repeating west and east column nodes: %" PRIu64 "\n", n_dup);
}
C->data = gmt_M_memory (GMT, C->data, C->npoints, struct SURFACE_DATA);
if (C->set_low == 1)
C->low_limit = C->data[kmin].z;
else if (C->set_low == 2 && C->low_limit > C->data[kmin].z) {
/* C->low_limit = data[kmin].z; */
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Your lower value is > than min data value.\n");
}
if (C->set_high == 1)
C->high_limit = C->data[kmax].z;
else if (C->set_high == 2 && C->high_limit < C->data[kmax].z) {
/* C->high_limit = data[kmax].z; */
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Your upper value is < than max data value.\n");
}
return (0);
}
GMT_LOCAL int surface_load_constraints (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, int transform) {
/* Deal with the constants or grids supplied via -L */
unsigned int i, j;
uint64_t ij;
double yy;
struct GMTAPI_CTRL *API = GMT->parent;
/* Load lower/upper limits, verify range, deplane, and rescale */
if (C->set_low > 0) {
if (C->set_low < 3) {
if ((C->Low = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, C->Grid)) == NULL) return (API->error);
for (ij = 0; ij < C->mxmy; ij++) C->Low->data[ij] = (gmt_grdfloat)C->low_limit;
}
else {
if ((C->Low = GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, C->low_file, NULL)) == NULL) return (API->error); /* Get header only */
if (C->Low->header->n_columns != C->Grid->header->n_columns || C->Low->header->n_rows != C->Grid->header->n_rows) {
GMT_Report (API, GMT_MSG_ERROR, "Lower limit file not of proper dimension!\n");
return (GMT_RUNTIME_ERROR);
}
if (GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, C->low_file, C->Low) == NULL) return (API->error);
}
if (transform) {
for (j = 0; j < C->Grid->header->n_rows; j++) {
yy = (double)(C->Grid->header->n_rows - j - 1);
for (i = 0; i < C->Grid->header->n_columns; i++) {
ij = gmt_M_ijp (C->Grid->header, j, i);
if (gmt_M_is_fnan (C->Low->data[ij])) continue;
C->Low->data[ij] -= (gmt_grdfloat)(C->plane_c0 + C->plane_c1 * i + C->plane_c2 * yy);
C->Low->data[ij] *= (gmt_grdfloat)C->r_z_scale;
}
}
}
C->constrained = true;
}
if (C->set_high > 0) {
if (C->set_high < 3) {
if ((C->High = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, C->Grid)) == NULL) return (API->error);
for (ij = 0; ij < C->mxmy; ij++) C->High->data[ij] = (gmt_grdfloat)C->high_limit;
}
else {
if ((C->High = GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, C->high_file, NULL)) == NULL) return (API->error); /* Get header only */
if (C->High->header->n_columns != C->Grid->header->n_columns || C->High->header->n_rows != C->Grid->header->n_rows) {
GMT_Report (API, GMT_MSG_ERROR, "Upper limit file not of proper dimension!\n");
return (GMT_RUNTIME_ERROR);
}
if (GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, C->high_file, C->High) == NULL) return (API->error);
}
if (transform) {
for (j = 0; j < C->Grid->header->n_rows; j++) {
yy = (double)(C->n_rows - j - 1);
for (i = 0; i < C->Grid->header->n_columns; i++) {
ij = gmt_M_ijp (C->Grid->header, j, i);
if (gmt_M_is_fnan (C->High->data[ij])) continue;
C->High->data[ij] -= (gmt_grdfloat)(C->plane_c0 + C->plane_c1 * i + C->plane_c2 * yy);
C->High->data[ij] *= (gmt_grdfloat)C->r_z_scale;
}
}
}
C->constrained = true;
}
return (0);
}
GMT_LOCAL int surface_write_grid (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, char *grdfile) {
/* Uses v.2.0 netCDF grd format - hence need to transpose original column grid to be GMT compatible. This will be rewritten, maybe */
uint64_t index, k;
int i, j, err;
gmt_grdfloat *u = C->Grid->data, *v2 = NULL;
if ((err = surface_load_constraints (GMT, C, false)) != 0) return (err); /* Reload constraints but this time do not transform data */
strcpy (C->Grid->header->title, "Data gridded with continuous surface splines in tension");
v2 = gmt_M_memory_aligned (GMT, NULL, C->Grid->header->size, gmt_grdfloat);
index = C->ij_sw_corner;
if (GMT->common.R.active[GSET]) { /* Pixel registration request. Reset limits to the original extents */
gmt_M_memcpy (C->Grid->header->wesn, C->wesn_orig, 4, double);
C->Grid->header->registration = GMT->common.R.registration;
/* Must reduce n_columns,n_rows by 1 to exclude the extra padding for pixel grids */
C->Grid->header->n_columns--; C->n_columns--;
C->Grid->header->n_rows--; C->n_rows--;
gmt_set_grddim (GMT, C->Grid->header); /* Reset all integer dimensions and xy_off */
}
for (i = 0; i < C->n_columns; i++, index += C->my) {
for (j = 0; j < C->n_rows; j++) {
k = gmt_M_ijp (C->Grid->header, j, i);
v2[k] = u[index + C->n_rows - j - 1];
if (C->set_low && !gmt_M_is_fnan (C->Low->data[k]) && v2[k] < C->Low->data[k]) v2[k] = C->Low->data[k];
if (C->set_high && !gmt_M_is_fnan (C->High->data[k]) && v2[k] > C->High->data[k]) v2[k] = C->High->data[k];
}
}
if (C->periodic) { /* Ensure periodicity of E-W boundaries */
for (j = 0; j < C->n_rows; j++) {
k = gmt_M_ijp (C->Grid->header, j, 0);
v2[k] = v2[k+C->n_columns-1] = (gmt_grdfloat)(0.5 * (v2[k] + v2[k+C->n_columns-1])); /* Set these to the same as their average */
}
}
gmt_M_free_aligned (GMT, C->Grid->data); /* Free original column-oriented grid */
C->Grid->data = v2; /* Hook in new scanline-oriented grid */
if (GMT_Write_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, grdfile, C->Grid) != GMT_NOERROR) {
return (GMT->parent->error);
}
if ((C->set_low > 0 && C->set_low < 3) && GMT_Destroy_Data (GMT->parent, &C->Low) != GMT_NOERROR) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to free C->Low\n");
}
if ((C->set_high > 0 && C->set_high < 3) && GMT_Destroy_Data (GMT->parent, &C->High) != GMT_NOERROR) {
GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to free C->High\n");
}
return (0);
}
GMT_LOCAL uint64_t surface_iterate (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, int mode) {
/* Main finite difference solver */
uint64_t ij, briggs_index, ij_v2, iteration_count = 0, ij_sw, ij_se;
unsigned int current_max_iterations = C->max_iterations * C->grid;
int i, j, k, kase;
int x_case, y_case, x_w_case, x_e_case, y_s_case, y_n_case;
unsigned char *iu = C->iu;
double current_limit = C->converge_limit / C->grid;
double change, max_change = 0.0, busum, sum_ij;
double *b = NULL;
gmt_grdfloat *u = C->Grid->data;
double x_0_const = 4.0 * (1.0 - C->boundary_tension) / (2.0 - C->boundary_tension);
double x_1_const = (3 * C->boundary_tension - 2.0) / (2.0 - C->boundary_tension);
double y_denom = 2 * C->l_epsilon * (1.0 - C->boundary_tension) + C->boundary_tension;
double y_0_const = 4 * C->l_epsilon * (1.0 - C->boundary_tension) / y_denom;
double y_1_const = (C->boundary_tension - 2 * C->l_epsilon * (1.0 - C->boundary_tension) ) / y_denom;
sprintf (C->format, "%%4ld\t%%c\t%%8" PRIu64 "\t%s\t%s\t%%10" PRIu64 "\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
do {
briggs_index = 0; /* Reset the constraint table stack pointer */
max_change = -1.0;
/* Fill in auxiliary boundary values (in new way) */
/* First set (1-T)d2[]/dn2 + Td[]/dn = 0 along edges */
for (i = 0; i < C->n_columns; i += C->grid) {
/* set BC on south side */
ij = C->ij_sw_corner + i * C->my;
/* u[ij - 1] = 2 * u[ij] - u[ij + grid]; */
u[ij - 1] = (gmt_grdfloat)(y_0_const * u[ij] + y_1_const * u[ij + C->grid]);
/* set BC on north side */
ij = C->ij_nw_corner + i * C->my;
/* u[ij + 1] = 2 * u[ij] - u[ij - grid]; */
u[ij + 1] = (gmt_grdfloat)(y_0_const * u[ij] + y_1_const * u[ij - C->grid]);
}
if (C->periodic) { /* Set periodic boundary conditions in longitude at west and east boundaries */
for (j = 0; j < C->n_rows; j += C->grid) {
ij_sw = C->ij_sw_corner + j;
ij_se = C->ij_se_corner + j;
u[ij_sw+C->offset[0][5]] = u[ij_se+C->offset[20][5]];
u[ij_se+C->offset[20][6]] = u[ij_sw+C->offset[0][6]];
u[ij_se] = u[ij_sw] = 0.5f * (u[ij_se] + u[ij_sw]); /* Set to average of east and west */
}
}
else { /* Regular natural BC */
for (j = 0; j < C->n_rows; j += C->grid) {
/* set BC on west side */
ij = C->ij_sw_corner + j;
/* u[ij - my] = 2 * u[ij] - u[ij + grid_east]; */
u[ij - C->my] = (gmt_grdfloat)(x_1_const * u[ij + C->grid_east] + x_0_const * u[ij]);
/* set BC on east side */
ij = C->ij_se_corner + j;
/* u[ij + my] = 2 * u[ij] - u[ij - grid_east]; */
u[ij + C->my] = (gmt_grdfloat)(x_1_const * u[ij - C->grid_east] + x_0_const * u[ij]);
}
}
/* Now set d2[]/dxdy = 0 at each of the 4 corners */
ij = C->ij_sw_corner;
u[ij - C->my - 1] = u[ij + C->grid_east - 1] + u[ij - C->my + C->grid] - u[ij + C->grid_east + C->grid];
ij = C->ij_nw_corner;
u[ij - C->my + 1] = u[ij + C->grid_east + 1] + u[ij - C->my - C->grid] - u[ij + C->grid_east - C->grid];
ij = C->ij_se_corner;
u[ij + C->my - 1] = u[ij - C->grid_east - 1] + u[ij + C->my + C->grid] - u[ij - C->grid_east + C->grid];
ij = C->ij_ne_corner;
u[ij + C->my + 1] = u[ij - C->grid_east + 1] + u[ij + C->my - C->grid] - u[ij - C->grid_east - C->grid];
/* Now set (1-T)dC/dn + Tdu/dn = 0 at each edge */
/* New experiment: only dC/dn = 0 */
x_w_case = 0;
x_e_case = C->block_nx - 1;
for (i = 0; i < C->n_columns; i += C->grid, x_w_case++, x_e_case--) {
if(x_w_case < 2)
x_case = x_w_case;
else if(x_e_case < 2)
x_case = 4 - x_e_case;
else
x_case = 2;
/* South side */
kase = x_case * 5;
ij = C->ij_sw_corner + i * C->my;
u[ij + C->offset[kase][11]] =
(gmt_grdfloat)(u[ij + C->offset[kase][0]] + C->eps_m2*(u[ij + C->offset[kase][1]] + u[ij + C->offset[kase][3]]
- u[ij + C->offset[kase][8]] - u[ij + C->offset[kase][10]])
+ C->two_plus_em2 * (u[ij + C->offset[kase][9]] - u[ij + C->offset[kase][2]]) );
/* + tense * C->eps_m2 * (u[ij + C->offset[kase][2]] - u[ij + C->offset[kase][9]]) / (1.0 - tense); */
/* North side */
kase = x_case * 5 + 4;
ij = C->ij_nw_corner + i * C->my;
u[ij + C->offset[kase][0]] =
-(gmt_grdfloat)(-u[ij + C->offset[kase][11]] + C->eps_m2 * (u[ij + C->offset[kase][1]] + u[ij + C->offset[kase][3]]
- u[ij + C->offset[kase][8]] - u[ij + C->offset[kase][10]])
+ C->two_plus_em2 * (u[ij + C->offset[kase][9]] - u[ij + C->offset[kase][2]]) );
/* - tense * C->eps_m2 * (u[ij + C->offset[kase][2]] - u[ij + C->offset[kase][9]]) / (1.0 - tense); */
}
y_s_case = 0;
y_n_case = C->block_ny - 1;
for (j = 0; j < C->n_rows; j += C->grid, y_s_case++, y_n_case--) {
if(y_s_case < 2)
y_case = y_s_case;
else if(y_n_case < 2)
y_case = 4 - y_n_case;
else
y_case = 2;
if (C->periodic) { /* Set periodic boundary conditions in longitude */
/* West side */
kase = y_case;
ij_sw = C->ij_sw_corner + j;
ij_se = C->ij_se_corner + j;
u[ij_sw+C->offset[kase][4]] = u[ij_se+C->offset[20+kase][4]];
/* East side */
kase = 20 + y_case;
u[ij_se + C->offset[kase][7]] = u[ij_sw+C->offset[y_case][7]];
}
else { /* Natural BCs */
/* West side */
kase = y_case;
ij = C->ij_sw_corner + j;
u[ij+C->offset[kase][4]] =
u[ij + C->offset[kase][7]] + (gmt_grdfloat)(C->eps_p2 * (u[ij + C->offset[kase][3]] + u[ij + C->offset[kase][10]]
-u[ij + C->offset[kase][1]] - u[ij + C->offset[kase][8]])
+ C->two_plus_ep2 * (u[ij + C->offset[kase][5]] - u[ij + C->offset[kase][6]]));
/* + tense * (u[ij + C->offset[kase][6]] - u[ij + C->offset[kase][5]]) / (1.0 - tense); */
/* East side */
kase = 20 + y_case;
ij = C->ij_se_corner + j;
u[ij + C->offset[kase][7]] =
- (gmt_grdfloat)(-u[ij + C->offset[kase][4]] + C->eps_p2 * (u[ij + C->offset[kase][3]] + u[ij + C->offset[kase][10]]
- u[ij + C->offset[kase][1]] - u[ij + C->offset[kase][8]])
+ C->two_plus_ep2 * (u[ij + C->offset[kase][5]] - u[ij + C->offset[kase][6]]) );
/* - tense * (u[ij + C->offset[kase][6]] - u[ij + C->offset[kase][5]]) / (1.0 - tense); */
}
}
/* That's it for the boundary points. Now loop over all data */
x_w_case = 0;
x_e_case = C->block_nx - 1;
for (i = 0; i < C->n_columns; i += C->grid, x_w_case++, x_e_case--) {
if (x_w_case < 2)
x_case = x_w_case;
else if (x_e_case < 2)
x_case = 4 - x_e_case;
else
x_case = 2;
y_s_case = 0;
y_n_case = C->block_ny - 1;
ij = C->ij_sw_corner + i * C->my;
for (j = 0; j < C->n_rows; j += C->grid, ij += C->grid, y_s_case++, y_n_case--) {
if (iu[ij] == SURFACE_IS_CONSTRAINED) continue; /* Point is fixed, nothing to do */
if (y_s_case < 2)
y_case = y_s_case;
else if (y_n_case < 2)
y_case = 4 - y_n_case;
else
y_case = 2;
kase = x_case * 5 + y_case;
sum_ij = 0.0;
if (iu[ij] == SURFACE_IS_UNCONSTRAINED) { /* Point is unconstrained */
for (k = 0; k < 12; k++)
sum_ij += (u[ij + C->offset[kase][k]] * C->coeff[0][k]);
}