#include <iostream.h>
#include <math.h>

struct point {
  double z, y;
};


struct rectangle {
  double z_size, y_size;
  double area;
  point centre;
};


double square (double x) {
  return x * x;
}


double cube (double x) {
  return x * x * x;
}


void calculate_moments (rectangle rectangles[], int n,
                        point &centroid,
                        double &Iz, double &Iy, double &Iyz) {

  double total_area, z_total, y_total;
  int i;

  // compute location of the centroid

  total_area = z_total = y_total = 0;

  for (i = 0; i < n; i++) {
    z_total += rectangles[i].centre.z * rectangles[i].area;
    y_total += rectangles[i].centre.y * rectangles[i].area;
    total_area += rectangles[i].area;
  }
  centroid.z = z_total / total_area;
  centroid.y = y_total / total_area;

  // compute Iz, Iy, and Iyz
  
  Iz = Iy = Iyz = 0;

  for (i = 0; i < n; i++) {
    Iz += ((rectangles[i].z_size * cube(rectangles[i].y_size)) / 12) +
          (rectangles[i].area * square(rectangles[i].centre.y - centroid.y));
    Iy += ((rectangles[i].y_size * cube(rectangles[i].z_size)) / 12) +
          (rectangles[i].area * square(rectangles[i].centre.z - centroid.z));
    Iyz += rectangles[i].area * (rectangles[i].centre.z - centroid.z) *
                                (rectangles[i].centre.y - centroid.y);
  }

}


int main () {

  const int max_rectangles = 40;
  rectangle rectangles[max_rectangles];
  double Iz, Iy, Iyz;
  int n = 0; // number of rectangles
  point centroid;

  for (;;) {

    cout << "\nEnter z and y size of rectangle (0, 0 to quit): ";
    cin >> rectangles[n].z_size >>  rectangles[n].y_size;

    if ((rectangles[n].z_size <= 0) ||
        (rectangles[n].y_size <= 0)) {
      break;
    }

    // compute area
    rectangles[n].area = rectangles[n].z_size * rectangles[n].y_size;

    cout << "Enter z and y co-ordinates of rectangle centre: ";
    cin >> rectangles[n].centre.z >>  rectangles[n].centre.y;
;
    if (++n == max_rectangles) {
      cout << "\n** Max nember of rectangles reached. **\n";
      break;
    }

  }

  // do the calculations
  calculate_moments (rectangles, n, centroid, Iz, Iy, Iyz);

  // output the results
  cout << "\nCentroid z = " << centroid.z
       << "Centroid y = " << centroid.y << endl
       << "\nIz = " << Iz << " Iy = " << Iy << " Iyz = " << Iyz << endl;


  return 0;

}