協同產品設計實習第1週任務

網際協同設計資料整合, 網際 2D 程式繪圖回顧

網際協同設計資料整合

利用 Github 倉儲中的組員間 git submodule 設定, 可以進行網際協同設計資料整合.

導入組員設計資料範例: http://mde.tw/2016fallcadp/blog/li-yong-git-submodule-gong-neng-jin-xing-zu-yuan-zi-liao-she-ji-zheng-he.html

導入組員程式範例: https://scrum-1.github.io/2016fallcadp_ag100/blog/dao-ru-zu-yuan-cheng-shi-ce-shi.html

網際 2D 程式繪圖回顧

利用 HTML5 Canvas 與 Brython, 可以進行網際 2D 設計模擬繪圖:

https://developer.mozilla.org/en-US/docs/Web/API/Canvas_API/Tutorial/Transformations

https://www.codeproject.com/Articles/598955/CoordinateplussystemplusinplusHTML-plusCanvas-cpl

http://blog.carbonfive.com/2011/03/31/taming-2d-transforms/

http://blog.carbonfive.com/2011/02/17/visualizing-skillsets-in-html5-canvas-part-1/

利用網際 Python3 程式進行 2D 繪圖:

以上為 canvas1, 座標軸 x 向右為正, y 向下為正.

以上為 canvas2, 座標軸 x 向右為正, y 向下為正.

以上為 canvas3, 座標軸 x 向右為正, y 向上為正.

W1 實習任務

請自行分組每班分為八組, 各組協調後每一位組員均採固定座位就坐, 請各組設法列出各組員座位圖後, 以全班協同方式直接在各組網誌上呈現各組員學號與座位圖, 以 a 班為例, 各組倉儲名稱分別為 2017springcd_ag1~2017springcd_ag8, 各組的網誌中均必須設法呈現全班的電腦教室座位圖. (問題: 有沒有辦法在每週上課後第1堂下課之前, 在各組網誌上呈現當週各學員的出席情況與座次表?)
各組必須準備一個隨身硬碟儲存下載的可攜系統, 並且分別下載 tiny2017_50MB.7z, 以及 tiny2017_1GB.7z, 各組期中簡報時, 必須說明如何從最基本的 50MB 系統逐一納入各類工具得到 1GB 的最終可攜系統. (問題: 各組員會不會從無到有, 自行打造此一可攜程式系統?)
每四組將分配一台電腦當作區網協同伺服主機, 請各組分別指派一名組員負責, 向助教報到, 以便學習如何搭建區域網路上的 Fossil SCM 協同主機, 並負責為各組員建立及管理相關帳號. (問題: 各組有沒有能力自行維護區域網路上的協同產品設計主機?)
本學期課程將會使用 Github, Bitbucket, Vimeo, Youtube, Fossil SCM (由各組自行建立) 與 Onshape 等系統, 請各組員確定已經利用學號作為代號, 擁用各系統的擷取帳號. (問題: 如何呈現階段性的設計實習成果?)
請各組員確認已經會在 Solvespace 與 Onshape 中完成四連桿機構的組立, 並且輸出 stl 格式檔案後, 輸入 V-rep 中. (課程終極目標: 本課程將嘗試從電腦輔助機械設計進入運動模擬, 納入簡單的機電整合與傳動, 最後期望每班協同模擬並列印組立出兩台四足行走機構)
本學期每週上課結束前, 各組均必須直接在 Github Page 中以 Reveal.js 格式, 完成各週的協同實習簡報檔, 其中包含各學員與各組任務執行進度與自評.

請根據上述網際繪圖程式架構, 以 Brython 完成下列四連桿機構的示意繪圖:

其中旋轉軸點為 A 座標為 (x1, y1), 旋轉軸端點為 B 座標為 (x2, y2), 第2連桿端點為 C 座標為 (x3, y3), 第3軸的固定端點則為 D 座標為 (x4, y4).

另外, 以 A 點為起點的主動旋轉桿長為 d1, BC 連桿長為 d2, CD 桿長度為 d3, AD 桿長為 d4, BC 連桿上與旋轉路徑目標點 E 對應的點為 F, 與 B 點的距離為 d5, 與 E 點垂直距離為 d6, 主動旋轉軸的逆時鐘旋轉角度則為 t.

此一平面四連桿的輸入為 x1, y1, x4, y4, d1, d2, d3, d5, d6, 以及 t, 輸出則為 E 點的運動路徑.

利用 Sympy 求解:

