Log In | Get Help   
Home My Page Projects Code Snippets Project Openings Mareframe
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files
[mareframe] Diff of /trunk/gadget/growthcalc.cc
[mareframe] / trunk / gadget / growthcalc.cc Repository:
ViewVC logotype

Diff of /trunk/gadget/growthcalc.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3, Thu Apr 30 12:52:16 2015 UTC revision 4, Thu Apr 30 17:32:47 2015 UTC
# Line 32  Line 32 
32    keeper->clearLast();    keeper->clearLast();
33  }  }
34    
35  void GrowthCalcA::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcA::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
36    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
37    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
38    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
39    
40    growthPar.Update(TimeInfo);    growthPar.Update(TimeInfo);
41    double tempL = TimeInfo->getTimeStepSize() * growthPar[0] *    double tempL = TimeInfo->getTimeStepSize() * growthPar[0] *
# Line 44  Line 44 
44        (growthPar[7] * Area->getTemperature(area, TimeInfo->getTime()) + growthPar[8]);        (growthPar[7] * Area->getTemperature(area, TimeInfo->getTime()) + growthPar[8]);
45    
46    int i;    int i;
47    for (i = 0; i < Lgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
48      Lgrowth[i] = tempL * pow(LgrpDiv->meanLength(i), growthPar[1]) * Fphi[i];      Lgrowth[i] = tempL * pow(LgrpDiv->meanLength(i), growthPar[1]) * Fphi[i];
49      if (Lgrowth[i] < 0.0)      if (Lgrowth[i] < 0.0)
50        Lgrowth[i] = 0.0;        Lgrowth[i] = 0.0;
# Line 111  Line 111 
111    }    }
112  }  }
113    
114  void GrowthCalcB::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcB::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
115    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
116    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
117    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
118    
119    int i, inarea = this->areaNum(area);    int i, inarea = this->areaNum(area);
120    for (i = 0; i < Lgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
121      Lgrowth[i] = (*lgrowth[inarea])[TimeInfo->getTime()][i];      Lgrowth[i] = (*lgrowth[inarea])[TimeInfo->getTime()][i];
122      Wgrowth[i] = (*wgrowth[inarea])[TimeInfo->getTime()][i];      Wgrowth[i] = (*wgrowth[inarea])[TimeInfo->getTime()][i];
123      if ((handle.getLogLevel() >= LOGWARN) && ((Lgrowth[i] < 0.0) || (Wgrowth[i] < 0.0)))      if ((handle.getLogLevel() >= LOGWARN) && ((Lgrowth[i] < 0.0) || (Wgrowth[i] < 0.0)))
# Line 192  Line 192 
192   * final form of the function is   * final form of the function is
193   * dw/dt = a0*exp(a1*T)*((w/a2)^a4 - (w/a3)^a5)   * dw/dt = a0*exp(a1*T)*((w/a2)^a4 - (w/a3)^a5)
194   * For no temperature dependency a1 = 0 */   * For no temperature dependency a1 = 0 */
195  void GrowthCalcC::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcC::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
196    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
197    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
198    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
199    
200    wgrowthPar.Update(TimeInfo);    wgrowthPar.Update(TimeInfo);
201    lgrowthPar.Update(TimeInfo);    lgrowthPar.Update(TimeInfo);
# Line 216  Line 216 
216    double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[0] *    double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[0] *
217        exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));        exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));
218    
219    for (i = 0; i < Wgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
220      if (numGrow[i].W < verysmall || isZero(tempW)) {      if (numGrow[i].W < verysmall || isZero(tempW)) {
221        Wgrowth[i] = 0.0;        Wgrowth[i] = 0.0;
222        Lgrowth[i] = 0.0;        Lgrowth[i] = 0.0;
# Line 308  Line 308 
308   * part is derived from the weight increase part by assuming a formula   * 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   * 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. */   * but instead the weight increases until it reaches the curve. */
311  void GrowthCalcD::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcD::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
312    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
313    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
314    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
315    
316    wgrowthPar.Update(TimeInfo);    wgrowthPar.Update(TimeInfo);
317    lgrowthPar.Update(TimeInfo);    lgrowthPar.Update(TimeInfo);
# Line 332  Line 332 
332    double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[2] *    double tempW = TimeInfo->getTimeStepSize() * wgrowthPar[2] *
333        exp(wgrowthPar[4] * Area->getTemperature(area, TimeInfo->getTime()) + wgrowthPar[5]);        exp(wgrowthPar[4] * Area->getTemperature(area, TimeInfo->getTime()) + wgrowthPar[5]);
334    
335    for (i = 0; i < Wgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
336      if (numGrow[i].W < verysmall) {      if (numGrow[i].W < verysmall) {
337        Wgrowth[i] = 0.0;        Wgrowth[i] = 0.0;
338        Lgrowth[i] = 0.0;        Lgrowth[i] = 0.0;
# Line 456  Line 456 
456   *       AreaEffect   *       AreaEffect
457   *       StepEffect   *       StepEffect
458   * Length increase is upgraded in the same way as earlier. */   * Length increase is upgraded in the same way as earlier. */
459  void GrowthCalcE::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcE::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
460    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
461    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
462    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
463    
464    wgrowthPar.Update(TimeInfo);    wgrowthPar.Update(TimeInfo);
465    lgrowthPar.Update(TimeInfo);    lgrowthPar.Update(TimeInfo);
# Line 482  Line 482 
482    double tempW = factor * TimeInfo->getTimeStepSize() * wgrowthPar[0] *    double tempW = factor * TimeInfo->getTimeStepSize() * wgrowthPar[0] *
483        exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));        exp(wgrowthPar[1] * Area->getTemperature(area, TimeInfo->getTime()));
484    
485    for (i = 0; i < Wgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
486      if (numGrow[i].W < verysmall || isZero(tempW)) {      if (numGrow[i].W < verysmall || isZero(tempW)) {
487        Wgrowth[i] = 0.0;        Wgrowth[i] = 0.0;
488        Lgrowth[i] = 0.0;        Lgrowth[i] = 0.0;
# Line 555  Line 555 
555      delete wgrowth[a];      delete wgrowth[a];
556  }  }
557    
558  void GrowthCalcF::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcF::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
559    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
560    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
561    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
562    
563    growthPar.Update(TimeInfo);    growthPar.Update(TimeInfo);
564    int i, t, inarea;    int i, t, inarea;
# Line 566  Line 566 
566    inarea = this->areaNum(area);    inarea = this->areaNum(area);
567    double kval = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());    double kval = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
568    
569    for (i = 0; i < Lgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
570      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * kval;      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * kval;
571      Wgrowth[i] = (*wgrowth[inarea])[t][i];      Wgrowth[i] = (*wgrowth[inarea])[t][i];
572      if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))      if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
# Line 619  Line 619 
619      delete wgrowth[a];      delete wgrowth[a];
620  }  }
621    
622  void GrowthCalcG::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcG::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
623    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
624    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
625    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
626    
627    //written by kgf 24/10 00    //written by kgf 24/10 00
628    //Gives linear growth (growthPar[0] == 0) or    //Gives linear growth (growthPar[0] == 0) or
# Line 637  Line 637 
637      handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is positive");      handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is positive");
638    
639    if (isZero(growthPar[0])) {    if (isZero(growthPar[0])) {
640      for (i = 0; i < Lgrowth.Size(); i++) {      for (i = 0; i < size; i++) {
641        Lgrowth[i] = kval;        Lgrowth[i] = kval;
642        Wgrowth[i] = (*wgrowth[inarea])[t][i];        Wgrowth[i] = (*wgrowth[inarea])[t][i];
643        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
644          handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");          handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
645      }      }
646    } else if (isEqual(growthPar[0], 1.0)) {    } else if (isEqual(growthPar[0], 1.0)) {
647      for (i = 0; i < Lgrowth.Size(); i++) {      for (i = 0; i < size; i++) {
648        Lgrowth[i] = kval * LgrpDiv->meanLength(i);        Lgrowth[i] = kval * LgrpDiv->meanLength(i);
649        Wgrowth[i] = (*wgrowth[inarea])[t][i];        Wgrowth[i] = (*wgrowth[inarea])[t][i];
650        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
651          handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");          handle.logMessage(LOGWARN, "Warning in growth calculation - weight growth parameter is negative");
652      }      }
653    } else {    } else {
654      for (i = 0; i < Lgrowth.Size(); i++) {      for (i = 0; i < size; i++) {
655        Lgrowth[i] = kval * pow(LgrpDiv->meanLength(i), growthPar[0]);        Lgrowth[i] = kval * pow(LgrpDiv->meanLength(i), growthPar[0]);
656        Wgrowth[i] = (*wgrowth[inarea])[t][i];        Wgrowth[i] = (*wgrowth[inarea])[t][i];
657        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))        if ((handle.getLogLevel() >= LOGWARN) && (Wgrowth[i] < 0.0))
# Line 681  Line 681 
681    
682  /* Simplified 2 parameter length based Von Bertalanffy growth function  /* Simplified 2 parameter length based Von Bertalanffy growth function
683   * compare with GrowthCalcC for the more complex weight based version */   * compare with GrowthCalcC for the more complex weight based version */
 void GrowthCalcH::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  
   const PopInfoVector& numGrow, const AreaClass* const Area,  
   const TimeClass* const TimeInfo, const DoubleVector& Fphi,  
   const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {  
   
   growthPar.Update(TimeInfo);  
   //JMB - first some error checking  
   if ((handle.getLogLevel() >= LOGWARN) && (growthPar.didChange(TimeInfo))) {  
     if (isZero(growthPar[1]) || isZero(growthPar[2]))  
       handle.logMessage(LOGWARN, "Warning in growth calculation - growth parameter is zero");  
     if (LgrpDiv->maxLength() > growthPar[0])  
       handle.logMessage(LOGWARN, "Warning in growth calculation - length greater than length infinity");  
   }  
   
   double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());  
   int i;  
   for (i = 0; i < Lgrowth.Size(); i++)  
     Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;  
 }  
   
   
   
684  void GrowthCalcH::calcGrowth(int area, double* Lgrowth, double* Wgrowth,  void GrowthCalcH::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
685    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
686    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
# Line 748  Line 726 
726    
727  /* Simplified 4 parameter Jones growth function  /* Simplified 4 parameter Jones growth function
728   * compare with GrowthCalcD for the more complex version */   * compare with GrowthCalcD for the more complex version */
729  void GrowthCalcI::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcI::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
730    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
731    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
732    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
733    
734    growthPar.Update(TimeInfo);    growthPar.Update(TimeInfo);
735    //JMB - first some error checking    //JMB - first some error checking
# Line 767  Line 745 
745        exp(growthPar[3] * Area->getTemperature(area, TimeInfo->getTime()));        exp(growthPar[3] * Area->getTemperature(area, TimeInfo->getTime()));
746    
747    int i;    int i;
748    for (i = 0; i < Wgrowth.Size(); i++) {    for (i = 0; i < size; i++) {
749      if (numGrow[i].W < verysmall) {      if (numGrow[i].W < verysmall) {
750        Wgrowth[i] = 0.0;        Wgrowth[i] = 0.0;
751        Lgrowth[i] = 0.0;        Lgrowth[i] = 0.0;
# Line 807  Line 785 
785  /* Simplified 2 parameter length based Von Bertalanffy growth function  /* Simplified 2 parameter length based Von Bertalanffy growth function
786   * compare with GrowthCalcC for the more complex weight based version   * compare with GrowthCalcC for the more complex weight based version
787   * with a non-zero value for t0 (compare to GrowthCalcH for simpler version */   * with a non-zero value for t0 (compare to GrowthCalcH for simpler version */
788  void GrowthCalcJ::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcJ::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
789    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
790    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
791    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
792    
793    growthPar.Update(TimeInfo);    growthPar.Update(TimeInfo);
794    //JMB - first some error checking    //JMB - first some error checking
# Line 823  Line 801 
801    
802    double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());    double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
803    int i;    int i;
804    for (i = 0; i < Lgrowth.Size(); i++)    for (i = 0; i < size; i++)
805      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;
806  }  }
807    
# Line 847  Line 825 
825  }  }
826    
827  /* Simplified length based Gompertz growth function */  /* Simplified length based Gompertz growth function */
828  void GrowthCalcK::calcGrowth(int area, DoubleVector& Lgrowth, DoubleVector& Wgrowth,  void GrowthCalcK::calcGrowth(int area, double* Lgrowth, double* Wgrowth,
829    const PopInfoVector& numGrow, const AreaClass* const Area,    const PopInfoVector& numGrow, const AreaClass* const Area,
830    const TimeClass* const TimeInfo, const DoubleVector& Fphi,    const TimeClass* const TimeInfo, const DoubleVector& Fphi,
831    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv) {                    const DoubleVector& MaxCon, const LengthGroupDivision* const LgrpDiv, int size) {
832    
833    growthPar.Update(TimeInfo);    growthPar.Update(TimeInfo);
834    //JMB - first some error checking    //JMB - first some error checking
# Line 863  Line 841 
841    
842    double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());    double mult = 1.0 - exp(-growthPar[1] * TimeInfo->getTimeStepSize());
843    int i;    int i;
844    for (i = 0; i < Lgrowth.Size(); i++)    for (i = 0; i < size; i++)
845      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;      Lgrowth[i] = (growthPar[0] - LgrpDiv->meanLength(i)) * mult;
846  }  }
Legend:
Removed from v.3  
changed lines
  Added in v.4
root@forge.cesga.es
 ViewVC Help
Powered by ViewVC 1.0.0  

Powered By FusionForge