fb:porticula NoPaste

Info / Hilfe Liste
Übersicht / Liste Neu
Datei hochladen Suche
Quellcode suchen
Dateidownload
estimate.c

Uploader:	ytwinky
Datum/Zeit:	13.06.2011 16:21:51
#include "stdio.h"

#include "stdlib.h"

#include "math.h"



#define PI 3.141592653589793238462643



double Factorial(int);



int main(int argc,char **argv)

{

    int i,j,n;

    double a,b,h,e;

    double len,s,theta;

    double x,y,xlast,ylast,dt;

    double fact1,fact2;



    printf("Enter axes length of ellipse: ");

    scanf("%lf %lf",&a,&b);

    e = sqrt(1 - b*b / (a*a));

    h = (a-b)*(a-b) / ((a+b)*(a+b));

    fprintf(stderr,"e: %lf\n",e);

   fprintf(stderr,"h: %lf\n",h);



   n = 10;

   for (j=0;j<7;j++) {

      len = 0;

      for (i=0;i<=n;i++) {

         theta = 0.5 * 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;

      }

      fprintf(stderr,"Numerical             : %.12f (%9d segments)\n",

            len,n);

      n *= 10;

   }



    n = 10;

   for (j=0;j<7;j++) {

      len = 0;

        dt = 0.5 * PI / n;

      for (i=0;i<=n;i++) {

            theta = 0.5 * PI * i / (double)n;

            len += sqrt(1 - e*e*sin(theta)*sin(theta)) * dt;

      }

        len *= a;

      fprintf(stderr,"Integral              : %.12f (%9d segments)\n",

         len,n);

      n *= 10;

   }





   if (a == b)

    fprintf(stderr,"Circle circumference  : %.12f\n",a*0.5*PI);



   len = PI * sqrt(2*(a*a + b*b) - 0.5*(a-b)*(a-b));

   fprintf(stderr,"Anonymous             : %.12f\n",0.25*len);



   len = 0.25 * PI * (a+b) * (3 * (1+h/4) + 1/(1-h/4));

   fprintf(stderr,"Hudson                : %.12f\n",0.25*len);



   len = PI * (3 * (a + b) - sqrt((a + 3*b)*(3*a + b)));

   fprintf(stderr,"Ramanujan 1           : %.12f\n",0.25*len);



   len = PI * (a + b) * (1 + 3*h / (10 + sqrt(4 - 3 * h)));

   fprintf(stderr,"Ramanujan 11          : %.12f\n",0.25*len);



   s = log(2.0) / log(PI/2);

   len = 4 * pow(pow(a,s) + pow(b,s),1.0/s);

   fprintf(stderr,"Necat                 : %.12f (s: %.15f)\n",0.25*len,s);



    if (e < 0.99)

        s = (3*PI - 8) / (8 - 2*PI);

    else

        s = log(2.0) / log(2.0 / (4 - PI));

    len = 4 * (a+b) - 2 * (4-PI) * a * b / pow(pow(a,s)/2 + pow(b,s)/2,1.0/s);

   fprintf(stderr,"Cantrell              : %.12f (s: %.15f)\n",0.25*len,s);



    for (n=2;n<70;n+=10) {

        len = 1;

        for (i=1;i<n;i++) {

            fact1 = Factorial(2*i);

            fact2 = Factorial(i);

            s = fact1;

         s /= pow(2.0,(double)i)*fact2;

         s /= pow(2.0,(double)i)*fact2;

         s *= s;

            s *= -pow(e,2.0*i) / (2.0 * i - 1);

            len += s;

        }

        len *= 2*PI;

        len *= a;

        fprintf(stderr,"Exact                 : %.12f (%02d terms)\n",0.25*len,n);

    }



    return(0);

}



/*

   Iterative definition of factorial

*/

double Factorial(int n)

{

   int i;

   double factor = 1;



   for (i=2;i<=n;i++)

      factor *= i;



   return(factor);

}