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

#define DISJOINT 1
#define INSIDE 2
#define OVERLAP 3

#define PI 3.141592653589793

typedef struct {
   double x;
   double y;
} point;

typedef struct {
   point center;
   double radius;
} circle;

void prnpoint ( point P )
{
   printf("(%10.6lf,%10.6lf )", P.x, P.y);
}

void initcircle ( circle *C, double c, double d, double r )
{
   (C -> center).x = c;
   (C -> center).y = d;
   C -> radius = r;
}

void gencircles ( circle *C1, circle *C2, circle *C3 )
{
   (C1->center).x = rand() % 21;
   (C1->center).y = rand() % 21;
   C1->radius = 5 + rand() % 16;

   do (C2->center).x = rand() % 21; while ((C2->center).x == (C1->center).x);
   do (C2->center).y = rand() % 21; while ((C2->center).y == (C1->center).y);
   C2->radius = 5 + rand() % 16;

   while (1) {
      (C3->center).x = rand() % 21;
      if ((C3->center).x == (C1->center).x) continue;
      if ((C3->center).x == (C2->center).x) continue;
      break;
   }
   while (1) {
      (C3->center).y = rand() % 21;
      if ((C3->center).y == (C1->center).y) continue;
      if ((C3->center).y == (C2->center).y) continue;
      break;
   }
   C3->radius = 5 + rand() % 16;

   /* Special cases */
   // initcircle(C1,0,4,20); initcircle(C2,10,20,16); initcircle(C3,16,9,15);  /* 5(a) */
   // initcircle(C1,16,14,18); initcircle(C2,17,1,10); initcircle(C3,19,20,18);  /* 5(b)(i) */
   // initcircle(C1,0,13,15); initcircle(C2,19,2,9); initcircle(C3,20,19,10);  /* 5(b)(ii) */
   // initcircle(C1,13,5,12); initcircle(C2,15,16,14); initcircle(C3,17,20,14);  /* 5(c)(i) */
   // initcircle(C1,0,13,20); initcircle(C2,12,17,11); initcircle(C3,19,19,14);  /* 5(c)(ii) */

   printf("+++ Input data for three circles\n");
   printf("    C1: %2d %2d %2d\n", (int)(C1->center).x, (int)(C1->center).y, (int)(C1->radius));
   printf("    C2: %2d %2d %2d\n", (int)(C2->center).x, (int)(C2->center).y, (int)(C2->radius));
   printf("    C3: %2d %2d %2d\n", (int)(C3->center).x, (int)(C3->center).y, (int)(C3->radius));
}

void drawcircle ( char grid[][64], circle C, char symbol )
{
   int x, y;
   int c, d, r;

   c = (C.center).x; d = (C.center).y; r = C.radius;
   for (x=c-r; x<=c+r; ++x) {
      y = round(sqrt(r*r-(x-c)*(x-c)));
      grid[d+y+20][x+20] = grid[d-y+20][x+20] = symbol;
   }
   for (y=d-r; y<=d+r; ++y) {
      x = round(sqrt(r*r-(y-d)*(y-d)));
      grid[y+20][c+x+20] = grid[y+20][c-x+20] = symbol;
   }
}

void drawcircles ( circle C1, circle C2, circle C3 )
{
   char grid[64][64];
   int x, y;
   FILE *fp;

   for (x=0; x<=60; ++x) {
      for (y=0; y<=60; ++y) grid[y][x] = ' ';
   }

   grid[20][20] = '+';
   for (x=0; x<=60; ++x) if (x != 20) grid[20][x] = '-';
   for (y=0; y<=60; ++y) if (y != 20) grid[y][20] = '|';

   drawcircle(grid,C1,'1');
   drawcircle(grid,C2,'2');
   drawcircle(grid,C3,'3');

   grid[(int)(C1.center).y+20][(int)(C1.center).x+20] = '1';
   grid[(int)(C2.center).y+20][(int)(C2.center).x+20] = '2';
   grid[(int)(C3.center).y+20][(int)(C3.center).x+20] = '3';

   fp = (FILE *)fopen("circle.txt", "w");
   for (y=60; y>=0; --y) {
      for (x=0; x<=60; ++x) fprintf(fp,"%c ", grid[y][x]);
      fprintf(fp,"\n");
   }
   fclose(fp);
}

