共通垂線(Skew lines)「ねじれの位置にある2直線間の最短距離」をsympyとVBAソルバーとFreeCADでやってみたい。

オリジナル

https://vigne-cla.com/5-4/


# SymPy Liveで(オリジナルの偏微分と同じ)最後に値を代入
????教えて下さい。
symbolの配列の使い方を教えて下さい。p0 p1 p2 を配列にできますか？
よろしくお願いします。


以下のソースコードを貼り付けて、Evaluateです。

https://live.sympy.org/

```python
from sympy import *
var('p0 p1 p2')
var('v0 v1 v2')
var('q0 q1 q2')
var('w0 w1 w2')
s = Symbol('s')
t = Symbol('t')
PQ2 = ( (q0+t*w0) - (p0+s*v0) )**2\
     +( (q1+t*w1) - (p1+s*v1) )**2\
     +( (q2+t*w2) - (p2+s*v2) )**2
dPQ2_ds = diff(PQ2, s)
dPQ2_dt = diff(PQ2, t)
ans = solve([dPQ2_ds, dPQ2_dt],[s,t])
s, t = ans[s], ans[t]
xs, ys, zs = p0+s*v0, p1+s*v1, p2+s*v2
xt, yt, zt = q0+t*w0, q1+t*w1, q2+t*w2
PQ2 = ( (q0+t*w0) - (p0+s*v0) )**2\
     +( (q1+t*w1) - (p1+s*v1) )**2\
     +( (q2+t*w2) - (p2+s*v2) )**2
PQ = PQ2**0.5
print("#",PQ)
#
mySubs={p0:-1,p1:1,p2:1,v0:1,v1:0,v2:-1,q0:0,q1:3,q2:2,w0:1,w1:2,w2:-1}
print("#",'On line 1 : (x, y, z) = ({}, {}, {})'.format(xs.subs(mySubs), ys.subs(mySubs), zs.subs(mySubs)))
print("#",'On line 2 : (x, y, z) = ({}, {}, {})'.format(xt.subs(mySubs), yt.subs(mySubs), zt.subs(mySubs)))
print("#",'min distance = {}'.format(PQ.subs(mySubs)))
# ((-p0 + q0 - v0*(-p0*v0*w1**2 - ..........   (以下省略)
# On line 1 : (x, y, z) = (-2, 1, 2)
# On line 2 : (x, y, z) = (-1, 1, 3)
# min distance = 1.41421356237310
```

# SymPy Liveで(外積と内積で最短距離だけ)
参考

https://www.geisya.or.jp/~mwm48961/linear_algebra/line_plane2.htm#six


1.3　SymPyによる内積の計算＜

https://ictsr4.com/l0120.html

```python
from sympy import *
var('p1_x,p1_y,p1_z')
var('p2_x,p2_y,p2_z')
var('v1_x,v1_y,v1_z')
var('v2_x,v2_y,v2_z')
p1 =Matrix([p1_x,p1_y,p1_z])
p2 =Matrix([p2_x,p2_y,p2_z])
v1 =Matrix([v1_x,v1_y,v1_z])
v2 =Matrix([v2_x,v2_y,v2_z])
d=abs((p2-p1).dot(v1.cross(v2))/v1.cross(v2).norm())
print("#",d)
mySubs={p1_x:-1,p1_y:1,p1_z:1,v1_x:1,v1_y:0,v1_z:-1, \
        p2_x: 0,p2_y:3,p2_z:2,v2_x:1,v2_y:2,v2_z:-1}
print("#",d.subs(mySubs),float(d.subs(mySubs)))
# Abs((p1_x - p2_x)*(v1_y*v2_z - v1_z*v2_y) -  (以下省略)
# sqrt(2) 1.4142135623730951
```

１行で書いてみました。値を上記の計算式のprint d(公式?)に代入です。

```python
from sympy import *
print("#",\
  abs((Matrix([0,3, 2]) -     Matrix([-1,1,1])).dot(Matrix([1,0,-1]).cross(Matrix([1,2,-1])))  \
      /Matrix([1,0,-1]).cross(Matrix([1,2,-1])).norm())                                        \
     )
# sqrt(2)
```


# SymPy Liveで(sympy3D)
勉強中。これが本命案です。
(20221027)対応済み
```python
from sympy import *
def myNejire(A,vA,B,vB):
    p1 = Matrix([ A.x, A.y, A.z])
    p2 = Matrix([ B.x, B.y, B.z])
    v1 = Matrix([vA.x,vA.y,vA.z])
    v2 = Matrix([vB.x,vB.y,vB.z])
    return abs((p2 - p1).dot(v1.cross(v2)) / v1.cross(v2).norm())
print("#",myNejire(Point(-1,1,1),Point(1,0,-1),Point(0,3,2),Point(1,2,-1)))
# sqrt(2)
```

