// 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>
//	"stdc.h"
//	"config.hh"
//	"../jtutil/util.hh"
//	"../jtutil/array.hh"
//	"../jtutil/linear_map.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 the following 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 a ghost zone 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_interp  objects.
// Every patch has (points to) 4  ghost_zone  objects, one for each of
// the patch's sides.
//
// A  patch_interp  object does the actual interpolation of data from
// within a patch (for filling in data in another patch's ghost zone).
// A  patch_interp  object points to the patch and patch_edge where it
// will be interpolating.
//
// For example, suppose we have two patches p and q with a common
// angular boundary.  Then the desired network of pointers looks like
// this (omitting the  patch_edge  objects for simplicity):
//
// +-----+                                                       +-----+
// |     | <--> p.interpatch_ghost_zone ---> q.patch_interp ---> |     |
// |  p  |                                                       |  q  |
// |     | <--- p.patch_interp <--- q.interpatch_ghost_zone <--> |     |
// +-----+                                                       +-----+
//
// Because of the mutual pointers, we can't easily construct (say) p's
// interpatch_ghost_zone until after q itself has been constructed, and
// vice versa.  Moreover, the  patch_interp::  constructor needs the
// adjacent-side  ghost_zone  objects to already exist, and it needs to
// know the iperp range of the interpolation region, which can only be
// computed from the adjacent-patch  interpatch_ghost_zone  object.
//
// The solution adopted here is to use a 3-phase algorithm, ultimately
// driven by the patch_system constructor:
// * The patch constructors themselves construct the  patch_edge  objects
//   and links them to/from the patches.
// * The patch_system constructor calls the appropriate functions
//	patch::create_mirror_symmetry_ghost_zone()
//	patch::create_periodic_symmetry_ghost_zone()
//	patch::create_interpatch_ghost_zone()
//   to construct the  ghost_zone  objects and link them to/from the
//   patches.
// * The patch_system constructor calls the functions
//	interpatch_ghost_zone::finish_setup()
//   to finish setting up the  interpatch_ghost_zone  objects, construct
//   the other patch's  patch_interp  objects, and finish linking the
//    interpatch_ghost_zone  objects to the  patch_interp  objects.
//

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

//
// 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_interp;
class patch_system;

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

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

	// each patch has a unique 0-origin small-integer patch number,
	// usually denoted  pn
	int patch_number() const { return patch_number_; }

	// each patch has a unique human-readable patch name for debugging etc
	const char* name() const { return name_; }	// typically "+z" etc

	// are we a +[xyz] or -[xyz] patch?
	bool is_plus() const { return is_plus_; }
	// ... values for the  is_plus_in  constructor argument
	static const bool patch_is_plus = true;
	static const bool patch_is_minus = false;

	// are we a (+/-) x or y or z patch?
	// ... 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'

	// are two patches really the same patch?
	// n.b. this does *not* compare any of the gridfn data!
	bool operator==(const patch& other_patch) const
		{ return patch_number() == other_patch.patch_number(); }

	// (rho,sigma,tau) coordinates as singleton coordinate sets
	local_coords::coords_set coords_set_rho() const
		{ return coords_set_rho_; } 
	local_coords::coords_set coords_set_sigma() const
		{ return coords_set_sigma_; }
	local_coords::coords_set coords_set_tau() const
		{ return coords_set_tau_; }

	// {rho,sigma} coordinate set
	local_coords::coords_set coords_set_rho_sigma() const
		{ return coords_set_rho() | coords_set_sigma(); }

	// (rho,sigma) coordinates 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) coordinates *****
	//
