SageMath로 Costa 곡면 그리기

양성덕
20210411

바로 앞에서 생각하기를 gr instance에 들어있는 자료에 대칭성을 적용하여 전체를 만들 수 있다고 생각했는데 쉬면서 생각해 보니 일단은 함수에 대칭성을 적용시켜서도 전체를 그릴 수 있음을 깨닫게 되었다. 문제는 이러면 여덟 조각을 전부 다 함수를 이용하여 그리는 것이기 때문에 시간이 많이 걸릴 것이라는 것과, 대칭성 때문에 향이 뒤집어지는 것에 대해 아무 것도 할 수 없다는 점이다. 하지만 일단은 그림 전체를 만들어 내는 것이 필요하다 싶어 이 아이디어대로 본다. 결과는, 역시 생각했던 대로 시간이 조금 걸리지만 어쨌거나 Costa 곡면 전체를 얻어 낼 수 있었다.

참고 문헌

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

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

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)

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

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

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

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

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

var('y')

gr1 = parametric_plot3d( (X,Y,Z) , (x, 1, xlen),(y, 1, ylen))

gr2 = parametric_plot3d( (X2,Y2,Z2) , (x, 1, xlen),(y, 1, ylen))

gr12 = parametric_plot3d( (mX,Y,Z) , (x, 1, xlen),(y, 1, ylen))
gr13 = parametric_plot3d( (mX,mY,Z) , (x, 1, xlen),(y, 1, ylen))
gr14 = parametric_plot3d( (X,mY,Z) , (x, 1, xlen),(y, 1, ylen))

gr22 = parametric_plot3d( (mX2,Y2,Z2) , (x, 1, xlen),(y, 1, ylen))
gr23 = parametric_plot3d( (mX2,mY2,Z2) , (x, 1, xlen),(y, 1, ylen))
gr24 = parametric_plot3d( (X2,mY2,Z2) , (x, 1, xlen),(y, 1, ylen))

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

gr.show(frame=False)