Ice Nanotube

python code NodeBox nanotube visualization

シンプルにしたが、頂点に○が書けなくなった。頂点・面・辺の描き順はいい解決策がない。

キラルなナノチューブも自在に描けると便利かな。

import svgwrite as sw
from math import *
import numpy as np


def rotmat(x,y,z):
    """
    3Dプリンタ用の造形言語やソフトウェアでよくある、x,y,zの順の回転を実現する。
    ただし、角度はラジアン。
    行列を横ベクトルの後ろから掛ける記法を採用しているので、行列の順はz,y,xとなる。
    """
    Rz = np.array([[cos(z), -sin(z), 0.0],
                   [sin(z), cos(z), 0.0],
                   [0,      0,      1.0]])
    Rx = np.array([[1,      0,      0],
                   [0.0, cos(x), -sin(x)],
                   [0.0, sin(x), cos(x)]])
    Ry = np.array([[cos(y), 0.0, sin(y)],
                   [0,      1,      0],
                   [-sin(y), 0.0, cos(y)]])
    return np.dot(Rx, np.dot(Ry, Rz))


def depth(a,c):
    """
    簡易なパースペクティブ変換
    """
    zoom = exp(a[2]/c)
    return np.array([a[0]*zoom,a[1]*zoom,a[2]])


def draw_poly(svg, vs, **attr):
    """
    SVG直打ちでポリゴンを描く
    """
    group = svg.add( svg.g( id="test" ) )
    p = []
    p.append(["M", vs[0][0], vs[0][1]])
    for v in vs[1:]:
        p.append(["L", v[0], v[1]])
    p.append(["Z"])
    path = sw.path.Path(d=p, **attr)
    group.add(path)


class Prim():
    """
    ポリゴンなどのPrimitiveのクラス。
    verticesは並進も回転も行わない。(最後にまとめて実施)
    """
    def __init__(self, typ, center, vertices=None, **kwargs):
        self.typ = typ
        self.center = center
        self.attr = kwargs
        self.vertices = vertices
        self.rotmat = np.identity(3)
    def rotate(self, x,y,z):
        R = rotmat(x,y,z)
        self.rotmat = np.dot(self.rotmat, R)
        self.center = np.dot(self.center, R)
    def translate(self, X):
        self.center += X
    def draw_svg(self, svg, d, offset=(0,0), zoom=1):
        if self.typ == "p":
            vs = []
            for v in prim.vertices:
                v = np.dot(v, prim.rotmat)
                v += prim.center
                v = depth(v, d)
                v *= zoom
                v += np.array([offset[0], offset[1], 0.0])
                vs.append(v)
            draw_poly(svg, vs, **self.attr)
        elif self.typ == "@":
            v = depth(prim.center, d)
            v *= zoom
            v += np.array([offset[0], offset[1], 0.0])
            svg.add(sw.shapes.Ellipse(center=v[:2], r=[6,6], **self.attr))
             
    


def hextube(r,n,m, **kwargs):
    """
    六角形でおおわれたチューブを描く。
    r: 半径
    n: 周上の六角形の数(チューブの太さ)
    m: 長手方向の六角形の数(チューブの太さ)
    """
    l = r*sin(pi/n)/sqrt(3.0)*2
    hexagons = []
    for ix in range(m):
        for ith in range(n):
            c = r*cos(2*pi*ith/n)
            s = r*sin(2*pi*ith/n)
            ch = r*cos(pi*(2*ith+1)/n)
            sh = r*sin(pi*(2*ith+1)/n)
            c1 = r*cos(pi*(2*ith+2)/n)
            s1 = r*sin(pi*(2*ith+2)/n)
            c1h = r*cos(pi*(2*ith+3)/n)
            s1h = r*sin(pi*(2*ith+3)/n)
            hexagon = np.array([
                (l*(ix*3), c, s),
                (l*(ix*3+1), c, s),
                (l*(ix*3+1.5), ch, sh),
                (l*(ix*3+1), c1, s1),
                (l*(ix*3), c1, s1),
                (l*(ix*3-0.5), ch, sh),])
            com = np.sum(hexagon, axis=0) / hexagon.shape[0]
            hexagon -= com
            hexagons.append(Prim("p", center=com, vertices=hexagon, **kwargs))
            hexagon = np.array([
                (l*(ix*3+1.5), ch, sh),
                (l*(ix*3+2.5), ch, sh),
                (l*(ix*3+3.0), c1, s1),
                (l*(ix*3+2.5), c1h, s1h),
                (l*(ix*3+1.5), c1h, s1h),
                (l*(ix*3+1.0), c1, s1),])
            com = np.sum(hexagon, axis=0) / hexagon.shape[0]
            hexagon -= com
            hexagons.append(Prim("p", center=com, vertices=hexagon, **kwargs))
    return hexagons

