ヤコビの楕円関数 ｓｎ　ｃｎ　ｄｎ

　ヤコビの楕円関数　ｓｎ　ｃｎ　ｄｎの計算プログラムです。
Mathematics Source Library　C & ASM　にあるプログラムをＤｅｌｐｈｉ用に変換したものです。
ｓｎ，ｃｎ，ｄｎの詳細については ヤコビの楕円関数を参照してください。

　離心率（ｋ）は１≧ｋ≧－１の範囲で指定します。
パラメタ（ｍ）は１≧ｍ≧０の範囲で指定し、角度(α)は±π/2の範囲です。
ｓｎさえ計算出来れば、ｃｎ，ｄｎ は ｓｎ から簡単に変換が可能です。

プログラム

//------------------------------
//  ヤコビの楕円関数 sn cn dn
//------------------------------
unit Main;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
  Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.ExtCtrls;

type
  TForm1 = class(TForm)
    Memo1: TMemo;
    snukbtn: TButton;
    u_Edit: TLabeledEdit;
    K_Edit: TLabeledEdit;
    kRadioButton: TRadioButton;
    mRadioButton: TRadioButton;
    aRadioButton: TRadioButton;
    m_Edit: TLabeledEdit;
    a_Edit: TLabeledEdit;
    procedure snukbtnClick(Sender: TObject);
  private
    { Private 宣言 }
    function jacobi_sn(u: double; arg: char; x: double): double;
    function Jacobi_cn(u: double; arg: char; x: double): double;
    function Jacobi_dn(u: double; arg: char; x: double): double;
    function Jacobi_am(u: double; arg: char; x: double): double;
  public
    { Public 宣言 }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses Math;


// 振幅関数
function TForm1.Jacobi_am(u: double; arg: char; x: double): double;
const
  KN = 30;
  EPSILON = 1E-14;
var
  a: array of double;
  g: array of double;
  c: array of double;
  two_n : double;
  phi   : double;
  k     : double;
  i, j, N : integer;
  half    : double;
begin
  if x = 0 then begin
    result := u;
    exit;
  end;
  case arg of
    'a' : k := sin(abs(x));
    'm' : k := sqrt(abs(x));
    else k := abs(x);
  end;
  if k = 1 then begin
    result := 2 * arctan(exp(u)) - pi / 2;
    exit;
  end;
  half := 1 / 2;
  N := 1;
  setlength(a, N);
  setlength(g, N);
  setlength(c, N);
  a[0] := 1;
  g[0] := sqrt(1 - k * k);
  c[0] := k;
  two_n := 1;
  for i := 0 to KN do begin
    if abs(a[i] - g[i]) < a[i] * EPSILON then break;
//    if a[i] - g[i] = 0 then break;
    two_n := two_n + two_n;
    inc(N);
    setlength(a, N);
    setlength(g, N);
    setlength(c, N);
    a[i + 1] := half * (a[i] + g[i]);
    g[i + 1] := sqrt(a[i] * g[i]);
    c[i + 1] := half * (a[i] - g[i]);
  end;
  phi := two_n * a[i] * u;
  for j := i downto 1 do
    phi := half * (phi + arcsin(c[j] * sin(phi) / a[j]));
  result := phi;
  memo1.Lines.Add('Landen Loop='+ intTostr(N - 1));
end;

// ヤコビ楕円関数
function TForm1.jacobi_sn(u: double; arg: char; x: double): double;
begin
  result := sin(Jacobi_am(u, arg, x));
end;

function TForm1.Jacobi_cn(u: double; arg: char; x: double): double;
begin
  result := cos(Jacobi_am(u, arg, x) );
//  result := sqrt(1 - jacobi_sn(u, arg, x) * jacobi_sn(u, arg, x));
end;

function TForm1.Jacobi_dn(u: double; arg: char; x: double): double;
var
  sn : double;
begin
  sn := sin(Jacobi_am(u, arg, x));
  case arg of
      'm' : begin
              result :=  sqrt(1.0 - x * sn * sn);
              exit;
            end;
      'a' : x := sin(x);
  end;
  x := x * sn;
  result := sqrt(1.0 - x * x );
end;

// sn(u,k),cn(u,k),dn(u,k)計算
procedure TForm1.snukbtnClick(Sender: TObject);
var
  x, u, snuk  : double;
  cnuk, dnuk  : double;
  ch          : integer;
  arg         : char;
begin
  x := 0;
  if kradiobutton.Checked then val(k_edit.Text, x, ch);
  if mradiobutton.Checked then val(m_edit.Text, x, ch);
  if aradiobutton.Checked then val(a_edit.Text, x, ch);
  if ch <> 0 then begin
    application.MessageBox('parameterに間違いがあります。','注意',0);
    exit;
  end;
  val(u_edit.Text, u, ch);
  if ch <> 0 then begin
    application.MessageBox('値(u)に間違いがあります。','注意',0);
    exit;
  end;
  arg := '0';
  if kradiobutton.Checked then begin
    if abs(x) > 1 then begin
      application.MessageBox('離心率の最大値は１です。','注意',0);
      exit;
    end;
    arg := 'k';
  end;
  if mradiobutton.Checked then begin
    if x > 1 then begin
      application.MessageBox('ｍの最大値は１です。','注意',0);
      exit;
    end;
    if x < 0 then begin
      application.MessageBox('ｍの最小値は０です。','注意',0);
      exit;
    end;
    arg := 'm';
  end;
  if aradiobutton.Checked then begin
    if abs(x) > pi / 2 then begin
      application.MessageBox(pchar('αの最大値は' + floattostrF(pi / 2, fffixed, 8,6) + 'です。'),'注意',0);
      exit;
    end;
    arg := 'a';
  end;
  memo1.Clear;
  snuk := jacobi_sn(u, arg, x);
  cnuk := jacobi_cn(u, arg, x);
  dnuk := jacobi_dn(u, arg, x);
  if arg = 'a' then arg := 'α';
  memo1.Lines.Add('sn(u,' + arg + ')=' + floatTostr(snuk));
  memo1.Lines.Add('');
  memo1.Lines.Add('cn(u,' + arg + ')=' + floatTostr(cnuk));
  memo1.Lines.Add('');
  memo1.Lines.Add('dn(u,' + arg + ')=' + floatTostr(dnuk));
end;

end.

　　　　　jacobi_elliptic_sn_cn_dn_functions.zip

各種プログラム計算例に戻る

　　　　　　最初に戻る