sphvol.c

Go to the documentation of this file.
00001 
00002 /****************************************************************************
00003 * MODULE:       R-Tree library 
00004 *              
00005 * AUTHOR(S):    Antonin Guttman - original code
00006 *               Daniel Green (green@superliminal.com) - major clean-up
00007 *                               and implementation of bounding spheres
00008 *               
00009 * PURPOSE:      Multidimensional index
00010 *
00011 * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *****************************************************************************/
00017 
00018 /*
00019  *                   SPHERE VOLUME
00020  *                   by Daniel Green
00021  *                   dgreen@superliminal.com
00022  *
00023  * Calculates and prints the volumes of the unit hyperspheres for
00024  * dimensions zero through the given value, or 9 by default.
00025  * Prints in the form of a C array of double called sphere_volumes.
00026  *
00027  * From formule in "Regular Polytopes" by H.S.M Coxeter, the volume
00028  * of a hypersphere of dimension d is:
00029  *        Pi^(d/2) / gamma(d/2 + 1)
00030  * 
00031  * This implementation works by first computing the log of the above
00032  * function and then returning the exp of that value in order to avoid
00033  * instabilities due to the huge values that the real gamma function
00034  * would return.
00035  *
00036  * Multiply the output volumes by R^n to get the volume of an n
00037  * dimensional sphere of radius R.
00038  */
00039 
00040 #include <stdlib.h>
00041 #include <stdio.h>
00042 #include <math.h>
00043 #include <grass/gis.h>
00044 
00045 static void print_volume(int dimension, double volume)
00046 {
00047     fprintf(stdout, "\t%.6f,  /* dimension %3d */\n", volume, dimension);
00048 }
00049 
00050 static double sphere_volume(double dimension)
00051 {
00052     double log_gamma, log_volume;
00053 
00054     log_gamma = gamma(dimension / 2.0 + 1);
00055     log_volume = dimension / 2.0 * log(M_PI) - log_gamma;
00056     return exp(log_volume);
00057 }
00058 
00059 extern int main(int argc, char *argv[])
00060 {
00061     int dim, max_dims = 9;
00062 
00063     if (2 == argc)
00064         max_dims = atoi(argv[1]);
00065 
00066     fprintf(stdout, "static const double sphere_volumes[] = {\n");
00067     for (dim = 0; dim < max_dims + 1; dim++)
00068         print_volume(dim, sphere_volume(dim));
00069     fprintf(stdout, "};\n");
00070     return 0;
00071 }