def squaretube(r,n,m, **kwargs):
    l = 2*r*sin(pi/n)
    squares = []
    for ix in range(m):
        for ith in range(n):
            c = r*cos(2*pi*ith/n)
            s = r*sin(2*pi*ith/n)
            c1 = r*cos(pi*(2*ith+2)/n)
            s1 = r*sin(pi*(2*ith+2)/n)
            square = np.array([
                (l*(ix), c, s),
                (l*(ix+1), c, s),
                (l*(ix+1), c1, s1),
                (l*(ix), c1, s1),])
            com = np.sum(square, axis=0) / square.shape[0]
            square -= com
            squares.append(Prim("p", center=com, vertices=square, **kwargs))
#    for ix in range(m+1):
#        for ith in range(n):
#            c = r*cos(2*pi*ith/n)
#            s = r*sin(2*pi*ith/n)
#            com = np.array((l*(ix), c, s))
#            squares.append(Prim("@", center=com, **kwargs))
    return squares



# ナノチューブの準備
r=12.8/2
n = 16
hexagons = hextube(r,n,9,
                   stroke_width=1, fill="#fff", stroke="#000", fill_opacity=0.7)
r = 2.75
n = 6
squares  = squaretube(r,n,17,
                      stroke_width=2, fill="#fff", stroke="#000", fill_opacity=0.7)
# 表示方法
theta = 1.3      #ANGLE
d = 50.0         #PERSPECTIVE
offsx = 100
offsy = 150

# SVGの準備
svg = sw.Drawing()

# すべてのプリミティブを原点中心で回転する
prims = hexagons + squares
for prim in prims:
    prim.rotate(theta, theta, 0)

# 回転したあとで、重心でソートし、遠い順に描く。
for prim in sorted(hexagons + squares, key=lambda x: x.center[2]):
    prim.draw_svg(svg, d, zoom=10, offset=(offsx, offsy))

# 標準出力
print(svg.tostring())

svgwriteを使ってpythonでsvgを出力するように変更。

ポリゴンを3次元座標変換してsvgにする汎用ライブラリが欲しい。

import svgwrite as sw
from math import *
import numpy as np

def hextube(r,n,m):
    l = r*sin(pi/n)/sqrt(3.0)*2
    hexagons = []
    for ix in range(m):
        for ith in range(n):
            c = r*cos(2*pi*ith/n)
            s = r*sin(2*pi*ith/n)
            ch = r*cos(pi*(2*ith+1)/n)
            sh = r*sin(pi*(2*ith+1)/n)
            c1 = r*cos(pi*(2*ith+2)/n)
            s1 = r*sin(pi*(2*ith+2)/n)
            c1h = r*cos(pi*(2*ith+3)/n)
            s1h = r*sin(pi*(2*ith+3)/n)
            hexagon = []
            hexagon.append((l*(ix*3), c, s))
            hexagon.append((l*(ix*3+1), c, s))
            hexagon.append((l*(ix*3+1.5), ch, sh))
            hexagon.append((l*(ix*3+1), c1, s1))
            hexagon.append((l*(ix*3), c1, s1))
            hexagon.append((l*(ix*3-0.5), ch, sh))
            hexagons.append(hexagon)
            hexagon = []
            hexagon.append((l*(ix*3+1.5), ch, sh))
            hexagon.append((l*(ix*3+2.5), ch, sh))
            hexagon.append((l*(ix*3+3.0), c1, s1))
            hexagon.append((l*(ix*3+2.5), c1h, s1h))
            hexagon.append((l*(ix*3+1.5), c1h, s1h))
            hexagon.append((l*(ix*3+1.0), c1, s1))
            hexagons.append(hexagon)
    return hexagons

