jw_cad 外部変形 － (742) gawkで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する －

外部変形は データのやり取りをテキストファイルで行うので プログラム言語は 自由に選ぶことができます。図形は機能的かつシンプルなため、数多くのユーザーに受け入れられています。

gawkで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する

その① 台形則で楕円積分

:gawkで楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する
@echo off
for /f %%n in ('gawk "/^#!/ { print NR }" %0') do (
  copy jwc_temp.txt myfiles > nul
  more +%%n %0 | gawk -f - myfiles > jwc_temp.txt
)
goto:eof

REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e

#!ここから more +n %0 の n 行目:最初の行は 0 行
function rad(x) { return x * PI / 180 }
function apa(o, p) { return (90 * (int( (o + 1) / 2) * 2 - 1) - p) * (-1) ^ (o - 1) }
function f(x, k) { return sqrt(1 - (k * sin(x)) ^ 2) }
function arc_l(r, p1, p2, w,  k, Ek, lp, n1, n2, n3) {
  #r  半径
  #p1 始角 ﾟ
  #p2 終角 ﾟ
  #w  扁平率
  #k  離心率
  #Ek 第二種完全楕円積分の値
  if (w > 1) {
    w = 1 / w
    r /= w
    p1 -= 90
    p2 -= 90
  }
  k = sqrt(1 - w * w)
  Ek = E(rad(90), k)
  while (p1 <  0  ) p1 += 360
  while (p1 >= 360) p1 -= 360
  while (p2 <  0  ) p2 += 360
  while (p2 >= 360) p2 -= 360
  if (p1 == p2) return r * 4 * Ek
  n1 = int(p1 / 90) + 1
  if (n1 > 4) n1 = 4
  lp = Ek * ( (n1 + 1) % 2) + E(rad(apa(n1, p1)), k) * (-1) ^ (n1 - 1)
  n2 = int(p2 / 90) + 1
  if (n2 > 4) n2 = 4
  if (n1 == n2) n1 += 4
  n3 = 2 * int( (n2 - 1) / 2) + 5
  lp += Ek * ( (n3 - n1) % 4) + E(rad(apa(n2, p2)), k) * (-1) ^ n2
  if (lp < 0) lp += 4 * Ek
  return r * lp
}
function E(p, k,  i, l) { #台形則
  N = 10000
  hh = 1.0 / N
  x = 0.0; l = 0.0
  for (i = 1; i <= N; i++) {
    s1 = f(x * p, k)
    x += hh
    s2 = f(x * p, k)
    l += (s1 + s2) * hh / 2.0
  }
  return l * p
}
BEGIN { CONVFMT = OFMT = "%.15g"; PI = atan2(1, 1) * 4 }
/^hs/ { for (i = 0; i <= 15; i++) hs[sprintf("%x", i)] = $(i + 2) }
/^lg[0-9a-f]/ { plg = substr($0, 3, 1); if (lg == "") lg = plg }
/^ci/ {
  p1 = p2 = 0; r = $4; w = 1
  if (NF == 8) {
    p1 = $5; p2 = $6; w = $7
    if (w > 1) { w = 1 / w; r /= w; p1 -= 90; p2 -= 90 }
  }
  if (w == 1) {
    lr = r * (p1 == p2 ? 360 : p2 - p1) * PI / 180
  } else {
    lr = arc_l(r, p1, p2, w)
  }
  printf "h#Lr = %.03f ", lr * hs[lg] / hs[plg]
}

その② シンプソン則で楕円積分

:楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する Simpson
@echo off
for /f %%n in ('gawk "/^#!/ { print NR }" %0') do (
  copy jwc_temp.txt myfiles > nul
  more +%%n %0 | gawk -f - myfiles > jwc_temp.txt
)
goto:eof

REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e

