#include "stdio.h" #include "stdlib.h" #include "math.h" #ifndef MIN #define MIN(x,y) (x < y ? x : y) #endif #ifndef MAX #define MAX(x,y) (x > y ? x : y) #endif /* The following is an attempt at evaluating various estimates of the circumference of an ellipse. */ double Bessel(double,double, double, int); double Exact(double,double,double,int); void ErrorString(char *,double,double); int main(int argc,char **argv) { long i,j,n; double a,b,h,e,p,d,dt,t; double s,m,tol,len,truelen,theta; double x,y,xlast=0,ylast=0; char errstr[64]; if (argc < 3) { fprintf(stderr,"Usage: %s axis1 axis2\n",argv[0]); exit(-1); } a = atof(argv[1]); b = atof(argv[2]); if (a < b) { a = atof(argv[2]); b = atof(argv[1]); } e = sqrt(1 - b*b / (a*a)); h = (a-b)*(a-b) / ((a+b)*(a+b)); fprintf(stderr,"e : %.10lf\n",e); fprintf(stderr,"h : %.10lf\n",h); fprintf(stderr,"b/a: %.10lf\n",b/a); /* High segment estimate, calculate truelen truelen = 0; n = 100000000L; for (i=0;i<=n;i++) { theta = 0.5 * M_PI * i / (double)n; x = a * cos(theta); y = b * sin(theta); if (i > 0) truelen += sqrt((x-xlast)*(x-xlast) + (y-ylast)*(y-ylast)); xlast = x; ylast = y; } truelen *= 4; */ // Revised courtesy of Charles Karney tol = pow(0.5,27), s = 0, m = 1; x = MAX(a,b); y = MIN(a,b); while (x - y > tol * y) { t = (x + y) / 2; y = sqrt(x * y); x = t; m *= 2; s += m * (x - y) * (x - y); } truelen = M_PI * ((a + b) * (a + b) - s) / (x + y); // Sum segments of one quadrant n = 10; for (j=0;j<8;j++) { len = 0; for (i=0;i<=n;i++) { theta = 0.5 * M_PI * i / (double)n; x = a * cos(theta); y = b * sin(theta); if (i > 0) len += sqrt((x-xlast)*(x-xlast) + (y-ylast)*(y-ylast)); xlast = x; ylast = y; } len *= 4; ErrorString(errstr,len,truelen); fprintf(stderr,"Numerical : %17.14lf %s (%9ld segments)\n",len,errstr,n); n *= 10; } /* Integral n = 10; for (j=0;j<8;j++) { len = 0; dt = 0.5 * M_PI / n; for (i=0;i<=n;i++) { theta = 0.5 * M_PI * i / (double)n; len += sqrt(1 - e*e*sin(theta)*sin(theta)) * dt; } len *= 4*a; ErrorString(errstr,len,truelen); fprintf(stderr,"Integral : %17.14lf %s (%9ld segments)\n",len,errstr,n); n *= 10; } */ // Revised courtesy of Charles Karney n = 10; for (j=0;j<8;j++) { len = 0; dt = 0.5 * M_PI / n; for (i=0;i 1e-3) { sprintf(s," Error: %8.6lf ",err); return; } err *= 1000; if (err > 1e-3) { sprintf(s," Error: %4.3lf x10(-3) ",err); return; } err *= 1000; if (err > 1e-3) { sprintf(s," Error: %4.3lf x10(-6) ",err); return; } err *= 1000; if (err > 1e-3) { sprintf(s," Error: %4.3lf x10(-9) ",err); return; } err *= 1000; sprintf(s," Error: %4.3lf x10(-12)",err); } /* Bessel formula */ double Bessel(double a,double b, double h, int n) { double len=1,sum; int i,j; for (i=1;i0;j--) { if (j > 1) sum *= (2*(j-1)-1); sum /= (2*j); } len += pow(h,(double)i) * sum * sum; } len *= M_PI * (a + b); return(len); } double Exact(double a,double b,double e,int n) { double s,len; double fact=1; int i; len = 1; for (i=1;i