/* $Header$ */

#include <assert.h>
#include <stdio.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "mpi.h"
#include <math.h>

#include "carpet.hh"


//#define DEBUG


#define round(x) (x<0?ceil((x)-0.5):floor((x)+0.5))

#define INITIAL_VAL 10000


//TODO:  overlap with finer levels still seems to be incorrect!


extern "C" void SphericalSurface_Setup (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_REAL const pi = 3.1415926535897932384626433832795028841971693993751;
  
  int group;
  cGroup groupinfo;
  cGroupDynamicData groupdata;
  
  int n;
  int ierr;
  
  
  
  if (nsurfaces == 0) return;
  
  
  
  group = CCTK_GroupIndex ("SphericalSurface::sf_radius");
  assert (group>=0);
  ierr = CCTK_GroupData (group, &groupinfo);
  assert (!ierr);
  ierr = CCTK_GroupDynamicData (cctkGH, group, &groupdata);
  assert (!ierr);
  
  
  
  for (n=0; n<nsurfaces; ++n) {
    
    if (!auto_res[n])
    {
      // set resolution according to given parameters
    
      /* internal consistency checks */
      assert (groupdata.dim == 2);
      assert (groupdata.gsh[0] >= ntheta[n]);
      assert (groupdata.gsh[1] >= nphi[n]);
      
      assert (ntheta[n] >= 3*nghoststheta[n] && ntheta[n] <= maxntheta);
      assert (nphi[n] >= 3*nghostsphi[n] && nphi[n] <= maxnphi);
      
      
      
      /* copy parameters into grid functions */
      sf_ntheta[n] = ntheta[n];
      sf_nphi[n] = nphi[n];
      sf_nghoststheta[n] = nghoststheta[n];
      sf_nghostsphi[n] = nghostsphi[n];
    }
    else
    {
      // we already have an approximate delta-spacing calculated in "SetupRes"
      // based on that, we calculate the number of necessary gridpoints...
      
      double theta_range = pi;
      double phi_range = 2*pi;
      
      if (symmetric_x[n])
         phi_range /= 2;
         
      if (symmetric_y[n])
         phi_range /= 2;
      
      if (symmetric_z[n])
         theta_range /= 2;
      
      sf_ntheta[n] = round(theta_range/sf_delta_theta_estimate[n]) + 2*nghoststheta[n];
      sf_nphi[n]   = round(phi_range/sf_delta_phi_estimate[n]) + 2*nghostsphi[n];
      
      // ...make number of theta-gridpoints odd...
      if (sf_ntheta[n] % 2 != 1)
        sf_ntheta[n] += 1;
      
      // ... and make inner phi-gridpoints devidable by 4 since some thorns require that (eg IsolatedHorizon)
      sf_nphi[n]   += 4 - (sf_nphi[n]-2*nghostsphi[n]) % 4 + 2*nghostsphi[n];
      
      // consistency check
      assert (sf_ntheta[n] >= 3*nghoststheta[n] && sf_ntheta[n] <= maxntheta);
      assert (sf_nphi[n] >= 3*nghostsphi[n] && sf_nphi[n] <= maxnphi);
      
      // check, if there are enough maxntheta,maxnphi gridpoints to remove symmetries
      // (some thorns require that)
      if (symmetric_z[n])
      {
         if (2*(sf_ntheta[n]-2*nghoststheta[n])+2*nghoststheta[n] > maxntheta)
	    CCTK_WARN(0, "maxntheta is not big enough to remove z-symmetry (some thorns may require that)!");
      }
      
      if (symmetric_x[n])
      {
         if (symmetric_y[n])
         {
            if (4*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
            {
               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
            }
         }
         else
         {
            if (2*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
            {
               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
            }
         }
      }
      else
      {
         if (symmetric_y[n])
         {
            if (2*(sf_nphi[n]-2*nghostsphi[n])+2*nghostsphi[n] > maxnphi)
            {
               CCTK_WARN(0, "maxnphi is not big enough to remove x-symmetry and y-symmetry (some thorns may require that)!");
            }
         }
      }

      
      sf_nghoststheta[n] = nghoststheta[n];
      sf_nghostsphi[n] = nghostsphi[n];
      
      if (verbose)
      {
         printf("SphericalSurface[%d]: setting sf_ntheta = %d\n", n, sf_ntheta[n]);
         printf("SphericalSurface[%d]: setting sf_nphi   = %d\n", n, sf_nphi[n]);
      }
    }
    
    if (verbose)
    {
       printf("SphericalSurface[%d]: finest refinement-level that fully contains this surface:  sf_maxreflevel = %d\n", n, sf_maxreflevel[n]);
       printf("SphericalSurface[%d]: finest refinement-level that overlaps with this surface :  sf_minreflevel = %d\n", n, sf_minreflevel[n]);
    }
    
    /* coordinates in the theta direction */
    /* avoid_sf_origin_theta = 1 */
    if (symmetric_z[n]) {
      
      /* upper hemisphere: z>=0, theta in (0, pi/2) */
      sf_delta_theta[n] = pi/2 / (sf_ntheta[n] - 2*nghoststheta[n] - 0.5);
      sf_origin_theta[n] = - (nghoststheta[n] - 0.5) * sf_delta_theta[n];
      
    } else {
      
      /* both hemispheres: theta in (0, pi) */
      sf_delta_theta[n] = pi / (sf_ntheta[n] - 2*nghoststheta[n]);
      sf_origin_theta[n] = - (nghoststheta[n] - 0.5) * sf_delta_theta[n];
      
    }
    
    
    
    /* coordinates in the phi direction */
    /* avoid_sf_origin_phi = 0 */
    if (symmetric_x[n]) {
      if (symmetric_y[n]) {
        
        /* one quadrant: x>=0, y>=0, phi in [0, pi/2] */
        assert (sf_nphi[n] - 2*nghostsphi[n] >= 1);
        sf_delta_phi[n] = pi/2 / (sf_nphi[n] - 2*nghostsphi[n] - 1);
        sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
        
      } else {
        
        /* two quadrants: x>=0, phi in [-pi/2, pi/2] */
        assert (sf_nphi[n] - 2*nghostsphi[n] >= 2);
        sf_delta_phi[n] = pi / (sf_nphi[n] - 2*nghostsphi[n] - 1);
        sf_origin_phi[n] = - pi/2 - nghostsphi[n] * sf_delta_phi[n];
        
      }
    } else {
      if (symmetric_y[n]) {
        
        /* two quadrants: y>=0, phi in [0, pi] */
        assert (sf_nphi[n] - 2*nghostsphi[n] >= 2);
        sf_delta_phi[n] = pi / (sf_nphi[n] - 2*nghostsphi[n] - 1);
        sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
        
      } else {
        
        /* all quadrants: phi in [0, 2pi) */
        assert (sf_nphi[n] - 2*nghostsphi[n] >= 4);
        sf_delta_phi[n] = 2*pi / (sf_nphi[n] - 2*nghostsphi[n]);
        sf_origin_phi[n] = - nghostsphi[n] * sf_delta_phi[n];
        
      }
    }
    
#ifdef DEBUG
      printf("SphericalSurface[%d]: sf_delta_theta = %.6f, sf_delta_phi = %.6f\n", n, sf_delta_theta[n], sf_delta_phi[n]);
#endif
    
    /* mark surface as uninitialised */
    sf_active[n] = 0;
    sf_valid[n] = 0;
    
  } /* for n */
}



extern "C" void SphericalSurface_SetupRes (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  static int first_call = 1;
  int n;
  
  //printf("test\n");
  
  for (n=0; n<nsurfaces; ++n) 
  {
    
    // try to find out with which Carpet::reflevel completely contains
    // the current surface and will intersect and and store the
    // result so that other thorns such as the kick-thorn
    // can decide which reflevel to use
     
    // set resolution according to Radius and Cartesian resolution
      
    CCTK_REAL my_radius;
      
    if (!set_spherical[n] && !set_elliptic[n])
    {
       // we have a problem...without a radius we cannot estimate a proper resolution
       // Here, we simply assume r=1 if radius[n] == 0.
       // If you want a better estimate of the resolution, set the radius parameter to an initial guess!
       
       if (radius[n] == 0)
          my_radius = 1;
       else
	  my_radius = radius[n];
    }
    else
    {
       my_radius = radius[n];
    }
      
#ifdef DEBUG
    printf("SphericalSurface: myradius = %f\n", my_radius);
#endif
    // get theta/phi delta-spacing by looking at Cartesian grid
      
    // first test, if we have an overlap with a certain reflevel
    // we do this by calculating a sphere based on the corners of the Cartesian local patch
    // on this processor.
    unsigned long l[2];
      
    l[0] = CCTK_GFINDEX3D(cctkGH, 0,             0,             0);
    l[1] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1);
      
    int is_overlapping = 0;
    double cart_radius = 0;    // radius of Cart. grid
    double cart_origin[3];     // origin of Cart. grid
      
    // calculate radius of Cartesian grid
    cart_origin[0] = (x[l[1]] - x[l[0]]) / 2.0;
    cart_origin[1] = (y[l[1]] - y[l[0]]) / 2.0;
    cart_origin[2] = (z[l[1]] - z[l[0]]) / 2.0;
      
    cart_radius = sqrt(cart_origin[0]*cart_origin[0] + cart_origin[1]*cart_origin[1] + cart_origin[2]*cart_origin[2]);
      
#ifdef DEBUG
    printf("SphericalSurface: cart_radius = %f, cart_origin = %f, %f, %f\n", cart_radius, cart_origin[0], cart_origin[1], cart_origin[2]);
#endif

    // set origin of Cartesian grid
    cart_origin[0] += x[l[0]];
    cart_origin[1] += y[l[0]];
    cart_origin[2] += z[l[0]];
      
    // norm of displacement vector between origin of Cart. grid and origin of SphericalSurface
    double dist = sqrt( (cart_origin[0]-origin_x[n]) * (cart_origin[0]-origin_x[n])
                      + (cart_origin[1]-origin_y[n]) * (cart_origin[1]-origin_y[n])
                      + (cart_origin[2]-origin_z[n]) * (cart_origin[2]-origin_z[n]));
      
    // does our spherical shell overlap with the processor's local grid patch? 
    if (my_radius+cart_radius >= dist &&  // Cart. grid and SphericalSurface overlap
        cart_radius+dist > my_radius  &&  // Cart. grid patch is not fully contained in this SphericalSurface
        -cart_radius+dist < my_radius)    // Cart. grid patch is not outside of the spherical shell
       is_overlapping = 1;
         
    // get resolution
    if (is_overlapping)
    {
       if (CCTK_Equals(auto_res_grid[n], "overlap"))
       {
          sf_delta_theta_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                                          + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
                                                                          + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
         
          sf_delta_phi_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                                        + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
                                                                        + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
       }
       
       // store reflevel
       sf_minreflevel[n] = Carpet::reflevel;
    }
    else if (first_call)
    {
       // set to very high value initially (otherwise would be zero
       // and reported as minimum value in following MPI communication)
       sf_delta_theta_estimate[n] = 10000;
       sf_delta_phi_estimate[n] = 10000; 
    }
    
    if (CCTK_Equals(auto_res_grid[n], "overlap"))
    {
      // find minimum value of "sf_delta_phi/theta" over all processors
      double dummy = 0;
      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
      sf_delta_theta_estimate[n] = dummy;
         
      dummy = 0;
      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
      sf_delta_phi_estimate[n] = dummy;
    }
    
    // find maximum value of "sf_minreflevel" over all processors
    int idummy;
    MPI_Allreduce(&sf_minreflevel[n], &idummy, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
    sf_minreflevel[n] = idummy;



    // find reflevel that completely contains this surface
    // FIXME: This algorithm does not take into account multiple
    //        disconnected refinement regions.
    
    
    //safety margin of 6 points around sphere
    CCTK_INT np = 6;
   
    CCTK_REAL max_radius = my_radius + ( sqrt(  np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                              + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1)
                                              + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ) ); 

#ifdef DEBUG						 
    /*printf("SphericalSurface: calculated max_radius = %.12f\n", max_radius);
    printf("SphericalSurface: buffer_radius = %.12f\n", sqrt(  np*np*CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                             + np*np*CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1)
                                                             + np*np*CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2) ));*/
#endif

   
    unsigned long k[8];
      
    k[0] = CCTK_GFINDEX3D(cctkGH, 0,             0,             0);
    k[1] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, 0,             0);
    k[2] = CCTK_GFINDEX3D(cctkGH, 0,             cctk_lsh[1]-1, 0);
    k[3] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, 0);
    k[4] = CCTK_GFINDEX3D(cctkGH, 0,             0,             cctk_lsh[2]-1);
    k[5] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, 0,             cctk_lsh[2]-1);
    k[6] = CCTK_GFINDEX3D(cctkGH, 0,             cctk_lsh[1]-1, cctk_lsh[2]-1);
    k[7] = CCTK_GFINDEX3D(cctkGH, cctk_lsh[0]-1, cctk_lsh[1]-1, cctk_lsh[2]-1);
   
    int i;
    int is_contained = 0;
    for (i=0; i < 8; i++)
    {
       if (max_radius < sqrt(x[k[i]]*x[k[i]] + y[k[i]]*y[k[i]] + z[k[i]]*z[k[i]]))
          is_contained = 1;
    }
   
    int max_rl = 0;
    if (is_contained)
    {
       if (CCTK_Equals(auto_res_grid[n], "fully contained"))
       {
          sf_delta_theta_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                                          + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
                                                                          + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
         
          sf_delta_phi_estimate[n] = auto_res_ratio[n] / my_radius * sqrt(CCTK_DELTA_SPACE(0)*CCTK_DELTA_SPACE(0) 
                                                                        + CCTK_DELTA_SPACE(1)*CCTK_DELTA_SPACE(1) 
                                                                        + CCTK_DELTA_SPACE(2)*CCTK_DELTA_SPACE(2));
       }
    
       max_rl = Carpet::reflevel; //(int) floor(log((double) (cctk_levfac[0])/log(2.0))+0.5)+1; 
    }
    else if (first_call)
    {
       // set to zero initially
       sf_delta_theta_estimate[n] = 0;
       sf_delta_phi_estimate[n] = 0; 
    }
    
    
    // find minimum value of "sf_maxreflevel" over all processors
    MPI_Allreduce(&max_rl, &idummy, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
    if (idummy != 0 || first_call)
	sf_maxreflevel[n] = idummy;
    
    // this condition can become true because of this approximate algorithm
    // which does not take into account disconnected ref. regions
    if (sf_maxreflevel[n] > sf_minreflevel[n])
        sf_maxreflevel[n] = sf_minreflevel[n];   // can still be incorrect!
	
    
    if (CCTK_Equals(auto_res_grid[n], "fully contained"))
    {
      // find maximum value of "sf_delta_phi/theta" over all processors
      double dummy = 0;
      MPI_Allreduce(&sf_delta_theta_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
      sf_delta_theta_estimate[n] = dummy;
         
      dummy = 0;
      MPI_Allreduce(&sf_delta_phi_estimate[n], &dummy, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
      sf_delta_phi_estimate[n] = dummy;
    }
    
#ifdef DEBUG
    printf("SphericalSurface: sf_maxreflevel[%d] = %d,  sf_minreflevel[%d] = %d \n", n, sf_maxreflevel[n], n, sf_minreflevel[n]);
    printf("SphericalSurface: sf_delta_theta[%d] = %.6f, sf_delta_phi[%d] = %.6f,  \n", n, sf_delta_theta_estimate[n], n, sf_delta_phi_estimate[n]);
#endif

  }
  
  first_call = 0;
  
}