def squaretube(r,n,m):
    l = 2*r*sin(pi/n)
    squares = []
    for ix in range(m):
        for ith in range(n):
            c = r*cos(2*pi*ith/n)
            s = r*sin(2*pi*ith/n)
            c1 = r*cos(pi*(2*ith+2)/n)
            s1 = r*sin(pi*(2*ith+2)/n)
            square = []
            square.append((l*(ix), c, s))
            square.append((l*(ix+1), c, s))
            square.append((l*(ix+1), c1, s1))
            square.append((l*(ix), c1, s1))
            squares.append(square)
    return squares


def outprod(a,b):
    return (a[1]*b[2]-a[2]*b[1],
            a[2]*b[0]-a[0]*b[2],
            a[0]*b[1]-a[1]*b[0])


def normal(polygon):
    #print polygon
    vec = [0.0] * 3
    for i in range(len(polygon)):
        op = outprod(polygon[i-1],polygon[i])
        for i in range(3):
            vec[i] += op[i]
    return vec


def rotatey(a,th):
    c = cos(th)
    s = sin(th)
    return (a[0]*c-a[2]*s,a[1],a[0]*s+a[2]*c)

def rotatex(a,th):
    c = cos(th)
    s = sin(th)
    return (a[0], a[1]*c-a[2]*s, a[1]*s+a[2]*c)

def depth(a,c):
    zoom = exp(a[2]/c)
    return (a[0]*zoom,a[1]*zoom,a[2])


#pos = []
#file = open("001.q.xyz","r")
#lin = file.readline()
#lin = file.readline()
#for lin in file:
#    elem,x,y,z = lin.split()
#    if elem == "O":
#        pos.append((float(x),float(y),float(z)))
#
#bonds = []
#for i in range(len(pos)):
#    xi,yi,zi = pos[i]
#    for j in range(i+1,len(pos)):
#        xj,yj,zj = pos[j]
#        d = (xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2
#        if d < 10:
#            bonds.append((pos[i],pos[j]))
            
# size(600,300)
r=12.8/2
n = 16
hexagons = hextube(r,n,9)
r = 2.75
n = 6
squares  = squaretube(r,n,17)

theta = 1.3      #ANGLE
d = 50.0         #PARSE
offsx = 100
offsy = 150

#backside hexagons
fill = (1,1,1,0.7)
stroke = "#000"
stroke_width = 1
svg = sw.Drawing()
group = svg.add( svg.g( id="test" ) )
for hexagon in hexagons:
    projected = []
    for vertex in hexagon:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        p = []
        p.append(["M", projected[0][0]*10+offsx, projected[0][1]*10+offsy])
        for vertex in projected:
            p.append(["L", vertex[0]*10+offsx, vertex[1]*10+offsy])
        p.append(["Z"])
        path = sw.path.Path(d=p, stroke_width=stroke_width, fill="#fff", stroke=stroke)
        group.add(path)
#strokewidth(2)
#for v0,v1 in bonds:
#    p0 = depth(rotatey(rotatex(v0,theta),theta),d)
#    p1 = depth(rotatey(rotatex(v1,theta),theta),d)
#    line(p0[0]*10+offsx, p0[1]*10+offsy,
#         p1[0]*10+offsx, p1[1]*10+offsy)
#fill(0)
#nostroke()
#for v0,v1 in bonds:
#    p0 = depth(rotatey(rotatex(v0,theta),theta),d)
#    p1 = depth(rotatey(rotatex(v1,theta),theta),d)
#    oval(p0[0]*10+offsx-4, p0[1]*10+offsy-4,8,8)
#    oval(p1[0]*10+offsx-4, p1[1]*10+offsy-4,8,8)


