/*
   SPHERE.I
   Routines for manipulating spherical coordinates

    Copyright (C) 2008  Thibaut Paumard <paumard@users.sourceforge.net>

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

   History:

   v. 0.2, Tue, 05 Feb 2008 14:13:50 +0100 Thibaut Paumard
      added functions: spldefault, sphere_area, sphere_mesh, splmk,
        sphere_test
      added keywords: pole, rotate (thus todo item is done)
      todo: rotate curves nicely
   
   v. 0.1, Mon, 28 Jan 2008 20:33:50 +0100 Thibaut Paumard
      initial release. Todo: support rotations in sphere_proj().
    
*/

extern sphere;
/* DOCUMENT sphere.i -- spherical geometry

    sphere2cart:  convert sperical coordinates to cartesian;
    cart2sphere:  convert cartesian coordinates to spherical;
    dotproduct:   scalar products of 3D vectors;
    crossproduct: cross product of 3D vectors;
    norm:         Euclidian norm;
    sphere_mesh:  return a spherical grid;
    sphere_area:  compute area of each cell in a spherical grid;
    sphere_proj:  compute projection from a sphere to a plane;
    splg:         display a projected spherical graph;
    splm:         display a projected spherical mesh on paper;
    splmk:        display projected spherical points on paper;
    spli:         display a projected spherical image on paper;
    spldefault:   set defaults for spl* functions;
    shisto:       compute spherical histograms;
    sphere_test:  test suite for this package.

   Spherical coordinates here are (r, theta, phi) where r is the
   distance from point M to center O, theta is the angle between the
   plane containing the z and x axes and the plane containing the z
   axis and OM,(in [-pi, pi]) and phi is the angle between Oz and OM
   (in [0, pi]). All angles are expressed in radians.

   Most routines are meant to deal with a sphere. r=1 is then implied
   and not included in the calling sequence. Coordinates (either
   spherical or cartesian) output in a single array by the routines,
   the last dimension being of length 3 (or 2 for 2D spherical
   coordinates). For instance, sphere2cart(r, theta, phi)(..,2)==y.

   When coordinates are expected as inputs, they may either be given
   in 2 or 3 individual array, or as a single array: for instance
   sphere2cart(r, theta, phi)==sphere2cart([r, theta, phi]).
    
   An important concept here is that of coordinate grid, represented
   by two arrays theta and phi which form a mesh. theta and phi have
   identical dimensionality. theta(i,) does not depend on i and
   phi(,j) does not depend on j. theta(,j) is monotonically increasing
   and phi(i,) is monotonic. The grid needs not be regular and may not
   cover the entire sphere.

   SEE ALSO:sphere2cart, cart2sphere, dotproduct, crossproduct, norm,
            sphere_proj, splg, splm, spli, shisto
 */

func sphere2cart(r, theta, phi) {
/* DOCUMENT x_y_z=sphere2cart(r, theta, phi)
         or x_y_z=sphere2cart(r_theta_phi)

    convert spherical coordinates to cartesian.
         
   SEE ALSO: sphere, cart2sphere
 */
  if (is_void(theta)) {
    theta=r(..,-1);
    phi=r(..,0);
    if (dimsof(r)(dimsof(r)(1))==3) r=r(..,1);
    else r=1;
  } else if (is_void(phi)) {
    phi=theta;
    theta=r;
    r=1;
  }
  x=r*cos(theta)*sin(phi);
  y=r*sin(theta)*sin(phi);
  z=r*cos(phi);
  return [x,y,z];
}

func cart2sphere(x, y, z) {
/* DOCUMENT r_theta_phi=cart2sphere(x, y, z)
         or r_theta_phi=cart2sphere(x_y_z)

    convert cartesian coordinates to spherical.
         
   SEE ALSO: sphere, sphere2cart
 */
  if (is_void(y)) {
    z=x(..,3);
    y=x(..,2);
    x=x(..,1);
  }
  r=sqrt(x^2+y^2+z^2);
  theta=atan(y,x);
  phi=atan(sqrt(x^2+y^2),z);
  return [r,theta,phi];
}

func dotproduct(v1,v2) {
/* DOCUMENT value=dotproduct(v1,v2)

    compute the dot or scalar product of pairs of vectors. V1 and V2
    must be conformable arrays with 3 as their last dimension. For
    instance, an NxPx3 array represents NxP vectors. VALUE will then
    be an NxP array containing the products.
     
   SEE ALSO: crossproduct, norm
 */
  return v1(..,1)*v2(..,1)+
         v1(..,2)*v2(..,2)+
         v1(..,3)*v2(..,3);
}

