jw_cad 外部変形 - (862) pythonで楕円弧の長さ(表示画面のスケール)を計測表示する -

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

 

pythonで楕円弧の長さ(表示画面のスケール)を計測表示する

:pythonで楕円弧の長さ(表示画面のスケール)を計測表示する Gauss-Legendre
@echo off
if not exist %~dpn0.py (
  for /f "delims=:" %%n in ('findstr /n "^#!" %0') do (
    more +%%n %0 > %~dpn0.py
  )
)
python %~dpn0.py > jwc_temp.txt
goto:eof

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

#!この次の行からプログラムを書いてください
# coding: shift_jis
# 円データのスケール補正
def sci_bai(s, p=1):
  if p != 1:
    i = 0
    for x in s[:3]:
      s[i] = x * p
      i += 1
  return s

def sgn(x):
  if x == None or x == 0:
    return 0
  elif x > 1:
    return 1
  else:
    return -1

def fx(x, k):
  return sqrt(1.0 - (k * sin(x)) ** 2)

def apa(o, p):
  return (90.0 * (2 * int( (o + 1) / 2) - 1) - p) * (-1) ** (o - 1)

def ee(p, g):
  if type(g) == int or type(g) == float:
    k = g
    n = 60
    e, h = gauss(n)
  else:
    k, n, e, h = g
  if k == 0 or n < 1:
    return p
  l = 0.0
  m = int(n / 2.0) + n % 2
  for i in range(1, int(n / 2.0) + 1):
    l += h[i] * fx(e[i] * p, k)
  if m == int(n / 2.0) + 1: #奇数個の場合
    l += h[m] * fx(e[m] * p, k) / 2.0
  return l * p

def gauss(n):
  m = int(n / 2.0) + n % 2
  ee = [0] * (m + 1)
  hh = [0] * (m + 1)
  e1 = n * (n + 1.0)
  for i in range(1, m + 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 range(2, n + 1):
      t1   = x0 * pk
      pkp1 = t1 - pkm1 - (t1 - pkm1) / k + t1
      pkm1 = pk
      pk   = pkp1
      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)
      ee[i] = xi
      hh[i] = wi
  return ee, hh

def arc_l(rc=1.0, p1=0.0, p2=0.0, w=1.0):
  if w == 1:
    p = p2 - p1
    while p <= 0: p += 360
    while p > 360: p -= 360
    return p * pi / 180.0 * rc
  else:
    if w > 1:
      w = 1.0 / w
      rc = rc / w
      p1 = p1 - 90.0
      p2 = p2 - 90.0
    n = int(12 * (1.5 / w)) #ガウス積分点の数
    if n > 300 : n = 300
    k = sqrt(1.0 - w ** 2)
    e, h = gauss(n)
    g = [k, n, e, h]
    ek = ee(pi / 2.0, g)
    # 弧長をガウス積分(第2種楕円積分)を利用して計算する
    while p1 < 0: p1 += 360
    while p1 >= 360: p1 -= 360
    while p2 < 0: p2 += 360
    while p2 >= 360: p2 -= 360
    if p1 == p2:
      lp = 4.0 * ek
    else:
      n1 = 1 + int(p1 / 90.0)
      if n1 > 4: n1 = 4
      lp = ek * ((n1 + 1.0) % 2) + ee(apa(n1, p1) * pi / 180, g) * (-1.0) ** (n1 - 1)
      n2 = 1 + int(p2 / 90.0)
      if n2 > 4: n2 = 4
      if n1 == n2: n1 += 4
      n3 = 5 + 2 * int((n2 - 1.0) / 2.0)
      lp += ek * ((n3 - n1) % 4) + ee(apa(n2, p2) * pi / 180, g) * (-1.0) ** n2
      if lp < 0.0: lp += 4.0 * ek
    lr = lp * rc
  return lr

import sys,re
from math import *

hs = {}
lg = ""

try :
  f = open("myfiles", "r")
  for line in f:
    F = line.split()
    if re.compile("^hs").search(line):
      i = 0
      for x in F[1:]:
        hs[str(hex(i)[2:])] = float(x)
        i += 1
    if re.compile("^lg").search(line):
      plg = line[2:3]
      if lg == "": lg = plg
    if re.compile("^ci").search(line):
      ci1 = F[1:]
      if len(ci1) == 3: ci1 = ci1 + [0, 360, 1, 0]
      ci1 = map(float, ci1)
      ci1 = sci_bai(ci1, hs[lg] / hs[plg])
      x1, y1, r1, p1, p2, w, d = ci1
      lr = arc_l(r1, p1, p2, w)
      print("h#円弧長 L = %.3f" % lr)
  f.close()
except :
  sys.exit()