#backside squares
fill = 1,1,1,0.7
stroke = "#000"
stroke_width = 2
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        p = []
        p.append(["M", projected[0][0]*10+offsx, projected[0][1]*10+offsy])
        for vertex in projected:
            p.append(["L", vertex[0]*10+offsx, vertex[1]*10+offsy])
        p.append("Z")
        path = sw.path.Path(d=p, stroke_width=stroke_width, fill="#fff", stroke=stroke, fill_opacity=0.7)
        group.add(path)
        
#backside atoms
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        for vertex in projected:
            u = sw.shapes.Ellipse(center=(vertex[0]*10+offsx, vertex[1]*10+offsy), r=(3,3),
                                 stroke_width=stroke_width, stroke=stroke, fill="#fff", fill_opacity=0.7)
            group.add(u)

#frontside squares
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        p = []
        p.append(["M", projected[0][0]*10+offsx, projected[0][1]*10+offsy])
        for vertex in projected:
            p.append(["L", vertex[0]*10+offsx, vertex[1]*10+offsy])
        p.append("Z")
        path = sw.path.Path(d=p, stroke_width=stroke_width, fill="#fff", stroke=stroke, fill_opacity=0.7)
        group.add(path)

#frontside atoms
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        for vertex in projected:
            u = sw.shapes.Ellipse(center=(vertex[0]*10+offsx, vertex[1]*10+offsy), r=(3,3),
                                  stroke_width=stroke_width, stroke=stroke, fill="#fff",
                                  fill_opacity=0.7)
            group.add(u)

#frontside hexagons
fill = (1,1,1,0.7)
stroke = "#000"
stroke_width = 1
for hexagon in hexagons:
    projected = []
    for vertex in hexagon:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        p = []
        p.append(["M", projected[0][0]*10+offsx, projected[0][1]*10+offsy])
        for vertex in projected:
            p.append(["L", vertex[0]*10+offsx, vertex[1]*10+offsy])
        p.append("Z")
        path = sw.path.Path(d=p, stroke_width=stroke_width, fill="#fff", stroke=stroke, fill_opacity=0.7)
        group.add(path)

print(svg.tostring())

NodeBox用の古いコード。numpyとSVGで書きなおしたいね。 def hextube(r,n,m): l = rsin(pi/n)/sqrt(3.0)2 print l hexagons = [] for ix in range(m): for ith in range(n): c = rcos(2piith/n) s = rsin(2piith/n) ch = rcos(pi(2ith+1)/n) sh = rsin(pi(2ith+1)/n) c1 = rcos(pi(2ith+2)/n) s1 = rsin(pi(2ith+2)/n) c1h = rcos(pi(2ith+3)/n) s1h = rsin(pi(2ith+3)/n) hexagon = [] hexagon.append((l(ix3), c, s)) hexagon.append((l(ix3+1), c, s)) hexagon.append((l(ix3+1.5), ch, sh)) hexagon.append((l(ix3+1), c1, s1)) hexagon.append((l(ix3), c1, s1)) hexagon.append((l(ix3-0.5), ch, sh)) hexagons.append(hexagon) hexagon = [] hexagon.append((l(ix3+1.5), ch, sh)) hexagon.append((l(ix3+2.5), ch, sh)) hexagon.append((l(ix3+3.0), c1, s1)) hexagon.append((l(ix3+2.5), c1h, s1h)) hexagon.append((l(ix3+1.5), c1h, s1h)) hexagon.append((l(ix3+1.0), c1, s1)) hexagons.append(hexagon) return hexagons

def squaretube(r,n,m):
    l = 2*r*sin(pi/n)
    print l
    squares = []
    for ix in range(m):
        for ith in range(n):
            c = r*cos(2*pi*ith/n)
            s = r*sin(2*pi*ith/n)
            c1 = r*cos(pi*(2*ith+2)/n)
            s1 = r*sin(pi*(2*ith+2)/n)
            square = []
            square.append((l*(ix), c, s))
            square.append((l*(ix+1), c, s))
            square.append((l*(ix+1), c1, s1))
            square.append((l*(ix), c1, s1))
            squares.append(square)
    return squares


def outprod(a,b):
    return (a[1]*b[2]-a[2]*b[1],
            a[2]*b[0]-a[0]*b[2],
            a[0]*b[1]-a[1]*b[0])


