/***************************** calcinit.c *********************************
                          Xinyu Deng (cindy)
                             Feb.25, 1994
***************************************************************************/
/**************************************************************************
modification: change output format for c1, c2, x0 y0
***************************************************************************/
/* this program calculate initial value of iteration */

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

/* constants definition */
#ifndef PI
#define PI       3.141592653589793238462643383280
#endif

#ifndef ERROR
#define ERROR    -1
#endif


typedef struct curve
{
  int type;       /* ellipse, rose */
  double c1, c2;
  double x0, y0;   /* curve offset center */
} curve;

/*** global variable declaration ***/
int n;       /* number of subdivision point on curve */
int m;       /* number of curves */
int mn;
curve **cs;
double delta_theta;
double *x, *y;
FILE *output1, *output2;

/* function prototype */
void allocate_xy(void);
void initellip(int);
void initellip2(int);
void initellip3(int);
void initellip4(int);
void initrose(int);

void main(int argc, char *argv[])
{
  int i, ans, c_type;
  char buf[100];
  char *str;

  if( argc < 3 )
  { 
     fprintf(stderr, "Usage: %s <output file1> <output file2>\n", argv[0]);
     exit( -1 );
  }

/* open output file */
   if( (output1 = fopen(argv[1], "w") ) == NULL )
   {
     fprintf(stderr,"%s: can't open %s for output1\n", argv[0], argv[1]);
     exit( -1 );
   }

/* open output file */
   if( (output2 = fopen(argv[2], "w") ) == NULL )
   {
     fprintf(stderr,"%s: can't open %s for output2\n", argv[0], argv[2]);
     exit( -1 );
   }

/* input parameters */
   printf("\nEnter number of curves to be used: ");
   fflush( stdin );
   scanf(" %d", &m);

   if( m < 1 )
   {
     fprintf(stderr, "Invalid number of curves: %d\n", m);
     exit( -1 );
   }

   printf("\nEnter number of points to be used: ");
   fflush( stdin );
   scanf(" %d", &n);

   if( n < 1 )
   {
     fprintf(stderr, "Invalid number of points: %d\n", n);
     exit( -1 );
   }
 
   mn = m*n;
   cs = (curve **) calloc(m, sizeof(curve *) );
   for( i=0; i<m; i++)
   {
     cs[i] = (curve *) calloc(1, sizeof( curve ) );
     ans = 0;
     do
     {
       printf("\nPlease choose a curve type for curve #%d:\n", i+1);
       printf("1. ellipse\n2.  long ellipse\n3. long circle\n4. dumbbell\n5. rose\n");
       printf("Enter your choice: ");
       fflush( stdin );
       scanf(" %d", &c_type);
       if( c_type !=1 && c_type != 2 && c_type != 3 && c_type != 4 && c_type != 5)
       {
          printf("Invalid curve type: %d, try again !!\n", c_type);
          ans = 0;
       }
       else
       {
          ans = 1;
       }
     }  while ( ans != 1 );
     cs[i]->type = c_type;
     printf("\nEnter curve parameters(c1, c2) for curve #%d: ", i+1);
     fflush( stdin );
     gets( buf );
     str = strtok( buf, " ,");
     cs[i]->c1 = atof( str );
     str = strtok( NULL, " ,");
     cs[i]->c2 = atof( str );

     printf("\nEnter curve center(x0, y0) for curve #%d: ", i+1);
     fflush( stdin );
     gets( buf );
     str = strtok(buf, " ,");
     cs[i]->x0 = atof( str );
     str = strtok( NULL, " ,");
     cs[i]->y0 = atof( str );
   }

/* calculate delta_theta */
   delta_theta = 2.0 * PI / (double) n;

/*allocate memory for  variable x and y */
   allocate_xy();


/* write input curve to output file */
   fprintf(output1, "k = %d\n", 0);
   fprintf(output1, "m = %d\n", m);
   fprintf(output1, "n = %d\n", n);
   fprintf(output2, "#k = %d\n", 0);
   fprintf(output2, "#m = %d\n", m);
   fprintf(output2, "#n = %d\n", n);
   for( i=0; i<m; i++ )
   {
     if ( cs[i]->type == 1 )
     {
        fprintf(output1, "curve #%d, Ellipse\n", i+1);
        fprintf(output2, "#curve #%d, Ellipse\n", i+1);

        fprintf(output1, "c1 = %lf\n", cs[i]->c1);
        fprintf(output2, "#c1 = %lf\n",cs[i]->c1);

        fprintf(output1, "c2 = %lf\n", cs[i]->c2);
        fprintf(output2, "#c2 = %lf\n",cs[i]->c2);

        fprintf(output1, "x = %lf\n", cs[i]->x0);
        fprintf(output2, "#x = %lf\n",cs[i]->x0);

        fprintf(output1, "y = %lf\n", cs[i]->y0);
        fprintf(output2, "#y = %lf\n",cs[i]->y0);

        fprintf(output1, "[\n");
        initellip(i); 
        fprintf(output1, "]\n");
        fprintf(output2, "\n");
     }
     else if (cs[i]->type == 2)
     {
        fprintf(output1, "curve #%d, Ellipse\n", i+1);
        fprintf(output2, "#curve #%d, Ellipse\n", i+1);

        fprintf(output1, "c1 = %lf\n", cs[i]->c1);
        fprintf(output2, "#c1 = %lf\n",cs[i]->c1);

        fprintf(output1, "c2 = %lf\n", cs[i]->c2);
        fprintf(output2, "#c2 = %lf\n",cs[i]->c2);

        fprintf(output1, "x = %lf\n", cs[i]->x0);
        fprintf(output2, "#x = %lf\n",cs[i]->x0);

        fprintf(output1, "y = %lf\n", cs[i]->y0);
        fprintf(output2, "#y = %lf\n",cs[i]->y0);

        fprintf(output1, "[\n");
        initellip2(i); 
        fprintf(output1, "]\n");
        fprintf(output2, "\n");
     }
     else if (cs[i]->type == 3)
     {
        fprintf(output1, "curve #%d, Ellipse\n", i+1);
        fprintf(output2, "#curve #%d, Ellipse\n", i+1);

        fprintf(output1, "c1 = %lf\n", cs[i]->c1);
        fprintf(output2, "#c1 = %lf\n",cs[i]->c1);

        fprintf(output1, "c2 = %lf\n", cs[i]->c2);
        fprintf(output2, "#c2 = %lf\n",cs[i]->c2);

        fprintf(output1, "x = %lf\n", cs[i]->x0);
        fprintf(output2, "#x = %lf\n",cs[i]->x0);

        fprintf(output1, "y = %lf\n", cs[i]->y0);
        fprintf(output2, "#y = %lf\n",cs[i]->y0);

        fprintf(output1, "[\n");
        initellip3(i); 
        fprintf(output1, "]\n");
        fprintf(output2, "\n");
     }
     else if (cs[i]->type == 4)
     {
        fprintf(output1, "curve #%d, Ellipse\n", i+1);
        fprintf(output2, "#curve #%d, Ellipse\n", i+1);

        fprintf(output1, "c1 = %lf\n", cs[i]->c1);
        fprintf(output2, "#c1 = %lf\n",cs[i]->c1);

        fprintf(output1, "c2 = %lf\n", cs[i]->c2);
        fprintf(output2, "#c2 = %lf\n",cs[i]->c2);

        fprintf(output1, "x = %lf\n", cs[i]->x0);
        fprintf(output2, "#x = %lf\n",cs[i]->x0);

        fprintf(output1, "y = %lf\n", cs[i]->y0);
        fprintf(output2, "#y = %lf\n",cs[i]->y0);

        fprintf(output1, "[\n");
        initellip4(i); 
        fprintf(output1, "]\n");
        fprintf(output2, "\n");
     }
     else
     {
        fprintf(output1, "curve #%d, Rose\n", i+1);
        fprintf(output2, "#curve #%d, Rose\n", i+1);

        fprintf(output1, "c1 = %lf\n", cs[i]->c1);
        fprintf(output2, "#c1 = %lf\n",cs[i]->c1);

        fprintf(output1, "c2 = %lf\n", cs[i]->c2);
        fprintf(output2, "#c2 = %lf\n",cs[i]->c2);

        fprintf(output1, "x = %lf\n", cs[i]->x0);
        fprintf(output2, "#x = %lf\n",cs[i]->x0);

        fprintf(output1, "y = %lf\n", cs[i]->y0);
        fprintf(output2, "#y = %lf\n",cs[i]->y0);

        fprintf(output1, "[\n");
        initrose(i);
        fprintf(output1, "]\n");
        fprintf(output2, "\n");
     }
   }
   fclose( output1 );
   fclose( output2 );

/* deallocate memory */
   for(i=0; i<m; i++)
   {
     if( cs[i] ) free( cs[i] );
   }
   if( cs ) free( cs );
   if( x ) free ( ( void * ) x );
   if( y ) free ( ( void * ) y );

}  /* end of main */

