// 
//  CoordSys.h
//  Qdos2 - Coordinate System Conversion Template Functions
//
//  Janet Barnes
//  Nov 1999
//

#ifndef QDOS_COORDSYS_H
#define QDOS_COORDSYS_H

#include <math.h>

namespace QDOS {

// Notes:
// 1. Range checking for angles - none at present, but sure there be 
// e.g check that phi in range [0,pi] and if not should a warning be 
// issued or an error thrown? 
// 2. Test for error conditions in atan2??

// include versions with const T& for output types - allows const objects
// with mutable data to be included e.g. matrix subarrays

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ arctan_yx
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- arctan_yx
template<class T1, class T2>
void arctan_yx(const T1& y_in,const T1& x_in, T2& rads_out)
{
// Cover for maths library routine atan2 - returns angle in radians in
// range [0.2pi)
// Requires that class T has a non-member function atan2 defined if T
// is not a double or convertible to a double
// NB atan2(0.,0) returns 0, atan2(y,0.) returns pi/2, atan2(0.,x) returns 0
    rads_out=atan2(y_in,x_in);
    T1 zero;
    zero=0;
    if (y_in < zero)
    {
        const double twopi=8.*std::atan2(1.,1.);
        rads_out+=twopi;
    }
} 

// const version for const objects with mutable data
template<class T1, class T2>
void arctan_yx(const T1& y_in,const T1& x_in, const T2& rads_out)
{
    rads_out=atan2(y_in,x_in);
    T1 zero;
    zero=0;
    if (y_in < zero)
    {
        const double twopi=8.*std::atan2(1.,1.);
        rads_out+=twopi;
    }
}     
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rad_to_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- rad_to_deg
template<class T1, class T2>
void rad_to_deg(const T1& rad_in,T2& deg_out)
{
// Convert radians to degrees
   T1 pi=4*atan((T1)1);       
   deg_out=180*rad_in/pi;
}

// const version
template<class T1, class T2>
void rad_to_deg(const T1& rad_in,const T2& deg_out)
{
// Convert radians to degrees
   T1 pi=4*atan((T1)1);       
   deg_out=180*rad_in/pi;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ deg_to_rad
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- deg_to_rad
template<class T1, class T2>
void deg_to_rad(const T1& deg_in,T2& rad_out)
{
// Convert degrees to radians
   double pi=4*std::atan(1.);       
   rad_out=pi*deg_in/180;
}

// const version
template<class T1, class T2>
void deg_to_rad(const T1& deg_in,const T2& rad_out)
{
// Convert degrees to radians
   double pi=4*std::atan(1.);       
   rad_out=pi*deg_in/180;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xy_to_rp
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xy_to_rp
template <class T1, class T2>
void xy_to_rp( const T1& xy_in, T2& rp_out)
{
// Transform two-dimensional Cartesian, or xy, coordinates to polar coordinates
// ((r),(p)hi)
// phi is a counterclockwise angular displacement in radians from the
// positive x-axis.
// r is the magnitude of the vector
// rp_out is an indexable class where rp_out[0]=(r), rp_out[1]=(p)hi
// Similarly, xy_in[0]=x, xy_in[1]=y
    xy_to_rp(xy_in[0],xy_in[1],rp_out[0],rp_out[1]);
} 

template <class T1,class T2>
void xy_to_rp( const T1& x_in, const T1& y_in, T2& r_out,T2& phi_out)
{
// Transform two-dimensional Cartesian, or xy, coordinates to polar coordinates
// ((r),(p)hi)
// phi is a counterclockwise angular displacement in radians from the
// positive x-axis.
    arctan_yx(y_in,x_in,phi_out); 
   r_out=sqrt(x_in*x_in+y_in*y_in);
} 

// const versions for const objects with mutable data
template <class T1, class T2>
void xy_to_rp( const T1& xy_in, const T2& rp_out)
{
    xy_to_rp(xy_in[0],xy_in[1],rp_out[0],rp_out[1]);
}    

template <class T1, class T2>
void xy_to_rp( const T1& x_in, const T1& y_in, const T2& r_out, const T2& phi_out)
{
    arctan_yx(y_in,x_in,phi_out);
    r_out=sqrt(x_in*x_in+y_in*y_in);
} 
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xy_to_rp_deg
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------ xy_to_rp_deg
template <class T1,class T2>
void xy_to_rp_deg( const T1& xy_in, T2& rp_out)
{
// Transform two-dimensional Cartesian, or xy, coordinates to polar coordinates
// ((r),(p)hi)
// phi is a counterclockwise angular displacement in degrees from the
// positive x-axis.
// r is the magnitude of the vector
// rp_out is an indexable class where rp_out[0]=(r), rp_out[1]=(p)hi
// Similarly, xy_in[0]=x, xy_in[1]=y
    xy_to_rp_deg(xy_in[0],xy_in[1],rp_out[0],rp_out[1]);
}    

template <class T1,class T2>
void xy_to_rp_deg( const T1& x_in, const T1& y_in, T2& r_out,T2& phi_out)
{
// Transform two-dimensional Cartesian, or xy, coordinates to polar coordinates
// ((r),(p)hi)
// phi is a counterclockwise angular displacement in degrees from the
// positive x-axis.
    T2 phi_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    r_out=sqrt(x_in*x_in+y_in*y_in);
}

// const versions
template <class T1,class T2>
void xy_to_rp_deg( const T1& xy_in, const T2& rp_out)
{
    xy_to_rp_deg(xy_in[0],xy_in[1],rp_out[0],rp_out[1]);
}    

template <class T1,class T2>
void xy_to_rp_deg( const T1& x_in, const T1& y_in, const T2& r_out,\
    const T2& phi_out)
{
    T2 phi_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    r_out=sqrt(x_in*x_in+y_in*y_in);
} 
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rp_to_xy
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- rp_to_xy
template <class T1,class T2>
void rp_to_xy( const T1& rp_in, T2& xy_out)
{
// Transform  polar coordinates ((r),(p)hi) to two-dimensional Cartesian,
// or xy, coordinates
// phi is a counterclockwise angular displacement in radians from the
// positive x-axis.
// rp_in is an indexable class where rp_in[0]=(r), rp_in[1]=(p)hi
// Similarly, xy_out[0]=x, xy_out[1]=y

    rp_to_xy(rp_in[0],rp_in[1],xy_out[0],xy_out[1]);
}  

template <class T1,class T2>
void rp_to_xy( const T1& r_in, const T1& phi_in, T2& x_out, T2& y_out)
{
// Transform  polar coordinates ((r),(p)hi) to two-dimensional Cartesian,
// or xy, coordinates
// phi is a counterclockwise angular displacements in radians from the
// positive x-axis.
    x_out=r_in*cos(phi_in);
    y_out=r_in*sin(phi_in);
}

// const versions
template <class T1,class T2>
void rp_to_xy( const T1& rp_in, const T2& xy_out)
{
    rp_to_xy(rp_in[0],rp_in[1],xy_out[0],xy_out[1]);
}  

template <class T1,class T2>
void rp_to_xy( const T1& r_in, const T1& phi_in, const T2& x_out, \
    const T2& y_out)
{
    x_out=r_in*cos(phi_in);
    y_out=r_in*sin(phi_in);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rp_deg_to_xy
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- rp_deg_to_xy
template <class T1,class T2>
void rp_deg_to_xy( const T1& rp_in, T2& xy_out)
{
// Transform  polar coordinates ((r),(p)hi) to two-dimensional Cartesian,
// or xy, coordinates
// phi is a counterclockwise angular displacement in degrees from the
// positive x-axis.
// rp_in is an indexable class where rp_in[0]=(r), rp_in[1]=(p)hi
// Similarly, xy_out[0]=x, xy_out[1]=y

    rp_deg_to_xy(rp_in[0],rp_in[1],xy_out[0],xy_out[1]);
}  

template <class T1,class T2>
void rp_deg_to_xy( const T1& r_in, const T1& phi_in, T2& x_out, T2& y_out)
{
// Transform  polar coordinates ((r),(p)hi) to two-dimensional Cartesian,
// or xy, coordinates
// phi is a counterclockwise angular displacements in degrees from the
// positive x-axis.
    T1 phi_rads;
    deg_to_rad(phi_in,phi_rads);
    x_out=r_in*cos(phi_rads);
    y_out=r_in*sin(phi_rads);
}

// const versions
template <class T1,class T2>
void rp_deg_to_xy( const T1& rp_in, const T2& xy_out)
{
    rp_deg_to_xy(rp_in[0],rp_in[1],xy_out[0],xy_out[1]);
}  

template <class T1,class T2>
void rp_deg_to_xy( const T1& r_in, const T1& phi_in, const T2& x_out, \
    const T2& y_out)
{
    T1 phi_rads;
    deg_to_rad(phi_in,phi_rads);
    x_out=r_in*cos(phi_rads);
    y_out=r_in*sin(phi_rads);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Rpz_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- Rpz_to_xyz
template <class T1,class T2>
void Rpz_to_xyz( const T1& Rpz_in, T2& xyz_out)
{
// Transform cylindrical coordinates ((R),(p)hi,(z)) to Cartesian, or
// xyz, coordinates
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in radians from the 
// positive x-axis, and z is the height above the x-y plane.
// Rpz_in is an indexable class where Rpz_in[0]=(R), Rpz_in[1]=(p)hi
// and Rpz_in[2]=(z).
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    Rpz_to_xyz(Rpz_in[0],Rpz_in[1],Rpz_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void Rpz_to_xyz( const T1& R_in,const T1& phi_in,const T1& z_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform cylindrical coordinates ((R),(p)hi,(z)) to Cartesian, or
// xyz, coordinates
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacements in radians from the 
// positive x-axis, and z is the height above the x-y plane.
  x_out=R_in*cos(phi_in);
  y_out=R_in*sin(phi_in);
  z_out=z_in;
} 

// const versions
template <class T1,class T2>
void Rpz_to_xyz( const T1& Rpz_in, const T2& xyz_out)
{
    Rpz_to_xyz(Rpz_in[0],Rpz_in[1],Rpz_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void Rpz_to_xyz( const T1& R_in, const T1& phi_in, const T1& z_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
  x_out=R_in*cos(phi_in);
  y_out=R_in*sin(phi_in);
  z_out=z_in;
} 
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Rpz_deg_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------- Rpz_deg_to_xyz
template <class T1,class T2>
void Rpz_deg_to_xyz( const T1& Rpz_in, T2& xyz_out)
{
// Transform cylindrical coordinates ((R),(p)hi,(z)) to Cartesian, or
// xyz, coordinates
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in degrees from the 
// positive x-axis, and z is the height above the x-y plane.
// Rpz_in is an indexable class where Rpz_in[0]=(R), Rpz_in[1]=(p)hi
// and Rpz_in[2]=(z).
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    Rpz_deg_to_xyz(Rpz_in[0],Rpz_in[1],Rpz_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void Rpz_deg_to_xyz( const T1& R_in,const T1& phi_in,const T1& z_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform cylindrical coordinates ((R),(p)hi,(z)) to Cartesian, or
// xyz, coordinates
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacements in degrees from the 
// positive x-axis, and z is the height above the x-y plane.
    T1 phi_rads;
    deg_to_rad(phi_in,phi_rads);   
    x_out=R_in*cos(phi_rads);
    y_out=R_in*sin(phi_rads);
    z_out=z_in;
} 

// const versions
template <class T1,class T2>
void Rpz_deg_to_xyz( const T1& Rpz_in, const T2& xyz_out)
{
    Rpz_deg_to_xyz(Rpz_in[0],Rpz_in[1],Rpz_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void Rpz_deg_to_xyz( const T1& R_in,const T1& phi_in,const T1& z_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
    T1 phi_rads;
    deg_to_rad(phi_in,phi_rads);   
    x_out=R_in*cos(phi_rads);
    y_out=R_in*sin(phi_rads);
    z_out=z_in;
} 
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_Rpz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- xyz_to_Rpz
template <class T1,class T2>
void xyz_to_Rpz( const T1& xyz_in, T2& Rpz_out)
{
// Transform Cartesian, or xyz, coordinates to cylindrical coordinates
// ((R),(p)hi,(z))
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in radians from the 
// positive x-axis, and z is the height above the x-y plane.
// Rpz_out is an indexable class where Rpz_out[0]=(R), Rpz_out[1]=(p)hi
// and Rpz_out[2]=(z).
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z
    xyz_to_Rpz(xyz_in[0],xyz_in[1],xyz_in[2],Rpz_out[0],Rpz_out[1],\
        Rpz_out[2]);
}    

template <class T1,class T2>
void xyz_to_Rpz( const T1& x_in,const T1& y_in,const T1& z_in,T2& R_out,\
    T2& phi_out,T2& z_out)
{
// Transform Cartesian, or xyz, coordinates to cylindrical coordinates
// ((R),(p)hi,(z))
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in radians from the 
// positive x-axis, and z is the height above the x-y plane.
    arctan_yx(y_in,x_in,phi_out);
    R_out=sqrt(x_in*x_in+y_in*y_in);
    z_out=z_in;
}

// const versions
template <class T1,class T2>
void xyz_to_Rpz( const T1& xyz_in, const T2& Rpz_out)
{
    xyz_to_Rpz(xyz_in[0],xyz_in[1],xyz_in[2],Rpz_out[0],Rpz_out[1],\
        Rpz_out[2]);
}    

template <class T1,class T2>
void xyz_to_Rpz( const T1& x_in,const T1& y_in,const T1& z_in,\
    const T2& R_out, const T2& phi_out, const T2& z_out)
{
    arctan_yx(y_in,x_in,phi_out);
    R_out=sqrt(x_in*x_in+y_in*y_in);
    z_out=z_in;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_Rpz_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------- xyz_to_Rpz_deg
template <class T1,class T2>
void xyz_to_Rpz_deg( const T1& xyz_in, T2& Rpz_out)
{
// Transform Cartesian, or xyz, coordinates to cylindrical coordinates
// ((R),(p)hi,(z))
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in degrees from the 
// positive x-axis, and z is the height above the x-y plane.
// Rpz_out is an indexable class where Rpz_out[0]=(R), Rpz_out[1]=(p)hi
// and Rpz_out[2]=(z).
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z
    xyz_to_Rpz_deg(xyz_in[0],xyz_in[1],xyz_in[2],Rpz_out[0],Rpz_out[1],\
        Rpz_out[2]);
}    

template <class T1,class T2>
void xyz_to_Rpz_deg( const T1& x_in,const T1& y_in,const T1& z_in,T2& R_out,\
    T2& phi_out,T2& z_out)
{
// Transform Cartesian, or xyz, coordinates to cylindrical coordinates
// ((R),(p)hi,(z))
// R is the distance from the origin to a point in the x-y plane, phi is a
// counterclockwise angular displacement in degrees from the 
// positive x-axis, and z is the height above the x-y plane.
    T2 phi_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    R_out=sqrt(x_in*x_in+y_in*y_in);
    z_out=z_in;
}  

// const versions
template <class T1,class T2>
void xyz_to_Rpz_deg( const T1& xyz_in, const T2& Rpz_out)
{
    xyz_to_Rpz_deg(xyz_in[0],xyz_in[1],xyz_in[2],Rpz_out[0],Rpz_out[1],\
        Rpz_out[2]);
}    

template <class T1,class T2>
void xyz_to_Rpz_deg( const T1& x_in,const T1& y_in,const T1& z_in,\
    const T2& R_out, const T2& phi_out, const T2& z_out)
{
    T2 phi_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    R_out=sqrt(x_in*x_in+y_in*y_in);
    z_out=z_in;
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_to_x
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rtp_to_x
template <class T1,class T2>
void rtp_to_x( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& x_out)
{
// Obtain cartesian x-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,pi), phi in range [0,2pi)
    T2 xy=r_in*sin(theta_in);
    x_out=xy*cos(phi_in);
}

// const version
template<class T1, class T2>
void rtp_to_x( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& x_out)
{
    T2 xy=r_in*sin(theta_in);
    x_out=xy*cos(phi_in);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_to_y
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rtp_to_y
template <class T1,class T2>
void rtp_to_y( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& y_out)
{
// Obtain cartesian y-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,pi), phi in range [0,2pi)
    T2 xy=r_in*sin(theta_in);
    y_out=xy*sin(phi_in);
}

// const version
template<class T1, class T2>
void rtp_to_y( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& y_out)
{
    T2 xy=r_in*sin(theta_in);
    y_out=xy*sin(phi_in);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_to_z
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rtp_to_z
template <class T1,class T2>
void rtp_to_z( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& z_out)
{
// Obtain cartesian z-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,pi), phi in range [0,2pi)
    z_out=r_in*cos(theta_in);
}

// const version
template<class T1,class T2>
void rtp_to_z( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& z_out)
{
   z_out=r_in*cos(theta_in);
}                                  
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- rtp_to_xyz
template <class T1,class T2>
void rtp_to_xyz( const T1& rtp_in, T2& xyz_out)
{
// Transform spherical coordinates ((r),(t)heta,(p)hi) to Cartesian, or
// xyz, coordinates
// Theta and phi are angular displacements in radians from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,pi), phi in range [0,2pi)
// rtp_in is an indexable class where rtp_in[0]=(r), rtp_in[1]=(t)heta
// and rtp_in[2]=(p)hi.
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    rtp_to_xyz(rtp_in[0],rtp_in[1],rtp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void rtp_to_xyz( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform spherical coordinates ((r),(t)heta,(p)hi) to Cartesian, or
// xyz, coordinates
// Theta and phi are angular displacements in radians from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,pi), phi in range [0,2pi)
    T2 xy=r_in*sin(theta_in);
    x_out=xy*cos(phi_in);
    y_out=xy*sin(phi_in);
    z_out=r_in*cos(theta_in);
}

// const versions
template <class T1,class T2>
void rtp_to_xyz( const T1& rtp_in, const T2& xyz_out)
{
    rtp_to_xyz(rtp_in[0],rtp_in[1],rtp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void rtp_to_xyz( const T1& r_in,const T1& theta_in,const T1& phi_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
    T2 xy=r_in*sin(theta_in);
    x_out=xy*cos(phi_in);
    y_out=xy*sin(phi_in);
    z_out=r_in*cos(theta_in);
}

// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_deg_to_x
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ---------------------------------------------------------------- rtp_deg_to_x
template <class T1,class T2>
void rtp_deg_to_x( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& x_out)
{
// Obtain cartesian x-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,180), phi in range [0,360)
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
}

// const version
template<class T1, class T2>
void rtp_deg_to_x( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& x_out)
{
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads); 
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_deg_to_y
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- rtp_deg_to_y
template <class T1,class T2>
void rtp_deg_to_y( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& y_out)
{
// Obtain cartesian y-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,180), phi in range [0,260)
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    y_out=xy*sin(phi_rads);
}

// const version
template<class T1, class T2>
void rtp_deg_to_y( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& y_out)
{
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    y_out=xy*sin(phi_rads);  
}
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_deg_to_z
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// -------------------------------------------------------------- rtp_deg_to_z
template <class T1,class T2>
void rtp_deg_to_z( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& z_out)
{
// Obtain cartesian z-coordinate from spherical coordinates 
// ((r),(t)heta,(p)hi where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis, and r is the magnitude of the radial distance.
// theta in range [0,180), phi in range [0,360)
    T1 theta_rads;
    deg_to_rad(theta_in,theta_rads);   
    z_out=r_in*cos(theta_rads);
}

// const version
template<class T1, class T2>
void rtp_deg_to_z( const T1& r_in,const T1& theta_in,const T1& phi_in,\
const T2& z_out)
{
    T1 theta_rads;
    deg_to_rad(theta_in,theta_rads);   
    z_out=r_in*cos(theta_rads);
}                              
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rtp_deg_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------- rtp_deg_to_xyz
template <class T1,class T2>
void rtp_deg_to_xyz( const T1& rtp_in, T2& xyz_out)
{
// Transform spherical coordinates ((r),(t)heta,(p)hi) to Cartesian, or
// xyz, coordinates
// Theta and phi are angular displacements in degrees from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,180), phi in range [0,360)
// rtp_in is an indexable class where rtp_in[0]=(r), rtp_in[1]=(t)heta
// and rtp_in[2]=(p)hi.
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    rtp_deg_to_xyz(rtp_in[0],rtp_in[1],rtp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void rtp_deg_to_xyz( const T1& r_in,const T1& theta_in,const T1& phi_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform spherical coordinates ((r),(t)heta,(p)hi) to Cartesian, or
// xyz, coordinates
// Theta and phi are angular displacements in degrees from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,180), phi in range [0,360)
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
    y_out=xy*sin(phi_rads);
    z_out=r_in*cos(theta_rads);
}

// const vesrions
template <class T1,class T2>
void rtp_deg_to_xyz( const T1& rtp_in, const T2& xyz_out)
{
    rtp_deg_to_xyz(rtp_in[0],rtp_in[1],rtp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}    

template <class T1,class T2>
void rtp_deg_to_xyz( const T1& r_in,const T1& theta_in,const T1& phi_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
    T1 theta_rads,phi_rads;
    deg_to_rad(theta_in,theta_rads);    
    T2 xy=r_in*sin(theta_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
    y_out=xy*sin(phi_rads);
    z_out=r_in*cos(theta_rads);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_r
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_r
template <class T1,class T2>
void xyz_to_r( const T1& x_in,const T1& y_in,const T1& z_in,T2& r_out)
{
// Obtain magnitude(r) from Cartesian, or xyz, coordinates
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    T1 sqr_z=z_in*z_in;
    r_out=sqrt(sqr_x+sqr_y+sqr_z);
}

// const versions
template <class T1,class T2>
void xyz_to_r( const T1& x_in,const T1& y_in,const T1& z_in, const T2& r_out)
{
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    T1 sqr_z=z_in*z_in;
    r_out=sqrt(sqr_x+sqr_y+sqr_z);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_t
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_t
template <class T1,class T2>
void xyz_to_t( const T1& x_in,const T1& y_in,const T1& z_in,T2& theta_out)
{
// Obtain theta(t) from Cartesian, or xyz, coordinates
// Theta_out is angular displacement in radians from the positive z-axis [0,pi)
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_out);
}

// const versions
template <class T1,class T2>
void xyz_to_t( const T1& x_in,const T1& y_in,const T1& z_in,const T2& theta_out)
{
// Obtain theta(t) from Cartesian, or xyz, coordinates
// Theta_out is angular displacement in radians from the positive z-axis [0,pi)
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_out);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_p
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_p
template <class T1,class T2>
void xyz_to_p( const T1& x_in,const T1& y_in,const T1& z_in,T2& phi_out)
{
// Obtain phi(p) from Cartesian, or xyz, coordinates
// phi_out is angular displacement in radians from the positive x-axis [0,2pi)
   arctan_yx(y_in,x_in,phi_out);   
}

// const versions
template <class T1,class T2>
void xyz_to_p( const T1& x_in,const T1& y_in,const T1& z_in,const T2& phi_out)
{
   arctan_yx(y_in,x_in,phi_out);   
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_rtp
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ------------------------------------------------------------- xyz_to_rtp
template <class T1,class T2>
void xyz_to_rtp( const T1& xyz_in, T2& rtp_out)
{
// Transform Cartesian, or xyz, coordinates to spherical coordinates
// ((r),(t)heta,(p)hi)
// Theta and phi are angular displacements in radians from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,pi), phi in range [0,2pi)
// rtp_out is an indexable class where rtp_out[0]=(r), rtp_out[1]=(t)heta
// and rtp_out[2]=(p)hi.
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z
    xyz_to_rtp(xyz_in[0],xyz_in[1],xyz_in[2],rtp_out[0],rtp_out[1],\
        rtp_out[2]);
}    

template <class T1,class T2>
void xyz_to_rtp( const T1& x_in,const T1& y_in,const T1& z_in,T2& r_out,\
    T2& theta_out,T2& phi_out)
{
// Transform Cartesian, or xyz, coordinates to spherical coordinates
// ((r),(t)heta,(p)hi)
// Theta and phi are angular displacements in radians from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,pi), phi in range [0,2pi)
    arctan_yx(y_in,x_in,phi_out);
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    T1 sqr_z=z_in*z_in;
    r_out=sqrt(sqr_x+sqr_y+sqr_z);
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_out);
}

// const versions
template <class T1,class T2>
void xyz_to_rtp( const T1& xyz_in, const T2& rtp_out)
{
    xyz_to_rtp(xyz_in[0],xyz_in[1],xyz_in[2],rtp_out[0],rtp_out[1],\
        rtp_out[2]);
}    

template <class T1,class T2>
void xyz_to_rtp( const T1& x_in,const T1& y_in,const T1& z_in,const T2& r_out,\
    const T2& theta_out,const T2& phi_out)
{
    arctan_yx(y_in,x_in,phi_out);
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    T1 sqr_z=z_in*z_in;
    r_out=sqrt(sqr_x+sqr_y+sqr_z);
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_out);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_t
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_t
template <class T1,class T2>
void xyz_to_t_deg( const T1& x_in,const T1& y_in,const T1& z_in,T2& theta_out)
{
// Obtain theta(t) from Cartesian, or xyz, coordinates
// Theta_out is angular displacement in degrees from the positive z-axis [0,180)
    T2 theta_rads;
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_rads);
    rad_to_deg(theta_rads,theta_out);
}

// const versions
template <class T1,class T2>
void xyz_to_t_deg( const T1& x_in,const T1& y_in,const T1& z_in,\
const T2& theta_out)
{
    T2 theta_rads;
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_rads);
    rad_to_deg(theta_rads,theta_out);
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_p_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_p_deg
template <class T1,class T2>
void xyz_to_p_deg( const T1& x_in,const T1& y_in,const T1& z_in,T2& phi_out)
{
// Obtain phi(p) from Cartesian, or xyz, coordinates
// phi_out is angular displacement in degrees from the positive x-axis [0,360)
   T2 phi_rads;
   arctan_yx(y_in,x_in,phi_rads);
   rad_to_deg(phi_rads,phi_out);   
}

// const versions
template <class T1,class T2>
void xyz_to_p_deg( const T1& x_in,const T1& y_in,const T1& z_in,const T2& phi_out)
{
  T2 phi_rads;
  arctan_yx(y_in,x_in,phi_rads);
  rad_to_deg(phi_rads,phi_out); 
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_rtp_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------- xyz_to_rtp_deg
template <class T1,class T2>
void xyz_to_rtp_deg( const T1& xyz_in, T2& rtp_out)
{
// Transform Cartesian, or xyz, coordinates to spherical coordinates
// ((r),(t)heta,(p)hi)
// Theta and phi are angular displacements in degrees from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,180), phi in range [0,360)
// rtp_out is an indexable class where rtp_out[0]=(r), rtp_out[1]=(t)heta
// and rtp_out[2]=(p)hi.
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z
    xyz_to_rtp_deg(xyz_in[0],xyz_in[1],xyz_in[2],rtp_out[0],rtp_out[1],\
        rtp_out[2]);
}    

template <class T1,class T2>
void xyz_to_rtp_deg( const T1& x_in,const T1& y_in,const T1& z_in,T2& r_out,\
    T2& theta_out,T2& phi_out)
{
// Transform Cartesian, or xyz, coordinates to spherical coordinates
// ((r),(t)heta,(p)hi)
// Theta and phi are angular displacements in degrees from the positive z-axis
// and positive x-axis, respectively
// theta in range [0,180), phi in range [0,360)
    T2 phi_rads,theta_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    r_out=sqrt(sqr_x+sqr_y+z_in*z_in);
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_rads);
    rad_to_deg(theta_rads,theta_out);
}

// const versions
template <class T1,class T2>
void xyz_to_rtp_deg( const T1& xyz_in, const T2& rtp_out)
{
    xyz_to_rtp_deg(xyz_in[0],xyz_in[1],xyz_in[2],rtp_out[0],rtp_out[1],\
        rtp_out[2]);
}    

template <class T1,class T2>
void xyz_to_rtp_deg( const T1& x_in,const T1& y_in,const T1& z_in,\
    const T2& r_out, const T2& theta_out, const T2& phi_out)
{
    T2 phi_rads,theta_rads;
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    r_out=sqrt(sqr_x+sqr_y+z_in*z_in);
    arctan_yx((T1)sqrt(sqr_x+sqr_y),z_in,theta_rads);
    rad_to_deg(theta_rads,theta_out);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_to_x
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rlp_to_x
template <class T1,class T2>
void rlp_to_x( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& x_out)
{
// Obtain cartesian x-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
    T2 xy=r_in*cos(lat_in);
    x_out=xy*cos(phi_in);
}

// const version
template<class T1, class T2>
void rlp_to_x( const T1& r_in,const T1& lat_in,const T1& phi_in,const T2& x_out)
{
    T2 xy=r_in*cos(lat_in);
    x_out=xy*cos(phi_in);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_to_y
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rlp_to_y
template <class T1,class T2>
void rlp_to_y( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& y_out)
{
// Obtain cartesian y-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
    T2 xy=r_in*cos(lat_in);
    y_out=xy*sin(phi_in);
}

// const version
template <class T1,class T2>
void rlp_to_y( const T1& r_in,const T1& lat_in,const T1& phi_in,const T2& y_out)
{
    T2 xy=r_in*cos(lat_in);
    y_out=xy*sin(phi_in);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_to_z
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rlp_to_z
template <class T1,class T2>
void rlp_to_z( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& z_out)
{
// Obtain cartesian z-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
    z_out=r_in*sin(lat_in);
}

// const version
template <class T1,class T2>
void rlp_to_z( const T1& r_in,const T1& lat_in,const T1& phi_in,const T2& z_out)
{ z_out=r_in*sin(lat_in);}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- rlp_to_xyz
template <class T1,class T2>
void rlp_to_xyz( const T1& rlp_in, T2& xyz_out)
{
// Transform coordinates expressed in terms of magnitude, latitude, and
// longitude (r,(l)atitude,(p)hi) to Cartesian, or xyz, coordinates
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
// rlp_in is an indexable class where rlp_in[0]=(r), rlp_in[1]=(l)atitude,
// and rlp_in[2]=(p)hi
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    rlp_to_xyz(rlp_in[0],rlp_in[1],rlp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}

template <class T1,class T2>
void rlp_to_xyz( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform coordinates expressed in terms of magnitude, latitude, and
// longitude (r,(l)atitude,(p)hi) to Cartesian, or xyz, coordinates
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
    T2 xy=r_in*cos(lat_in);
    x_out=xy*cos(phi_in);
    y_out=xy*sin(phi_in);
    z_out=r_in*sin(lat_in);

}

// const versions
template <class T1,class T2>
void rlp_to_xyz( const T1& rlp_in, const T2& xyz_out)
{
    rlp_to_xyz(rlp_in[0],rlp_in[1],rlp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}

template <class T1,class T2>
void rlp_to_xyz( const T1& r_in,const T1& lat_in,const T1& phi_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
    T2 xy=r_in*cos(lat_in);
    x_out=xy*cos(phi_in);
    y_out=xy*sin(phi_in);
    z_out=r_in*sin(lat_in);

}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_deg_to_x
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- rlp_deg_to_x
template <class T1,class T2>
void rlp_deg_to_x( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& x_out)
{
// Obtain cartesian x-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
}

// const version
template <class T1,class T2>
void rlp_deg_to_x( const T1& r_in,const T1& lat_in,const T1& phi_in,const T2& x_out)
{
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);  
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_deg_to_y
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- rlp_deg_to_y
template <class T1,class T2>
void rlp_deg_to_y( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& y_out)
{
// Obtain cartesian y-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);
    y_out=xy*sin(phi_rads);
}

// const version
template <class T1,class T2>
void rlp_deg_to_y( const T1& r_in,const T1& lat_in,const T1& phi_in,\
const T2& y_out)
{
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);
    y_out=xy*sin(phi_rads);
 
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_to_z
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------------- rlp_to_z
template <class T1,class T2>
void rlp_deg_to_z( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& z_out)
{
// Obtain cartesian z-coordinate from coordinates expressed in terms of 
// magnitude, latitude, and longitude (r,(l)atitude,(p)hi) 
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
    T1 lat_rads;
    deg_to_rad(lat_in,lat_rads);
    z_out=r_in*sin(lat_rads);
}

// const version
template <class T1,class T2>
void rlp_deg_to_z( const T1& r_in,const T1& lat_in,const T1& phi_in,\
const T2& z_out)
{  
    T1 lat_rads;
    deg_to_rad(lat_in,lat_rads);
    z_out=r_in*sin(lat_rads);
}

// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rlp_deg_to_xyz
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- rlp_deg_to_xyz
template <class T1,class T2>
void rlp_deg_to_xyz( const T1& rlp_in, T2& xyz_out)
{
// Transform coordinates expressed in terms of magnitude, latitude, and
// longitude (r,(l)atitude,(p)hi) to Cartesian, or xyz, coordinates
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
// rlp_in is an indexable class where rlp_in[0]=(r), rlp_in[1]=(l)atitude,
// and rlp_in[2]=(p)hi
// Similarly, xyz_out[0]=x, xyz_out[1]=y, xyz_out[2]=z
    rlp_deg_to_xyz(rlp_in[0],rlp_in[1],rlp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}

template <class T1,class T2>
void rlp_deg_to_xyz( const T1& r_in,const T1& lat_in,const T1& phi_in,T2& x_out,\
    T2& y_out,T2& z_out)
{
// Transform coordinates expressed in terms of magnitude, latitude, and
// longitude (r,(l)atitude,(p)hi) to Cartesian, or xyz, coordinates
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
    y_out=xy*sin(phi_rads);
    z_out=r_in*sin(lat_rads);
}

// const versions
template <class T1,class T2>
void rlp_deg_to_xyz( const T1& rlp_in, const T2& xyz_out)
{
    rlp_deg_to_xyz(rlp_in[0],rlp_in[1],rlp_in[2],xyz_out[0],xyz_out[1],\
        xyz_out[2]);
}

template <class T1,class T2>
void rlp_deg_to_xyz( const T1& r_in,const T1& lat_in,const T1& phi_in,\
    const T2& x_out, const T2& y_out, const T2& z_out)
{
    T1 lat_rads,phi_rads;
    deg_to_rad(lat_in,lat_rads);
    T2 xy=r_in*cos(lat_rads);
    deg_to_rad(phi_in,phi_rads);
    x_out=xy*cos(phi_rads);
    y_out=xy*sin(phi_rads);
    z_out=r_in*sin(lat_rads);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_l
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_l
template <class T1,class T2>
void xyz_to_l( const T1& x_in,const T1& y_in,const T1& z_in,T2& lat_out)
{
// Obtain latitude(l) { [pi/2,-pi/2) }, from Cartesian, or xyz, coordinates
// latitude = pi/2 - theta, where theta is angular displacement in radians from 
// the positive z-axis
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    lat_out=atan2(z_in,sqrt(sqr_x+sqr_y)); 
}

// const versions
template <class T1,class T2>
void xyz_to_l( const T1& x_in,const T1& y_in,const T1& z_in,const T2& lat_out)
{
   T1 sqr_x=x_in*x_in;
   T1 sqr_y=y_in*y_in;
   lat_out=atan2(z_in,sqrt(sqr_x+sqr_y)); 
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_rlp
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_rlp
template <class T1,class T2>
void xyz_to_rlp( const T1& xyz_in, T2& rlp_out)
{
// Transform Cartesian, or xyz, coordinates to coordinates expressed in terms
// of magnitude, latitude, and longitude (r,(l)atitude,(p)hi)
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
// rlp_out is an indexable class where rlp_out[0]=(r), rlp_out[1]=(l)atitude,
// and rlp_out[2]=(p)hi
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z

    xyz_to_rlp(xyz_in[0],xyz_in[1],xyz_in[2],rlp_out[0],rlp_out[1],\
        rlp_out[2]);
}

template <class T1,class T2>
void xyz_to_rlp(const T1& x_in, const T1& y_in, const T1& z_in, T2& r_out,\
    T2& lat_out, T2& phi_out)
{
// Transform Cartesian, or xyz, coordinates to coordinates expressed in terms
// of magnitude, latitude, and longitude (r,(l)atitude,(p)hi)
// latitude = pi/2 - theta, where theta is the angular distance in radians 
// from the positive z axis, phi is the angular distance in radians from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [pi/2,-pi/2), phi in range [0,2pi)
    T1 x_sqr=x_in*x_in;
    T1 y_sqr=y_in*y_in;
    r_out=sqrt(x_sqr+y_sqr+z_in*z_in);
    lat_out=atan2(z_in,sqrt(x_sqr+y_sqr));
    arctan_yx(y_in,x_in,phi_out);
}

// const versions
template <class T1,class T2>
void xyz_to_rlp( const T1& xyz_in, const T2& rlp_out)
{
    xyz_to_rlp(xyz_in[0],xyz_in[1],xyz_in[2],rlp_out[0],rlp_out[1],\
        rlp_out[2]);
}

template <class T1,class T2>
void xyz_to_rlp(const T1& x_in, const T1& y_in, const T1& z_in,\
    const T2& r_out, const T2& lat_out, const T2& phi_out)
{
    T1 x_sqr=x_in*x_in;
    T1 y_sqr=y_in*y_in;
    r_out=sqrt(x_sqr+y_sqr+z_in*z_in);
    lat_out=atan2(z_in,sqrt(x_sqr+y_sqr));
    arctan_yx(y_in,x_in,phi_out);
}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_l_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// --------------------------------------------------------------- xyz_to_l_deg
template <class T1,class T2>
void xyz_to_l_deg( const T1& x_in,const T1& y_in,const T1& z_in,T2& lat_out)
{
// Obtain latitude(l) { [90,-90) }, from Cartesian, or xyz, coordinates
// latitude = pi/2 - theta, where theta is angular displacement in degrees from 
// the positive z-axis
    T2 lat_rads;
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    lat_rads=atan2(z_in,sqrt(sqr_x+sqr_y));
    rad_to_deg(lat_rads,lat_out); 
}

// const versions
template <class T1,class T2>
void xyz_to_l_deg( const T1& x_in,const T1& y_in,const T1& z_in,const T2& lat_out)
{
    T2 lat_rads;
    T1 sqr_x=x_in*x_in;
    T1 sqr_y=y_in*y_in;
    lat_rads=atan2(z_in,sqrt(sqr_x+sqr_y));
    rad_to_deg(lat_rads,lat_out); 

}
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xyz_to_rlp_deg
// ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
// ----------------------------------------------------------- xyz_to_rlp_deg
template <class T1,class T2>
void xyz_to_rlp_deg( const T1& xyz_in, T2& rlp_out)
{
// Transform Cartesian, or xyz, coordinates to coordinates expressed in terms
// of magnitude, latitude, and longitude (r,(l)atitude,(p)hi)
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
// rlp_out is an indexable class where rlp_out[0]=(r), rlp_out[1]=(l)atitude,
// and rlp_out[2]=(p)hi
// Similarly, xyz_in[0]=x, xyz_in[1]=y, xyz_in[2]=z

    xyz_to_rlp_deg(xyz_in[0],xyz_in[1],xyz_in[2],rlp_out[0],rlp_out[1],\
        rlp_out[2]);
}

template <class T1,class T2>
void xyz_to_rlp_deg(const T1& x_in, const T1& y_in, const T1& z_in, T2& r_out,\
    T2& lat_out, T2& phi_out)
{
// Transform Cartesian, or xyz, coordinates to coordinates expressed in terms
// of magnitude, latitude, and longitude (r,(l)atitude,(p)hi)
// latitude = 90 - theta, where theta is the angular distance in degrees 
// from the positive z axis, phi is the angular distance in degrees from the
// positive x axis and is equal to the longitude, and r is the magnitude of 
// the radial distance.
// latitude in range [90,-90), phi in range [0,360)
    T1 x_sqr=x_in*x_in;
    T1 y_sqr=y_in*y_in;
    r_out=sqrt(x_sqr+y_sqr+z_in*z_in);
    T2 lat_rads,phi_rads;
    lat_rads=atan2(z_in,sqrt(x_sqr+y_sqr));
    rad_to_deg(lat_rads,lat_out);
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
}

// const versions
template <class T1,class T2>
void xyz_to_rlp_deg( const T1& xyz_in,  const T2& rlp_out)
{
    xyz_to_rlp_deg(xyz_in[0],xyz_in[1],xyz_in[2],rlp_out[0],rlp_out[1],\
        rlp_out[2]);
}

template <class T1,class T2>
void xyz_to_rlp_deg(const T1& x_in, const T1& y_in, const T1& z_in, \
    const T2& r_out, const T2& lat_out, const T2& phi_out)
{
    T1 x_sqr=x_in*x_in;
    T1 y_sqr=y_in*y_in;
    r_out=sqrt(x_sqr+y_sqr+z_in*z_in);
    T2 lat_rads,phi_rads;
    lat_rads=atan2(z_in,sqrt(x_sqr+y_sqr));
    rad_to_deg(lat_rads,lat_out);
    arctan_yx(y_in,x_in,phi_rads);
    rad_to_deg(phi_rads,phi_out);
}

}  // end namespace QDOS
#endif
