Annotation of /trunk/gadget/growthcalc.cc

Parent Directory | Revision Log

---

Revision 4 - (view) (download)

1 :	agomez	1	#include "growthcalc.h"
2 :			#include "errorhandler.h"
3 :			#include "areatime.h"
4 :			#include "keeper.h"
5 :			#include "readfunc.h"
6 :			#include "readword.h"
7 :			#include "gadget.h"
8 :			#include "global.h"
9 :
10 :			// ********************************************************
11 :			// Functions for GrowthCalcBase
12 :			// ********************************************************
13 :			GrowthCalcBase::GrowthCalcBase(const IntVector& Areas) : LivesOnAreas(Areas) {
14 :			}
15 :
16 :			// ********************************************************
17 :			// Functions for GrowthCalcA
18 :			// ********************************************************
19 :			GrowthCalcA::GrowthCalcA(CommentStream& infile, const IntVector& Areas,
20 :			const TimeClass* const TimeInfo, Keeper* const keeper)
21 :			: GrowthCalcBase(Areas), numGrowthConstants(9) {
22 :
23 :			keeper->addString("growthcalcA");
24 :			growthPar.resize(numGrowthConstants, keeper);
25 :
26 :			char text[MaxStrLength];
27 :			strncpy(text, "", MaxStrLength);
28 :			infile >> text >> ws;
29 :			if (strcasecmp(text, "growthparameters") != 0)
30 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
31 :			growthPar.read(infile, TimeInfo, keeper);
32 :			keeper->clearLast();
33 :			}
34 :
35 :	ulcessvp	4	void GrowthCalcA::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
36 :			const PopInfoVector& numGrow, const AreaClass* const Area,
37 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
38 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
39 :	agomez	1
40 :			growthPar.Update(TimeInfo);
41 :			double tempL = TimeInfo->getTimeStepSize() * growthPar[0] *
42 :			(growthPar[2] * Area->getTemperature(area, TimeInfo->getTime()) + growthPar[3]);
43 :			double tempW = TimeInfo->getTimeStepSize() * growthPar[4] *
44 :			(growthPar[7] * Area->getTemperature(area, TimeInfo->getTime()) + growthPar[8]);
45 :
46 :			int i;
47 :	ulcessvp	4	for (i = 0; i < size; i++) {
48 :	agomez	1	Lgrowth[i] = tempL * pow(LgrpDiv->meanLength(i), growthPar[1]) * Fphi[i];
49 :			if (Lgrowth[i] < 0.0)
50 :			Lgrowth[i] = 0.0;
51 :
52 :			if (numGrow[i].W < verysmall || isZero(tempW))
53 :			Wgrowth[i] = 0.0;
54 :			else
55 :			Wgrowth[i] = tempW * pow(numGrow[i].W, growthPar[5]) * (Fphi[i] - growthPar[6]);
56 :			}
57 :			}
58 :
59 :			// ********************************************************
60 :			// Functions for GrowthCalcB
61 :			// ********************************************************
62 :			GrowthCalcB::GrowthCalcB(CommentStream& infile, const IntVector& Areas,
63 :			const TimeClass* const TimeInfo, Keeper* const keeper,
64 :			const AreaClass* const Area, const CharPtrVector& lenindex)
65 :			: GrowthCalcBase(Areas) {
66 :
67 :			char datafilename[MaxStrLength];
68 :			strncpy(datafilename, "", MaxStrLength);
69 :			ifstream datafile;
70 :			CommentStream subdata(datafile);
71 :
72 :			int i;
73 :			for (i = 0; i < Areas.Size(); i++) {
74 :			lgrowth.resize(new FormulaMatrix(TimeInfo->numTotalSteps() + 1, lenindex.Size(), 0.0));
75 :			wgrowth.resize(new FormulaMatrix(TimeInfo->numTotalSteps() + 1, lenindex.Size(), 0.0));
76 :			}
77 :
78 :			keeper->addString("growthcalcB");
79 :
80 :			readWordAndValue(infile, "lengthgrowthfile", datafilename);
81 :			datafile.open(datafilename, ios::in);
82 :			handle.checkIfFailure(datafile, datafilename);
83 :			handle.Open(datafilename);
84 :			readGrowthAmounts(subdata, TimeInfo, Area, lgrowth, lenindex, Areas);
85 :			handle.Close();
86 :			datafile.close();
87 :			datafile.clear();
88 :
89 :			readWordAndValue(infile, "weightgrowthfile", datafilename);
90 :			datafile.open(datafilename, ios::in);
91 :			handle.checkIfFailure(datafile, datafilename);
92 :			handle.Open(datafilename);
93 :			readGrowthAmounts(subdata, TimeInfo, Area, wgrowth, lenindex, Areas);
94 :			handle.Close();
95 :			datafile.close();
96 :			datafile.clear();
97 :
98 :			for (i = 0; i < Areas.Size(); i++) {
99 :			(*lgrowth[i]).Inform(keeper);
100 :			(*wgrowth[i]).Inform(keeper);
101 :			}
102 :
103 :			keeper->clearLast();
104 :			}
105 :
106 :			GrowthCalcB::~GrowthCalcB() {
107 :			int a;
108 :			for (a = 0; a < lgrowth.Size(); a++) {
109 :			delete lgrowth[a];
110 :			delete wgrowth[a];
111 :			}
112 :			}
113 :
114 :	ulcessvp	4	void GrowthCalcB::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
115 :			const PopInfoVector& numGrow, const AreaClass* const Area,
116 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
117 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
118 :	agomez	1
119 :			int i, inarea = this->areaNum(area);
120 :	ulcessvp	4	for (i = 0; i < size; i++) {
121 :	agomez	1	Lgrowth[i] = (*lgrowth[inarea])[TimeInfo->getTime()][i];
122 :			Wgrowth[i] = (*wgrowth[inarea])[TimeInfo->getTime()][i];
123 :			if ((handle.getLogLevel() >= LOGWARN) && ((Lgrowth[i] < 0.0) || (Wgrowth[i] < 0.0)))
124 :			handle.logMessage(LOGWARN, "Warning in growth calculation - negative growth parameter");
125 :			}
126 :			}
127 :
128 :			// ********************************************************
129 :			// Functions for GrowthCalcC
130 :			// ********************************************************
131 :			GrowthCalcC::GrowthCalcC(CommentStream& infile, const IntVector& Areas,
132 :			const TimeClass* const TimeInfo, const LengthGroupDivision* const LgrpDiv,
133 :			Keeper* const keeper, const char* refWeightFile)
134 :			: GrowthCalcBase(Areas), numWeightGrowthConstants(6), numLengthGrowthConstants(9) {
135 :
136 :			keeper->addString("growthcalcC");
137 :			wgrowthPar.resize(numWeightGrowthConstants, keeper);
138 :			lgrowthPar.resize(numLengthGrowthConstants, keeper);
139 :			refWeight.resize(LgrpDiv->numLengthGroups(), 0.0);
140 :
141 :			int i, j;
142 :			char text[MaxStrLength];
143 :			strncpy(text, "", MaxStrLength);
144 :			infile >> text >> ws;
145 :			if (strcasecmp(text, "wgrowthparameters") != 0)
146 :			handle.logFileUnexpected(LOGFAIL, "wgrowthparameters", text);
147 :			wgrowthPar.read(infile, TimeInfo, keeper);
148 :
149 :			infile >> text >> ws;
150 :			if (strcasecmp(text, "lgrowthparameters") != 0)
151 :			handle.logFileUnexpected(LOGFAIL, "lgrowthparameters", text);
152 :			lgrowthPar.read(infile, TimeInfo, keeper);
153 :
154 :			//read information on reference weights.
155 :			keeper->addString("referenceweights");
156 :			ifstream subfile;
157 :			subfile.open(refWeightFile, ios::in);
158 :			handle.checkIfFailure(subfile, refWeightFile);
159 :			handle.Open(refWeightFile);
160 :			CommentStream subcomment(subfile);
161 :			DoubleMatrix tmpRefW;
162 :			readRefWeights(subcomment, tmpRefW);
163 :			handle.Close();
164 :			subfile.close();
165 :			subfile.clear();
166 :
167 :			//Interpolate the reference weights. First there are some error checks.
168 :			if (LgrpDiv->meanLength(0) < tmpRefW[0][0] ||
169 :			LgrpDiv->meanLength(LgrpDiv->numLengthGroups() - 1) > tmpRefW[tmpRefW.Nrow() - 1][0])
170 :			handle.logFileMessage(LOGFAIL, "lengths for reference weights must span the range of growth lengths");
171 :
172 :			double ratio, tmplen;
173 :			int pos = 0;
174 :			for (j = 0; j < LgrpDiv->numLengthGroups(); j++) {
175 :			tmplen = LgrpDiv->meanLength(j);
176 :			for (i = pos; i < tmpRefW.Nrow() - 1; i++) {
177 :			if (((tmplen > tmpRefW[i][0]) || (isEqual(tmplen, tmpRefW[i][0]))) &&
178 :			((tmplen < tmpRefW[i + 1][0]) || (isEqual(tmplen, tmpRefW[i + 1][0])))) {
179 :
180 :			ratio = (tmplen - tmpRefW[i][0]) / (tmpRefW[i + 1][0] - tmpRefW[i][0]);
181 :			refWeight[j] = tmpRefW[i][1] + ratio * (tmpRefW[i + 1][1] - tmpRefW[i][1]);
182 :			pos = i;
183 :			}
184 :			}
185 :			}
186 :			keeper->clearLast();
187 :			keeper->clearLast();
188 :			}
189 :
190 :			/* Von Bertalanffy growth function. dw/dt = a*w^n - b*w^m;
191 :			* As a generalisation a and b are made temperature dependent so the
192 :			* final form of the function is
193 :			* dw/dt = a0*exp(a1*T)*((w/a2)^a4 - (w/a3)^a5)
194 :			* For no temperature dependency a1 = 0 */
195 :	ulcessvp	4	void GrowthCalcC::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
196 :			const PopInfoVector& numGrow, const AreaClass* const Area,
197 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
198 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
199 :	agomez	1
200 :			wgrowthPar.Update(TimeInfo);
201 :			lgrowthPar.Update(TimeInfo);
202 :
203 :			//JMB - first some error checking
204 :			if (handle.getLogLevel() >= LOGWARN) {
205 :			if (isZero(wgrowthPar[2]) || isZero(wgrowthPar[3]))
206 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is zero");
207 :			if (isZero(lgrowthPar[6]) || isZero(lgrowthPar[7]))
208 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is zero");
209 :			if (lgrowthPar[5] < 0.0)
210 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is negative");
211 :			}
212 :
213 :			int i;
214 :			double fx;
215 :			double ratio = lgrowthPar[0] + lgrowthPar[8] * (lgrowthPar[1] + lgrowthPar[2] * lgrowthPar[8]);
216 :			double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[0] *
217 :			exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));
218 :
219 :	ulcessvp	4	for (i = 0; i < size; i++) {
220 :	agomez	1	if (numGrow[i].W < verysmall || isZero(tempW)) {
221 :			Wgrowth[i] = 0.0;
222 :			Lgrowth[i] = 0.0;
223 :			} else {
224 :			Wgrowth[i] = tempW * (pow(numGrow[i].W / wgrowthPar[2], wgrowthPar[4]) -
225 :			pow(numGrow[i].W / wgrowthPar[3], wgrowthPar[5]));
226 :
227 :			if (Wgrowth[i] < verysmall) {
228 :			Wgrowth[i] = 0.0;
229 :			Lgrowth[i] = 0.0;
230 :			} else {
231 :			fx = lgrowthPar[3] + (lgrowthPar[4] * (numGrow[i].W - ratio * refWeight[i]) / numGrow[i].W);
232 :			if (fx > lgrowthPar[5])
233 :			fx = lgrowthPar[5];
234 :			if (fx < verysmall)
235 :			Lgrowth[i] = 0.0;
236 :			else
237 :			Lgrowth[i] = fx * Wgrowth[i] / (lgrowthPar[6] * lgrowthPar[7] *
238 :			pow(LgrpDiv->meanLength(i), lgrowthPar[7] - 1.0));
239 :			}
240 :			}
241 :			}
242 :			}
243 :
244 :			// ********************************************************
245 :			// Functions for GrowthCalcD
246 :			// ********************************************************
247 :			GrowthCalcD::GrowthCalcD(CommentStream& infile, const IntVector& Areas,
248 :			const TimeClass* const TimeInfo, const LengthGroupDivision* const LgrpDiv,
249 :			Keeper* const keeper, const char* refWeightFile)
250 :			: GrowthCalcBase(Areas), numWeightGrowthConstants(6), numLengthGrowthConstants(8) {
251 :
252 :			keeper->addString("growthcalcD");
253 :			wgrowthPar.resize(numWeightGrowthConstants, keeper);
254 :			lgrowthPar.resize(numLengthGrowthConstants, keeper);
255 :			refWeight.resize(LgrpDiv->numLengthGroups(), 0.0);
256 :
257 :			int i, j;
258 :			char text[MaxStrLength];
259 :			strncpy(text, "", MaxStrLength);
260 :			infile >> text >> ws;
261 :			if (strcasecmp(text, "wgrowthparameters") != 0)
262 :			handle.logFileUnexpected(LOGFAIL, "wgrowthparameters", text);
263 :			wgrowthPar.read(infile, TimeInfo, keeper);
264 :
265 :			infile >> text >> ws;
266 :			if (strcasecmp(text, "lgrowthparameters") != 0)
267 :			handle.logFileUnexpected(LOGFAIL, "lgrowthparameters", text);
268 :			lgrowthPar.read(infile, TimeInfo, keeper);
269 :
270 :			//read information on reference weights.
271 :			keeper->addString("referenceweights");
272 :			ifstream subfile;
273 :			subfile.open(refWeightFile, ios::in);
274 :			handle.checkIfFailure(subfile, refWeightFile);
275 :			handle.Open(refWeightFile);
276 :			CommentStream subcomment(subfile);
277 :			DoubleMatrix tmpRefW;
278 :			readRefWeights(subcomment, tmpRefW);
279 :			handle.Close();
280 :			subfile.close();
281 :			subfile.clear();
282 :
283 :			//Interpolate the reference weights. First there are some error checks.
284 :			if (LgrpDiv->meanLength(0) < tmpRefW[0][0] ||
285 :			LgrpDiv->meanLength(LgrpDiv->numLengthGroups() - 1) > tmpRefW[tmpRefW.Nrow() - 1][0])
286 :			handle.logFileMessage(LOGFAIL, "lengths for reference weights must span the range of growth lengths");
287 :
288 :			double ratio, tmplen;
289 :			int pos = 0;
290 :			for (j = 0; j < LgrpDiv->numLengthGroups(); j++) {
291 :			tmplen = LgrpDiv->meanLength(j);
292 :			for (i = pos; i < tmpRefW.Nrow() - 1; i++) {
293 :			if (((tmplen > tmpRefW[i][0]) || (isEqual(tmplen, tmpRefW[i][0]))) &&
294 :			((tmplen < tmpRefW[i + 1][0]) || (isEqual(tmplen, tmpRefW[i + 1][0])))) {
295 :
296 :			ratio = (tmplen - tmpRefW[i][0]) / (tmpRefW[i + 1][0] - tmpRefW[i][0]);
297 :			refWeight[j] = tmpRefW[i][1] + ratio * (tmpRefW[i + 1][1] - tmpRefW[i][1]);
298 :			pos = i;
299 :			}
300 :			}
301 :			}
302 :			keeper->clearLast();
303 :			keeper->clearLast();
304 :			}
305 :
306 :			/* Growth function from Jones 1978. Found from experiment in captivity.
307 :			* Jones formula only applies to weight increase. The length increase
308 :			* part is derived from the weight increase part by assuming a formula
309 :			* w = a*l^b. If the weight is below the curve no length increase takes place
310 :			* but instead the weight increases until it reaches the curve. */
311 :	ulcessvp	4	void GrowthCalcD::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
312 :			const PopInfoVector& numGrow, const AreaClass* const Area,
313 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
314 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
315 :	agomez	1
316 :			wgrowthPar.Update(TimeInfo);
317 :			lgrowthPar.Update(TimeInfo);
318 :
319 :			//JMB - first some error checking
320 :			if (handle.getLogLevel() >= LOGWARN) {
321 :			if (isZero(wgrowthPar[0]))
322 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is zero");
323 :			if (isZero(lgrowthPar[6]) || isZero(lgrowthPar[7]))
324 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is zero");
325 :			if (lgrowthPar[5] < 0.0)
326 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is negative");
327 :			}
328 :
329 :			int i;
330 :			double ratio, fx;
331 :			double tempC = TimeInfo->getTimeStepSize() / wgrowthPar[0];
332 :			double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[2] *
333 :			exp(wgrowthPar[4] * Area->getTemperature(area, TimeInfo->getTime()) + wgrowthPar[5]);
334 :
335 :	ulcessvp	4	for (i = 0; i < size; i++) {
336 :	agomez	1	if (numGrow[i].W < verysmall) {
337 :			Wgrowth[i] = 0.0;
338 :			Lgrowth[i] = 0.0;
339 :			} else {
340 :			Wgrowth[i] = Fphi[i] * MaxCon[i] * tempC / (pow(numGrow[i].W, wgrowthPar[1])) -
341 :			tempW * pow(numGrow[i].W, wgrowthPar[3]);
342 :
343 :			if (Wgrowth[i] < verysmall) {
344 :			Wgrowth[i] = 0.0;
345 :			Lgrowth[i] = 0.0;
346 :			} else {
347 :			ratio = lgrowthPar[0] + Fphi[i] * (lgrowthPar[1] + lgrowthPar[2] * Fphi[i]);
348 :			fx = lgrowthPar[3] + (lgrowthPar[4] * (numGrow[i].W - ratio * refWeight[i]) / numGrow[i].W);
349 :			if (fx > lgrowthPar[5])
350 :			fx = lgrowthPar[5];
351 :			if (fx < verysmall)
352 :			Lgrowth[i] = 0.0;
353 :			else
354 :			Lgrowth[i] = fx * Wgrowth[i] / (lgrowthPar[6] * lgrowthPar[7] *
355 :			pow(LgrpDiv->meanLength(i), lgrowthPar[7] - 1.0));
356 :			}
357 :			}
358 :			}
359 :			}
360 :
361 :			// ********************************************************
362 :			// Functions for GrowthCalcE
363 :			// ********************************************************
364 :			GrowthCalcE::GrowthCalcE(CommentStream& infile, const IntVector& Areas,
365 :			const TimeClass* const TimeInfo, const LengthGroupDivision* const LgrpDiv,
366 :			Keeper* const keeper, const char* refWeightFile)
367 :			: GrowthCalcBase(Areas), numWeightGrowthConstants(6), numLengthGrowthConstants(9) {
368 :
369 :			keeper->addString("growthcalcE");
370 :			wgrowthPar.resize(numWeightGrowthConstants, keeper);
371 :			lgrowthPar.resize(numLengthGrowthConstants, keeper);
372 :			refWeight.resize(LgrpDiv->numLengthGroups(), 0.0);
373 :			yearEffect.resize(TimeInfo->getLastYear() - TimeInfo->getFirstYear() + 1, keeper);
374 :			stepEffect.resize(TimeInfo->numSteps(), keeper);
375 :			areaEffect.resize(Areas.Size(), keeper);
376 :
377 :			int i, j;
378 :			char text[MaxStrLength];
379 :			strncpy(text, "", MaxStrLength);
380 :			infile >> text >> ws;
381 :			if (strcasecmp(text, "wgrowthparameters") != 0)
382 :			handle.logFileUnexpected(LOGFAIL, "wgrowthparameters", text);
383 :			wgrowthPar.read(infile, TimeInfo, keeper);
384 :
385 :			infile >> text >> ws;
386 :			if (strcasecmp(text, "lgrowthparameters") != 0)
387 :			handle.logFileUnexpected(LOGFAIL, "lgrowthparameters", text);
388 :			lgrowthPar.read(infile, TimeInfo, keeper);
389 :
390 :			infile >> text >> ws;
391 :			if (strcasecmp(text, "yeareffect") != 0)
392 :			handle.logFileUnexpected(LOGFAIL, "yeareffect", text);
393 :			for (i = 0; i < yearEffect.Size(); i++)
394 :			if (!(infile >> yearEffect[i]))
395 :			handle.logFileMessage(LOGFAIL, "invalid format of yeareffect vector");
396 :			yearEffect.Inform(keeper);
397 :
398 :			infile >> text >> ws;
399 :			if (strcasecmp(text, "stepeffect") != 0)
400 :			handle.logFileUnexpected(LOGFAIL, "stepeffect", text);
401 :			for (i = 0; i < stepEffect.Size(); i++)
402 :			if (!(infile >> stepEffect[i]))
403 :			handle.logFileMessage(LOGFAIL, "invalid format of stepeffect vector");
404 :			stepEffect.Inform(keeper);
405 :
406 :			infile >> text >> ws;
407 :			if (strcasecmp(text, "areaeffect") != 0)
408 :			handle.logFileUnexpected(LOGFAIL, "areaeffect", text);
409 :			for (i = 0; i < areaEffect.Size(); i++)
410 :			if (!(infile >> areaEffect[i]))
411 :			handle.logFileMessage(LOGFAIL, "invalid format of areaeffect vector");
412 :			areaEffect.Inform(keeper);
413 :
414 :			//read information on reference weights.
415 :			keeper->addString("referenceweights");
416 :			ifstream subfile;
417 :			subfile.open(refWeightFile, ios::in);
418 :			handle.checkIfFailure(subfile, refWeightFile);
419 :			handle.Open(refWeightFile);
420 :			CommentStream subcomment(subfile);
421 :			DoubleMatrix tmpRefW;
422 :			readRefWeights(subcomment, tmpRefW);
423 :			handle.Close();
424 :			subfile.close();
425 :			subfile.clear();
426 :
427 :			//Interpolate the reference weights.
428 :			if (LgrpDiv->meanLength(0) < tmpRefW[0][0] ||
429 :			LgrpDiv->meanLength(LgrpDiv->numLengthGroups() - 1) > tmpRefW[tmpRefW.Nrow() - 1][0])
430 :			handle.logFileMessage(LOGFAIL, "lengths for reference weights must span the range of growth lengths");
431 :
432 :			double ratio, tmplen;
433 :			int pos = 0;
434 :			for (j = 0; j < LgrpDiv->numLengthGroups(); j++) {
435 :			tmplen = LgrpDiv->meanLength(j);
436 :			for (i = pos; i < tmpRefW.Nrow() - 1; i++) {
437 :			if (((tmplen > tmpRefW[i][0]) || (isEqual(tmplen, tmpRefW[i][0]))) &&
438 :			((tmplen < tmpRefW[i + 1][0]) || (isEqual(tmplen, tmpRefW[i + 1][0])))) {
439 :
440 :			ratio = (tmplen - tmpRefW[i][0]) / (tmpRefW[i + 1][0] - tmpRefW[i][0]);
441 :			refWeight[j] = tmpRefW[i][1] + ratio * (tmpRefW[i + 1][1] - tmpRefW[i][1]);
442 :			pos = i;
443 :			}
444 :			}
445 :			}
446 :			keeper->clearLast();
447 :			keeper->clearLast();
448 :			}
449 :
450 :			/* Growthfunction to be tested for capelin.
451 :			* The weight growth here is given by the formula
452 :			* dw/dt = a0*factor(year)*factor(area)*factor(Step)*w^a1 - a2*w^a3
453 :			* This is Von Bertanlanffy equation with possibility to make the growth
454 :			* year, area and Step dependent. 3 vectors were added to the class
455 :			* YearEffect
456 :			* AreaEffect
457 :			* StepEffect
458 :			* Length increase is upgraded in the same way as earlier. */
459 :	ulcessvp	4	void GrowthCalcE::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
460 :			const PopInfoVector& numGrow, const AreaClass* const Area,
461 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
462 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
463 :	agomez	1
464 :			wgrowthPar.Update(TimeInfo);
465 :			lgrowthPar.Update(TimeInfo);
466 :
467 :			//JMB - first some error checking
468 :			if (handle.getLogLevel() >= LOGWARN) {
469 :			if (isZero(wgrowthPar[2]) || isZero(wgrowthPar[3]))
470 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is zero");
471 :			if (isZero(lgrowthPar[6]) || isZero(lgrowthPar[7]))
472 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is zero");
473 :			if (lgrowthPar[5] < 0.0)
474 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is negative");
475 :			}
476 :
477 :			int i;
478 :			double fx;
479 :			double factor = yearEffect[TimeInfo->getYear() - TimeInfo->getFirstYear()] *
480 :			stepEffect[TimeInfo->getStep() - 1] * areaEffect[this->areaNum(area)];
481 :			double ratio = lgrowthPar[0] + lgrowthPar[8] * (lgrowthPar[1] + lgrowthPar[2] * lgrowthPar[8]);
482 :			double tempW = factor * TimeInfo->getTimeStepSize() * wgrowthPar[0] *
483 :			exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));
484 :
485 :	ulcessvp	4	for (i = 0; i < size; i++) {
486 :	agomez	1	if (numGrow[i].W < verysmall || isZero(tempW)) {
487 :			Wgrowth[i] = 0.0;
488 :			Lgrowth[i] = 0.0;
489 :			} else {
490 :			Wgrowth[i] = tempW * (pow(numGrow[i].W / wgrowthPar[2], wgrowthPar[4]) -
491 :			pow(numGrow[i].W / wgrowthPar[3], wgrowthPar[5]));
492 :
493 :			if (Wgrowth[i] < verysmall) {
494 :			Wgrowth[i] = 0.0;
495 :			Lgrowth[i] = 0.0;
496 :			} else {
497 :			fx = lgrowthPar[3] + (lgrowthPar[4] * (numGrow[i].W - ratio * refWeight[i]) / numGrow[i].W);
498 :			if (fx > lgrowthPar[5])
499 :			fx = lgrowthPar[5];
500 :			if (fx < verysmall)
501 :			Lgrowth[i] = 0.0;
502 :			else
503 :			Lgrowth[i] = fx * Wgrowth[i] / (lgrowthPar[6] * lgrowthPar[7] *
504 :			pow(LgrpDiv->meanLength(i), lgrowthPar[7] - 1.0));
505 :			}
506 :			}
507 :			}
508 :			}
509 :
510 :			// ********************************************************
511 :			// Functions for GrowthCalcF
512 :			// ********************************************************
513 :			/* The Norwegian Growthfunction Equations 6 and 7 on page 7 in
514 :			* Description of a multispecies model for the Barents Sea.
515 :			* parameter # 0 corresponds to C4 in the equation etc. */
516 :			GrowthCalcF::GrowthCalcF(CommentStream& infile, const IntVector& Areas,
517 :			const TimeClass* const TimeInfo, Keeper* const keeper,
518 :			const AreaClass* const Area, const CharPtrVector& lenindex)
519 :			: GrowthCalcBase(Areas), numGrowthConstants(2) {
520 :
521 :			keeper->addString("growthcalcF");
522 :			growthPar.resize(numGrowthConstants, keeper);
523 :
524 :			int i;
525 :			char text[MaxStrLength];
526 :			strncpy(text, "", MaxStrLength);
527 :			infile >> text >> ws;
528 :			if (strcasecmp(text, "growthparameters") != 0)
529 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
530 :			growthPar.read(infile, TimeInfo, keeper);
531 :
532 :			for (i = 0; i < Areas.Size(); i++)
533 :			wgrowth.resize(new FormulaMatrix(TimeInfo->numTotalSteps() + 1, lenindex.Size(), 0.0));
534 :
535 :			ifstream datafile;
536 :			CommentStream subdata(datafile);
537 :			readWordAndValue(infile, "weightgrowthfile", text);
538 :			datafile.open(text, ios::in);
539 :			handle.checkIfFailure(datafile, text);
540 :			handle.Open(text);
541 :			readGrowthAmounts(subdata, TimeInfo, Area, wgrowth, lenindex, Areas);
542 :			handle.Close();
543 :			datafile.close();
544 :			datafile.clear();
545 :
546 :			for (i = 0; i < Areas.Size(); i++)
547 :			(*wgrowth[i]).Inform(keeper);
548 :
549 :			keeper->clearLast();
550 :			}
551 :
552 :			GrowthCalcF::~GrowthCalcF() {
553 :			int a;
554 :			for (a = 0; a < wgrowth.Size(); a++)
555 :			delete wgrowth[a];
556 :			}
557 :
558 :	ulcessvp	4	void GrowthCalcF::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
559 :			const PopInfoVector& numGrow, const AreaClass* const Area,
560 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
561 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
562 :	agomez	1
563 :			growthPar.Update(TimeInfo);
564 :			int i, t, inarea;
565 :			t = TimeInfo->getTime();
566 :			inarea = this->areaNum(area);
567 :			double kval = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
568 :
569 :	ulcessvp	4	for (i = 0; i < size; i++) {
570 :	agomez	1	Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * kval;
571 :			Wgrowth[i] = (*wgrowth[inarea])[t][i];
572 :			if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
573 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
574 :			}
575 :			}
576 :
577 :			// ********************************************************
578 :			// Functions for GrowthCalcG
579 :			// ********************************************************
580 :			GrowthCalcG::GrowthCalcG(CommentStream& infile, const IntVector& Areas,
581 :			const TimeClass* const TimeInfo, Keeper* const keeper,
582 :			const AreaClass* const Area, const CharPtrVector& lenindex)
583 :			: GrowthCalcBase(Areas), numGrowthConstants(2) {
584 :
585 :			keeper->addString("growthcalcG");
586 :			growthPar.resize(numGrowthConstants, keeper);
587 :
588 :			int i;
589 :			char text[MaxStrLength];
590 :			strncpy(text, "", MaxStrLength);
591 :			infile >> text >> ws;
592 :			if (strcasecmp(text, "growthparameters") != 0)
593 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
594 :			growthPar.read(infile, TimeInfo, keeper);
595 :
596 :			for (i = 0; i < Areas.Size(); i++)
597 :			wgrowth.resize(new FormulaMatrix(TimeInfo->numTotalSteps() + 1, lenindex.Size(), 0.0));
598 :
599 :			ifstream datafile;
600 :			CommentStream subdata(datafile);
601 :			readWordAndValue(infile, "weightgrowthfile", text);
602 :			datafile.open(text, ios::in);
603 :			handle.checkIfFailure(datafile, text);
604 :			handle.Open(text);
605 :			readGrowthAmounts(subdata, TimeInfo, Area, wgrowth, lenindex, Areas);
606 :			handle.Close();
607 :			datafile.close();
608 :			datafile.clear();
609 :
610 :			for (i = 0; i < Areas.Size(); i++)
611 :			(*wgrowth[i]).Inform(keeper);
612 :
613 :			keeper->clearLast();
614 :			}
615 :
616 :			GrowthCalcG::~GrowthCalcG() {
617 :			int a;
618 :			for (a = 0; a < wgrowth.Size(); a++)
619 :			delete wgrowth[a];
620 :			}
621 :
622 :	ulcessvp	4	void GrowthCalcG::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
623 :			const PopInfoVector& numGrow, const AreaClass* const Area,
624 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
625 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
626 :	agomez	1
627 :			//written by kgf 24/10 00
628 :			//Gives linear growth (growthPar[0] == 0) or
629 :			//growth decreasing with length (growthPar[0] < 0)
630 :			growthPar.Update(TimeInfo);
631 :			int i, t, inarea;
632 :			t = TimeInfo->getTime();
633 :			inarea = this->areaNum(area);
634 :			double kval = growthPar[1] * TimeInfo->getTimeStepSize();
635 :
636 :			if ((handle.getLogLevel() >= LOGWARN) && (growthPar[0] > 0.0))
637 :			handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is positive");
638 :
639 :			if (isZero(growthPar[0])) {
640 :	ulcessvp	4	for (i = 0; i < size; i++) {
641 :	agomez	1	Lgrowth[i] = kval;
642 :			Wgrowth[i] = (*wgrowth[inarea])[t][i];
643 :			if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
644 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
645 :			}
646 :			} else if (isEqual(growthPar[0], 1.0)) {
647 :	ulcessvp	4	for (i = 0; i < size; i++) {
648 :	agomez	1	Lgrowth[i] = kval * LgrpDiv->meanLength(i);
649 :			Wgrowth[i] = (*wgrowth[inarea])[t][i];
650 :			if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
651 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
652 :			}
653 :			} else {
654 :	ulcessvp	4	for (i = 0; i < size; i++) {
655 :	agomez	1	Lgrowth[i] = kval * pow(LgrpDiv->meanLength(i), growthPar[0]);
656 :			Wgrowth[i] = (*wgrowth[inarea])[t][i];
657 :			if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
658 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
659 :			}
660 :			}
661 :			}
662 :
663 :			// ********************************************************
664 :			// Functions for GrowthCalcH
665 :			// ********************************************************
666 :			GrowthCalcH::GrowthCalcH(CommentStream& infile, const IntVector& Areas,
667 :			const TimeClass* const TimeInfo, Keeper* const keeper)
668 :			: GrowthCalcBase(Areas), numGrowthConstants(4) {
669 :
670 :			keeper->addString("growthcalcH");
671 :			growthPar.resize(numGrowthConstants, keeper);
672 :
673 :			char text[MaxStrLength];
674 :			strncpy(text, "", MaxStrLength);
675 :			infile >> text >> ws;
676 :			if (strcasecmp(text, "growthparameters") != 0)
677 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
678 :			growthPar.read(infile, TimeInfo, keeper);
679 :			keeper->clearLast();
680 :			}
681 :
682 :			/* Simplified 2 parameter length based Von Bertalanffy growth function
683 :			* compare with GrowthCalcC for the more complex weight based version */
684 :	ulcessvp	2	void GrowthCalcH::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
685 :			const PopInfoVector& numGrow, const AreaClass* const Area,
686 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
687 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
688 :
689 :			growthPar.Update(TimeInfo);
690 :			//JMB - first some error checking
691 :			if ((handle.getLogLevel() >= LOGWARN) && (growthPar.didChange(TimeInfo))) {
692 :			if (isZero(growthPar[1]) || isZero(growthPar[2]))
693 :			handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is zero");
694 :			if (LgrpDiv->maxLength() > growthPar[0])
695 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length greater than length infinity");
696 :			}
697 :
698 :			double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
699 :			int i;
700 :			for (i = 0; i < size; i++)
701 :			Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;
702 :			}
703 :
704 :
705 :
706 :
707 :
708 :	agomez	1	// ********************************************************
709 :			// Functions for GrowthCalcI
710 :			// ********************************************************
711 :			GrowthCalcI::GrowthCalcI(CommentStream& infile, const IntVector& Areas,
712 :			const TimeClass* const TimeInfo, Keeper* const keeper)
713 :			: GrowthCalcBase(Areas), numGrowthConstants(6) {
714 :
715 :			keeper->addString("growthcalcI");
716 :			growthPar.resize(numGrowthConstants, keeper);
717 :
718 :			char text[MaxStrLength];
719 :			strncpy(text, "", MaxStrLength);
720 :			infile >> text >> ws;
721 :			if (strcasecmp(text, "growthparameters") != 0)
722 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
723 :			growthPar.read(infile, TimeInfo, keeper);
724 :			keeper->clearLast();
725 :			}
726 :
727 :			/* Simplified 4 parameter Jones growth function
728 :			* compare with GrowthCalcD for the more complex version */
729 :	ulcessvp	4	void GrowthCalcI::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
730 :			const PopInfoVector& numGrow, const AreaClass* const Area,
731 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
732 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
733 :	agomez	1
734 :			growthPar.Update(TimeInfo);
735 :			//JMB - first some error checking
736 :			if ((handle.getLogLevel() >= LOGWARN) && (growthPar.didChange(TimeInfo))) {
737 :			if (isZero(growthPar[0]) || isZero(growthPar[1]))
738 :			handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is zero");
739 :			if (isZero(growthPar[4]) || isZero(growthPar[5]))
740 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length growth parameter is zero");
741 :			}
742 :
743 :			double tempC = TimeInfo->getTimeStepSize() * growthPar[0];
744 :			double tempW = TimeInfo->getTimeStepSize() * growthPar[1] *
745 :			exp(growthPar[3] * Area->getTemperature(area, TimeInfo->getTime()));
746 :
747 :			int i;
748 :	ulcessvp	4	for (i = 0; i < size; i++) {
749 :	agomez	1	if (numGrow[i].W < verysmall) {
750 :			Wgrowth[i] = 0.0;
751 :			Lgrowth[i] = 0.0;
752 :			} else {
753 :			Wgrowth[i] = tempC * Fphi[i] * MaxCon[i] - tempW * pow(numGrow[i].W, growthPar[2]);
754 :
755 :			if (Wgrowth[i] < verysmall) {
756 :			Wgrowth[i] = 0.0;
757 :			Lgrowth[i] = 0.0;
758 :			} else {
759 :			Lgrowth[i] = Wgrowth[i] / (growthPar[4] * growthPar[5] *
760 :			pow(LgrpDiv->meanLength(i), growthPar[5] - 1.0));
761 :			}
762 :			}
763 :			}
764 :			}
765 :
766 :			// ********************************************************
767 :			// Functions for GrowthCalcJ
768 :			// ********************************************************
769 :			GrowthCalcJ::GrowthCalcJ(CommentStream& infile, const IntVector& Areas,
770 :			const TimeClass* const TimeInfo, Keeper* const keeper)
771 :			: GrowthCalcBase(Areas), numGrowthConstants(5) {
772 :
773 :			keeper->addString("GrowthCalcJ");
774 :			growthPar.resize(numGrowthConstants, keeper);
775 :
776 :			char text[MaxStrLength];
777 :			strncpy(text, "", MaxStrLength);
778 :			infile >> text >> ws;
779 :			if (strcasecmp(text, "growthparameters") != 0)
780 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
781 :			growthPar.read(infile, TimeInfo, keeper);
782 :			keeper->clearLast();
783 :			}
784 :
785 :			/* Simplified 2 parameter length based Von Bertalanffy growth function
786 :			* compare with GrowthCalcC for the more complex weight based version
787 :			* with a non-zero value for t0 (compare to GrowthCalcH for simpler version */
788 :	ulcessvp	4	void GrowthCalcJ::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
789 :			const PopInfoVector& numGrow, const AreaClass* const Area,
790 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
791 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
792 :	agomez	1
793 :			growthPar.Update(TimeInfo);
794 :			//JMB - first some error checking
795 :			if ((handle.getLogLevel() >= LOGWARN) && (growthPar.didChange(TimeInfo))) {
796 :			if (isZero(growthPar[1]) || isZero(growthPar[2]))
797 :			handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is zero");
798 :			if (LgrpDiv->maxLength() > growthPar[0])
799 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length greater than length infinity");
800 :			}
801 :
802 :			double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
803 :			int i;
804 :	ulcessvp	4	for (i = 0; i < size; i++)
805 :	agomez	1	Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;
806 :			}
807 :
808 :			// ********************************************************
809 :			// Functions for GrowthCalcK
810 :			// ********************************************************
811 :			GrowthCalcK::GrowthCalcK(CommentStream& infile, const IntVector& Areas,
812 :			const TimeClass* const TimeInfo, Keeper* const keeper)
813 :			: GrowthCalcBase(Areas), numGrowthConstants(5) {
814 :
815 :			keeper->addString("GrowthCalcK");
816 :			growthPar.resize(numGrowthConstants, keeper);
817 :
818 :			char text[MaxStrLength];
819 :			strncpy(text, "", MaxStrLength);
820 :			infile >> text >> ws;
821 :			if (strcasecmp(text, "growthparameters") != 0)
822 :			handle.logFileUnexpected(LOGFAIL, "growthparameters", text);
823 :			growthPar.read(infile, TimeInfo, keeper);
824 :			keeper->clearLast();
825 :			}
826 :
827 :			/* Simplified length based Gompertz growth function */
828 :	ulcessvp	4	void GrowthCalcK::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
829 :			const PopInfoVector& numGrow, const AreaClass* const Area,
830 :			const TimeClass* const TimeInfo, const DoubleVector& Fphi,
831 :			const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
832 :	agomez	1
833 :			growthPar.Update(TimeInfo);
834 :			//JMB - first some error checking
835 :			if ((handle.getLogLevel() >= LOGWARN) && (growthPar.didChange(TimeInfo))) {
836 :			if (isZero(growthPar[1]) || isZero(growthPar[2]))
837 :			handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is zero");
838 :			if (LgrpDiv->maxLength() > growthPar[0])
839 :			handle.logMessage(LOGWARN, "Warning in growth calculation - length greater than length infinity");
840 :			}
841 :
842 :			double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
843 :			int i;
844 :	ulcessvp	4	for (i = 0; i < size; i++)
845 :	agomez	1	Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;
846 :			}

---

root@forge.cesga.es	ViewVC Help
Powered by ViewVC 1.0.0