////////////////////////
// Evaluate Jacobi elliptic functions sn cn dn
// by means library routines.

// Excellent reference for the math and concepts:
//  http://www.und.edu/instruct/schwalm/MAA_Presentation_10-02/handout.pdf

#include <iostream>             /* for cout */
#include <iomanip>              /* for setw() */
#include <boost/math/special_functions/jacobi_elliptic.hpp>
using namespace std;

int main(){
  int halfpts = 200;    // number of points in each half of x-axis
  double xmax = 10;     // plot abscissa from -xmax to +xmax
  double kk = 0.95;     // the modulus

// column headers:
  cout << setw(10) << "x"
       << setw(15) << "sn(x)"
       << setw(15) << "cn(x)"
       << setw(15) << "dn(x)"
       << setw(15) << "check"
       << endl;

// main loop
// Note total number of points is 1+2*halfpts
  for (int ii = -halfpts; ii <= halfpts; ii++){

// Write this as double(halfpts)
// to make it absolutely clear to the compiler
// that we want floating-point math, not integer math
    double xx;          // the argument
    xx = ii * xmax / double(halfpts);

    double sn, cn, dn;

// The library routine calculates all three (sn, cn, dn)
// but returns the answers in an inelegant way:
// the sn result gets returned as the function value
// while cn and dn get returned via pointers in the argument list:
    sn = boost::math::jacobi_elliptic(kk, xx, &cn, &dn);

    cout << setw(10) << xx
         << setw(15) << sn
         << setw(15) << cn
         << setw(15) << dn
         << setw(15) << sn*sn + cn*cn - 1.
         << endl;
  }
}