CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
interpolate.cpp
Go to the documentation of this file.
1
11
12/* ===============================================================================================================================
13
14 This file is part of CoastalME, the Coastal Modelling Environment.
15
16 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
19
20 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
21
22===============================================================================================================================*/
23#include <assert.h>
24
25#include <cfloat>
26
27#include <vector>
28using std::vector;
29
30#include "cme.h"
31#include "simulation.h"
32
33//===============================================================================================================================
36//===============================================================================================================================
37double CSimulation::dGetInterpolatedValue(vector<double> const* pVdXdata, vector<double> const* pVdYdata, double dX, bool bExtrapolate)
38{
39 int const size = static_cast<int>(pVdXdata->size());
40
41 int i = 0; // Find left end of interval for interpolation
42
43 if (dX >= pVdXdata->at(size - 2)) // Special case: beyond right end
44 {
45 i = size - 2;
46 }
47
48 else
49 {
50 while (dX > pVdXdata->at(i + 1))
51 i++;
52 }
53
54 double const dXL = pVdXdata->at(i);
55 double dYL = pVdYdata->at(i);
56 double const dXR = pVdXdata->at(i + 1);
57 double dYR = pVdYdata->at(i + 1); // Points on either side (unless beyond ends)
58
59 if (! bExtrapolate) // If beyond ends of array and not extrapolating
60 {
61 if (dX < dXL)
62 dYR = dYL;
63
64 if (dX > dXR)
65 dYL = dYR;
66 }
67
68 double const ddYdX = (dYR - dYL) / (dXR - dXL); // Gradient
69
70 return (dYL + ddYdX * (dX - dXL)); // Linear interpolation
71}
72
73//===============================================================================================================================
76//===============================================================================================================================
77double CSimulation::dGetInterpolatedValue(vector<int> const* pVnXdata, vector<double> const* pVdYdata, int nX, bool bExtrapolate)
78{
79 unsigned int const nSize = static_cast<unsigned int>(pVnXdata->size());
80
81 int i = 0; // Find left end of interval for interpolation
82
83 if (nX >= pVnXdata->at(nSize - 2)) // Special case: beyond right end
84 {
85 i = nSize - 2;
86 }
87
88 else
89 {
90 while (nX > pVnXdata->at(i + 1))
91 i++;
92 }
93
94 int const nXL = pVnXdata->at(i);
95 int const nXR = pVnXdata->at(i + 1);
96
97 double dYL = pVdYdata->at(i);
98 double dYR = pVdYdata->at(i + 1); // Points on either side (unless beyond ends)
99
100 if (! bExtrapolate) // If beyond ends of array and not extrapolating
101 {
102 if (nX < nXL)
103 dYR = dYL;
104
105 if (nX > nXR)
106 dYL = dYR;
107 }
108
109 double const ddYdX = (dYR - dYL) / static_cast<double>(nXR - nXL); // Gradient
110
111 return dYL + ddYdX * static_cast<double>(nX - nXL); // Linear interpolation
112}
113
114//===============================================================================================================================
116//===============================================================================================================================
117int CSimulation::nFindIndex(vector<double> const* pVdX, double const dValueIn)
118{
119 double dLastValue = DBL_MAX;
120 int nIndexFound = 0;
121
122 for (unsigned int i = 0; i < pVdX->size(); ++i)
123 {
124 double const dThisValue = tAbs(dValueIn - pVdX->at(i));
125
126 if (dThisValue <= dLastValue)
127 {
128 dLastValue = dThisValue;
129 nIndexFound = i;
130 }
131 }
132
133 return nIndexFound;
134}
135
136//===============================================================================================================================
138//===============================================================================================================================
139vector<double> CSimulation::VdInterpolateCShoreProfileOutput(vector<double> const* pVdX, vector<double> const* pVdY, vector<double> const* pVdXNew)
140{
141 int const nXSize = static_cast<int>(pVdX->size());
142 int const nXNewSize = static_cast<int>(pVdXNew->size());
143
144 // assert(nXSize > 0);
145 // assert(nXNewSize > 0);
146
147 double dX;
148 double dY;
149 vector<double> VdYNew(nXNewSize, 0.0);
150
151 for (int i = 0; i < nXNewSize; ++i)
152 {
153 int const idx = nFindIndex(pVdX, pVdXNew->at(i));
154
155 if (pVdX->at(idx) > pVdXNew->at(i))
156 {
157 if (idx > 0)
158 {
159 dX = pVdX->at(idx) - pVdX->at(idx - 1);
160 dY = pVdY->at(idx) - pVdY->at(idx - 1);
161 }
162
163 else
164 {
165 dX = pVdX->at(idx + 1) - pVdX->at(idx);
166 dY = pVdY->at(idx + 1) - pVdY->at(idx);
167 }
168 }
169
170 else
171 {
172 if (idx < nXSize - 1)
173 {
174 dX = pVdX->at(idx + 1) - pVdX->at(idx);
175 dY = pVdY->at(idx + 1) - pVdY->at(idx);
176 }
177
178 else
179 {
180 dX = pVdX->at(idx) - pVdX->at(idx - 1);
181 dY = pVdY->at(idx) - pVdY->at(idx - 1);
182 }
183 }
184
185 // Safety check: this crashes (divide by zero) if there are identical consecutive values in pVdX, and thus if dX becomes 0. To prevent this, if dX is near zero, set to a small non-zero number
186 if (bFPIsEqual(dX, 0.0, TOLERANCE))
187 dX = 1e-10;
188
189 double const dM = dY / dX;
190 double const dB = pVdY->at(idx) - pVdX->at(idx) * dM;
191
192 // VdYNew[i] = (pVdXNew->at(i) * dM) + dB;
193 VdYNew[nXNewSize - 1 - i] = (pVdXNew->at(i) * dM) + dB;
194 }
195
196 return VdYNew;
197}
vector< double > VdInterpolateCShoreProfileOutput(vector< double > const *, vector< double > const *, vector< double > const *)
Returns a linearly interpolated vector of doubles, to make CShore profile output compatible with CME....
int nFindIndex(vector< double > const *, double const)
This is used by VdInterpolateCShoreProfileOutput, it returns the index of the value in pVdX which is ...
double dGetInterpolatedValue(vector< double > const *, vector< double > const *, double, bool)
This file contains global definitions for CoastalME.
double const TOLERANCE
Definition cme.h:823
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1303
T tAbs(T a)
Definition cme.h:1277
Contains CSimulation definitions.