# SymPy Liveで(ソルバー)最短距離
勉強中

# vbaで(ソルバー)最短距離

????教えて下さい。
行数を少なくする方法を教えて下さい。
よろしくお願いします。



```vb
Const p0 = -1
Const p1 = 1
Const p2 = 1
Const v0 = 1
Const v1 = 0
Const v2 = -1
Const q0 = 0
Const q1 = 3
Const q2 = 2
Const w0 = 1
Const w1 = 2
Const w2 = -1
Function myR1C1toA1(i, j)
    myR1C1toA1 = Application.ConvertFormula("R" & i & "C" & j, xlR1C1, xlA1)
End Function
Function myKyotuSuisen(p0, p1, p2, v0, v1, v2, q0, q1, q2, w0, w1, w2)
     my0 = myR1C1toA1(2, 1) & "+" & myR1C1toA1(2, 7) & "*" & myR1C1toA1(2, 4) & "-" & myR1C1toA1(1, 1) & "-" & myR1C1toA1(1, 7) & "*" & myR1C1toA1(1, 4)
     my1 = myR1C1toA1(2, 2) & "+" & myR1C1toA1(2, 7) & "*" & myR1C1toA1(2, 5) & "-" & myR1C1toA1(1, 2) & "-" & myR1C1toA1(1, 7) & "*" & myR1C1toA1(1, 5)
     my2 = myR1C1toA1(2, 3) & "+" & myR1C1toA1(2, 7) & "*" & myR1C1toA1(2, 6) & "-" & myR1C1toA1(1, 3) & "-" & myR1C1toA1(1, 7) & "*" & myR1C1toA1(1, 6)
     Range(myR1C1toA1(3, 1)).Formula = "= (" & my0 & ")^ 2 + (" & my1 & ")^ 2 +(" & my2 & ")^ 2"
     Dim ws As Worksheet: Set ws = ActiveSheet
     SolverReset
     SolverOk setCell:=ws.Range(myR1C1toA1(3, 1)), _
                   MaxMinVal:=2, _
                   ByChange:=ws.Range(myR1C1toA1(1, 7) & ":" & myR1C1toA1(2, 7)), _
                   EngineDesc:="GRG Nonlinear"
     SolverOptions AssumeNonNeg:=False
     SolverSolve UserFinish:=True
End Function
Sub myFormat()
    Range("A1:G5").Select
    Selection.NumberFormatLocal = "0.000_ "
    Range("A6").Select
End Sub
Sub aaa_Main()
    Dim myTS As Variant
    ActiveSheet.Cells.Clear
    Cells(1, 1) = p0
    Cells(1, 2) = p1
    Cells(1, 3) = p2
    Cells(1, 4) = v0
    Cells(1, 5) = v1
    Cells(1, 6) = v2
    Cells(2, 1) = q0
    Cells(2, 2) = q1
    Cells(2, 3) = q2
    Cells(2, 4) = w0
    Cells(2, 5) = w1
    Cells(2, 6) = w2
    myTS = myKyotuSuisen(p0, p1, p1, v0, v1, v2, q0, q1, q2, w0, w1, w2)
    Range(myR1C1toA1(4, 1)).Formula = "=" & myR1C1toA1(1, 1) & " + " & myR1C1toA1(1, 7) & " * " & myR1C1toA1(1, 4)
    Range(myR1C1toA1(4, 2)).Formula = "=" & myR1C1toA1(1, 2) & " + " & myR1C1toA1(1, 7) & " * " & myR1C1toA1(1, 5)
    Range(myR1C1toA1(4, 3)).Formula = "=" & myR1C1toA1(1, 3) & " + " & myR1C1toA1(1, 7) & " * " & myR1C1toA1(1, 6)
    Range(myR1C1toA1(5, 1)).Formula = "=" & myR1C1toA1(2, 1) & " + " & myR1C1toA1(2, 7) & " * " & myR1C1toA1(2, 4)
    Range(myR1C1toA1(5, 2)).Formula = "=" & myR1C1toA1(2, 2) & " + " & myR1C1toA1(2, 7) & " * " & myR1C1toA1(2, 5)
    Range(myR1C1toA1(5, 3)).Formula = "=" & myR1C1toA1(2, 3) & " + " & myR1C1toA1(2, 7) & " * " & myR1C1toA1(2, 6)
    myFormat
    Range(myR1C1toA1(6, 1)).Formula = "=sqrt(" & myR1C1toA1(3, 1) & ")"
End Sub
'-1.000  1.000   1.000   1.000   0.000   -1.000  -1.000
' 0.000  3.000   2.000   1.000   2.000   -1.000  -1.000
' 2.000
'-2.000  1.000   2.000
'-1.000  1.000   3.000
' 1.414213562
```

