2018/07/16
微小角分割全周長計算64ビットだと間違った答えを出すのを修正(値の０初期化漏れ)
2018/07/07
周長計算の角度の処理修正 同じ角度を入れた場合と、開始角と終了角の値が近い場合の演算エラー修正
2018/06/12
計算方法が分り辛かったので説明を追加しました。
2018/03/31
プログラムのランデン変換時の配列の大きさを最適化し、計算精度を上げました。
楕円弧長の計算
楕円弧長の計算をします。
楕円弧長は、円と違って単純に計算できません、計算する方法としては、分割積分をするか、楕円積分になります。
楕円弧の面積は、円の面積に対して、半径比に比例するので、円での面積を求め、半径比で計算れば、簡単に求めるとが出来ますので、此処では取り上げません 楕円弧部各種計算 を参照してください。
簡単に楕円弧長を計算する方法は、微小角に分割して計算する方法です。
分割数を大きくすれば、精度は上げられますが、計算上の有効桁数により、限度があります。
楕円弧長の計算には第二種積分を使用します。
不完全積分には E(X,K) で表される積分で範囲をX=ＳＩＮ(α)で表される積分方法もあります。
この場合、αの最大値はπ/2となります。
α=sin-1x
X=1の時 π/2 です。
楕円積分は、上図の様に短軸を基準にして長軸方向の距離分積分をしたものですが、角度が π／２ を超える場合は、９０°＋ α' ９０°×２＋ α' の様に分割して計算します。
積分範囲を角度で指定する場合、角度αを単純に X=sinα
で換算すると、問題があります。
確かに、楕円長軸上の積分なので、計算は成り立ちますが、積分範囲 0から45°の積分と 0から135°の積分の値が同じになってしまっては積分公式を無視する事になります。
あくまでも、楕円積分式の積分の範囲は±π／２ 又は±90°で X=sinαとなります。
ランデン変換による場合は、９０°毎に分割しなくても計算出来る場合がありますが、計算の精度は９０°迄が高いので、90度に分割して計算したほうが高い精度の結果を得ることができます。
特に64ビットでコンパイルすると、データーが８バイトになってしまう為、有効桁数が最悪小数点以下11桁程度になってしまう事があるので注意が必要です。
実際の計算では、有効桁数が小数点以下11桁あれば、問題ないと思われますが、高精度を必要とする場合は、多倍長演算が必要です。
E(α,K)の楕円積分式で微小角分割積分をすると、誤差が大きくなります、楕円上の座標x,yを計算し距離により、積分したほうが精度が高くなります。
分割数としては、360°を100万分割すると、小数点以下８桁以上の精度を得ることができます。
楕円積分を使用する場合、楕円弧長を計算するのには使いづらいので、通常の楕円に対する角度が使用できるようにします。
例えば、楕円の開始点あるいは、終了点を角度θで入力する場合、直線の傾き角がθである楕円の中心を通る直線と楕円の交点から、楕円積分の角度α(円変換時の角度)を求めて計算します。
まず楕円と傾きθの楕円の中心を通る直線の交点を求めます。
楕円の中心を通る傾きθの直線をy=cxとするとc=tan(θ)
y=cx を 楕円式に代入して x を求めると
x=√(1/(1/a2+c2/b2))
これで上図の交点1の座標が求まります。
次に、長径aを半径とする円の座標にyの値を変換しますxの値は変わりません。
y'=y*a/b
これで交点2の座標が求まります。
x と y' の値から、楕円積分用の角度を求めます。
φ=arctan2(x,y')
問題はここからです。
長径基準の角度になっているので、短軸基準の角度にして積分計算を行う必要があります。
積分計算上角度を90°単位で分割しておこないます。
上図左の場合、X=1-cos(φ)が楕円積分の範囲です。
1 はsin(90°)で、cos(φ)は sin(α)=sin(90°-φ)
α = 90°-φ となります。
90°の楕円積分値から、α = 90°-φ の角度の積分値を差し引けば第１象限の楕円積分値となります。
この様にして、第１象限から第４象限迄の計算式を作成します。
奇数象限と、偶数象限の２つの計算式が出来ます。
積分角度 φ (0°～360°) 楕円積分値 A, A' (0～2π) 離心率 k=√(1 - a2 / b2)
I = φ÷90°
I = 0, 2 α = (I+1) x
90° - φ 奇数象限
1, 3 α = φ - I x
90° 偶数象限
A = E(90°,K) 90°積分値
A' =
E( α ,K) α 積分値
I = 0, 2
c = (I + 1) x A - A' 奇数象限
I = 1, 3
c = I x A + A' 偶数象限
cの値に楕円の長径a
を乗ずれば0°からの楕円弧長になります。
第二種楕円積分プログラム
楕円積分E(α,K)なので、角度は楕円に対する角度ではなく、積分用の角度です。(短軸基準です)
計算精度にExtended １０バイトを指定していますが、これなら３６０°の積分でもそれなりの精度が出ますが、Double ８バイトにすると、可成り精度が悪くなります。
ここのプログラムは、K=1 の時の計算は、角度 < 90°で 90°は計算プログラムの中の自然対数Ln(X)の Xの値がゼロになるので計算できません。
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)
Button1: TButton;
Edit1: TEdit;
LabeledEdit1: TLabeledEdit;
LabeledEdit2: TLabeledEdit;
Label1: TLabel;
Edit2: TEdit;
Label2: TLabel;
Edit3: TEdit;
Label3: TLabel;
procedure Button1Click(Sender: TObject);
private
{ Private 宣言 }
function second_imperfect_elliptic_integral(Q, K: Extended): Extended;
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses Math;
//-------------------------------
// 第１種第２種不完全積分
// FBnKn[0] 第１種不完全積分解
// EBnKn[0] 第２種不完全積分解
// ランデン変換
// Q 角度 deg
// K 離心率
//-------------------------------
function TForm1.second_imperfect_elliptic_integral(Q, K: Extended): Extended;
var
I, MI : integer;
Kn : array of Extended; // Kn
knd : array of Extended; // k'n
N2SKn : array of Extended; // 2/(1+Kn)
T2SKnd : array of Extended; // П2/(1+Kn')
BRad : array of Extended; // β(rad)
SinBn : array of Extended; // sinβn
FBnKn : array of Extended; // F(βn,kn)
EBnKn : array of Extended; // E(βn,kn)
LnD : Extended;
begin
// K = 0 の時は円なのでpiの角度値rad。
if K = 0 then begin
result := pi * Q / 180;
Edit2.Text := floatTostr(result);
exit;
end;
// 1 > K > 0 の時
MI := 0;
LND := 0;
setlength(kn, MI + 1); // kn
kn[0] := K;
while LND <> Kn[MI] do begin
LnD := Kn[MI];
inc(MI);
setlength(kn, MI + 1); // kn
kn[MI] := 2 * sqrt(Kn[MI - 1])/(1 + Kn[MI - 1]);
end;
Dec(MI);
setlength(kn, MI + 1); // kn
setlength(knd, MI + 1); // k'n
setlength(N2SKn, MI + 1); // 2/(1+Kn)
setlength(T2SKnd, MI + 1); // П2/(1+Kn')
setlength(BRad, MI + 1); // β(rad)
setlength(SinBn, MI + 1); // sinβn
setlength(FBnKn, MI + 1); // F(βn,kn)
setlength(EBnKn, MI + 1); // E(βn,kn)
knd[0] := 1;
for I := 1 to MI do
knd[I] := (1 - kn[I - 1]) / (1 + Kn[I - 1]);
for I := 0 to MI do
N2SKn[I] := 2 / (1 + kn[I]);
T2SKnd[MI] := N2SKn[MI];
for I := MI - 1 downto 0 do
T2SKnd[I] := T2SKnd[I + 1] * N2SKn[I];
BRad[0] := Q / 180 * pi;
for I := 1 to MI do
BRad[I] := (arcsin(kn[I - 1] * sin(Brad[I - 1])) + Brad[I - 1]) / 2;
for I := 0 to MI do
SinBn[I] := sin(Brad[I]);
LnD := LN((1 + sin(Brad[MI])) / cos(Brad[MI]));
for I := 0 to MI do
FBnKn[I] := T2SKnd[I] * LnD;
EBnKn[MI] := sin(Brad[MI]);
for I := MI - 1 downto 0 do
EBnKn[I] := (2 * EBnKn[I + 1] + 2 * knd[I + 1] * FBnKn[I + 1]) / (1 + knd[I + 1]) - kn[I] * SinBn[I];
result := EBnKn[0];
Edit2.Text := floatTostr(result);
end;
//--------------------------------------
// 角度分割計算とランデン変換による計算
//--------------------------------------
procedure TForm1.Button1Click(Sender: TObject);
var
I, N : integer;
K, Q, dq, Wq : Extended;
x1, y1, x2, y2 : Extended;
S, ds, b : Extended;
dx, dy : Extended;
rnd : Extended;
begin
val(LabeledEdit1.Text, K, I);
if I <> 0 then begin
application.MessageBox('Kの値に間違いがあります。','離心率K',0);
exit;
end;
if abs(K) > 1 then begin
application.MessageBox('Kの値は１迄です。','離心率K',0);
exit;
end;
val(LabeledEdit2.Text, Q, I);
if I <> 0 then begin
application.MessageBox('αの値に間違いがあります。','積分角α',0);
exit;
end;
if abs(Q) > 360 then begin
application.MessageBox('αの値が大きすぎます。','積分角α',0);
exit;
end;
// ランデン変換による計算
rnd := second_imperfect_elliptic_integral(Q, K);
// 角度分割による計算
b := sqrt(1 - K * K); // 半径b計算 半径aは１で省略
N := round(Q * 10000); // 分割数設定
N := abs(N);
if N = 0 then N := 1; // 零での除算防止
Q := Q / 180 * pi; // 角度radに変換
dq := Q / N; // Δα計算
S := 0; // 積分値クリア
x1 := cos(0) * b; // x1=cos()*bは短軸基準 長軸基準はx1=cos()
y1 := sin(0); // y1=sin() は短軸基準 長軸基準はy1=sin()*b
for I := 1 to N do begin // 指定角α迄積分計算
wq := I * dq;
x2 := cos(wq) * b; // x2=cos()*bは短軸基準 長軸基準はx2=cos()
y2 := sin(wq); // y2=sin() は短軸基準 長軸基準はy2=sin()*b
dx := x2 - x1;
dy := y2 - y1;
ds := sqrt(dx * dx + dy * dy);
S := S + ds;
x1 := x2;
y1 := y2;
end;
Edit1.Text := floattostr(S);
Edit3.Text := floattostr(abs(rnd)-S);
end;
end.
楕円弧長計算プログラム
楕円弧長なので、角度は楕円に対する角度です。
(長軸基準です)
楕円積分は、９０°単位で最大値π/2として分割して計算しています。
短半径bが非常に小さい場合は、Kの値が１になるので、１になってもエラーにならない様にしています。
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)
Image1: TImage;
Panel1: TPanel;
arEdit: TLabeledEdit;
brEdit: TLabeledEdit;
sttEdit: TLabeledEdit;
endEdit: TLabeledEdit;
Button1: TButton;
procedure FormCreate(Sender: TObject);
procedure Button1Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
procedure ImageClear;
function draw_ellipse: boolean;
function EllipseInput(var ar, br, sttQ, endQ: double): boolean;
procedure ellipse_calc(ar, br, q, magnification: double; var x, y: double);
procedure Linedlaw(x1, y1, x2, y2: double; Xsift, Ysift: integer; LF: boolean);
procedure divisionpartition;
function second_imperfect_elliptic_integral(Q, K: double): double;
procedure elliptic_integral;
function circle_calc(ar, br, q: double): double;
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
uses System.Math;
//-----------------------------
// 表示画像消去
//-----------------------------
procedure TForm1.ImageClear;
begin
image1.Canvas.Brush.Style := bsSolid;
image1.Canvas.Brush.Color := clWhite;
image1.Canvas.FillRect(Rect(0, 0, image1.Width, image1.Height));
end;
//--------------------------
// ２点間線引き
//--------------------------
procedure TForm1.Linedlaw(x1, y1, x2, y2: double; Xsift, Ysift: Integer; LF: boolean);
var
xp1, yp1, xp2, yp2: integer;
begin
if LF then begin
xp1 := round(x1) + Xsift;
yp1 := image1.Height - round(y1) - Ysift;
image1.Canvas.MoveTo(xp1, yp1);
end;
xp2 := round(x2) + Xsift;
yp2 := image1.Height - round(y2) - Ysift;
image1.Canvas.LineTo(xp2, yp2);
end;
//--------------------
// 楕円の値入力処理
//--------------------
function TForm1.EllipseInput(var ar, br, sttQ, endQ: double): boolean;
var
ch : integer;
begin
result := false;
val(arEdit.Text, ar, ch);
if ch <> 0 then begin
application.MessageBox('半径aは数値ではありません。', '半径a', 0);
exit;
end;
if ar <= 0 then begin
application.MessageBox('半径aがゼロ 又はゼロ以下です。', '半径a', 0);
exit;
end;
val(brEdit.Text, br, ch);
if ch <> 0 then begin
application.MessageBox('半径bは数値ではありません。', '半径b', 0);
exit;
end;
if br <= 0 then begin
application.MessageBox('半径bがゼロ 又はゼロ以下です。', '半径b', 0);
exit;
end;
val(sttEdit.Text, sttQ, ch);
if ch <> 0 then begin
application.MessageBox('開始角は数値ではありません。', '開始角', 0);
exit;
end;
val(endEdit.Text, endQ, ch);
if ch <> 0 then begin
application.MessageBox('終了角は数値ではありません。', '終了角', 0);
exit;
end;
// 角度が0～360°の範囲に無かったら修正
Repeat
if sttQ > 360 then sttQ := sttQ - 360;
if sttQ < 0 then sttQ := SttQ + 360;
if endQ > 360 then endQ := endQ - 360;
if endQ < 0 then endQ := endQ + 360;
Until (sttQ = 0) and (sttQ < 360) and (endQ > 0) and (endQ = 360);
// 入力角度に長径を基準とした角度θをを使用しています。
// 楕円積分は短径を基準にsinαで表される距離に対する積分となります。
// 入力された角度をそのまま使用する事は出来ません。
// 楕円の中心位置を通る直線 y=cx の cを c= tanθ として楕円との交点を求め
// 交点yの値を円変換時の値y'=y*a/bに変換、φ=arctan2(y',x)として
// 円上の角度に変換して、x=cos(φ)*a y=sin(φ)*b θ=arctan(x,y)のθが入力値と
// 一致するようにします。
// 楕円積分時は更に、長径基準φを短径基準αに変換して計算します。
sttQ := circle_calc(ar, br, sttQ);
endQ := circle_calc(ar, br, endQ);
// 円上の角度に変換することで
// 角度がマイナスに成ることがあるので再度修正
if endQ < 0 then endQ := endQ + 360;
if sttQ < 0 then sttQ := SttQ + 360;
result := true;
end;
//-----------------------
// 楕円作図
//-----------------------
function TForm1.draw_ellipse: boolean;
const
divideN = 2000; // 作図用分割数
draw_margin = 40; // 作図余白
var
ar, br : double;
sttQ, endQ : double;
maxr : double;
I : integer;
dq : double;
x1, y1 : double;
x2, y2 : double;
Xsift, Ysift : integer;
magnification : double;
sttN, endN : integer;
// x1,y1-x2,y2 線引き
procedure DlawCalc;
begin
ellipse_calc(ar, br, dq * I, magnification, x2, y2);
if I = sttN + 1 then Linedlaw(x1, y1, x2, y2, Xsift, Ysift, True)
else Linedlaw(x1, y1, x2, y2, Xsift, Ysift, False);
x1 := x2;
y1 := y2;
end;
begin
result := true;
if not EllipseInput(ar, br, sttQ, endQ) then begin
result := false;
exit;
end;
// 画像消去
ImageClear;
// 作図設定
Xsift := image1.Width div 2;
Ysift := image1.Height div 2;
sttQ := sttQ / 180 * pi;
endQ := endQ / 180 * pi;
// 作図計算ピッチ計算
dq := 2 * pi / divideN;
// 大きい方の半径選択
maxr := ar;
if br > ar then maxr := br;
// 作図倍率設定
magnification := (Xsift - draw_margin) / maxr;
// 楕円作図
image1.Canvas.Pen.Style := psSolid;
image1.Canvas.Pen.Color := clBlue;
image1.Canvas.Pen.Width := 1;
sttN := 0;
// 最初の角度
ellipse_calc(ar, br, 0, magnification, x1, y1);
// 最後まで
for I := 1 to divideN do DlawCalc;
// 楕円弧作図
image1.Canvas.Pen.Color := clRed;
image1.Canvas.Pen.Width := 3;
sttN := round(sttQ / dq);
endN := trunc(endQ / dq);
// 最初の角度
ellipse_calc(ar, br, sttQ, magnification, x1, y1);
// 開始角より終了角が大きい場合
if endQ >= sttQ then begin
ellipse_calc(ar, br, sttN * dq, magnification, x1, y1);
for I := sttN + 1 to endN do DlawCalc;
end
// 開始角より終了角が小さい場合
else begin
// 開始角より360°迄作図
for I := sttN + 1 to divideN do DlawCalc;
// 0°(=360°)から終了角迄作図
for I := 1 to endN do DlawCalc;
end;
// 最後の角度
ellipse_calc(ar, br, endQ, magnification, x2, y2);
Linedlaw(x1, y1, x2, y2, Xsift, Ysift, True);
// 中心線作図
image1.Canvas.Pen.Width := 1;
image1.Canvas.Pen.Style := psDashDot;
image1.Canvas.Pen.Color := clGreen;
Linedlaw(-50, 0, 50, 0, Xsift, Ysift, True);
Linedlaw( 0, -50, 0, 50, Xsift, Ysift, True);
end;
//-------------------------------------------------------------
// 円変換時の 角度計算
// 積分角度 q 角度
// ar , br 半径
// x^2/a^2 + y^2/b^2 = 1 と 角度θの直線y = cx の連立方程式の
// の解法による円変換時の角度の計算
//-------------------------------------------------------------
function TForm1.circle_calc(ar, br, q: double): double;
var
c, ga, sq, cb: double;
x, y : double;
begin
if (q = 0) or (q = 90) or (q = 180) or (q = 270) or (q = 360) then begin
result := q;
exit;
end;
// 楕円と直線の交点の計算
sq := pi / 180 * q;
c := tan(sq); // y = cx 傾斜c計算
cb := c * c / br / br; // y^2/b^2 のyに代入
ga := 1 / ar / ar; //
x := sqrt(1 / (ga + cb)); // xの値
// xの値がプラスしか計算されないので角度により符号設定
if (q > 90) and (q < 270) then x := - abs(x);
y := c * x * ar / br; // 円上のyの値
// 角度に変換
result := arctan2(y, x) / pi * 180;
end;
//-------------------------------------------------------------
// 楕円上の座標計算
// q 角度
// ar , br 半径
// magnification 倍率
// x, y 座標
// q 角度は中心に対するxy座標の角度ではないので注意が必要です。
//-------------------------------------------------------------
procedure TForm1.ellipse_calc(ar, br, q, magnification: double; var x, y: double);
begin
x := ar * cos(q) * magnification;
y := br * sin(q) * magnification;
end;
//---------------------------------
// 周長 楕円弧長の分割計算
//---------------------------------
procedure TForm1.divisionpartition;
const
divideN = 1440000; // 全周計算分割数
divide4 = 4000; // 指定範囲分割数基礎数
PI2 = 2 * pi; // 2π
K180 = 180; // 180°
K360 = 360; // 360°
var
ar, br : double; // 半径a , b
sttQ, endQ : double; // 積分開始角、積分終了角 deg単位
I : integer; //
dq : double; // 微小角Δq
x1, y1 : double; // 始点座標
x2, y2 : double; // 終点座標
sttN : integer; // 分割数
Lng, All : double; // 楕円弧長 周長
Q : double; // 座標計算角度
basQ, topQ : double; // 積分開始角、積分終了角 rad単位
// 座標の計算と微小角に対する距離の計算と加算
procedure LngCalc;
var
xh2, yh2 : double;
begin
Q := basQ + dq * I; // 座標計算角度
if Q >= PI2 then Q := Q - PI2;
ellipse_calc(ar, br, Q, 1, x2, y2); // 座標計算
xh2 := x2 - x1;
yh2 := y2 - y1;
Lng := Lng + sqrt(xh2 * xh2 + yh2 * yh2);
x1 := x2;
y1 := y2;
end;
begin
// 各値の入力値変換
if not EllipseInput(ar, br, sttQ, endQ) then exit;
// 全周計算ピッチ計算
basQ := 0; // 2018/07/16追加
dq := PI2 / divideN;
Lng := 0;
// 全周長計算
ellipse_calc(ar, br, 0, 1, x1, y1);
for I := 1 to divideN do LngCalc;
All := Lng; // 全周長
// 指定範囲計算
basQ := sttQ / K180 * pi;
topQ := endQ / K180 * pi;
// 分割数と微小角度Δqの計算
if endQ >= sttQ then begin
sttN := round((endQ - sttq) * divide4); // 分割数
if sttN = 0 then sttn := 1;
dq := (topQ - basQ) / sttN; // Δq
end
else begin
sttN := round((endQ + K360 - sttQ) * divide4); // 分割数
if sttN = 0 then sttn := 1;
dq := (topQ + PI2 - basQ) / sttN; // Δq
end;
Lng := 0;
// 最初の座標
ellipse_calc(ar, br, basQ, 1, x1, y1);
for I := 1 to sttn do begin
LngCalc;
end;
Image1.Canvas.TextOut(20, 2,'分割計算 楕円弧長= ' + floatTostr(Lng));
Image1.Canvas.TextOut(20,17,'分割計算 周 長 = ' + floatTostr(ALL));
end;
//-------------------------------
// 第１種第２種楕円積分
// FBnKn[0] 第１種積分解
// EBnKn[0] 第２種積分解
// ランデン変換
// Q 角度 deg
// K 離心率
//-------------------------------
function TForm1.second_imperfect_elliptic_integral(Q, K: double): double;
var
I, MI : integer;
Kn : array of double; // Kn
knd : array of double; // k'n
N2SKn : array of double; // 2/(1+Kn)
T2SKnd : array of double; // П2/(1+Kn')
BRad : array of double; // β(rad)
SinBn : array of double; // sinβn
FBnKn : array of double; // F(βn,kn)
EBnKn : array of double; // E(βn,kn)
LnD : double;
begin
// K = 0 の時は円なのでpiの角度値radをかえします。
if K = 0 then begin
result := pi * Q / 180;
exit;
end;
// この時は１をかえします。
if (K = 1) and (Q = 90) then begin
result := 1;
exit;
end;
// (1 > K > 0 and Q <= 90) or (K = 1 and Q < 90)の時
MI := 0;
LND := 0;
setlength(kn, MI + 1); // kn
kn[0] := K;
while LND <> Kn[MI] do begin
LnD := Kn[MI];
inc(MI);
setlength(kn, MI + 1); // kn
kn[MI] := 2 * sqrt(Kn[MI - 1])/(1 + Kn[MI - 1]);
end;
Dec(MI);
setlength(kn, MI + 1); // kn
setlength(knd, MI + 1); // k'n
setlength(N2SKn, MI + 1); // 2/(1+Kn)
setlength(T2SKnd, MI + 1); // П2/(1+Kn')
setlength(BRad, MI + 1); // β(rad)
setlength(SinBn, MI + 1); // sinβn
setlength(FBnKn, MI + 1); // F(βn,kn)
setlength(EBnKn, MI + 1); // E(βn,kn)
kn[0] := K;
for I := 1 to MI do
kn[I] := 2 * sqrt(Kn[I - 1])/(1 + Kn[I - 1]);
knd[0] := 1;
for I := 1 to MI do
knd[I] := (1 - kn[I - 1]) / (1 + Kn[I - 1]);
for I := 0 to MI do
N2SKn[I] := 2 / (1 + kn[I]);
T2SKnd[MI] := N2SKn[MI];
for I := MI - 1 downto 0 do
T2SKnd[I] := T2SKnd[I + 1] * N2SKn[I];
BRad[0] := Q / 180 * pi;
for I := 1 to MI do
BRad[I] := (arcsin(kn[I - 1] * sin(Brad[I - 1])) + Brad[I - 1]) / 2;
for I := 0 to MI do
SinBn[I] := sin(Brad[I]);
LnD := LN((1 + sin(Brad[MI])) / cos(Brad[MI]));
for I := 0 to MI do
FBnKn[I] := T2SKnd[I] * LnD;
EBnKn[MI] := sin(Brad[MI]);
for I := MI - 1 downto 0 do
EBnKn[I] := (2 * EBnKn[I + 1] + 2 * knd[I + 1] * FBnKn[I + 1]) / (1 + knd[I + 1]) - kn[I] * SinBn[I];
result := EBnKn[0];
end;
//--------------------------------------------------
// 楕円積分による周長計算
// ランデン変換による計算をしていますが90°位までが
// 精度が良いので90°に分けて計算しています。
// 楕円積分は短径が基準となっていますが、角度は長径が
// 基準となっているので楕円弧長計算時変換しています。
//--------------------------------------------------
procedure TForm1.elliptic_integral;
var
ar, br : double;
sttQ, endQ : double;
K, A, L : double;
sttL, endL : double;
IND : integer;
// 角度による楕円弧長計算
function Elliptical_arc(q : double; I: integer): double;
var
sttendR : double;
begin
result := 0;
case I of
0, 2 : begin
sttendR := (I + 1) * 90 - q;
result := (I + 1) * A - ar * second_imperfect_elliptic_integral(sttendR, K);
end;
1, 3 : begin
sttendR := q - (90 * I);
result := I * A + ar * second_imperfect_elliptic_integral(sttendR, K);
end;
end;
end;
begin
// 各値の入力値変換
if not EllipseInput(ar, br, sttQ, endQ) then exit
if (sttQ = 360) and (endQ <> 360) then sttQ := 0;
// 半径Brの方が大きかったら入れ替え 角度を90°変更
if br > ar then begin
k := ar;
ar := br;
br := k;
if (sttq <> 0) or (endQ <> 360) then begin
sttQ := sttQ + 90;
endQ := endQ + 90;
if sttQ >= 360 then sttQ := sttQ - 360;
if endQ >= 360 then endQ := endQ - 360;
end;
end;
sttL := 0;
endL := 0;
K := sqrt(1 - br * br / ar / ar); // k 離心率
// 角度90°円弧長計算
A := ar * second_imperfect_elliptic_integral(90, K); // 半径 x 第２種不完全積分
// 範囲計算
for IND := 3 downto 0 do begin
// 開始円弧長計算
if (IND * 90 <= sttQ) and (sttQ < (IND + 1) * 90) then
sttL := Elliptical_arc(sttQ, IND);
// 終了円弧長計算
if (IND * 90 < endQ) and (endQ <= (IND + 1) * 90) then
endL := Elliptical_arc(endQ, IND);
end;
// 開始角と終了角の大きさで計算方法変更
if endQ >= sttQ then L := endL - sttL
else L := 4 * A - sttL + endL;
if endQ = sttQ then L := 0;
Image1.Canvas.TextOut(20, image1.Height - 32,'不完全積分 楕円弧長= ' + floatTostr(L));
Image1.Canvas.TextOut(20, image1.Height - 17,'不完全積分 周 長 = ' + floatTostr(4 * A));
end;
//----------------------
// 周長 楕円弧長計算
//----------------------
procedure TForm1.Button1Click(Sender: TObject);
begin
if not draw_ellipse then exit; // 作図
divisionpartition; // 分割計算
elliptic_integral; // 不完全楕円積分計算
end;
//------------
// 初期設定
//------------
procedure TForm1.FormCreate(Sender: TObject);
begin
ClientHeight := 329;
ClientWidth := 489;
Image1.Top := 0;
Image1.Left := 0;
Image1.Height := ClientHeight;
Image1.Width := ClientHeight;
Image1.Canvas.Font.Height := 15;
Panel1.Top := 8;
Panel1.Left := 352;
Panel1.Height := 313;
Panel1.Width := 113;
end;
end.
EllipsePerimeter.zip