|
1 #====================================================================== |
|
2 # P L O T _ A C C E L E R A T I O N _ R E S I D U A L S . P L |
|
3 # doc: Tue May 17 21:40:08 2016 |
|
4 # dlm: Wed May 18 19:43:18 2016 |
|
5 # (c) 2016 A.M. Thurnherr |
|
6 # uE-Info: 46 37 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # HISTORY: |
|
10 # May 17, 2016: - created from [plot_attitude_biases.pl] |
|
11 # May 18, 2016: - made work |
|
12 |
|
13 require "$ANTS/libGMT.pl"; |
|
14 |
|
15 sub plot_acceleration_residuals($) |
|
16 { |
|
17 my($pfn) = @_; # plot file name |
|
18 |
|
19 my($xmin) = -1.5; |
|
20 my($xmax) = 1.5; |
|
21 my($ymin) = -0.03; |
|
22 my($ymax) = 0.03; |
|
23 my($bin_size) = 0.1; # m/s^3 |
|
24 |
|
25 my($min_thin) = 200; |
|
26 my($min_fat) = 500; |
|
27 my($btstrp_ndraw) = 20; |
|
28 my($excluded_surf_layer) = 150; |
|
29 my($hist_height) = 0.015; |
|
30 |
|
31 my($xo1) = -0.03; |
|
32 my($xo2) = -0.01; |
|
33 my($xo3) = 0.01; |
|
34 my($xo4) = 0.03; |
|
35 |
|
36 #------------------------------------------------------- |
|
37 # Bin-Average & Create Histogram from Beampair Residuals |
|
38 # - use 0.1 m/s^3 |
|
39 # - also calculate mean/rms pitch |
|
40 #------------------------------------------------------- |
|
41 |
|
42 my(@pHistDC,@rHistDC,@pSumDC,@rSumDC,$pHistDC,$rHistDC,$pSumDC,$rSumDC); |
|
43 my(@pHistUC,@rHistUC,@pSumUC,@rSumUC,$pHistUC,$rHistUC,$pSumUC,$rSumUC); |
|
44 my(@pValsDC,@rValsDC,@pValsUC,@rValsUC,$mode); |
|
45 my(@w_ttSumDC,@w_ttSumUC); |
|
46 for (my($e)=$firstGoodEns; $e<=$lastGoodEns; $e++) { |
|
47 next unless numberp($LADCP{ENSEMBLE}[$e]->{CTD_DEPTH}); |
|
48 my(@bindepth) = calc_binDepths($e); |
|
49 |
|
50 for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) { |
|
51 next if ($bindepth[$bin] <= $excluded_surf_layer); |
|
52 next unless ($bin+1>=$outGrid_firstBin && $bin+1<=$outGrid_lastBin); |
|
53 next unless numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin]) && |
|
54 numberp($LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin]); |
|
55 |
|
56 my($hi) = int(($CTD{W_tt}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}]-$xmin)/$bin_size); |
|
57 next unless ($hi >= 0); |
|
58 my($bi) = $bindepth[$bin]/$opt_o; |
|
59 |
|
60 if ($e < $LADCP_atbottom) { # DOWNCAST |
|
61 $w_ttSumDC += $CTD{W_tt}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}]; $pHistDC++; $rHistDC++; |
|
62 $pSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $DNCAST{MEDIAN_W}[$bi]; |
|
63 $rSumDC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $DNCAST{MEDIAN_W}[$bi]; |
|
64 push(@{$pValsDC[$hi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $DNCAST{MEDIAN_W}[$bi]); |
|
65 push(@{$rValsDC[$hi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $DNCAST{MEDIAN_W}[$bi]); |
|
66 $pSumDC[$hi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $DNCAST{MEDIAN_W}[$bi]; $pHistDC[$hi]++; |
|
67 $rSumDC[$hi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $DNCAST{MEDIAN_W}[$bi]; $rHistDC[$hi]++; |
|
68 $mode = $pHistDC[$hi] if ($pHistDC[$hi] > $mode); |
|
69 $mode = $rHistDC[$hi] if ($rHistDC[$hi] > $mode); |
|
70 } else { # UPCAST |
|
71 $w_ttSumUC += $CTD{W_tt}[$LADCP{ENSEMBLE}[$e]->{CTD_SCAN}]; $pHistUC++; $rHistUC++; |
|
72 $pSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $UPCAST{MEDIAN_W}[$bi]; |
|
73 $rSumUC += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $UPCAST{MEDIAN_W}[$bi]; |
|
74 push(@{$pValsUC[$hi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $UPCAST{MEDIAN_W}[$bi]); |
|
75 push(@{$rValsUC[$hi]},$LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $UPCAST{MEDIAN_W}[$bi]); |
|
76 $pSumUC[$hi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W12}[$bin] - $UPCAST{MEDIAN_W}[$bi]; $pHistUC[$hi]++; |
|
77 $rSumUC[$hi] += $LADCP{ENSEMBLE}[$e]->{SSCORRECTED_OCEAN_W34}[$bin] - $UPCAST{MEDIAN_W}[$bi]; $rHistUC[$hi]++; |
|
78 $mode = $pHistUC[$hi] if ($pHistUC[$hi] > $mode); |
|
79 $mode = $rHistUC[$hi] if ($rHistUC[$hi] > $mode); |
|
80 } |
|
81 } |
|
82 } |
|
83 |
|
84 #---------- |
|
85 # Plot Data |
|
86 #---------- |
|
87 |
|
88 my($R) = "-R$xmin/$xmax/$ymin/$ymax"; # plot frame |
|
89 GMT_begin($pfn,'-JX10/10',$R,"-P -Bg0.5f0.1a0.5:'Acceleration Derivative (w\@-tt\@-) [m s\@+-3\@+]':/g0.01f0.001a0.01:'Beam-Plane Residual Vertical Velocity [m/s]':WeSn"); |
|
90 |
|
91 # ZERO LINE |
|
92 GMT_psxy('-W4,CornflowerBlue'); |
|
93 print(GMT "$xmin 0\n$xmax 0\n"); |
|
94 |
|
95 # HISTOGRAMS |
|
96 GMT_psxy('-W2,coral'); |
|
97 for (my($i)=0; $i<@pHistDC; $i++) { |
|
98 next unless ($pHistDC[$i] >= $min_thin); |
|
99 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", |
|
100 $i*$bin_size+$xmin+$xo2-($bin_size/2),$ymin, |
|
101 $i*$bin_size+$xmin+$xo2-($bin_size/2),$ymin+$hist_height*$pHistDC[$i]/$mode, |
|
102 $i*$bin_size+$xmin+$xo2+($bin_size/2),$ymin+$hist_height*$pHistDC[$i]/$mode, |
|
103 $i*$bin_size+$xmin+$xo2+($bin_size/2),$ymin, |
|
104 $i*$bin_size+$xmin+$xo2+($bin_size/2),$ymin+$hist_height*$pHistDC[$i]/$mode); |
|
105 } |
|
106 GMT_psxy('-W2,SeaGreen'); |
|
107 for (my($i)=0; $i<@pHistUC; $i++) { |
|
108 next unless ($pHistUC[$i] >= $min_thin); |
|
109 printf(GMT ">\n%f %f\n%f %f\n%f %f\n>\n%f %f\n%f %f\n", |
|
110 $i*$bin_size+$xmin+$xo3-($bin_size/2),$ymin, |
|
111 $i*$bin_size+$xmin+$xo3-($bin_size/2),$ymin+$hist_height*$pHistUC[$i]/$mode, |
|
112 $i*$bin_size+$xmin+$xo3+($bin_size/2),$ymin+$hist_height*$pHistUC[$i]/$mode, |
|
113 $i*$bin_size+$xmin+$xo3+($bin_size/2),$ymin, |
|
114 $i*$bin_size+$xmin+$xo3+($bin_size/2),$ymin+$hist_height*$pHistUC[$i]/$mode); |
|
115 } |
|
116 |
|
117 # DC PITCH |
|
118 GMT_psxy('-Ey0.2/2,coral'); |
|
119 for (my($i)=0; $i<@pHistDC; $i++) { |
|
120 next unless ($pHistDC[$i] >= $min_fat); |
|
121 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits |
|
122 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
123 } |
|
124 GMT_psxy('-Ey0.2/1,coral'); # dc pitch |
|
125 for (my($i)=0; $i<@pHistDC; $i++) { |
|
126 next unless ($pHistDC[$i]>=$min_thin && $pHistDC[$i]<$min_fat); |
|
127 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsDC[$i]}); # 95% bootstrap conf limits |
|
128 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo1,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
129 } |
|
130 GMT_psxy('-Sc0.15 -Gcoral'); |
|
131 for (my($i)=0; $i<@pHistDC; $i++) { |
|
132 next unless ($pHistDC[$i] >= $min_thin); |
|
133 printf(GMT "%f %f\n",$i*$bin_size+$xmin+$xo1,$pSumDC[$i]/$pHistDC[$i]); # errorbar center symbol |
|
134 } |
|
135 if ($pHistDC) { |
|
136 GMT_psxy('-W1,coral,8_2:0'); # averages (lines) |
|
137 printf(GMT ">\n%f %f\n%f %f\n",$xmin,$pSumDC/$pHistDC,$xmax,$pSumDC/$pHistDC); # bias |
|
138 } |
|
139 |
|
140 # DC ROLL |
|
141 GMT_psxy('-Ey0.2/2,coral'); |
|
142 for (my($i)=0; $i<@rHistDC; $i++) { |
|
143 next unless ($rHistDC[$i] >= $min_fat); |
|
144 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]}); |
|
145 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo2,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
146 } |
|
147 GMT_psxy('-Ey0.2/1,coral'); |
|
148 for (my($i)=0; $i<@rHistDC; $i++) { |
|
149 next unless ($rHistDC[$i]>=$min_thin && $rHistDC[$i]<$min_fat); |
|
150 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsDC[$i]}); |
|
151 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo2,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
152 } |
|
153 GMT_psxy('-Sx0.25 -W2,coral'); |
|
154 for (my($i)=0; $i<@rHistDC; $i++) { |
|
155 next unless ($rHistDC[$i] >= $min_thin); |
|
156 printf(GMT "%f %f\n",$i*$bin_size+$xmin+$xo2,$rSumDC[$i]/$rHistDC[$i]); |
|
157 } |
|
158 if ($rHistDC) { |
|
159 GMT_psxy('-W1,coral,2_2:0'); |
|
160 printf(GMT ">\n%f %f\n%f %f\n",$xmin,$rSumDC/$rHistDC,$xmax,$rSumDC/$rHistDC); |
|
161 } |
|
162 |
|
163 # UC PITCH |
|
164 GMT_psxy('-Ey0.2/2,SeaGreen'); |
|
165 for (my($i)=0; $i<@pHistUC; $i++) { |
|
166 next unless ($pHistUC[$i] >= $min_fat); |
|
167 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]}); |
|
168 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
169 } |
|
170 GMT_psxy('-Ey0.2/1,SeaGreen'); |
|
171 for (my($i)=0; $i<@pHistUC; $i++) { |
|
172 next unless ($pHistUC[$i]>=$min_thin && $pHistUC[$i]<$min_fat); |
|
173 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$pValsUC[$i]}); |
|
174 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo3,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
175 } |
|
176 GMT_psxy('-Sc0.15 -GSeaGreen'); |
|
177 for (my($i)=0; $i<@pHistUC; $i++) { |
|
178 next unless ($pHistUC[$i] >= $min_thin); |
|
179 printf(GMT "%f %f\n",$i*$bin_size+$xmin+$xo3,$pSumUC[$i]/$pHistUC[$i]); |
|
180 } |
|
181 if ($pHistUC) { |
|
182 GMT_psxy('-W1,SeaGreen,8_2:0'); |
|
183 printf(GMT ">\n%f %f\n%f %f\n",$xmin,$pSumUC/$pHistUC,$xmax,$pSumUC/$pHistUC); |
|
184 } |
|
185 |
|
186 # UC ROLL |
|
187 GMT_psxy('-Ey0.2/2,SeaGreen'); |
|
188 for (my($i)=0; $i<@rHistUC; $i++) { |
|
189 next unless ($rHistUC[$i] >= $min_fat); |
|
190 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]}); |
|
191 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo4,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
192 } |
|
193 GMT_psxy('-Ey0.2/1,SeaGreen'); |
|
194 for (my($i)=0; $i<@rHistUC; $i++) { |
|
195 next unless ($rHistUC[$i]>=$min_thin && $rHistUC[$i]<$min_fat); |
|
196 my($minLim,$maxLim) = &bootstrap($btstrp_ndraw,0.95,\&avg,@{$rValsUC[$i]}); |
|
197 printf(GMT "%f %f %f\n",$i*$bin_size+$xmin+$xo4,($maxLim+$minLim)/2,($maxLim-$minLim)/2); |
|
198 } |
|
199 GMT_psxy('-Sx0.25 -W2,SeaGreen'); |
|
200 for (my($i)=0; $i<@rHistUC; $i++) { |
|
201 next unless ($rHistUC[$i] >= $min_thin); |
|
202 printf(GMT "%f %f\n",$i*$bin_size+$xmin+$xo4,$rSumUC[$i]/$rHistUC[$i]); |
|
203 } |
|
204 if ($rHistUC) { |
|
205 GMT_psxy('-W1,SeaGreen,2_2:0'); |
|
206 printf(GMT ">\n%f %f\n%f %f\n",$xmin,$rSumUC/$rHistUC,$xmax,$rSumUC/$rHistUC); |
|
207 } |
|
208 |
|
209 # ANNOTATIONS |
|
210 GMT_unitcoords(); # LABELS |
|
211 GMT_pstext('-F+f9,Helvetica,orange+jTR -N -Gwhite'); |
|
212 print(GMT "0.99 0.99 V$VERSION\n"); |
|
213 GMT_pstext('-F+f14,Helvetica,blue+jBL -N'); # profile id |
|
214 print(GMT "0.0 1.03 $P{out_basename} $P{run_label}\n"); |
|
215 |
|
216 GMT_setR($R); # FINISH PLOT |
|
217 GMT_end(); |
|
218 } |
|
219 |
|
220 1; # return true on require |