Rev 639 | Rev 667 | Go to most recent revision | Only display areas with differences | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 639 | Rev 664 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp': |
2 | #include "global.h"
|
2 | #include "global.h"
|
3 | #include "math.h"
|
3 | #include "math.h"
|
4 | #include "tree.h"
|
4 | #include "tree.h"
|
5 | 5 | ||
6 | #include "grid.h"
|
6 | #include "grid.h"
|
7 | 7 | ||
8 | #include "stamp.h"
|
8 | #include "stamp.h"
|
9 | #include "species.h"
|
9 | #include "species.h"
|
10 | #include "resourceunit.h"
|
10 | #include "resourceunit.h"
|
11 | #include "model.h"
|
11 | #include "model.h"
|
12 | #include "snag.h"
|
12 | #include "snag.h"
|
13 | 13 | ||
14 | // static varaibles
|
14 | // static varaibles
|
15 | FloatGrid *Tree::mGrid = 0; |
15 | FloatGrid *Tree::mGrid = 0; |
16 | HeightGrid *Tree::mHeightGrid = 0; |
16 | HeightGrid *Tree::mHeightGrid = 0; |
17 | int Tree::m_statPrint=0; |
17 | int Tree::m_statPrint=0; |
18 | int Tree::m_statAboveZ=0; |
18 | int Tree::m_statAboveZ=0; |
19 | int Tree::m_statCreated=0; |
19 | int Tree::m_statCreated=0; |
20 | int Tree::m_nextId=0; |
20 | int Tree::m_nextId=0; |
21 | 21 | ||
22 | 22 | ||
23 | 23 | ||
24 | /** get distance and direction between two points.
|
24 | /** get distance and direction between two points.
|
25 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */
|
25 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */
|
26 | float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle) |
26 | float dist_and_direction(const QPointF &PStart, const QPointF &PEnd, float &rAngle) |
27 | {
|
27 | {
|
28 | float dx = PEnd.x() - PStart.x(); |
28 | float dx = PEnd.x() - PStart.x(); |
29 | float dy = PEnd.y() - PStart.y(); |
29 | float dy = PEnd.y() - PStart.y(); |
30 | float d = sqrt(dx*dx + dy*dy); |
30 | float d = sqrt(dx*dx + dy*dy); |
31 | // direction:
|
31 | // direction:
|
32 | rAngle = atan2(dx, dy); |
32 | rAngle = atan2(dx, dy); |
33 | return d; |
33 | return d; |
34 | }
|
34 | }
|
35 | 35 | ||
36 | 36 | ||
37 | // lifecycle
|
37 | // lifecycle
|
38 | Tree::Tree() |
38 | Tree::Tree() |
39 | {
|
39 | {
|
40 | mDbh = mHeight = 0; |
40 | mDbh = mHeight = 0; |
41 | mRU = 0; mSpecies = 0; |
41 | mRU = 0; mSpecies = 0; |
42 | mFlags = mAge = 0; |
42 | mFlags = mAge = 0; |
43 | mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.; |
43 | mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.; |
44 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
44 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
45 | mLightResponse = 0.; |
45 | mLightResponse = 0.; |
46 | mId = m_nextId++; |
46 | mId = m_nextId++; |
47 | m_statCreated++; |
47 | m_statCreated++; |
48 | }
|
48 | }
|
49 | 49 | ||
50 | float Tree::crownRadius() const |
50 | float Tree::crownRadius() const |
51 | {
|
51 | {
|
52 | Q_ASSERT(mStamp!=0); |
52 | Q_ASSERT(mStamp!=0); |
53 | return mStamp->crownRadius(); |
53 | return mStamp->crownRadius(); |
54 | }
|
54 | }
|
55 | 55 | ||
56 | float Tree::biomassBranch() const |
56 | float Tree::biomassBranch() const |
57 | {
|
57 | {
|
58 | return mSpecies->biomassBranch(mDbh); |
58 | return mSpecies->biomassBranch(mDbh); |
59 | }
|
59 | }
|
60 | 60 | ||
61 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
61 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
62 | {
|
62 | {
|
63 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
63 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
64 | }
|
64 | }
|
65 | 65 | ||
66 | /// dumps some core variables of a tree to a string.
|
66 | /// dumps some core variables of a tree to a string.
|
67 | QString Tree::dump() |
67 | QString Tree::dump() |
68 | {
|
68 | {
|
69 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
69 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
70 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
70 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
71 | .arg(position().x()).arg(position().y()) |
71 | .arg(position().x()).arg(position().y()) |
72 | .arg(mRU->index()).arg(mLRI); |
72 | .arg(mRU->index()).arg(mLRI); |
73 | return result; |
73 | return result; |
74 | }
|
74 | }
|
75 | 75 | ||
76 | void Tree::dumpList(DebugList &rTargetList) |
76 | void Tree::dumpList(DebugList &rTargetList) |
77 | {
|
77 | {
|
78 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
78 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
79 | << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea; |
79 | << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea; |
80 | }
|
80 | }
|
81 | 81 | ||
82 | void Tree::setup() |
82 | void Tree::setup() |
83 | {
|
83 | {
|
84 | if (mDbh<=0 || mHeight<=0) |
84 | if (mDbh<=0 || mHeight<=0) |
85 | return; |
85 | return; |
86 | // check stamp
|
86 | // check stamp
|
87 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
87 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
88 | mStamp = species()->stamp(mDbh, mHeight); |
88 | mStamp = species()->stamp(mDbh, mHeight); |
89 | if (!mStamp) { |
89 | if (!mStamp) { |
90 | throw IException("Tree::setup() with invalid stamp!"); |
90 | throw IException("Tree::setup() with invalid stamp!"); |
91 | }
|
91 | }
|
92 | 92 | ||
93 | mFoliageMass = species()->biomassFoliage(mDbh); |
93 | mFoliageMass = species()->biomassFoliage(mDbh); |
94 | mCoarseRootMass = species()->biomassRoot(mDbh); // coarse root (allometry) |
94 | mCoarseRootMass = species()->biomassRoot(mDbh); // coarse root (allometry) |
95 | mFineRootMass = mFoliageMass * species()->finerootFoliageRatio(); // fine root (size defined by finerootFoliageRatio) |
95 | mFineRootMass = mFoliageMass * species()->finerootFoliageRatio(); // fine root (size defined by finerootFoliageRatio) |
96 | mWoodyMass = species()->biomassWoody(mDbh); |
96 | mWoodyMass = species()->biomassWoody(mDbh); |
97 | 97 | ||
98 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
|
98 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
|
99 | mLeafArea = mFoliageMass * species()->specificLeafArea(); |
99 | mLeafArea = mFoliageMass * species()->specificLeafArea(); |
100 | mOpacity = 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea()); |
100 | mOpacity = 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea()); |
101 | mNPPReserve = (1+species()->finerootFoliageRatio())*mFoliageMass; // initial value |
101 | mNPPReserve = (1+species()->finerootFoliageRatio())*mFoliageMass; // initial value |
102 | mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment |
102 | mDbhDelta = 0.1; // initial value: used in growth() to estimate diameter increment |
103 | 103 | ||
104 | }
|
104 | }
|
105 | 105 | ||
106 | void Tree::setAge(const int age, const float treeheight) |
106 | void Tree::setAge(const int age, const float treeheight) |
107 | {
|
107 | {
|
108 | mAge = age; |
108 | mAge = age; |
109 | if (age==0) { |
109 | if (age==0) { |
110 | // estimate age using the tree height
|
110 | // estimate age using the tree height
|
111 | mAge = mSpecies->estimateAge(treeheight); |
111 | mAge = mSpecies->estimateAge(treeheight); |
112 | }
|
112 | }
|
113 | }
|
113 | }
|
114 | 114 | ||
115 | //////////////////////////////////////////////////
|
115 | //////////////////////////////////////////////////
|
116 | //// Light functions (Pattern-stuff)
|
116 | //// Light functions (Pattern-stuff)
|
117 | //////////////////////////////////////////////////
|
117 | //////////////////////////////////////////////////
|
118 | 118 | ||
119 | #define NOFULLDBG
|
119 | #define NOFULLDBG
|
120 | //#define NOFULLOPT
|
120 | //#define NOFULLOPT
|
121 | 121 | ||
122 | 122 | ||
123 | void Tree::applyLIP() |
123 | void Tree::applyLIP() |
124 | {
|
124 | {
|
125 | if (!mStamp) |
125 | if (!mStamp) |
126 | return; |
126 | return; |
127 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
127 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
128 | QPoint pos = mPositionIndex; |
128 | QPoint pos = mPositionIndex; |
129 | int offset = mStamp->offset(); |
129 | int offset = mStamp->offset(); |
130 | pos-=QPoint(offset, offset); |
130 | pos-=QPoint(offset, offset); |
131 | 131 | ||
132 | float local_dom; // height of Z* on the current position |
132 | float local_dom; // height of Z* on the current position |
133 | int x,y; |
133 | int x,y; |
134 | float value, z, z_zstar; |
134 | float value, z, z_zstar; |
135 | int gr_stamp = mStamp->size(); |
135 | int gr_stamp = mStamp->size(); |
136 | int grid_x, grid_y; |
136 | int grid_x, grid_y; |
137 | float *grid_value_ptr; |
137 | float *grid_value_ptr; |
138 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
138 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
139 | // this should not happen because of the buffer
|
139 | // this should not happen because of the buffer
|
140 | return; |
140 | return; |
141 | }
|
141 | }
|
142 | grid_y = pos.y(); |
142 | grid_y = pos.y(); |
143 | for (y=0;y<gr_stamp; ++y) { |
143 | for (y=0;y<gr_stamp; ++y) { |
144 | 144 | ||
145 | grid_value_ptr = mGrid->ptr(pos.x(), grid_y); |
145 | grid_value_ptr = mGrid->ptr(pos.x(), grid_y); |
146 | grid_x = pos.x(); |
146 | grid_x = pos.x(); |
147 | for (x=0;x<gr_stamp;++x) { |
147 | for (x=0;x<gr_stamp;++x) { |
148 | // suppose there is no stamping outside
|
148 | // suppose there is no stamping outside
|
149 | 149 | ||
150 | local_dom = mHeightGrid->valueAtIndex(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height; |
150 | local_dom = mHeightGrid->valueAtIndex(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height; |
151 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
151 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
152 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
152 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
153 | value = (*mStamp)(x,y); // stampvalue |
153 | value = (*mStamp)(x,y); // stampvalue |
154 | value = 1. - value*mOpacity * z_zstar; // calculated value |
154 | value = 1. - value*mOpacity * z_zstar; // calculated value |
155 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
155 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
156 | value = qMax(value, 0.02f); // limit value |
156 | value = qMax(value, 0.02f); // limit value |
157 | 157 | ||
158 | *grid_value_ptr++ *= value; |
158 | *grid_value_ptr++ *= value; |
159 | 159 | ||
160 | grid_x++; |
160 | grid_x++; |
161 | }
|
161 | }
|
162 | grid_y++; |
162 | grid_y++; |
163 | }
|
163 | }
|
164 | 164 | ||
165 | m_statPrint++; // count # of stamp applications... |
165 | m_statPrint++; // count # of stamp applications... |
166 | }
|
166 | }
|
167 | 167 | ||
168 | /// helper function for gluing the edges together
|
168 | /// helper function for gluing the edges together
|
169 | /// index: index at grid
|
169 | /// index: index at grid
|
170 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)
|
170 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)
|
171 | /// buffer: size of buffer around simulation area (in pixels)
|
171 | /// buffer: size of buffer around simulation area (in pixels)
|
172 | inline int torusIndex(int index, int count, int buffer, int ru_index) |
172 | inline int torusIndex(int index, int count, int buffer, int ru_index) |
173 | {
|
173 | {
|
174 | return buffer + ru_index + (index-buffer+count)%count; |
174 | return buffer + ru_index + (index-buffer+count)%count; |
175 | }
|
175 | }
|
176 | 176 | ||
177 | 177 | ||
178 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.
|
178 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.
|
179 | */
|
179 | */
|
180 | void Tree::applyLIP_torus() |
180 | void Tree::applyLIP_torus() |
181 | {
|
181 | {
|
182 | if (!mStamp) |
182 | if (!mStamp) |
183 | return; |
183 | return; |
184 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
184 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
185 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
185 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
186 | QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
186 | QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
187 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
187 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
188 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index |
188 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index |
189 | 189 | ||
190 | int offset = mStamp->offset(); |
190 | int offset = mStamp->offset(); |
191 | pos-=QPoint(offset, offset); |
191 | pos-=QPoint(offset, offset); |
192 | 192 | ||
193 | float local_dom; // height of Z* on the current position |
193 | float local_dom; // height of Z* on the current position |
194 | int x,y; |
194 | int x,y; |
195 | float value; |
195 | float value; |
196 | int gr_stamp = mStamp->size(); |
196 | int gr_stamp = mStamp->size(); |
197 | int grid_x, grid_y; |
197 | int grid_x, grid_y; |
198 | float *grid_value; |
198 | float *grid_value; |
199 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
199 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
200 | // todo: in this case we should use another algorithm!!! necessary????
|
200 | // todo: in this case we should use another algorithm!!! necessary????
|
201 | return; |
201 | return; |
202 | }
|
202 | }
|
203 | float z, z_zstar; |
203 | float z, z_zstar; |
204 | int xt, yt; // wraparound coordinates |
204 | int xt, yt; // wraparound coordinates |
205 | for (y=0;y<gr_stamp; ++y) { |
205 | for (y=0;y<gr_stamp; ++y) { |
206 | grid_y = pos.y() + y; |
206 | grid_y = pos.y() + y; |
207 | yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m |
207 | yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m |
208 | for (x=0;x<gr_stamp;++x) { |
208 | for (x=0;x<gr_stamp;++x) { |
209 | // suppose there is no stamping outside
|
209 | // suppose there is no stamping outside
|
210 | grid_x = pos.x() + x; |
210 | grid_x = pos.x() + x; |
211 | xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x()); |
211 | xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x()); |
212 | 212 | ||
213 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height; |
213 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height; |
214 | 214 | ||
215 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
215 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
216 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
216 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
217 | value = (*mStamp)(x,y); // stampvalue |
217 | value = (*mStamp)(x,y); // stampvalue |
218 | value = 1. - value*mOpacity * z_zstar; // calculated value |
218 | value = 1. - value*mOpacity * z_zstar; // calculated value |
219 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
219 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
220 | value = qMax(value, 0.02f); // limit value |
220 | value = qMax(value, 0.02f); // limit value |
221 | 221 | ||
222 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
222 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
223 | *grid_value *= value; |
223 | *grid_value *= value; |
224 | }
|
224 | }
|
225 | }
|
225 | }
|
226 | 226 | ||
227 | m_statPrint++; // count # of stamp applications... |
227 | m_statPrint++; // count # of stamp applications... |
228 | }
|
228 | }
|
229 | 229 | ||
230 | /** heightGrid()
|
230 | /** heightGrid()
|
231 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.
|
231 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.
|
232 | */
|
232 | */
|
233 | void Tree::heightGrid() |
233 | void Tree::heightGrid() |
234 | {
|
234 | {
|
235 | 235 | ||
236 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
236 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
237 | 237 | ||
238 | // count trees that are on height-grid cells (used for stockable area)
|
238 | // count trees that are on height-grid cells (used for stockable area)
|
239 | mHeightGrid->valueAtIndex(p).increaseCount(); |
239 | mHeightGrid->valueAtIndex(p).increaseCount(); |
240 | if (mHeight > mHeightGrid->valueAtIndex(p).height) |
240 | if (mHeight > mHeightGrid->valueAtIndex(p).height) |
241 | mHeightGrid->valueAtIndex(p).height=mHeight; |
241 | mHeightGrid->valueAtIndex(p).height=mHeight; |
242 | 242 | ||
243 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
243 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
244 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
244 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
245 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
245 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
246 | if (index_eastwest - r < 0) { // east |
246 | if (index_eastwest - r < 0) { // east |
247 | mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight); |
247 | mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight); |
248 | }
|
248 | }
|
249 | if (index_eastwest + r >= cPxPerHeight) { // west |
249 | if (index_eastwest + r >= cPxPerHeight) { // west |
250 | mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight); |
250 | mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight); |
251 | }
|
251 | }
|
252 | if (index_northsouth - r < 0) { // south |
252 | if (index_northsouth - r < 0) { // south |
253 | mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight); |
253 | mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight); |
254 | }
|
254 | }
|
255 | if (index_northsouth + r >= cPxPerHeight) { // north |
255 | if (index_northsouth + r >= cPxPerHeight) { // north |
256 | mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight); |
256 | mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight); |
257 | }
|
257 | }
|
258 | 258 | ||
259 | 259 | ||
260 | // without spread of the height grid
|
260 | // without spread of the height grid
|
261 | 261 | ||
262 | // // height of Z*
|
262 | // // height of Z*
|
263 | // const float cellsize = mHeightGrid->cellsize();
|
263 | // const float cellsize = mHeightGrid->cellsize();
|
264 | //
|
264 | //
|
265 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
265 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
266 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
266 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
267 | // int dist[9];
|
267 | // int dist[9];
|
268 | // dist[3] = index_northsouth * 2 + 1; // south
|
268 | // dist[3] = index_northsouth * 2 + 1; // south
|
269 | // dist[1] = index_eastwest * 2 + 1; // west
|
269 | // dist[1] = index_eastwest * 2 + 1; // west
|
270 | // dist[5] = 10 - dist[3]; // north
|
270 | // dist[5] = 10 - dist[3]; // north
|
271 | // dist[7] = 10 - dist[1]; // east
|
271 | // dist[7] = 10 - dist[1]; // east
|
272 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
272 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
273 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
273 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
274 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
274 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
275 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
275 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
276 | // dist[4] = 0; // center cell
|
276 | // dist[4] = 0; // center cell
|
277 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
277 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
278 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
278 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
279 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
279 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
280 | // */
|
280 | // */
|
281 | //
|
281 | //
|
282 | //
|
282 | //
|
283 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
283 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
284 | // int ix, iy;
|
284 | // int ix, iy;
|
285 | // int ring;
|
285 | // int ring;
|
286 | // float hdom;
|
286 | // float hdom;
|
287 | //
|
287 | //
|
288 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
288 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
289 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
289 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
290 | // ring = qMax(abs(ix), abs(iy));
|
290 | // ring = qMax(abs(ix), abs(iy));
|
291 | // QPoint pos(ix+p.x(), iy+p.y());
|
291 | // QPoint pos(ix+p.x(), iy+p.y());
|
292 | // if (mHeightGrid->isIndexValid(pos)) {
|
292 | // if (mHeightGrid->isIndexValid(pos)) {
|
293 | // float &rHGrid = mHeightGrid->valueAtIndex(pos).height;
|
293 | // float &rHGrid = mHeightGrid->valueAtIndex(pos).height;
|
294 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
294 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
295 | // continue;
|
295 | // continue;
|
296 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
296 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
297 | // hdom = mHeight - dist[direction];
|
297 | // hdom = mHeight - dist[direction];
|
298 | // if (ring>1)
|
298 | // if (ring>1)
|
299 | // hdom -= (ring-1)*10;
|
299 | // hdom -= (ring-1)*10;
|
300 | //
|
300 | //
|
301 | // rHGrid = qMax(rHGrid, hdom); // write value
|
301 | // rHGrid = qMax(rHGrid, hdom); // write value
|
302 | // } // is valid
|
302 | // } // is valid
|
303 | // } // for (y)
|
303 | // } // for (y)
|
304 | }
|
304 | }
|
305 | 305 | ||
306 | void Tree::heightGrid_torus() |
306 | void Tree::heightGrid_torus() |
307 | {
|
307 | {
|
308 | // height of Z*
|
308 | // height of Z*
|
309 | 309 | ||
310 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
310 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
311 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels) |
311 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels) |
312 | p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m |
312 | p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m |
313 | p.setY((p.y()-bufferOffset)%10 + bufferOffset); |
313 | p.setY((p.y()-bufferOffset)%10 + bufferOffset); |
314 | 314 | ||
315 | 315 | ||
316 | // torus coordinates: ru_offset = coords of lower left corner of 1ha patch
|
316 | // torus coordinates: ru_offset = coords of lower left corner of 1ha patch
|
317 | QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y()); |
317 | QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y()); |
318 | 318 | ||
319 | // count trees that are on height-grid cells (used for stockable area)
|
319 | // count trees that are on height-grid cells (used for stockable area)
|
320 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
320 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
321 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
321 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
322 | v.increaseCount(); |
322 | v.increaseCount(); |
323 | v.height = qMax(v.height, mHeight); |
323 | v.height = qMax(v.height, mHeight); |
324 | 324 | ||
325 | 325 | ||
326 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
326 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
327 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
327 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
328 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
328 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
329 | if (index_eastwest - r < 0) { // east |
329 | if (index_eastwest - r < 0) { // east |
330 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()), |
330 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()), |
331 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
331 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
332 | v.height = qMax(v.height, mHeight); |
332 | v.height = qMax(v.height, mHeight); |
333 | }
|
333 | }
|
334 | if (index_eastwest + r >= cPxPerHeight) { // west |
334 | if (index_eastwest + r >= cPxPerHeight) { // west |
335 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()), |
335 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()), |
336 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
336 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
337 | v.height = qMax(v.height, mHeight); |
337 | v.height = qMax(v.height, mHeight); |
338 | }
|
338 | }
|
339 | if (index_northsouth - r < 0) { // south |
339 | if (index_northsouth - r < 0) { // south |
340 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
340 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
341 | torusIndex(p.y()-1,10,bufferOffset,ru_offset.y())); |
341 | torusIndex(p.y()-1,10,bufferOffset,ru_offset.y())); |
342 | v.height = qMax(v.height, mHeight); |
342 | v.height = qMax(v.height, mHeight); |
343 | }
|
343 | }
|
344 | if (index_northsouth + r >= cPxPerHeight) { // north |
344 | if (index_northsouth + r >= cPxPerHeight) { // north |
345 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
345 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
346 | torusIndex(p.y()+1,10,bufferOffset,ru_offset.y())); |
346 | torusIndex(p.y()+1,10,bufferOffset,ru_offset.y())); |
347 | v.height = qMax(v.height, mHeight); |
347 | v.height = qMax(v.height, mHeight); |
348 | }
|
348 | }
|
349 | 349 | ||
350 | 350 | ||
351 | 351 | ||
352 | 352 | ||
353 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
353 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
354 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
354 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
355 | // int dist[9];
|
355 | // int dist[9];
|
356 | // dist[3] = index_northsouth * 2 + 1; // south
|
356 | // dist[3] = index_northsouth * 2 + 1; // south
|
357 | // dist[1] = index_eastwest * 2 + 1; // west
|
357 | // dist[1] = index_eastwest * 2 + 1; // west
|
358 | // dist[5] = 10 - dist[3]; // north
|
358 | // dist[5] = 10 - dist[3]; // north
|
359 | // dist[7] = 10 - dist[1]; // east
|
359 | // dist[7] = 10 - dist[1]; // east
|
360 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
360 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
361 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
361 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
362 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
362 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
363 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
363 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
364 | // dist[4] = 0; // center cell
|
364 | // dist[4] = 0; // center cell
|
365 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
365 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
366 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
366 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
367 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
367 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
368 | // */
|
368 | // */
|
369 | //
|
369 | //
|
370 | //
|
370 | //
|
371 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
371 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
372 | // int ix, iy;
|
372 | // int ix, iy;
|
373 | // int ring;
|
373 | // int ring;
|
374 | // float hdom;
|
374 | // float hdom;
|
375 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
375 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
376 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
376 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
377 | // ring = qMax(abs(ix), abs(iy));
|
377 | // ring = qMax(abs(ix), abs(iy));
|
378 | // QPoint pos(ix+p.x(), iy+p.y());
|
378 | // QPoint pos(ix+p.x(), iy+p.y());
|
379 | // QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()),
|
379 | // QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()),
|
380 | // torusIndex(pos.y(),10,bufferOffset,ru_offset.y()));
|
380 | // torusIndex(pos.y(),10,bufferOffset,ru_offset.y()));
|
381 | // if (mHeightGrid->isIndexValid(p_torus)) {
|
381 | // if (mHeightGrid->isIndexValid(p_torus)) {
|
382 | // float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height;
|
382 | // float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height;
|
383 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
383 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
384 | // continue;
|
384 | // continue;
|
385 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
385 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
386 | // hdom = mHeight - dist[direction];
|
386 | // hdom = mHeight - dist[direction];
|
387 | // if (ring>1)
|
387 | // if (ring>1)
|
388 | // hdom -= (ring-1)*10;
|
388 | // hdom -= (ring-1)*10;
|
389 | //
|
389 | //
|
390 | // rHGrid = qMax(rHGrid, hdom); // write value
|
390 | // rHGrid = qMax(rHGrid, hdom); // write value
|
391 | // } // is valid
|
391 | // } // is valid
|
392 | // } // for (y)
|
392 | // } // for (y)
|
393 | }
|
393 | }
|
394 | 394 | ||
395 | 395 | ||
396 | 396 | ||
397 | void Tree::readLIF() |
397 | void Tree::readLIF() |
398 | {
|
398 | {
|
399 | if (!mStamp) |
399 | if (!mStamp) |
400 | return; |
400 | return; |
401 | const Stamp *reader = mStamp->reader(); |
401 | const Stamp *reader = mStamp->reader(); |
402 | if (!reader) |
402 | if (!reader) |
403 | return; |
403 | return; |
404 | QPoint pos_reader = mPositionIndex; |
404 | QPoint pos_reader = mPositionIndex; |
405 | 405 | ||
406 | int offset_reader = reader->offset(); |
406 | int offset_reader = reader->offset(); |
407 | int offset_writer = mStamp->offset(); |
407 | int offset_writer = mStamp->offset(); |
408 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
408 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
409 | 409 | ||
410 | pos_reader-=QPoint(offset_reader, offset_reader); |
410 | pos_reader-=QPoint(offset_reader, offset_reader); |
411 | 411 | ||
412 | float local_dom; |
412 | float local_dom; |
413 | 413 | ||
414 | int x,y; |
414 | int x,y; |
415 | double sum=0.; |
415 | double sum=0.; |
416 | double value, own_value; |
416 | double value, own_value; |
417 | float *grid_value; |
417 | float *grid_value; |
418 | float z, z_zstar; |
418 | float z, z_zstar; |
419 | int reader_size = reader->size(); |
419 | int reader_size = reader->size(); |
420 | int rx = pos_reader.x(); |
420 | int rx = pos_reader.x(); |
421 | int ry = pos_reader.y(); |
421 | int ry = pos_reader.y(); |
422 | for (y=0;y<reader_size; ++y, ++ry) { |
422 | for (y=0;y<reader_size; ++y, ++ry) { |
423 | grid_value = mGrid->ptr(rx, ry); |
423 | grid_value = mGrid->ptr(rx, ry); |
424 | for (x=0;x<reader_size;++x) { |
424 | for (x=0;x<reader_size;++x) { |
425 | 425 | ||
426 | local_dom = mHeightGrid->valueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
426 | local_dom = mHeightGrid->valueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
427 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
427 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
428 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
428 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
429 | 429 | ||
430 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
430 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
431 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
431 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
432 | own_value = qMax(own_value, 0.02); |
432 | own_value = qMax(own_value, 0.02); |
433 | value = *grid_value++ / own_value; // remove impact of focal tree |
433 | value = *grid_value++ / own_value; // remove impact of focal tree |
434 | 434 | ||
435 | //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset);
|
435 | //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset);
|
436 | //if (value>0.)
|
436 | //if (value>0.)
|
437 | sum += value * (*reader)(x,y); |
437 | sum += value * (*reader)(x,y); |
438 | 438 | ||
439 | }
|
439 | }
|
440 | }
|
440 | }
|
441 | mLRI = sum; |
441 | mLRI = sum; |
442 | // LRI correction...
|
442 | // LRI correction...
|
443 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
443 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
444 | if (hrel<1.) |
444 | if (hrel<1.) |
445 | mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel); |
445 | mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel); |
446 | 446 | ||
447 | // read dominance field...
|
447 | // read dominance field...
|
448 | // this applies only if some trees are potentially *higher* than the dominant height grid
|
448 | // this applies only if some trees are potentially *higher* than the dominant height grid
|
449 | //if (dom_height < m_Height) {
|
449 | //if (dom_height < m_Height) {
|
450 | // if tree is higher than Z*, the dominance height
|
450 | // if tree is higher than Z*, the dominance height
|
451 | // a part of the crown is in "full light".
|
451 | // a part of the crown is in "full light".
|
452 | // m_statAboveZ++;
|
452 | // m_statAboveZ++;
|
453 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height;
|
453 | // mImpact = 1. - (1. - mImpact)*dom_height/m_Height;
|
454 | //}
|
454 | //}
|
455 | 455 | ||
456 | if (mLRI > 1.) |
456 | if (mLRI > 1.) |
457 | mLRI = 1.; |
457 | mLRI = 1.; |
458 | 458 | ||
459 | // Finally, add LRI of this Tree to the ResourceUnit!
|
459 | // Finally, add LRI of this Tree to the ResourceUnit!
|
460 | mRU->addWLA(mLeafArea, mLRI); |
460 | mRU->addWLA(mLeafArea, mLRI); |
461 | 461 | ||
462 | 462 | ||
463 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
463 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
464 | //mRU->addWLA(mLRI*mLeafArea, mLeafArea);
|
464 | //mRU->addWLA(mLRI*mLeafArea, mLeafArea);
|
465 | }
|
465 | }
|
466 | 466 | ||
467 | /// Torus version of read stamp (glued edges)
|
467 | /// Torus version of read stamp (glued edges)
|
468 | void Tree::readLIF_torus() |
468 | void Tree::readLIF_torus() |
469 | {
|
469 | {
|
470 | if (!mStamp) |
470 | if (!mStamp) |
471 | return; |
471 | return; |
472 | const Stamp *reader = mStamp->reader(); |
472 | const Stamp *reader = mStamp->reader(); |
473 | if (!reader) |
473 | if (!reader) |
474 | return; |
474 | return; |
475 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
475 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
476 | 476 | ||
477 | QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
477 | QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
478 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
478 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
479 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index |
479 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index |
480 | 480 | ||
481 | int offset_reader = reader->offset(); |
481 | int offset_reader = reader->offset(); |
482 | int offset_writer = mStamp->offset(); |
482 | int offset_writer = mStamp->offset(); |
483 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
483 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
484 | 484 | ||
485 | pos_reader-=QPoint(offset_reader, offset_reader); |
485 | pos_reader-=QPoint(offset_reader, offset_reader); |
486 | 486 | ||
487 | float local_dom; |
487 | float local_dom; |
488 | 488 | ||
489 | int x,y; |
489 | int x,y; |
490 | double sum=0.; |
490 | double sum=0.; |
491 | double value, own_value; |
491 | double value, own_value; |
492 | float *grid_value; |
492 | float *grid_value; |
493 | float z, z_zstar; |
493 | float z, z_zstar; |
494 | int reader_size = reader->size(); |
494 | int reader_size = reader->size(); |
495 | int rx = pos_reader.x(); |
495 | int rx = pos_reader.x(); |
496 | int ry = pos_reader.y(); |
496 | int ry = pos_reader.y(); |
497 | int xt, yt; // wrapped coords |
497 | int xt, yt; // wrapped coords |
498 | 498 | ||
499 | for (y=0;y<reader_size; ++y) { |
499 | for (y=0;y<reader_size; ++y) { |
500 | yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y()); |
500 | yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y()); |
501 | for (x=0;x<reader_size;++x) { |
501 | for (x=0;x<reader_size;++x) { |
502 | xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x()); |
502 | xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x()); |
503 | grid_value = mGrid->ptr(xt,yt); |
503 | grid_value = mGrid->ptr(xt,yt); |
504 | 504 | ||
505 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
505 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
506 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
506 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45° line) |
507 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
507 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
508 | 508 | ||
509 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
509 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
510 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
510 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
511 | own_value = qMax(own_value, 0.02); |
511 | own_value = qMax(own_value, 0.02); |
512 | value = *grid_value++ / own_value; // remove impact of focal tree |
512 | value = *grid_value++ / own_value; // remove impact of focal tree |
513 | 513 | ||
514 | // debug for one tree in HJA
|
514 | // debug for one tree in HJA
|
515 | //if (id()==178020)
|
515 | //if (id()==178020)
|
516 | // qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y);
|
516 | // qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y);
|
517 | //if (_isnan(value))
|
517 | //if (_isnan(value))
|
518 | // qDebug() << "isnan" << id();
|
518 | // qDebug() << "isnan" << id();
|
519 | if (value * (*reader)(x,y)>1.) |
519 | if (value * (*reader)(x,y)>1.) |
520 | qDebug() << "LIFTorus: value>1."; |
520 | qDebug() << "LIFTorus: value>1."; |
521 | sum += value * (*reader)(x,y); |
521 | sum += value * (*reader)(x,y); |
522 | 522 | ||
523 | //} // isIndexValid
|
523 | //} // isIndexValid
|
524 | }
|
524 | }
|
525 | }
|
525 | }
|
526 | mLRI = sum; |
526 | mLRI = sum; |
527 | 527 | ||
528 | // LRI correction...
|
528 | // LRI correction...
|
529 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
529 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
530 | if (hrel<1.) |
530 | if (hrel<1.) |
531 | mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel); |
531 | mLRI = species()->speciesSet()->LRIcorrection(mLRI, hrel); |
532 | 532 | ||
533 | 533 | ||
534 | if (isnan(mLRI)) { |
534 | if (isnan(mLRI)) { |
535 | qDebug() << "LRI invalid (nan)!" << id(); |
535 | qDebug() << "LRI invalid (nan)!" << id(); |
536 | mLRI=0.; |
536 | mLRI=0.; |
537 | //qDebug() << reader->dump();
|
537 | //qDebug() << reader->dump();
|
538 | }
|
538 | }
|
539 | if (mLRI > 1.) |
539 | if (mLRI > 1.) |
540 | mLRI = 1.; |
540 | mLRI = 1.; |
541 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
541 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
542 | 542 | ||
543 | // Finally, add LRI of this Tree to the ResourceUnit!
|
543 | // Finally, add LRI of this Tree to the ResourceUnit!
|
544 | mRU->addWLA(mLeafArea, mLRI); |
544 | mRU->addWLA(mLeafArea, mLRI); |
545 | }
|
545 | }
|
546 | 546 | ||
547 | 547 | ||
548 | void Tree::resetStatistics() |
548 | void Tree::resetStatistics() |
549 | {
|
549 | {
|
550 | m_statPrint=0; |
550 | m_statPrint=0; |
551 | m_statCreated=0; |
551 | m_statCreated=0; |
552 | m_statAboveZ=0; |
552 | m_statAboveZ=0; |
553 | m_nextId=1; |
553 | m_nextId=1; |
554 | }
|
554 | }
|
555 | 555 | ||
556 | void Tree::calcLightResponse() |
556 | void Tree::calcLightResponse() |
557 | {
|
557 | {
|
558 | // calculate a light response from lri:
|
558 | // calculate a light response from lri:
|
559 | // http://iland.boku.ac.at/individual+tree+light+availability
|
559 | // http://iland.boku.ac.at/individual+tree+light+availability
|
560 | double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.); // Eq. (3) |
560 | double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.); // Eq. (3) |
561 | mLightResponse = mSpecies->lightResponse(lri); // Eq. (4) |
561 | mLightResponse = mSpecies->lightResponse(lri); // Eq. (4) |
562 | mRU->addLR(mLeafArea, mLightResponse); |
562 | mRU->addLR(mLeafArea, mLightResponse); |
563 | 563 | ||
564 | }
|
564 | }
|
565 | 565 | ||
566 | //////////////////////////////////////////////////
|
566 | //////////////////////////////////////////////////
|
567 | //// Growth Functions
|
567 | //// Growth Functions
|
568 | //////////////////////////////////////////////////
|
568 | //////////////////////////////////////////////////
|
569 | 569 | ||
570 | /** grow() is the main function of the yearly tree growth.
|
570 | /** grow() is the main function of the yearly tree growth.
|
571 | The main steps are:
|
571 | The main steps are:
|
572 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability
|
572 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability
|
573 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation
|
573 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation
|
574 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???)
|
574 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???)
|
575 | Further activties: * the age of the tree is increased
|
575 | Further activties: * the age of the tree is increased
|
576 | * the mortality sub routine is executed
|
576 | * the mortality sub routine is executed
|
577 | * seeds are produced */
|
577 | * seeds are produced */
|
578 | void Tree::grow() |
578 | void Tree::grow() |
579 | {
|
579 | {
|
580 | TreeGrowthData d;
|
580 | TreeGrowthData d;
|
581 | mAge++; // increase age |
581 | mAge++; // increase age |
582 | // step 1: get "interception area" of the tree individual [m2]
|
582 | // step 1: get "interception area" of the tree individual [m2]
|
583 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert)
|
583 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert)
|
584 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
584 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
585 | 585 | ||
586 | // step 2: calculate GPP of the tree based
|
586 | // step 2: calculate GPP of the tree based
|
587 | // (1) get the amount of GPP for a "unit area" of the tree species
|
587 | // (1) get the amount of GPP for a "unit area" of the tree species
|
588 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
588 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
589 | // (2) GPP (without aging-effect) in kg Biomass / year
|
589 | // (2) GPP (without aging-effect) in kg Biomass / year
|
590 | double raw_gpp = raw_gpp_per_area * effective_area; |
590 | double raw_gpp = raw_gpp_per_area * effective_area; |
591 | 591 | ||
592 | // apply aging according to the state of the individuum
|
592 | // apply aging according to the state of the individuum
|
593 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
593 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
594 | mRU->addTreeAging(mLeafArea, aging_factor); |
594 | mRU->addTreeAging(mLeafArea, aging_factor); |
595 | double gpp = raw_gpp * aging_factor; // |
595 | double gpp = raw_gpp * aging_factor; // |
596 | d.NPP = gpp * cAutotrophicRespiration; // respiration loss (0.47), cf. Waring et al 1998. |
596 | d.NPP = gpp * cAutotrophicRespiration; // respiration loss (0.47), cf. Waring et al 1998. |
597 | 597 | ||
598 | //DBGMODE(
|
598 | //DBGMODE(
|
599 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
599 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
600 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
600 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
601 | dumpList(out); // add tree headers |
601 | dumpList(out); // add tree headers |
602 | out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
602 | out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
603 | }
|
603 | }
|
604 | //); // DBGMODE()
|
604 | //); // DBGMODE()
|
605 | if (d.NPP>0.) |
605 | if (d.NPP>0.) |
606 | partitioning(d); // split npp to compartments and grow (diameter, height) |
606 | partitioning(d); // split npp to compartments and grow (diameter, height) |
607 | 607 | ||
608 | // mortality
|
608 | // mortality
|
609 | if (Model::settings().mortalityEnabled) |
609 | if (Model::settings().mortalityEnabled) |
610 | mortality(d); |
610 | mortality(d); |
611 | 611 | ||
612 | mStressIndex = d.stress_index; |
612 | mStressIndex = d.stress_index; |
613 | 613 | ||
614 | if (!isDead()) |
614 | if (!isDead()) |
615 | mRU->resourceUnitSpecies(species()).statistics().add(this, &d); |
615 | mRU->resourceUnitSpecies(species()).statistics().add(this, &d); |
616 | 616 | ||
617 | // regeneration
|
617 | // regeneration
|
618 | mSpecies->seedProduction(mAge, mHeight, mPositionIndex); |
618 | mSpecies->seedProduction(mAge, mHeight, mPositionIndex); |
619 | 619 | ||
620 | }
|
620 | }
|
621 | 621 | ||
622 | /** partitioning of this years assimilates (NPP) to biomass compartments.
|
622 | /** partitioning of this years assimilates (NPP) to biomass compartments.
|
623 | Conceptionally, the algorithm is based on Duursma, 2007.
|
623 | Conceptionally, the algorithm is based on Duursma, 2007.
|
624 | @sa http://iland.boku.ac.at/allocation */
|
624 | @sa http://iland.boku.ac.at/allocation */
|
625 | inline void Tree::partitioning(TreeGrowthData &d) |
625 | inline void Tree::partitioning(TreeGrowthData &d) |
626 | {
|
626 | {
|
627 | double npp = d.NPP; |
627 | double npp = d.NPP; |
628 | // add content of reserve pool
|
628 | // add content of reserve pool
|
629 | npp += mNPPReserve; |
629 | npp += mNPPReserve; |
630 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
630 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
631 | const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio()); |
631 | const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio()); |
632 | double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100% |
632 | double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100% |
633 | 633 | ||
634 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
634 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
635 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
635 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
636 | // turnover rates
|
636 | // turnover rates
|
637 | const double &to_fol = species()->turnoverLeaf(); |
637 | const double &to_fol = species()->turnoverLeaf(); |
638 | const double &to_root = species()->turnoverRoot(); |
638 | const double &to_root = species()->turnoverRoot(); |
639 | // the turnover rate of wood depends on the size of the reserve pool:
|
639 | // the turnover rate of wood depends on the size of the reserve pool:
|
640 | 640 | ||
641 | 641 | ||
642 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
642 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
643 | 643 | ||
644 | apct_root = rus.prod3PG().rootFraction(); |
644 | apct_root = rus.prod3PG().rootFraction(); |
645 | d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots |
645 | d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots |
646 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage) |
646 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage) |
647 | 647 | ||
648 | // Duursma 2007, Eq. (20)
|
648 | // Duursma 2007, Eq. (20)
|
649 | apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf ); |
649 | apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf ); |
650 | if (apct_wood<0) |
650 | if (apct_wood<0) |
651 | apct_wood = 0.; |
651 | apct_wood = 0.; |
652 | apct_foliage = 1. - apct_root - apct_wood; |
652 | apct_foliage = 1. - apct_root - apct_wood; |
653 | 653 | ||
654 | 654 | ||
655 | //DBGMODE(
|
655 | //DBGMODE(
|
656 | if (apct_foliage<0 || apct_wood<0) |
656 | if (apct_foliage<0 || apct_wood<0) |
657 | qDebug() << "transfer to foliage or wood < 0"; |
657 | qDebug() << "transfer to foliage or wood < 0"; |
658 | if (npp<0) |
658 | if (npp<0) |
659 | qDebug() << "NPP < 0"; |
659 | qDebug() << "NPP < 0"; |
660 | // );
|
660 | // );
|
661 | 661 | ||
662 | // Change of biomass compartments
|
662 | // Change of biomass compartments
|
663 | double sen_root = mFineRootMass * to_root; |
663 | double sen_root = mFineRootMass * to_root; |
664 | double sen_foliage = mFoliageMass * to_fol; |
664 | double sen_foliage = mFoliageMass * to_fol; |
665 | if (ru()->snag()) |
665 | if (ru()->snag()) |
666 | ru()->snag()->addTurnoverLitter(this->species(), sen_foliage, sen_root); |
666 | ru()->snag()->addTurnoverLitter(this->species(), sen_foliage, sen_root); |
667 | 667 | ||
668 | // Roots
|
668 | // Roots
|
669 | // http://iland.boku.ac.at/allocation#belowground_NPP
|
669 | // http://iland.boku.ac.at/allocation#belowground_NPP
|
670 | mFineRootMass -= sen_root; // reduce only fine root pool |
670 | mFineRootMass -= sen_root; // reduce only fine root pool |
671 | double delta_root = apct_root * npp; |
671 | double delta_root = apct_root * npp; |
672 | // 1st, refill the fine root pool
|
672 | // 1st, refill the fine root pool
|
673 | double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass; |
673 | double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass; |
674 | if (fineroot_miss>0.){ |
674 | if (fineroot_miss>0.){ |
675 | double delta_fineroot = qMin(fineroot_miss, delta_root); |
675 | double delta_fineroot = qMin(fineroot_miss, delta_root); |
676 | mFineRootMass += delta_fineroot; |
676 | mFineRootMass += delta_fineroot; |
677 | delta_root -= delta_fineroot; |
677 | delta_root -= delta_fineroot; |
678 | }
|
678 | }
|
679 | // 2nd, the rest of NPP allocated to roots go to coarse roots
|
679 | // 2nd, the rest of NPP allocated to roots go to coarse roots
|
680 | double max_coarse_root = species()->biomassRoot(mDbh); |
680 | double max_coarse_root = species()->biomassRoot(mDbh); |
681 | mCoarseRootMass += delta_root; |
681 | mCoarseRootMass += delta_root; |
682 | if (mCoarseRootMass > max_coarse_root) { |
682 | if (mCoarseRootMass > max_coarse_root) { |
683 | // if the coarse root pool exceeds the value given by the allometry, then the
|
683 | // if the coarse root pool exceeds the value given by the allometry, then the
|
684 | // surplus is accounted as turnover
|
684 | // surplus is accounted as turnover
|
685 | if (ru()->snag()) |
685 | if (ru()->snag()) |
686 | ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root); |
686 | ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root); |
687 | 687 | ||
688 | mCoarseRootMass = max_coarse_root; |
688 | mCoarseRootMass = max_coarse_root; |
689 | }
|
689 | }
|
690 | 690 | ||
691 | // Foliage
|
691 | // Foliage
|
692 | double delta_foliage = apct_foliage * npp - sen_foliage; |
692 | double delta_foliage = apct_foliage * npp - sen_foliage; |
693 | mFoliageMass += delta_foliage; |
693 | mFoliageMass += delta_foliage; |
694 | if (isnan(mFoliageMass)) |
694 | if (isnan(mFoliageMass)) |
695 | qDebug() << "foliage mass invalid!"; |
695 | qDebug() << "foliage mass invalid!"; |
696 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
696 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
697 | 697 | ||
698 | mLeafArea = mFoliageMass * species()->specificLeafArea(); // update leaf area |
698 | mLeafArea = mFoliageMass * species()->specificLeafArea(); // update leaf area |
699 | 699 | ||
700 | // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild,
|
700 | // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild,
|
701 | // foliage_mass_allo: simply higher chance for stress
|
701 | // foliage_mass_allo: simply higher chance for stress
|
702 | // note: npp = NPP + reserve (see above)
|
702 | // note: npp = NPP + reserve (see above)
|
703 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.); |
703 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.); |
704 | 704 | ||
705 | // Woody compartments
|
705 | // Woody compartments
|
706 | // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth
|
706 | // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth
|
707 | // (1) transfer to reserve pool
|
707 | // (1) transfer to reserve pool
|
708 | double gross_woody = apct_wood * npp; |
708 | double gross_woody = apct_wood * npp; |
709 | double to_reserve = qMin(reserve_size, gross_woody); |
709 | double to_reserve = qMin(reserve_size, gross_woody); |
710 | mNPPReserve = to_reserve; |
710 | mNPPReserve = to_reserve; |
711 | double net_woody = gross_woody - to_reserve; |
711 | double net_woody = gross_woody - to_reserve; |
712 | double net_stem = 0.; |
712 | double net_stem = 0.; |
713 | mDbhDelta = 0.; |
713 | mDbhDelta = 0.; |
714 | 714 | ||
715 | 715 | ||
716 | if (net_woody > 0.) { |
716 | if (net_woody > 0.) { |
717 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter)
|
717 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter)
|
718 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
718 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
719 | d.NPP_stem = net_stem; |
719 | d.NPP_stem = net_stem; |
720 | mWoodyMass += net_woody; |
720 | mWoodyMass += net_woody; |
721 | // (3) growth of diameter and height baseed on net stem increment
|
721 | // (3) growth of diameter and height baseed on net stem increment
|
722 | grow_diameter(d); |
722 | grow_diameter(d); |
723 | }
|
723 | }
|
724 | 724 | ||
725 | //DBGMODE(
|
725 | //DBGMODE(
|
726 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
726 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
727 | && isDebugging() ) { |
727 | && isDebugging() ) { |
728 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
728 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
729 | dumpList(out); // add tree headers |
729 | dumpList(out); // add tree headers |
730 | out << npp << apct_foliage << apct_wood << apct_root |
730 | out << npp << apct_foliage << apct_wood << apct_root |
731 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index; |
731 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index; |
732 | }
|
732 | }
|
733 | 733 | ||
734 | //); // DBGMODE()
|
734 | //); // DBGMODE()
|
735 | DBGMODE(
|
735 | DBGMODE(
|
736 | if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000 |
736 | if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000 |
737 | || mNPPReserve>4000.) { |
737 | || mNPPReserve>4000.) { |
738 | qDebug() << "Tree:partitioning: invalid or unlikely pools."; |
738 | qDebug() << "Tree:partitioning: invalid or unlikely pools."; |
739 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
739 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
740 | DebugList dbg; dumpList(dbg); |
740 | DebugList dbg; dumpList(dbg); |
741 | qDebug() << dbg; |
741 | qDebug() << dbg; |
742 | } ); |
742 | } ); |
743 | 743 | ||
744 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump()
|
744 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump()
|
745 | + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8")
|
745 | + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8")
|
746 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root)
|
746 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root)
|
747 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/
|
747 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/
|
748 | 748 | ||
749 | }
|
749 | }
|
750 | 750 | ||
751 | 751 | ||
752 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp).
|
752 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp).
|
753 | Refer to XXX for equations and variables.
|
753 | Refer to XXX for equations and variables.
|
754 | This function updates the dbh and height of the tree.
|
754 | This function updates the dbh and height of the tree.
|
755 | The equations are based on dbh in meters! */
|
755 | The equations are based on dbh in meters! */
|
756 | inline void Tree::grow_diameter(TreeGrowthData &d) |
756 | inline void Tree::grow_diameter(TreeGrowthData &d) |
757 | {
|
757 | {
|
758 | // determine dh-ratio of increment
|
758 | // determine dh-ratio of increment
|
759 | // height increment is a function of light competition:
|
759 | // height increment is a function of light competition:
|
760 | double hd_growth = relative_height_growth(); // hd of height growth |
760 | double hd_growth = relative_height_growth(); // hd of height growth |
761 | double d_m = mDbh / 100.; // current diameter in [m] |
761 | double d_m = mDbh / 100.; // current diameter in [m] |
762 | double net_stem_npp = d.NPP_stem; |
762 | double net_stem_npp = d.NPP_stem; |
763 | 763 | ||
764 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
764 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
765 | 765 | ||
766 | const double mass_factor = species()->volumeFactor() * species()->density(); |
766 | const double mass_factor = species()->volumeFactor() * species()->density(); |
767 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
767 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
768 | 768 | ||
769 | // factor is in diameter increment per NPP [m/kg]
|
769 | // factor is in diameter increment per NPP [m/kg]
|
770 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
770 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
771 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
771 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
772 | 772 | ||
773 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9)
|
773 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9)
|
774 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
774 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
775 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
775 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
776 | 776 | ||
777 | // the final increment is then:
|
777 | // the final increment is then:
|
778 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
778 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
779 | double res_final = 0.; |
779 | double res_final = 0.; |
780 | if (fabs(stem_residual) > 1.) { |
780 | if (fabs(stem_residual) > 1.) { |
781 | 781 | ||
782 | // calculate final residual in stem
|
782 | // calculate final residual in stem
|
783 | res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
783 | res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
784 | if (fabs(res_final)>1.) { |
784 | if (fabs(res_final)>1.) { |
785 | // for large errors in stem biomass due to errors in diameter increment (> 1kg), we solve the increment iteratively.
|
785 | // for large errors in stem biomass due to errors in diameter increment (> 1kg), we solve the increment iteratively.
|
786 | // first, increase increment with constant step until we overestimate the first time
|
786 | // first, increase increment with constant step until we overestimate the first time
|
787 | // then,
|
787 | // then,
|
788 | d_increment = 0.02; // start with 2cm increment |
788 | d_increment = 0.02; // start with 2cm increment |
789 | bool reached_error = false; |
789 | bool reached_error = false; |
790 | double step=0.01; // step-width 1cm |
790 | double step=0.01; // step-width 1cm |
791 | double est_stem; |
791 | double est_stem; |
792 | do { |
792 | do { |
793 | est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment |
793 | est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment |
794 | stem_residual = est_stem - (stem_mass + net_stem_npp); |
794 | stem_residual = est_stem - (stem_mass + net_stem_npp); |
795 | 795 | ||
796 | if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg |
796 | if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg |
797 | break; |
797 | break; |
798 | if (stem_residual > 0.) { |
798 | if (stem_residual > 0.) { |
799 | d_increment -= step; |
799 | d_increment -= step; |
800 | reached_error=true; |
800 | reached_error=true; |
801 | } else { |
801 | } else { |
802 | d_increment += step; |
802 | d_increment += step; |
803 | }
|
803 | }
|
804 | if (reached_error) |
804 | if (reached_error) |
805 | step /= 2.; |
805 | step /= 2.; |
806 | } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm |
806 | } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm |
807 | }
|
807 | }
|
808 | }
|
808 | }
|
809 | 809 | ||
810 | if (d_increment<0.f) |
810 | if (d_increment<0.f) |
811 | qDebug() << "Tree::grow_diameter: d_inc < 0."; |
811 | qDebug() << "Tree::grow_diameter: d_inc < 0."; |
812 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
812 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
813 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
813 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
814 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
814 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
815 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
815 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
816 | 816 | ||
817 | //DBGMODE(
|
817 | //DBGMODE(
|
818 | // do not calculate res_final twice if already done
|
818 | // do not calculate res_final twice if already done
|
819 | DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
819 | DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
820 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump()); |
820 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump()); |
821 | 821 | ||
822 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
822 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
823 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
823 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
824 | dumpList(out); // add tree headers |
824 | dumpList(out); // add tree headers |
825 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
825 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
826 | }
|
826 | }
|
827 | 827 | ||
828 | //); // DBGMODE()
|
828 | //); // DBGMODE()
|
829 | 829 | ||
830 | d_increment = qMax(d_increment, 0.); |
830 | d_increment = qMax(d_increment, 0.); |
831 | 831 | ||
832 | // update state variables
|
832 | // update state variables
|
833 | mDbh += d_increment*100; // convert from [m] to [cm] |
833 | mDbh += d_increment*100; // convert from [m] to [cm] |
834 | mDbhDelta = d_increment*100; // save for next year's growth |
834 | mDbhDelta = d_increment*100; // save for next year's growth |
835 | mHeight += d_increment * hd_growth; |
835 | mHeight += d_increment * hd_growth; |
836 | 836 | ||
837 | // update state of LIP stamp and opacity
|
837 | // update state of LIP stamp and opacity
|
838 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
838 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
839 | // calculate the CrownFactor which reflects the opacity of the crown
|
839 | // calculate the CrownFactor which reflects the opacity of the crown
|
840 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
840 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
841 | mOpacity = 1. - exp(-k * mLeafArea / mStamp->crownArea()); |
841 | mOpacity = 1. - exp(-k * mLeafArea / mStamp->crownArea()); |
842 | 842 | ||
843 | }
|
843 | }
|
844 | 844 | ||
845 | 845 | ||
846 | /// return the HD ratio of this year's increment based on the light status.
|
846 | /// return the HD ratio of this year's increment based on the light status.
|
847 | inline double Tree::relative_height_growth() |
847 | inline double Tree::relative_height_growth() |
848 | {
|
848 | {
|
849 | double hd_low, hd_high; |
849 | double hd_low, hd_high; |
850 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
850 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
851 | 851 | ||
852 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
852 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
853 | DBG_IF_X(hd_low < 20 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high)); |
853 | DBG_IF_X(hd_low < 20 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high)); |
854 | 854 | ||
855 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees)
|
855 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees)
|
856 | // use the corrected LRI (see tracker#11)
|
856 | // use the corrected LRI (see tracker#11)
|
857 | double lri = limit(mLRI * mRU->LRImodifier(),0.,1.); |
857 | double lri = limit(mLRI * mRU->LRImodifier(),0.,1.); |
858 | double hd_ratio = hd_high - (hd_high-hd_low)*lri; |
858 | double hd_ratio = hd_high - (hd_high-hd_low)*lri; |
859 | return hd_ratio; |
859 | return hd_ratio; |
860 | }
|
860 | }
|
861 | 861 | ||
862 | /** This function is called if a tree dies.
|
862 | /** This function is called if a tree dies.
|
863 | @sa ResourceUnit::cleanTreeList(), remove() */
|
863 | @sa ResourceUnit::cleanTreeList(), remove() */
|
864 | void Tree::die(TreeGrowthData *d) |
864 | void Tree::die(TreeGrowthData *d) |
865 | {
|
865 | {
|
866 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
866 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
- | 867 | mRU->treeDied(); |
|
867 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
868 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
868 | rus.statisticsDead().add(this, d); // add tree to statistics |
869 | rus.statisticsDead().add(this, d); // add tree to statistics |
869 | if (ru()->snag()) |
870 | if (ru()->snag()) |
870 | ru()->snag()->addMortality(this); |
871 | ru()->snag()->addMortality(this); |
871 | }
|
872 | }
|
872 | 873 | ||
873 | void Tree::remove(double removeFoliage, double removeBranch, double removeStem ) |
874 | void Tree::remove(double removeFoliage, double removeBranch, double removeStem ) |
874 | {
|
875 | {
|
875 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
876 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
- | 877 | mRU->treeDied(); |
|
876 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
878 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
877 | rus.statisticsMgmt().add(this, 0); |
879 | rus.statisticsMgmt().add(this, 0); |
878 | if (ru()->snag()) |
880 | if (ru()->snag()) |
879 | ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage); |
881 | ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage); |
880 | }
|
882 | }
|
881 | 883 | ||
882 | void Tree::mortality(TreeGrowthData &d) |
884 | void Tree::mortality(TreeGrowthData &d) |
883 | {
|
885 | {
|
884 | // death if leaf area is 0
|
886 | // death if leaf area is 0
|
885 | if (mFoliageMass<0.00001) |
887 | if (mFoliageMass<0.00001) |
886 | die(); |
888 | die(); |
887 | 889 | ||
888 | double p_death, p_stress, p_intrinsic; |
890 | double p_death, p_stress, p_intrinsic; |
889 | p_intrinsic = species()->deathProb_intrinsic(); |
891 | p_intrinsic = species()->deathProb_intrinsic(); |
890 | p_stress = species()->deathProb_stress(d.stress_index); |
892 | p_stress = species()->deathProb_stress(d.stress_index); |
891 | p_death = p_intrinsic + p_stress; |
893 | p_death = p_intrinsic + p_stress; |
892 | double p = drandom(ru()->randomGenerator()); //0..1 |
894 | double p = drandom(ru()->randomGenerator()); //0..1 |
893 | if (p<p_death) { |
895 | if (p<p_death) { |
894 | // die...
|
896 | // die...
|
895 | die(); |
897 | die(); |
896 | }
|
898 | }
|
897 | }
|
899 | }
|
898 | 900 | ||
899 | //////////////////////////////////////////////////
|
901 | //////////////////////////////////////////////////
|
900 | //// value functions
|
902 | //// value functions
|
901 | //////////////////////////////////////////////////
|
903 | //////////////////////////////////////////////////
|
902 | 904 | ||
903 | double Tree::volume() const |
905 | double Tree::volume() const |
904 | {
|
906 | {
|
905 | /// @see Species::volumeFactor() for details
|
907 | /// @see Species::volumeFactor() for details
|
906 | const double volume_factor = species()->volumeFactor(); |
908 | const double volume_factor = species()->volumeFactor(); |
907 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
909 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
908 | return volume; |
910 | return volume; |
909 | }
|
911 | }
|
910 | 912 | ||
911 | /// return the basal area in m2
|
913 | /// return the basal area in m2
|
912 | double Tree::basalArea() const |
914 | double Tree::basalArea() const |
913 | {
|
915 | {
|
914 | // A = r^2 * pi = d/2*pi; from cm->m: d/200
|
916 | // A = r^2 * pi = d/2*pi; from cm->m: d/200
|
915 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
917 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
916 | return b; |
918 | return b; |
917 | }
|
919 | }
|
918 | 920 |