/*--------------------------------------------------------------------------
     initellip() - this function calculate initial value of ellipse
---------------------------------------------------------------------------*/
void initellip(int i)
{
  int i1;
  for (i1=i*n+1; i1<=(i+1)*n; i1++)
  {
    x[i1] = cs[i]->c1*cos(delta_theta*(i1-1))+cs[i]->x0;
    y[i1] = cs[i]->c2*sin(delta_theta*(i1-1))+cs[i]->y0;
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
  }
  fprintf(output1, "!%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
  fprintf(output2, "%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);

}

/*--------------------------------------------------------------------------
     initellip2() - this function calculate initial value of ellipse
---------------------------------------------------------------------------*/
void initellip2(int i)
{
  int i1;
  for (i1=i*n+1; i1<=(i+1)*n; i1++)
  {
    x[i1] = cs[i]->c1*cos(delta_theta*(i1-1)-0.35*sin(2*delta_theta*(i1-1)))
            +cs[i]->x0;
    y[i1] = cs[i]->c2*sin(delta_theta*(i1-1)-0.35*sin(2*delta_theta*(i1-1)))
            +cs[i]->y0;
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
  }
  fprintf(output1, "!%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
  fprintf(output2, "%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
}


/*--------------------------------------------------------------------------

     initellip4() - this function calculate initial value of dumbbell 
---------------------------------------------------------------------------*/
void initellip4(int i)
{
  int i1,i2,nn;
  double b,b1,r,theta;
    b= cs[i]->c1;
    r= cs[i]->c2;
    b1= b+r*cos(1.0*PI/24.0);
/*    
   b1= b+r*sqrt(3.0)/2.0;
*/
   i2= i*n;
   nn=n/6;
  for (i1=i2+1; i1<=i2+2*nn; i1++)
   {
    theta = -23.0*PI/24.0 + (double)(i1-i2)*23.0/24.0*PI/(double)(nn) ;
/*
    theta = -5.0*PI/6.0 + (double)(i1-i2)*5.0/6.0*PI/(double)(nn) ;
*/
    x[i1] = b1+ r*cos(theta );
    y[i1] =  r*sin(theta); 
    printf("i1= %d\n", i1);
    printf(" %20.15f      %20.15f      \n", x[i1], y[i1], theta);
    fprintf(output1, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
    fprintf(output2, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
   }

     i2=i*n+2*nn;
  for (i1=i2+1;i1<=i2+nn; i1++)
   { 
    x[i1] = b-(double)(i1-i2)*b*2.0/(double)nn;
/*
    y[i1] = 0.5*r;
*/
    y[i1] = r*sin(1.0*PI/24.0);
    fprintf(output1, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
    fprintf(output2, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
   }

   i2= i*n+3*nn;
  for (i1=i2+1; i1<=i2+2*nn; i1++)
   {
    theta = PI/24.0 + (double)(i1-i2)*23.0/24.0*PI/(double)(nn) ;
/*
    theta = PI/6.0 + (double)(i1-i2)*5.0/6.0*PI/(double)(nn) ;
*/
    x[i1] =-b1+ r*cos(theta );
    y[i1] =  r*sin(theta); 
    printf("i1= %d\n", i1);
    printf(" %20.15f      %20.15f      \n", x[i1], y[i1], theta);
    fprintf(output1, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
    fprintf(output2, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
   }

     i2=i*n+5*nn;
  for (i1=i2+1;i1<=i2+nn; i1++)
   { 
    x[i1] = -b+(double)(i1-i2)*b*2.0/(double)nn;
/*
    y[i1] =-0.5*r;
*/
    y[i1] = -r*sin(1.0*PI/24.0);

    fprintf(output1, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
    fprintf(output2, " %20.15f      %20.15f      %20.15f\n", x[i1], y[i1], theta);
  }

  fprintf(output1, "!%20.15f      %20.15f      %20.15f\n", x[i*n+1], y[i*n+1], theta);
  fprintf(output2, "%20.15f      %20.15f      %20.15f\n", x[i*n+1], y[i*n+1], theta);
}


/*--------------------------------------------------------------------------
     initellip3() - this function calculate initial value of long circle
---------------------------------------------------------------------------*/
void initellip3(int i)
{
  int i1,i2;
  double b,r;
    b= cs[i]->c1;
    r= cs[i]->c2;

   i2= i*n;
  for (i1=i2+1; i1<=i2+n/5; i1++)
   { 
    x[i1] = b+ r*cos(-PI/2.0+(double)(i1-i2)*5.0*PI/(double)n );
    y[i1] =  r*sin(-PI/2.0 +(double) (i1-i2)*5.0*PI/(double)n); 
    printf("i1= %d\n", i1);
    printf(" %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }

   i2=i*n+n/5;
  for (i1=i2+1;i1<=i2+n/10; i1++)
   { 
    x[i1] = b-(i1-i2)*10.0/(double)n;
    y[i1] = r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }

   i2= i*n+ 3*n/10;
  for (i1=i2+1;i1<=i2+n/10; i1++)
   { 
    x[i1] = b-1.0-(i1-i2)*2.0*(b-1)*10.0/n;
    y[i1] = r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }

    i2 =i1+ 2*n/5;
  for (i1=i2+1; i1<=i2+n/10; i1++)
   { 
    x[i1] = -b+1.0 -(i1-i2)*10.0/n;
    y[i1] =  r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }


   i2= i*n+ n/2;
  for (i1=i2+1;i1<=i2+n/5; i1++)
   { 
    x[i1] =-b+ r*cos(PI/2.0+(double)((i1-i2)*5.0*PI/n) );
    y[i1] =    r*sin(PI/2.0+(double)((i1-i2)*5.0*PI/n) );
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }

    i2 =i1+ 7*n/10;
  for (i1=i2+1; i1<=i2+n/10; i1++)
   { 
    x[i1] = -b +(i1-i2)*10.0/n;
    y[i1] = -r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }


   i2= i*n+ 4*n/5;
  for (i1=i2+1;i1<=i2+n/10; i1++)
   { 
    x[i1] =-b+1.0+(i1-i2)*2.0*(b-1.0)*10.0/n;
    y[i1] =-r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }
  
  i2 =i1+ 9*n/10;
  for (i1=i2+1; i1<=i2+n/10; i1++)
   { 
    x[i1] = b-1.0 +(i1-i2)*10.0/n;
    y[i1] = -r; 
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
   }
  
  fprintf(output1, "!%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
  fprintf(output2, "%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
}


/*--------------------------------------------------------------------------
     initrose() - this function calculate initial value of ellipse
---------------------------------------------------------------------------*/
void initrose(int i)
{
  int i1;
  for (i1=i*n+1; i1<=(i+1)*n; i1++)
  {
    x[i1] = (2+cs[i]->c1*sin(3*delta_theta*(i1-1)))
                 *cos(delta_theta*(i1-1))+cs[i]->x0;
    y[i1] = (2+cs[i]->c2*sin(3*delta_theta*(i1-1)))
                 *sin(delta_theta*(i1-1))+cs[i]->y0;
    fprintf(output1, " %20.15f      %20.15f      \n", x[i1], y[i1]);
    fprintf(output2, " %20.15f      %20.15f      \n", x[i1], y[i1]);
  }
  fprintf(output1, "!%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
  fprintf(output2, "%20.15f      %20.15f      \n", x[i*n+1], y[i*n+1]);
}

/*---------------------------------------------------------------------------
        allocate_xy() - this function allocates memory for x and y
----------------------------------------------------------------------------*/
void allocate_xy()
{
  x = (double *) calloc(mn+1, sizeof( double ) );
  if ( ! x )
  {
    fprintf(stderr, "Out of memory...\n");
    exit ( -1 );
  }

  y = (double *) calloc(mn+1, sizeof( double ) );
  if ( ! y )
  {
    fprintf(stderr, "Out of memory...\n");
    exit ( -1 );
  }
}