--- trunk/gadget/grow.cc	2014/02/10 17:09:07	1
+++ trunk/gadget/grow.cc	2015/04/29 12:55:30	2
@@ -2,31 +2,135 @@
 #include "grower.h"
 #include "agebandmatrix.h"
 #include "gadget.h"
+#include "doublematrix.h"
+#include "matrix.h"
 
 /* Update the agebandmatrix to reflect the calculated growth  */
 /* Lgrowth contains the ratio of each length group that grows */
 /* by a certain number of length groups, and Wgrowth contains */
 /* the weight increase for each entry in Lgrowth              */
 
+
+  void AgeBandMatrix::Grow(const Matrix& Lgrowth, const Matrix& Wgrowth) {
+    int i, lgrp, grow, maxlgrp;
+    double num, wt, tmp;
+//    ofstream outfile;
+//      outfile.open("aaa", ios::out);
+//    _Lgrowth.Print(outfile);
+//    int a;
+//    cout << "1-----------------" << endl;
+//    cin >> a;
+
+
+    maxlgrp = Lgrowth.Nrow();
+
+
+  //   double** lGrowth = ToMatrix(Lgrowth, maxlgrp, aux);
+  //   double** wGrowth = ToMatrix(Wgrowth, maxlgrp, aux);
+
+
+    for (i = 0; i < nrow; i++) {
+  	int maxCol = v[i]->maxCol();
+  	int minCol = v[i]->minCol();
+      //the part that grows to or above the highest length group
+      num = 0.0;
+      wt = 0.0;
+
+      for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
+        for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
+          tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
+          num += tmp;
+          wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
+        }
+      }
+//      cout << "wt1:" << wt << endl;
+
+      lgrp = maxCol - 1;
+      if (isZero(num) || (wt < verysmall)) {
+        (*v[i])[lgrp].setToZero();
+      } else {
+        (*v[i])[lgrp].W = wt / num;
+        (*v[i])[lgrp].N = num;
+      }
+
+
+      //the central diagonal part of the length division
+      for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
+        num = 0.0;
+        wt = 0.0;
+        for (grow = 0; grow < maxlgrp; grow++) {
+          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
+          num += tmp;
+          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
+        }
+//        cout << "wt2:" << wt << endl;
+        if (isZero(num) || (wt < verysmall)) {
+          (*v[i])[lgrp].setToZero();
+        } else {
+          (*v[i])[lgrp].W = wt / num;
+          (*v[i])[lgrp].N = num;
+        }
+      }
+
+
+      //the lowest part of the length division
+      for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
+        num = 0.0;
+        wt = 0.0;
+        for (grow = 0; grow <= lgrp - minCol; grow++) {
+          tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
+          num += tmp;
+          wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
+        }
+
+        if (isZero(num) || (wt < verysmall)) {
+          (*v[i])[lgrp].setToZero();
+        } else {
+          (*v[i])[lgrp].W = wt / num;
+          (*v[i])[lgrp].N = num;
+        }
+      }
+//      cout << "wt3:" << wt << endl;
+    }
+  }
+
+
+
+
+
+
+
+
+
+  /********************************************************************************/
+
 /* JMB changed to deal with very small weights a bit better   */
 void AgeBandMatrix::Grow(const DoubleMatrix& Lgrowth, const DoubleMatrix& Wgrowth) {
   int i, lgrp, grow, maxlgrp;
   double num, wt, tmp;
 
   maxlgrp = Lgrowth.Nrow();
+
+//   double** lGrowth = ToMatrix(Lgrowth, maxlgrp, aux);
+//   double** wGrowth = ToMatrix(Wgrowth, maxlgrp, aux);
+
+
   for (i = 0; i < nrow; i++) {
+	int maxCol = v[i]->maxCol();
+	int minCol = v[i]->minCol();
     //the part that grows to or above the highest length group
     num = 0.0;
     wt = 0.0;
-    for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {
-      for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {
+
+    for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
+      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
         tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
         num += tmp;
         wt += tmp * (Wgrowth[grow][lgrp] + (*v[i])[lgrp].W);
       }
     }
 
-    lgrp = v[i]->maxCol() - 1;
+    lgrp = maxCol - 1;
     if (isZero(num) || (wt < verysmall)) {
       (*v[i])[lgrp].setToZero();
     } else {
@@ -34,8 +138,9 @@
       (*v[i])[lgrp].N = num;
     }
 
+
     //the central diagonal part of the length division
-    for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {
+    for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
       num = 0.0;
       wt = 0.0;
       for (grow = 0; grow < maxlgrp; grow++) {
@@ -53,10 +158,10 @@
     }
 
     //the lowest part of the length division
-    for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {
+    for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
       num = 0.0;
       wt = 0.0;
-      for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {
+      for (grow = 0; grow <= lgrp - minCol; grow++) {
         tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
         num += tmp;
         wt += tmp * (Wgrowth[grow][lgrp - grow] + (*v[i])[lgrp - grow].W);
@@ -78,14 +183,17 @@
   int i, lgrp, grow, maxlgrp, age;
   double num, wt, matnum, tmp, ratio;
 
+
   maxlgrp = Lgrowth.Nrow();
   for (i = 0; i < nrow; i++) {
+	  int maxCol = v[i]->maxCol();
+	    int minCol = v[i]->minCol();
     age = i + minage;
     num = 0.0;
     wt = 0.0;
     matnum = 0.0;
-    for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {
-      for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {
+    for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
+      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
         ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);
         tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
         matnum += (tmp * ratio);
@@ -94,7 +202,7 @@
       }
     }
 
-    lgrp = v[i]->maxCol() - 1;
+    lgrp = maxCol - 1;
     if (isZero(num) || (wt < verysmall)) {
       //no fish grow to this length cell
       (*v[i])[lgrp].setToZero();
@@ -114,7 +222,8 @@
       Mat->storeMatureStock(area, age, lgrp, matnum, wt / num);
     }
 
-    for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {
+
+    for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
       num = 0.0;
       wt = 0.0;
       matnum = 0.0;
@@ -146,11 +255,11 @@
       }
     }
 
-    for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {
+    for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
       num = 0.0;
       wt = 0.0;
       matnum = 0.0;
-      for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {
+      for (grow = 0; grow <= lgrp - minCol; grow++) {
         ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);
         tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
         matnum += (tmp * ratio);
@@ -186,13 +295,16 @@
   double num;
 
   maxlgrp = Lgrowth.Nrow();
+
   for (i = 0; i < nrow; i++) {
+	  int maxCol = v[i]->maxCol();
+	    int minCol = v[i]->minCol();
     num = 0.0;
-    for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--)
-      for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++)
+    for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--)
+      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++)
         num += (Lgrowth[grow][lgrp] * (*v[i])[lgrp].N);
 
-    lgrp = v[i]->maxCol() - 1;
+    lgrp = maxCol - 1;
     if (isZero(num)) {
       (*v[i])[lgrp].setToZero();
     } else {
@@ -200,7 +312,7 @@
       (*v[i])[lgrp].W = Weight[lgrp];
     }
 
-    for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {
+    for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
       num = 0.0;
       for (grow = 0; grow < maxlgrp; grow++)
         num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);
@@ -213,9 +325,9 @@
       }
     }
 
-    for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {
+    for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
       num = 0.0;
-      for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++)
+      for (grow = 0; grow <= lgrp - minCol; grow++)
         num += (Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N);
 
       if (isZero(num)) {
@@ -236,12 +348,15 @@
   double num, matnum, tmp, ratio;
 
   maxlgrp = Lgrowth.Nrow();
+
   for (i = 0; i < nrow; i++) {
+	  int maxCol = v[i]->maxCol();
+	    int minCol = v[i]->minCol();
     age = i + minage;
     num = 0.0;
     matnum = 0.0;
-    for (lgrp = v[i]->maxCol() - 1; lgrp >= v[i]->maxCol() - maxlgrp; lgrp--) {
-      for (grow = v[i]->maxCol() - lgrp - 1; grow < maxlgrp; grow++) {
+    for (lgrp = maxCol - 1; lgrp >= maxCol - maxlgrp; lgrp--) {
+      for (grow = maxCol - lgrp - 1; grow < maxlgrp; grow++) {
         ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp].W);
         tmp = Lgrowth[grow][lgrp] * (*v[i])[lgrp].N;
         matnum += (tmp * ratio);
@@ -249,7 +364,7 @@
       }
     }
 
-    lgrp = v[i]->maxCol() - 1;
+    lgrp = maxCol - 1;
     if (isZero(num)) {
       //no fish grow to this length cell
       (*v[i])[lgrp].setToZero();
@@ -269,7 +384,7 @@
       Mat->storeMatureStock(area, age, lgrp, matnum, Weight[lgrp]);
     }
 
-    for (lgrp = v[i]->maxCol() - 2; lgrp >= v[i]->minCol() + maxlgrp - 1; lgrp--) {
+    for (lgrp = maxCol - 2; lgrp >= minCol + maxlgrp - 1; lgrp--) {
       num = 0.0;
       matnum = 0.0;
       for (grow = 0; grow < maxlgrp; grow++) {
@@ -299,10 +414,10 @@
       }
     }
 
-    for (lgrp = v[i]->minCol() + maxlgrp - 2; lgrp >= v[i]->minCol(); lgrp--) {
+    for (lgrp = minCol + maxlgrp - 2; lgrp >= minCol; lgrp--) {
       num = 0.0;
       matnum = 0.0;
-      for (grow = 0; grow <= lgrp - v[i]->minCol(); grow++) {
+      for (grow = 0; grow <= lgrp - minCol; grow++) {
         ratio = Mat->calcMaturation(age, lgrp, grow, (*v[i])[lgrp - grow].W);
         tmp = Lgrowth[grow][lgrp - grow] * (*v[i])[lgrp - grow].N;
         matnum += (tmp * ratio);