path: root/src/patch/patch.hh

// patch.hh -- describes a coordinate/grid patch
// $Id$
//
// ***** how patch boundaries are handled *****
// patch - abstract base class to describe a coordinate/grid patch
//
// z_patch - derived class for a +/- z patch
// x_patch - derived class for a +/- x patch
// y_patch - derived class for a +/- y patch
//

//
// prerequisites:
//	<stdio.h>
//	<assert.h>
//	<math.h>
//	"jt/util++.hh"
//	"jt/array.hh"
//	"jt/linear_map.hh"
//	fp.hh
//	coords.hh
//	grid.hh
//	fd_grid.hh
//

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// ***** how patch boundaries are handled *****
//

//
// Basically, we handle patch boundaries using the usual "ghost zone"
// technique, interpolating values from neighboring patches as necessary.
//
// In more detail, we use 5 interrelated types of objects to handle
// patch boundaries:
//
// A  patch_edge  object represents the basic geometry of a min/max
// rho/sigma side of a patch, i.e. it provides which-side-am-I predicates,
// coordinate conversions between (perp,par) and (rho,sigma), etc.
// Every patch has (points to) 4  patch_edge  objects, one for each of
// the patch's sides.
//
// A  ghost_zone  object describes a patch's ghost zone, and knows how
// to fill in gridfns there based on either the patch system's symmetry
// or interpolation from a neighboring patch.  ghost_zone is an abstract
// base class, from which we derive two classes:
// * A  symmetry_ghost_zone  object describes an patch border which
//   is a (discrete) symmetry of spacetime, either mirror-image or
//   periodic.  Such an object knows how to fill in ghost-zone gridfn
//   data from the "other side" of the symmetry.
// * An  interpatch_ghost_zone  object describes a ghost zone which
//   overlaps another patch.  Such an object knows how to get ghost
//   zone gridfn data from the other patch.  More accurately, it gets
//   the data by asking (calling) the appropriate one of the other
//   patch's  patch_frontier objects.
// Every patch has (points to) 4  ghost_zone  objects, one for each of
// the patch's sides.
//
// A  patch_frontier  object represents the part of a patch *inside*
// its nominal grid from which data is interpolated to compute gridfns
// in another patch's ghost zone.  A  patch_frontier  object stores
// information about the interpolator (which interpolator to use, any
// options for it) and any cached data for the interpolator.  Every
// patch has (points to) between 0 and 4  patch_frontier  objects,
// one for each of the patch's sides where there's an adjacent patch
// and thus from which interpolation is done.
//

