A..1 easy_elliptic.hpp

/*
 * easy_elliptic.hpp --- boost を簡単に使って、K(k),sn(u;k),cn(u;k),dn(u;k)計算
 *  よく考えてみたら、agm.cpp の内容を移せば、boost に依存しないように出来る。
 *  つまり全く自前の計算で済ませられる。これは驚いた。
 */

#include <boost/math/special_functions.hpp>

static double K(double k)
{
  return boost::math::ellint_1(k);
}

static double F(double k, double phi)
{
  return boost::math::ellint_1(k, phi);
}

static double sn(double k, double x)
{
  return boost::math::jacobi_sn(k, x);
}

static double cn(double k, double x)
{
  return boost::math::jacobi_cn(k, x);
}

static double dn(double k, double x)
{
  return boost::math::jacobi_dn(k, x);
}

// 複素関数として sn()
static std::complex<double> jacobi_sn(double k, std::complex<double> z)
{
  double kp,m,x,y,s,c,d,s1,c1,d1,den;
  x = z.real(); y = z.imag(); kp = sqrt(1 - k * k); m = k * k;
#ifdef OLD
  s  = boost::math::jacobi_sn(k,  x);
  c  = boost::math::jacobi_cn(k,  x);
  d  = boost::math::jacobi_dn(k,  x);
  s1 = boost::math::jacobi_sn(kp, y);
  c1 = boost::math::jacobi_cn(kp, y);
  d1 = boost::math::jacobi_dn(kp, y);
#else
  s =  boost::math::jacobi_elliptic(k, x, &c, &d);
  s1 = boost::math::jacobi_elliptic(kp, y, &c1, &d1);
#endif
  den = c1 * c1 + m * s * s * s1 * s1; // denominator: 分母
  return std::complex<double>(s * d1 / den, c * d * s1 * c1 / den);
}

// 複素数版 Jacobi の sn(), cn(), dn()
static std::complex<double> jacobi_elliptic(double k,
					    std::complex<double> z,
					    std::complex<double> *cn,
					    std::complex<double> *dn)
{
  double kp,m,x,y,s,c,d,s1,c1,d1,den;
  std::complex<double> sn;
  x = z.real(); y = z.imag(); kp = sqrt(1 - k * k); m = k * k;
#ifdef OLD
  s  = boost::math::jacobi_sn(k, x);
  c  = boost::math::jacobi_cn(k, x);
  d  = boost::math::jacobi_dn(k, x);
  s1 = boost::math::jacobi_sn(kp, y);
  c1 = boost::math::jacobi_cn(kp, y);
  d1 = boost::math::jacobi_dn(kp, y);
#else
  s =  boost::math::jacobi_elliptic(k, x, &c, &d);
  s1 = boost::math::jacobi_elliptic(kp, y, &c1, &d1);
#endif  
  den = c1 * c1 + m * s * s * s1 * s1; // denominator: 分母
  sn  = std::complex<double>(s * d1 / den, c * d * s1 * c1 / den);
  *cn = std::complex<double>(c * c1 / den, - s * d * s1 * d1 / den);
  *dn = std::complex<double>(d * c1 * d1 / den, - m * s * c * s1 / den);
  return sn;
}

// まだ十分テストしていないけれど
static double my_jacobi_elliptic(double k, double u, double *cn, double *dn)
{
  int n, maxn = 10, N, success;
  double kp, a[maxn+1], b[maxn+1], c[maxn+1], phi[maxn+1], eps, myK, pi;

  success = 0;
  pi = 4.0 * atan(1.0);
  // AGM
  eps = 1e-15;
  kp = sqrt(1 - k * k);
  a[0] = 1.0; b[0] = kp; c[0] = k;
  for (n = 1; n <= maxn; n++) {
    a[n] = (a[n-1] + b[n-1]) / 2;
    b[n] = sqrt(a[n-1] * b[n-1]);
    c[n] = (a[n-1] - b[n-1]) / 2;
    if (c[n] < eps * a[n]) {
      success = 1;
      break;
    }
  }
  if (!success) {
    std::cerr << "収束しませんでした。" << std::endl;
  }
  N = n;
  // 完全楕円積分
  myK = pi / (2 * a[N]); // これは返さない...
  // 楕円関数の amplitude の計算
  phi[N] = (1 << n) * a[N] * u; // φN
  for (n = N; n >= 1; n--)
    phi[n-1] = (phi[n] + asin(c[n] * sin(phi[n]) / a[n])) / 2;
  //
  *cn = cos(phi[0]);
  *dn = cos(phi[0]) / cos(phi[1]-phi[0]);
  return sin(phi[0]);
}