spherical.h

Go to the documentation of this file.

/** @file spherical.h
 */
/**
# Spherical coordinates
 
The default radius of the sphere is set to one. */
 
double Radius = 1.;
 
/**
Rather than the classical [mathematical
convention](http://en.wikipedia.org/wiki/Spherical_coordinate_system),
we use the geographic convention: *x* is the longitude within
\f$[-180:180]\f$ degrees and *y* is the latitude within \f$[-90:90]\f$
degrees. \f$\Delta\f$ is the characteristic cell length expressed in the
same units as *Radius*. */
 
macro VARIABLES (Point point = point, int _ig = ig, int _jg = jg, int _kg = kg)
{
  VARIABLES (point, _ig, _jg, _kg);
  x = x < -180. ? x + 360. : x > 180. ? x - 360. : x;
  Delta_x = Delta_y = Delta;
  Delta *= Radius*pi/180.;
}
 
/**
On trees we need to define consistent refinement functions. The
default restriction functions (averages) are already consistent
(i.e. they preserve volume and length integrals).
 
Mathematically we have
\f[
cm = fm.y = \cos(\phi)
\f]
\f[
fm.x = 1
\f]
with \f$\phi\f$ the latitude in radians. This is the definition we use for
the length scale factors *fm*.
 
For the volume scale factor *cm*, more care needs to be taken to
guarantee discrete volume conservation i.e.
\f[
\sum_{children}cm_{child} = 4cm_{parent}
\f]
This is achieved by computing the exact area of a cell i.e.
\f[
\Delta^2 \int^{\phi + d\phi/2}_{\phi - d \phi/2} \cos\phi d\phi =
\Delta^2  (\sin(\phi + d\phi/2) - \sin(\phi - d\phi/2)) = \Delta^2 cm
\f]
*/
 
#if TREE
static void refine_cm_lonlat (Point point, scalar cm)
{
  double phi = y*pi/180., dphi = Delta/(2.*Radius);
  double sinphi = sin(phi);
  fine(cm,0,0) = fine(cm,1,0) = (sinphi - sin(phi - dphi))/dphi;
  fine(cm,0,1) = fine(cm,1,1) = (sin(phi + dphi) - sinphi)/dphi;
}
 
static void refine_face_x_lonlat (Point point, scalar fm)
{
  if (!is_refined(neighbor(-1)))
    fine(fm,0,0) = fine(fm,0,1) = 1.;
  if (!is_refined(neighbor(1)) && neighbor(1).neighbors)
    fine(fm,2,0) = fine(fm,2,1) = 1.;
  fine(fm,1,0) = fine(fm,1,1) = 1.;
}
 
static void refine_face_y_lonlat (Point point, scalar fm)
{
  double phi = y*pi/180., dphi = Delta/(2.*Radius);
  if (!is_refined(neighbor(0,-1)))
    fine(fm,0,0) = fine(fm,1,0) = cos(phi - dphi);
  if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors)
    fine(fm,0,2) = fine(fm,1,2) = cos(phi + dphi);
  fine(fm,0,1) = fine(fm,1,1) = cos(phi);
}
#endif
 
/** @brief Event: metric (i = 0) */
void event_metric (void) {
 
  /**
  We initialise the scale factors, taking care to first allocate the
  fields if they are still constant. */
 
  if (is_constant(cm)) {
    scalar * l = list_copy (all);
    cm = {0} /* new scalar */;
    free (all);
    all = list_concat ({cm}, l);
    free (l);
  }
  scalar cmv = cm;
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double phi = y*pi/180., dphi = Delta/(2.*Radius);
    cmv[] = (sin(phi + dphi) - sin(phi - dphi))/(2.*dphi);
  }
 
  if (is_constant(fm.x)) {
    scalar * l = list_copy (all);
    fm = {0} /* new vector */;
    free (all);
    all = list_concat ((scalar *){fm}, l);
    free (l);
  }
  vector fmv = fm;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    fmv.x[] = 1.;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    fmv.y[] = cos(y*pi/180.);
 
  /**
  We set our refinement/prolongation functions on trees. */
 
#if TREE
  cm.refine = cm.prolongation = refine_cm_lonlat;
  fm.x.prolongation = refine_face_x_lonlat;
  fm.y.prolongation = refine_face_y_lonlat;
#endif
}