//
// There are some unobvious complications to how we set up ghost zone
// gridfn data.  For example, notice that (for example) if we have octant
// symmetry, then data near the +z axis with x and y both negative, comes
// from the +z patch with *two* mirror-image symmetry reflections, i.e.
// a single application of the mirror symmetry doesn't suffice.
//
// Because of this, and because we need data in the "corners" of the
// ghost zones, we actually use a 3-phase algorithm to fill in data in
// the ghost zones:
// * We first fill in data at all the non-corner points in all the
//   symmetry borders, by using the symmetries to get this data from
//   patch nominal-grids.
// * We then fill in data in all the interpatch borders in all the
//   patches, by interpatch interpolating as described above.
// * Finally, we fill in data at all the corner points in all the
//   symmetry borders, by using the symmetries to get this data from
//   patch nominal-grids or ghost zones.
//
// For example, consider the +z patch in an octant patch system.
// The following illustration is looking down the z axis; for purposes
// of illuatration, we take the two patch coordinates to be (x,y).
//
//                   #                                                   **
//   <s-y>  <s-y>    #     i+y    i+y    i+y    i+y    i+y    i+y    i+y
//  (-2,7) (-1,7)  (0,7)  (1,7)  (2,7)  (3,7)  (4,7)  (5,7)  (6,7)  (7,7)
//    s-x    s-x     #                                              *i+x
//                   #                                            **
//   <s-y>  <s-y>    #     i+y    i+y    i+y    i+y    i+y    i+y
//  (-2,6) (-1,6)  (0,6)  (1,6)  (2,6)  (3,6)  (4,6)  (5,6)  (6,6)  (7,6)
//    s-x    s-x     #                                       *i+x    i+x
//                   #                                     **
//                   #                                   **
//  (-2,5) (-1,5)   2,5)--(1,5)--(2,5)--(3,5)--(4,5)--(5,5)  (6,5)  (7,5)   
//    s-x    s-x     #                                  |     i+x    i+x
//                   #                                  |
//                   #                                  |
//  (-2,4) (-1,4)  (0,4)  (1,4)  (2,4)  (3,4)  (4,4)  (5,4)  (6,4)  (7,4)
//    s-x    s-x     #                                  |     i+x    i+x
//                   #                                  |
//                   #                                  |
//  (-2,3) (-1,3)  (0,3)  (1,3)  (2,3)  (3,3)  (4,3)  (5,3)  (6,3)  (7,3)
//    s-x    s-x     #                                  |     i+x    i+x
//                   #                                  |
//                   #                                  |
//  (-2,2) (-1,2)  (0,2)  (1,2)  (2,2)  (3,2)  (4,2)  (5,2)  (6,2)  (7,2)
//    s-x    s-x     #                                  |     i+x    i+x
//                   #                                  |
//                   #                                  |
//  (-2,1) (-1,1)  (0,1)  (1,1)  (2,1)  (3,1)  (4,1)  (5,1)  (6,1)  (7,1)
//    s-x    s-x     #                                  |     i+x    i+x
//                   #                                  |
//                   #                                  |
// #(-2,0)#(-1,0)##(0,0)##(1,0)##(2,0)##(3,0)##(4,0)##(5,0)##(6,0)##(7,0)##
//    s-x    s-x     #                                        i+x    i+x
//                   #
//   <s-y>  <s-y>   s-y    s-y    s-y    s-y    s-y    s-y    s-y    s-y 
//  (-2,-1)(-1,-1) (0,-1) (1,-1) (2,-1) (3,-1) (4,-1) (5,-1) (6,-1) (7,1)
//   <s-x>  <s-x>   s-x    s-x    s-x    s-x    s-x    s-x    s-x    s-x
//                   #
//   <s-y>  <s-y>   s-y    s-y    s-y    s-y    s-y    s-y    s-y    s-y 
//  (-2,-2)(-1,-2) (0,-2) (1,-2) (2,-2) (3,-2) (4,-2) (5,-2) (6,-2) (7,0)
//   <s-x>  <s-x>   s-x    s-x    s-x    s-x    s-x    s-x    s-x    s-x
//                   #
//
// For this example,
// * The xz plane and yz plane are marked with ### lines
// * The nominal +z patch is ([0,5],[0,5]), i.e. 0 <= x,y <= 5; its
//   boundary lines are shown with single lines --- and | .
// * The +z patch's ghost zones are
//	-x: ([-2,-1],[-2, 7])
//	+x: ([ 6, 7],[-2, 7])
//	-y: ([-2, 7],[-2,-1])
//	+y: ([-2, 7],[ 6, 7])
// * Points marked with an "s-x" below or an "s-y" above are computed via
//   mirroring across the -x boundary (yz plane) or -y boundary (xz plane)
//   respectively, as described by the +z patch's -x or -y symmetry_ghost_zone
//   object respectively.  This is done in phase 1 of our 3-phase algorithm.
// * The +z patch's frontiers are
//	+x: ([ 3,4],[-2,7])
//	+y: ([-2,7],[ 3,4])
//   Note that in both cases the frontier includes the points computed
//   by symmetry (in phase 1 of our 3-phase algorithm) on the adjacent
//   edges! There are no -x or -y frontiers, since no interpolation is
//   needed across those boundaries of this patch.
// * Ghost-zone points marked with an "i+x" below or an "i+y" below are
//   computed via interpatch interpolation from the appropriate other patch,
//   as described by the +z patch's +x or +y  interpatch_ghost_zone  object
//   respectively.  This is done in phase 2 of our 3-phase algorithm.
// * The diagonal ** line shows the boundary between the +x and +y borders;
//   points there may be interpolated via either of the two possible
//    interpatch_ghost_zone  objects.
//	[At present we require all points in a given ghost zone
//	to be interpolated from the same neighboring patch and
//	 patch_frontier  object, so must arbitrarily choose one
//	of the two neighbors for the diagonal points.  In theory
//	it would be better to take the average of the two neighbors,
//	but in practice this doesn't matter for horizon finding.
//	(For doing time evolution in my multipatch scheme it
//	_does_ matter...)]
// * Points marked with an "<s-x>" below or an "<s-y>" above are computed
//   via mirroring across the -x boundary (yz plane) or -y boundary
//   (xz plane) respectively, as described by the +z patch's -x or -y
//    symmetry_ghost_zone  object respectively.  This is done in phase 3
//   of our 3-phase algorithm.
//

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// patch - abstract base class to describe a generic coordinate/grid patch
//

