剧情提要：
[机器小伟]在[工程师阿伟]的陪同下进入了[九转金丹]之第五转的修炼。
这次要研究的是[推理与证明]。
正剧开始：
星历2016年04月24日 12:16:07, 银河系厄尔斯星球中华帝国江南行省。
[工程师阿伟]正在和[机器小伟]一起研究[推理与证明]。
<span style="font-size:18px;">>>>
23.5 13.865424623862042
#海伦-秦九韶公式
def HQFormula(a, b, c):
p = (a+b+c)/2;
S = math.sqrt(p*(p-a)*(p-b)*(p-c));
return S;
#例3
def tmp():
#a, b, c = 3, 4, 5
a, b, c = 3, 4, 5
S1= a*b/2
S2 = a*c/2
S3 = b*c/2
S4 = HQFormula((a*a+b*b)**0.5, (a*a+c*c)**0.5, (b*b+c*c)**0.5);
print(S1+S2+S3, S4);</span>
<span style="font-size:18px;"> if (1) {
var array = new Array();
array[0] = [1];
array[1] = [1,1];
for (var i = 2; i < 10; i++) {
array[i] = [1];
for (var j = 1; j < i; j++) {
array[i].push(array[i-1][j-1]+array[i-1][j]);
}
array[i].push(1);
}
var width = 600, height = 400;
var x = width/2, y = 30;
var measure = 0, s = '';
var len = array.length;
for (var i = 0; i < len; i++) {
s = array[i].join(' ');
measure = plot.measureText(s);
plot.fillText(s, x-measure/2, y, measure);
y += 30;
}
}</span>
关于求二面角，小伟找到了这个公式：
<span style="font-size:18px;">>>>
cos = 0.9218142082600806 角度是： 22.80723538641745 度
#[数] dihedral angle;
#二面角余弦
def dihedral():
xyz = [-1.2999760,0.0173840,-0.7162670, 1.0107690,1.5229620,-0.0945670, 0.0932470,-1.0086090,1.6265070, -1.2280340,-1.4934410,1.8123150, 0.0932470,-1.0086090,1.6265070, 1.0734260,-1.5230970,2.5155820]
#六个点
ax1 = xyz[0]
ay1 = xyz[1]
az1 = xyz[2]
bx2 = xyz[3]
by2 = xyz[4]
bz2 = xyz[5]
cx3 = xyz[6]
cy3 = xyz[7]
cz3 = xyz[8]
dx1 = xyz[9]
dy1 = xyz[10]
dz1 = xyz[11]
ex2 = xyz[12]
ey2 = xyz[13]
ez2 = xyz[14]
fx3 = xyz[15]
fy3 = xyz[16]
fz3 = xyz[17]
#面一法线
nx = ((cz3-az1)/(cy3-ay1)-(bz2-az1)/(by2-ay1))/((bx2-ax1)/(by2-ay1)-(cx3-ax1)/(cy3-ay1))
ny = ((cz3-az1)/(cx3-ax1)-(bz2-az1)/(bx2-ax1))/((by2-ay1)/(bx2-ax1)-(cy3-ay1)/(cx3-ax1))
nz = 1
#面二法线
mx = ((fz3-dz1)/(fy3-dy1)-(ez2-dz1)/(ey2-dy1))/((ex2-dx1)/(ey2-dy1)-(fx3-dx1)/(fy3-dy1))
my = ((fz3-dz1)/(fx3-dx1)-(ez2-dz1)/(ex2-dx1))/((ey2-dy1)/(ex2-dx1)-(fy3-dy1)/(fx3-dx1))
mz = 1
cosAngle = (nx*mx+ny*my+nz*mz)/((math.sqrt(nx**2+ny**2+nz**2))*(math.sqrt(mx**2+my**2+mz**2)));
print('cos = ', cosAngle, '角度是：', 180/math.pi*math.acos(cosAngle), '度');
</span>
但是在验证的过程中小伟发现这个算法好像有问题：
比如对于下面这个例子：
求一下VAB和ABC的二面角
<span style="font-size:18px;">>>>
[2, 5, 2, 5, 0, 1, 0, 0, 0, 1, 1, 5, 5, 0, 1, 0, 0, 0]
cos = 0.4911436350228293 角度是： 60.584222864016645 度
#[数] dihedral angle;
#二面角余弦
#暂时只能算不垂直或平行于xy, xz, yz任一平面的两平面的二面角
#公式寻找中...
def dihedral(points):
#points格式是 []*18, 分六个点，每个点x, y, z坐标排列
'''
xyz = [-1.2999760,0.0173840,-0.7162670, 1.0107690,1.5229620,-0.0945670, 0.0932470,-1.0086090,1.6265070, -1.2280340,-1.4934410,1.8123150, 0.0932470,-1.0086090,1.6265070, 1.0734260,-1.5230970,2.5155820]
'''
if (len(points) != 18):
return 'inf';
xyz = points;
#六个点
ax1 = xyz[0]
ay1 = xyz[1]
az1 = xyz[2]
bx2 = xyz[3]
by2 = xyz[4]
bz2 = xyz[5]
cx3 = xyz[6]
cy3 = xyz[7]
cz3 = xyz[8]
dx1 = xyz[9]
dy1 = xyz[10]
dz1 = xyz[11]
ex2 = xyz[12]
ey2 = xyz[13]
ez2 = xyz[14]
fx3 = xyz[15]
fy3 = xyz[16]
fz3 = xyz[17]
#面一法线
nx = ((cz3-az1)/(cy3-ay1)-(bz2-az1)/(by2-ay1))/((bx2-ax1)/(by2-ay1)-(cx3-ax1)/(cy3-ay1))
ny = ((cz3-az1)/(cx3-ax1)-(bz2-az1)/(bx2-ax1))/((by2-ay1)/(bx2-ax1)-(cy3-ay1)/(cx3-ax1))
nz = 1
#面二法线
mx = ((fz3-dz1)/(fy3-dy1)-(ez2-dz1)/(ey2-dy1))/((ex2-dx1)/(ey2-dy1)-(fx3-dx1)/(fy3-dy1))
my = ((fz3-dz1)/(fx3-dx1)-(ez2-dz1)/(ex2-dx1))/((ey2-dy1)/(ex2-dx1)-(fy3-dy1)/(fx3-dx1))
mz = 1
cosAngle = (nx*mx+ny*my+nz*mz)/((math.sqrt(nx**2+ny**2+nz**2))*(math.sqrt(mx**2+my**2+mz**2)));
print('cos = ', cosAngle, '角度是：', 180/math.pi*math.acos(cosAngle), '度');
#求二面角
def tmp():
V = [2, 5, 2]
A = [5, 0, 1]
B = [0, 0, 0]
C = [1, 1, 5]
points = V+ A + B+ C+A+B;
#print(len(points));
print(points);
dihedral(points);</span>
再求一下VBC与ABC的，或是VAB与VBC的
<span style="font-size:18px;">>>> [2, 5, 2, 0, 0, 0, 1, 1, 5, 1, 1, 5, 5, 0, 1, 0, 0, 0] cos = -0.2558139534883721 角度是： 104.82182106205012 度 >>> ================================ RESTART ================================ >>> [2, 5, 2, 5, 0, 1, 1, 1, 5, 1, 1, 5, 5, 0, 1, 0, 0, 0] cos = -0.19218663979154185 角度是： 101.08042157684446 度</span>
怎么都会是100多度呢，看着不像啊，但小伟也说不准是对是错。
先放着吧。
<span style="font-size:18px;"> if (1) {
var r = 20;
config.setSector(1,1,1,1);
config.graphPaper2D(0, 0, r);
config.axis2D(0, 0,190);
//坐标轴设定
var scaleX = 4*r, scaleY = 4*r;
var spaceX = 2, spaceY = 2;
var xS = -10, xE = 10;
var yS = -10, yE = 10;
config.axisSpacing(xS, xE, spaceX, scaleX, 'X');
config.axisSpacing(yS, yE, spaceY, scaleY, 'Y');
var array = [[2, 5, 2], [5, 0, 1], [0, 0, 0], [1, 1, 5]];
//array = shape.xyzSort(array);
var size = array.length;
var array2D = [];
for (var i = 0; i < size; i++) {
array2D.push(shape.point3D(array[i][0], array[i][1], array[i][2]));
}
//去除重复点
var pointArray = removeDuplicatedPoint(array2D);
//无重复的点的数量
var points = pointArray.length;
//得到距离阵列
//格式为[[点1序号，点2序号, 距离值], ...]
var distanceArray = distanceSort(pointArray);
//边的数量
var edges = distanceArray.length;
//存放需要连通的边
var linkedArray = [];
//连通的边的数量
var links = 0;
//每个顶点相关的边的集合
var edgeOfVertex = [];
for (var i = 0; i < points; i++) {
//获得顶点相关的边的集合
edgeOfVertex = [];
for (var j = 0; j < edges; j++) {
if (distanceArray[j][0] == i ||
distanceArray[j][1] == i) {
edgeOfVertex.push(distanceArray[j]);
}
}
//根据起始点寻找最短长度的两条边
edgeOfVertex.sort(function(a, b) {
return a[2] - b[2];
});
var choice = 4;
if (edgeOfVertex.length > choice) {
edgeOfVertex = edgeOfVertex.slice(0, choice);
}
linkedArray = linkedArray.concat(edgeOfVertex);
}
//document.write(linkedArray.join(' , ')+'<br/>');
linkedArray = removeDuplicatedPoint(linkedArray);
links = linkedArray.length;
//document.write(linkedArray.join(' , ')+'<br/>');
var startPoint, endPoint, x1, y1, x2, y2;
//比例缩放
var scale = 40;
for (var i = 0; i < links; i++) {
startPoint = linkedArray[i][0];
endPoint = linkedArray[i][1];
x1 = pointArray[startPoint][0];
y1 = pointArray[startPoint][1];
x2 = pointArray[endPoint][0];
y2 = pointArray[endPoint][1];
shape.vectorDraw([[x1,y1], [x2, y2]], 'red', scale);
}
shape.pointDraw(pointArray, 'blue', scale, 1, 'VABC');
plot.setFillStyle('blue');
plot.fillText('向量图', -270, -170, 300);
}</span>
本节到此结束，欲知后事如何，请看下回分解。
时间： 2024-10-22 06:22:02