1#!/usr/bin/env bash
2#
3# Test script that exercise the various options in x2sys.
4# We generate a grid and some fake tracks and sample the
5# grid, then add various systematic errors to the tracks
6# and finally try to solve for and undo the errors.
7
8# 1. Make a gmt surface grid with a Mexican hat bump in the middle
9
10delete=0	# Set to 0 for debug where files are not removed
11
12gmt grdmath -R-4/4/-4/4 -I0.1 0 0 CDIST DUP DUP MUL NEG 4 DIV EXP EXCH 3 MUL COS MUL = hat.nc
13
14# 2. Create 3 fake tracks
15cat << EOF | gmt mapproject -G+uc | gmt sample1d -Fl -T2 -I0.1 | gmt grdtrack -Ghat.nc > trackA.xydz
16-3	-3
173	3
183	1
19-2	1
20EOF
21cat << EOF | gmt mapproject -G+uc | gmt sample1d -Fl -T2 -I0.1 | gmt grdtrack -Ghat.nc > trackB.xydz
223	-1
230	2
240	-2
25-2	0
26EOF
27# THrow in a wrench by scaling the grid by 1.1 before sampling track C:
28gmt grdmath hat.nc 1.1 MUL = $$.nc
29cat << EOF | gmt mapproject -G+uc | gmt sample1d -Fl -T2 -I0.1 | gmt grdtrack -G$$.nc > trackC.xydz
30-3	-1
312	-1
322	3
33-2	-1
34EOF
35rm -f $$.nc
36
37# Create a *.def file and initialize the tag
38
39cat << EOF > xydz.def
40# This file applies to plain 4-column ASCII files with no header record.
41#
42#---------------------------------------------------------------------
43#ASCII		# The input file is ASCII
44#SKIP 0		# The number of header records to skip
45#---------------------------------------------------------------------
46#name	intype	NaN-proxy?	NaN-proxy	scale	offset	oformat
47x	a	N		0		1	0	-
48y	a	N		0		1	0	-
49d	a	N		0		1	0	-
50z	a	N		0		1	0	-
51EOF
52
53gmt x2sys_init FAKE -Dxydz -V -F -R-5/5/-5/5 -I1 -Cc
54rm -f xydz.def
55
56gmt x2sys_cross -TFAKE track[ABC].xydz -Qe -V -2 > fake_COE_orig.txt
57
58gmt makecpt -T-1/1 -Cpolar -Z > hat.cpt
59PS=x2sys_1.ps
60gmt grdview hat.nc -Chat.cpt -JX6i -P -B1g1 -BWSne -Qs -Wc1p -K -X1.25i > $PS
61gmt psxy -R -J -O -K trackA.xydz -W2p,red >> $PS
62gmt psxy -R -J -O -K trackB.xydz -W2p,blue >> $PS
63gmt psxy -R -J -O -K trackC.xydz -W2p,black >> $PS
64gmt psxy -R -J -O -K fake_COE_orig.txt -Sc0.15 -W0.5p >> $PS
65head -1 trackA.xydz | gmt psxy -R -J -O -K -Sc0.1i -Gred >> $PS
66head -1 trackB.xydz | gmt psxy -R -J -O -K -Sc0.1i -Gblue >> $PS
67head -1 trackC.xydz | gmt psxy -R -J -O -K -Sc0.1i -Gblack >> $PS
68cut -f3,4 trackA.xydz | gmt psxy -R0/16/-1/1 -JX6i/2i -Bx5f1 -By0.2 -BWSne+t"X2SYS 1\072 3 tracks" -O -K -Y6.5i -W2p,red >> $PS
69cut -f3,4 trackB.xydz | gmt psxy -R -J -O -K -W2p,blue >> $PS
70cut -f3,4 trackC.xydz | gmt psxy -R -J -O -W2p,black >> $PS
71gv $PS &
72
73#---------------------------------------------------------
74# Test 1: Add some constant shifts and try to resolve them
75#---------------------------------------------------------
76
77gmt math trackA.xydz -C3 0.2 ADD -Ca = trackAc.xydz
78gmt math trackB.xydz -C3 0.1 SUB -Ca = trackBc.xydz
79gmt math trackC.xydz -C3 0.05 ADD -Ca = trackCc.xydz
80
81# Obtain all external crossovers
82
83gmt x2sys_cross -TFAKE track[ABC]c.xydz -Qe -V -2 > fake_COE_constant.txt
84gmt x2sys_list -TFAKE -Cz fake_COE_constant.txt -Fnc > COE.txt
85
86# Make track plots for shifted tracks
87
88PS=x2sys_2.ps
89cut -f3,4 trackA.xydz | gmt psxy -R0/16/-1.5/1.5 -JX6i/2.25i -X1.25i -P -Y7.25i -Bx5f1 -By0.2g10 -BWSne+t"X2SYS 2\072 Constants added" -K -W0.25p,red,- > $PS
90cut -f3,4 trackAc.xydz | gmt psxy -R -J -O -K -W2p,red >> $PS
91cut -f3,4 trackB.xydz | gmt psxy -R -J -O -K -W0.25p,blue,- >> $PS
92cut -f3,4 trackBc.xydz | gmt psxy -R -J -O -K -W2p,blue >> $PS
93cut -f3,4 trackC.xydz | gmt psxy -R -J -O -K -W0.25p,black,- >> $PS
94cut -f3,4 trackCc.xydz | gmt psxy -R -J -O -K -W2p,black >> $PS
95
96cut -f3,4 trackAc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,red >> $PS
97gmt x2sys_list -TFAKE -Cz fake_COE_constant.txt -Fdc -StrackA | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
98cut -f3,4 trackBc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,blue >> $PS
99gmt x2sys_list -TFAKE -Cz fake_COE_constant.txt -Fdc -StrackB | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
100cut -f3,4 trackCc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,black >> $PS
101gmt x2sys_list -TFAKE -Cz fake_COE_constant.txt -Fdc -StrackC | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
102gmt psxy -R -J -O -T >> $PS
103gv $PS &
104
105# Solve for constants
106gmt x2sys_solve COE.txt -TFAKE -Cz -Ec -V > $X2SYS_HOME/FAKE/corr_const.lis
107A=$(grep trackA corr.lis | cut -f3)
108B=$(grep trackB corr.lis | cut -f3)
109C=$(grep trackC corr.lis | cut -f3)
110
111# Correct tracks
112gmt x2sys_datalist -TFAKE -Lcorr_const.lis trackAc.xydz > trackAcc.xydz
113gmt x2sys_datalist -TFAKE -Lcorr_const.lis trackBc.xydz > trackBcc.xydz
114gmt x2sys_datalist -TFAKE -Lcorr_const.lis trackCc.xydz > trackCcc.xydz
115
116# Make track plots for corrected tracks
117gmt x2sys_cross -TFAKE track[ABC]cc.xydz -Qe -V -2 > fake_COE_constant_corr.txt
118
119PS=x2sys_3.ps
120cut -f3,4 trackA.xydz | gmt psxy -R0/16/-1.5/1.5 -JX6i/2.25i -X1.25i -P -Y7.25i -Bx5f1 -By0.2g10 -BWSne+t"X2SYS 3\072 Constants resolved" -K -W0.25p,red,- > $PS
121cut -f3,4 trackAcc.xydz | gmt psxy -R -J -O -K -W2p,red >> $PS
122cut -f3,4 trackB.xydz | gmt psxy -R -J -O -K -W0.25p,blue,- >> $PS
123cut -f3,4 trackBcc.xydz | gmt psxy -R -J -O -K -W2p,blue >> $PS
124cut -f3,4 trackC.xydz | gmt psxy -R -J -O -K -W0.25p,black,- >> $PS
125cut -f3,4 trackCcc.xydz | gmt psxy -R -J -O -K -W2p,black >> $PS
126
127cut -f3,4 trackAcc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,red >> $PS
128gmt x2sys_list -TFAKE -Cz fake_COE_constant_corr.txt -Fdc -StrackA | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
129cut -f3,4 trackBcc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,blue >> $PS
130gmt x2sys_list -TFAKE -Cz fake_COE_constant_corr.txt -Fdc -StrackB | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
131cut -f3,4 trackCcc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,black >> $PS
132gmt x2sys_list -TFAKE -Cz fake_COE_constant_corr.txt -Fdc -StrackC | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
133gmt psxy -R -J -O -T >> $PS
134gv $PS &
135#---------------------------------------------------------
136# Test 2: Add some linear drifts and try to resolve them
137#---------------------------------------------------------
138
139awk '{printf "%s\t%s\t%s\t%g\n", $1, $2, $3, $4 + 0.2 + 0.04*$3}' trackA.xydz > trackAd.xydz
140awk '{printf "%s\t%s\t%s\t%g\n", $1, $2, $3, $4 - 0.1 + 0.03*$3}' trackB.xydz > trackBd.xydz
141awk '{printf "%s\t%s\t%s\t%g\n", $1, $2, $3, $4 + 0.05 - 0.025*$3}' trackC.xydz > trackCd.xydz
142
143# Obtain all external crossovers
144
145gmt x2sys_cross -TFAKE track[ABC]d.xydz -Qe -V -2 > fake_COE_drift.txt
146gmt x2sys_list -TFAKE -Cz fake_COE_drift.txt -Fndc > COE.txt
147
148# Make track plots for shifted tracks
149
150PS=x2sys_4.ps
151cut -f3,4 trackA.xydz | gmt psxy -R0/16/-1.5/1.5 -JX6i/2.25i -X1.25i -P -Y7.25i -Bx5f1 -By0.2g10 -BWSne+t"X2SYS 4\072 Drifts added" -K -W0.25p,red,- > $PS
152cut -f3,4 trackAd.xydz | gmt psxy -R -J -O -K -W2p,red >> $PS
153cut -f3,4 trackB.xydz | gmt psxy -R -J -O -K -W0.25p,blue,- >> $PS
154cut -f3,4 trackBd.xydz | gmt psxy -R -J -O -K -W2p,blue >> $PS
155cut -f3,4 trackC.xydz | gmt psxy -R -J -O -K -W0.25p,black,- >> $PS
156cut -f3,4 trackCd.xydz | gmt psxy -R -J -O -K -W2p,black >> $PS
157
158cut -f3,4 trackAd.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,red >> $PS
159gmt x2sys_list -TFAKE -Cz fake_COE_drift.txt -Fdc -StrackA | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
160cut -f3,4 trackBd.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,blue >> $PS
161gmt x2sys_list -TFAKE -Cz fake_COE_drift.txt -Fdc -StrackB | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
162cut -f3,4 trackCd.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,black >> $PS
163gmt x2sys_list -TFAKE -Cz fake_COE_drift.txt -Fdc -StrackC | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
164gmt psxy -R -J -O -T >> $PS
165gv $PS &
166
167# Solve for trends
168gmt x2sys_solve COE.txt -TFAKE -Cz -Ed -V > $X2SYS_HOME/FAKE/corr_trend.lis
169
170# Correct tracks
171gmt x2sys_datalist -TFAKE -Lcorr_trend.lis trackAd.xydz > trackAdc.xydz
172gmt x2sys_datalist -TFAKE -Lcorr_trend.lis trackBd.xydz > trackBdc.xydz
173gmt x2sys_datalist -TFAKE -Lcorr_trend.lis trackCd.xydz > trackCdc.xydz
174
175# Make track plots for corrected tracks
176gmt x2sys_cross -TFAKE track[ABC]dc.xydz -Qe -V -2 > fake_COE_drift_corr.txt
177
178PS=x2sys_5.ps
179cut -f3,4 trackA.xydz | gmt psxy -R0/16/-1.5/1.5 -JX6i/2.25i -X1.25i -P -Y7.25i -Bx5f1 -By0.2g10 -BWSne+t"X2SYS 5\072 Drifts resolved" -K -W0.25p,red,- > $PS
180cut -f3,4 trackAdc.xydz | gmt psxy -R -J -O -K -W2p,red >> $PS
181cut -f3,4 trackB.xydz | gmt psxy -R -J -O -K -W0.25p,blue,- >> $PS
182cut -f3,4 trackBdc.xydz | gmt psxy -R -J -O -K -W2p,blue >> $PS
183cut -f3,4 trackC.xydz | gmt psxy -R -J -O -K -W0.25p,black,- >> $PS
184cut -f3,4 trackCdc.xydz | gmt psxy -R -J -O -K -W2p,black >> $PS
185
186cut -f3,4 trackAdc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,red >> $PS
187gmt x2sys_list -TFAKE -Cz fake_COE_drift_corr.txt -Fdc -StrackA | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
188cut -f3,4 trackBdc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,blue >> $PS
189gmt x2sys_list -TFAKE -Cz fake_COE_drift_corr.txt -Fdc -StrackB | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
190cut -f3,4 trackCdc.xydz | gmt psxy -R -J -O -K -Bx5f1 -By0.2g10 -BWSne -Y-2.25i -W1p,black >> $PS
191gmt x2sys_list -TFAKE -Cz fake_COE_drift_corr.txt -Fdc -StrackC | gmt psxy -R -J -O -K -Sc0.05 -Gblack >> $PS
192gmt psxy -R -J -O -T >> $PS
193gv $PS &
194if [ $delete -eq 1 ]; then
195	rm -f hat.nc hat.cpt track[ABC]*.xydz fake_COE_*.txt COE.txt corr_const.lis corr_trend.lis
196fi
197