/*
 *  find_center.c
 *  copyright 2002 Scott Hotton
 *
 *  This program finds the "phyllotactic center" for a finite set of points.  
 *  It does this by implementing gradient descent via Euler's method on a 
 *  gradient flow.  The gradient is formed from the "sum of squares" of the 
 *  divergence angles made by the points.
 */
#include<math.h>
#include<stdio.h>

main()
{
 int j, J=0, count=0;
 double x[501], y[501], sum_x=0.0, sum_y=0.0, dist_sq;
 double theta[501], Dx_theta[501], Dy_theta[501];      /* assumes a maximum */
 double delta[500], Dx_delta[500], Dy_delta[500];      /* of 500 data points */
 double sum_delta;
 double sum_Dx_delta, sum_delta_Dx_delta;
 double sum_Dy_delta, sum_delta_Dy_delta;
 double Dx_S, Dy_S, epsilon=0.00001, tolerance=0.00000000001;

 for( j=1 ; j<=500 ; ++j){             /* reads in coordinates of each point */
    if( scanf("%lf",x+j)==EOF ) break;  
    sum_x+=x[j];                      /* ignores everything after 500 points */
    if( scanf("%lf",y+j)==EOF ){      /* checks there is a y-coordinate for  */
      fprintf(stderr,"last line incomplete\n");  /* for each x-coordinate    */
      exit(1);
    }
    sum_y+=y[j];
    J=j;
 }

 x[0]=sum_x/J;            /* use center of mass as starting estimate */
 y[0]=sum_y/J;

 printf("Starting estimate for center: (%lf, %lf)\n",x[0],y[0]);

 do{

   for( j=1 ; j<=J ; j++){
      theta[j]=atan2(y[j]-y[0],x[j]-x[0]);
      dist_sq=(x[j]-x[0])*(x[j]-x[0])+(y[j]-y[0])*(y[j]-y[0]);
      Dx_theta[j]=(y[j]-y[0])/dist_sq;
      Dy_theta[j]=(x[0]-x[j])/dist_sq;
   }
   sum_delta=0.0;
   sum_Dx_delta=sum_Dy_delta=sum_delta_Dx_delta=sum_delta_Dy_delta=0.0;
   for( j=1 ; j<J ; j++){
      delta[j]=theta[j+1]-theta[j];
 
       /* put delta[j] into the interval [-pi,pi] */

      while( delta[j]<-M_PI ) delta[j]+=2*M_PI;
      while( M_PI<delta[j]  ) delta[j]-=2*M_PI;
    
      sum_delta += delta[j];

      Dx_delta[j] = Dx_theta[j+1]-Dx_theta[j];  /* intermediate terms in the */
      sum_Dx_delta += Dx_delta[j];              /* computation of the grad */
      sum_delta_Dx_delta += delta[j]*Dx_delta[j];

      Dy_delta[j]=Dy_theta[j+1]-Dy_theta[j];
      sum_Dy_delta+=Dy_delta[j];
      sum_delta_Dy_delta+=delta[j]*Dy_delta[j];
   }

   Dx_S = (sum_delta_Dx_delta) - (sum_delta)*(sum_Dx_delta)/(J-1);
   Dy_S = (sum_delta_Dy_delta) - (sum_delta)*(sum_Dy_delta)/(J-1);
   x[0] -= epsilon*Dx_S;
   y[0] -= epsilon*Dy_S;  /* update estimate for the center */

   ++count;                            /* output estimate for center every */
   if( count%1000000==0 ){             /* million iterations */
     printf("Estimate for center: (%lf, %lf)",x[0],y[0]);
     printf(" after %d iterations\n",count);
   }

  /* run this do-loop while the change to the estimate is large enough */

 }while( fabs(epsilon*Dx_S)>tolerance || fabs(epsilon*Dy_S)>tolerance );

 printf("Final estimate for center: (%lf, %lf)\n",x[0],y[0]);
}