/* Tests the sphere slice volume functions. */
/* Last edited on 2023-03-31 03:51:33 by stolfi */

/* Must define _GNU_SOURCE in order to get {asprintf} */
#define _GNU_SOURCE
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>

#include <affirm.h>
#include <bool.h>
#include <jsmath.h>

#include <tsp_lib.h>

#define PI M_PI
/* An obscure mathematical constant. */

#define DMAX 6
/* Max sphere dimension. */

#define NX 4
/* Number of positional slices in each half-axis. */

#define NW 6
/* Number of angular slices between 0 and PI/2. */

int main (int argc, char **argv);

void print_header(FILE *wr, int imin, int imax, double zmax);
void print_dashes(FILE *wr, int imin, int imax);
/* Print header and dashes for dimension, total volume, and columns
  ranging from {imin} to {imax} inclusive, where the last column has
  parameter {zmax}. */

int main (int argc, char **argv)
  { 
    int d;
    
    /* Slices by angle: */
    { int imin = -1;
      int imax = NW;
      double wmax = M_PI/2;
      fprintf(stderr, "Volume fraction between lat = 0 and lat = z:\n");
      print_header(stderr, imin, imax, wmax);
      print_dashes(stderr, imin, imax);
      for (d = 1; d <= DMAX; d++)
        { fprintf(stderr, "%3d", d); 
          double v = ball_vol(d);
          fprintf(stderr, "  %8.5f ", v); 
          int i;
          for (i = imin; i <= imax; i++)
            { double w = wmax*((double)i)/((double)imax);
              double f = ball_zone_vol_frac_ang(d, w);
              fprintf(stderr, " %8.5f", f); 
              /* Checking: */
              if (d <= 4)
                { double ff;
                  if (d == 1)
                    { ff = sin(w)/2; }
                  else if (d == 2)
                    { ff = (w + sin(w)*cos(w))/M_PI; }
                  else if (d == 3)
                    { ff = sin(w)*(3 - sin(w)*sin(w))/4; }
                  else if (d == 4)
                    { ff = (w + 2*sin(2*w)/3 + sin(4*w)/12)/M_PI; }
                  else { assert(FALSE); ff = 0.0; }
                  affirm(fabs(ff - f) <= 1.0e-6, "ball_zone_vol_frac_ang error");
                }
            }
          fprintf(stderr, "\n");
        }
      fprintf(stderr, "\n");
    }
    
    /* Slices by position: */
    { int imin = -NX;
      int imax = NX;
      double xmax = 1.0;
      fprintf(stderr, "Volume fraction between pos = -1 and pos = z:\n");
      print_header(stderr, imin, imax, xmax);
      print_dashes(stderr, imin, imax);
      for (d = 0; d <= DMAX; d++)
        { fprintf(stderr, "%3d", d); 
          double v = ball_vol(d);
          fprintf(stderr, "  %8.5f ", v); 
          int i;
          for (i = imin; i <= imax; i++)
            { double u = xmax*((double)i)/((double)imax);
              double f = ball_cap_vol_frac_pos(d, u);
              fprintf(stderr, " %8.5f", f); 
            }
          fprintf(stderr, "\n");
        }
      fprintf(stderr, "\n");
    }
            
    return 0;
  }

void print_header(FILE *wr, int imin, int imax, double zmax)
  {
    fprintf(stderr, "%3s", "dim"); 
    fprintf(stderr, "  %-8s ", "tot.vol."); 
    int i;
    for (i = imin; i <= imax; i++)
      { double z = zmax*((double)i)/((double)imax);
        fprintf(stderr, " %8.5f", z); 
      }
    fprintf(stderr, "\n");
  }

void print_dashes(FILE *wr, int imin, int imax)
  {
    fprintf(stderr, "%3s", "---"); 
    fprintf(stderr, "  %-8s ", "--------"); 
    int i;
    for (i = imin; i <= imax; i++)
      { fprintf(stderr, " %-8s", "--------"); }
    fprintf(stderr, "\n");
  }