sascha@925: package de.intevation.flys.artifacts.math;
sascha@925: 
sascha@925: import de.intevation.flys.artifacts.model.WKms;
sascha@925: import de.intevation.flys.artifacts.model.WKmsImpl;
sascha@925: 
sascha@925: import java.util.Arrays;
sascha@925: 
sascha@925: public abstract class WKmsOperation
sascha@925: {
sascha@925:     public static final double EPSILON = 1e-6;
sascha@925: 
sascha@925:     public static final class KmW 
sascha@925:     implements                Comparable<KmW>
sascha@925:     {
sascha@925:         protected double km;
sascha@925:         protected double w;
sascha@925: 
sascha@925:         public KmW(double km, double w) {
sascha@925:             this.km = km;
sascha@925:             this.w  = w;
sascha@925:         }
sascha@925: 
sascha@925:         public int compareTo(KmW other) {
sascha@925:             return km < other.km
sascha@925:                 ? -1
sascha@925:                 : km > other.km ? +1 : 0;
sascha@925:         }
sascha@925: 
sascha@925:         public boolean kmEquals(KmW other) {
sascha@925:             return Math.abs(km - other.km) < EPSILON;
sascha@925:         }
sascha@925: 
sascha@925:         public double subtract(KmW other) {
sascha@925:             return w - other.w;
sascha@925:         }
sascha@925:     } // class KmW
sascha@925: 
sascha@925:     public static final WKmsOperation SUBTRACTION = new WKmsOperation() {
sascha@925: 
sascha@925:         @Override
sascha@925:         public WKms operate(WKms a, WKms b) {
sascha@925:             return subtract(a, b);
sascha@925:         }
sascha@925:     };
sascha@925: 
sascha@925:     protected WKmsOperation() {
sascha@925:     }
sascha@925: 
sascha@925:     public abstract WKms operate(WKms a, WKms b);
sascha@925: 
felix@1148:     /**
felix@1148:      * Subtract two series from each other, interpolate values
felix@1148:      * missing in one series in the other.
felix@1148:      */
sascha@925:     public static WKms subtract(WKms minuend, WKms subtrahend) {
sascha@925: 
sascha@925:         int M = minuend   .size();
sascha@925:         int S = subtrahend.size();
sascha@925: 
sascha@925:         // Don't subtract empty sets
sascha@925:         if (M < 1 || S < 1) {
sascha@925:             return new WKmsImpl();
sascha@925:         }
sascha@925: 
sascha@925:         KmW [] ms = new KmW[M];
sascha@925:         KmW [] ss = new KmW[S];
sascha@925: 
sascha@925:         for (int i = 0; i < M; ++i) {
sascha@925:             ms[i] = new KmW(minuend.getKm(i), minuend.getW(i));
sascha@925:         }
sascha@925: 
sascha@925:         for (int i = 0; i < S; ++i) {
sascha@925:             ss[i] = new KmW(subtrahend.getKm(i), subtrahend.getW(i));
sascha@925:         }
sascha@925: 
sascha@925:         Arrays.sort(ms);
sascha@925:         Arrays.sort(ss);
sascha@925: 
sascha@925:         // no overlap -> empty result set
sascha@925:         if (ms[0].km > ss[S-1].km || ss[0].km > ms[M-1].km) {
sascha@925:             return new WKmsImpl();
sascha@925:         }
sascha@925: 
sascha@925:         WKmsImpl result = new WKmsImpl();
sascha@925: 
sascha@925:         int mi = 0;
sascha@925:         int si = 0;
sascha@925: 
sascha@925:         OUT: while (mi < M && si < S) {
sascha@925:             KmW m = ms[mi];
sascha@925:             KmW s = ss[si];
sascha@925: 
sascha@925:             if (m.km + EPSILON < s.km) {
sascha@925:                 // minuend is before subtrahend
sascha@925: 
sascha@925:                 while (ms[mi].km + EPSILON < s.km) {
sascha@925:                     if (++mi >= M) {
sascha@925:                         break OUT;
sascha@925:                     }
sascha@925:                 }
sascha@925: 
sascha@925:                 if (ms[mi].km + EPSILON > s.km) {
sascha@925:                     double mw = Linear.linear(
sascha@925:                         s.km,
sascha@925:                         ms[mi-1].km, ms[mi].km,
sascha@925:                         ms[mi-1].w,  ms[mi].w);
sascha@925:                     result.add(s.km, mw - s.w);
sascha@925:                     ++si;
sascha@925:                 }
sascha@925:                 else { // s.km == ms[mi].km
sascha@925:                     result.add(s.km, ms[mi].subtract(s));
sascha@925:                     ++mi;
sascha@925:                     ++si;
sascha@925:                 }
sascha@925:             }
sascha@925:             else if (m.km > s.km + EPSILON) {
sascha@925:                 // subtrahend is before minuend
sascha@925: 
sascha@925:                 while (m.km > ss[si].km + EPSILON) {
sascha@925:                     if (++si >= S) {
sascha@925:                         break OUT;
sascha@925:                     }
sascha@925:                 }
sascha@925: 
sascha@925:                 if (ss[si].km + EPSILON > m.km) {
sascha@925:                     double sw = Linear.linear(
sascha@925:                         m.km,
sascha@925:                         ss[si-1].km, ss[si].km,
sascha@925:                         ss[si-1].w,  ss[si].w);
sascha@925:                     result.add(m.km, m.w - sw);
sascha@925:                 }
sascha@925:                 else { // ss[si].km == m.km
sascha@925:                     result.add(m.km, m.subtract(ss[si]));
sascha@925:                     ++mi;
sascha@925:                     ++si;
sascha@925:                 }
sascha@925:             }
sascha@925:             else { // m.km == s.km
sascha@925:                 result.add(s.km, m.subtract(s));
sascha@925:                 ++mi;
sascha@925:                 ++si;
sascha@925:             }
sascha@925:         }
sascha@925: 
sascha@925:         return result;
sascha@925:     }
sascha@925: }
sascha@925: // vim:set ts=4 sw=4 si et sta sts=4 fenc=utf-8 :