FC2ブログ
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

# 楕円体の体積 Volume of a Ellipsoid
# 引数 半軸 半軸 半軸 ($SemiaxisA, $SemiaxisB, $SemiaxisC)
# 戻り値 楕円体の体積 (@Ellipsoid)
sub ELLIPSOID{
my ($SemiaxisA, $SemiaxisB, $SemiaxisC) = @_;
my @Ellipsoid = ();
my $Pi = atan2(1, 1) * 4;
my $X = 0;
my $K = 0;

# 半軸の確認
if(($SemiaxisA <= 0) || ($SemiaxisB <= 0) || ($SemiaxisC <= 0) || ($SemiaxisA < $SemiaxisB) || ($SemiaxisA < $SemiaxisC) || ($SemiaxisB < $SemiaxisC)){
return 0;
}

# 第一種不完全楕円積分 第二種不完全楕円積分で使用する値
$X = sqrt(1 - (($SemiaxisC * $SemiaxisC) / ($SemiaxisA * $SemiaxisA)));
$K = sqrt(1 - (($SemiaxisC * $SemiaxisC) / ($SemiaxisB * $SemiaxisB))) / $X;

# 楕円体の体積 Volume of a Ellipsoid
# [0] 楕円体の体積 [1] 楕円体の表面積
$Ellipsoid [0] = (4 * $Pi * $SemiaxisA * $SemiaxisB * $SemiaxisC) / 3;
$Ellipsoid [1] = (2 * $Pi) * (($SemiaxisC * $SemiaxisC) + ($SemiaxisA * $SemiaxisB * $X * &COMPLETEELLIPTICINTEGRALSECOND($X, $K)) + ((($SemiaxisB * ($SemiaxisC * $SemiaxisC)) / ($SemiaxisA * $X)) * &INCOMPLETEELLIPTICINTEGRALFIRST($X ,$K)));

return @Ellipsoid ;
}

# 第一種不完全楕円積分 Incomplete Elliptic Integral of the First Kind
# 引数 値 値 ($X, $K)
# 戻り値 第一種不完全楕円積分 ($IncompleteEllipticIntegralOfTheFirstKind)
sub INCOMPLETEELLIPTICINTEGRALFIRST{
my ($X, $K) = @_;
my $IncompleteEllipticIntegralOfTheFirstKind = 0;
my $SimpsonsRule = 0;
my $a = 0;
my $b = $X;
my $N = 1000;
my $h = ($b - $a) / $N;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;

# 値の確認
if(($X < -1) || (1 < $X) || (($K * $X) < -1) || (1 < ($K * $X))){
return 0;
}

# シンプソンの公式 Simpson's Rule
$SumA = &INTEGRANDFIRST($K, $a);
$SumB = &INTEGRANDFIRST($K, $b);

# 分割数は適当
for(my $i = 1; $i < $N; $i++){
my $T = $a + ($i * $h);
my $num = ($i % 2 == 0 ? 2 : 4);
my $tmp = &INTEGRANDFIRST($K, $T);

$Sum += $num * $tmp;
}

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 第一種不完全楕円積分 Incomplete Elliptic Integral of the First Kind
$IncompleteEllipticIntegralOfTheFirstKind = $SimpsonsRule;

return $IncompleteEllipticIntegralOfTheFirstKind;
}

# 第二種不完全楕円積分 Complete Elliptic Integral of the Second Kind
# 引数 値 値 ($X, $K)
# 戻り値 第二種不完全楕円積分 ($CompleteEllipticIntegralOfTheSecondKind)
sub COMPLETEELLIPTICINTEGRALSECOND{
my ($X, $K) = @_;
my $CompleteEllipticIntegralOfTheSecondKind = 0;
my $SimpsonsRule = 0;
my $a = 0;
my $b = $X;
my $N = 1000;
my $h = ($b - $a) / $N;
my $Sum = 0;
my $SumA = 0;
my $SumB = 0;

# 値の確認
if(($X < -1) || (1 < $X) || (($K * $X) < -1) || (1 < ($K * $X))){
return 0;
}

# シンプソンの公式 Simpson's Rule
$SumA = &INTEGRANDSECOND($K, $a);
$SumB = &INTEGRANDSECOND($K, $b);

for(my $i = 1; $i < $N; $i++){
my $T = $a + ($i * $h);
my $num = ($i % 2 == 0 ? 2 : 4);
my $tmp = &INTEGRANDSECOND($K, $T);

$Sum += $num * $tmp;
}

# シンプソンの公式 Simpson's Rule
$SimpsonsRule = ($SumA + $Sum + $SumB) * ($h / 3);

# 第二種不完全楕円積分 Complete Elliptic Integral of the Second Kind
$CompleteEllipticIntegralOfTheSecondKind = $SimpsonsRule;

return $CompleteEllipticIntegralOfTheSecondKind;
}

# 被積分関数 Integrand
# 引数 変数 変数 ($K, $T)
# 戻り値 被積分関数 ($Function)
sub INTEGRANDFIRST{
my ($K, $T) = @_;
my $Function = 0;

# 被積分関数 Integrand
$Function = 1 / sqrt((1 - ($T * $T)) * (1 - (($K * $K) * ($T * $T))));

return $Function;
}

# 被積分関数 Integrand
# 引数 変数 変数 ($K, $T)
# 戻り値 被積分関数 ($Function)
sub INTEGRANDSECOND{
my ($K, $T) = @_;
my $Function = 0;

# 被積分関数 Integrand
$Function = sqrt((1 - (($K * $K) * ($T * $T))) / (1 - ($T * $T)));

return $Function;
}


参考URL
楕円体の体積 - 高精度計算サイト
楕円体 - Wikipedia
第1種不完全楕円積分 F(x,k) - 高精度計算サイト
第2種不完全楕円積分 E(x,k) - 高精度計算サイト
オンライン コンパイラ/インタプリタ
テクニカル分析
プロフィール

Author:雨宮
Firefoxを使用しているので気づかなかったけど、IE6でソースコードを上手くコピーできない

5/3
携帯用ならIE6でもソースコードをコピーできる
携帯用

検索フォーム


あわせて読みたいブログパーツ
一寸先は闇 RSS

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。