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