56 strTexture =
"coarse";
59 double dStillToErodeOnPolygon = dErosionTargetOnPolygon;
74 LogStream <<
m_ulIter <<
": " <<
ERR <<
"while eroding unconsolidated " + strTexture +
" sediment on polygon " << pPolygon->
nGetPolygonCoastID() <<
", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile <<
") for depth of closure = " <<
m_dDepthOfClosure << endl;
89 int const nUpCoastPartProfileLen = nIndex + 1;
94 vector<CGeom2DIPoint> PtiVUpCoastPartProfileCell;
98 for (
int n = nUpCoastPartProfileLen - 1; n >= 0; n--)
101 PtiVUpCoastPartProfileCell.push_back(Pti);
104 int const nUpCoastProfileCoastPoint = pUpCoastProfile->
nGetCoastPoint();
105 int const nDownCoastProfileCoastPoint = pDownCoastProfile->
nGetCoastPoint();
106 int const nXUpCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetX();
107 int const nYUpCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetY();
111 vector<int> nVCoastPoint;
113 if (nDownCoastProfileCoastPoint ==
m_VCoast[nCoast].nGetCoastlineSize() - 1)
116 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
118 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
119 nVCoastPoint.push_back(nCoastPoint);
125 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
127 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
128 nVCoastPoint.push_back(nCoastPoint);
132 double const dAllTargetPerProfile = dErosionTargetOnPolygon / nCoastSegLen;
135 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen,
m_Rand[1]);
138 for (
int n = 0; n < nCoastSegLen; n++)
141 int const nCoastPoint = nVCoastPoint[n];
144 int const nCoastX = PtiCoastPoint.
nGetX();
145 int const nCoastY = PtiCoastPoint.
nGetY();
161 int const nXOffset = nCoastX - PtiVUpCoastPartProfileCell.back().nGetX();
162 int const nYOffset = nCoastY - PtiVUpCoastPartProfileCell.back().nGetY();
165 vector<CGeom2DIPoint> VPtiParProfile;
167 for (
int m = 0; m < nUpCoastPartProfileLen; m++)
170 CGeom2DIPoint const PtiTmp(PtiVUpCoastPartProfileCell[m].nGetX() + nXOffset, PtiVUpCoastPartProfileCell[m].nGetY() + nYOffset);
171 VPtiParProfile.push_back(PtiTmp);
175 int nParProfEndX = VPtiParProfile[0].nGetX();
176 int nParProfEndY = VPtiParProfile[0].nGetY();
185 VPtiParProfile[0].SetX(nParProfEndX);
186 VPtiParProfile[0].SetY(nParProfEndY);
189 bool bHitEdge =
false;
190 bool bEndProfile =
false;
191 bool bZeroGradient =
false;
192 bool bEnoughEroded =
false;
195 int nInlandOffset = -1;
197 double const dParProfCoastElev =
m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetSedimentTopElev();
198 double const dParProfEndElev =
m_pRasterGrid->m_Cell[nParProfEndX][nParProfEndY].dGetSedimentTopElev();
200 vector<double> VdParProfileDeanElev;
203 vector<int> VnParProfLenEachOffset;
204 vector<double> VdAmountEachOffset;
205 vector<vector<CGeom2DIPoint>> VVPtiParProfileEachOffset;
206 vector<vector<double>> VVdParProfileDeanElevEachOffset;
214 if (nInlandOffset > 0)
219 LogStream <<
m_ulIter <<
": reached end of up-coast profile " << nUpCoastProfile <<
" during down-coast erosion of unconsolidated " + strTexture +
" sediment for coast " << nCoast <<
" polygon " << pPolygon->
nGetPolygonCoastID() <<
" (nCoastPoint = " << nCoastPoint <<
" nInlandOffset = " << nInlandOffset <<
")" << endl;
229 int const nXUpCoastStartOffset = PtiUpCoastTmp.
nGetX() - nXUpCoastProfileExistingCoastPoint;
230 int const nYUpCoastStartOffset = PtiUpCoastTmp.
nGetY() - nYUpCoastProfileExistingCoastPoint;
231 int const nXUpCoastThisStart = nCoastX - nXUpCoastStartOffset;
232 int const nYUpCoastThisStart = nCoastY - nYUpCoastStartOffset;
249 nXParNew = nXUpCoastThisStart + nXOffset,
250 nYParNew = nYUpCoastThisStart + nYOffset;
262 if ((VPtiParProfile.back().nGetX() != nXParNew) || (VPtiParProfile.back().nGetY() != nYParNew))
266 VPtiParProfile.push_back(PtiTmp);
274 VPtiParProfile.push_back(PtiTmp);
278 nParProfLen =
static_cast<int>(VPtiParProfile.size());
300 double const dElevDiff = dParProfCoastElev - dParProfEndElev;
308 bZeroGradient =
true;
313 if (dParProfileLen <= 0)
317 double const dParProfA = dElevDiff / pow(dParProfileLen,
DEAN_POWER);
318 VdParProfileDeanElev.resize(nParProfLen, 0);
319 double const dInc = dParProfileLen / (nParProfLen - 1);
322 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfCoastElev, dParProfA,
false, 0, 0);
324 vector<double> dVParProfileNow(nParProfLen, 0);
325 vector<bool> bVProfileValid(nParProfLen,
true);
327 for (
int m = 0; m < nParProfLen; m++)
329 int const nX = VPtiParProfile[nParProfLen - m - 1].nGetX();
330 int const nY = VPtiParProfile[nParProfLen - m - 1].nGetY();
338 bVProfileValid[m] =
false;
344 bVProfileValid[m] =
false;
346 dVParProfileNow[m] =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
350 double const dParProfTotDiff =
dSubtractProfiles(&dVParProfileNow, &VdParProfileDeanElev, &bVProfileValid);
389 if (dParProfTotDiff > dAllTargetPerProfile)
393 bEnoughEroded =
true;
401 LogStream <<
m_ulIter <<
": leaving loop because nInlandOffset (" << nInlandOffset <<
") >= MIN_INLAND_OFFSET_UNCONS_EROSION) and dParProfTotDiff = " << dParProfTotDiff << endl;
407 VdAmountEachOffset.push_back(dParProfTotDiff);
408 VnParProfLenEachOffset.push_back(nParProfLen);
409 VVPtiParProfileEachOffset.push_back(VPtiParProfile);
410 VVdParProfileDeanElevEachOffset.push_back(VdParProfileDeanElev);
415 if (bHitEdge || bEndProfile || bZeroGradient)
419 double dStillToErodeOnProfile = dAllTargetPerProfile;
424 int nOffsetForLargestPossible = -1;
425 double dLargestPossibleErosion = 0;
427 for (
unsigned int nn = 0; nn < VdAmountEachOffset.size(); nn++)
429 if (VdAmountEachOffset[nn] > dLargestPossibleErosion)
431 dLargestPossibleErosion = VdAmountEachOffset[nn];
432 nOffsetForLargestPossible = nn;
437 if (nOffsetForLargestPossible < 0)
441 nInlandOffset = nOffsetForLargestPossible;
442 dStillToErodeOnProfile = dLargestPossibleErosion;
443 nParProfLen = VnParProfLenEachOffset[nInlandOffset];
444 VPtiParProfile = VVPtiParProfileEachOffset[nInlandOffset];
445 VdParProfileDeanElev = VVdParProfileDeanElevEachOffset[nInlandOffset];
451 int const nRet =
nDoParallelProfileUnconsErosion(pPolygon, nCoastPoint, nCoastX, nCoastY, nTexture, nInlandOffset, nParProfLen, &VPtiParProfile, &VdParProfileDeanElev, dStillToErodeOnProfile, dStillToErodeOnPolygon, dEroded);
466 if (dStillToErodeOnPolygon <= 0)
502int CSimulation::nDoParallelProfileUnconsErosion(
CGeomCoastPolygon* pPolygon,
int const nCoastPoint,
int const nCoastX,
int const nCoastY,
int const nTexture,
int const nInlandOffset,
int const nParProfLen, vector<CGeom2DIPoint>
const* pVPtiParProfile, vector<double>
const* pVdParProfileDeanElev,
double& dStillToErodeOnProfile,
double& dStillToErodeOnPolygon,
double& dTotEroded)
504 for (
int nDistSeawardFromNewCoast = 0; nDistSeawardFromNewCoast < nParProfLen; nDistSeawardFromNewCoast++)
507 if (dStillToErodeOnPolygon <= 0)
510 LogStream <<
m_ulIter <<
": AAA in polygon " << pPolygon->
nGetPolygonCoastID() <<
" at coast point " << nCoastPoint <<
" nInlandOffset = " << nInlandOffset <<
", leaving loop because enough erosion for polygon, dStillToErodeOnPolygon = " << dStillToErodeOnPolygon <<
" dStillToErodeOnProfile = " << dStillToErodeOnProfile << endl;
516 if (dStillToErodeOnProfile <= 0)
519 LogStream <<
m_ulIter <<
": BBB in polygon " << pPolygon->
nGetPolygonCoastID() <<
" at coast point " << nCoastPoint <<
" nInlandOffset = " << nInlandOffset <<
", leaving loop because enough erosion for profile, dStillToErodeOnPolygon = " << dStillToErodeOnPolygon <<
" dStillToErodeOnProfile = " << dStillToErodeOnProfile << endl;
524 CGeom2DIPoint PtiTmp = pVPtiParProfile->at(nParProfLen - nDistSeawardFromNewCoast - 1);
544 if (!
m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
547 double const dThisElevNow =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
552 double const dElevDiff = dThisElevNow - pVdParProfileDeanElev->at(nDistSeawardFromNewCoast);
554 if ((dElevDiff > 0) && (
m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()))
560 int const nThisLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
569 double const dToErode =
tMin(dElevDiff, dStillToErodeOnProfile, dStillToErodeOnPolygon);
580 dTotEroded += dRemoved;
581 dStillToErodeOnProfile -= dRemoved;
582 dStillToErodeOnPolygon -= dRemoved;
610 else if ((dElevDiff < 0) && (
m_pRasterGrid->m_Cell[nX][nY].bIsInContiguousSea()))
617 double dTotToDeposit =
tMin(-dElevDiff, dTotEroded);
619 int const nTopLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopLayerAboveBasement();
625 if (dTotToDeposit > 0)
627 dTotToDeposit =
tMin(dTotToDeposit, dTotEroded);
631 double const dSandNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
632 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dTotToDeposit);
637 dTotEroded -= dTotToDeposit;
639 dStillToErodeOnProfile += dTotToDeposit;
640 dStillToErodeOnPolygon += dTotToDeposit;
649 double const dCoarseNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
650 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dTotToDeposit);
655 dTotEroded -= dTotToDeposit;
657 dStillToErodeOnProfile += dTotToDeposit;
658 dStillToErodeOnPolygon += dTotToDeposit;
673 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dTotToDeposit);
785 strTexture =
"coarse";
805 LogStream <<
m_ulIter <<
": " <<
ERR <<
"while depositing " + strTexture +
" unconsolidated sediment for coast " << nCoast <<
" polygon " << pPolygon->
nGetPolygonCoastID() <<
", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile <<
") for depth of closure = " <<
m_dDepthOfClosure << endl;
820 int const nUpCoastDeanLen = nIndex + 1;
831 int const nUpCoastProfileCoastPoint = pUpCoastProfile->
nGetCoastPoint();
832 int const nDownCoastProfileCoastPoint = pDownCoastProfile->
nGetCoastPoint();
833 int const nXUpCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetX();
834 int const nYUpCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nUpCoastProfileCoastPoint)->nGetY();
838 vector<int> nVCoastPoint;
840 if (nDownCoastProfileCoastPoint ==
m_VCoast[nCoast].nGetCoastlineSize() - 1)
843 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
845 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
846 nVCoastPoint.push_back(nCoastPoint);
852 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
854 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
855 nVCoastPoint.push_back(nCoastPoint);
858 double dStillToDepositOnPoly = dTargetToDepositOnPoly;
859 double dTargetToDepositOnProfile = dTargetToDepositOnPoly / nCoastSegLen;
860 double dStillToDepositOnProfile;
863 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen,
m_Rand[1]);
904 for (
int n = 0; n < nCoastSegLen; n++)
907 int const nCoastPoint = nVCoastPoint[n];
909 int const nCoastX = PtiCoastPoint.
nGetX();
910 int const nCoastY = PtiCoastPoint.
nGetY();
913 int const nXOffset = nCoastX - nXUpCoastProfileExistingCoastPoint;
914 int const nYOffset = nCoastY - nYUpCoastProfileExistingCoastPoint;
915 int nSeawardOffset = -1;
917 vector<CGeom2DIPoint> PtiVParProfile;
918 vector<double> VdParProfileDeanElev;
927 int nParProfLen = nUpCoastDeanLen + nSeawardOffset;
939 PtiVParProfile.resize(0);
941 for (
int m = 0; m < nParProfLen; m++)
945 PtiVParProfile.push_back(PtiTmp);
949 int nSeaEndX = PtiVParProfile.back().nGetX();
950 int nSeaEndY = PtiVParProfile.back().nGetY();
959 PtiVParProfile.back().SetX(nSeaEndX);
960 PtiVParProfile.back().SetY(nSeaEndY);
963 double const dParProfEndElev =
m_pRasterGrid->m_Cell[nSeaEndX][nSeaEndY].dGetSedimentTopElev();
969 double const dParProfLen =
dGetDistanceBetween(&PtiVParProfile.front(), &PtiVParProfile.back());
972 double dParProfDeanLen = dParProfLen - (nSeawardOffset *
m_dCellSide);
975 if (dParProfDeanLen <= 0)
979 double const dParProfA = (dParProfStartElev - dParProfEndElev) / pow(dParProfDeanLen,
DEAN_POWER);
981 nParProfLen =
static_cast<int>(PtiVParProfile.size());
982 VdParProfileDeanElev.resize(nParProfLen, 0);
988 double const dInc = dParProfDeanLen / (nParProfLen - nSeawardOffset - 2);
991 double const dCoastElev =
m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetSedimentTopElev();
994 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfStartElev, dParProfA,
true, nSeawardOffset, dCoastElev);
996 double dParProfTotDiff = 0;
998 for (
int m = 0; m < nParProfLen; m++)
1001 int nX = PtiTmp.
nGetX();
1002 int nY = PtiTmp.
nGetY();
1019 if (!
m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1021 double const dTmpElev =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
1022 double const dDiff = VdParProfileDeanElev[m] - dTmpElev;
1024 dParProfTotDiff += dDiff;
1075 if (dParProfTotDiff >= dTargetToDepositOnProfile)
1086 double dDepositedOnProfile = 0;
1090 dStillToDepositOnProfile = dTargetToDepositOnProfile;
1092 for (
unsigned int nSeawardFromCoast = 0; nSeawardFromCoast < PtiVParProfile.size(); nSeawardFromCoast++)
1112 int nX = PtiTmp.
nGetX();
1113 int nY = PtiTmp.
nGetY();
1130 if (!
m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1133 double const dThisElevNow =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
1138 double const dElevDiff = VdParProfileDeanElev[nSeawardFromCoast] - dThisElevNow;
1143 bool bDeposited =
false;
1144 double dToDepositHere = 0;
1149 dToDepositHere =
tMin(dElevDiff, dStillToDepositOnProfile, dStillToDepositOnPoly);
1153 int const nTopLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopLayerAboveBasement();
1165 double const dSandNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
1167 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dToDepositHere);
1172 dDepositedOnProfile += dToDepositHere;
1173 dDepositedOnPoly += dToDepositHere;
1178 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1180 dStillToDepositOnPoly -= dToDepositHere;
1181 dStillToDepositOnProfile -= dToDepositHere;
1193 double const dCoarseNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
1195 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dToDepositHere);
1205 dDepositedOnProfile += dToDepositHere;
1206 dDepositedOnPoly += dToDepositHere;
1211 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1213 dStillToDepositOnPoly -= dToDepositHere;
1214 dStillToDepositOnProfile -= dToDepositHere;
1256 m_pRasterGrid->m_Cell[nX][nY].SetPotentialBeachErosion(-dElevDiff);
1261 int const nThisLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1272 double dSandRemoved = 0;
1276 dDepositedOnProfile -= dSandRemoved;
1277 dStillToDepositOnProfile += dSandRemoved;
1280 dDepositedOnPoly -= dSandRemoved;
1281 dStillToDepositOnPoly += dSandRemoved;
1299 double dCoarseRemoved = 0;
1303 dDepositedOnProfile -= dCoarseRemoved;
1304 dStillToDepositOnProfile += dCoarseRemoved;
1307 dDepositedOnPoly -= dCoarseRemoved;
1308 dStillToDepositOnPoly += dCoarseRemoved;
1336 if (dTargetToDepositOnPoly > 0)
1349 LogStream <<
m_ulIter <<
": " <<
ERR <<
"while depositing beach for coast " << nCoast <<
" polygon " << pPolygon->
nGetPolygonCoastID() <<
", could not find the seaward end point of the up-coast profile (" << nUpCoastProfile <<
") for depth of closure = " <<
m_dDepthOfClosure << endl;
1364 int const nDownCoastDeanLen = nIndex1 + 1;
1375 int const nXDownCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nDownCoastProfileCoastPoint)->nGetX();
1376 int const nYDownCoastProfileExistingCoastPoint =
m_VCoast[nCoast].pPtiGetCellMarkedAsCoastline(nDownCoastProfileCoastPoint)->nGetY();
1381 nVCoastPoint.resize(0);
1383 if (nUpCoastProfileCoastPoint == 0)
1386 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint + 1;
1388 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint <= nDownCoastProfileCoastPoint; nCoastPoint++)
1389 nVCoastPoint.push_back(nCoastPoint);
1395 nCoastSegLen = nDownCoastProfileCoastPoint - nUpCoastProfileCoastPoint;
1397 for (
int nCoastPoint = nUpCoastProfileCoastPoint; nCoastPoint < nDownCoastProfileCoastPoint; nCoastPoint++)
1398 nVCoastPoint.push_back(nCoastPoint);
1402 shuffle(nVCoastPoint.begin(), nVCoastPoint.begin() + nCoastSegLen,
m_Rand[1]);
1405 dTargetToDepositOnProfile = dStillToDepositOnPoly / nCoastSegLen;
1408 bool bEnoughDepositedOnPolygon =
false;
1411 for (
int n = 0; n < nCoastSegLen; n++)
1414 if (bEnoughDepositedOnPolygon)
1418 int const nCoastPoint = nVCoastPoint[n];
1420 int const nCoastX = PtiCoastPoint.
nGetX();
1421 int const nCoastY = PtiCoastPoint.
nGetY();
1424 int const nXOffset = nCoastX - nXDownCoastProfileExistingCoastPoint;
1425 int const nYOffset = nCoastY - nYDownCoastProfileExistingCoastPoint;
1426 int nSeawardOffset = -1;
1428 vector<CGeom2DIPoint> PtiVParProfile;
1429 vector<double> VdParProfileDeanElev;
1438 int nParProfLen = nDownCoastDeanLen + nSeawardOffset;
1450 PtiVParProfile.resize(0);
1452 for (
int m = 0; m < nParProfLen; m++)
1456 PtiVParProfile.push_back(PtiTmp);
1460 int nSeaEndX = PtiVParProfile.back().nGetX();
1461 int nSeaEndY = PtiVParProfile.back().nGetY();
1469 PtiVParProfile.back().SetX(nSeaEndX);
1470 PtiVParProfile.back().SetY(nSeaEndY);
1473 double const dParProfEndElev =
m_pRasterGrid->m_Cell[nSeaEndX][nSeaEndY].dGetSedimentTopElev();
1479 double const dParProfLen =
dGetDistanceBetween(&PtiVParProfile.front(), &PtiVParProfile.back());
1482 double dParProfDeanLen = dParProfLen - (nSeawardOffset *
m_dCellSide);
1485 if (dParProfDeanLen <= 0)
1486 dParProfDeanLen = 1;
1489 double const dParProfA = (dParProfStartElev - dParProfEndElev) / pow(dParProfDeanLen,
DEAN_POWER);
1491 nParProfLen =
static_cast<int>(PtiVParProfile.size());
1492 VdParProfileDeanElev.resize(nParProfLen, 0);
1494 double const dInc = dParProfDeanLen / (nParProfLen - nSeawardOffset - 2);
1497 double const dCoastElev =
m_pRasterGrid->m_Cell[nCoastX][nCoastY].dGetSedimentTopElev();
1500 CalcDeanProfile(&VdParProfileDeanElev, dInc, dParProfStartElev, dParProfA,
true, nSeawardOffset, dCoastElev);
1502 double dParProfTotDiff = 0;
1504 for (
int m = 0; m < nParProfLen; m++)
1507 int nX = PtiTmp.
nGetX();
1508 int nY = PtiTmp.
nGetY();
1525 if (!
m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1527 double const dTmpElev =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
1528 double const dDiff = VdParProfileDeanElev[m] - dTmpElev;
1530 dParProfTotDiff += dDiff;
1581 if (dParProfTotDiff >= dTargetToDepositOnProfile)
1592 double dDepositedOnProfile = 0;
1594 dStillToDepositOnProfile = dTargetToDepositOnProfile;
1596 for (
unsigned int nSeawardFromCoast = 0; nSeawardFromCoast < PtiVParProfile.size(); nSeawardFromCoast++)
1603 bEnoughDepositedOnPolygon =
true;
1618 int nX = PtiTmp.
nGetX();
1619 int nY = PtiTmp.
nGetY();
1636 if (!
m_pRasterGrid->m_Cell[nX][nY].bBeachErosionOrDepositionThisIter())
1639 double const dThisElevNow =
m_pRasterGrid->m_Cell[nX][nY].dGetSedimentTopElev();
1644 double const dElevDiff = VdParProfileDeanElev[nSeawardFromCoast] - dThisElevNow;
1649 bool bDeposited =
false;
1650 double dToDepositHere = 0;
1655 dToDepositHere =
tMin(dElevDiff, dStillToDepositOnProfile, dStillToDepositOnPoly);
1659 int const nTopLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopLayerAboveBasement();
1671 double const dSandNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetSandDepth();
1673 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetSandDepth(dSandNow + dToDepositHere);
1678 dDepositedOnProfile += dToDepositHere;
1679 dDepositedOnPoly += dToDepositHere;
1684 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1686 dStillToDepositOnPoly -= dToDepositHere;
1687 dStillToDepositOnProfile -= dToDepositHere;
1699 double const dCoarseNow =
m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->dGetCoarseDepth();
1701 m_pRasterGrid->m_Cell[nX][nY].pGetLayerAboveBasement(nTopLayer)->pGetUnconsolidatedSediment()->SetCoarseDepth(dCoarseNow + dToDepositHere);
1706 if (dDepositedOnPoly > dTargetToDepositOnPoly)
1711 dDepositedOnProfile += dToDepositHere;
1712 dDepositedOnPoly += dToDepositHere;
1717 m_pRasterGrid->m_Cell[nX][nY].IncrBeachDeposition(dToDepositHere);
1719 dStillToDepositOnPoly -= dToDepositHere;
1720 dStillToDepositOnProfile -= dToDepositHere;
1756 m_pRasterGrid->m_Cell[nX][nY].SetPotentialBeachErosion(-dElevDiff);
1761 int const nThisLayer =
m_pRasterGrid->m_Cell[nX][nY].nGetTopNonZeroLayerAboveBasement();
1772 double dSandRemoved = 0;
1776 dDepositedOnProfile -= dSandRemoved;
1777 dStillToDepositOnProfile += dSandRemoved;
1780 dDepositedOnPoly -= dSandRemoved;
1781 dStillToDepositOnPoly += dSandRemoved;
1799 double dCoarseRemoved = 0;
1803 dDepositedOnProfile -= dCoarseRemoved;
1804 dStillToDepositOnProfile += dCoarseRemoved;
1807 dDepositedOnPoly -= dCoarseRemoved;
1808 dStillToDepositOnPoly += dCoarseRemoved;