//
// There are 3 types of patches, z, x, and y.  Each type uses two of
// (mu,nu,phi) as its angular coordinates (rho,sigma); the remaining
// "unused" one of (mu,nu,phi) is tau.
//
//	z patch ==> (rho,sigma) = (mu,nu)    tau = phi
//	x patch ==> (rho,sigma) = (nu,phi)   tau = mu
//	y patch ==> (rho,sigma) = (mu,phi)   tau = nu
//

// forward declarations
class patch_edge;
class ghost_zone;
class symmetry_ghost_zone;
class interpatch_ghost_zone;
class patch_frontier;
class patch_system;

//
// const qualifiers refer to the gridfn values
//
class	patch
	: public fd_grid
	{
public:
	//
	// ***** patch system, type, and coordinate metadata *****
	//

	// to which patch system do we belong?
	patch_system& my_patch_system() const
		{ return my_patch_system_; }

	// patch number in our patch system
	int patch_number() const { return patch_number_; }

	// human-readable patch name for debugging etc
	// ... should be unique among all patches
	const char *name() const { return name_; }	// typically "+z" etc

	// patch type
	// ... are we a +[xyz] or -[xyz] patch?
	bool is_plus() const { return is_plus_; }
	// ... are we (+/-) x or y or z?
	// ... n.b. type is `char' because this is handy for both
	//	    switch() and human-readable printing
	char ctype() const { return ctype_; }		// 'z' or 'x' or 'y'

	// (rho,sigma,tau) coordinates as singleton coordinate sets
	local_coords::coords_set coords_set_rho() const
		{ return coords_set_rho_; } 
	local_coords::coords_set coord_set_sigma() const
		{ return coords_set_sigma_; }
	local_coords::coords_set coord_set_tau() const
		{ return coords_set_tau_; }
	// ... as human-readable character strings
	//     (for labelling output files etc)
	virtual const char *name_of_rho() const = 0;
	virtual const char *name_of_sigma() const = 0;


	//
	// ***** (rho,sigma,tau) coordinate conversions *****
	//

	// convert (rho,sigma) --> tau
	virtual fp tau_of_rho_sigma(fp rho, fp sigma) const = 0;

	// convert (rho,sigma) --> (mu,nu,phi)
	virtual fp mu_of_rho_sigma(fp rho, fp sigma) const = 0;
	virtual fp nu_of_rho_sigma(fp rho, fp sigma) const = 0;
	virtual fp phi_of_rho_sigma(fp rho, fp sigma) const = 0;

	// convert (rho,sigma) <--> usual polar spherical (theta,phi)
	virtual void theta_phi_of_rho_sigma(fp rho, fp sigma,
					    fp& ps_theta, fp& ps_phi)
		const = 0;
	virtual void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi,
					    fp& rho, fp& sigma)
		const = 0;

	// convert (r,rho,sigma) <--> (x,y,z)
	virtual void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma,
					fp& x, fp& y, fp& z)
		const = 0;
	virtual void r_rho_sigma_of_xyz(fp x, fp y, fp z,
					fp& r, fp& rho, fp& sigma)
		const = 0;
	void global_xyz_of_r_rho_sigma(fp r, fp rho, fp sigma,
				       fp& x, fp& y, fp& z)
		const;
	void r_rho_sigma_of_global_xyz(fp x, fp y, fp z,
				       fp& r, fp& rho, fp& sigma)
		const;

	// plotting coordinates (dpx,dpy)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - dphi"
	//     (used for labelling output files)
	virtual const char *name_of_dpx() const = 0;
	virtual const char *name_of_dpy() const = 0;
	// ... (irho,isimga) --> (px,py)
	virtual fp dpx_of_irho_isigma(int irho, int isigma) const = 0;
	virtual fp dpy_of_irho_isigma(int irho, int isigma) const = 0;


	//
	// ***** patch edges ****
	//
	const patch_edge& min_rho_patch_edge() const
		{ return min_rho_patch_edge_; }
	const patch_edge& max_rho_patch_edge() const
		{ return max_rho_patch_edge_; }
	const patch_edge& min_sigma_patch_edge() const
		{ return min_sigma_patch_edge_; }
	const patch_edge& max_sigma_patch_edge() const
		{ return max_sigma_patch_edge_; }


	//
	// ***** ghost zones *****
	//
	const ghost_zone& min_rho_ghost_zone() const
		{ return *min_rho_ghost_zone_; }
	const ghost_zone& max_rho_ghost_zone() const
		{ return *max_rho_ghost_zone_; }
	const ghost_zone& min_sigma_ghost_zone() const
		{ return *min_sigma_ghost_zone_; }
	const ghost_zone& max_sigma_ghost_zone() const
		{ return *max_sigma_ghost_zone_; }


	//
	// ***** set up patch border/frontier subobjects *****
	//

	// set up a mirror-symmetry patch border
	// ... this is a friend of class symmetry_ghost_zone
	void setup_mirror_sym_ghost_zone(const patch_edge& edge_in);

	// set up a periodic-BC patch border
	// ... this is a friend of class symmetry_ghost_zone
	void setup_periodic_BC_ghost_zone
	   (const patch_edge& my_edge_in, const patch_edge& symmetry_edge_in,
	    int my_edge_sample_ipar, int symmetry_edge_sample_ipar,
	    bool ipar_map_sign);	// one of cpm_map::{plus,minus}_map

	// set up an interpatch patch border,
	// but don't link it to the other patch's frontier
	// (which doesn't exist yet)
	// ... this is a friend of class interpatch_ghost_zone
	void setup_interpatch_ghost_zone(const patch_edge& edge_in,
				     const patch_edge& other_edge_in);

	// set up *another* patch's frontier from which we need
	// to interpolate data, and set up our border to link to it
	void setup_other_patch_frontier
		(const patch_edge& my_edge_in,
		 const patch_edge& other_edge_in,
		 const interpolator_pars& interpolator_pars_in);

