Below are some samples of the code I used. I'm two lazy to include everything.

plot3d([r*cos(theta), r*sin(theta), sqrt(r)*cos(theta/2)], r=0..1, theta=0..4*Pi, color= sqrt(r)*sin(theta/2),grid=[20,50]);Since the functions involved in the second set of plots are singular, I used the "view" option to truncate values. For example the plot of the real part of the ℘-function

was given by

P := (x,y) -> WeierstrassP(x + y*I,1,0);For example, the frames of the second animation of an elliptic curve were generated using the function "Gr" below. (Some tweaking had to be done to get the graphs to come out right, which explains some of the strange constants.)

plot3d(Re(P(x,y)), x= 0..8, y=0..8, view=-4..4, grid=[50,50]);

P := (x,y) -> WeierstrassP(x + y*I,1,0);The code for the plot of the the j-function is given below. Unfortunately, Maple doesn't what this function is, so it has to expressed in terms of theta functions (to simplify matters, I'm dropping some constant factors). Also I had to use a cutoff function "trun" so that Maple could handle the plot properly.

PP := (x,y) -> WeierstrassPPrime(x + y*I,1,0);

IR := (theta, z) -> cos(theta)*Re(z) + sin(theta)*Im(z);

Gr1 := theta -> plot3d([Re(P(x,y)), Im(P(x,y)), 0.3*IR(theta,PP(x,y))] , x=0..3.74,y=0..3.74, color=[sin(2*Pi*x/3.74), 0.5, sin(2*Pi*y/3.74)], view=[-1..1,-1..1, -1..1], grid=[35,35]);

L1 := theta -> spacecurve([Re(P(x,1.854)), Im(P(x,1.854)), 0.3*IR(theta,PP(x,1.854))] , x=0..3.74, color=red, thickness=2);

L2 := theta -> spacecurve([Re(P(x,0)), Im(P(x,0)), 0.3*IR(theta, PP(x,0))] , x=1..2.7, color=red, thickness=2);

L3 := theta->spacecurve([Re(P(1.854,y)), Im(P(1.854,y)), 0.3*IR(theta,PP(1.854,y))] , y=0..3.74, color=yellow, thickness=2);

Gr := theta-> display({Gr1(theta) , L1(theta),L2(theta),L3(theta)});

lambda := proc (tau)

local q;

q := exp(I*Pi*tau);

return JacobiTheta2(0, q)^4/JacobiTheta3(0, q)^4

end proc;

J := proc (tau)

local lam;

lam := lambda(tau);

return (1-lam+lam^2)^3/(lam^2*(1-lam)^2)

end proc;

trun := x -> 1000*arctan((1/1000)*x); plot3d(0, x = -1.5 .. 1.5, y = 0.001 .. 0.9, color = trun(Re(evalf(J(x+I*y)))), grid = [180, 180], style = patchnogrid);

#include "colors.inc"The degenerate quartic in the same section:

#include "functions.inc"

#declare T = function {f_torus(x,y,z,0.8,0.18)}

background{White}

camera {

location <0, 4, 4>

look_at 0

angle 36

}

light_source {<-100,200,100> colour rgb 1}

isosurface {

function { T(x-0.7,y,z)*T(x+0.7,y,z) - 0.02}

max_gradient 2

accuracy 0.001

contained_by{sphere{0,2}}

pigment {Green}

finish {phong 1}

}

#include "colors.inc"

camera {

location <0, .1, -65>

look_at 0

angle 30

}

background { color Gray50}

light_source { <300, 300, -1000> White }

#declare HT =difference {

torus {

4, 1

rotate -90*x

pigment { Blue }

}

box {<-5, -5, -1>, <5, 0, 1> }

}

#declare Cyl = cylinder {

<0,8,0>, <0,-8,0>, 1

pigment{Blue}

}

#declare Cyl2 = blob {

threshold .5

cylinder {

<0,8,0>, <0,-8,0>, 2, 1

pigment{Blue}

}

sphere {<1,0,0>, 3, 1 pigment {Blue}}

sphere {<1,-4,0>, 3, 1 pigment {Blue}}

sphere {<1,4,0>, 3, 1 pigment {Blue}}

}

#declare Cyl3 = blob {

threshold .5

cylinder {

<0,4,0>, <0,-4,0>, 2, 1

pigment{Green}

}

sphere {<-1,0,0>, 3, 1 pigment {Green}}

sphere {<-1,-4,0>, 3, 1 pigment {Green}}

sphere {<-1,4,0>, 3, 1 pigment {Green}}

}

union {

object {HT translate 8*y}

object {HT

rotate 180*x

translate -8*y

}

object{Cyl2 translate x*4 }

object{Cyl translate -x*4 }

object{Cyl3 translate x*9.2 }

}

*Purchasesd with funding from the NSF