from sympy import *
'''
已知四連桿四個關鍵點座標分別為 A (x1, y1), B (x2, y2), C (x3, y3) 與 D (x4, y4)
且 E (x5, y5) 點相關參考x 座標距離為 d5, 而 y座標距離為 d6, 以及輸入角度逆時鐘轉 t 度
以 (x1, y1), (x4, y4), d1, d2, d3, d5, d6 及 t 等 10 個參數作為輸入, 
求 E 點座標 (x5, y5)
假設 AB 連桿長度為 d1, BC 連桿長度為 d2, CD 連桿長度為 d3, AD 距離為 d4
'''
x1, x2, x3, x4, x5 = symbols('x1 x2 x3 x4 x5')
y1, y2, y3, y4, y5 = symbols('y1 y2 y3 y4 y5')
d1, d2, d3, d4, d5, d6, t, t3 = symbols('d1 d2 d3 d4 d5 d6 t t3')
ah, bh, aj, dj, bd, hj, dk, bk = symbols('ah bh aj dj bd hj dk bk')
# angle daj defined as daj
daj, adj, bad, bcd, bdc, bdk = symbols('daj adj bad bcd bdc bdk')
# degree factor
degree, pi = symbols('degree pi')
degree = pi/180.0
# 假設 B 點的絕對 y 座標方向投影點為 H
d1 = sqrt((x1-x2)**2+(y1-y2)**2)
#print(d1)
d2 = sqrt((x2-x3)**2+(y2-y3)**2)
d3 = sqrt((x3-x4)**2+(y3-y4)**2)
d4 = sqrt((x1-x4)**2+(y1-y4)**2)
ah = d1*cos(t)
bh = sqrt(d1**2 - ah**2)
aj = Abs(x4-x1)
dj = Abs(y4-y1)
dk = aj - ah
bk =  bh - dj
t3 = bdc + bdk
# for daj, dj**2 = d4**2+aj**2 -2*d4*aj*cos(daj)
pos = 1
if pos == 1:
    daj = solve(-dj**2+d4**2+aj**2 -2*d4*aj*cos(daj), daj)[0]
else:
    daj = solve(-dj**2+d4**2+aj**2 -2*d4*aj*cos(daj), daj)[1]
#print(daj)
# for adj, aj**2=d4**2+dj**2-2*d4*aj*cos(adj)
if pos == 1:
    adj = solve(-aj**2+d4**2+dj**2-2*d4*aj*cos(adj), adj)[0]
else:
    adj = solve(-aj**2+d4**2+dj**2-2*d4*aj*cos(adj), adj)[0]
#print(adj)
bad = t*degree - daj
# according triangle tad find bd
#bd**2 = d1**2+d4**2-2*d1*d4*cos(bad)
if pos == 1:
    bd = solve(-bd**2+d1**2+d4**2-2*d1*d4*cos(bad), bd)[0]
else:
    bd = solve(-bd**2+d1**2+d4**2-2*d1*d4*cos(bad), bd)[1]
print(bd)

if pos == 1:
    bcd = solve(-bd**2+d2**2+d3**2-2*d2*d3*cos(bcd), bcd)[0]
else:
    bcd = solve(-bd**2+d2**2+d3**2-2*d2*d3*cos(bcd), bcd)[1]

if pos == 1:
    bdk = solve(-bk**2+bd**2+dk**2-2*bd*dk*cos(bdk), bdk)[0]
else:
    bdk = solve(-bk**2+bd**2+dk**2-2*bd*dk*cos(bdk), bdk)[1]

if pos == 1:
    bdc = solve(-d2**2+d3**2+bd**2-2*d3*bd*cos(bdc), bdc)[0]
else:
    bdc = solve(-d2**2+d3**2+bd**2-2*d3*bd*cos(bdc), bdc)[1]
print(t3)

PLAP 三角形符號運算

已知 Point a, Lengeh ac, Angle cab, Point b, 求 c 點座標.

利用 Solvespace 繪製 2D 圖:

利用 Jupyterhub 執行網際符號 (Symbolic) 與數值 (Numerical) 分析運算:

Sympy 1.0 手冊.pdf

PLAP 方程式推導

已知 Point a, Length ac, Angle bac 與 Point b, 求 c 點座標.

Sympy 符號運算程式碼:

#PLAP
from sympy import symbols, sqrt, solve, cos, sin, Abs

