#include <vcl.h>
#include <math.h>
#include <stdio.h>

#include "Data_Gen.h"
#include "ListConst.h"
#include "ProcExt.h"
#include "message.h"

//---------------------------------------------------------
// Вычислить геометрические параметры сечения тонкостнного профиля по точкам
int Get_geom_param_XY(double *X_, double* Y_, double* thick, double hG, int n_seg, GEOM_PARAM_SECT* param_sect) {

  double lk;
  double Ar, dA[20];
  double Sx, Sy, Xc, Yc;
  double Ix, Iy, Ixy;
  double x_k[20], y_k[20];
  double x_[20], y_[20];
  int i, k;
  double alfa, Iu, Iv;
  double Itor;
  double omega[20];
  double Iwx, Iwy, ax, ay, aX, aY;
  double W_, Iw;

  Ar = 0;
  Sx = Sy = 0;
  Itor = 0;
  for (k=1; k<=n_seg; k++) {
    lk = sqrt((X_[k] - X_[k-1])*(X_[k] - X_[k-1]) + (Y_[k] - Y_[k-1])*(Y_[k] - Y_[k-1]));
    dA[k] = lk*thick[k];
    Ar += dA[k];
    Sx += 0.5*(Y_[k] + Y_[k-1])*dA[k];
    Sy += 0.5*(X_[k] + X_[k-1])*dA[k];
    Itor += lk*pow(thick[k], 3.0)/3;
  }
  // Координаты ц.т.
  Xc = Sy/Ar;  Yc = Sx/Ar;
  for (k=0; k<=n_seg; k++) {
    x_k[k] = X_[k] - Xc;
    y_k[k] = Y_[k] - Yc;
  }
  Ix = Iy = Ixy = 0;
  for (k=1; k<=n_seg; k++) {
    Ix += (y_k[k-1]*y_k[k-1] + y_k[k-1]*y_k[k] + y_k[k]*y_k[k])*dA[k]/3.0;
    Iy += (x_k[k-1]*x_k[k-1] + x_k[k-1]*x_k[k] + x_k[k]*x_k[k])*dA[k]/3.0;
    Ixy += (x_k[k-1]*y_k[k-1] + x_k[k]*y_k[k])*dA[k]/3.0 + (x_k[k-1]*y_k[k] + x_k[k]*y_k[k-1])*dA[k]/6.0;
  }
  alfa = -0.5*atan(2*Ixy/(Ix - Iy));
  Iu = 0.5*(Ix + Iy) + 0.5*(Ix - Iy)*cos(2*alfa) - Ixy*sin(2*alfa);
  Iv = Ix + Iy - Iu;
  // Секториальные координаты
  for (i=0; i<=n_seg; i++) {
     omega[i] = 0;
     for (k=1; k<=i; k++) {
       omega[i] += x_k[k-1]*y_k[k] - x_k[k]*y_k[k-1];
     }
  }
  // Секториальные центробежные моменты инерции
  Iwx = Iwy = 0;
  for (k=1; k<=n_seg; k++) {
    Iwx += (omega[k-1]*x_k[k-1] + omega[k]*x_k[k])*dA[k]/3 + (omega[k-1]*x_k[k] + omega[k]*x_k[k-1])*dA[k]/6;
    Iwy += (omega[k-1]*y_k[k-1] + omega[k]*y_k[k])*dA[k]/3 + (omega[k-1]*y_k[k] + omega[k]*y_k[k-1])*dA[k]/6;
  }
  // Координаты центра сдвига
  aX = Xc + (Ixy*Iwx - Iy*Iwy)/(Ixy*Ixy - Ix*Iy);
  aY = Yc - (Ixy*Iwy - Ix*Iwx)/(Ixy*Ixy - Ix*Iy);
  // Координаты центра сдвига в главных центральных осях
  ax = (aX - Xc)*cos(alfa) + (aY - Yc)*sin(alfa) ;
  ay = -(aX - Xc)*sin(alfa) + (aY - Yc)*cos(alfa) ;
  // Координаты узловых точек в главных центральных осях

  double x0_max = -1e6;
  double x0_min = 1e6;
  double y0_max = -1e6;
  double y0_min = 1e6;

  for (k=0; k<=n_seg; k++) {
    x_[k] = x_k[k]*cos(alfa) + y_k[k]*sin(alfa);
    y_[k] = -x_k[k]*sin(alfa) + y_k[k]*cos(alfa);
    x0_max = MAX(x0_max, x_[k]);
    y0_max = MAX(y0_max, y_[k]);
    x0_min = MIN(x0_min, x_[k]);
    y0_min = MIN(y0_min, y_[k]);
  }
  // Секториальные координаты относительно центра сдвига
  for (i=0; i<=n_seg; i++) {
     omega[i] = 0;
     for (k=1; k<=i; k++) {
       omega[i] += (X_[k-1] - aX)*(Y_[k] - aY) - (X_[k] - aX)*(Y_[k-1] - aY);
     }
  }
  // Нормализованные секториальные координаты
  W_ = 0;
  for (k=1; k<=n_seg; k++) {
    W_ += 0.5*(omega[k-1] + omega[k])*dA[k];
  }
  for (k=0; k<=n_seg; k++) {
    omega[k] -= W_/Ar;
  }
  // Cекториальный момент инерции
  Iw = 0;
  for (k=1; k<=n_seg; k++) {
    Iw += (omega[k-1]*omega[k-1] + omega[k-1]*omega[k] + omega[k]*omega[k])*dA[k]/3;
  }
  param_sect->Ar = Round(Ar, 1);
  param_sect->Iyy = Round(Ix, 0);   // относительно горизонтальной оси
  param_sect->Izz = Round(Iy, 0);
  param_sect->Iyz = Round(Ixy, 0);
  param_sect->Itor = Round(Itor, 0);
  param_sect->Zct = Round(hG/2 - Yc, 2);
  param_sect->Yct = Round(Xc, 2);
  param_sect->Ze = Round(hG/2 - aY, 2);
  param_sect->Ye = Round(aX, 2);
  param_sect->alfa = alfa;
  param_sect->I1 = Round(Iu, 0);
  param_sect->I2 = Round(Iv, 0);
  param_sect->Iw = Round(Iw, 0);
  if ((ABS(alfa)) < PI_/4) {
    param_sect->Wu_max = ABS(Iu/y0_max);
    param_sect->Wu_min = ABS(Iu/y0_min);
    param_sect->Wv_max = ABS(Iv/x0_max);
    param_sect->Wv_min = ABS(Iv/x0_min);
  }
  else {
    param_sect->Wu_max = ABS(Iu/x0_max);
    param_sect->Wu_min = ABS(Iu/x0_min);
    param_sect->Wv_max = ABS(Iv/y0_max);
    param_sect->Wv_min = ABS(Iv/y0_min);
  }

}
