For XLL 2

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()