mschaefer@9432: /** Copyright (C) 2017 by Bundesanstalt für Gewässerkunde
gernotbelger@9313:  * Software engineering by
gernotbelger@9313:  *  Björnsen Beratende Ingenieure GmbH
gernotbelger@9313:  *  Dr. Schumacher Ingenieurbüro für Wasser und Umwelt
gernotbelger@9313:  *
gernotbelger@9313:  * This file is Free Software under the GNU AGPL (>=v3)
gernotbelger@9313:  * and comes with ABSOLUTELY NO WARRANTY! Check out the
gernotbelger@9313:  * documentation coming with Dive4Elements River for details.
gernotbelger@9313:  */
mschaefer@9432: 
gernotbelger@9313: package org.dive4elements.river.artifacts.bundu.bezugswst;
gernotbelger@9313: 
gernotbelger@9313: import java.util.ArrayList;
mschaefer@9450: import java.util.Calendar;
gernotbelger@9313: import java.util.List;
gernotbelger@9313: 
gernotbelger@9313: import org.dive4elements.artifacts.CallContext;
mschaefer@9432: import org.dive4elements.river.artifacts.access.FixRealizingAccess;
gernotbelger@9313: import org.dive4elements.river.artifacts.bundu.BUNDUArtifact;
gernotbelger@9318: import org.dive4elements.river.artifacts.bundu.BunduResultType;
mschaefer@9446: import org.dive4elements.river.artifacts.common.AbstractResultType;
gernotbelger@9313: import org.dive4elements.river.artifacts.common.GeneralResultType;
gernotbelger@9313: import org.dive4elements.river.artifacts.common.ResultRow;
mschaefer@9432: import org.dive4elements.river.artifacts.model.Calculation;
mschaefer@9432: import org.dive4elements.river.artifacts.model.Calculation.Problem;
gernotbelger@9313: import org.dive4elements.river.artifacts.model.CalculationResult;
mschaefer@9432: import org.dive4elements.river.artifacts.model.WQKms;
mschaefer@9432: import org.dive4elements.river.artifacts.model.fixings.FixRealizingCalculation;
mschaefer@9432: import org.dive4elements.river.artifacts.model.fixings.FixRealizingResult;
gernotbelger@9499: import org.dive4elements.river.artifacts.model.river.RiverInfoProvider;
gernotbelger@9313: import org.dive4elements.river.artifacts.resources.Resources;
mschaefer@9432: import org.dive4elements.river.artifacts.sinfo.tkhstate.BedHeightsFinder;
mschaefer@9432: import org.dive4elements.river.artifacts.sinfo.tkhstate.WinfoArtifactWrapper;
gernotbelger@9313: import org.dive4elements.river.artifacts.sinfo.util.CalculationUtils;
gernotbelger@9313: import org.dive4elements.river.artifacts.sinfo.util.RiverInfo;
gernotbelger@9323: import org.dive4elements.river.artifacts.sinfo.util.WstInfo;
mschaefer@9528: import org.dive4elements.river.artifacts.states.WaterlevelData;
mschaefer@9432: import org.dive4elements.river.exports.WaterlevelDescriptionBuilder;
mschaefer@9444: import org.dive4elements.river.model.BedHeightValueType;
gernotbelger@9313: import org.dive4elements.river.model.River;
mschaefer@9451: import org.dive4elements.river.utils.DoubleUtil;
gernotbelger@9313: 
gernotbelger@9313: class BezugswstCalculation {
gernotbelger@9313: 
mschaefer@9432:     // private static Logger log = Logger.getLogger(BezugswstCalculation.class);
gernotbelger@9313: 
mschaefer@9446:     /**
mschaefer@9446:      * Additional depth (m) to compute the excavation volume
mschaefer@9446:      */
mschaefer@9446:     private final static double EXCAVATION_DEPTH = 0.2; // REMARK Sollte von außen einstellbar sein
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Excavation costs (euro) per cubic meter
mschaefer@9446:      */
mschaefer@9446:     private final static double EXPENSE_PER_CBM = 12.0; // REMARK Sollte von außen einstellbar sein
mschaefer@9446: 
gernotbelger@9313:     private final CallContext context;
gernotbelger@9313: 
mschaefer@9446:     private final List<ResultRow> rows;
mschaefer@9446: 
mschaefer@9450:     private Double missKmFrom;
mschaefer@9450: 
mschaefer@9450:     private Double missKmTo;
mschaefer@9450: 
gernotbelger@9313:     public BezugswstCalculation(final CallContext context) {
gernotbelger@9313:         this.context = context;
mschaefer@9446:         this.rows = new ArrayList<>();
gernotbelger@9313:     }
gernotbelger@9313: 
mschaefer@9432:     /**
mschaefer@9432:      * Calculates the result rows of a bundu bzws workflow
mschaefer@9432:      */
gernotbelger@9313:     public CalculationResult calculate(final BUNDUArtifact bunduartifact) {
gernotbelger@9313: 
mschaefer@9432:         // Get input data
gernotbelger@9313:         final String user = CalculationUtils.findArtifactUser(this.context, bunduartifact);
gernotbelger@9313:         final BunduAccess access = new BunduAccess(bunduartifact);
gernotbelger@9313:         final River river = access.getRiver();
gernotbelger@9313:         final RiverInfo riverInfo = new RiverInfo(river);
mschaefer@9432:         final String calcModeLabel = Resources.getMsg(this.context.getMeta(), "bundu_bezugswst");
mschaefer@9432:         final boolean preprocessing = access.getPreprocessing();
gernotbelger@9323:         final int startYear = access.getStartYear();
gernotbelger@9323:         final int endYear = access.getBezugsJahr();
gernotbelger@9323:         final Integer ud = access.getUd();
mschaefer@9450:         this.missKmFrom = access.getMissingVolFrom();
mschaefer@9450:         this.missKmTo = access.getMissingVolTo();
gernotbelger@9330: 
gernotbelger@9495:         final BezugswstCalculationResults results = new BezugswstCalculationResults(calcModeLabel, user, riverInfo, access.getRange(),
gernotbelger@9495:                 access.isCalculateMissingVolume());
gernotbelger@9330: 
mschaefer@9432:         final Calculation problems = new Calculation();
gernotbelger@9331: 
mschaefer@9432:         // Calculate the wspl for the selected river range as in fixa awspl
gernotbelger@9465:         bunduartifact.addStringData("wq_isq", "true");
mschaefer@9432:         final WinfoArtifactWrapper winfo = new WinfoArtifactWrapper(bunduartifact);
mschaefer@9528:         final RiverInfoProvider riverInfoProvider = RiverInfoProvider.forRange(this.context, river, access.getRange());
mschaefer@9438:         final FixRealizingResult fixResult = calculateWspl(bunduartifact, problems);
mschaefer@9438:         if (fixResult == null)
mschaefer@9432:             return new CalculationResult(results, problems);
mschaefer@9438: 
mschaefer@9438:         final WQKms wqkms = fixResult.getWQKms()[0];
mschaefer@9528:         // We have no wst year as the wst is created by a calculation; we do not need it though
mschaefer@9528:         final int wspYear = -1;
mschaefer@9528:         // Remark: showAllGauges true for Fixierungsanalyse, false for WInfo, so true here as well
mschaefer@9528:         final boolean showAllGauges = true;
mschaefer@9528:         final WaterlevelData waterlevel = new WaterlevelData(wqkms, wspYear, showAllGauges, true);
mschaefer@9528:         final RiverInfoProvider riverInfoProvider2 = riverInfoProvider.forWaterlevel(waterlevel);
mschaefer@9528:         final WstInfo wstInfo = new WstInfo(wqkms.getName(), 0, riverInfoProvider2.getReferenceGauge(), true);
gernotbelger@9331: 
gernotbelger@9573:         // Fetch the bed levels of the selected sounding
mschaefer@9480:         final Integer bedHeightId = access.getBedHeightID();
mschaefer@9480:         final BedHeightsFinder bedHeightsFinder = (bedHeightId != null) ? BedHeightsFinder.forId(problems, bedHeightId, access.getRange())
mschaefer@9480:                 : BedHeightsFinder.NullFinder();
gernotbelger@9332: 
mschaefer@9432:         // Fetch the river channel data
mschaefer@9432:         final ChannelFinder channelFinder = ChannelFinder.loadValues(problems, river, access.getBezugsJahr());
mschaefer@9432:         if (channelFinder == null)
mschaefer@9432:             return new CalculationResult(results, problems);
gernotbelger@9332: 
mschaefer@9446:         // Compute the result rows
mschaefer@9432:         for (int i = 0; i <= wqkms.size() - 1; i++) {
mschaefer@9528:             this.rows.add(createRow(wqkms.getKm(i), wqkms.getW(i), wqkms.getQ(i), bedHeightsFinder, channelFinder, riverInfoProvider2, wstInfo));
mschaefer@9446:         }
mschaefer@9446: 
mschaefer@9446:         // Compute the missing volumes
gernotbelger@9495:         if (access.isCalculateMissingVolume()) {
mschaefer@9480:             if ((bedHeightsFinder == null) || bedHeightsFinder.isNull())
gernotbelger@9465:                 return new CalculationResult(results, problems);
mschaefer@9450:             computeMissingVolumes(problems);
mschaefer@9451:             final BedQualityCalculator bqCalculator = computeDensities(problems, bunduartifact, access, river);
mschaefer@9451:             computeMissingMasses(problems, bqCalculator);
gernotbelger@9330:         }
gernotbelger@9330: 
mschaefer@9432:         // Add the result to the results collection
mschaefer@9432:         final WaterlevelDescriptionBuilder descBuilder = new WaterlevelDescriptionBuilder(winfo, this.context);
mschaefer@9432:         final String qtext = descBuilder.getMetadataQ();
mschaefer@9446:         final BezugswstMainCalculationResult result = new BezugswstMainCalculationResult("bundu-bzws", this.rows, bedHeightsFinder.getInfo(), wstInfo,
mschaefer@9450:                 access.getFunction(), preprocessing, startYear, endYear, ud, qtext, wqkms, this.missKmFrom, this.missKmTo);
mschaefer@9432:         results.addResult(result, problems);
mschaefer@9432: 
mschaefer@9446:         // Create the missing volume results
mschaefer@9446:         if (access.getMissingVolFrom() != null) {
mschaefer@9446:             final String title1 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result1");
mschaefer@9446:             final BezugswstMissVolCalculationResult1 r1 = new BezugswstMissVolCalculationResult1(title1, this.rows);
mschaefer@9446:             results.addResult(r1, null);
mschaefer@9446: 
mschaefer@9446:             final String title2 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result2");
mschaefer@9450:             final List<ResultRow> rows2 = copyMissRows();
mschaefer@9450:             final BezugswstMissVolCalculationResult2 r2 = new BezugswstMissVolCalculationResult2(title2, rows2);
mschaefer@9446:             results.addResult(r2, null);
mschaefer@9446: 
mschaefer@9446:             final String title3 = Resources.getMsg(this.context.getMeta(), "bundu.export.csv.title.bezugswst.result3");
mschaefer@9446:             final List<ResultRow> totalRows = new ArrayList<>();
mschaefer@9450:             totalRows.add(createTotalsRow(problems));
mschaefer@9446:             final BezugswstMissVolCalculationResult3 r3 = new BezugswstMissVolCalculationResult3(title3, totalRows);
mschaefer@9446:             results.addResult(r3, null);
mschaefer@9446:         }
mschaefer@9432: 
mschaefer@9432:         return new CalculationResult(results, problems);
gernotbelger@9313:     }
gernotbelger@9313: 
gernotbelger@9313:     /**
mschaefer@9432:      * Calculates a w-q-longitudinal section for a river range and Q specified in an artifact
gernotbelger@9313:      */
mschaefer@9438:     private FixRealizingResult calculateWspl(final BUNDUArtifact bundu, final Calculation problems) {
gernotbelger@9313: 
mschaefer@9432:         final FixRealizingAccess access = new FixRealizingAccess(bundu);
mschaefer@9432:         final FixRealizingCalculation calc = new FixRealizingCalculation(access);
mschaefer@9432: 
gernotbelger@9479:         final CalculationResult res = calc.calculate(this.context.getMeta());
mschaefer@9432: 
mschaefer@9432:         final FixRealizingResult fixRes = (FixRealizingResult) res.getData();
mschaefer@9432: 
mschaefer@9432:         final List<Problem> problems2 = res.getReport().getProblems();
mschaefer@9432:         if (problems2 != null) {
mschaefer@9432:             for (final Problem problem : problems2) {
mschaefer@9432:                 problems.addProblem(problem);
mschaefer@9432:             }
mschaefer@9432:         }
mschaefer@9438:         return fixRes;
mschaefer@9432:     }
mschaefer@9432: 
mschaefer@9432:     /**
mschaefer@9432:      * Create a result row for a station
mschaefer@9432:      */
mschaefer@9432:     private ResultRow createRow(final double station, final double w, final double q, final BedHeightsFinder bedHeightsFinder,
mschaefer@9432:             final ChannelFinder channelFinder, final RiverInfoProvider riverInfoProv, final WstInfo wstInfo) {
mschaefer@9432: 
mschaefer@9446:         // Set W and Q
mschaefer@9432:         final ResultRow row = ResultRow.create();
mschaefer@9432:         row.putValue(GeneralResultType.station, station);
mschaefer@9432:         row.putValue(BunduResultType.bezugswst, w);
mschaefer@9432:         row.putValue(GeneralResultType.dischargeQwithUnit, q);
mschaefer@9432:         row.putValue(GeneralResultType.waterlevelLabel, wstInfo.getLabel());
mschaefer@9432:         row.putValue(GeneralResultType.gaugeLabel, riverInfoProv.findGauge(station));
mschaefer@9432:         row.putValue(GeneralResultType.location, riverInfoProv.getLocation(station));
mschaefer@9480:         if (bedHeightsFinder.getInfo() != null)
mschaefer@9480:             row.putValue(BunduResultType.sounding, bedHeightsFinder.getInfo().getDescription());
mschaefer@9480:         else
mschaefer@9480:             row.putValue(BunduResultType.sounding, "");
mschaefer@9446: 
mschaefer@9446:         // Set bed and channel bottom height
mschaefer@9432:         final double msh = bedHeightsFinder.getMeanBedHeight(station);
mschaefer@9535:         row.putValue(BunduResultType.heightMeanBed, msh);
mschaefer@9432:         if (!Double.isNaN(w) && !Double.isNaN(msh))
mschaefer@9535:             row.putValue(BunduResultType.flowdepthMeanBed, w - msh);
mschaefer@9432:         else
mschaefer@9535:             row.putValue(BunduResultType.flowdepthMeanBed, Double.NaN);
gernotbelger@9465: 
mschaefer@9432:         final double channelDepth = channelFinder.getDepth(station);
mschaefer@9432:         row.putValue(BunduResultType.channelDepth, channelDepth);
mschaefer@9446:         double channelHeight;
mschaefer@9432:         if (!Double.isNaN(w) && !Double.isNaN(channelDepth))
mschaefer@9505:             channelHeight = w - channelDepth;
mschaefer@9432:         else
mschaefer@9446:             channelHeight = Double.NaN;
mschaefer@9446:         row.putValue(BunduResultType.channelLowerEdge, channelHeight);
mschaefer@9446:         final double channelWidth = channelFinder.getWidth(station);
mschaefer@9446:         row.putValue(BunduResultType.channelWidth, channelWidth);
mschaefer@9450:         if (!Double.isNaN(channelHeight)) {
mschaefer@9450:             if (msh > channelHeight + 0.001)
mschaefer@9450:                 row.putValue(BunduResultType.missDepthMeanBed, msh - channelHeight);
mschaefer@9450:             else
mschaefer@9450:                 row.putValue(BunduResultType.missDepthMeanBed, 0.0);
mschaefer@9450:         }
mschaefer@9446: 
mschaefer@9446:         // Set field heights and missing heights
mschaefer@9444:         final List<Double> fieldHeights = new ArrayList<>();
mschaefer@9444:         final List<Double> fieldDepths = new ArrayList<>();
mschaefer@9446:         final List<Double> fieldMissDepths = new ArrayList<>();
mschaefer@9446:         final List<Double> fieldMissWidths = new ArrayList<>();
mschaefer@9450:         final List<Double> fieldNulls = new ArrayList<>();
mschaefer@9446:         int missFieldCnt = 0;
mschaefer@9444:         for (int i = BedHeightValueType.FIELD_FIRST_INDEX; i <= BedHeightValueType.FIELD_LAST_INDEX; i++) {
mschaefer@9444:             final double h = bedHeightsFinder.getFieldHeight(station, i);
mschaefer@9444:             fieldHeights.add(Double.valueOf(h));
mschaefer@9446:             fieldDepths.add(Double.valueOf(w - h));
mschaefer@9446:             if (h > channelHeight + 0.001) {
mschaefer@9446:                 missFieldCnt++;
mschaefer@9446:                 fieldMissDepths.add(Double.valueOf(h - channelHeight));
mschaefer@9446:                 fieldMissWidths.add(Double.valueOf(channelWidth / BedHeightValueType.FIELD_LAST_INDEX));
gernotbelger@9448:             } else {
mschaefer@9446:                 fieldMissDepths.add(Double.valueOf(0.0));
mschaefer@9446:                 fieldMissWidths.add(Double.valueOf(0.0));
mschaefer@9446:             }
mschaefer@9450:             fieldNulls.add(Double.NaN);
mschaefer@9444:         }
mschaefer@9450:         if (isKmInMissingVolumeRange(station)) {
mschaefer@9450:             row.putValue(BunduResultType.missDepthFields, fieldMissDepths);
mschaefer@9450:             row.putValue(BunduResultType.missWidthFields, fieldMissWidths);
mschaefer@9450:             row.putValue(BunduResultType.hasMissingDepth, (missFieldCnt >= 1));
gernotbelger@9465:         } else {
mschaefer@9450:             row.putValue(BunduResultType.missDepthFields, fieldNulls);
mschaefer@9450:             row.putValue(BunduResultType.missWidthFields, fieldNulls);
mschaefer@9450:             row.putValue(BunduResultType.hasMissingDepth, null);
mschaefer@9450:         }
mschaefer@9450:         row.putValue(BunduResultType.missVolumeFields, fieldNulls);
mschaefer@9450:         row.putValue(BunduResultType.missMassFields, fieldNulls);
mschaefer@9444:         row.putValue(BunduResultType.bedHeightFields, fieldHeights);
mschaefer@9444:         row.putValue(BunduResultType.depthFields, fieldDepths);
mschaefer@9446: 
mschaefer@9446:         // Preset the missing volume fields with NaN
mschaefer@9446:         row.putValue(BunduResultType.excavationCosts, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.excavationVolume, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missVolumeMeanBed, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missMassMeanBed, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missVolumeTotal, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missMassTotal, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.density, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missStationRangeFrom, Double.NaN);
mschaefer@9446:         row.putValue(BunduResultType.missStationRangeTo, Double.NaN);
mschaefer@9446: 
mschaefer@9432:         return row;
mschaefer@9432:     }
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Computes the missing volumes in a km range
mschaefer@9446:      */
mschaefer@9450:     private void computeMissingVolumes(final Calculation problems) {
mschaefer@9446:         // Search start km
mschaefer@9446:         int first = -1;
mschaefer@9446:         for (int j = 0; j <= this.rows.size() - 1; j++) {
mschaefer@9450:             if (isKmInMissingVolumeRange(this.rows.get(j).getDoubleValue(GeneralResultType.station))) {
mschaefer@9446:                 first = j;
mschaefer@9446:                 break;
mschaefer@9446:             }
mschaefer@9446:         }
mschaefer@9446:         if (first < 0)
mschaefer@9446:             return;
mschaefer@9535:         // Calculate all kms in missing volume calc range
mschaefer@9535:         double km;
mschaefer@9446:         int last = this.rows.size() - 1;
mschaefer@9450:         for (int i = first; i <= this.rows.size() - 1; i++) {
mschaefer@9535:             km = this.rows.get(i).getDoubleValue(GeneralResultType.station);
mschaefer@9535:             if (!isKmInMissingVolumeRange(km))
mschaefer@9446:                 break;
mschaefer@9535:             if (km > this.missKmTo.doubleValue() - 0.0001)
mschaefer@9446:                 last = i;
mschaefer@9446:             final List<Double> areas = new ArrayList<>();
mschaefer@9446:             final List<Double> volumes = new ArrayList<>();
mschaefer@9446:             double vTotal = 0.0;
mschaefer@9446:             double vExcav = 0.0;
mschaefer@9446:             for (int j = BedHeightValueType.FIELD_FIRST_INDEX; j <= BedHeightValueType.FIELD_LAST_INDEX; j++) {
mschaefer@9446:                 if (getFieldValue(i, BunduResultType.missDepthFields, j) > 0.0001) {
mschaefer@9446:                     computeMissingVolume(volumes, areas, i, first, last, j);
mschaefer@9446:                     vTotal += volumes.get(j - 1);
mschaefer@9446:                     vExcav += volumes.get(j - 1) + areas.get(j - 1) * EXCAVATION_DEPTH;
gernotbelger@9448:                 } else {
mschaefer@9446:                     volumes.add(Double.valueOf(0.0));
mschaefer@9446:                     areas.add(Double.valueOf(0.0));
mschaefer@9446:                 }
mschaefer@9446:             }
mschaefer@9450:             final double[] meanBedVolumeArea = computeMeanBedMissingAreaAndVolume(i, first, last);
mschaefer@9450:             this.rows.get(i).putValue(BunduResultType.missVolumeMeanBed, meanBedVolumeArea[0]);
mschaefer@9450:             this.rows.get(i).putValue(BunduResultType.missAreaMeanBed, meanBedVolumeArea[1]);
mschaefer@9446:             this.rows.get(i).putValue(BunduResultType.missVolumeFields, volumes);
mschaefer@9446:             this.rows.get(i).putValue(BunduResultType.missAreaFields, areas);
mschaefer@9446:             this.rows.get(i).putValue(BunduResultType.missVolumeTotal, vTotal);
mschaefer@9446:             this.rows.get(i).putValue(BunduResultType.excavationVolume, vExcav);
mschaefer@9450:             this.rows.get(i).putValue(BunduResultType.excavationCosts, vExcav * EXPENSE_PER_CBM);
mschaefer@9446:         }
mschaefer@9446:     }
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Computes the missing volume of a field of a km row
mschaefer@9446:      */
mschaefer@9446:     private void computeMissingVolume(final List<Double> volumes, final List<Double> areas, final int current, final int first, final int last,
mschaefer@9446:             final int field) {
mschaefer@9450: 
mschaefer@9446:         final double areaCurr = missingArea(current, first, last, field);
mschaefer@9446:         final double areaPrev = missingArea(current - 1, first, last, field);
mschaefer@9446:         final double areaNext = missingArea(current + 1, first, last, field);
mschaefer@9446:         final double kmCurr = missingKm(current);
mschaefer@9446:         final double kmPrev = missingKm(current - 1);
mschaefer@9446:         final double kmNext = missingKm(current + 1);
mschaefer@9450:         final double area1 = Double.isNaN(kmPrev) ? 0.0 : 0.5 * (areaCurr + areaPrev);
mschaefer@9450:         final double area2 = Double.isNaN(kmNext) ? 0.0 : 0.5 * (areaCurr + areaNext);
mschaefer@9450:         final double volume = Double.valueOf((Math.abs(kmCurr - kmPrev) * 500 * area1) + (Math.abs(kmNext - kmCurr) * 500 * area2));
mschaefer@9450:         volumes.add(volume);
mschaefer@9450:         if (!Double.isNaN(volume))
mschaefer@9446:             areas.add(Double.valueOf(area1 + area2));
mschaefer@9450:         else
mschaefer@9450:             areas.add(Double.NaN);
mschaefer@9446:     }
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Gets the missing area of a field and a row if in range, otherwise 0.0
mschaefer@9446:      */
mschaefer@9446:     private double missingArea(final int rowIndex, final int first, final int last, final int fieldIndex) {
mschaefer@9446:         if ((first <= rowIndex) && (rowIndex <= last))
gernotbelger@9448:             return getFieldValue(rowIndex, BunduResultType.missDepthFields, fieldIndex) * getFieldValue(rowIndex, BunduResultType.missWidthFields, fieldIndex);
mschaefer@9446:         else
mschaefer@9446:             return 0.0;
mschaefer@9446:     }
mschaefer@9446: 
mschaefer@9446:     /**
gernotbelger@9573:      * Computes the missing area and volume of the mean bed level of a km row
mschaefer@9450:      */
mschaefer@9450:     private double[] computeMeanBedMissingAreaAndVolume(final int current, final int first, final int last) {
mschaefer@9450: 
mschaefer@9450:         final double areaCurr = meanBedMissingArea(current, first, last);
mschaefer@9535:         if (areaCurr < 0.0001)
mschaefer@9535:             return new double[] { 0.0, 0.0 };
mschaefer@9450:         final double areaPrev = meanBedMissingArea(current - 1, first, last);
mschaefer@9450:         final double areaNext = meanBedMissingArea(current + 1, first, last);
mschaefer@9450:         final double kmCurr = missingKm(current);
mschaefer@9450:         final double kmPrev = missingKm(current - 1);
mschaefer@9450:         final double kmNext = missingKm(current + 1);
mschaefer@9450:         final double area1 = Double.isNaN(kmPrev) ? 0.0 : 0.5 * (areaCurr + areaPrev);
mschaefer@9450:         final double area2 = Double.isNaN(kmNext) ? 0.0 : 0.5 * (areaCurr + areaNext);
mschaefer@9450:         final double volume = Double.valueOf((Math.abs(kmCurr - kmPrev) * 500 * area1) + (Math.abs(kmNext - kmCurr) * 500 * area2));
mschaefer@9450:         final double area = Double.isNaN(volume) ? Double.NaN : Double.valueOf(area1 + area2);
mschaefer@9450:         return new double[] { volume, area };
mschaefer@9450:     }
mschaefer@9450: 
mschaefer@9450:     /**
gernotbelger@9573:      * Gets the missing area of the mean bed level and a row if in range, otherwise 0.0
mschaefer@9450:      */
mschaefer@9450:     private double meanBedMissingArea(final int rowIndex, final int first, final int last) {
mschaefer@9450:         if ((first <= rowIndex) && (rowIndex <= last)) {
mschaefer@9450:             final double dh = this.rows.get(rowIndex).getDoubleValue(BunduResultType.channelDepth)
mschaefer@9535:                     - this.rows.get(rowIndex).getDoubleValue(BunduResultType.flowdepthMeanBed);
mschaefer@9450:             if (dh > 0.0)
mschaefer@9450:                 return dh * this.rows.get(rowIndex).getDoubleValue(BunduResultType.channelWidth);
mschaefer@9450:             return 0.0;
mschaefer@9450:         }
mschaefer@9450:         return 0.0;
mschaefer@9450:     }
mschaefer@9450: 
mschaefer@9450:     /**
mschaefer@9446:      * Gets the km of a row if within range, otherwise NaN
mschaefer@9446:      */
mschaefer@9446:     private double missingKm(final int rowIndex) {
mschaefer@9446:         if ((0 <= rowIndex) && (rowIndex <= this.rows.size() - 1))
mschaefer@9446:             return this.rows.get(rowIndex).getDoubleValue(GeneralResultType.station);
mschaefer@9450:         return Double.NaN;
mschaefer@9450:     }
mschaefer@9450: 
mschaefer@9450:     /**
mschaefer@9451:      * Create a density calculator and compute the densities of the missing volume km range and time period according to the
mschaefer@9451:      * reference year
mschaefer@9451:      */
mschaefer@9451:     private BedQualityCalculator computeDensities(final Calculation problems, final BUNDUArtifact bunduartifact, final BunduAccess access, final River river) {
mschaefer@9451:         final BedQualityCalculator bqCalculator = new BedQualityCalculator(this.context, bunduartifact);
mschaefer@9451:         // REMARK 10km tolerance at start and end to enable interpolation there
mschaefer@9451:         final double[] kms = DoubleUtil.explode(access.getMissingVolFrom().doubleValue() - 10.0, access.getMissingVolTo().doubleValue() + 10.0,
mschaefer@9451:                 access.getStep().doubleValue() / 1000);
mschaefer@9451:         final Calendar endDay = Calendar.getInstance();
mschaefer@9451:         endDay.set(access.getBezugsJahr().intValue(), 11, 31);
mschaefer@9451:         final Calendar startDay = Calendar.getInstance();
mschaefer@9451:         // TODO Spezialregelung für den Rhein (bis 1999, 2000 bis 2009, ab 2010)
mschaefer@9451:         startDay.set(endDay.get(Calendar.YEAR) - 20, 0, 1);
mschaefer@9451:         bqCalculator.execute(problems, river, kms, startDay.getTime(), endDay.getTime());
mschaefer@9451:         return bqCalculator;
mschaefer@9451:     }
mschaefer@9451: 
mschaefer@9451:     /**
mschaefer@9450:      * Computes the missing masses
mschaefer@9450:      */
mschaefer@9451:     private void computeMissingMasses(final Calculation problems, final BedQualityCalculator densityFinder) {
mschaefer@9450:         for (final ResultRow row : this.rows) {
mschaefer@9450:             @SuppressWarnings("unchecked")
mschaefer@9450:             final List<Double> volumes = (List<Double>) row.getValue(BunduResultType.missVolumeFields);
mschaefer@9450:             if ((volumes == null) || Double.isNaN(volumes.get(0)))
mschaefer@9450:                 continue;
mschaefer@9451:             final double density = getDensity(row.getDoubleValue(GeneralResultType.station), densityFinder);
mschaefer@9450:             final List<Double> masses = new ArrayList<>();
mschaefer@9450:             double kmTotal = 0.0;
mschaefer@9450:             for (int j = BedHeightValueType.FIELD_FIRST_INDEX; j <= BedHeightValueType.FIELD_LAST_INDEX; j++) {
mschaefer@9450:                 final double m = volumes.get(j - 1) * density;
mschaefer@9450:                 masses.add(m);
mschaefer@9450:                 if (!Double.isNaN(m))
mschaefer@9450:                     kmTotal += m;
mschaefer@9450:             }
mschaefer@9450:             row.putValue(BunduResultType.density, density);
mschaefer@9450:             row.putValue(BunduResultType.missMassFields, masses);
mschaefer@9450:             row.putValue(BunduResultType.missMassTotal, kmTotal);
mschaefer@9450:             row.putValue(BunduResultType.missMassMeanBed, row.getDoubleValue(BunduResultType.missVolumeMeanBed) * density);
mschaefer@9450:         }
mschaefer@9446:     }
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Gets a value of one of the field list types of a row
mschaefer@9446:      *
mschaefer@9446:      * @param rowIndex
mschaefer@9446:      * @param type
mschaefer@9446:      * @param fieldIndex
mschaefer@9446:      *            1-based field index
mschaefer@9446:      */
mschaefer@9446:     private double getFieldValue(final int rowIndex, final AbstractResultType type, final int fieldIndex) {
mschaefer@9446:         @SuppressWarnings("unchecked")
mschaefer@9446:         final List<Double> values = (List<Double>) this.rows.get(rowIndex).getValue(type);
mschaefer@9446:         return values.get(fieldIndex - 1);
mschaefer@9446:     }
mschaefer@9446: 
mschaefer@9446:     /**
mschaefer@9446:      * Computes the volume and mass total of all rows with missing volumes
mschaefer@9446:      */
mschaefer@9450:     private ResultRow createTotalsRow(final Calculation problems) {
mschaefer@9446:         // Search start km
mschaefer@9446:         double vTotal = 0.0;
mschaefer@9446:         double mTotal = 0.0;
mschaefer@9450:         double eTotal = 0.0;
mschaefer@9450:         double cTotal = 0.0;
mschaefer@9446:         for (final ResultRow row : this.rows) {
mschaefer@9450:             final double volume = row.getDoubleValue(BunduResultType.missVolumeTotal);
mschaefer@9450:             final double mass = row.getDoubleValue(BunduResultType.missMassTotal);
mschaefer@9450:             final double excavation = row.getDoubleValue(BunduResultType.excavationVolume);
mschaefer@9450:             final double costs = row.getDoubleValue(BunduResultType.excavationCosts);
mschaefer@9446:             if (!Double.isNaN(volume) && !Double.isNaN(mass)) {
mschaefer@9446:                 vTotal += volume;
mschaefer@9446:                 mTotal += mass;
mschaefer@9450:                 eTotal += excavation;
mschaefer@9450:                 cTotal += costs;
mschaefer@9446:             }
mschaefer@9446:         }
mschaefer@9446:         final ResultRow sumRow = ResultRow.create();
mschaefer@9450:         sumRow.putValue(BunduResultType.missStationRangeFrom, Double.valueOf(this.missKmFrom));
mschaefer@9450:         sumRow.putValue(BunduResultType.missStationRangeTo, Double.valueOf(this.missKmTo));
mschaefer@9446:         sumRow.putValue(BunduResultType.missVolumeTotal, vTotal);
mschaefer@9446:         sumRow.putValue(BunduResultType.missMassTotal, mTotal);
mschaefer@9480:         sumRow.putValue(BunduResultType.excavationVolumeTotal, eTotal);
mschaefer@9480:         sumRow.putValue(BunduResultType.excavationCostsTotal, cTotal);
mschaefer@9446:         return sumRow;
mschaefer@9446:     }
mschaefer@9450: 
mschaefer@9450:     /**
mschaefer@9450:      * Copies the rows of the missing volume calculation range into a new list
mschaefer@9450:      */
mschaefer@9450:     private List<ResultRow> copyMissRows() {
mschaefer@9450:         final List<ResultRow> missRows = new ArrayList<>();
mschaefer@9450:         for (final ResultRow row : this.rows) {
mschaefer@9450:             final double km = row.getDoubleValue(GeneralResultType.station);
mschaefer@9450:             if (isKmInMissingVolumeRange(km))
mschaefer@9450:                 missRows.add(row);
mschaefer@9450:         }
mschaefer@9450:         return missRows;
mschaefer@9450:     }
mschaefer@9450: 
mschaefer@9451:     /**
mschaefer@9451:      * Gets the density of a km from the densities calculation
mschaefer@9451:      */
mschaefer@9451:     private double getDensity(final double km, final BedQualityCalculator densityFinder) {
mschaefer@9451:         return densityFinder.getDensity(km);
mschaefer@9450:     }
mschaefer@9450: 
mschaefer@9450:     /**
mschaefer@9450:      * Checks whether a km lies in the missing volume calculation range
mschaefer@9450:      */
mschaefer@9450:     private boolean isKmInMissingVolumeRange(final double km) {
mschaefer@9450:         if ((this.missKmFrom == null) || (this.missKmTo == null))
mschaefer@9450:             return false;
mschaefer@9450:         return (this.missKmFrom.doubleValue() - 0.0001 < km) && (km < this.missKmTo.doubleValue() + 0.0001);
mschaefer@9450:     }
gernotbelger@9313: }