func crossproduct(v1,v2) {
/* DOCUMENT vector=crossproduct(v1,v2)

    compute the cross product of pairs of vectors. V1 and V2 must be
    conformable arrays with 3 as their last dimension. For instance,
    an NxPx3 array represents NxP vectors. VECTOR will then be also an
    NxPx3 array containing the products.
     
   SEE ALSO: dotproduct, norm
 */
  return [v1(..,2)*v2(..,3)-v1(..,3)*v2(..,2),
          v1(..,3)*v2(..,1)-v1(..,1)*v2(..,3),
          v1(..,1)*v2(..,2)-v1(..,2)*v2(..,1)];
}

func norm(v) {
/* DOCUMENT value=norm(v)

    compute the Euclidian norm of one of more 3D vectors. V must be an
    array with a last dimension of 3. For instance, an NxPx3 array
    represents NxP vectors. VALUE will then be an NxP array containing
    the norms.
     
   SEE ALSO: dotproduct, crossproduct
 */
  return sqrt(dotproduct(v,v));
}

func sphere_proj(theta,phi,&z,proj=,rotate=, pole=) {
/* DOCUMENT x_y=sphere_proj(theta, phi)
         or x_y=sphere_proj([theta, phi])
   
    project a 2D object from a sphere onto a plane. The inputs are the
    spherical coordinates of the points to project (without R; see
    sphere). The output in a N1xN2x...xNnx2 array giving the x and y
    2D cartesian coordinates of each projected point.

   KEYWORDS:
   
    proj= a string selecting the of projection. One of "azimuthal
          equidistant" (or "aed"; default), "azimuthal equal area" (or
          "aea"), "aitoff", or "hammer".
    
   SEE ALSO: splg, splm, spli
 */
  extern _sphere_default_rot, _sphere_default_proj, _sphere_default_pole;

  // keywords init
  if (is_void(proj)) proj=_sphere_default_proj;
  if (is_void(rotate)) rotate=_sphere_default_rot;
  if (is_void(pole)) pole=_sphere_default_pole;

  // get possibly get THETA and PHI from single THETA_PHI
  _sphere_get_theta_phi,theta,phi;

  // pole
  if (pole=="south") {
    phi=pi-phi;
    if (numberof(theta)>=2) theta=-theta(0:1:-1,..);
    else theta=-theta;
    if (numberof(phi)>=2) phi=phi(0:1:-1,..);
    if (numberof(z)>=2) z=z(0:1:-1,..);
  }
  
  // rotation (ROTATE keyword)
  theta=theta-rotate;
  if (numberof(theta)==1) {
    theta=(theta+pi)%(2*pi)-pi;
  } else {
    mnt0=theta(min);
    theta+=(mnt0+pi)%(2*pi)-pi-mnt0;
    if (mnt0<-pi) theta+=2*pi;
    
    if (allof(theta>2*pi)) theta-=2*pi;
    if (dimsof(theta)(2) >= 4 && anyof(theta(2:-2,..)>pi)) {
      br=min(where2(theta>(pi))(1));
      theta=transpose(_(transpose(theta(br-1:-1,..))-2*pi,
                        transpose(theta(1:br-1,..))));
      phi=transpose(_(transpose(phi(br-1:-1,..)),
                      transpose(phi(1:br-1,..))));
      if (!is_void(z)) z=transpose(_(transpose(z(br-1:,..)),
                                     transpose(z(1:br-2,..))));
    }
  }
  //
  if (anyof(proj==["aitoff","hammer"])) {
    theta=theta/2.;
    around=[0,pi/2];
    cart=sphere2cart(theta,phi);
    e1=sphere2cart(around);
    e1=e1/norm(e1);
    e2=[0.,0.,1.]-e1*dotproduct([0.,0.,1.],e1);
    e2=e2/norm(e2);
    e3=crossproduct(e1,e2);
    e3=e3/norm(e3);
    phi=acos(dotproduct(cart,e1));
    x=dotproduct(cart,e2);
    y=dotproduct(cart,e3);
    theta=atan(y,x);
  }
  if (anyof(proj==["aed","azimuthal equidistant","aitoff"])) {
    x=(phi)*cos(theta+pi/2);
    y=(phi)*sin(theta+pi/2);
  } else if (anyof(proj==["aea","azimuthal equal area","hammer"])) {
    rho=sin(phi/2);
    x=rho*cos(theta+pi/2);
    y=rho*sin(theta+pi/2);
  } else error,"unkown projection: "+proj;
  if (anyof(proj==["aitoff","hammer"])) x=2*x;
  return [x,y];
}

