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 }