double getangle ( circle C, point P )
{
   return atan2(P.y-(C.center).y, P.x-(C.center).x);
}

int relpos ( circle C1, circle C2 )
{
   double c1, c2, d1, d2, r1, r2;
   double D, R;

   c1 = (C1.center).x; d1 = (C1.center).y; r1 = C1.radius;
   c2 = (C2.center).x; d2 = (C2.center).y; r2 = C2.radius;
   D = (c1-c2)*(c1-c2) + (d1-d2)*(d1-d2);
   R = (r1+r2)*(r1+r2);
   if (D >= R) return DISJOINT;
   R = (r1-r2)*(r1-r2);
   if (D <= R) return INSIDE;
   return OVERLAP;
}

void prnrelpos ( int pos )
{
   switch (pos) {
      case DISJOINT: printf("DISJOINT"); break;
      case INSIDE: printf("INSIDE"); break;
      case OVERLAP: printf("OVERLAP"); break;
      default: printf("UNKNOWN");
   }
}

double distance ( point P1, point P2 )
{
   double x1, y1, x2, y2;

   x1 = P1.x; y1 = P1.y;
   x2 = P2.x; y2 = P2.y;
   return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
}

double trianglearea ( point P1, point P2, point P3 )
{
   double a, b, c, s;

   a = distance(P2,P3); b = distance(P1,P3); c = distance(P1,P2);
   s = (a + b + c) / 2;
   return sqrt(s*(s-a)*(s-b)*(s-c));
}

double circlearea ( circle C )
{
   double r;

   r = C.radius;
   return PI * r * r;
}

double segmentarea ( circle C, point P1, point P2 )
{
   double c, d, r, x1, y1, x2, y2;
   double theta1, theta2, theta;
   double tarea, sarea;

   c = (C.center).x; d = (C.center).y; r = C.radius;
   x1 = P1.x; y1 = P1.y;
   x2 = P2.x; y2 = P2.y;

   theta1 = atan2(y1-d,x1-c); theta2 = atan2(y2-d,x2-c);
   theta = theta2 - theta1; if (theta < 0) theta += 2 * PI;

   sarea = theta * r * r / 2;
   tarea = trianglearea(C.center,P1,P2);

   return (theta <= PI) ? sarea - tarea : sarea + tarea;
}

void findintpts ( circle C1, circle C2, point *P1, point *P2 )
{
   double c1, d1, c2, d2, r1, r2;
   double u, v, U, V, W, D;
   double x, y;

   c1 = (C1.center).x; d1 = (C1.center).y; r1 = C1.radius;
   c2 = (C2.center).x; d2 = (C2.center).y; r2 = C2.radius;

   u = (c2-c1)/(d1-d2);
   v = ((c1*c1-c2*c2)+(d1*d1-d2*d2)+(r2*r2-r1*r1))/(2*(d1-d2));

   U = u*u+1;
   V = 2*(u*(v-d1)-c1);
   W = c1*c1+(v-d1)*(v-d1)-r1*r1;

   D = V*V-4*U*W;
   if (D < 0) {
      fprintf(stderr, "*** Error in findintpts: Non-intersecting circles\n");
      D = 0;
   }
   x = (-V+sqrt(D))/(2*U); y = u*x+v; P1->x = x; P1->y = y;
   x = (-V-sqrt(D))/(2*U); y = u*x+v; P2->x = x; P2->y = y;
}