extern splg, splm, splmk;
/* DOCUMENT splg, theta, phi
         or splm, theta, phi
         or splmk, theta, phi

    equivalent to plg, plm and plmk for data mapped on a sphere. The
    sphere is projected on a plane using sphere_proj() and the result
    in plotted using pl* routines.

   KEYWORDS:

    Each routine accepts keywords to pass through to both
    sphere_proj() and to the correspong planar plotting routine.
     
   SEE ALSO: plg, plm, plmk, splf, sphere_proj, sphere
 */

func splg(theta, phi, proj=, pole=, rotate=,
          legend=, hide=,
          type=, width=, color=, closed=, smooth=,
          marks=, marker=, mspace=, mphase=,
          rays=, arrowl=, arroww=, rspace=, rphase=
          ) {
  xy=sphere_proj(theta,phi,proj=proj,pole=pole,rotate=rotate);
  plg,xy(..,2),xy(..,1),legend=legend, hide=hide,
    type=type, width=width, color=color, closed=closed, smooth=smooth,
    marks=marks, marker=marker, mspace=mspace, mphase=mphase,
    rays=rays, arrowl=arrowl, arroww=arroww, rspace=rspace, rphase=rphase;
}

func splm(theta, phi, proj=, pole=, rotate=,
          legend=, hide=, type=, width=, color=, region=
          ) {
  xy=sphere_proj(theta,phi,proj=proj,pole=pole,rotate=rotate);
  plm,xy(..,2),xy(..,1),legend=legend, hide=hide,
    type=type, width=width, color=color,
    region=region;
}

func splmk(theta, phi, proj=, pole=, rotate=,
           marker=, msize=, width=, color=
          ) {
  xy=sphere_proj(theta,phi,proj=proj,pole=pole,rotate=rotate);
  plmk,xy(..,2),xy(..,1), marker=marker, msize=msize, width=width, color=color;
}

func spli(z, theta, phi, proj=, rotate=,
          legend=, hide=, top=, cmin=, cmax=,
          edges=, ecolor=, ewidth=) {
/* DOCUMENT spli, im
         or splf, im, theta, phi

    splf and spli are actually the same function. It has two names
    because it is meant to be the spherical equivalent to both pli and
    plf. As usual (see sphere), THETA and PHI can be specified in a
    single parameter as [theta, phi]. THETA and PHI default to
    regularly-spaced arrays which will map the image IM to cover the
    entire sphere. The sphere is projected using sphere_proj() and the
    projected image is drawn using plf.

   KEYWORDS:

    all the keywords accepted by sphere_proj and plf (and thus by
    pli).

   SEE ALSO: sphere, sphere_proj, pli, plf
 */
  if (is_void(theta)) {
    grid=dimsof(z)(-1:);
    theta=array(span(-pi,pi,grid(1)+1),grid(2)+1);
    phi=transpose(array(span(0,pi,grid(2)+1),grid(1)+1));
  }
  xy=sphere_proj(theta,phi,z,proj=proj,rotate=rotate);
  plf,z,xy(..,2),xy(..,1),
    legend=legend, hide=hide, top=top, cmin=cmin, cmax=cmax,
    edges=edges, ecolor=ecolor, ewidth=ewidth;
}
splf=spli;

