SageMath로 Costa surface 그리기

양성덕
20210410

거의 20년 전 내가 Mathematica를 잘 알아야겠다고 결심하게 된 계기 중 하나는 Weber가 Costa 곡면을 그릴 때 사용한 프로그램이 Mathematica임을 알게 된 것이었다. 그 후로 정말 여러 방법으로 Costa 곡면을 그려 보았다.

그런데 최근 들어 내 교육과 연구에 SageMath를 활용하기로 작정하면서 제일 먼저 든 궁금증이 Costa 곡면을 SageMath로 그릴 수 있는가 하는 것이었다. 그걸 그릴 수 있다면 내가 연구하는 많은 부분을 학생들에게 소개하는 데에 참으로 좋은 동기가 될 것으로 생각한다.

오늘은 토요일인데 집에서 이것 저것 하다가 문득 구글에서 "How to draw Costa surface with sagemath"를 검색하여 보았는데 2019년에 인터넷에 올려진 코드가 있더라! ([1] 참조) 그걸 돌려 보고 나서는 거의 20년 전에 참조했던 Weber의 Mathematica 파일에 있는 알고리즘을 거의 그대로 활용하여 SageMath로 Costa 곡면을 그려보았다. 물론 Weber와 Mathematica의 조합만큼 좋은 작품이 나올 리가 없겠지만, 어쨌거나 수학적으로 가장 핵심적인 부분은 그럭저럭 이야기할 수 있는 정도는 된 것 같다.

나중에, 좀더 SageMath를 잘 다루게 되면 다시 그려 올려놓을 생각이다. 많이 부족한 그림이지만 학생들이 직접 따라 그려보면서 이것저것 실험해 보기를 희망하며 부족한 그림 올려놓는다.

참고로, 혹시 그림이 까맣게 보인다면 Naver Whale, Mac의 사파리로 봐 보기 바란다. 나는 Chrome, Firefox로는 그림이 까맣게 보이는데 Safari로는 잘 보인다. Chrome으로 잘 안 보이는 것이 Naver Whale로는 잘 보이는 게 이번이 두 번째 경우인데 국산제품이 이렇게 우수한 성능을 보이는 데에 조금 뿌듯함을 느낀다. 그러나 Chrome이 최근 보안 정책을 강화하면서 안 돌아가는 프로그램들이 조금 늘어난 것 같은데 설마 구글이 할 수가 없어서 못 보여 주겠는가, 보안 정책 때문 아닌가 하고 생각해 본다.

Kilom691 의 코드

[1] https://commons.wikimedia.org/wiki/File:Costa_minimal_surface.stl

var('u, v')

e1 = CBF(1/2).elliptic_p(CBF(I))
ii = CBF(I)

def fx(u, v):
    z = CBF(u) + ii*CBF(v)
    z0 = CBF(z).elliptic_zeta(ii)
    z1 = CBF(z - 1/2).elliptic_zeta(ii)
    zi = CBF(z - ii/2).elliptic_zeta(ii)
    return float(1/2*(-z0+pi*CBF(u)+pi**2/4/e1+pi/(2*e1)*(z1-zi)).real_part())

def fy(u, v):
    z  = CBF(u) + ii*CBF(v)
    z0 = CBF(z).elliptic_zeta(ii)
    z1 = CBF(z-1/2).elliptic_zeta(ii)
    zi = CBF(z-ii/2).elliptic_zeta(ii)
    return float(1/2*(-ii*z0+pi*CBF(v)+pi**2/4/e1-pi/(2*e1)*(ii*z1-ii*zi)).real_part())

def fz(u, v):
    ell = (CBF(u)+ii*CBF(v)).elliptic_p(ii)
    return float(sqrt(CBF(2*pi))/4 * abs((ell-e1) / (ell+e1)).log())

P = parametric_plot3d((fx, fy, fz), (u, 0, 1), (v, 0, 1), xmin=-10, xmax=10, zmin=-3, zmax=3)

P.show()

나의 코드

Mathias Weber가 작성한 Mathematica Code의 내용을 sagemath로 나타내 본다. 실제 해 보니 최소한 수학적으로 필요한 함수는 sagemath에 다 들어 있다고 보인다.

[2] Matthias Weber, Costa's Minimal Surface, Mathematica Notebook File

var('rho')

rho = ( gamma(3/4)/sqrt(2)/gamma(5/4) ).n()
rho

0.955977594972250

Phi1(w) = 2*rho*I*sqrt(w)*hypergeometric((1/4,3/2,), (5/4,), w^2)
Phi2(w) = -2/(3*rho)*I*w*sqrt(w)*hypergeometric((1/2,3/4,), (7/4,), w^2)

omega1(z) = -1/2*( Phi1(z) - Phi2(z))
omega2(z) =  I/2*( Phi1(z) + Phi2(z))
omega3(z) =  1/2*log(-1+z) - 1/2*log(1+z)

F(z) = ( omega1(z).real(), omega2(z).real(), omega3(z).real() ) 

Mathematica의 계산결과와 비교하면서 이 식과 계산이 맞음을 확인한다.

F(1.1)

(-1.43407999876598, 0.000000000000000, -1.52226121886171)

F(I).n()

(0.762019208769613, -0.762019208769614, 0.000000000000000)

p(w) = sqrt( 1 + exp(w) )

var('y')

y

points = [(p(x+I*y).real(), p(x+I*y).imag()) for x in srange(-1,3,0.1) for y in srange(0*pi, pi, 0.05*pi, include_endpoint=True)] 

list_plot(points, aspect_ratio=1)

참고로, 다음과 같이 하면 안 그려진다. 이유가 뭘까?

# parametric_plot3d( (p(x+I*y).real(), p(x+I*y).imag(), 0), (x, -1,3), (y, 0, pi))

gr1 = [ F(p(x+I*y)).n() 
      for x in srange(-1,3,0.1, include_endpoint=True) 
      for y in srange(0.05*pi, pi, 0.05*pi, include_endpoint=True)
     ]

a = point3d(gr1)
a.show(frame=False, axes=True)

def sym2 (point3d): x,y,z = point3d; return (y,x,-z)

def sym12(point3d): x,y,z = point3d; return (-x,y,z)
def sym13(point3d): x,y,z = point3d; return (-x,-y,z)
def sym14(point3d): x,y,z = point3d; return (x,-y,z)

def sym22(point3d): x,y,z = point3d; return (-x,y,z)
def sym23(point3d): x,y,z = point3d; return (-x,-y,z)
def sym24(point3d): x,y,z = point3d; return (x,-y,z)

gr2  = list( map( sym2,   gr1 ) )

gr12 = list( map( sym12,  gr1 ) )
gr13 = list( map( sym13,  gr1 ) )
gr14 = list( map( sym14,  gr1 ) )

gr22 = list( map( sym22,  gr2 ) )
gr23 = list( map( sym23,  gr2 ) )
gr24 = list( map( sym24,  gr2 ) )

gr = gr1 + gr2 + gr12 + gr13 + gr14 + gr22 + gr23 + gr24

point3d(gr, frame=False)