Mercurial > dive4elements > river

comparison flys-artifacts/src/main/java/de/intevation/flys/artifacts/model/WstValueTable.java @ 1937:f07d64d5cbe1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
'W auf freier Strecke' calculation. Fetch corresponding Qs for given Ws from the WST model flys-artifacts/trunk@3318 c6561f87-3c4e-4783-a992-168aeb5c3f6f
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 25 Nov 2011 15:50:31 +0000
parents 9bec7d2f8c35
children 1d991c91285b
comparison
equal deleted inserted replaced
1936:0ad05cb691fc 1937:f07d64d5cbe1
17 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 17 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
18 18
19 import org.apache.commons.math.ArgumentOutsideDomainException; 19 import org.apache.commons.math.ArgumentOutsideDomainException;
20 20
21 import org.apache.commons.math.exception.MathIllegalArgumentException; 21 import org.apache.commons.math.exception.MathIllegalArgumentException;
22
23 import gnu.trove.TDoubleArrayList;
22 24
23 /** 25 /**
24 * W, Q and km data from database 'wsts' spiced with interpolation algorithms. 26 * W, Q and km data from database 'wsts' spiced with interpolation algorithms.
25 */ 27 */
26 public class WstValueTable 28 public class WstValueTable
360 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]); 362 ow = Linear.weight(weight, other.ws[index-1], other.ws[index]);
361 } 363 }
362 364
363 return Linear.weight(kmWeight, tw, ow); 365 return Linear.weight(kmWeight, tw, ow);
364 } 366 }
367
368 public double [] findQsForW(double w, WstValueTable table) {
369
370 TDoubleArrayList qs = new TDoubleArrayList();
371
372 if (ws.length > 0 && Math.abs(ws[0]-w) < 0.000001) {
373 double q = table.getQIndex(0, km);
374 if (!Double.isNaN(q)) {
375 qs.add(q);
376 }
377 }
378
379 for (int i = 1; i < ws.length; ++i) {
380 double w2 = ws[i];
381 if (Double.isNaN(w2)) {
382 continue;
383 }
384 if (Math.abs(w2-w) < 0.000001) {
385 double q = table.getQIndex(i, km);
386 if (!Double.isNaN(q)) {
387 qs.add(q);
388 }
389 continue;
390 }
391 double w1 = ws[i-1];
392 if (Double.isNaN(w1)) {
393 continue;
394 }
395
396 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
397 continue;
398 }
399
400 double q1 = table.getQIndex(i-1, km);
401 double q2 = table.getQIndex(i, km);
402 if (Double.isNaN(q1) || Double.isNaN(q2)) {
403 continue;
404 }
405
406 double q = Linear.linear(w, w1, w2, q1, q2);
407 qs.add(q);
408 }
409
410 return qs.toNativeArray();
411 }
412
413 public double [] findQsForW(
414 Row other,
415 double w,
416 double km,
417 WstValueTable table
418 ) {
419 TDoubleArrayList qs = new TDoubleArrayList();
420
421 double factor = Linear.factor(km, this.km, other.km);
422
423 if (ws.length > 0) {
424 double wt = Linear.weight(factor, ws[0], other.ws[0]);
425 if (!Double.isNaN(wt)) {
426 double q = table.getQIndex(0, km);
427 if (!Double.isNaN(q)) {
428 qs.add(q);
429 }
430 }
431 }
432
433 for (int i = 1; i < ws.length; ++i) {
434 double w2 = Linear.weight(factor, ws[i], other.ws[i]);
435 if (Double.isNaN(w2)) {
436 continue;
437 }
438 if (Math.abs(w2-w) < 0.000001) {
439 double q = table.getQIndex(i, km);
440 if (!Double.isNaN(q)) {
441 qs.add(q);
442 }
443 continue;
444 }
445 double w1 = Linear.weight(factor, ws[i-1], other.ws[i-1]);
446 if (Double.isNaN(w1)) {
447 continue;
448 }
449
450 if (w < Math.min(w1, w2) || w > Math.max(w1, w2)) {
451 continue;
452 }
453
454 double q1 = table.getQIndex(i-1, km);
455 double q2 = table.getQIndex(i, km);
456 if (Double.isNaN(q1) || Double.isNaN(q2)) {
457 continue;
458 }
459
460 double q = Linear.linear(w, w1, w2, q1, q2);
461 qs.add(q);
462 }
463
464 return qs.toNativeArray();
465 }
365 } // class Row 466 } // class Row
366 467
367 /** Rows in table. */ 468 /** Rows in table. */
368 protected List<Row> rows; 469 protected List<Row> rows;
369 470
666 int rowIndex = Collections.binarySearch(rows, new Row(km)); 767 int rowIndex = Collections.binarySearch(rows, new Row(km));
667 768
668 if (rowIndex < 0) { 769 if (rowIndex < 0) {
669 rowIndex = -rowIndex -1; 770 rowIndex = -rowIndex -1;
670 } 771 }
671 if (rowIndex >= rows.size()) 772 if (rowIndex >= rows.size()) {
672 rowIndex = rows.size() -1; 773 rowIndex = rows.size() -1;
774 }
673 775
674 Row startRow = rows.get(rowIndex); 776 Row startRow = rows.get(rowIndex);
675 777
676 for (int col = 0; col < columns.length; col++) { 778 for (int col = 0; col < columns.length; col++) {
677 qs[col] = columns[col].getQRangeTree().findQ(km); 779 qs[col] = columns[col].getQRangeTree().findQ(km);
701 } 803 }
702 804
703 return new double [][] {qs, ws}; 805 return new double [][] {qs, ws};
704 } 806 }
705 807
808 public double [] findQsForW(double km, double w) {
809
810
811 int rowIndex = Collections.binarySearch(rows, new Row(km));
812
813 if (rowIndex >= 0) {
814 return rows.get(rowIndex).findQsForW(w, this);
815 }
816
817 rowIndex = -rowIndex - 1;
818
819 if (rowIndex < 1 || rowIndex >= rows.size()) {
820 // do not extrapolate
821 return new double[0];
822 }
823
824
825 // Needs bilinear interpolation
826 Row r1 = rows.get(rowIndex-1);
827 Row r2 = rows.get(rowIndex);
828
829 return r1.findQsForW(r2, w, km, this);
830 }
831
706 public QPosition getQPosition(double km, double q) { 832 public QPosition getQPosition(double km, double q) {
707 return getQPosition(km, q, new QPosition()); 833 return getQPosition(km, q, new QPosition());
708 } 834 }
709 835
710 public QPosition getQPosition(double km, double q, QPosition qPosition) { 836 public QPosition getQPosition(double km, double q, QPosition qPosition) {
http://dive4elements.wald.intevation.org