SageMath로 Costa 곡면 그리기

양성덕
20210411

목표

Mathias Weber가 작성한 Mathematica Code의 내용을 충실히 따르자면 한 조각을 만든 후 그것에 대칭성을 적용하여 전체를 만든다.
여기서는 한 조각을 만들어 본다.

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

var('rho')

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

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

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

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

x, y = var('x, y', domain='real')

xlist=  [-3.5,-3.0,-2.5, -2.0,-1.5, -1.0,-0.5,-0.3,-0.2,-0.1,-0.05,-0.03,0.0,0.03,0.05,0.1,0.2,0.4,0.6,0.8,1.0,1.2,1.5,2.0 ,2.5,3.0 ,3.5,4.0,5.0,6.0 ];
xlen= len(xlist);


ylist= list(pi*vector([0.001,0.05,0.08,0.16,0.24,0.32,0.4,0.48,0.56,0.64,0.72, 0.8,0.88,0.95,0.98,0.99,1.0]));
ylen=len(ylist);


def X(x):
    if x<0.5 or x>xlen + 0.4:
        return -1
    else:
        return xlist[ int(round(x.real()) -1)]
    
def Y(y):
    if y<0.5 or y>ylen + 0.4:
        return -1
    else:
        return ylist[ int(round(y.real()) -1)]  

Weber가 사용한 위 함수 X, Y는 여기서는 문제를 일으킨다. 여러 가지로 시도해 보다가 결국 다음과 같이 넘어갔다.

f = [lambda x, y: omega1(z).substitute(z = p( xlist[ int(round(x) -1)] + I*ylist[ int(round(y) -1)] )).real(),
     lambda x, y: omega2(z).substitute(z = p( xlist[ int(round(x) -1)] + I*ylist[ int(round(y) -1)] )).real(),
     lambda x, y: omega3(z).substitute(z = p( xlist[ int(round(x) -1)] + I*ylist[ int(round(y) -1)] )).real()
    ]

gr1 = parametric_plot3d( f , (x, 1, xlen),(y, 1, ylen))

gr1.show(frame=False)

OK, 여러 시간 쓴 끝에 한 조각 만들어 내는 데는 성공했다. 자, 이제 이걸 어떻게 변환하나?

아이디어: gr1은 class instance다.

type(gr1)

<class 'sage.plot.plot3d.parametric_surface.ParametricSurface'>

이제 python class의 instance를 분해하는 Python의 일반론을 따라서 이걸 분해하고 자료를 바꾸면 될 것이다.