/* This is "rod6.c", written to accompany Introduction to
   Scientific Programming by Joseph L. Zachary. */

/* This program solves the rod-stacking problem described in
   Chapter 13 of the text. */


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



/* Returns the square of "x". */

double sq (double x)
{
  return(x*x);
}



/* Given the lengths of three sides of a triangle, returns
   the cosine of the angle between the first two. */

double cosine (double adj1, double adj2, double opp)
{
  return((sq(adj1) + sq(adj2) - sq(opp)) / (2*adj1*adj2));
}



/* Returns the distance between the points (x1,y1)
   and (x2,y2). */

double distance (double x1, double y1,
                 double x2, double y2)
{
  return(sqrt(sq(x1-x2) + sq(y1-y2)));
}



/* This function returns the y-coordinate of a circular
   rod of radius "rTop" that is suported by two other rods.
   The left supporting rod has radius "rLeft" and its center
   is at (xLeft,yLeft).  The right supporting rod has
   radius "rRight" and its center is at (xRight,yRight). */

double threeRodY (double xLeft, double yLeft,
                  double xRight, double yRight,
                  double rLeft, double rRight, double rTop)
{
  double a, b, c, d, e;   /* Sides of the two triangles from */
                          /* Figure 13.3 of the book. */

  double cosA, sinA,      /* Cosine and sine of angle alpha */
                          /* from Figure 13.3 of the book. */

         cosB, sinB;      /* Cosine and sine of angle beta */
                          /* from Figure 13.3 of the book. */


  /* Calculate the lengths of the sides of the triangles. */

  a = rLeft + rTop;
  b = rRight + rTop;
  c = distance(xLeft, yLeft, xRight, yRight);
  d = xRight - xLeft;
  e = yRight - yLeft;


  /* Calculate the cosines and sines of alpha and beta. */

  cosA = cosine(a, c, b);
  sinA = sqrt(1 - sq(cosA));
  cosB = d/c;
  sinB = e/c;


  /* Compute and return the y-coordinate of the supported rod. */

  return(yLeft + a * (sinA*cosB + sinB*cosA));

}



/* This function returns the x-coordinate of a circular
   rod of radius "rTop" that is suported by two other rods.
   The left supporting rod has radius "rLeft" and its center
   is at (xLeft,yLeft).  The right supporting rod has
   radius "rRight" and its center is at (xRight,yRight). */

double threeRodX (double xLeft, double yLeft,
                  double xRight, double yRight,
                  double rLeft, double rRight, double rTop)
{
  double a, b, c, d, e;   /* Sides of the two triangles from */
                          /* Figure 13.3 of the book. */
  
  double cosA, sinA,      /* Cosine and sine of angle alpha */
                          /* from Figure 13.3 of the book. */

         cosB, sinB;      /* Cosine and sine of angle alpha */
                          /* from Figure 13.3 of the book. */


  /* Calculate the lengths of the sides of the triangles. */

  a = rLeft + rTop;
  b = rRight + rTop;
  c = distance(xLeft, yLeft, xRight, yRight);
  d = xRight - xLeft;
  e = yRight - yLeft;


  /* Calculate the cosines and sines of alpha and beta. */
  
  cosA = cosine(a, c, b);
  sinA = sqrt(1 - sq(cosA));
  cosB = d/c;
  sinB = e/c;


  /* Compute and return the x-coordinate of the supported rod. */

  return(xLeft + a * (cosA*cosB - sinA*sinB));

}



/* Solves the rod-stacking problem described in Chapter 13 of the
   book.  The known center and coordinates are built into this
   function.  It calculates and displays the center coordinates
   of the five supported rods. */

void main (void)
{

  double x1, x2, x3, x4, x5, x6, x7, x8;   /* Center x-coordinates. */
  double y1, y2, y3, y4, y5, y6, y7, y8;   /* Center y-coordinates. */
  double r1, r2, r3, r4, r5, r6, r7, r8;   /* Radii. */


  /* Set up the known values. */

  x1 = 3; y1 = 3; r1 = 3;
  x2 = 20; y2 = 4; r2 = 4;
  x3 = 38; y3 = 2; r3 = 2;
  r4 = 9; r5 = 8; r6 = 3; r7 = 5; r8 = 4;


  /* Calculate the center coordinates of the supported rods. */

  x4 = threeRodX(x1, y1, x2, y2, r1, r2, r4);
  y4 = threeRodY(x1, y1, x2, y2, r1, r2, r4);

  x5 = threeRodX(x2, y2, x3, y3, r2, r3, r5);
  y5 = threeRodY(x2, y2, x3, y3, r2, r3, r5);

  x6 = threeRodX(x4, y4, x5, y5, r4, r5, r6);
  y6 = threeRodY(x4, y4, x5, y5, r4, r5, r6);

  x7 = threeRodX(x6, y6, x5, y5, r6, r5, r7);
  y7 = threeRodY(x6, y6, x5, y5, r6, r5, r7);

  x8 = threeRodX(x4, y4, x7, y7, r4, r7, r8);
  y8 = threeRodY(x4, y4, x7, y7, r4, r7, r8);


  /* Display results. */

  printf("Rod    X       Y\n");
  printf(" 4 %7.2f %7.2f\n", x4, y4);
  printf(" 5 %7.2f %7.2f\n", x5, y5);
  printf(" 6 %7.2f %7.2f\n", x6, y6);
  printf(" 7 %7.2f %7.2f\n", x7, y7);
  printf(" 8 %7.2f %7.2f\n", x8, y8);

}