Skip to content

Instantly share code, notes, and snippets.

@lxmfly123

lxmfly123/solve2.py

Last active March 25, 2020 15:02
Show Gist options
  • Save lxmfly123/9603123c6edea077ca8493bb3d4d37bc to your computer and use it in GitHub Desktop.
Save lxmfly123/9603123c6edea077ca8493bb3d4d37bc to your computer and use it in GitHub Desktop.
Download ZIP
For XLL 2
Raw
solve2.py
from scipy.optimize import fsolve
import scipy
import math
import numpy
import sys
C = []
S = []
keyLength = 0.3166
keyLength_b = 0.5483
keyLength_c = 0.2414
keyLength_d = 0.3981
A_1 = 0.75
A_2 = 0.75
_lambda_1 = 10
_lambda_2 = 20
M = 50
N = 50
def getStartAtom():
return [0, 0, 0];
# 令 y = 0，求此边缘曲线上的前 N 个原子
def getAtomsOnEdgeCurve():
C.append([])
n = 0
while n < N:
C[0].append(getAtomOnEdgeCurve(n))
n = n + 1;
return C[0];
# 令 y = 0，求此边缘曲线上的第 n 个原子
def getAtomOnEdgeCurve(n):
if n == 0:
return getStartAtom()
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
y,
math.pow(x - C[0][n - 1][0], 2) + math.pow(y - C[0][n - 1][1], 2) + math.pow(z - C[0][n - 1][2], 2) - math.pow(keyLength, 2)
]
return roundToZero(fsolve(eqs, [C[0][n - 1][0] + 1, 0, 0]).tolist())
# 工具函数，当坐标值应为零时，将求得的极小浮点数转化为零。
def roundToZero(atom):
def toZero(num):
if math.fabs(num) < 0.000000000001:
return 0
return num
return list(map(toZero, atom))
# 根据前行后序号原子和前行同序号原子求得对应原子（需先求出首行所有原子）
def getNthAtomOnMthCurve(n, m):
if m % 2 == 0 and n == 0:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m - 2][n][0], 2) + math.pow(y - C[m-2][n][1], 2) + math.pow(z - C[m - 2][n][2], 2) - math.pow(keyLength_b, 2),
math.pow(x - C[m - 1][n][0], 2) + math.pow(y - C[m - 1][n][1], 2) + math.pow(z - C[m - 1][n][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
];
return fsolve(eqs,[C[m - 2][n][0], C[m - 2][n][1] + keyLength_b, C[m - 1][n][2]]).tolist()
elif m % 2 == 0 and n == N - 1:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m - 2][n][0], 2) + math.pow(y - C[m-2][n][1], 2) + math.pow(z - C[m - 2][n][2], 2) - math.pow(keyLength_b, 2),
math.pow(x - C[m - 1][n - 1][0], 2) + math.pow(y - C[m - 1][n-1][1], 2) + math.pow(z - C[m - 1][n - 1][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
];
return fsolve(eqs,[C[m - 2][n][0], C[m - 2][n][1] + keyLength_b, C[m - 2][n][2]]).tolist()
elif m % 2 == 0:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m - 1][n][0], 2) + math.pow(y - C[m-1][n][1], 2) + math.pow(z - C[m - 1][n][2], 2) - math.pow(keyLength, 2),
math.pow(x - C[m - 1][n - 1][0], 2) + math.pow(y - C[m - 1][n-1][1], 2) + math.pow(z - C[m - 1][n - 1][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
]
return fsolve(eqs,[(C[m - 1][n][0] + C[m - 1][n - 1][0]) / 2, C[m - 1][n][1] + keyLength, C[m - 1][n][2]]).tolist()
else:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m-1][n + 1][0], 2) + math.pow(y - C[m-1][n + 1][1], 2) + math.pow(z - C[m-1][n + 1][2], 2) - math.pow(keyLength, 2),
math.pow(x - C[m - 1][n][0], 2) + math.pow(y - C[m - 1][n][1], 2) + math.pow(z - C[m - 1][n][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
];
return fsolve(eqs,[(C[m - 1][n][0] + C[m - 1][n + 1][0]) / 2, C[m - 1][n][1] + keyLength, (C[m - 1][n][2] + C[m-1][n+1][2]) /2]).tolist()
def getAtomsOnMthCurve(m):
if m == 0:
getAtomsOnEdgeCurve()
else:
r = []
n = 0;
while n < N - m % 2:
r.append(getNthAtomOnMthCurve(n, m))
n = n + 1
C.append(r)
return C[m];
# 根据本行前序号原子和前行同序号原子求得对应原子（需先求出首行所有原子）
def getNthAtomOnMthCurve2(n, m):
if m % 2 == 1:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m][n - 1][0], 2) + math.pow(y - C[m][n - 1][1], 2) + math.pow(z - C[m][n - 1][2], 2) - math.pow(keyLength, 2),
math.pow(x - C[m - 1][n][0], 2) + math.pow(y - C[m - 1][n][1], 2) + math.pow(z - C[m - 1][n][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
];
print(m, n)
return fsolve(eqs,[(C[m - 1][n][0] + C[m - 1][n + 1][0]) / 2, C[m - 1][n][1] + keyLength, C[m - 1][n][2]]).tolist()
else:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
math.pow(x - C[m][n - 1][0], 2) + math.pow(y - C[m][n - 1][1], 2) + math.pow(z - C[m][n - 1][2], 2) - math.pow(keyLength, 2),
math.pow(x - C[m - 1][n-1][0], 2) + math.pow(y - C[m - 1][n-1][1], 2) + math.pow(z - C[m - 1][n-1][2], 2) - math.pow(keyLength, 2),
A_1 * scipy.sin(2 * math.pi / _lambda_1 * y) + A_2 * scipy.sin(2 * math.pi / _lambda_2 * (x / 2 + math.sqrt(3) / 2 * y)) - z,
];
print(m, n)
return fsolve(eqs,[C[m][n-1][0] + keyLength, C[m][n-1][1], C[m][n-1][2]]).tolist()
def getAtomsOnMthCurve2(m):
if m == 0:
getAtomsOnEdgeCurve()
else:
C.append([])
n = 0;
while n < N - m % 2:
if n == 0:
C[m].append(getNthAtomOnMthCurve(n, m))
else:
C[m].append(getNthAtomOnMthCurve2(n, m))
n = n + 1
return C[m];
def scaleCByTen():
global C
C = list(map(lambda row:
list(map(lambda atom:
[atom[0] * 10, atom[1] * 10, atom[2] * 10],
row)),
C))
def scaleSByTen():
global S
S = list(map(lambda row:
list(map(lambda pair:
list(map(lambda atom:
[atom[0] * 10, atom[1] * 10, atom[2] * 10],
pair)),
row)),
S))
def testMthAtoms(m):
testedM = m
n = 0
while n < N - 1 - m % 2:
print(math.sqrt(math.fabs(math.pow(C[testedM][n][0]-C[testedM][n+1][0], 2) + math.pow(C[testedM][n][1]-C[testedM][n+1][1], 2) + math.pow(C[testedM][n][2]-C[testedM][n+1][2], 2) - math.pow(keyLength, 2))))
n = n + 1
def testMthAtoms2(m):
testedM = m
n = 0
while n < N - 1:
print(math.sqrt(math.fabs(math.pow(C[testedM][n][0]-C[testedM - 1][n+1][0], 2) + math.pow(C[testedM][n][1]-C[testedM-1][n+1][1], 2) + math.pow(C[testedM][n][2]-C[testedM - 1][n+1][2], 2) - math.pow(keyLength, 2))))
n = n + 1
def getCornorCAtoms(S_m, S_n):
# print(S_m, S_n)
if S_m % 2 == 1 and S_n == 0:
return [C[S_m - 1][S_n], C[S_m][S_n], C[S_m + 1][S_n]]
if S_m % 2 == 1 and S_n == N - 1:
return [C[S_m][S_n - 1], C[S_m - 1][S_n], C[S_m + 1][S_n]]
if S_m % 2 == 1:
return [C[S_m][S_n - 1], C[S_m][S_n], C[S_m + 1][S_n]]
return [C[S_m][S_n], C[S_m][S_n + 1], C[S_m + 1][S_n]]
def composeEqs2(S_m, S_n):
cornerCAtoms = getCornorCAtoms(S_m, S_n)
if S_m % 2 == 1 and S_n == 0:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
(x - cornerCAtoms[0][0])**2 + (y - cornerCAtoms[0][1])**2 + (z - cornerCAtoms[0][2])**2 - keyLength_d ** 2,
(x - cornerCAtoms[1][0])**2 + (y - cornerCAtoms[1][1])**2 + (z - cornerCAtoms[1][2])**2 - keyLength_c ** 2,
(x - cornerCAtoms[2][0])**2 + (y - cornerCAtoms[2][1])**2 + (z - cornerCAtoms[2][2])**2 - keyLength_c ** 2,
]
return eqs;
if S_m % 2 == 1 and S_n == N - 1:
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
(x - cornerCAtoms[0][0])**2 + (y - cornerCAtoms[0][1])**2 + (z - cornerCAtoms[0][2])**2 - keyLength_c ** 2,
(x - cornerCAtoms[1][0])**2 + (y - cornerCAtoms[1][1])**2 + (z - cornerCAtoms[1][2])**2 - keyLength_d ** 2,
(x - cornerCAtoms[2][0])**2 + (y - cornerCAtoms[2][1])**2 + (z - cornerCAtoms[2][2])**2 - keyLength_c ** 2,
]
return eqs;
# if S_m % 2 == 1:
# def eqs(i):
# x, y, z = i[0], i[1], i[2]
# return [
# (x - cornerCAtoms[0][0])**2 + (y - cornerCAtoms[0][1])**2 + (z - cornerCAtoms[0][2])**2 - keyLength_c ** 2,
# (x - cornerCAtoms[1][0])**2 + (y - cornerCAtoms[1][1])**2 + (z - cornerCAtoms[1][2])**2 - keyLength_c ** 2,
# (x - cornerCAtoms[2][0])**2 + (y - cornerCAtoms[2][1])**2 + (z - cornerCAtoms[2][2])**2 - keyLength_c ** 2,
# ]
# return eqs;
def eqs(i):
x, y, z = i[0], i[1], i[2]
return [
(x - cornerCAtoms[0][0])**2 + (y - cornerCAtoms[0][1])**2 + (z - cornerCAtoms[0][2])**2 - keyLength_c ** 2,
(x - cornerCAtoms[1][0])**2 + (y - cornerCAtoms[1][1])**2 + (z - cornerCAtoms[1][2])**2 - keyLength_c ** 2,
(x - cornerCAtoms[2][0])**2 + (y - cornerCAtoms[2][1])**2 + (z - cornerCAtoms[2][2])**2 - keyLength_c ** 2,
]
return eqs;
def getNthSAtomAtMth(S_n, S_m):
cornerCAtoms = getCornorCAtoms(S_m, S_n)
if S_m % 2 == 1 and S_n == 0:
return [
fsolve(composeEqs2(S_m, S_n),[cornerCAtoms[0][0], cornerCAtoms[0][1] + keyLength_d, cornerCAtoms[0][2] + 10]).tolist(),
fsolve(composeEqs2(S_m, S_n),[cornerCAtoms[0][0], cornerCAtoms[0][1] + keyLength_d, cornerCAtoms[0][2] - 10]).tolist()
]
else:
return [
fsolve(composeEqs2(S_m, S_n),[cornerCAtoms[2][0], (cornerCAtoms[0][1] + cornerCAtoms[2][1]) / 2, cornerCAtoms[0][2] + 10]).tolist(),
fsolve(composeEqs2(S_m, S_n),[cornerCAtoms[2][0], (cornerCAtoms[0][1] + cornerCAtoms[2][1]) / 2, cornerCAtoms[0][2] - 10]).tolist()
]
if __name__ == "__main__":
m = 0
maxM = M
while m < maxM:
getAtomsOnMthCurve(m)
# getAtomsOnMthCurve2(m)
m = m + 1
# testMthAtoms(9)
f = open("C:/Users/truef/Desktop/新建文件夹/r.txt", "w+", encoding="utf-8")
print(C, file = f,)
f.close()
# scaleByTen()
# print(C[maxM - 1])
# print('\n\n\n')
m = 0
maxM = M
while m < maxM - 1:
n = 0
S.append([])
while n < N - ((m + 1) % 2):
S[m].append(getNthSAtomAtMth(n, m))
n = n + 1
# getAtomsOnMthCurve2(m)
m = m + 1
# scaleCByTen()
# scaleSByTen()
f = open("C:/Users/truef/Desktop/新建文件夹/r.txt", "w+", encoding="utf-8")
print(C, file = f,)
f.close()
f = open("C:/Users/truef/Desktop/新建文件夹/r2.txt", "w+", encoding="utf-8")
print(S, file = f)
f.close()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment