ベクトルの垂直「2023 福島大学前期 食農学類【3】」をsympyとFreeCADでやってみたい。

・FreeCADで、図を多く取り入れわかりやすく?　できませんでした。
　作図も下手で申し訳ありません。
　datum planeでできませんでした。 planeでしました。
　planeの隅の4点の座標計算を探しています。

**・２平面の交差？直線の方程式です。
・内積の垂直条件について、
　???　零ベクトル(非ゼロベクトル)が、0個か,1個か,2個か,わかりませんでした。???**


オリジナル

https://kamelink.com/blog/2023/08/31/9-1-23-hukusimadai-syokunou-3

上と同じです。大学入試数学問題集成 様＞テキスト 食農学類【3】

https://mathexamtest.web.fc2.com/2023/202310141/2023101410100mj.html#top-0117

公式HP
https://nyushi.adb.fukushima-u.ac.jp/problem/Files/2023/03/R5_riko_shokuno_zenki.pdf#page=5
＜過去問題（学類入学試験）
https://nyushi.adb.fukushima-u.ac.jp/kakomon_gaku.html

# 内積の垂直条件
???　零ベクトル(非ゼロベクトル)が、0個か,1個か,2個か,わかりませんでした。???

https://qiita.com/mrrclb48z/items/9642ba4caebf5d9c9e72#chatgpt-35%E3%81%A7

# sympyで(数学入試問題 様の方法で)
勉強中。[余弦定理](https://qiita.com/mrrclb48z/items/9642ba4caebf5d9c9e72)。


# sympyで(内積で)
・今後、内積の垂直条件の時は、長さ>0の、判定文が必要になりました。
　いずれ、新作？myDotユーザー定義関数(判定&dotの計算)が必要になるかもしれません。
```python
from sympy import *
var('s'  ,real=True)
var('x,y',real=True)
O=Point( 0, 0,  0)
A=Point(40, 0,  0)
B=Point( 0, s,  0)
H=Point(40,30,120)
#(1)
print("#(1)|HA|>0",Point(H).distance(A)>0)
print("#(1)|HB|>0",Point(H).distance(B)>0)
eq1  =Eq((A-H).dot(B-H),0)
s_ans=solve(eq1,s)[0]
#(2) 
B    =B.subs({s:s_ans})
P    =Point(x,y,0)
print("#(2)|HB|>0",Point(H).distance(B)>0)
print("#(2)|HP|>0",Point(H).distance(P)>0)
eq2  =Eq(( solve( Eq((B-H).dot(P-H),0),y)[0] -y)*12,0)
# 
print("#(1)",O.distance(H))
print("#(2)",s_ans)
print("#(3)",eq2.lhs,"= 0,z = 0")

#(1)|HA|>0 True
#(1)|HB|>0 True
#(2)|HB|>0 True
#(2)|HP|>0 True

#(1) 130
#(2) 510
#(3) x - 12*y - 40 = 0,z = 0
```

# sympyで(結局、planeとplaneの交差で十分でした。)
？点Aが、点Pの軌跡上でした？
点Aが、∠BHA = 90◦,z=0
点Pが、∠BHP = 90◦,z=0 上の軌跡

平面を求めるだけなら、点Hを通る、法線ベクトル(H-B)の平面でした。
？(3)で内積の計算は、どこへいきましたか。？

**途中、空間上の2点を通る直線で、表現しました。**
**結局、最終行。planeとplaneの交差で十分でした。**

class sympy.geometry.plane.Plane(p1, a=None, b=None, **kwargs)
https://docs.sympy.org/latest/modules/geometry/plane.html#sympy.geometry.plane.Plane
equation(x=None, y=None, z=None)
https://docs.sympy.org/latest/modules/geometry/plane.html#sympy.geometry.plane.Plane.equation
intersection(o)
https://docs.sympy.org/latest/modules/geometry/plane.html#sympy.geometry.plane.Plane.intersection

```python
from sympy import *
var('s'    ,real=True)
var('x,y,z',real=True)
O=Point( 0, 0,  0)
A=Point(40, 0,  0)
B=Point( 0, s,  0)
H=Point(40,30,120)
#(1)
print("#(1)|HA|>0",Point(H).distance(A)>0)
print("#(1)|HB|>0",Point(H).distance(B)>0)
eq1  =Eq((A-H).dot(B-H),0)
s_ans=solve(eq1,s)[0]
#(2) 
B    =B.subs({s:s_ans})
P    =Point(x,y,0)
print("#(2)|HB|>0",Point(H).distance(B)>0)
print("#(2)|HP|>0",Point(H).distance(P)>0)
# 
print("#(1)",O.distance(H))
print("#(2)",s_ans)

#(1)|HA|>0 True
#(1)|HB|>0 True
#(2)|HB|>0 True
#(2)|HP|>0 True

#(1) 130
#(2) 510
#(3) 
myPlane0=Plane(O,normal_vector=(0,0,1) )
myPlane1=Plane(H,normal_vector=(H-B) )
myPlane2=(myPlane1.equation()/40).subs({z:P.z})
print("#(3A)",myPlane0.intersection(myPlane1)[0])
print("#(3A)",myPlane0.intersection(myPlane1)[0].equation())
# 
eq2  =Eq(( solve( Eq((B-H).dot(P-H),0),y)[0] -y)*12,0)
myLine=Line(eq2.lhs) 
print("#(3B)","x=",myLine.p1.x*3,"+",(myLine.p1.x-myLine.p2.x)*(-12),"*t",",",
              "y=",myLine.p1.y*3,"+",(myLine.p1.y-myLine.p2.y)*(-12),"*t",",",
              "z=0 (t:実数)"
     )
# 
print("#(3C)",  "(x-",myLine.p1.x*3,")/",(myLine.p1.x-myLine.p2.x)*(-12),
              "= (y-",myLine.p1.y*3,")/",(myLine.p1.y-myLine.p2.y)*(-12),",",  
              "  z=0"
     )

#(1)|HA|>0 True
#(1)|HB|>0 True
#(2)|HB|>0 True
#(2)|HP|>0 True

#(1) 130
#(2) 510

#(3A) Line3D(Point3D(40, 0, 0), Point3D(520, 40, 0))
#(3A) (-x + 12*y + 40, z)
#(3B) x= 0 + 12 *t , y= -10 + 1 *t , z=0 (t:実数)
#(3C) (x- 0 )/ 12 = (y- -10 )/ 1 ,   z=0
```

# FreeCADのマクロで作図
平面作成のマクロ＜document
https://wiki.freecad.org/Part_Plane
```python
# 
import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
### 3D作図
from sympy import *
var('s'    ,real=True)
var('x,y,z',real=True)
O=Point( 0, 0,  0)
A=Point(40, 0,  0)
B=Point( 0, s,  0)
H=Point(40,30,120)
#(1)
print("#(1)|HA|>0",Point(H).distance(A)>0)
print("#(1)|HB|>0",Point(H).distance(B)>0)
eq1  =Eq((A-H).dot(B-H),0)
s_ans=solve(eq1,s)[0]
#(2) 
B    =B.subs({s:s_ans})
P    =Point(x,y,0)
print("#(2)|HB|>0",Point(H).distance(B)>0)
print("#(2)|HP|>0",Point(H).distance(P)>0)
# 
print("#(1)",O.distance(H))
print("#(2)",s_ans)

#(1)|HA|>0 True
#(1)|HB|>0 True
#(2)|HB|>0 True
#(2)|HP|>0 True

#(1) 130
#(2) 510
#(3) 
myPlane0=Plane(O,normal_vector=(0,0,1) )
myPlane1=Plane(H,normal_vector=(H-B) )
myPlane2=(myPlane1.equation()/40).subs({z:P.z})
print("#(3A)",myPlane0.intersection(myPlane1)[0])
print("#(3A)",myPlane0.intersection(myPlane1)[0].equation())
# 
eq2  =Eq(( solve( Eq((B-H).dot(P-H),0),y)[0] -y)*12,0)
myLine=Line(eq2.lhs) 
print("#(3B)","x=",myLine.p1.x*3,"+",(myLine.p1.x-myLine.p2.x)*(-12),"*t",",",
              "y=",myLine.p1.y*3,"+",(myLine.p1.y-myLine.p2.y)*(-12),"*t",",",
              "z=0 (t:実数)"
     )
# 
print("#(3C)",  "(x-",myLine.p1.x*3,")/",(myLine.p1.x-myLine.p2.x)*(-12),
              "= (y-",myLine.p1.y*3,")/",(myLine.p1.y-myLine.p2.y)*(-12),",",  
              "  z=0"
     )

#(1)|HA|>0 True
#(1)|HB|>0 True
#(2)|HB|>0 True
#(2)|HP|>0 True

#(1) 130
#(2) 510

#(3A) Line3D(Point3D(40, 0, 0), Point3D(520, 40, 0))
#(3A) (-x + 12*y + 40, z)
#(3B) x= 0 + 12 *t , y= -10 + 1 *t , z=0 (t:実数)
#(3C) (x- 0 )/ 12 = (y- -10 )/ 1 ,   z=0
##########################################################################################################
### 3D作図
def myXYZ2Txt(A):
    return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
def myTxtXYZ(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 = '40.0 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)
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return 0
def myTxtXYZ_S(*xy_tx):
    for i in range(1,int(len(xy_tx)/2)+1):
        myTxtXYZ(xy_tx[2*i-2],xy_tx[2*i-1]+myXYZ2Txt(xy_tx[2*i-2]) ) 
    return
def myLine_C(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    myLine(args[i],args[0])
    return
#################################################################################
doc = App.activeDocument()
myTxtXYZ_S( A,"A",B,"B",H,"H" )
myLine_S  ( O,H   )
myLine_C  ( A,H,B )
#################################################################################
var('myX,myY,myZ')
def myDegree(myX,myY,myZ):
    r_x=0
    r_y=0
    r_z=0
    if myX==0:
       r_z=0    
    else:
       r_z=atan( myZ/myX ) *180/pi    
    if myY==0:
       r_z=90    
    else:
       r_z=atan( myZ/myY ) *180/pi    

    r_x=float(r_x)
    r_y=float(r_y)
    r_z=float(r_z)
    return r_x,r_y,r_z 
plane = doc.addObject("Part::Plane", "myPlane2")
plane.Length = 200*2
plane.Width = 400*2
plane.Placement = App.Placement(App.Vector(0-H.x*2, 0-H.y*2, 0-H.z*2), App.Rotation( 0, 0,myDegree(0,H.y,H.z)[2] ))
plane = doc.addObject("Part::Plane", "myPlane3")
#
(plx,ply,plz)             =(0-H.x*2 ,0-400 ,0)
(plane_Length,plane_Width)=(200*2,400*3)
plane.Length = plane_Length
plane.Width =  plane_Width
plane.Placement = App.Placement(App.Vector(plx,ply,plz), App.Rotation( 0, 0, 0 ))
##############################################################################
def my_Sphere_name(myCenter,myr,myname):
    sphere = doc.addObject("Part::Sphere", myname)
    sphere.Radius = float(myr)
    sphere.Angle1 = -90
    sphere.Angle2 =  90
    sphere.Angle3 = 360
    sphere.Placement = App.Placement(App.Vector(float(myCenter.x),float(myCenter.y),float(myCenter.z)), App.Rotation(0, 90, 0))
    #FreeCADGui.getDocument('Unnamed').getObject(myname).Transparency = 95
    FreeCADGui.getDocument('Unnamed').getObject(myname).Transparency = 80

my_Sphere_name((A+B)/2, A.distance(B)/2  ,"mySphere_AB")
#################################################################################
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")
```
isometric方向?です。
CADの操作で、transparencyを80にしました。

![0png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/321c4f4c-cb83-8e7d-d7c6-5156d40be2a6.png)

側面から?です。
**側面から見ても、文字が読めるいい方法ないですか？**現在、文字は上向きです。
2方向作図の、1方向非表示ですか。自動的に非表示は、可能ですか？

![1png.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/109698/7355fc47-0cae-3e5e-a877-05bad44172e8.png)

# sympyの実行環境
①私の環境は,pycharmです。
②よく聞くのは、Jupyterです。
③web上で、上記のソースを「SymPy Live shell」に、コピー貼り付けでもできました。
黒背景の右上に、マウスを移動すると、コピーマークが発生します。
???　タブレット環境で、コピー貼り付けが実行できませんでした。???

https://live.sympy.org/

# 参考

以下、いつもの?おすすめです。

https://qiita.com/mrrclb48z/items/75904cdc6e84ecef4d03#%E5%8F%82%E8%80%83wolframalpha-%E9%96%8B%E7%99%BA%E7%89%A9%E8%AA%9E