1 #====================================================================== |
|
2 # P L O T _ A T T I T U D E _ B I A S E S _ W . P L |
|
3 # doc: Sun May 15 16:08:59 2016 |
|
4 # dlm: Tue May 24 16:38:31 2016 |
|
5 # (c) 2016 A.M. Thurnherr |
|
6 # uE-Info: 15 31 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # HISTORY: |
|
10 # May 15, 2016: - created from [plot_mean_residuals.pl] |
|
11 # May 16, 2016: - continued |
|
12 # May 17, 2016: - renamed from [plot_attitude_biases.pl] |
|
13 # May 18, 2016: - added version |
|
14 # - expunged $realLastGoodEns |
|
15 # May 19, 2016: - added notes about wrong beam plane |
|
16 # May 24, 2016: - calc_binDepths() -> binDepths() |
|
17 |
|
18 # IMPORTANT NOTE: |
|
19 # - the variables prefixed with p/r refer to beam-pairs 1,2 and 3,4 respectively, |
|
20 # i.e. the p variables correspond to the roll plane and the r variables |
|
21 # correspond to the pitch plane |
|
22 |
|
23 require "$ANTS/libGMT.pl"; |
|
24 |
|
25 sub plot_attitude_biases_w($) |
|
26 { |
|
27 my($pfn) = @_; # plot file name |
|
28 |
|
29 my($xmin) = -round($opt_t); # full pitch range |
|
30 my($xmax) = round($opt_t); |
|
31 my($ymin) = -0.05; |
|
32 my($ymax) = 0.05; |
|
33 |
|
34 my($min_thin) = 200; |
|
35 my($min_fat) = 500; |
|
36 my($btstrp_ndraw) = 20; |
|
37 my($excluded_surf_layer) = 150; |
|
38 |
|
39 #-------------------------------------------------------- |
|
40 # Bin-Average & Create Histogram from Beampair Velocities |
|
41 # - use 1-degree bins for simplicity |
|
42 # - use gimbal pitch (not measured pitch) |
|
43 #-------------------------------------------------------- |
|
44 |
|
45 my(@pHistDC,@rHistDC,@pSumDC,@rSumDC,$pHistDC,$rHistDC,$pSumDC,$rSumDC); |
|
46 my(@pHistUC,@rHistUC,@pSumUC,@rSumUC,$pHistUC,$rHistUC,$pSumUC,$rSumUC); |
|
47 my(@pValsDC,@rValsDC,@pValsUC,@rValsUC,$mode); |
|
48 for (my($e)=$firstGoodEns; $e<=$lastGoodEns; $e++) { |
|
49 next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}); |
|
50 my(@bindepth) = binDepths($e); |
|
51 for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { |
|
52 next if ($bindepth[$bin] <= $excluded_surf_layer); |
|
53 next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin); |
|
54 next unless numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]) && |
|
55 numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]); |
|
56 my($pi) = int($LADCP{ENSEMBLE}[$e]->{GIMBAL_PITCH}+$opt_t); # pitch/roll indices |
|
57 my($ri) = int($LADCP{ENSEMBLE}[$e]->{ROLL}+$opt_t); |
|
58 if ($e < $LADCP_atbottom) { # downcast |
|
59 push(@{$pValsDC[$pi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]); |
|
60 push(@{$rValsDC[$ri]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]); |
|
61 $pSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistDC++; |
|
62 $rSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistDC++; |
|
63 $pSumDC[$pi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistDC[$pi]++; |
|
64 $rSumDC[$ri] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistDC[$ri]++; |
|
65 $mode = $pHistDC[$pi] if ($pHistDC[$pi] > $mode); |
|
66 $mode = $rHistDC[$ri] if ($rHistDC[$ri] > $mode); |
|
67 } else { # upcast |
|
68 push(@{$pValsUC[$pi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]); |
|
69 push(@{$rValsUC[$ri]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]); |
|
70 $pSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistUC++; |
|
71 $rSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistUC++; |
|
72 $pSumUC[$pi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]; $pHistUC[$pi]++; |
|
73 $rSumUC[$ri] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]; $rHistUC[$ri]++; |
|
74 $mode = $pHistUC[$pi] if ($pHistUC[$pi] > $mode); |
|
75 $mode = $rHistUC[$ri] if ($rHistUC[$ri] > $mode); |
|
76 } |
|
77 } |
|
78 } |
|
79 |
|
80 #---------- |
|
81 # Plot Data |
|
82 #---------- |
|
83 |
|
84 my($R) = "-R$xmin/$xmax/$ymin/$ymax"; # begin plot |
|
85 GMT_begin($pfn,'-JX10/10',$R,"-P -Bg5f1a5:'Pitch/Roll [\260]':/g0.01f0.001a0.01:'Beam-Plane Vertical Velocity [m/s]':WeSn"); |
|
86 |
|
87 # ZERO LINE |
|
88 GMT_psxy('-W4,CornflowerBlue'); |
|
89 print(GMT "$xmin 0\n$xmax 0\n"); |
|
90 |
|
91 # DC BEAMS 1,2 |
|
92 GMT_psxy('-Ey0.2/2,coral'); |
|
93 for (my($i)=0; $i<2*round($opt_t); $i++) { # error bars |
|
94 next unless ($pHistDC[$i] >= $min_fat); |
|
95 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits |
|
96 printf(GMT "%f %f %f\n",$i-round($opt_t)-0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
97 } |
|
98 GMT_psxy('-Ey0.2/1,coral'); |
|
99 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
100 next unless ($pHistDC[$i]>=$min_thin && $pHistDC[$i]<$min_fat); |
|
101 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits |
|
102 printf(GMT "%f %f %f\n",$i-round($opt_t)-0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
103 } |
|
104 GMT_psxy('-Sc0.15 -Gcoral'); |
|
105 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
106 next unless ($pHistDC[$i] >= $min_thin); |
|
107 printf(GMT "%f %f\n",$i-round($opt_t)-0.3,$pSumDC[$i]/$pHistDC[$i]); # errorbar center symbol |
|
108 printf(GMT "%f %f\n",$i-round($opt_t)-0.3,$ymin+0.03*$pHistDC[$i]/$mode); # histogram symbol |
|
109 } |
|
110 if ($pHistDC) { |
|
111 GMT_psxy('-W1,coral,8_2:0'); # average bias (horizontal line); |
|
112 printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$pSumDC/$pHistDC,$opt_t,$pSumDC/$pHistDC) |
|
113 } |
|
114 GMT_psxy('-W2,coral,8_2:0'); |
|
115 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
116 next unless ($pHistDC[$i] >= $min_thin); |
|
117 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", # histogram bar |
|
118 $i-round($opt_t)-0.3-0.5,$ymin, |
|
119 $i-round($opt_t)-0.3-0.5,$ymin+0.03*$pHistDC[$i]/$mode, |
|
120 $i-round($opt_t)-0.3+0.5,$ymin+0.03*$pHistDC[$i]/$mode, |
|
121 $i-round($opt_t)-0.3+0.5,$ymin, |
|
122 $i-round($opt_t)-0.3+0.5,$ymin+0.03*$pHistDC[$i]/$mode); |
|
123 } |
|
124 |
|
125 # DC BEAMS 3,4 |
|
126 GMT_psxy('-Ey0.2/2,coral'); |
|
127 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
128 next unless ($rHistDC[$i] >= $min_fat); |
|
129 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]}); |
|
130 printf(GMT "%f %f %f\n",$i-round($opt_t)-0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
131 } |
|
132 GMT_psxy('-Ey0.2/1,coral'); |
|
133 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
134 next unless ($rHistDC[$i]>=$min_thin && $rHistDC[$i]<$min_fat); |
|
135 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]}); |
|
136 printf(GMT "%f %f %f\n",$i-round($opt_t)-0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
137 } |
|
138 GMT_psxy('-Sx0.25 -W2,coral'); |
|
139 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
140 next unless ($rHistDC[$i] >= $min_thin); |
|
141 printf(GMT "%f %f\n",$i-round($opt_t)-0.1,$rSumDC[$i]/$rHistDC[$i]); |
|
142 printf(GMT "%f %f\n",$i-round($opt_t)-0.1,$ymin+0.03*$rHistDC[$i]/$mode); |
|
143 } |
|
144 if ($rHistDC) { |
|
145 GMT_psxy('-W1,coral,2_2:0'); |
|
146 printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$rSumDC/$rHistDC,$opt_t,$rSumDC/$rHistDC); |
|
147 } |
|
148 GMT_psxy('-W2,coral,2_2:0'); |
|
149 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
150 next unless ($rHistDC[$i] >= $min_thin); |
|
151 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", |
|
152 $i-round($opt_t)-0.1-0.5,$ymin, |
|
153 $i-round($opt_t)-0.1-0.5,$ymin+0.03*$rHistDC[$i]/$mode, |
|
154 $i-round($opt_t)-0.1+0.5,$ymin+0.03*$rHistDC[$i]/$mode, |
|
155 $i-round($opt_t)-0.1+0.5,$ymin, |
|
156 $i-round($opt_t)-0.1+0.5,$ymin+0.03*$rHistDC[$i]/$mode); |
|
157 } |
|
158 |
|
159 # UC BEAMS 1,2 |
|
160 GMT_psxy('-Ey0.2/2,SeaGreen'); |
|
161 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
162 next unless ($pHistUC[$i] >= $min_fat); |
|
163 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]}); |
|
164 printf(GMT "%f %f %f\n",$i-round($opt_t)+0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
165 } |
|
166 GMT_psxy('-Ey0.2/1,SeaGreen'); |
|
167 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
168 next unless ($pHistUC[$i]>=$min_thin && $pHistUC[$i]<$min_fat); |
|
169 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]}); |
|
170 printf(GMT "%f %f %f\n",$i-round($opt_t)+0.1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
171 } |
|
172 GMT_psxy('-Sc0.15 -GSeaGreen'); |
|
173 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
174 next unless ($pHistUC[$i] >= $min_thin); |
|
175 printf(GMT "%f %f\n",$i-round($opt_t)+0.1,$pSumUC[$i]/$pHistUC[$i]); |
|
176 printf(GMT "%f %f\n",$i-round($opt_t)+0.1,$ymin+0.03*$pHistUC[$i]/$mode); |
|
177 } |
|
178 if ($pHistUC) { |
|
179 GMT_psxy('-W1,SeaGreen,8_2:0'); |
|
180 printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$pSumUC/$pHistUC,$opt_t,$pSumUC/$pHistUC); |
|
181 } |
|
182 GMT_psxy('-W2,SeaGreen,8_2:0'); |
|
183 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
184 next unless ($pHistUC[$i] >= $min_thin); |
|
185 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", |
|
186 $i-round($opt_t)+0.1-0.5,$ymin, |
|
187 $i-round($opt_t)+0.1-0.5,$ymin+0.03*$pHistUC[$i]/$mode, |
|
188 $i-round($opt_t)+0.1+0.5,$ymin+0.03*$pHistUC[$i]/$mode, |
|
189 $i-round($opt_t)+0.1+0.5,$ymin, |
|
190 $i-round($opt_t)+0.1+0.5,$ymin+0.03*$pHistUC[$i]/$mode); |
|
191 } |
|
192 |
|
193 # UC BEAMS 3,4 |
|
194 GMT_psxy('-Ey0.2/2,SeaGreen'); |
|
195 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
196 next unless ($rHistUC[$i] >= $min_fat); |
|
197 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]}); |
|
198 printf(GMT "%f %f %f\n",$i-round($opt_t)+0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
199 } |
|
200 GMT_psxy('-Ey0.2/1,SeaGreen'); |
|
201 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
202 next unless ($rHistUC[$i]>=$min_thin && $rHistUC[$i]<$min_fat); |
|
203 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]}); |
|
204 printf(GMT "%f %f %f\n",$i-round($opt_t)+0.3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
205 } |
|
206 GMT_psxy('-Sx0.25 -W2,SeaGreen'); |
|
207 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
208 next unless ($rHistUC[$i] >= $min_thin); |
|
209 printf(GMT "%f %f\n",$i-round($opt_t)+0.3,$rSumUC[$i]/$rHistUC[$i]); |
|
210 printf(GMT "%f %f\n",$i-round($opt_t)+0.3,$ymin+0.03*$rHistUC[$i]/$mode); |
|
211 } |
|
212 if ($rHistUC) { |
|
213 GMT_psxy('-W1,SeaGreen,2_2:0'); |
|
214 printf(GMT ">\n%f %f\n%f %f\n",-$opt_t,$rSumUC/$rHistUC,$opt_t,$rSumUC/$rHistUC); |
|
215 } |
|
216 GMT_psxy('-W2,SeaGreen,2_2:0'); |
|
217 for (my($i)=0; $i<2*round($opt_t); $i++) { |
|
218 next unless ($rHistUC[$i] >= $min_thin); |
|
219 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", |
|
220 $i-round($opt_t)+0.3-0.5,$ymin, |
|
221 $i-round($opt_t)+0.3-0.5,$ymin+0.03*$rHistUC[$i]/$mode, |
|
222 $i-round($opt_t)+0.3+0.5,$ymin+0.03*$rHistUC[$i]/$mode, |
|
223 $i-round($opt_t)+0.3+0.5,$ymin, |
|
224 $i-round($opt_t)+0.3+0.5,$ymin+0.03*$rHistUC[$i]/$mode); |
|
225 } |
|
226 |
|
227 GMT_unitcoords(); # LABELS |
|
228 GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite'); |
|
229 print(GMT "0.99 0.01 V$VERSION\n"); |
|
230 |
|
231 GMT_pstext('-F+f14,Helvetica,blue+jBL -N'); # profile id |
|
232 print(GMT "0.0 1.03 $P{out_basename} $P{run_label}\n"); |
|
233 |
|
234 GMT_setR($R); # FINISH PLOT |
|
235 GMT_end(); |
|
236 } |
|
237 |
|
238 1; # return true on require |
|