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

float Atrap(float n, float m, float x) {
    float A;
    A=(n+m*x/2)*x;
    return A;
}
float Sfun(float n, float m, float x) {
    float S;
    S=x*(3*n*n+x*(3*n*m+x*m*m))/6.;
    return S;
}
float Ifun(float n, float m, float x) {
    float I;
    I=(4*n*n*n*x+6*m*pow(n*x,2)+4*n*x*pow(m*x,2)+x*pow(m*x,3))/12.;
    return I;
}
float Ixyfun(float n, float m, float x) {
    float I;
    I=(6*pow(n*x,2)+8*n*m*pow(x,3)+3*m*m*pow(x,4))/24.;
    return I;
}
void presek(float x[], float y[], int N,
                  float *xT, float *yT, float *A, float *Sx, float *Sy,
                  float *Ix, float *Iy, float *Ixy) {
    int i;
    float n,m;
    *A=0.0; *Sx=0.0; *Sy=0.0;
    *Ix=0.0; *Iy=0.0; *Ixy=0.0;
    for(i=0; i<N-1;i++) {
        if(x[i+1] != x[i]) {
            m=(y[i+1]-y[i])/(x[i+1]-x[i]);
            n=y[i]-m*x[i];
            *A  = *A - Atrap(n,m,x[i+1]) + Atrap(n,m,x[i]);
            *Sx = *Sx - Sfun(n,m,x[i+1]) + Sfun(n,m,x[i]);
            *Ix = *Ix - Ifun(n,m,x[i+1]) + Ifun(n,m,x[i]);
            *Ixy = *Ixy - Ixyfun(n,m,x[i+1]) + Ixyfun(n,m,x[i]);
        }
        if(y[i+1] != y[i]) {
            m=(x[i+1]-x[i])/(y[i+1]-y[i]);
            n=x[i]-m*y[i];
            *Sy = *Sy + Sfun(n,m,y[i+1]) - Sfun(n,m,y[i]);
            *Iy = *Iy + Ifun(n,m,y[i+1]) - Ifun(n,m,y[i]);
        }
    }
    *xT=*Sy/(*A);
    *yT=*Sx/(*A);
}
void Itez(float xT, float yT, float A, float Sx, float Sy,
          float Ix, float Iy, float Ixy,
          float *Jx, float *Jy, float *Jxy) {
    *Jx = Ix - 2*yT*Sx + yT*yT*A;  
    *Jy = Iy - 2*xT*Sy + xT*xT*A;  
    *Jxy = Ixy -xT*Sx - yT*Sy +xT*yT*A; 
}          
void Imax(float Jx, float Jy, float Jxy, 
          float *alph, float *Jmax, float *Jmin) {
    *alph = atan2(2*Jxy,Jy-Jx);
    *Jmax= (Jx+Jy)/2+sqrt(pow((Jy-Jx)/2,2)+pow(Jxy,2));                
    *Jmin= (Jx+Jy)/2-sqrt(pow((Jy-Jx)/2,2)+pow(Jxy,2));  
}              
int main()
{
  float x[12]={0,9,9,10,10.732,11,10.732,10,9,6,6,0};
  float y[12]={0,0,2,2.268,3,4,5,5.732,6,6,3,0};
  float xT,yT,A,Sx,Sy,Ix,Iy,Ixy,Jx,Jy,Jxy,alph,Jmax,Jmin;
  int n=12;
  presek(x,y,n,&xT,&yT,&A,&Sx,&Sy,&Ix,&Iy,&Ixy);
  printf("A, xT, yT = %f, %f, %f\n",A,xT,yT);
  printf("Sx, Sy = %f, %f\n",Sx,Sy);
  printf("Ix, Iy, Ixy = %f, %f, %f\n",Ix,Iy,Ixy);
  Itez(xT,yT,A,Sx,Sy,Ix,Iy,Ixy,&Jx,&Jy,&Jxy);
  printf("Jx, Jy, Jxy = %f, %f, %f\n",Jx,Jy,Jxy);
  Imax(Jx, Jy, Jxy, &alph, &Jmax, &Jmin);
  printf("Jmax, Jmin, alpha = %f, %f, %f\n",Jmax,Jmin,alph*180/3.14159265);
  system("PAUSE");	
  return 0;
}