// divergence_theorem.c  numerically check Gauss Dirvergence Theorem
// first cube, well with unequal sides, normals, n, are thus known
// integral over volume (del dot F(x,y,z)) dx dy dz =
// surface integral (F dot n) ds  summed for entire surface, 6 for cube
//
// cube x=0..1.0  y=0..1.5  z=0..2.0
// x, y, z direction vector <i, j, k>
// two sides each  xy 1.0 by 1.5  n=<0,0,1> at z=2.0 and <0,0,-1> at z=0
//                 xz 1.0 by 2.0  n=<0,1,0> at y=1.5 and <0,-1,0> at y=0 
//                 yz 1.5 by 2.0  n=<1,0,0> at x=1.0 and <-1,0,0) at x=0 
//
// F(x,y,z) = sin(x) i + cos(0.5*y) j + (z^2+1)/7 k 
//
// we could numerically evaluate div F(x,y,z) and use that
// in the numerical volume integration, yet for a first cut
// we use  div F(x,y,z) = cos(x) - 0.5*sin(0.5*y) + 2z/7 
// by computing the Jacobian and evaluating the determinant.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
static void F(double x, double y, double z, double ff[]);
static void divF(double x, double y, double z, double gf[]);

int main(int argc, char* argv[])
{
  double volume, xy0, xy1, xz0, xz1, yz0, yz1, surface;
  double xmin = 0.0;
  double xmax = 1.0;
  double x, dx;
  int nx = 10;
  double ymin = 0.0;
  double ymax = 1.5;
  double y, dy;
  int ny = 10;
  double zmin = 0.0;
  double zmax = 2.0;
  double z, dz;
  int nz = 10;
  int i, j, k, m;
  double gf[3], ff[3]; // x i, y j, z k 

  printf("divergence_theorem.c running,  numerical check\n");
  printf("In a 3D volume, F is a vector function, e.g.\n");
  printf("F(x,y,z) = sin(x) i + cos(0.5*y) j + (z^2+1)/7 k \n");
  printf("The volume chosen is a cube with dimensions:\n");
  printf("xmin=%f, xmax=%f\n", xmin, xmax);
  printf("ymin=%f, ymax=%f\n", ymin, ymax);
  printf("zmin=%f, zmax=%f\n", zmin, zmax);
  printf("volume integral((del dot F)dv)  should equal\n");
  printf("the sum of 6 surface integral((F dot normal))ds \n");
  printf("\n");

  for(m=0; m<3; m++) // three accuracies
  {
    dx = (xmax-xmin)/(double)nx;
    dy = (ymax-ymin)/(double)ny;
    dz = (zmax-zmin)/(double)nz;
    printf("nx=%d, dx=%f, ny=%d, dy=%f, nz=%d, dz=%f\n",
           nx, dx, ny, dy, nz, dz);

    // integral of volume
    volume = 0.0;
    for(i=0; i<nx; i++)
    {
      x = xmin + dx*(double)i;
      for(j=0; j<ny; j++)
      {
        y = ymin + dy*(double)j;
        for(k=0; k<nz; k++)
        {
          z = zmin + dz*(double)k;
          divF(x+dx/2.0, y+dy/2.0, z+dz/2.0, gf); // center of volume
          volume += gf[0]*dx*dy*dz;
        }
      }
    }
    printf("integral grad F(x,y,z) dx dy dz = %f \n", volume);

    // area integrals
    xy0 = 0.0; // XY planes
    xy1 = 0.0;
    for(i=0; i<nx; i++)
    {
      x = xmin + dx*(double)i;
      for(j=0; j<ny; j++)
      {
        y = ymin + dy*(double)j;
        F(x+dx/2.0, y+dy/2.0, zmin, ff);
        xy0 += -ff[2]*dx*dy; // n = <0,0,-1>
        F(x+dx/2.0, y+dy/2.0, zmax, ff); 
        xy1 +=  ff[2]*dx*dy; // n = <0,0,1>
      }
    }
    printf("surface xy0=%f, xy1=%f\n", xy0, xy1);

    xz0 = 0.0; // XZ planes
    xz1 = 0.0;
    for(i=0; i<nx; i++)
    {
      x = xmin + dx*(double)i;
      for(k=0; k<nz; k++)
      {
        z = zmin + dz*(double)k;
        F(x+dx/2.0, ymin, z+dz/2.0, ff);
        xz0 += -ff[1]*dx*dz; // n = <0,-1,0>
        F(x+dx/2.0, ymax, z+dz/2.0, ff); 
       xz1 +=  ff[1]*dx*dz; // n = <0,1,0>
      }
    }
    printf("surface xz0=%f, xz1=%f\n", xz0, xz1);

    yz0 = 0.0; // YZ planes
    yz1 = 0.0;
    for(j=0; j<ny; j++)
    {
      y = ymin + dy*(double)j;
      for(k=0; k<nz; k++)
      {
        z = zmin + dz*(double)k;
        F(xmin, y+dy/2.0, z+dz/2.0, ff);
        yz0 += -ff[0]*dy*dz; // n = <-1,0.0>
        F(xmax, y+dy/2.0, z+dz/2.0, ff); 
        yz1 +=  ff[0]*dy*dz; // n = <1,0,0>
      }
    }
    printf("surface yz0=%f, yz1=%f\n", yz0, yz1);
    surface = xy0 + xy1 + xz0 + xz1 + yz0 + yz1;
    printf("volume=%f, surface=%f, err=%e\n\n",
           volume, surface, volume-surface);
    nx *= 2; 
    ny *= 2; 
    nz *= 2; 
  } // end m loop
  printf("\ndivergence_theorem.c finished\n");
  return 0;
} // end main

static void F(double x, double y, double z, double ff[])
{
  ff[0] = sin(x);
  ff[1] = cos(0.5*y);
  ff[2] = (z*z+1.0)/7.0;
} // end F

static void divF(double x, double y, double z, double gf[])
{
  gf[0] = cos(x) -0.5*sin(0.5*y) + 2.0*z/7.0;
} // end divF

// end divergence_theron,c