//------------------------------------------------------------------------------
//	Jacobi の楕円関数に関連するメソッド
//		作成者：三上直樹，2018/01/20 作成，©三上直樹 2018
//------------------------------------------------------------------------------

using System;
using System.Windows.Forms;

using static System.Math;

namespace IIR_Design
{
	public static class Jacobi
	{
		// K/K' から母数（modulus）を求める
		public static double Modulus(double u)
		{
			double q = Exp(-PI*u);
			double a = 1.0;
			double b = 1.0;
			double c = 1.0;
			double d = q;
			int n;
			const int COUNT = 15;
			for (n=0; n<COUNT; n++)
			{
				a = a + 2.0*c*d;
				c = c*d*d;
				b = b + c;
				d = d*q;
				if (c < 1.0E-16) break;
			}
			if (n >= COUNT) MessageBox.Show("Modulus: 収束しません．");
			a = b/a;
			return 4.0*Sqrt(q)*a*a;
		}

		// 第１種完全楕円積分
		// 　算術幾何平均による方法を利用
		//	 山内他編：“電子計算機のための数値計算法Ⅲ”，p.258，培風館，1972 参照
		//	 1-x*x に対する値を計算
		public static double EllipC(double x)
		{
			double a = 1.0;
			double b = Sqrt(x);

			int n;
			const int COUNT = 20;
			for (n=0; n<COUNT; n++)
			{
				double at = (a + b)/2.0;
				b = Sqrt(a*b);
				a = at;
				if ((a - b)/a < 1.0E-8) break;
			}
			if (n >= COUNT)
			{
				MessageBox.Show("EllipC: 収束しません．");
				return 0;
			}
			else
				return (PI/2.0)/a;
		}

		// 第１種不完全楕円積分
		//	 ただし，kc = sqrt(1-k*k)
		public static double EllipI(double x, double kc)
		{
			double a = 1.0;
			double b = kc;
			double y = 1.0/x;
			double i = 0;

			int n;
			const int COUNT = 15;
			for (n=0; n<COUNT; n++)
			{
				double bt = a*b;
				a = a + b;
				b = 2.0*Sqrt(bt);
				y = y - bt/y;
				if (y == 0.0) y = Sqrt(bt)*1.0E-10;
				if (Abs(a - b) < a*1.0E-10) break;
				i = 2*i;
				if (y < 0.0) i++;
			}
			if (n >= COUNT)
			{
				MessageBox.Show("EllipI: 収束しません．");
				return 0;
			}
			if (y < 0.0) i++;
			return ( (Atan(a/y) + PI*i)/a );
		}

		// Jacobi の楕円関数，sn(u,kc), cn(u,kc), dn(u,kc) の計算
		//	 ただし，kc = 1-k*k
		public static bool SnCnDn(double u, double kc, out double sn, out double cn, out double dn)
		{
			sn = 0.0;
			cn = 1.0;
			dn = 1.0;
			if (u == 0.0) return true;

			double a = 1.0;
			double b = Sqrt(kc);
			const int COUNT = 16;
			double[] aa = new double[COUNT];
			double[] bb = new double[COUNT];
			int n;
			for (n=0; n<COUNT; n++)
			{
				aa[n] = a;
				bb[n] = b;
				double at = (a + b)/2.0;
				b = Sqrt(a*b);
				a = at;
				if ((a-b)/a < 1.0E-8) break;
			}
			if (n >= COUNT)
			{
				MessageBox.Show("SnCnDn: 収束しません．");
				return false;
			}

			double c = a/Tan(u*a);
			double d = 1.0;
			for (int k=n; k>=0; k--)
			{
				double e = c*c/a;
				c = c*d;
				a = aa[k];
				d = (e + bb[k])/(e + a);
			}

			sn = 1.0/Sqrt(1.0 + c*c);
			cn = sn*c;
			dn = d;
			return true;
		}
	}
}
