Mercurial > dive4elements > river
comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 1939:2730d17df021
Added the implementation of the 'Bezugslinienverfahren'.
flys-artifacts/trunk@3320 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 25 Nov 2011 18:52:11 +0000 |
parents | 1d991c91285b |
children | 2898b1ff6013 |
comparison
equal
deleted
inserted
replaced
1938:1d991c91285b | 1939:2730d17df021 |
---|---|
87 this.weight = weight; | 87 this.weight = weight; |
88 return this; | 88 return this; |
89 } | 89 } |
90 | 90 |
91 } // class Position | 91 } // class Position |
92 | |
93 public static final class SplineFunction { | |
94 | |
95 public PolynomialSplineFunction spline; | |
96 public double [] splineQs; | |
97 public double [] splineWs; | |
98 | |
99 public SplineFunction( | |
100 PolynomialSplineFunction spline, | |
101 double [] splineQs, | |
102 double [] splineWs | |
103 ) { | |
104 this.spline = spline; | |
105 this.splineQs = splineQs; | |
106 this.splineWs = splineWs; | |
107 } | |
108 | |
109 public double [][] sample( | |
110 int numSamples, | |
111 double km, | |
112 Calculation errors | |
113 ) { | |
114 double minQ = getQMin(); | |
115 double maxQ = getQMax(); | |
116 | |
117 double [] outWs = new double[numSamples]; | |
118 double [] outQs = new double[numSamples]; | |
119 | |
120 Arrays.fill(outWs, Double.NaN); | |
121 Arrays.fill(outQs, Double.NaN); | |
122 | |
123 double stepWidth = (maxQ - minQ)/numSamples; | |
124 | |
125 try { | |
126 double q = minQ; | |
127 for (int i = 0; i < outWs.length; ++i, q += stepWidth) { | |
128 outWs[i] = spline.value(outQs[i] = q); | |
129 } | |
130 } | |
131 catch (ArgumentOutsideDomainException aode) { | |
132 if (errors != null) { | |
133 // TODO: I18N | |
134 errors.addProblem(km, "spline interpolation failed"); | |
135 } | |
136 log.error("spline interpolation failed.", aode); | |
137 } | |
138 | |
139 return new double [][] { outWs, outQs }; | |
140 } | |
141 | |
142 public double getQMin() { | |
143 return Math.min(splineQs[0], splineQs[splineQs.length-1]); | |
144 } | |
145 | |
146 public double getQMax() { | |
147 return Math.max(splineQs[0], splineQs[splineQs.length-1]); | |
148 } | |
149 | |
150 /** Constructs a continues index between the columns to Qs. */ | |
151 public PolynomialSplineFunction createIndexQRelation() { | |
152 | |
153 double [] indices = new double[splineQs.length]; | |
154 for (int i = 0; i < indices.length; ++i) { | |
155 indices[i] = i; | |
156 } | |
157 | |
158 try { | |
159 SplineInterpolator interpolator = new SplineInterpolator(); | |
160 return interpolator.interpolate(indices, splineQs); | |
161 } | |
162 catch (MathIllegalArgumentException miae) { | |
163 // Ignore me! | |
164 } | |
165 return null; | |
166 } | |
167 } // class SplineFunction | |
92 | 168 |
93 /** | 169 /** |
94 * A row, typically a position where measurements were taken. | 170 * A row, typically a position where measurements were taken. |
95 */ | 171 */ |
96 public static final class Row | 172 public static final class Row |
164 ows[i] = Double.NaN; | 240 ows[i] = Double.NaN; |
165 } | 241 } |
166 } | 242 } |
167 } | 243 } |
168 | 244 |
169 public static final class SplineFunction { | |
170 | |
171 public PolynomialSplineFunction spline; | |
172 public double [] splineQs; | |
173 public double [] splineWs; | |
174 | |
175 public SplineFunction( | |
176 PolynomialSplineFunction spline, | |
177 double [] splineQs, | |
178 double [] splineWs | |
179 ) { | |
180 this.spline = spline; | |
181 this.splineQs = splineQs; | |
182 this.splineWs = splineWs; | |
183 } | |
184 | |
185 public double [][] sample( | |
186 int numSamples, | |
187 double km, | |
188 Calculation errors | |
189 ) { | |
190 double q1 = splineQs[0], qN = splineQs[splineQs.length-1]; | |
191 double minQ = Math.min(q1, qN); | |
192 double maxQ = Math.max(q1, qN); | |
193 | |
194 double [] outWs = new double[numSamples]; | |
195 double [] outQs = new double[numSamples]; | |
196 | |
197 Arrays.fill(outWs, Double.NaN); | |
198 Arrays.fill(outQs, Double.NaN); | |
199 | |
200 double stepWidth = (maxQ - minQ)/numSamples; | |
201 | |
202 try { | |
203 double q = minQ; | |
204 for (int i = 0; i < outWs.length; ++i, q += stepWidth) { | |
205 outWs[i] = spline.value(outQs[i] = q); | |
206 } | |
207 } | |
208 catch (ArgumentOutsideDomainException aode) { | |
209 if (errors != null) { | |
210 // TODO: I18N | |
211 errors.addProblem(km, "spline interpolation failed"); | |
212 } | |
213 log.error("spline interpolation failed.", aode); | |
214 } | |
215 | |
216 return new double [][] { outWs, outQs }; | |
217 } | |
218 } // class SplineFunction | |
219 | 245 |
220 public SplineFunction createSpline( | 246 public SplineFunction createSpline( |
221 WstValueTable table, | 247 WstValueTable table, |
222 Calculation errors | 248 Calculation errors |
223 ) { | 249 ) { |
819 return new double [][] {qs, ws}; | 845 return new double [][] {qs, ws}; |
820 } | 846 } |
821 | 847 |
822 public double [] findQsForW(double km, double w) { | 848 public double [] findQsForW(double km, double w) { |
823 | 849 |
824 | |
825 int rowIndex = Collections.binarySearch(rows, new Row(km)); | 850 int rowIndex = Collections.binarySearch(rows, new Row(km)); |
826 | 851 |
827 if (rowIndex >= 0) { | 852 if (rowIndex >= 0) { |
828 return rows.get(rowIndex).findQsForW(w, this); | 853 return rows.get(rowIndex).findQsForW(w, this); |
829 } | 854 } |
830 | 855 |
831 rowIndex = -rowIndex - 1; | 856 rowIndex = -rowIndex - 1; |
832 | 857 |
833 if (rowIndex < 1 || rowIndex >= rows.size()) { | 858 if (rowIndex < 1 || rowIndex >= rows.size()) { |
834 // do not extrapolate | 859 // Do not extrapolate. |
835 return new double[0]; | 860 return new double[0]; |
836 } | 861 } |
837 | 862 |
838 | 863 // Needs bilinear interpolation. |
839 // Needs bilinear interpolation | |
840 Row r1 = rows.get(rowIndex-1); | 864 Row r1 = rows.get(rowIndex-1); |
841 Row r2 = rows.get(rowIndex); | 865 Row r2 = rows.get(rowIndex); |
842 | 866 |
843 return r1.findQsForW(r2, w, km, this); | 867 return r1.findQsForW(r2, w, km, this); |
868 } | |
869 | |
870 protected SplineFunction createSpline(double km, Calculation errors) { | |
871 | |
872 int rowIndex = Collections.binarySearch(rows, new Row(km)); | |
873 | |
874 if (rowIndex >= 0) { | |
875 SplineFunction sf = rows.get(rowIndex).createSpline(this, errors); | |
876 if (sf == null && errors != null) { | |
877 // TODO: I18N | |
878 errors.addProblem(km, "cannot create w/q relation"); | |
879 } | |
880 return sf; | |
881 } | |
882 | |
883 rowIndex = -rowIndex - 1; | |
884 | |
885 if (rowIndex < 1 || rowIndex >= rows.size()) { | |
886 // Do not extrapolate. | |
887 if (errors != null) { | |
888 // TODO: I18N | |
889 errors.addProblem(km, "km not found"); | |
890 } | |
891 return null; | |
892 } | |
893 | |
894 // Needs bilinear interpolation. | |
895 Row r1 = rows.get(rowIndex-1); | |
896 Row r2 = rows.get(rowIndex); | |
897 | |
898 SplineFunction sf = r1.createSpline(r2, km, this, errors); | |
899 if (sf == null && errors != null) { | |
900 // TODO: I18N | |
901 errors.addProblem(km, "cannot create w/q relation"); | |
902 } | |
903 | |
904 return sf; | |
905 } | |
906 | |
907 /** 'Bezugslinienverfahren' */ | |
908 public double [][] relateWs( | |
909 double km1, | |
910 double km2, | |
911 int numSamples, | |
912 Calculation errors | |
913 ) { | |
914 SplineFunction sf1 = createSpline(km1, errors); | |
915 if (sf1 == null) { | |
916 return new double[2][0]; | |
917 } | |
918 | |
919 SplineFunction sf2 = createSpline(km1, errors); | |
920 if (sf2 == null) { | |
921 return new double[2][0]; | |
922 } | |
923 | |
924 PolynomialSplineFunction iQ1 = sf1.createIndexQRelation(); | |
925 if (iQ1 == null) { | |
926 if (errors != null) { | |
927 // TODO: I18N | |
928 errors.addProblem(km1, "cannot create index/q relation"); | |
929 } | |
930 return new double[2][0]; | |
931 } | |
932 | |
933 PolynomialSplineFunction iQ2 = sf1.createIndexQRelation(); | |
934 if (iQ2 == null) { | |
935 if (errors != null) { | |
936 // TODO: I18N | |
937 errors.addProblem(km2, "cannot create index/q relation"); | |
938 } | |
939 return new double[2][0]; | |
940 } | |
941 | |
942 int N = sf1.splineQs.length; | |
943 double stepWidth = N/(double)numSamples; | |
944 | |
945 PolynomialSplineFunction qW1 = sf1.spline; | |
946 PolynomialSplineFunction qW2 = sf2.spline; | |
947 | |
948 double [] ws1 = new double[numSamples]; | |
949 double [] ws2 = new double[numSamples]; | |
950 | |
951 Arrays.fill(ws1, Double.NaN); | |
952 Arrays.fill(ws2, Double.NaN); | |
953 | |
954 boolean hadErrors = false; | |
955 | |
956 double p = 0d; | |
957 for (int i = 0; i < numSamples; ++i, p += stepWidth) { | |
958 try { | |
959 double q1 = iQ1.value(p); | |
960 double w1 = qW1.value(q1); | |
961 double q2 = iQ2.value(p); | |
962 double w2 = qW2.value(q2); | |
963 ws1[i] = w1; | |
964 ws2[i] = w2; | |
965 } | |
966 catch (ArgumentOutsideDomainException aode) { | |
967 if (!hadErrors) { | |
968 // XXX: I'm not sure if this really can happen | |
969 // and if we should report this more than once. | |
970 hadErrors = true; | |
971 if (errors != null) { | |
972 // TODO: I18N | |
973 log.debug("W~W failed", aode); | |
974 errors.addProblem("W~W failed"); | |
975 } | |
976 } | |
977 } | |
978 } | |
979 | |
980 return new double [][] { ws1, ws2 }; | |
844 } | 981 } |
845 | 982 |
846 public QPosition getQPosition(double km, double q) { | 983 public QPosition getQPosition(double km, double q) { |
847 return getQPosition(km, q, new QPosition()); | 984 return getQPosition(km, q, new QPosition()); |
848 } | 985 } |