private:
	// set up a patch frontier
	// ... used (only) by  setup_other_patch_frontier()
	// ... this is a friend of class patch_frontier
	const patch_frontier& setup_patch_frontier
		(const patch_edge& edge_in,
		 const interpolator_pars& interpolator_pars_in);


protected:
	//
	// Note we can't set up the ghost zone and frontier info
	// in the constructor:  This info depends on knowing our
	// neighbouring patches, which may not exist yet.  Instead,
	// we set up this info later with separate setup_*() functions.
	//
	patch(patch_system &my_patch_system_in, int patch_number_in,
	      const char *name_in, bool is_plus_in, char ctype_in,
	      local_coords::coords_set coords_set_rho_in,
				       coords_set_sigma_in,
				       coords_set_tau_in,
	      const grid_arrays::array_pars& grid_array_pars,
	      const grid::grid_pars& grid_pars_in,
	      int N_gridfns_in);

	// destructor must be virtual to allow destruction
	// of derived classes via ptr/ref to this class
	virtual ~patch();


private:
	// we forbid copying and passing by value
	// by declaring the copy constructor and assignment operator
	// private, but never defining them
	patch(const patch& rhs);
	patch& operator=(const patch& rhs);


private:
	//
	// ***** data members *****
	//

	// type/coordinate metadata
	patch_system &my_patch_system_;
	const int patch_number_;
	const char *name_;
	const bool is_plus_;
	const char ctype_;
	const local_coords::coords_set coords_set_rho_,
				       coords_set_sigma_,
				       coords_set_tau_;

	// edges
	const patch_edge& min_rho_patch_edge_,
			& max_rho_patch_edge_,
			& min_sigma_patch_edge_,
			& max_sigma_patch_edge_;

	// ghost zones
	const ghost_zone *min_rho_ghost_zone_,
			 *max_rho_ghost_zone_;
	const ghost_zone *min_sigma_ghost_zone_,
			 *max_sigma_ghost_zone_;

	// frontiers (NULL pointers if no frontier for this border)
	const patch_frontier *min_rho_patch_frontier_,
			     *max_rho_patch_frontier_;
	const patch_frontier *min_sigma_patch_frontier_,
			     *max_sigma_patch_frontier_;
	};

//*****************************************************************************
//*****************************************************************************
//*****************************************************************************