func shisto(theta,phi,grid_theta,grid_phi) {
/* DOCUMENT hist=shisto(theta, phi, grid_theta, grid_phi)
         or hist=shisto(theta_phi, , grid_theta_phi)
         ...

    compute the histogram HIST of the positions (THETA, PHI) falling
    into each cell of the mesh defined in spherical coordinates by
    GRID_THETA and GRID_PHI. As usual (see sphere), theta and phi may
    be specified in a single array THETA_PHI=[THETA, PHI] but beware
    of the two commas before GRID_THETA in this case.

    GRID_THETA and GRID_PHI (which may be specified ans
    GRID_THETA_PHI=[GRID_THETA, GRID_PHI]) must be either monotonical
    vectors (1-dimensional arrays) or scalars. Setting GRID_THETA
    (resp. GRID_PHI) to the scalar N is equivalent to setting it to
    span(-pi, pi, N+1) (resp. span(0, pi, N+1). If only GRID_THETA is
    specified and is a scalar, then GRID_PHI is set to
    (floor(GRID_THETA/2)).

    Finally, GRID_THETA and GRID_PHI may be 2-dimensional and define a
    mesh, but this mesh must be a grid: it is necessary that for any
    index N, GRID_THETA(,N)==GRID_THETA(,1) and conversely
    GRID_PHI(N,)==GRID_PHI(1,).

   SEE ALSO: histo, plh, pli
 */
  if (is_void(grid_theta)) grid_theta=_sphere_default_mesh;
  if (is_void(grid_phi)) {
    if (numberof(grid_theta)==1) grid_phi=floor(grid_theta/2);
    else {
      grid_phi=grid_theta(..,2);
      grid_theta=grid_theta(..,1);
    }
  }
  if (numberof(grid_theta)==1) grid_theta=span(-pi,pi,grid_theta+1);
  else if (dimsof(grid_theta)(1)==2) grid_theta=grid_theta(,1);
  if (numberof(grid_phi)==1) grid_phi=span(0,pi,grid_phi+1);
  else if (dimsof(grid_phi)(1)==2) grid_phi=grid_phi(1,);
  _sphere_get_theta_phi,theta,phi;
  nx=dimsof(grid_theta)(2)-1;
  ny=dimsof(grid_phi)(2)-1;
  out=array(long,nx,ny);
  for (x=1;x<=nx;x++) {
    xmin=grid_theta(min:x:x+1);
    xmax=grid_theta(max:x:x+1);
    ind=where(theta>=xmin & theta<xmax);
    if (numberof(ind)) {
      for (y=1;y<=ny;y++) {
        ymin=grid_phi(min:y:y+1);
        ymax=grid_phi(max:y:y+1);
        out(x,y)=sum(phi(ind)>=ymin & phi(ind)<ymax);
      }
    }
  }
  return out;
}

func sphere_mesh(theta,phi) {
/* DOCUMENT mesh_theta_phi=sphere_mesh(theta_spec, phi_spec)
         or mesh_theta_phi=sphere_mesh([theta_spec, phi_spec])

    construct a regular spherical mesh. Returns a NxMx2 suitable for
    input to splf. THETA_SPEC and PHI_SPEC are either scalars or
    vectors. If they are both scalars, then N=THETA_SPEC, M=PHI_SPEC
    and the mesh is regular and covers the entire sphere. If they are
    both vectors, N=numberof(THETA_SPEC), M=numberof(PHI_SPEC) and the
    mesh is built by replicating M times THETA_SPEC and N times
    PHI_SPEC. So the output mesh may be only semi-regular and may not
    cover the entire sphere.

    If THETA_SPEC is not specified, it defaults to 36. If PHI_SPEC is
    not specified and THETA is scalar, PHI_SPEC defaults to
    THETA_SPEC/2. Therefore, the default mesh covers the entire sphere
    with parallels and meridians, and there is one prrallel/meridian
    every 10°.
         
   SEE ALSO: sphere, splm
 */
  if (is_void(theta)) theta=_sphere_default_mesh;
  if (is_void(phi)) {
    if (numberof(theta)==1) phi=long(theta)/2;
    else { phi=theta(..,2); theta=theta(..,1); }
  }
  if (numberof(theta)==1) theta=span(-pi,pi,theta+1);
  if (numberof(phi)==1) phi=span(0,pi,phi+1);
  nt=numberof(theta);
  theta=array(theta,numberof(phi));
  phi=transpose(array(phi,nt));
  return [theta, phi];
}

func sphere_area(mth,mph) {
/* DOCUMENT area=spere_area(grid_theta, grid_phi)
         or area=spere_area(grid_theta_phi)

     compute exact area of each cell in the mesh defined by grid_theta
     and grid_phi

   SEE ALSO: sphere, sphere_mesh
 */
  _sphere_get_theta_phi,mth,mph;
  return (mth(2:,:-1)-mth(:-1,:-1))*(cos(mph(:-1,:-1))-cos(mph(:-1,2:)));
}

func _sphere_get_theta_phi(&theta,&phi) {
  extern _sphere_default_mesh;
  if (is_void(theta)) theta=_sphere_default_mesh;
  if (is_void(phi)) {
    phi=theta(..,2);
    theta=theta(..,1);
  }
}

extern _sphere_default_mesh, _sphere_default_rot,
  _sphere_default_proj, _sphere_default_pole;