double intarea2 ( circle C1, circle C2 )
{
   double c1, r1, c2, r2;
   double A;
   int pos;
   point P1, P2, P;

   pos = relpos(C1,C2);
   if (pos == DISJOINT) return 0;

   c1 = (C1.center).x; r1 = C1.radius;
   c2 = (C2.center).x; r2 = C2.radius;

   if (pos == INSIDE)
      return (r1 <= r2) ? circlearea(C1) : circlearea(C2);

   findintpts(C1,C2,&P1,&P2);
   if (P1.y > P2.y) { P = P1; P1 = P2; P2 = P; }
   if (c1 < c2) {
      A = segmentarea(C1,P1,P2);
      A += segmentarea(C2,P2,P1);
   } else {
      A = segmentarea(C2,P1,P2);
      A += segmentarea(C1,P2,P1);
   }
   return A;
}

int incircle ( circle C, point P )
{
   double c, d, r, x, y;

   c = (C.center).x; d = (C.center).y; r = C.radius;
   x = P.x; y = P.y;

   return (((x-c)*(x-c)+(y-d)*(y-d)) < r*r) ? 1 : 0;
}

double intarea3 ( circle C1, circle C2, circle C3 )
{
   circle C;
   double r1, r2, r3;
   double theta1, theta2, theta3, theta4;
   double A;
   int pos12, pos13, pos23, in5, in6;
   point Q1, Q2, Q3, Q4, Q5, Q6, Q, P1, P2, P3, P4;

   /* First sort the circles by the x-coordinates. We actually require
      C1 to have the smallest x-coordinate. */
   if ((C1.center).x > (C2.center).x) { C = C1; C1 = C2; C2 = C; }
   if ((C2.center).x > (C3.center).x) { C = C2; C2 = C3; C3 = C; }
   if ((C1.center).x > (C2.center).x) { C = C1; C1 = C2; C2 = C; }

   r1 = C1.radius; r2 = C2.radius; r3 = C3.radius;
   pos12 = relpos(C1,C2); pos13 = relpos(C1,C3); pos23 = relpos(C2,C3);

   /* First handle the simple cases */
   if ((pos12==DISJOINT) || (pos13==DISJOINT) || (pos23==DISJOINT)) return 0;

   if (pos12==INSIDE) return (r1 < r2) ? intarea2(C1,C3) : intarea2(C2,C3);
   if (pos13==INSIDE) return (r1 < r3) ? intarea2(C1,C2) : intarea2(C2,C3);
   if (pos23==INSIDE) return (r2 < r3) ? intarea2(C1,C2) : intarea2(C1,C3);

   /* Here each pair of the given circles is overlapping */
   findintpts(C1,C2,&Q1,&Q2);
   if (Q1.y > Q2.y) { Q = Q1; Q1 = Q2; Q2 = Q; }
   theta1 = getangle(C1,Q1); theta2 = getangle(C1,Q2);

   findintpts(C1,C3,&Q3,&Q4);
   if (Q3.y > Q4.y) { Q = Q3; Q3 = Q4; Q4 = Q; }
   theta3 = getangle(C1,Q3); theta4 = getangle(C1,Q4);

   if (theta2 < theta1) theta2 += 2 * PI;
   if (theta4 < theta3) theta4 += 2 * PI;
   if (((theta2-theta3) > 2*PI) || ((theta2-theta4) > 2*PI)) {
      theta3 += 2*PI; theta4 += 2*PI;
   }
   if (((theta4-theta1) > 2*PI) || ((theta4-theta2) > 2*PI)) {
      theta1 += 2*PI; theta2 += 2*PI;
   }

   findintpts(C2,C3,&Q5,&Q6);

   /* Case of partially overlapping arcs */
   if ( ((theta3 < theta1) && (theta1 < theta4) && (theta4 < theta2)) ||
        ((theta1 < theta3) && (theta3 < theta2) && (theta2 < theta4)) ) {
      printf("--- Case 5(a)\n");
      A = 0;
      P3 = (incircle(C1,Q5)) ? Q5 : Q6;
      if (theta3 < theta1) {
         P1 = Q1; P2 = Q4;
         A += segmentarea(C1,P1,P2);
         A += segmentarea(C3,P2,P3);
         A += segmentarea(C2,P3,P1);
      } else {
         P1 = Q3; P2 = Q2;
         A += segmentarea(C1,P1,P2);
         A += segmentarea(C2,P2,P3);
         A += segmentarea(C3,P3,P1);
      }
      A += trianglearea(P1,P2,P3);
      return A;
   }

   /* Case of disjoint arcs */
   if ((theta2 < theta3) || (theta4 < theta1)) {
      in5 = incircle(C1,Q5); in6 = incircle(C1,Q6);
      if ((in5) && (in6)) {
         printf("--- Case 5(b)(i)\n");
         return intarea2(C2,C3);
      } else if ((!in5) && (!in6)) {
         printf("--- Case 5(b)(ii)\n");
         return 0;
      } else {
         printf("--- Case 5(b)(UNKNOWN)\n");
         return -1;
      }
   }

   /* Case of one arc inside other */
   in5 = incircle(C1,Q5); in6 = incircle(C1,Q6);
   if ((!in5) && (!in6)) {
      printf("--- Case 5(c)(i)\n");
      return (theta2 < theta4) ? intarea2(C1,C2) : intarea2(C1,C3);
   } else if ((in5) && (in6)) {
      printf("--- Case 5(c)(ii)\n");
      A = 0;
      if (Q5.y > Q6.y) {
         P3 = Q5; P4 = Q6;
      } else {
         P3 = Q6; P4 = Q5;
      }
      if (theta2 < theta4) {
         P1 = Q1; P2 = Q2;
         A += segmentarea(C1,P1,P2);
         A += segmentarea(C2,P2,P3);
         A += segmentarea(C3,P3,P4);
         A += segmentarea(C2,P4,P1);
      } else {
         P1 = Q3; P2 = Q4;
         A += segmentarea(C1,P1,P2);
         A += segmentarea(C3,P2,P3);
         A += segmentarea(C2,P3,P4);
         A += segmentarea(C3,P4,P1);
      }
      A += trianglearea(P1,P2,P3);
      A += trianglearea(P1,P3,P4);
      return A;
   } else {
      printf("--- Case 5(c)(UNKNOWN)\n");
      return -1;
   }
}