# inputs
ax, ay, bx, by, bac, ac = symbols('ax ay bx by bac ac')
# intermediate variables
ab, dab = symbols('ab dab')
ad, bd = symbols('ad bd')
# outputs
cx, cy = symbols('cx cy')
# 從 a, b 點座標求 ab, ad 與 bd
ab = sqrt((ax-bx)**2+(ay-by)**2)
ad = Abs(bx-ax)
bd = Abs(by-ay)
data = solve(-bd**2+ad**2+ab**2-2*ad*ab*cos(dab), dab)
# 第1組解
dab = data[0]
cx = ax+ac*cos(dab+bac)
cy = ay+ac*sin(dab+bac)
print("cx=", cx, "cy=", cy)
# 第二組解
dab = data[1]
cx = ax+ac*cos(dab+bac)
cy = ay+ac*sin(dab+bac)
print("cx=", cx, "cy=", cy)

PLLP 方程式推導

已知 Point a, Length ac, Length cb 與 Point b, 求 c 點座標.

#PLLP
from sympy import symbols, sqrt, solve, cos, sin, Abs

# inputs
ax, ay, bx, by, ac, cb = symbols('ax ay bx by ac cb')
# intermediate variables
ab, dab, bac, degree= symbols('ab dab bac degree')
ad, bd = symbols('ad bd')
# outputs
cx, cy = symbols('cx cy')
# 從 a, b 點座標求 ab
ab = sqrt((ax-bx)**2+(ay-by)**2)
#從三角形已知三邊長, 求解 cx, cy
data = solve([ac**2-((ax-cx)**2+(ay-cy)**2), cb**2-((cx-bx)**2+(cy-by)**2)], [cx, cy])
# 第1組解
print("cx = ", data[0][0])
print("cy = ", data[0][1])
# 第2組解
print("cx = ", data[1][0])
print("cy = ", data[1][1])

數值分析驗證程式碼:

from math import pi, cos, sin, sqrt, acos

radian = 180/pi
degree = pi/180

