jw_cad 外部変形 - (1086) clispで多角形の面積を計算する -

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

 

clispで多角形の面積を計算する

:clispで多角形の面積を計算する
@echo off
for /f "delims=:" %%n in ('findstr /n "^#!" %0') do (
  more +%%n %0 | clisp -q > nul
)
goto:eof

REM #jww
REM #1-%d 点を指示してください
REM #99#%d
REM #e

#!この次の行からプログラムを書いてください
;ユーザ定義定数
;文字列 str の old に一致する部分を new で置き換える
(defun gsub (old new str)
(let (s e)
  (if (search old str)
      (progn
        (setq s (search old str))
        (setq e (+ s (length old)))
        (setq str (concatenate 'string (subseq str 0 s) new (subseq str e)))
        (gsub old new str)
      )
      str
  )
))

( ;文字列 "x" を 数値 xd0 (倍精度 double-float) に変換
  defun to_f (x)
    (if (equal x nil) (setq x "0"))
    (if (listp x)
        (mapcar #'to_f x)
        (progn
          (if (typep x 'double-float)
              x
              (progn
                (if (stringp x) nil (setq x (write-to-string x)))
                (if (search "e" x)
                    (progn
                      (setq s (search "e" x))
                      (read-from-string
                        (concatenate 'string (subseq x 0 s) "d" (subseq x (+ 1 s) (length x)))
                      )
                    ) ;progn (search "e" x)
                    (if (search "d" x)
                        x
                        (read-from-string
                          (concatenate 'string x "d0")
                        )
                    )
                )
              ) ;progn (typep x 'double-float)
          )
        ) ;progn (listp x)
    )
)

(defun area (hp) ;多角形の面積
(let (n ax sx sy jx jy
      i x y gx gy ix iy xmin xmax ymin ymax)
  (setq hp (to_f hp))
  (setq n (- (length hp) 1))
  (setq ax 0
        sx 0
        sy 0
        jx 0
        jy 0)
  (if (equal (elt hp 0) (elt hp n))
    (setq n (- n 1))
    (progn
      (setq hp (append hp '( (0 0))))
      (setf (elt hp (+ n 1)) (elt hp 0))
    )
  )
  (setq x (loop for i to n collect 0))
  (setq y (loop for i to n collect 0))
  (loop for i to n
    do
    (setq x (elt (elt hp i) 0))
    (setq y (elt (elt hp i) 1))
    (setq x1 (elt (elt hp i) 0))
    (setq y1 (elt (elt hp i) 1))
    (setq x2 (elt (elt hp (+ i 1)) 0))
    (setq y2 (elt (elt hp (+ i 1)) 1))
;    a += 0.5 * (x1 - x2) * (y1 + y2) #左回りが正
;    sx += ( (x1 - x2) * (y1 * (2 * x1 + x2) + y2 * (2 * x2 + x1)) / 6.0)
;    sy += ( (y2 - y1) * (x1 * (2 * y1 + y2) + x2 * (2 * y2 + y1)) / 6.0)
;    ix += ( (x1 - x2) * (y1 * (2 * x1 + x2) ** 2 + y2 * (2 * x2 + x1) ** 2 + (x1 - x2) ** 2 * (y1 + y2) / 2) / 18.0)
;    iy += ( (y2 - y1) * (x1 * (2 * y1 + y2) ** 2 + x2 * (2 * y2 + y1) ** 2 + (y1 - y2) ** 2 * (x1 + x2) / 2) / 18.0)
    (setq ax (+ ax (* 0.5d0 (- x1 x2) (+ y1 y2))))
    (setq sx (+ sx (/ (* (- x1 x2) (+ (* y1 (+ (* 2 x1) x2)) (* y2 (+ (* 2 x2) x1)))) 6)))
    (setq sy (+ sy (/ (* (- y2 y1) (+ (* x1 (+ (* 2 y1) y2)) (* x2 (+ (* 2 y2) y1)))) 6)))
    (setq jx (+ jx (/ (* (- x1 x2) (+ (* y1 (pow (+ (* 2 x1) x2) 2)) (* y2 (pow (+ (* 2 x2) x1) 2)) (/ (* (pow (- x1 x2) 2) (+ y1 y2)) 2))) 18)))
    (setq jy (+ jy (/ (* (- y2 y1) (+ (* x1 (pow (+ (* 2 y1) y2) 2)) (* x2 (pow (+ (* 2 y2) y1) 2)) (/ (* (pow (- y1 y2) 2) (+ x1 x2)) 2))) 18)))
  )
  (setq gx (/ sx ax)) ;図心
  (setq gy (/ sy ax))
  (setq ix (- jx (* ax gx gx))) ;図心まわりの2次モーメント(Y軸) Iy-y
  (setq iy (- jy (* ax gy gy))) ;図心まわりの2次モーメント(X軸) Ix-x
  (setq xmin (- gx (min x))) ;図心からX軸方向の両側の最縁端距離
  (setq xmax (- gx (max x)))
  (setq ymin (- gy (min y))) ;図心からY軸方向の両側の最縁端距離
  (setq ymax (- gy (max y)))
  (list ax (list gx gy) (list ix iy) (list xmin xmax ymin ymax))
)); public :area

;本文 / JW_CAD 外部変形
( ;jwc_temp.txt から 入力
  with-open-file (f "jwc_temp.txt" :direction :input)
  (loop for line = (read-line f nil) while line do
    (setq a (regexp:regexp-split "\\s\\+" line))
    ;指示点データ
    (if (regexp:match "^hp[1-9][0-9]\\?\\(ln\\|ci\\|ch\\)\\?-\\?" line)
        (progn
          (setq hpn (gsub "hp" "" (car a)))
          (setq hpn (gsub "ln" "" hpn))
          (setq hpn (gsub "ci" "" hpn))
          (setq hpn (gsub "ch" "" hpn))
          (setq hpn (gsub "#" "" hpn))
          (setq hpn (parse-integer (gsub "-" "" hpn)))
          (if (>= (- hpn (length hp)) 0)
              (loop for i to (- hpn (length hp)) do
                (setq hp (append hp '( (0 0))))
              )
          )
          (setf (elt hp hpn) (cdr a))
        )
    )
  )
)

( ;多角形の面積 a : (ax) (gx gy) (ix iy) (xmin xmax ymin ymax)
  setq a (area (hp#n 1 hpn))
)

( ;jwc_temp.txt へ 出力
  with-open-file (f "jwc_temp.txt" :direction :output)
  (format f "h#Area = ~,1f Iy = ~,1f Ix = ~,1f~%"
    (elt a 0) ; 面積
    (elt (elt a 2) 0) ; 図心まわりの2次モーメント(Y軸) Iy-y
    (elt (elt a 2) 1) ; 図心まわりの2次モーメント(X軸) Ix-x
  )
)

 

 

○計算式について

線データで構成された閉鎖形の面積はつぎのように計算します。

台形を2つの三角形に分けて考えるとイメージしやすい。(△1-2-p1, ▽p1-p2-2)

面積

   ax = ∫ x0 y0 dA = ∑ ∫|x1,x2| f(x) dx
   ay = ∫ x0 y0 dA = ∑ ∫|y1,y2| g(y) dy
   ただし、正負の符号を考えて 線データは以下のようにする。
   -f(x) = y1 + (y2 - y1) / (x2 - x1) * (x - x1)
    g(y) = x1 + (x2 - x1) / (y2 - y1) * (y - y1)
  左図 で 線とX軸から成る台形の面積 ax を考えると
    ax = (x1 - x2) * (y1 + y2) / 2
       = (x1 * y1 + x1 * y2 - x2 * y1 - x2 * y2) / 2
      線とY軸から成る台形の面積 ay は
    ay =-(y1 - y2) * (x1 + x2) / 2
       =-(y1 * x1 + y1 * x2 - y2 * x1 - y2 * x2) / 2
  右図 で 線とX軸から成る台形を考えたときの三角形の面積 ax は
    ax = (x1 * y1 + x1 * y2 - x2 * y1 - x2 * y2) / 2
       + (x2 * y2 + x2 * y3 - x3 * y2 - x3 * y3) / 2
       + (x3 * y3 + x3 * y1 - x1 * y3 - x1 * y1) / 2
       = ∑(x[i] * y[i+1] - x[i+1] * y[i]) / 2
      線とY軸から成る台形を考えたときの三角形の面積 ay は
    ay = ax = ∑(x[i] * y[i+1] - x[i+1] * y[i]) / 2
  となります。辺の数だけ繰り返して求めます。
  → 左回りを正としています。

 

断面1次モーメント

   sx = ∫ x1 y0 dA = ∑ ∫|x1,x2| x f(x) dx
   sy = ∫ x0 y1 dA = ∑ ∫|y1,y2| y g(y) dy
   ただし、正負の符号を考えて 線データは以下のようにする。
   -f(x) = y1 + (y2 - y1) / (x2 - x1) * (x - x1)
    g(y) = x1 + (x2 - x1) / (y2 - y1) * (y - y1)
 左図 で 線とX軸から成る台形の1次モーメント sx
         線とY軸から成る台形の1次モーメント sy は
    sx = (x1 - x2) * (y1 * (2 * x1 + x2) + y2 * (2 * x2 + x1)) / 6
       = (x1 * y2 - x2 * y1) * (x1 + x2) / 6 + (x1 ** 2 * y1 - x2 ** 2 * y2) / 3
    sy =-(y1 - y2) * (x1 * (2 * y1 + y2) + x2 * (2 * y2 + y1)) / 6
       =-(y1 * x2 - y2 * x1) * (y1 + y2) / 6 + (y1 ** 2 * x1 - y2 ** 2 * x2) / 3
 となります。
  右図 で 線とX軸から成る台形を考えたときの三角形の1次モーメント sx は
    sx = (x1 * y2 - x2 * y1) * (x1 + x2) / 6 + (x1 ** 2 * y1 - x2 ** 2 * y2) / 3
       + (x2 * y3 - x3 * y2) * (x2 + x3) / 6 + (x2 ** 2 * y2 - x3 ** 2 * y3) / 3
       + (x3 * y1 - x1 * y3) * (x3 + x1) / 6 + (x3 ** 2 * y3 - x1 ** 2 * y1) / 3
       = ∑(x[i] * y[i+1] - x[i+1] * y[i]) * (x[i] + x[i+1]) / 6
      線とY軸から成る台形を考えたときの三角形の1次モーメント sy は
    sy = ∑(x[i] * y[i+1] - x[i+1] * y[i]) * (y[i] + y[i+1]) / 6
 と整理できます。
 図心 gx, gy は 以下のように求めます。
   gx = sx / ax
   gy = sy / ay
  → 左回りが正となるように補正しています。

 

断面2次モーメント

   ix = ∫ x2 y0 dA = ∑ ∫|x1,x2| x * x * f(x) dx
   iy = ∫ x0 y2 dA = ∑ ∫|y1,y2| y * y * g(y) dy
   ただし、正負の符号を考えて 線データは以下のようにする。
   -f(x) = y1 + (y2 - y1) / (x2 - x1) * (x - x1)
    g(y) = x1 + (x2 - x1) / (y2 - y1) * (y - y1)
 左図 で 線とX軸から成る台形の2次モーメント ix
         線とY軸から成る台形の2次モーメント iy は
    ix = (x1 - x2) * (y1 * (2 * x1 + x2) ** 2 + y2 * (2 * x2 + x1) ** 2 + (x1 - x2) ** 2 * (y1 + y2) / 2) / 18
       = (x1 * y2 - x2 * y1) * (x1 ** 2 + x1 * x2 + x2 ** 2) / 12 + (x1 ** 3 * y1 - x2 ** 3 * y2) / 4
    iy =-(y1 - y2) * (x1 * (2 * y1 + y2) ** 2 + x2 * (2 * y2 + y1) ** 2 + (y1 - y2) ** 2 * (x1 + x2) / 2) / 18
       =-(y1 * x2 - y2 * x1) * (y1 ** 2 + y1 * y2 + y2 ** 2) / 12 + (y1 ** 3 * x1 - y2 ** 3 * x2) / 4
 となります。
  右図 で 線とX軸から成る台形を考えたときの三角形の2次モーメント ix は
    ix = (x1 * y2 - x2 * y1) * (x1 ** 2 + x1 * x2 + x2 ** 2) / 12 + (x1 ** 3 * y1 - x2 ** 3 * y2) / 4
       + (x2 * y3 - x3 * y2) * (x2 ** 2 + x2 * x3 + x3 ** 2) / 12 + (x2 ** 3 * y2 - x3 ** 3 * y3) / 4
       + (x3 * y1 - x1 * y3) * (x3 ** 2 + x3 * x1 + x1 ** 2) / 12 + (x3 ** 3 * y3 - x1 ** 3 * y1) / 4
       = ∑(x[i] * y[i+1] - x[i+1] * y[i]) * ((x[i] + x[i+1]) ** 2 - x[i] * x[i+1]) / 12
      線とY軸から成る台形を考えたときの三角形の2次モーメント iy は
    iy = ∑(x[i] * y[i+1] - x[i+1] * y[i]) * ((y[i] + y[i+1]) ** 2 - y[i] * y[i+1]) / 12
 と整理できます。
 図心回りの2次モーメント jxc, jyc は 以下のように求めます。
   jyc = ix - ax * gx ** 2
   jxc = iy - ay * gy ** 2
  → 左回りが正となるように補正しています。