#!
function rad(x) { return x * PI / 180 }
function apa(o, p) { return (90 * (int( (o + 1) / 2) * 2 - 1) - p) * (-1) ^ (o - 1) }
function f(x, k) { return sqrt(1 - (k * sin(x)) ^ 2) }
function arc_l(r, p1, p2, w,  k, Ek, lp, n1, n2, n3) {
  #r  半径
  #p1 始角 ﾟ
  #p2 終角 ﾟ
  #w  扁平率
  #k  離心率
  #Ek 第二種完全楕円積分の値
  if (w > 1) {
    w = 1 / w
    r /= w
    p1 -= 90
    p2 -= 90
  }
  k = sqrt(1 - w * w)
  Ek = E(rad(90), k)
  while (p1 <  0  ) p1 += 360
  while (p1 >= 360) p1 -= 360
  while (p2 <  0  ) p2 += 360
  while (p2 >= 360) p2 -= 360
  if (p1 == p2) return r * 4 * Ek
  n1 = int(p1 / 90) + 1
  if (n1 > 4) n1 = 4
  lp = Ek * ( (n1 + 1) % 2) + E(rad(apa(n1, p1)), k) * (-1) ^ (n1 - 1)
  n2 = int(p2 / 90) + 1
  if (n2 > 4) n2 = 4
  if (n1 == n2) n1 += 4
  n3 = 2 * int( (n2 - 1) / 2) + 5
  lp += Ek * ( (n3 - n1) % 4) + E(rad(apa(n2, p2)), k) * (-1) ^ n2
  if (lp < 0) lp += 4 * Ek
  return r * lp
}
function E(p, k,  i, l) { #Simpson
  N = 1000
  x = H = p / (N * 2)
  i = s1 = 0
  s2 = f(x, k)
  while (i != N-1) {
    i++; 
    x += H; 
    s1 += f(x,k); 
    x += H; 
    s2 += f(x,k) 
  }
  return H / 3 * (f(0, k) + f(p, k) + 4 * s2 + 2 * s1)
}
BEGIN { CONVFMT = OFMT = "%.15g"; PI = atan2(1, 1) * 4 }
/^hs/ { for (i = 0; i <= 15; i++) hs[sprintf("%x", i)] = $(i + 2) }
/^lg[0-9a-f]/ { plg = substr($0, 3, 1); if (lg == "") lg = plg }
/^ci/ {
  CONVFMT = OFMT = "%.15g"
  r = $4; p1 = p2 = 0; w = 1
  if (NF == 8) {
    p1 = $5; p2 = $6; w = $7
    if (w > 1) { w = 1 / w; r /= w; p1 -= 90; p2 -= 90 }
  }
  if (w == 1) {
    lr = r * (p1 == p2 ? 360 : p2 - p1) * PI / 180
  } else {
    lr = arc_l(r, p1, p2, w)
  }
  printf "h#Lr = %.03f ", lr * hs[lg] / hs[plg]
}

その③ ガウスールジャンドル則で楕円積分

:楕円弧の長さ(表示画面のｽｹｰﾙ)を計測表示する Gauss
@echo off
for /f %%n in ('gawk "/^#!/ { print NR }" %0') do (
  copy jwc_temp.txt myfiles > nul
  more +%%n %0 | gawk -f - myfiles > jwc_temp.txt
)
goto:eof

REM #jww
REM #1ci 長さを計測する円・円弧を指示してください
REM #e

