[原创]：马尔法蒂问题的lisp求解

highflybird

所谓的马尔法蒂问题是：
在一个已知三角形内画三个圆，每个圆与其它两个圆以及三角形的两边相切。
这个问题有几何解，即用纯几何法来完成。
我这里编写了一个lisp程序求解。
加载程序，在CAD下运行mal 即可。

2. 
   
3. ;;*****************************************************************************
   
4. ;;马尔法蒂问题lisp求解---------------------------------------------------------
   
5. (defun C:Mal (/           pa        pb   pc          a    b    c         p    xa   xb        xc
   
6.               ya   yb        yc   aga  agb  agc  ang         ta   tb   tc        ja
   
7.               jb   jc        ha   hb          hc   vpa  vpb         vpc  cen  ra        rb
   
8.               rc   la        lb   lc          cena cenb cenc
   
9.              )
   
10. ;;(defun C:Mal ()
    
11.   (graphscr)
    
12.   (setq oldmode (getvar "osmode"))
    
13.   (setq oce (getvar "cmdecho"))
    
14.   (setvar "cmdecho" 0)
    
15.   ;;输入数据
    
16.   (command ".ucs" "W")
    
17.   (setq pa (getpoint "请输入第一点：\n"))
    
18.   (setq pb (getpoint "请输入第二点：\n"))
    
19.   (setq pc (getpoint "请输入第三点：\n"))
    
20.   (command ".ucs" "O" pa)
    
21.   (setq        pa (trans pa 0 1)
    
22.         pb (trans pb 0 1)
    
23.         pc (trans pc 0 1)
    
24.   )
    
25.   ;;求三角形边长和半周长  
    
26.   (setq        a (distance pb pc)
    
27.         b (distance pc pa)
    
28.         c (distance pa pb)
    
29.         p (/ (+ a b c) 2)
    
30.   )
    
31.   ;;计算切线长
    
32.   (setq        xa (sqrt (abs (- (* p p) (* p a))))
    
33.         xb (sqrt (abs (- (* p p) (* p b))))
    
34.         xc (sqrt (abs (- (* p p) (* p c))))
    
35.   )
    
36.   (setq        ya (sqrt (* p a))
    
37.         yb (sqrt (* p b))
    
38.         yc (sqrt (* p c))
    
39.   )
    
40.   (setq        aga (angle '(0 0) (list xa ya))
    
41.         agb (angle '(0 0) (list xb yb))
    
42.         agc (angle '(0 0) (list xc yc))
    
43.   )
    
44.   (setq ang (/ (+ aga agb agc) 2))
    
45.   (setq ta (* p (* (sin (- ang aga))) (* (sin (- ang aga)))))
    
46.   (setq tb (* p (* (sin (- ang agb))) (* (sin (- ang agb)))))
    
47.   (setq tc (* p (* (sin (- ang agc))) (* (sin (- ang agc)))))
    
48.   ;;***************************************************************************
    
49.   ;;求三角形内心---------------------------------------------------------------
    
50.   (defun cen_incir (pa pb pc)
    
51.     (setq jc (angle pa pb)
    
52.           ja (angle pb pc)
    
53.           jb (angle pc pa)
    
54.     )
    
55.     (setq ha (/ (+ jb jc pi) 2)
    
56.           hb (/ (+ jc ja pi) 2)
    
57.           hc (/ (+ ja jb pi) 2)
    
58.     )
    
59.     (setq vpa (polar pa ha p)
    
60.           vpb (polar pb hb p)
    
61.           vpc (polar pc hc p)
    
62.     )
    
63.     (inters pa vpa pb vpb nil)
    
64.   )
    
65.   (setq cen (cen_incir pa pb pc))
    
66.   ;;***************************************************************************
    
67.   ;;求每个圆的半径，圆心位置---------------------------------------------------
    
68.   (defun tan (x)
    
69.     (/ (sin x) (cos x))
    
70.   )
    
71.   ;;定义正切函数
    
72.   (if (> 1e-16
    
73.          (abs (* (sin (- jb jc)) (sin (- jc ja)) (sin (- ja jb))))
    
74.       )
    
75.     (progn
    
76.       (alert "你输入的三点在一条直线上,请重新输入!")
    
77.       (command ".UCS" "P")
    
78.       (command ".UCS" "P")
    
79.       (setvar "osmode" oldmode)
    
80.       (setvar "cmdecho" oce)
    
81.       (princ)
    
82.     )
    
83.     ;;判断输入的三点是否在同一条直线上
    
84.     (progn
    
85.       (setq ra (* ta (/ 1 (abs (tan (/ (- jb jc) 2))))))
    
86.       (setq rb (* tb (/ 1 (abs (tan (/ (- jc ja) 2))))))
    
87.       (setq rc (* tc (/ 1 (abs (tan (/ (- ja jb) 2))))))
    
88.       (if (or nil (< (abs ra) 1e-8) (< (abs ra) 1e-8) (< (abs ra) 1e-8))
    
89.         (progn
    
90.           (princ "你输入的三点在一条直线上,请重新输入!")
    
91.           (command ".UCS" "P")
    
92.           (command ".UCS" "P")
    
93.           (setvar "cmdecho" oce)
    
94.           (setvar "osmode" oldmode)
    
95.           (princ)
    
96.         )
    
97.         ;;判断圆的半径是否为零
    
98.         (progn
    
99.           (defun gougu (x y)
    
100.             (sqrt (+ (* x x) (* y y)))
     
101.           )
     
102.           (setq        la (gougu ta ra)
     
103.                 lb (gougu tb rb)
     
104.                 lc (gougu tc rc)
     
105.           )
     
106.           (setq        cena (polar pa (angle pa cen) la)
     
107.                 cenb (polar pb (angle pb cen) lb)
     
108.                 cenc (polar pc (angle pc cen) lc)
     
109.           )
     
110.           ;;*******************************************************************
     
111.           ;;画圆---------------------------------------------------------------
     
112.           (setvar "osmode" 0)
     
113.           (command ".line" pa pb pc "C")
     
114.           (command ".CIRCLE" cena ra)
     
115.           (command ".CIRCLE" cenb rb)
     
116.           (command ".CIRCLE" cenc rc)
     
117.           (command ".UCS" "P")
     
118.           (command ".UCS" "P")
     
119.           (setvar "osmode" oldmode)
     
120.           (setvar "cmdecho" oce)
     
121.           (princ)
     
122.         )
     
123.       )
     
124.     )
     
125.   )
     
126. )

<br /> ;;*****************************************************************************<br /> ;;马尔法蒂问题lisp求解---------------------------------------------------------<br /> (defun C:Mal (/           pa        pb   pc          a    b    c         p    xa   xb        xc<br />               ya   yb        yc   aga  agb  agc  ang         ta   tb   tc        ja<br />               jb   jc        ha   hb          hc   vpa  vpb         vpc  cen  ra        rb<br />               rc   la        lb   lc          cena cenb cenc<br />              )<br /> ;;(defun C:Mal ()<br />   (graphscr)<br />   (setq oldmode (getvar "osmode"))<br />   (setq oce (getvar "cmdecho"))<br />   (setvar "cmdecho" 0)<br />   ;;输入数据<br />   (command ".ucs" "W")<br />   (setq pa (getpoint "请输入第一点：\n"))<br />   (setq pb (getpoint "请输入第二点：\n"))<br />   (setq pc (getpoint "请输入第三点：\n"))<br />   (command ".ucs" "O" pa)<br />   (setq        pa (trans pa 0 1)<br />         pb (trans pb 0 1)<br />         pc (trans pc 0 1)<br />   )<br />   ;;求三角形边长和半周长  <br />   (setq        a (distance pb pc)<br />         b (distance pc pa)<br />         c (distance pa pb)<br />         p (/ (+ a b c) 2)<br />   )<br />   ;;计算切线长<br />   (setq        xa (sqrt (abs (- (* p p) (* p a))))<br />         xb (sqrt (abs (- (* p p) (* p b))))<br />         xc (sqrt (abs (- (* p p) (* p c))))<br />   )<br />   (setq        ya (sqrt (* p a))<br />         yb (sqrt (* p b))<br />         yc (sqrt (* p c))<br />   )<br />   (setq        aga (angle '(0 0) (list xa ya))<br />         agb (angle '(0 0) (list xb yb))<br />         agc (angle '(0 0) (list xc yc))<br />   )<br />   (setq ang (/ (+ aga agb agc) 2))<br />   (setq ta (* p (* (sin (- ang aga))) (* (sin (- ang aga)))))<br />   (setq tb (* p (* (sin (- ang agb))) (* (sin (- ang agb)))))<br />   (setq tc (* p (* (sin (- ang agc))) (* (sin (- ang agc)))))<br />   ;;***************************************************************************<br />   ;;求三角形内心---------------------------------------------------------------<br />   (defun cen_incir (pa pb pc)<br />     (setq jc (angle pa pb)<br />           ja (angle pb pc)<br />           jb (angle pc pa)<br />     )<br />     (setq ha (/ (+ jb jc pi) 2)<br />           hb (/ (+ jc ja pi) 2)<br />           hc (/ (+ ja jb pi) 2)<br />     )<br />     (setq vpa (polar pa ha p)<br />           vpb (polar pb hb p)<br />           vpc (polar pc hc p)<br />     )<br />     (inters pa vpa pb vpb nil)<br />   )<br />   (setq cen (cen_incir pa pb pc))<br />   ;;***************************************************************************<br />   ;;求每个圆的半径，圆心位置---------------------------------------------------<br />   (defun tan (x)<br />     (/ (sin x) (cos x))<br />   )<br />   ;;定义正切函数<br />   (if (> 1e-16<br />          (abs (* (sin (- jb jc)) (sin (- jc ja)) (sin (- ja jb))))<br />       )<br />     (progn<br />       (alert "你输入的三点在一条直线上,请重新输入!")<br />       (command ".UCS" "P")<br />       (command ".UCS" "P")<br />       (setvar "osmode" oldmode)<br />       (setvar "cmdecho" oce)<br />       (princ)<br />     )<br />     ;;判断输入的三点是否在同一条直线上<br />     (progn<br />       (setq ra (* ta (/ 1 (abs (tan (/ (- jb jc) 2))))))<br />       (setq rb (* tb (/ 1 (abs (tan (/ (- jc ja) 2))))))<br />       (setq rc (* tc (/ 1 (abs (tan (/ (- ja jb) 2))))))<br />       (if (or nil (< (abs ra) 1e-8) (< (abs ra) 1e-8) (< (abs ra) 1e-8))<br />         (progn<br />           (princ "你输入的三点在一条直线上,请重新输入!")<br />           (command ".UCS" "P")<br />           (command ".UCS" "P")<br />           (setvar "cmdecho" oce)<br />           (setvar "osmode" oldmode)<br />           (princ)<br />         )<br />         ;;判断圆的半径是否为零<br />         (progn<br />           (defun gougu (x y)<br />             (sqrt (+ (* x x) (* y y)))<br />           )<br />           (setq        la (gougu ta ra)<br />                 lb (gougu tb rb)<br />                 lc (gougu tc rc)<br />           )<br />           (setq        cena (polar pa (angle pa cen) la)<br />                 cenb (polar pb (angle pb cen) lb)<br />                 cenc (polar pc (angle pc cen) lc)<br />           )<br />           ;;*******************************************************************<br />           ;;画圆---------------------------------------------------------------<br />           (setvar "osmode" 0)<br />           (command ".line" pa pb pc "C")<br />           (command ".CIRCLE" cena ra)<br />           (command ".CIRCLE" cenb rb)<br />           (command ".CIRCLE" cenc rc)<br />           (command ".UCS" "P")<br />           (command ".UCS" "P")<br />           (setvar "osmode" oldmode)<br />           (setvar "cmdecho" oce)<br />           (princ)<br />         )<br />       )<br />     )<br />   )<br /> )