int main ()
{
   circle C1, C2, C3;
   point P1, P2;

   srand((unsigned int)time(NULL));

   gencircles(&C1,&C2,&C3);
   drawcircles(C1,C2,C3);

   printf("\n+++ Relative positions of\n");
   printf("    C1,C2: "); prnrelpos(relpos(C1,C2)); printf("\n");
   printf("    C1,C3: "); prnrelpos(relpos(C1,C3)); printf("\n");
   printf("    C2,C3: "); prnrelpos(relpos(C2,C3)); printf("\n");

   printf("\n+++ Intersection points\n");
   if (relpos(C1,C2) == OVERLAP) {
      findintpts(C1,C2,&P1,&P2);
      printf("    Between C1 and C2: ");
      prnpoint(P1); printf(" and ");
      prnpoint(P2); printf("\n");
   }
   if (relpos(C1,C3) == OVERLAP) {
      findintpts(C1,C3,&P1,&P2);
      printf("    Between C1 and C3: ");
      prnpoint(P1); printf(" and ");
      prnpoint(P2); printf("\n");
   }
   if (relpos(C2,C3) == OVERLAP) {
      findintpts(C2,C3,&P1,&P2);
      printf("    Between C2 and C3: ");
      prnpoint(P1); printf(" and ");
      prnpoint(P2); printf("\n");
   }

   printf("\n+++ Circle areas\n");
   printf("    area(C1) = %11.6lf\n", circlearea(C1));
   printf("    area(C2) = %11.6lf\n", circlearea(C2));
   printf("    area(C3) = %11.6lf\n", circlearea(C3));

   printf("\n+++ Intersection areas of two circles\n");
   printf("    area(C1,C2) = %11.6lf\n", intarea2(C1,C2));
   printf("    area(C1,C3) = %11.6lf\n", intarea2(C1,C3));
   printf("    area(C2,C3) = %11.6lf\n", intarea2(C2,C3));

   printf("\n+++ Intersection area of three circles\n");
   printf("    area(C1,C2,C3) = %11.6lf\n", intarea3(C1,C2,C3));

   exit(0);
}