public:

	// 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) <--> local (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 fp   rho_of_xyz(fp x, fp y, fp z) const = 0;
	virtual fp sigma_of_xyz(fp x, fp y, fp z) const = 0;

	// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
	//                         with respect to the local coordinate system
	virtual void xyzcos_of_rho_sigma(fp rho, fp sigma,
					 fp& xcos, fp& ycos, fp& zcos)
		const = 0;

	// partial (rho,sigma) / partial (x,y,z)
	virtual fp partial_rho_wrt_x(fp x, fp y, fp z) const = 0;
	virtual fp partial_rho_wrt_y(fp x, fp y, fp z) const = 0;
	virtual fp partial_rho_wrt_z(fp x, fp y, fp z) const = 0;
	virtual fp partial_sigma_wrt_x(fp x, fp y, fp z) const = 0;
	virtual fp partial_sigma_wrt_y(fp x, fp y, fp z) const = 0;
	virtual fp partial_sigma_wrt_z(fp x, fp y, fp z) const = 0;

	// partial^2 (rho,sigma) / partial (xx,xy,xz,yy,yz)
	virtual fp partial2_rho_wrt_xx(fp x, fp y, fp z) const = 0;
	virtual fp partial2_rho_wrt_xy(fp x, fp y, fp z) const = 0;
	virtual fp partial2_rho_wrt_xz(fp x, fp y, fp z) const = 0;
	virtual fp partial2_rho_wrt_yy(fp x, fp y, fp z) const = 0;
	virtual fp partial2_rho_wrt_yz(fp x, fp y, fp z) const = 0;
	virtual fp partial2_rho_wrt_zz(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_xx(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_xy(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_xz(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_yy(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_yz(fp x, fp y, fp z) const = 0;
	virtual fp partial2_sigma_wrt_zz(fp x, fp y, fp z) const = 0;

	// plotting coordinates (dpx,dpy)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - drho = 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_rho_sigma(fp rho, fp sigma) const = 0;
	virtual fp dpy_of_rho_sigma(fp rho, fp sigma) const = 0;


	//
	// ***** patch edges ****
	//
public:
	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_; }
	const patch_edge& minmax_ang_patch_edge(bool want_min, bool want_rho)
		const
		{
		return want_min ? (want_rho ? min_rho_patch_edge()
					    : min_sigma_patch_edge())
				: (want_rho ? max_rho_patch_edge()
					    : max_sigma_patch_edge());
		}

	// find which patch edge is adjacent to neighboring patch q,
	// or error_exit() if it's not actually a neighboring patch
	// ... computation done using only (rho,sigma) coordinate sets
	//     and min/max dang bounds ==> ok to use in setting up ghost zones
	// ... N_overlap_points = number of grid points (grid spacings
	//     in the perpendicular direction) these patches' nominal grids
	//     overlap,
	//     ... if this is nonzero, then these patches must have
	//         the *same* grid spacing in the perpendicular direction
	//     ... e.g. delta_dang = 5, this patch max_dang = 50,
	//         other patch min_dang = 40 ==> N_overlap_points = 3
	//		p   p   p   p   p
	//			q   q   q   q   q
	const patch_edge& edge_adjacent_to_patch(const patch& q,
						 int N_overlap_points = 0)
		const;


	//
	// ***** ghost zones *****
	//
public:
	ghost_zone& min_rho_ghost_zone() const
		{
		assert(min_rho_ghost_zone_ != NULL);
		return *min_rho_ghost_zone_;
		}
	ghost_zone& max_rho_ghost_zone() const
		{
		assert(max_rho_ghost_zone_ != NULL);
		return *max_rho_ghost_zone_;
		}
	ghost_zone& min_sigma_ghost_zone() const
		{
		assert(min_sigma_ghost_zone_ != NULL);
		return *min_sigma_ghost_zone_;
		}
	ghost_zone& max_sigma_ghost_zone() const
		{
		assert(max_sigma_ghost_zone_ != NULL);
		return *max_sigma_ghost_zone_;
		}

	// ... these fns are defined explicitly (unlike for grid:: stuff)
	//     because they're actually used by Jacobian code in ../gr/
	ghost_zone& minmax_rho_ghost_zone(bool want_min)
		const
		{
		return want_min ? min_rho_ghost_zone()
				: max_rho_ghost_zone();
		}
	ghost_zone& minmax_sigma_ghost_zone(bool want_min)
		const
		{
		return want_min ? min_sigma_ghost_zone()
				: max_sigma_ghost_zone();
		}

	ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho)
		const
		{
		return want_rho ? minmax_rho_ghost_zone(want_min)
				: minmax_sigma_ghost_zone(want_min);
		}

	ghost_zone& ghost_zone_on_edge(const patch_edge &e) const;

	// which of the two ghost zones at a specifieid corner,
	// contains a specified point?
	ghost_zone& corner_ghost_zone_containing_point
		(bool rho_is_min, bool sigma_is_min,	// specifies corner
		 int irho, int isigma)			// specifies point
		const;

	// get ghost zone on specified edge,
	// assert() that it is indeed interpatch,
	// static_cast<> to interpatch_ghost_zone ptr
	interpatch_ghost_zone& interpatch_ghost_zone_on_edge
		(const patch_edge &e)
		const;


	//
	// ***** set up ghost zones
	//
public:

	// assert() that this ghost zone hasn't been set up yet,
	// then set it up as mirror-symmetry
	void create_mirror_symmetry_ghost_zone(const patch_edge& edge);

	// assert() that this ghost zone hasn't been set up yet,
	// then set it up as periodic-symmetry
	void create_periodic_symmetry_ghost_zone
		(const patch_edge& my_edge, const patch_edge& symmetry_edge,
		 bool ipar_map_is_plus);

	// assert() that this ghost zone hasn't been set up yet,
	// then set it up as interpatch
	// ... this only sets up ghost zone in skeletal form; use
	//     interpatch_ghost_zone::finish_setup()  to complete
	//     the setup process
	void create_interpatch_ghost_zone
		(const patch_edge& my_edge, const patch_edge& other_edge,
		 int N_overlap_points);

	// assert() that all ghost zones
	// are fully setup
	void assert_all_ghost_zones_fully_setup() const;

private:
	// helper function for setup_*_ghost_zone():
	// assert() that ghost zone pointer on specified edge is NULL
	// (i.e. that we haven't already setup this ghost zone),
	// then assign new value to it
	void set_ghost_zone(const patch_edge& edge, ghost_zone* gzp);


	//
	// ***** constructor, destructor, et al *****
	//
protected:
	// ... used only from derived classes
	// ... doesn't set up ghost zone info, since this depends on
	//     knowing our neighbouring patches, which might not exist yet
	// ... saves a pointer to name_in[], so this should have a
	//     lifetime at least as long as that of this object
	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,
	      local_coords::coords_set coords_set_sigma_in,
	      local_coords::coords_set coords_set_tau_in,
	      const grid_arrays::grid_array_pars& grid_array_pars_in,
	      const grid::grid_pars& grid_pars_in);
public:
	// 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);


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

	// 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_;
	const patch_edge& max_rho_patch_edge_;
	const patch_edge& min_sigma_patch_edge_;
	const patch_edge& max_sigma_patch_edge_;

	// ghost zones
	// ... pointers are set to NULL by ctor,
	//     reset to non-NULL by set_ghost_zone(), which is called by
	//	  create_mirror_symmetry_ghost_zone()
	//	  create_periodic_symmetry_ghost_zone()
	//	  create_interpatch_ghost_zone()
	ghost_zone* min_rho_ghost_zone_;
	ghost_zone* max_rho_ghost_zone_;
	ghost_zone* min_sigma_ghost_zone_;
	ghost_zone* max_sigma_ghost_zone_;
	};

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

//
// 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); }
	fp rho_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::mu_of_yz(y,z); }
	fp sigma_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::nu_of_xz(x,z); }

	// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
	//                         with respect to the local coordinate system
	void xyzcos_of_rho_sigma(fp rho, fp sigma,
				 fp& xcos, fp& ycos, fp& zcos)
		const
		{ local_coords::xyzcos_of_mu_nu(rho,sigma, xcos,ycos,zcos); }

	// partial (rho,sigma) / partial (x,y,z)
	fp partial_rho_wrt_x(fp x, fp y, fp z) const { return 0.0; }
	fp partial_rho_wrt_y(fp x, fp y, fp z) const
		{ return local_coords::partial_mu_wrt_y(y,z); }
	fp partial_rho_wrt_z(fp x, fp y, fp z) const
		{ return local_coords::partial_mu_wrt_z(y,z); }
	fp partial_sigma_wrt_x(fp x, fp y, fp z) const
		{ return local_coords::partial_nu_wrt_x(x,z); }
	fp partial_sigma_wrt_y(fp x, fp y, fp z) const { return 0.0; }
	fp partial_sigma_wrt_z(fp x, fp y, fp z) const
		{ return local_coords::partial_nu_wrt_z(x,z); }

	// partial^2 (rho,sigma) / partial (xx,xy,xz,yy,yz)
	fp partial2_rho_wrt_xx(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_xy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_xz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_yy(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_yy(y,z); }
	fp partial2_rho_wrt_yz(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_yz(y,z); }
	fp partial2_rho_wrt_zz(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_zz(y,z); }
	fp partial2_sigma_wrt_xx(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_xx(x,z); }
	fp partial2_sigma_wrt_xy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_xz(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_xz(x,z); }
	fp partial2_sigma_wrt_yy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_yz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_zz(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_zz(x,z); }

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

	// 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::grid_array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_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); }
	fp rho_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::nu_of_xz(x, z); }
	fp sigma_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::phi_of_xy(x, y); }

	// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
	//                         with respect to the local coordinate system
	void xyzcos_of_rho_sigma(fp rho, fp sigma,
				 fp& xcos, fp& ycos, fp& zcos)
		const
		{ local_coords::xyzcos_of_nu_phi(rho,sigma, xcos,ycos,zcos); }

	// partial (rho,sigma) / partial (x,y,z)
	fp partial_rho_wrt_x(fp x, fp y, fp z) const
		{ return local_coords::partial_nu_wrt_x(x,z); }
	fp partial_rho_wrt_y(fp x, fp y, fp z) const { return 0.0; }
	fp partial_rho_wrt_z(fp x, fp y, fp z) const
		{ return local_coords::partial_nu_wrt_z(x,z); }
	fp partial_sigma_wrt_x(fp x, fp y, fp z) const
		{ return local_coords::partial_phi_wrt_x(x,y); }
	fp partial_sigma_wrt_y(fp x, fp y, fp z) const
		{ return local_coords::partial_phi_wrt_y(x,y); }
	fp partial_sigma_wrt_z(fp x, fp y, fp z) const { return 0.0; }

	// partial^2 (rho,sigma) / partial (xx,xy,xz,yy,yz)
	fp partial2_rho_wrt_xx(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_xx(x,z); }
	fp partial2_rho_wrt_xy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_xz(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_xz(x,z); }
	fp partial2_rho_wrt_yy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_yz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_zz(fp x, fp y, fp z) const
		{ return local_coords::partial2_nu_wrt_zz(x,z); }
	fp partial2_sigma_wrt_xx(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_xx(x,y); }
	fp partial2_sigma_wrt_xy(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_xy(x,y); }
	fp partial2_sigma_wrt_xz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_yy(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_yy(x,y); }
	fp partial2_sigma_wrt_yz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_zz(fp x, fp y, fp z) const { return 0.0; }

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

	// 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::grid_array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_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); }
	fp rho_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::mu_of_yz(y, z); }
	fp sigma_of_xyz(fp x, fp y, fp z) const
		{ return local_coords::phi_of_xy(x, y); }

	// convert (rho,sigma) --> direction cosines (xcos,ycos,zcos)
	//                         with respect to the local coordinate system
	void xyzcos_of_rho_sigma(fp rho, fp sigma,
				 fp& xcos, fp& ycos, fp& zcos)
		const
		{ local_coords::xyzcos_of_mu_phi(rho,sigma, xcos,ycos,zcos); }

	// partial (rho,sigma) / partial (x,y,z)
	fp partial_rho_wrt_x(fp x, fp y, fp z) const { return 0.0; }
	fp partial_rho_wrt_y(fp x, fp y, fp z) const
		{ return local_coords::partial_mu_wrt_y(y,z); }
	fp partial_rho_wrt_z(fp x, fp y, fp z) const
		{ return local_coords::partial_mu_wrt_z(y,z); }
	fp partial_sigma_wrt_x(fp x, fp y, fp z) const
		{ return local_coords::partial_phi_wrt_x(x,y); }
	fp partial_sigma_wrt_y(fp x, fp y, fp z) const
		{ return local_coords::partial_phi_wrt_y(x,y); }
	fp partial_sigma_wrt_z(fp x, fp y, fp z) const { return 0.0; }

	// partial^2 (rho,sigma) / partial (xx,xy,xz,yy,yz)
	fp partial2_rho_wrt_xx(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_xy(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_xz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_rho_wrt_yy(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_yy(y,z); }
	fp partial2_rho_wrt_yz(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_yz(y,z); }
	fp partial2_rho_wrt_zz(fp x, fp y, fp z) const
		{ return local_coords::partial2_mu_wrt_zz(y,z); }
	fp partial2_sigma_wrt_xx(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_xx(x,y); }
	fp partial2_sigma_wrt_xy(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_xy(x,y); }
	fp partial2_sigma_wrt_xz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_yy(fp x, fp y, fp z) const
		{ return local_coords::partial2_phi_wrt_yy(x,y); }
	fp partial2_sigma_wrt_yz(fp x, fp y, fp z) const { return 0.0; }
	fp partial2_sigma_wrt_zz(fp x, fp y, fp z) const { return 0.0; }

	// plotting coordinates (px,py)
	// ... character string describing how (dpx,dpy) are
	//     defined in terms of (mu,nu,phi), eg "90 - drho = 90 - dphi"
	//     (used for labelling output files)
	const char* name_of_dpx() const
		{
		return is_plus() ? "90 - dsigma = 90 - dphi"
				 : "90 + dsigma = 90 + dphi";
		}
	const char* name_of_dpy() const { return "drho = dmu"; }
	// ... (rho,simga) --> (px,py)
	fp dpx_of_rho_sigma(fp rho, fp sigma) const
		{
		const fp dsigma = jtutil::degrees_of_radians(sigma);
		return is_plus() ? 90.0 - dsigma : 90.0 + dsigma;
		}
	fp dpy_of_rho_sigma(fp rho, fp sigma) const
		{ return jtutil::degrees_of_radians(rho); }


	// 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::grid_array_pars& grid_array_pars_in,
		const grid::grid_pars& grid_pars_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);
	};