sphcoords.h

Go to the documentation of this file.

/* This file is part of the Tomographer project, which is distributed under the
 * terms of the MIT license.
 *
 * The MIT License (MIT)
 *
 * Copyright (c) 2016 ETH Zurich, Institute for Theoretical Physics, Philippe Faist
 * Copyright (c) 2017 Caltech, Institute for Quantum Information and Matter, Philippe Faist
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef SPHCOORDS_H
#define SPHCOORDS_H

#include <cmath>

#include <Eigen/Eigen>

#include <tomographer/tools/cxxutil.h> // tomographer_assert()


namespace Tomographer {
namespace MathTools {
namespace SphCoords {


template<typename Der1, typename Der2>
inline void cart_to_sph(Eigen::MatrixBase<Der2>& rtheta, const Eigen::MatrixBase<Der1>& cart)
{
  { using namespace Eigen; EIGEN_STATIC_ASSERT_LVALUE(Der2); }
  typedef typename Eigen::MatrixBase<Der1>::Scalar Scalar;

  tomographer_assert(cart.cols() == 1 && rtheta.cols() == 1);
  tomographer_assert(cart.rows() == rtheta.rows());

  const Eigen::Index ds = cart.rows()-1; // dimension of the sphere

  // see http://people.sc.fsu.edu/~jburkardt/cpp_src/hypersphere_properties/hypersphere_properties.cpp
  // and http://en.wikipedia.org/wiki/N-sphere#Spherical_coordinates
  //
  // Remember that artcot(x/y) for us can be replaced by atan2(y, x)

  //
  // R coordinate -- rtheta(0)
  //
  rtheta(0) = sqrt(cart.transpose()*cart);

  //
  // theta coordinates -- rtheta(1..ds)
  //

  // REMEMBER: .block() arguments are (OFFSET, SIZE) and not (STARTOFFSET, ENDOFFSET)
  //           .segment() arguments are (STARTPOS, LEN) and not (STARTPOS, ENDPOS)

  // all except last
  rtheta.segment(1,ds-1).setZero();

  Eigen::Index i;

  Scalar val = pow(cart(ds),2);
  for (i = ds-1; i >= 1; --i) {
    val = val + pow(cart(i),2);
    rtheta(i) = val;
  }

  for (i = 1; i < ds; i++) {
    rtheta(i) = atan2(sqrt(rtheta(i)), cart(i-1));
  }

  // last angle, theta(ds-1) == rtheta(ds) :
  val = sqrt(pow(cart(ds), 2) + pow(cart(ds-1), 2)) + cart(ds-1);

  rtheta(ds) = 2.0 * atan2(cart(ds), val);
}


template<typename Der1, typename Der2>
inline void sphsurf_to_cart(Eigen::MatrixBase<Der2>& cart, const Eigen::MatrixBase<Der1>& theta,
                            const typename Eigen::MatrixBase<Der1>::Scalar R = 1.0)
{
  { using namespace Eigen; EIGEN_STATIC_ASSERT_LVALUE(Der2); }
  // same as sph_to_cart, except that only the angles are given, and R is given
  // separately, defaulting to 1.

  //  typedef typename Eigen::MatrixBase<Der1>::Scalar Scalar;

  tomographer_assert(cart.cols() == 1 && theta.cols() == 1);
  tomographer_assert(cart.rows() == theta.rows() + 1);

  const Eigen::Index ds = theta.rows(); // dimension of the sphere

  // see http://people.sc.fsu.edu/~jburkardt/cpp_src/hypersphere_properties/hypersphere_properties.cpp

  cart.setConstant(R); // R coordinate
  
  Eigen::Index i;
  for (i = 0; i < ds; ++i) {
    cart(i) *= cos(theta(i));
    cart.segment(i+1, ds+1 - (i+1)) *= sin(theta(i));
  }
}

template<typename Der1, typename Der2>
inline void sph_to_cart(Eigen::MatrixBase<Der2>& cart, const Eigen::MatrixBase<Der1>& rtheta)
{
  { using namespace Eigen; EIGEN_STATIC_ASSERT_LVALUE(Der2); }
  //  typedef typename Eigen::MatrixBase<Der1>::Scalar Scalar;

  tomographer_assert(cart.cols() == 1 && rtheta.cols() == 1);
  tomographer_assert(cart.rows() == rtheta.rows());

  const Eigen::Index ds = rtheta.rows()-1; // dimension of the sphere

  // see http://people.sc.fsu.edu/~jburkardt/cpp_src/hypersphere_properties/hypersphere_properties.cpp

  sphsurf_to_cart(cart, rtheta.segment(1,ds), rtheta(0));
}



template<typename Der1>
inline typename Eigen::MatrixBase<Der1>::Scalar cart_to_sph_jacobian(const Eigen::MatrixBase<Der1>& rtheta)
{
  const Eigen::Index ds = rtheta.rows()-1; // dimension of the sphere
  typename Eigen::MatrixBase<Der1>::Scalar jac = pow(rtheta(0), (int)ds); // r^{n-1}

  Eigen::Index i;
  for (i = 0; i < ds-1; ++i) {
    jac *= pow(sin(rtheta(1+i)), (int)(ds-1-i));
  }

  return jac;
}

template<typename Der1>
inline typename Eigen::MatrixBase<Der1>::Scalar surf_sph_jacobian(const Eigen::MatrixBase<Der1>& theta)
{
  const Eigen::Index ds = theta.rows();

  typename Eigen::MatrixBase<Der1>::Scalar jac = 1;

  Eigen::Index i;
  for (i = 0; i < ds-1; ++i) {
    jac *= pow(sin(theta(i)), (int)(ds-1-i));
  }

  return jac;
}


template<typename Der1, typename Der2>
inline void sphsurf_diffjac(Eigen::ArrayBase<Der1> & dxdtheta, const Eigen::MatrixBase<Der2>& theta)
{
  //
  // For this calculation & the 2nd derivatives, see [Drafts&Calculations, Vol. V,
  // 13.02.2015 (~60%)]
  //

  { using namespace Eigen; EIGEN_STATIC_ASSERT_LVALUE(Der1); }
  typedef typename Eigen::ArrayBase<Der1>::Scalar Scalar;

  const Eigen::Index ds = theta.rows();
  const Eigen::Index n = ds + 1;

  tomographer_assert(theta.cols() == 1);
  tomographer_assert(dxdtheta.rows() == (int)n);
  tomographer_assert(dxdtheta.cols() == (int)ds);

  Eigen::Array<Scalar, Der2::RowsAtCompileTime, 1> sintheta(theta.rows());
  sintheta = theta.array().sin();
  Eigen::Array<Scalar, Der2::RowsAtCompileTime, 1> costheta(theta.rows());
  costheta = theta.array().cos();

  Eigen::Index i, k, mm;
  for (i = 0; i < ds; ++i) {
    for (k = 0; k < n; ++k) {
      //std::cout << "k,i = "<< k<< ", "<< i << "\n";
      Scalar val = 0.0;
      // pick the right case
      if (i > k) {
        //std::cout << "[situation i>k]\n";
        val = 0.0;
      } else if (i == k && k < n - 1) {
        //std::cout << "[situation i == k && k < n-1]\n";
        val = -1;
        for (mm = 0; mm <= i; ++mm) {
          val *= sintheta(mm);
        }
      } else if (i == k && k == n-1) {
        //std::cout << "[situation i == k && k == n-1]\n";
        val = costheta(i);
        for (mm = 0; mm < i; ++mm) {
          val *= sintheta(mm);
        }
      } else if (i < k && k < n-1) {
        //std::cout << "[situation i < k && k < n-1]\n";
        val = costheta(i) * costheta(k);
        for (mm = 0; mm < k; ++mm) {
          if (mm == i) {
            continue;
          }
          val *= sintheta(mm);
        }
      } else { // i < k && k == n-1
        //std::cout << "[else: {situation i < k && k == n-1}]\n";
        val = costheta(i);
        for (mm = 0; mm < n-1; ++mm) {
          if (mm == i) {
            continue;
          }
          val *= sintheta(mm);
        }
      }
      dxdtheta(k,i) = val;
      //std::cout << "dxdtheta(k="<<k<<",i="<<i<<") = "<< val << "\n";
    }
  }
}

template<typename Der1, typename Der2>
inline void sphsurf_diffjac2(Eigen::ArrayBase<Der1> & ddxddtheta, const Eigen::MatrixBase<Der2>& theta)
{
  //
  // For this calculation & the 2nd derivatives, see [Drafts&Calculations, Vol. V,
  // 13.02.2015 (~60%)]
  //

  { using namespace Eigen; EIGEN_STATIC_ASSERT_LVALUE(Der1); }
  typedef typename Eigen::ArrayBase<Der1>::Scalar Scalar;

  const Eigen::Index ds = theta.rows();
  const Eigen::Index n = ds + 1;

  tomographer_assert(theta.cols() == 1);
  tomographer_assert(ddxddtheta.rows() == (int)n);
  tomographer_assert(ddxddtheta.cols() == (int)(ds*ds));

  Eigen::Array<Scalar, Der2::RowsAtCompileTime, 1> sintheta(theta.rows());
  sintheta = theta.array().sin();
  Eigen::Array<Scalar, Der2::RowsAtCompileTime, 1> costheta(theta.rows());
  costheta = theta.array().cos();

  Eigen::Index i, j, k, mm;
  for (k = 0; k < n; ++k) {
    for (i = 0; i < ds; ++i) {
      for (j = 0; j <= i; ++j) {
        //std::cout << "k,i,j = "<< k<< ", "<< i << ", " << j << "\n";
        Scalar val = 0.0;
        if (i > k) {
          val = 0.0;
        } else if (i == k && k < n-1) {
          val = -costheta(j);
          for (mm = 0; mm <= i; ++mm) {
            if (mm == j) {
              continue;
            }
            val *= sintheta(mm);
          }
        } else if (i == k && k == n-1) {
          if (j == i) {
            val = -1.0;
            for (mm = 0; mm <= i; ++mm) {
              val *= sintheta(mm);
            }
          } else { // j < i
            val = costheta(i)*costheta(j);
            for (mm = 0; mm < i; ++mm) { // i not included
              if (mm == j) {
                continue;
              }
              val *= sintheta(mm);
            }
          }
        } else if (i < k && k < n-1) {
          if (j == i) {
            val = -costheta(k);
            for (mm = 0; mm < k; ++mm) {
              val *= sintheta(mm);
            }
          } else { // j < i
            val = costheta(j)*costheta(i)*costheta(k);
            for (mm = 0; mm < k; ++mm) {
              if (mm == j || mm == i) {
                continue;
              }
              val *= sintheta(mm);
            }
          }
        } else { // i < k && k == n-1
          if (j == i) {
            val = -1.0;
            for (mm = 0; mm < n-1; ++mm) {
              val *= sintheta(mm);
            }
          } else { // j < i
            val = costheta(i)*costheta(j);
            for (mm = 0; mm < n-1; ++mm) {
              if (mm == i || mm == j) {
                continue;
              }
              val *= sintheta(mm);
            }
          }
        }
        // set the calculated value
        ddxddtheta(k, i+ds*j) = val;
        ddxddtheta(k, j+ds*i) = val; // symmetric
      }
    }
  }
}


} // namespace SphCoords
} // namespace MathTools
} // namespace Tomographer


#endif