---

dumbell.c

Go to the documentation of this file.

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

typedef struct {
  double a, b, c, d, e, f;
} poly6_t;

#define POLY(p, x) (p.a + (x)*(p.b + (x)*(p.c + (x)*(p.d + (x)*(p.e + (x)*p.f)))))
#define POLY1(p, x) (p.b + (x)*(2.*p.c + (x)*(3.*p.d + (x)*(4.*p.e + (x)*5.*p.f))))
#define POLY2(p, x) (2.*p.c + (x)*(6.*p.d + (x)*(12.*p.e + (x)*20.*p.f)))


#define GAUSSTOL 1.e-10

static int gauss(double A[6][6], double B[6], double C[6])
{
  int i, j, k, i1, n = 6;
  double max, temp;
  
  for (j = 0; j <= n - 2; j++) {
    i1 = j;
    max = fabs(A[j][j]);
    for (i = j + 1; i < n; i++) {
      if (fabs(A[i][j]) >= max) {
        i1 = i;
        max = fabs(A[i][j]);
      }
    }
    for (k = j; k < n; k++) {
      max = A[j][k]; A[j][k] = A[i1][k]; A[i1][k] = max;
    }
    max = B[j]; B[j] = B[i1]; B[i1] = max;
    for (i = j + 1; i < n; i++) {
      if (fabs(A[j][j]) < GAUSSTOL)
        return 1;
      else
        temp = A[i][j] / A[j][j];
      for (k = j; k < n; k++)
        A[i][k] -= temp * A[j][k];
      B[i] -= temp * B[j];
    }
  }

  if (fabs(A[n-1][n-1]) < GAUSSTOL)
    return 1;
  C[n - 1] = B[n - 1] / A[n - 1][n - 1];
  for (i = n - 2; i >= 0; i--) {
    C[i] = B[i];
    for (k = i + 1; k < n; k++)
      C[i] -= A[i][k] * C[k];
    if (fabs(A[i][i]) < GAUSSTOL)
      return 1;
    C[i] /= A[i][i];
  }
  return 0;
}


void polycircle(double x1, double y1, double dy1, double ddy1,
                double x2, double y2, double dy2, double ddy2,
                poly6_t * poly)
{
  double A[6][6], B[6], C[6];
  B[0] = y1;
  A[0][0] = 1.0; A[0][1] = x1; A[0][2] = x1*x1; A[0][3] = x1*x1*x1; 
  A[0][4] = x1*x1*x1*x1; A[0][5] = x1*x1*x1*x1*x1;

  B[1] = dy1;
  A[1][0] = 0.0; A[1][1] = 1.0; A[1][2] = 2.*x1; A[1][3] = 3.*x1*x1; 
  A[1][4] = 4.*x1*x1*x1; A[1][5] = 5.*x1*x1*x1*x1;

  B[2] = ddy1;
  A[2][0] = 0.0; A[2][1] = 0.0; A[2][2] = 2.; A[2][3] = 6.*x1; 
  A[2][4] = 12.*x1*x1; A[2][5] = 20.*x1*x1*x1;

  B[3] = y2;
  A[3][0] = 1.0; A[3][1] = x2; A[3][2] = x2*x2; A[3][3] = x2*x2*x2; 
  A[3][4] = x2*x2*x2*x2; A[3][5] = x2*x2*x2*x2*x2;

  B[4] = dy2;
  A[4][0] = 0.0; A[4][1] = 1.0; A[4][2] = 2.*x2; A[4][3] = 3.*x2*x2; 
  A[4][4] = 4.*x2*x2*x2; A[4][5] = 5.*x2*x2*x2*x2;

  B[5] = ddy2;
  A[5][0] = 0.0; A[5][1] = 0.0; A[5][2] = 2.; A[5][3] = 6.*x2; 
  A[5][4] = 12.*x2*x2; A[5][5] = 20.*x2*x2*x2;

  if (gauss(A, B, C)) {
    fprintf(stderr, "gauss: no solution\n");
    exit(1);
  }
  poly->a = C[0]; poly->b = C[1]; poly->c = C[2]; poly->d = C[3];
  poly->e = C[4]; poly->f = C[5];
}



int main(int argc, char * argv[])
{
  int i, n = 200, n1;
  double radius = atof(argv[1]);
  double height = atof(argv[2]);
  double width = atof(argv[3]);
  double theta = height/radius, t, x;
  double xc1 = radius + width/2. + 0.5, xc2 = -radius - width/2. + 0.5;
  poly6_t poly;

  for (i = 0; i < n; i++) {
    t = (double)i/(double)(n - 1)*(PI - theta);
    printf("%g %g\n", xc1 + radius*cos(t), radius*sin(t));
  }
  polycircle(xc1 + radius*cos(PI - theta), 
             radius*sin(PI - theta),
             1./tan(theta),
             -20./radius,
             xc2 + radius*cos(theta), 
             radius*sin(theta),
             -1./tan(theta),
             -20./radius,
             &poly);
  n1 = 50;
  for (i = 1; i < n1 - 1; i++) {
    x = xc1 + radius*cos(PI - theta) + (double)i/(double)(n1 - 1)*
      (xc2 - xc1 + 2.*radius*cos(theta));
    printf("%g %g\n", x, POLY(poly, x));
  }
  for (i = 0; i < n; i++) {
    t = theta + (double)i/(double)(n - 1)*(PI - theta);
    printf("%g %g\n", xc2 + radius*cos(t), radius*sin(t));
  }
  return 0;
}

---