# FreeCADのマクロで
上のsympy計算に、FreeCAD作図を追加しました。
```python
# 上のsympy計算に、FreeCAD作図を追加しました。
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
from sympy import *
var('p0 p1 p2')
var('v0 v1 v2')
var('q0 q1 q2')
var('w0 w1 w2')
s = Symbol('s')
t = Symbol('t')
PQ2 = ( (q0+t*w0) - (p0+s*v0) )**2\
     +( (q1+t*w1) - (p1+s*v1) )**2\
     +( (q2+t*w2) - (p2+s*v2) )**2
dPQ2_ds = diff(PQ2, s)
dPQ2_dt = diff(PQ2, t)
ans = solve([dPQ2_ds, dPQ2_dt],[s,t])
s, t = ans[s], ans[t]
xs, ys, zs = p0+s*v0, p1+s*v1, p2+s*v2
xt, yt, zt = q0+t*w0, q1+t*w1, q2+t*w2
PQ2 = ( (q0+t*w0) - (p0+s*v0) )**2\
     +( (q1+t*w1) - (p1+s*v1) )**2\
     +( (q2+t*w2) - (p2+s*v2) )**2
PQ = PQ2**0.5
print("#",PQ)
mySubs={p0:-1,p1:1,p2:1,v0:1,v1:0,v2:-1,q0:0,q1:3,q2:2,w0:1,w1:2,w2:-1}
print("#",'On line 1 : (x, y, z) = ({}, {}, {})'.format(xs.subs(mySubs), ys.subs(mySubs), zs.subs(mySubs)))
print("#",'On line 2 : (x, y, z) = ({}, {}, {})'.format(xt.subs(mySubs), yt.subs(mySubs), zt.subs(mySubs)))
print("#",'min distance = {}'.format(PQ.subs(mySubs)))
# ((-p0 + q0 - v0*(-p0*v0*w1**2 - ..........   (以下省略)
# On line 1 : (x, y, z) = (-2, 1, 2)
# On line 2 : (x, y, z) = (-1, 1, 3)
# min distance = 1.41421356237310
##########################################################################################################
# 3D作図
def myPoint3D_subs(myA,mySubs):
    return Point3D(myA.x.subs(mySubs),myA.y.subs(mySubs),myA.z.subs(mySubs))
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxt_XYZ(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=float(A.z)
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.2 mm'
def myLine(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),float(A.z)
    Bx,By,Bz=float(B.x),float(B.y),float(B.z)
    pl = FreeCAD.Placement()
    pl.Rotation.Q = (0.4247081540122249, 0.17592004639554645, 0.33985110062924484, 0.8204732460821097)
    pl.Base = FreeCAD.Vector(-3.9166066876399563, -2.1670824762243774, 1.7495260956243028)
    points = [FreeCAD.Vector(Ax,Ay,Az), FreeCAD.Vector(Bx,By,Bz)]
    line = Draft.make_wire(points, placement=pl, closed=False, face=True, support=None)
    Draft.autogroup(line)
myL=-3
SPA=myPoint3D_subs(Point3D(p0       ,p1       ,p2       ),mySubs)
SPB=myPoint3D_subs(Point3D(p0+v0*myL,p1+v1*myL,p2+v2*myL),mySubs)
SQA=myPoint3D_subs(Point3D(q0       ,q1       ,q2       ),mySubs)
SQB=myPoint3D_subs(Point3D(q0+w0*myL,q1+w1*myL,q2+w2*myL),mySubs)
LP =  Point3D     (xs.subs(mySubs), ys.subs(mySubs), zs.subs(mySubs))
LQ =  Point3D     (xt.subs(mySubs), yt.subs(mySubs), zt.subs(mySubs))

myTxt_XYZ(SPA,"SPA"+myXYZ2Txt(SPA))
myTxt_XYZ(SPB,"SPB"+myXYZ2Txt(SPB))
myTxt_XYZ(SQA,"SQA"+myXYZ2Txt(SQA))
myTxt_XYZ(SQB,"SQB"+myXYZ2Txt(SQB))
myTxt_XYZ(LP ,"LP" +myXYZ2Txt(LP ))
myTxt_XYZ(LQ ,"LQ" +myXYZ2Txt(LQ ))
myLine(SPA,SPB)
myLine(SQA,SQB)
myLine(LP ,LQ )
doc = App.activeDocument()
App.ActiveDocument.addObject("App::Origin", "Origin")
App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")
```
ねじれの位置に見えますか？
![png1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/e2c7e61c-7e11-902c-cdf6-17bb10598684.png)

# 参考

https://ja.wikipedia.org/wiki/%E3%81%AD%E3%81%98%E3%82%8C%E3%81%AE%E4%BD%8D%E7%BD%AE

https://en.wikipedia.org/wiki/Skew_lines

https://www.wolframalpha.com/widgets/gallery/view.jsp?id=a7c7f3e02b7060d65396c823de265b0f&reportprob=1