#PLAP
def plap(ax, ay, ac, bac, bx, by, pos):
    if pos == 0:
        cx= ac*cos(bac - acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ax 
        cy= ac*sin(bac - acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ay
    else:
        cx= ac*cos(bac + acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ax 
        cy= ac*sin(bac + acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ay
    return cx, cy

#PLLP
def pllp(ax, ay, ac, cb, bx, by, pos):
    if pos == 0:
        cx =  -((ay - by)*(-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 - sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx)) + (ac**2 - ax**2 - ay**2 + bx**2 + by**2 - cb**2)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))/(2*(ax - bx)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
        cy =  (-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(-ax + bx))/(2*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
    else:
        cx =  -((ay - by)*(-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx)) + (ac**2 - ax**2 - ay**2 + bx**2 + by**2 - cb**2)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))/(2*(ax - bx)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
        cy =  (-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx))/(2*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
    return cx, cy

ax = -60
ay = 0
bx = 0
by = 0
bac = 50*degree
ac = 30
cd = 50
db = 60
ce = 50
ed = 50
cx, cy = plap(ax, ay, ac, bac, bx, by, 0)
print("cx=", cx, "cy=", cy)
dx, dy = pllp(cx, cy, cd, db, bx, by, 0)
print("dx=", dx, "dy=", dy)
ex, ey = pllp(cx, cy, ce, ed, dx, dy, 0)
print("ex=", ex, "ey=", ey)

結果:

cx= -40.716371709403816 cy= 22.98133329356934

dx= -6.698073034033397 dy= 59.62495968661744

ex= -55.44153371488418 ey= 70.76385733649067

上述符號運算式推導之方程式, 透過數值分析運算與下列圖解法進行比對驗證:

數值分析對應之 Solvespace 圖檔.7z (或 fourbar_numerical_solvespace.7z)

曲柄滑塊機構

曲柄滑塊機構:

曲柄滑塊數值分析驗證:

sol 0 符合上述曲柄滑塊圖中的 e 點座標.

PLPP 方程式推導

利用 sympy 推導 PLPP 方程式:

#PLPP
from sympy import symbols, sqrt, solve

# inputs
bx, by, be, cx, cy, dx, dy = symbols('bx by be cx cy dx dy')
# intermediate variables
cd, m= symbols('cd m')
# outputs
ex, ey = symbols('ex ey')
# e on line cd
cd = sqrt((cx-dx)**2+(cy-dy)**2)
m = (dx-cx)/(dy-cy)
data = solve([be-sqrt((bx-ex)**2+(by-ey)**2), ex-cx-m*(ey-cy)] ,  [ex, ey])
print(data)

PLPP 數值分析驗證:

#plpp numerical
from math import sqrt

def plpp(bx, by, be, cx, cy, dx, dy, sol):
    data = [(((cx - dx)*(bx*cx*cy - bx*cx*dy - bx*cy*dx + bx*dx*dy + by*cy**2 - 2*by*cy*dy + by*dy**2 + cx**2*dy - cx*cy*dx - cx*dx*dy + cy*dx**2 + (-cy + dy)*sqrt(be**2*cx**2 - 2*be**2*cx*dx + be**2*cy**2 - 2*be**2*cy*dy + be**2*dx**2 + be**2*dy**2 - bx**2*cy**2 + 2*bx**2*cy*dy - bx**2*dy**2 + 2*bx*by*cx*cy - 2*bx*by*cx*dy - 2*bx*by*cy*dx + 2*bx*by*dx*dy - 2*bx*cx*cy*dy + 2*bx*cx*dy**2 + 2*bx*cy**2*dx - 2*bx*cy*dx*dy - by**2*cx**2 + 2*by**2*cx*dx - by**2*dx**2 + 2*by*cx**2*dy - 2*by*cx*cy*dx - 2*by*cx*dx*dy + 2*by*cy*dx**2 - cx**2*dy**2 + 2*cx*cy*dx*dy - cy**2*dx**2)) - (cx*dy - cy*dx)*(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2))/((cy - dy)*(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2)), (bx*cx*cy - bx*cx*dy - bx*cy*dx + bx*dx*dy + by*cy**2 - 2*by*cy*dy + by*dy**2 + cx**2*dy - cx*cy*dx - cx*dx*dy + cy*dx**2 + (-cy + dy)*sqrt(be**2*cx**2 - 2*be**2*cx*dx + be**2*cy**2 - 2*be**2*cy*dy + be**2*dx**2 + be**2*dy**2 - bx**2*cy**2 + 2*bx**2*cy*dy - bx**2*dy**2 + 2*bx*by*cx*cy - 2*bx*by*cx*dy - 2*bx*by*cy*dx + 2*bx*by*dx*dy - 2*bx*cx*cy*dy + 2*bx*cx*dy**2 + 2*bx*cy**2*dx - 2*bx*cy*dx*dy - by**2*cx**2 + 2*by**2*cx*dx - by**2*dx**2 + 2*by*cx**2*dy - 2*by*cx*cy*dx - 2*by*cx*dx*dy + 2*by*cy*dx**2 - cx**2*dy**2 + 2*cx*cy*dx*dy - cy**2*dx**2))/(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2)), (((cx - dx)*(bx*cx*cy - bx*cx*dy - bx*cy*dx + bx*dx*dy + by*cy**2 - 2*by*cy*dy + by*dy**2 + cx**2*dy - cx*cy*dx - cx*dx*dy + cy*dx**2 + (cy - dy)*sqrt(be**2*cx**2 - 2*be**2*cx*dx + be**2*cy**2 - 2*be**2*cy*dy + be**2*dx**2 + be**2*dy**2 - bx**2*cy**2 + 2*bx**2*cy*dy - bx**2*dy**2 + 2*bx*by*cx*cy - 2*bx*by*cx*dy - 2*bx*by*cy*dx + 2*bx*by*dx*dy - 2*bx*cx*cy*dy + 2*bx*cx*dy**2 + 2*bx*cy**2*dx - 2*bx*cy*dx*dy - by**2*cx**2 + 2*by**2*cx*dx - by**2*dx**2 + 2*by*cx**2*dy - 2*by*cx*cy*dx - 2*by*cx*dx*dy + 2*by*cy*dx**2 - cx**2*dy**2 + 2*cx*cy*dx*dy - cy**2*dx**2)) - (cx*dy - cy*dx)*(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2))/((cy - dy)*(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2)), (bx*cx*cy - bx*cx*dy - bx*cy*dx + bx*dx*dy + by*cy**2 - 2*by*cy*dy + by*dy**2 + cx**2*dy - cx*cy*dx - cx*dx*dy + cy*dx**2 + (cy - dy)*sqrt(be**2*cx**2 - 2*be**2*cx*dx + be**2*cy**2 - 2*be**2*cy*dy + be**2*dx**2 + be**2*dy**2 - bx**2*cy**2 + 2*bx**2*cy*dy - bx**2*dy**2 + 2*bx*by*cx*cy - 2*bx*by*cx*dy - 2*bx*by*cy*dx + 2*bx*by*dx*dy - 2*bx*cx*cy*dy + 2*bx*cx*dy**2 + 2*bx*cy**2*dx - 2*bx*cy*dx*dy - by**2*cx**2 + 2*by**2*cx*dx - by**2*dx**2 + 2*by*cx**2*dy - 2*by*cx*cy*dx - 2*by*cx*dx*dy + 2*by*cy*dx**2 - cx**2*dy**2 + 2*cx*cy*dx*dy - cy**2*dx**2))/(cx**2 - 2*cx*dx + cy**2 - 2*cy*dy + dx**2 + dy**2))]
    return data[sol]
for i in range(2):
    ex, ey = plpp(bx=-34.71, by=28.18, be=40, cx=0, cy=12, dx=10, dy=30, sol=i)
    print("sol", i, ": ex=", ex, "ey=", ey)

Jansen 多連桿機構運動模擬

下列多連桿機構 k 點座標可以利用上述 PLAP 與 PLLP 方程式疊代求得:

數值分析求得:

kx= -30.806349547073083 ky= -84.02289073817713

Jansen 多連桿機構數值分析驗證程式碼:

from math import pi, cos, sin, sqrt, acos

radian = 180/pi
degree = pi/180

#PLAP
def plap(ax, ay, ac, bac, bx, by, ccw):
    if ccw == 1:
        cx= ac*cos(bac - acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ax 
        cy= ac*sin(bac - acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ay
    else:
        cx= ac*cos(bac + acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ax 
        cy= ac*sin(bac + acos((ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 + abs(ax - bx)**2 - abs(ay - by)**2)/(2*sqrt(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2)*abs(ax - bx)))) + ay
    return cx, cy

#PLLP
def pllp(ax, ay, ac, cb, bx, by, cw):
    if cw == 1:
        cx =  -((ay - by)*(-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 - sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx)) + (ac**2 - ax**2 - ay**2 + bx**2 + by**2 - cb**2)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))/(2*(ax - bx)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
        cy =  (-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(-ax + bx))/(2*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
    else:
        cx =  -((ay - by)*(-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx)) + (ac**2 - ax**2 - ay**2 + bx**2 + by**2 - cb**2)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))/(2*(ax - bx)*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
        cy =  (-ac**2*ay + ac**2*by + ax**2*ay + ax**2*by - 2*ax*ay*bx - 2*ax*bx*by + ay**3 - ay**2*by + ay*bx**2 - ay*by**2 + ay*cb**2 + bx**2*by + by**3 - by*cb**2 + sqrt((-ac**2 + 2*ac*cb + ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2 - cb**2)*(ac**2 + 2*ac*cb - ax**2 + 2*ax*bx - ay**2 + 2*ay*by - bx**2 - by**2 + cb**2))*(ax - bx))/(2*(ax**2 - 2*ax*bx + ay**2 - 2*ay*by + bx**2 + by**2))
    return cx, cy

ax = -38
ay = 0
# b 為原點
bx = 0
by = 0
cx = 0
cy = 7.8
# m 為配合 PLAP 新增固定點
mx = 30
my = 7.8
# dcm ccw 方向角度
dcm = 30*degree
cd = 15
# 三角形 dcm 為 ccw plap d=(a, cd, dcm, m)
dx, dy = plap(cx, cy, cd, dcm, mx, my, ccw=1)
print("dx=", dx, "dy=", dy)
# 三角形 aed 為 cw pllp e=(a, ae, ed, d)
ae = 41.5
ed = 50
ex, ey = pllp(ax, ay, ae, ed, dx, dy, cw=1)
print("ex=", ex, "ey=", ey)
# 三角形 afe 為 cw pllp f=(a, af, fe, e)
af = 40.1
fe = 55.8
fx, fy = pllp(ax, ay, af, fe, ex, ey, cw=1)
print("fx=", fx, "fy=", fy)
# 三角形 dha 為 cw pllp h=(d, dh, ha, a)
dh = 61.9
ha = 39.3
hx, hy = pllp(dx, dy, dh, ha, ax, ay, cw=1)
print("hx=", hx, "hy=", hy)
# 三角形 hgf 為 cw pllp g=(h, hg, gf, f)
hg = 36.7
gf = 39.4
gx, gy = pllp(hx, hy, hg, gf, fx, fy, cw=1)
print("gx=", gx, "gy=", gy)
# 三角形 hkg 為 cw pllp k=(h, hk, kg, g)
hk = 49
kg = 65.7
kx, ky = pllp(hx, hy, hk, kg, gx, gy, cw=1)
print("kx=", kx, "ky=", ky)