def normal(polygon):
    #print polygon
    vec = [0.0] * 3
    for i in range(len(polygon)):
        op = outprod(polygon[i-1],polygon[i])
        for i in range(3):
            vec[i] += op[i]
    return vec


def rotatey(a,th):
    c = cos(th)
    s = sin(th)
    return (a[0]*c-a[2]*s,a[1],a[0]*s+a[2]*c)

def rotatex(a,th):
    c = cos(th)
    s = sin(th)
    return (a[0], a[1]*c-a[2]*s, a[1]*s+a[2]*c)

def depth(a,c):
    zoom = exp(a[2]/c)
    return (a[0]*zoom,a[1]*zoom,a[2])


#pos = []
#file = open("001.q.xyz","r")
#lin = file.readline()
#lin = file.readline()
#for lin in file:
#    elem,x,y,z = lin.split()
#    if elem == "O":
#        pos.append((float(x),float(y),float(z)))
#
#bonds = []
#for i in range(len(pos)):
#    xi,yi,zi = pos[i]
#    for j in range(i+1,len(pos)):
#        xj,yj,zj = pos[j]
#        d = (xi-xj)**2 + (yi-yj)**2 + (zi-zj)**2
#        if d < 10:
#            bonds.append((pos[i],pos[j]))
            
size(600,300)
from math import *
r=12.8/2
n = 16
hexagons = hextube(r,n,9)
r = 2.75
n = 6
squares  = squaretube(r,n,17)

theta = 1.3      #ANGLE
d = 50.0         #PARSE
offsx = 100
offsy = 150

#backside hexagons
fill(1,1,1,0.7)
stroke(0)
strokewidth(1)
for hexagon in hexagons:
    projected = []
    for vertex in hexagon:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        beginpath(projected[0][0]*10+offsx, projected[0][1]*10+offsy)
        for vertex in projected:
            lineto(vertex[0]*10+offsx, vertex[1]*10+offsy)
        endpath()
#strokewidth(2)
#for v0,v1 in bonds:
#    p0 = depth(rotatey(rotatex(v0,theta),theta),d)
#    p1 = depth(rotatey(rotatex(v1,theta),theta),d)
#    line(p0[0]*10+offsx, p0[1]*10+offsy,
#         p1[0]*10+offsx, p1[1]*10+offsy)
#fill(0)
#nostroke()
#for v0,v1 in bonds:
#    p0 = depth(rotatey(rotatex(v0,theta),theta),d)
#    p1 = depth(rotatey(rotatex(v1,theta),theta),d)
#    oval(p0[0]*10+offsx-4, p0[1]*10+offsy-4,8,8)
#    oval(p1[0]*10+offsx-4, p1[1]*10+offsy-4,8,8)


#backside squares
fill(1,1,1,0.7)
stroke(0)
strokewidth(2)
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        beginpath(projected[0][0]*10+offsx, projected[0][1]*10+offsy)
        for vertex in projected:
            lineto(vertex[0]*10+offsx, vertex[1]*10+offsy)
        endpath()
#backside atoms
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] > 0:
        for vertex in projected:
            oval(vertex[0]*10+offsx-3, vertex[1]*10+offsy-3,6,6)
#frontside squares
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        beginpath(projected[0][0]*10+offsx, projected[0][1]*10+offsy)
        for vertex in projected:
            lineto(vertex[0]*10+offsx, vertex[1]*10+offsy)
        endpath()
#frontside atoms
for square in squares:
    projected = []
    for vertex in square:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        for vertex in projected:
            oval(vertex[0]*10+offsx-3, vertex[1]*10+offsy-3,6,6)

#frontside hexagons
fill(1,1,1,0.7)
stroke(0)
strokewidth(1)
for hexagon in hexagons:
    projected = []
    for vertex in hexagon:
        projected.append(depth(rotatey(rotatex(vertex,theta),theta),d))
    n = normal(projected)
    if n[2] < 0:
        beginpath(projected[0][0]*10+offsx, projected[0][1]*10+offsy)
        for vertex in projected:
            lineto(vertex[0]*10+offsx, vertex[1]*10+offsy)
        endpath()

Linked from


Edit