//
// This class describes a +/- z patch.  It doesn't define any new
// functions not already present in  class patch ; it "just" defines
// non-virtual versions of all the pure virtual functions defined there;
//
//	z patch ==> (rho,sigma) = (mu,nu)    tau = phi
//
class	z_patch
	: public patch
	{
public:
	// human-readable names of (rho,sigma)
	const char *name_of_rho() const { return "mu"; }
	const char *name_of_sigma() const { return "nu"; }

	// convert (rho,sigma) --> tau
	fp tau_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::phi_of_mu_nu(rho, sigma); }

	// convert (rho,sigma) --> (mu,nu,phi)
	fp mu_of_rho_sigma(fp rho, fp sigma) const { return rho; }
	fp nu_of_rho_sigma(fp rho, fp sigma) const { return sigma; }
	fp phi_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::phi_of_mu_nu(rho, sigma); }

	// convert (rho,sigma) <--> usual polar spherical (theta,phi)
	void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi)
		const
		{
		local_coords::theta_phi_of_mu_nu(rho, sigma, ps_theta, ps_phi);
		}
	void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma)
		const
		{
		local_coords::mu_nu_of_theta_phi(ps_theta, ps_phi, rho, sigma);
		}

	// convert (r,rho,sigma) <--> (x,y,z)
	void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z)
		const
		{ local_coords::xyz_of_r_mu_nu(r, rho, sigma, x, y, z); }
	void r_rho_sigma_of_xyz(fp x, fp y, fp z, fp& r, fp& rho, fp& sigma)
		const
		{ local_coords::r_mu_nu_of_xyz(x, y, z, r, rho, sigma); }

	// plotting coordinates (px,py)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - dphi"
	//     (used for labelling output files)
	const char *name_of_dpx() const { return "dnu"; }
	const char *name_of_dpy() const
		{ return is_plus() ? "dmu" : "dmu - 180"; }
	// ... (irho,isimga) --> (px,py)
	fp dpx_of_irho_isigma(int irho, int isigma) const
		{ return jtutil::degrees_of_radians(sigma_of_isigma(isigma)); }
	fp dpy_of_irho_isigma(int irho, int sigma) const
		{
		fp dmu = jtutil::degrees_of_radians(rho_of_irho(irho));
		return is_plus() ? dmu : dmu - 180.0;
		}

	// constructor, destructor
	z_patch(patch_system &my_patch_system_in, int patch_number_in,
		const char *name_in, bool is_plus_in,
		const grid_arrays::array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_in,
		int N_gridfns_in);
	~z_patch() { }

private:
        // we forbid copying and passing by value
        // by declaring the copy constructor and assignment operator
        // private, but never defining them
	z_patch(const z_patch& rhs);
	z_patch& operator=(const z_patch& rhs);
	};

//*****************************************************************************