#!
BEGIN { CONVFMT = OFMT = "%.15g"; PI = atan2(1, 1) * 4 }
function abs(x) { return x > 0 ? x : -x }
function rad(x) { return x * PI / 180 }
function apa(o, p) { return (90 * (int( (o + 1) / 2) * 2 - 1) - p) * (-1) ^ (o - 1) }
function fx(x, k) { return sqrt(1 - (k * sin(x)) ^ 2) }
function arc_l(r, p1, p2, w,  k, Ek, lp, n1, n2, n3) {
  #r  半径
  #p1 始角 ﾟ
  #p2 終角 ﾟ
  #w  扁平率
  #k  離心率
  #Ek 第二種完全楕円積分の値
  if (w > 1) {
    w = 1 / w
    r /= w
    p1 -= 90
    p2 -= 90
  }
  k = sqrt(1 - w * w)
  Ek = E(rad(90), k)
  while (p1 <  0  ) p1 += 360
  while (p1 >= 360) p1 -= 360
  while (p2 <  0  ) p2 += 360
  while (p2 >= 360) p2 -= 360
  if (p1 == p2) return r * 4 * Ek
  n1 = int(p1 / 90) + 1
  if (n1 > 4) n1 = 4
  lp = Ek * ( (n1 + 1) % 2) + E(rad(apa(n1, p1)), k) * (-1) ^ (n1 - 1)
  n2 = int(p2 / 90) + 1
  if (n2 > 4) n2 = 4
  if (n1 == n2) n1 += 4
  n3 = 2 * int( (n2 - 1) / 2) + 5
  lp += Ek * ( (n3 - n1) % 4) + E(rad(apa(n2, p2)), k) * (-1) ^ n2
  if (lp < 0) lp += 4 * Ek
  return r * lp
}
function E(p, k,  i, l, n, e, a, s_n) {
  if (k == "") {
    k = p
    p = PI / 2
  }
  if (PREC > 99) {
    if (k < 0 || 1 < k) { printf("h#E p, %s 計算できませんでした\n", k); return }
    if (k > 1-1e-30) { return sin(p) }
    n = 2; e = 0; a = 1; s_n = p
    c = cos(p); s = sin(p)
    j = 0
    while (abs(a * s_n) > 1e-60) {
      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) { break }
    }
    return e
  } else {
    l = 0
    if (N > 1) for (i = 1; i <= int(N / 2); i++) l += __h__[i] * fx(__e__[i] * p, k)
    if (N % 2 == 1) l += __h__[M] * fx(__e__[M] * p, k) / 2 #奇数個の場合
    return l * p
  }
}
function gear(bpa,  ww) {
  if (bpa < 0 || N > 0) { #140429
    if (N == 0) N = int(abs(bpa))
  } else {
    ww = bpa; if (bpa > 1) ww = 1 / bpa
    if (ww == 0 || ww == "") ww = 0.1
    N = int(12 * (1.5 / ww))
    if (N > 300) N = 300
  }
  eh(N)
  return bpa
}
function eh(N) { # 積分点デ－タ
  M  = int(N / 2) + N % 2
  E1 = N * (N + 1)
  for (I = 1; I <= M; I++) { #DO 1
    T   = (4 * I - 1) * PI / (4 * N + 2)
    X0  = (1. - (1. - 1. / N) / (8. * N * N)) * cos(T)
    PKM1= 1.
    PK  = X0
    for (K = 2; K <= N; K++) { #DO 3
      T1  = X0 * PK
      PKP1= T1 - PKM1 - (T1 - PKM1) / K + T1
      PKM1= PK
      PK  = PKP1
    } # 3
    DEN = 1. - X0 * X0
    D1  = N * (PKM1 - X0 * PK)
    DPN = D1 / DEN
    D2PN= (2. * X0 * DPN - E1 * PK) / DEN
    D3PN= (4. * X0 * D2PN + (2. - E1) * DPN) / DEN
    D4PN= (6. * X0 * D3PN + (6. - E1) * D2PN) / DEN
    U   = PK / DPN
    V   = D2PN / DPN
    H   = -U * (1. + .5 * U * (V + U * (V * V - D3PN / (3. * DPN))))
    P   = PK + H * (DPN + .5 * H * (D2PN + H / 3. * (D3PN + .25 * H * D4PN)))
    DP  = DPN + H * (D2PN + .5 * H * (D3PN + H * D4PN / 3.))
    H   = H - P / DP
    XI  = -(X0 + H)
    FX  = D1 - H * E1 * (PK + .5 * H * (DPN + H / 3. * (D2PN + .25 * H * (D3PN + .2 * H * D4PN))))
    WI  = 2. * (1. - XI * XI) / (FX * FX)
    __e__[I] = XI
    __h__[I] = WI
  } # 1
}
BEGIN { CONVFMT = OFMT = "%.15g"; PI = atan2(1, 1) * 4 }
/^hs/ { for (i = 0; i <= 15; i++) hs[sprintf("%x", i)] = $(i + 2) }
/^lg[0-9a-f]/ { plg = substr($0, 3, 1); if (lg == "") lg = plg }
/^ci/ {
  p1 = p2 = 0; r = $4; w = 1
  if (NF == 8) {
    p1 = $5; p2 = $6; w = $7
    if (w > 1) { w = 1 / w; r /= w; p1 -= 90; p2 -= 90 }
  }
  if (w == 1) {
    lr = r * (p1 == p2 ? 360 : p2 - p1) * PI / 180
  } else {
    lr = arc_l(r, p1, p2, gear(w))
  }
  printf "h#Lr = %.03f ", lr * hs[lg] / hs[plg]
}