func spldefault(void,theta_phi=,proj=,rotate=,pole=) {
/* DOCUMENT spldefault, keyword=value

    set defaults for sherical plots

   KEYWORDS:
    theta_phi= [theta, phi] default coordinate grid
    proj, rotate, pole: see sphere_proj

   SEE ALSO: pldefault, sphere, sphere_proj
 */
  extern _sphere_default_mesh, _sphere_default_rot,
    _sphere_default_pole, _sphere_default_proj;
  if (!is_void(pole)) _sphere_default_pole=pole;
  if (!is_void(theta_phi)) _sphere_default_mesh=theta_phi;
  if (!is_void(rotate)) _sphere_default_rot=rotate;
  if (!is_void(proj)) _sphere_default_proj=proj;
}

spldefault,theta_phi=sphere_mesh(36),rotate=0.,proj="aed",pole="north";

func sphere_test {
/* DOCUMENT sphere_test

    Test suite for sphere.i

   SEE ALSO: sphere
 */
  write,format="Old default projection:              %s\n",
    _sphere_default_proj;
  write,format="Old default pole:                    %s\n",
    _sphere_default_pole;
  write,format="Old default rotation angle:          %g x pi\n",
    _sphere_default_rot/pi;
  write,format="Shape of old default grid: %i x %i\n",
    dimsof(_sphere_default_mesh)(2)-1, dimsof(_sphere_default_mesh)(3)-1;

  default_proj=_sphere_default_proj;
  default_pole=_sphere_default_pole;
  default_rot =_sphere_default_rot ;
  default_mesh=_sphere_default_mesh;

  write,format="%s\n","Running spldefault...";
  spldefault,theta_phi=sphere_mesh(36),rotate=-pi/2,proj="hammer",pole="north";

  write,format="New default projection:              %s\n",
    _sphere_default_proj;
  write,format="New default pole:                    %s\n",
    _sphere_default_pole;
  write,format="New default rotation angle:          %g x pi\n",
    _sphere_default_rot/pi;
  write,format="Shape of new default grid: %i x %i\n",
    dimsof(_sphere_default_mesh)(2)-1, dimsof(_sphere_default_mesh)(3)-1;
  
  nmc=long(1e5);
  write,format="%s %i %s\n","Generating",nmc,"uniformly distributed point on the sphere";
  thetas=random(nmc)*2*pi-pi;
  phis=asin(random(nmc)*2-1)+pi/2;
  write,format="%s\n","Running shisto...";
  hist=shisto(thetas,phis);
  write,format="Number of random points per cell has min %i and max %i\n",min(hist),max(hist);
  write,format="%s\n","Running sphere_area...";
  area=sphere_area();
  write,format="Cell area has min %g and max %g\n",min(area),max(area);
  write,format="Density of random points has mean %g and RMS %g\n",avg(hist/area),reform(hist/area,numberof(area))(rms);


  write, format="%s\n","";
  write, format="%s\n",
    "Will now plot 6 times (once per projection and pole) the following:";
  write, format="%s\n",
    "  * density map of random points;";
  write, format="%s\n",
    "  * spherical coordinates grid mesh;";
  write, format="%s\n",
    "  * three point with (theta, phi)=(0,0), (0,pi/2), (pi/2, pi/2);";
  write, format="%s\n",
    "  * curve defined as theta=phi.";
  write, format="%s\n",
    "The curve appears to connect the North pole to the South pole in a straight";
  write, format="%s\n",
    "line, this is a limitation of the sphere_proj(rotate=...) function.";
  proj=["aed","aed","aitoff","aea","aea","hammer"];
  pole=["north","south","north","north","south","north"];
  thetas=phis=span(0,pi,181);
  for (win=1;win<=6;win++) {
    window, win, style="nobox.gs";
    limits,square=1;
    limits;
    write, format="%s\n","";
    write,format="Opening window %i\n",win;
    write,format="Setting projection type to %s and pole to %s\n",proj(win),pole(win);
    spldefault,proj=proj(win),pole=pole(win);
    write,format="New default projection:              %s\n",
      _sphere_default_proj;
    write,format="New default pole:                    %s\n",
      _sphere_default_pole;
    write,format="%s\n","Running spli...";
    spli,hist/area;
    write,format="%s\n","Running splm...";
    splm;
    write,format="%s\n","Running splmk...";
    splmk,0,0,marker=0;
    splmk,0,pi/2,marker=1;
    splmk,pi/2,pi/2,marker=2;
    splg,thetas,phis,color="red";
  }

  _sphere_default_proj=_default_proj;
  _sphere_default_pole=_default_pole;
  _sphere_default_rot =_default_rot ;
  _sphere_default_mesh=_default_mesh;
}