外部変形は データのやり取りをテキストファイルで行うので プログラム言語は 自由に選ぶことができます。図形は機能的かつシンプルなため、数多くのユーザーに受け入れられています。
rubyで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する
その① 台形則
:rubyで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する 台形則
@echo off
ruby -x %0 jwc_temp.txt
goto:eof
REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e
#!ruby -Ks -an -i.bak
def rad(x) #度をラジアンに変換
x ? x * PI / 180.0 : 0
end
def apa(o, q) #角度補正
return (90.0 * (2 * ( (o + 1) / 2).to_i - 1) - q) * (-1) ** (o - 1)
end
def arc_l(r, p1, p2, w) #楕円弧の長さ
#r, p1, p2, w 半径、始角 ﾟ、終角 ﾟ、扁平率
if w > 1
w = 1.0 / w
r /= w
p1 -= 90
p2 -= 90
end
k = sqrt(1.0 - w * w) #離心率
ek = e(k) #第二種完全楕円積分
p1 += 360 while p1 < 0
p1 -= 360 while p1 >= 360
p2 += 360 while p2 < 0
p2 -= 360 while p2 >= 360
return r * 4.0 * ek if p1 == p2
n1 = 1 + (p1 / 90.0).to_i
n1 = 4 if n1 > 4
lp = ek * ( (n1 + 1) % 2) + e(rad(apa(n1, p1)), k) * (-1.0) ** (n1 - 1)
n2 = 1 + (p2 / 90.0).to_i
n2 = 4 if n2 > 4
n1 += 4 if n1 == n2
n3 = 5 + 2 * ( (n2 - 1) / 2.0).to_i
lp += ek * ( (n3 - n1) % 4) + e(rad(apa(n2, p2)), k) * (-1.0) ** n2
lp += 4.0 * ek if lp < 0.0
return r * lp
end
def e(p, k = nil) #台形則 : 第二種不完全楕円積分
def fx(x, k)
return sqrt(1.0 - (k * sin(x)) ** 2)
end
if k == nil #第二種完全楕円積分の補正
k = p
p = PI / 2
end
n = 10000
h = p / n
s1 = s2 = 0.0
for i in 0..n-1
s1 += fx(x = h * i, k)
s2 += fx(x + h, k)
end
return h / 2.0 * (s1 + s2)
end
BEGIN { include Math }
case $_
when /^hs/
hs = {}
for i in 0..15
hs[ sprintf("%x", i) ] = $F[i + 1].to_f
end
when /^lg[0-9a-f]/
plg = $_[2, 1].to_s
$lg = plg if $lg == nil
when /^ci/
p = hs[$lg] / hs[plg]
x, y, r = $F[1..3].map { |z| z.to_f * p }
p1, p2, w, d = 0, 0, 1, 0
if $F.size == 8
p1, p2, w, d = $F[4..7].map { |z| z.to_f }
end
if w == 1
lr = r * rad(p1 == p2 ? 360 : p2 - p1)
else
lr = arc_l(r, p1, p2, w)
end
puts "h#円弧長 L = %.3f" % lr
end
__END__
その② Simpson則
:rubyで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する Simpson
@echo off
ruby -x %0 jwc_temp.txt
goto:eof
REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e
#!ruby -Ks -an -i.bak
def rad(x) #度をラジアンに変換
x ? x * PI / 180.0 : 0
end
def apa(o, q) #角度補正
return (90.0 * (2 * ( (o + 1) / 2).to_i - 1) - q) * (-1) ** (o - 1)
end
def arc_l(r, p1, p2, w) #楕円弧の長さ
#r, p1, p2, w 半径、始角 ﾟ、終角 ﾟ、扁平率
if w > 1
w = 1.0 / w
r /= w
p1 -= 90
p2 -= 90
end
k = sqrt(1.0 - w * w) #離心率
ek = e(k) #第二種完全楕円積分
p1 += 360 while p1 < 0
p1 -= 360 while p1 >= 360
p2 += 360 while p2 < 0
p2 -= 360 while p2 >= 360
return r * 4.0 * ek if p1 == p2
n1 = 1 + (p1 / 90.0).to_i
n1 = 4 if n1 > 4
lp = ek * ( (n1 + 1) % 2) + e(rad(apa(n1, p1)), k) * (-1.0) ** (n1 - 1)
n2 = 1 + (p2 / 90.0).to_i
n2 = 4 if n2 > 4
n1 += 4 if n1 == n2
n3 = 5 + 2 * ( (n2 - 1) / 2.0).to_i
lp += ek * ( (n3 - n1) % 4) + e(rad(apa(n2, p2)), k) * (-1.0) ** n2
lp += 4.0 * ek if lp < 0.0
return r * lp
end
def e(p, k = nil) #Simpson : 第二種不完全楕円積分
def fx(x, k)
return sqrt(1.0 - (k * sin(x)) ** 2)
end
if k == nil #第二種完全楕円積分の補正
k = p
p = PI / 2
end
n = 1000
x = h = p / (n * 2.0)
s1 = 0
s2 = fx(x, k)
for i in 1..n-1
s1 += fx(x += h, k)
s2 += fx(x += h, k)
end
return h / 3.0 * (fx(0, k) + fx(p, k) + 4.0 * s2 + 2.0 * s1)
end
BEGIN { include Math }
case $_
when /^hs/
hs = {}
for i in 0..15
hs[ sprintf("%x", i) ] = $F[i + 1].to_f
end
when /^lg[0-9a-f]/
plg = $_[2, 1].to_s
$lg = plg if $lg == nil
when /^ci/
p = hs[$lg] / hs[plg]
x, y, r = $F[1..3].map { |z| z.to_f * p }
p1, p2, w, d = 0, 0, 1, 0
if $F.size == 8
p1, p2, w, d = $F[4..7].map { |z| z.to_f }
end
if w == 1
lr = r * rad(p1 == p2 ? 360 : p2 - p1)
else
lr = arc_l(r, p1, p2, w)
end
puts "h#円弧長 L = %.3f" % lr
end
__END__
その③ Gauss-Legendre
:rubyで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する Gauss-Legendre
@echo off
ruby -x %0 jwc_temp.txt
goto:eof
REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e
#!ruby -Ks -an -i.bak
def rad(x) #度をラジアンに変換
x ? x * PI / 180.0 : 0
end
def apa(o, q) #角度補正
return (90.0 * (2 * ( (o + 1) / 2).to_i - 1) - q) * (-1) ** (o - 1)
end
def arc_l(r, p1, p2, w) #楕円弧の長さ
#r, p1, p2, w 半径、始角 ﾟ、終角 ﾟ、扁平率
if w > 1
w = 1.0 / w
r /= w
p1 -= 90
p2 -= 90
end
k = sqrt(1.0 - w * w) #離心率
ek = e(k) #第二種完全楕円積分
p1 += 360 while p1 < 0
p1 -= 360 while p1 >= 360
p2 += 360 while p2 < 0
p2 -= 360 while p2 >= 360
return r * 4.0 * ek if p1 == p2
n1 = 1 + (p1 / 90.0).to_i
n1 = 4 if n1 > 4
lp = ek * ( (n1 + 1) % 2) + e(rad(apa(n1, p1)), k) * (-1.0) ** (n1 - 1)
n2 = 1 + (p2 / 90.0).to_i
n2 = 4 if n2 > 4
n1 += 4 if n1 == n2
n3 = 5 + 2 * ( (n2 - 1) / 2.0).to_i
lp += ek * ( (n3 - n1) % 4) + e(rad(apa(n2, p2)), k) * (-1.0) ** n2
lp += 4.0 * ek if lp < 0.0
return r * lp
end
def e(p, k = nil) #Gauss-Legendre : 第二種不完全楕円積分
def fx(x, k)
return sqrt(1.0 - (k * sin(x)) ** 2)
end
if k == nil #第二種完全楕円積分の補正
k = p
p = PI / 2
end
l = 0
if $n > 1
for i in 1..($n / 2).to_i
l += $h[i] * fx($e[i] * p, k)
end
end
if ($n % 2) == 1
m = ($n / 2).to_i + $n % 2
l += $h[m] * fx($e[m] * p, k) / 2.0 #奇数個の場合
end
return l * p
end
def gear(bpa)
ww = bpa
if bpa > 1 then ww = 1.0 / bpa end
if ww == 0 || ww == nil then ww = 0.1 end
$n = (12 * (1.5 / ww)).to_i
$n = 300 if $n > 300
eh($n)
return bpa
end
def eh(n) # 積分点デ－タ
$e =
$h =
m = (n / 2).to_i + n % 2
e1 = n * (n + 1.0)
for i in 1..m do #DO 1
t = (4.0 * i - 1.0) * PI / (4.0 * n + 2.0)
x0 = (1.0 - (1.0 - 1.0 / n) / (8.0 * n * n)) * cos(t)
pkm1= 1.0
pk = x0
for k in 2..n do #DO 3
t1 = x0 * pk
pkp1 = t1 - pkm1 - (t1 - pkm1) / k + t1
pkm1 = pk
pk = pkp1
end #DO 3
den = 1.0 - x0 * x0
d1 = n * (pkm1 - x0 * pk)
dpn = d1 / den
d2pn = (2.0 * x0 * dpn - e1 * pk) / den
d3pn = (4.0 * x0 * d2pn + (2.0 - e1) * dpn) / den
d4pn = (6.0 * x0 * d3pn + (6.0 - e1) * d2pn) / den
u = pk / dpn
v = d2pn / dpn
h = -1.0 * u * (1.0 + 0.5 * u * (v + u * (v * v - d3pn / (3.0 * dpn))))
pp = pk + h * (dpn + 0.5 * h * (d2pn + h / 3.0 * (d3pn + 0.25 * h * d4pn)))
dp = dpn + h * (d2pn + 0.5 * h * (d3pn + h * d4pn / 3.0))
h -= pp / dp
xi = x0 + h
fx = d1 - h * e1 * (pk + 0.5 * h * (dpn + h / 3.0 * (d2pn + 0.25 * h * (d3pn + 0.2 * h * d4pn))))
wi = 2.0 * (1.0 - xi * xi) / (fx * fx)
$e[i] = xi
$h[i] = wi
end #DO 1
end
BEGIN { include Math }
case $_
when /^hs/
hs = {}
for i in 0..15
hs[ sprintf("%x", i) ] = $F[i + 1].to_f
end
when /^lg[0-9a-f]/
plg = $_[2, 1].to_s
$lg = plg if $lg == nil
when /^ci/
p = hs[$lg] / hs[plg]
x, y, r = $F[1..3].map { |z| z.to_f * p }
p1, p2, w, d = 0, 0, 1, 0
if $F.size == 8
p1, p2, w, d = $F[4..7].map { |z| z.to_f }
end
if w == 1
lr = r * rad(p1 == p2 ? 360 : p2 - p1)
else
lr = arc_l(r, p1, p2, gear(w))
end
puts "h#円弧長 L = %.3f" % lr
end
__END__
計算精度は Gauss-Legendre による楕円積分が良好です。楕円積分は数値積分のほか級数解も紹介されていましたがいまはリンクぎれのようです。
rubyに書き直したメモがあるので級数解による不完全楕円積分の例を紹介します。
@echo off
ruby -x %~f0
pause
goto:eof
#!ruby -Ks
=begin
○楕円積分 数値計算ノート(リンクは切れています)
http://www.ab.cyberhome.ne.jp/~naka-kevin/Elliptic-Integral.html
より 抜粋
楕円積分 (∫φ_0は0からφまでの定積分を表しています)
第1種楕円積分 E1(φ,k)=∫φ_0 1/√(1-k^2 sin^2φ) dφ=φ+(1/2) k^2 S_2(φ)+(1/2)(3/4) k^4 S_4(φ)+ …
第2種楕円積分 E2(φ,k)=∫φ_0 √(1-k^2 sin^2φ) dφ=φ+(-1/2) k^2 S_2(φ)+(-1/2)(3/4) k^4 S_4(φ)+(-1/2)(3/4)(5/6) k^6 S_6(φ)+ …
(ここで∫φ_0 sin^nφ dφ=S_n(φ)と略記)
S_n(φ)は漸化式 S_n(φ) =( (n-1)S_n-2(φ) -sin^(n-1)φcosφ)/nで計算出来る
S_0(φ) =φ
S_2(φ) =(1/2)S_0(φ)-sinφcosφ/2=φ/2-sinφcosφ/2
S_4(φ) =(3/4)S_2(φ)-sin^3φcosφ/4=(3/4)(φ/2-sinφcosφ/2)-sin^3φcosφ/4=(3/8)φ-(1/8)sinφcosφ(2sin^2φ+3)
…
○上記 ホームページの解説により
無限級数解も実用的であることがわかりました。
処理速度、計算精度ともに ガウスールジャンドルの数値積分 程度 と思われます。
計算誤差の問題はないようです。倍精度なら倍精度で返せるようです。
=end
include Math
def E1(φ, k=nil) #第１種楕円積分
if k == nil
k = φ
φ = PI / 2
end
if k < 0 or 1 < k then return print "h#E1 φ, %s 計算できませんでした\n" % k end
if k > 1-1e-15 then return log( (1.0 + sin(φ)) / (1.0 - sin(φ))) / 2.0 end
n, f, a, s_n = 2.0, 0.0, 1.0, φ
c, s = cos(φ), sin(φ)
j = 0
while (a * s_n).abs > 1e-20
f += a * s_n
a *= (n - 1) / n * k * k
s_n = ( (n - 1) * s_n - s ** (n - 1) * c) / n
n += 2
j += 1
if j > 30000 then break end
end
return f
end
def E2(φ, k=nil) #第２種楕円積分
if k == nil
k = φ
φ = PI / 2
end
if k < 0 or 1 < k then return print "h#E2 φ, %s 計算できませんでした\n" % k end
if k > 1-1e-15 then return sin(φ) end
n, e, a, s_n = 2.0, 0.0, 1.0, φ
c, s = cos(φ), sin(φ)
j = 0
while (a * s_n).abs > 1e-20
e += a * s_n
a *= (n - 3) / n * k * k
s_n = ( (n - 1) * s_n - s ** (n - 1) * c) / n
n += 2
j += 1
if j > 30000 then break end
end
return e
end
φ= PI/6
k = sqrt(1-0.5**2)
p E1(k) #第１種完全楕円積分
p E1(φ, k) #第１種不完全楕円積分
p E2(k) #第２種完全楕円積分
p E2(φ, k) #第２種不完全楕円積分
def arc_l(r,p1,p2,w)
ek = E2(PI/2,k=sqrt(1-w*w))
if p1 == p2 then p2 += 360 end
lp = ( (p2-p1).abs/360.0).to_i*4.0*ek
p1 += 360 while p1 < 0
p1 -= 360 while p1 >= 360
p2 += 360 while p2 < 0
p2 -= 360 while p2 >= 360
if p1 == p2; return r*lp; end
n1 = 1+(p1/90.0).to_i
if n1 > 4; n1 = 4; end
lp += ek*( (n1+1.0).remainder(2))+E2(apa(n1,p1)*PI/180,k)*(-1.0)**(n1-1)
n2 = 1+(p2/90.0).to_i
if n2 > 4; n2 = 4; end
if n1 == n2; n1 += 4; end
n3 = 5+2*( (n2-1.0)/2.0).to_i
lp += ek*( (n3-n1).remainder(4))+E2(apa(n2,p2)*PI/180,k)*(-1.0)**n2
if lp < 0.0; lp += 4.0 * ek; end
return r*lp
end
def apa(o,q)
return (90.0*(2*( (o+1)/2).to_i-1)-q)*(-1)**(o-1)
end
p arc_l(半径=50, 始角=30, 終角=120, 扁平比=0.5)