//
// This class describes a +/- x patch.  It doesn't define any new
// functions not already present in  class patch ; it "just" defines
// non-virtual versions of all the pure virtual functions defined there;
//
//	x patch ==> (rho,sigma) = (nu,phi)   tau = mu
//
class	x_patch
	: public patch
	{
public:
	// human-readable names of (rho,sigma)
	const char *name_of_rho() const { return "nu"; }
	const char *name_of_sigma() const { return "phi"; }

	// convert (rho,sigma) --> tau
	fp tau_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::mu_of_nu_phi(rho, sigma); }

	// convert (rho,sigma) --> (mu,nu,phi)
	fp nu_of_rho_sigma(fp rho, fp sigma) const { return rho; }
	fp phi_of_rho_sigma(fp rho, fp sigma) const { return sigma; }
	fp mu_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::mu_of_nu_phi(rho, sigma); }

	// convert (rho,sigma) <--> usual polar spherical (theta,phi)
	void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi)
		const
		{
		local_coords::theta_phi_of_nu_phi(rho, sigma, ps_theta, ps_phi);
		}
	void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma)
		const
		{
		local_coords::nu_phi_of_theta_phi(ps_theta, ps_phi, rho, sigma);
		}

	// convert (r,rho,sigma) <--> (x,y,z)
	void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z)
		const
		{ local_coords::xyz_of_r_nu_phi(r, rho, sigma, x, y, z); }
	void r_rho_sigma_of_xyz(fp x, fp y, fp z, fp& r, fp& rho, fp& sigma)
		const
		{ local_coords::r_nu_phi_of_xyz(x, y, z, r, rho, sigma); }

	// plotting coordinates (px,py)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - dphi"
	//     (used for labelling output files)
	const char *name_of_dpx() const { return "dnu"; }
	const char *name_of_dpy() const { return "dphi"; }
	// ... (irho,isimga) --> (px,py)
	fp dpx_of_rho_sigma(int irho, int isigma) const
		{ return jtutil::degrees_of_radians(rho_of_irho(irho); }
	fp dpy_of_rho_sigma(int irho, int isigma) const { return sigma; }
		{ return jtutil::degrees_of_radians(sigma_of_isigma(isigma); }

	// constructor, destructor
	x_patch(patch_system &my_patch_system_in, int patch_number_in,
		const char *name_in, bool is_plus_in,
		const grid_arrays::array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_in,
		int N_gridfns_in);
	~x_patch() { }

private:
        // we forbid copying and passing by value
        // by declaring the copy constructor and assignment operator
        // private, but never defining them
	x_patch(const x_patch& rhs);
	x_patch& operator=(const x_patch& rhs);
	};

//*****************************************************************************

//
// This class describes a +/- y patch.  It doesn't define any new
// functions not already present in  class patch ; it "just" defines
// non-virtual versions of all the pure virtual functions defined there;
//
//	y patch ==> (rho,sigma) = (mu,phi)   tau = nu
//
class	y_patch
	: public patch
	{
public:
	// human-readable names of (rho,sigma)
	const char *name_of_rho() const { return "mu"; }
	const char *name_of_sigma() const { return "phi"; }

	// convert (rho,sigma) --> tau
	fp tau_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::nu_of_mu_phi(rho, sigma); }

	// convert (rho,sigma) --> (mu,nu,phi)
	fp mu_of_rho_sigma(fp rho, fp sigma) const { return rho; }
	fp phi_of_rho_sigma(fp rho, fp sigma) const { return sigma; }
	fp nu_of_rho_sigma(fp rho, fp sigma) const
		{ return local_coords::nu_of_mu_phi(rho, sigma); }

	// convert (rho,sigma) <--> usual polar spherical (theta,phi)
	void theta_phi_of_rho_sigma(fp rho, fp sigma, fp& ps_theta, fp& ps_phi)
		const
		{
		local_coords::theta_phi_of_mu_phi(rho, sigma, ps_theta, ps_phi);
		}
	void rho_sigma_of_theta_phi(fp ps_theta, fp ps_phi, fp& rho, fp& sigma)
		const
		{
		local_coords::mu_phi_of_theta_phi(ps_theta, ps_phi, rho, sigma);
		}

	// convert (r,rho,sigma) <--> (x,y,z)
	void xyz_of_r_rho_sigma(fp r, fp rho, fp sigma, fp& x, fp& y, fp& z)
		const
		{ local_coords::xyz_of_r_mu_phi(r, rho, sigma, x, y, z); }
	void r_rho_sigma_of_xyz(fp x, fp y, fp z, fp& r, fp& rho, fp& sigma)
		const
		{ local_coords::r_mu_phi_of_xyz(x, y, z, r, rho, sigma); }

	// plotting coordinates (px,py)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - dphi"
	//     (used for labelling output files)
	const char *name_of_dpx() const
		{
		// assume known value for max_dsigma()
		// ==> we can use a constant string
		assert( fuzzy<fp>::EQ(max_dsigma(), 90.0) );
		return "90 - dphi";
		}
	const char *name_of_dpy() const { return "dmu"; }
	// ... (rho,simga) --> (px,py)
	fp dpx_of_irho_isigma(int irho, int isigma) const
		{ jtutil::degrees_of_radians(sigma_of_isigma(isigma); }
	fp dpy_of_irho_isigma(int irho, int isigma) const
		{ return jtutil::degrees_of_radians(rho_of_irho(irho); }


	// constructor, destructor
	y_patch(patch_system &my_patch_system_in, int patch_number_in,
		const char *name_in, bool is_plus_in,
		const grid_arrays::array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_in,
		int N_gridfns_in);
	~y_patch() { }

private:
        // we forbid copying and passing by value
        // by declaring the copy constructor and assignment operator
        // private, but never defining them
	y_patch(const y_patch& rhs);
	y_patch& operator=(const y_patch& rhs);
	};