Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/Calculation4.java @ 655:913b52064449

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactored version of "Berechnung 4" flys-artifacts/trunk@2053 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Sun, 05 Jun 2011 18:24:46 +0000
parents
children baea7981477a
comparison
equal deleted inserted replaced
654:bbc966c81809 655:913b52064449
1 package de.intevation.flys.artifacts.model;
2
3 import java.util.List;
4 import java.util.Arrays;
5 import java.util.Comparator;
6 import java.util.Collections;
7
8 import de.intevation.flys.utils.DoubleUtil;
9
10 import de.intevation.flys.model.River;
11 import de.intevation.flys.model.Gauge;
12
13 import de.intevation.flys.artifacts.model.WstValueTable.QPosition;
14
15 import de.intevation.flys.artifacts.math.Function;
16 import de.intevation.flys.artifacts.math.Linear;
17 import de.intevation.flys.artifacts.math.Identity;
18 import de.intevation.flys.artifacts.math.BackJumpCorrector;
19
20 import org.apache.log4j.Logger;
21
22 public class Calculation4
23 {
24 private static Logger logger = Logger.getLogger(Calculation4.class);
25
26 public static final double MINIMAL_STEP_WIDTH = 1e-5;
27
28 public static final Comparator<Segment> REF_CMP =
29 new Comparator<Segment>() {
30 public int compare(Segment a, Segment b) {
31 double d = a.referencePoint - b.referencePoint;
32 if (d < 0d) return -1;
33 return d > 0d ? +1 : 0;
34 }
35 };
36
37 protected List<Segment> segments;
38
39 protected boolean isQ;
40
41 public Calculation4() {
42 }
43
44 public Calculation4(List<Segment> segments, River river, boolean isQ) {
45
46 this.segments = segments;
47 this.isQ = isQ;
48
49 int numResults = 1;
50
51 // assign reference points
52 for (Segment segment: segments) {
53 Gauge gauge = river.determineGauge(
54 segment.getFrom(), segment.getTo());
55
56 segment.setReferencePoint(gauge != null
57 ? gauge.getStation().doubleValue()
58 : 0.5*(segment.getFrom() + segment.getTo()));
59
60 double [] values = segment.values;
61
62 if (numResults == -1) {
63 numResults = values.length;
64 }
65 else if (numResults != values.length) {
66 throw new IllegalArgumentException("wrong length of values");
67 }
68
69 // convert to Q if needed
70 if (!isQ && gauge != null) {
71 double [][] table = new DischargeTables(
72 river.getName(), gauge.getName()).getFirstTable();
73
74 // need the original values for naming
75 segment.backup();
76
77 for (int i = 0; i < values.length; ++i) {
78 values[i] = DischargeTables.getQForW(table, values[i]);
79 }
80 }
81 } // for all segments
82
83 Collections.sort(segments, REF_CMP);
84 }
85
86 public WQKms [] calculate(
87 WstValueTable table,
88 double from, double to, double step
89 ) {
90 boolean debug = logger.isDebugEnabled();
91
92 if (debug) {
93 logger.debug(
94 "calculate from " + from + " to " + to + " step " + step);
95 logger.debug("# segments: " + segments.size());
96 }
97
98 if (segments.isEmpty()) {
99 logger.debug("no segments found");
100 return new WQKms[0];
101 }
102
103 int numResults = segments.get(0).values.length;
104
105 if (numResults < 1) {
106 logger.debug("no values given");
107 return new WQKms[0];
108 }
109
110
111 WQKms [] results = new WQKms[numResults];
112 if (debug) {
113 logger.debug("after last segment -> gleichwertig");
114 }
115 for (int i = 0; i < results.length; ++i) {
116 results[i] = new WQKms();
117 }
118
119 if (Math.abs(step) < MINIMAL_STEP_WIDTH) {
120 step = MINIMAL_STEP_WIDTH;
121 }
122
123 if (from > to) {
124 step = -step;
125 }
126
127 QPosition [] qPositions = new QPosition[numResults];
128
129 Function [] functions = new Function[numResults];
130
131 double [] out = new double[2];
132
133 Segment sentinel = new Segment(Double.MAX_VALUE);
134 Segment s1 = sentinel, s2 = sentinel;
135
136 for (double pos = from;
137 from < to ? pos <= to : pos >= to;
138 pos = DoubleUtil.round(pos + step)
139 ) {
140 if (pos < s1.referencePoint || pos > s2.referencePoint) {
141 if (debug) {
142 logger.debug("need to find new interval for " + pos);
143 }
144 // find new interval
145 if (pos <= segments.get(0).referencePoint) {
146 // before first segment -> "gleichwertig"
147 if (debug) {
148 logger.debug("before first segment -> gleichwertig");
149 }
150 Segment first = segments.get(0);
151 double [] values = first.values;
152 double refPos = first.referencePoint;
153 for (int i = 0; i < qPositions.length; ++i) {
154 qPositions[i] = table.getQPosition(
155 refPos, values[i]);
156 }
157 sentinel.setReferencePoint(-Double.MAX_VALUE);
158 s1 = sentinel;
159 s2 = segments.get(0);
160 Arrays.fill(functions, Identity.IDENTITY);
161 }
162 else if (pos >= segments.get(segments.size()-1).referencePoint) {
163 // after last segment -> "gleichwertig"
164 if (debug) {
165 logger.debug("after last segment -> gleichwertig");
166 }
167 Segment last = segments.get(segments.size()-1);
168 double [] values = last.values;
169 double refPos = last.referencePoint;
170 for (int i = 0; i < qPositions.length; ++i) {
171 qPositions[i] = table.getQPosition(
172 refPos, values[i]);
173 }
174 sentinel.setReferencePoint(Double.MAX_VALUE);
175 s1 = last;
176 s2 = sentinel;
177 Arrays.fill(functions, Identity.IDENTITY);
178 }
179 else { // "ungleichwertig"
180 // find matching interval
181 if (debug) {
182 logger.debug("in segments -> ungleichwertig");
183 }
184 s1 = s2 = null;
185 for (int i = 1, N = segments.size(); i < N; ++i) {
186 Segment si1 = segments.get(i-1);
187 Segment si = segments.get(i);
188 if (debug) {
189 logger.debug("check " + pos + " in " +
190 si1.referencePoint + " - " + si.referencePoint);
191 }
192 if (pos >= si1.referencePoint
193 && pos <= si. referencePoint) {
194 s1 = si1;
195 s2 = si;
196 break;
197 }
198 }
199
200 if (s1 == null) {
201 throw new IllegalStateException("no interval found");
202 }
203
204 Segment anchor, free;
205
206 if (from > to) { anchor = s1; free = s2; }
207 else { anchor = s2; free = s1; }
208
209 // build transforms based on "gleichwertiger" phase
210 for (int i = 0; i < qPositions.length; ++i) {
211 QPosition qi = table.getQPosition(
212 anchor.referencePoint,
213 anchor.values[i]);
214
215 if ((qPositions[i] = qi) == null) {
216 // TODO: error report
217 functions[i] = Identity.IDENTITY;
218 }
219 else {
220 double qA = table.getQ(qi, anchor.referencePoint);
221 double qF = table.getQ(qi, free .referencePoint);
222
223 functions[i] = Double.isNaN(qA) || Double.isNaN(qF)
224 ? Identity.IDENTITY
225 : new Linear(
226 qA, qF,
227 anchor.values[i], free.values[i]);
228 }
229 } // build transforms
230 } // "ungleichwertiges" interval
231 } // find matching interval
232
233 for (int i = 0; i < qPositions.length; ++i) {
234 QPosition qPosition = qPositions[i];
235
236 if (qPosition == null) {
237 continue;
238 }
239
240 if (table.interpolate(pos, out, qPosition, functions[i])) {
241 results[i].add(out[0], out[1], pos);
242 }
243 else {
244 // TODO: error report
245 }
246 }
247 }
248
249 // Backjump correction
250 for (int i = 0; i < results.length; ++i) {
251 BackJumpCorrector bjc = new BackJumpCorrector();
252
253 double [] ws = results[i].getWs();
254 double [] kms = results[i].getKms();
255
256 if (bjc.doCorrection(kms, ws)) {
257 results[i] = new WQCKms(results[i], bjc.getCorrected());
258 }
259 // TODO: error report
260 }
261
262 // name the curves
263 for (int i = 0; i < results.length; ++i) {
264 results[i].setName(createName(i));
265 }
266
267 return results;
268 }
269
270 protected String createName(int index) {
271 // TODO: i18n
272 StringBuilder sb = new StringBuilder(isQ ? "Q" : "W");
273 sb.append(" benutzerdefiniert (");
274 for (int i = 0, N = segments.size(); i < N; ++i) {
275 if (i > 0) {
276 sb.append("; ");
277 }
278 Segment segment = segments.get(i);
279 sb.append((segment.backup != null
280 ? segment.backup
281 : segment.values)[index]);
282 }
283 sb.append(')');
284 return sb.toString();
285 }
286 }
287 // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf8 :
http://dive4elements.wald.intevation.org