/* $Header$ */

#include <assert.h>

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



void SphericalSurfaceInfo_Setup (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_REAL const pi = 3.1415926535897932384626433832795028841971693993751;
  int n;
  
  
  
  for (n=0; n<nsurfaces; ++n) {
    
    /* internal consistency checks */
    assert (ntheta[n] >= 3*nghoststheta && ntheta[n] <= maxntheta);
    assert (nphi[n] >= 3*nghostsphi && nphi[n] <= maxnphi);
    
    
    
    /* 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 / (ntheta[n] - 2*nghoststheta - 0.5);
      sf_origin_theta[n] = - (nghoststheta - 0.5) * sf_delta_theta[n];
      
    } else {
      
      /* both hemispheres: theta in (0, pi) */
      sf_delta_theta[n] = pi / (ntheta[n] - 2*nghoststheta);
      sf_origin_theta[n] = - (nghoststheta - 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 (nphi[n] - 2*nghostsphi >= 1);
        sf_delta_phi[n] = pi/2 / (nphi[n] - 2*nghostsphi - 1);
        sf_origin_phi[n] = - nghostsphi * sf_delta_phi[n];
        
      } else {
        
        /* two quadrants: x>=0, phi in [-pi/2, pi/2] */
        assert (nphi[n] - 2*nghostsphi >= 2);
        sf_delta_phi[n] = pi / (nphi[n] - 2*nghostsphi - 1);
        sf_origin_phi[n] = - pi/2 - nghostsphi * sf_delta_phi[n];
        
      }
    } else {
      if (symmetric_y[n]) {
        
        /* two quadrants: y>=0, phi in [0, pi] */
        assert (nphi[n] - 2*nghostsphi >= 2);
        sf_delta_phi[n] = pi / (nphi[n] - 2*nghostsphi - 1);
        sf_origin_phi[n] = - nghostsphi * sf_delta_phi[n];
        
      } else {
        
        /* all quadrants: phi in [0, 2pi) */
        assert (nphi[n] - 2*nghostsphi >= 4);
        sf_delta_phi[n] = 2*pi / (nphi[n] - 2*nghostsphi);
        sf_origin_phi[n] = - nghostsphi * sf_delta_phi[n];
        
      }
    }
    
  }
}