Stereographic duoprism with Asymptote

DuoprismStereo_ASY.md

DuoprismStereo_ASYMPTOTE.asy

settings.render = 4;
settings.outformat = "eps";
import tube;
size(200,0);

currentprojection=orthographic(4,4,4);
currentlight = light(gray(0.85), ambient=black, specularfactor=3,
                     (100,100,100), specular=gray(0.9), viewport=false);
currentlight.background = rgb("363940ff");

// files to be saved -----------------------------------------------------------
string[] files = {
	"DP000", "DP001", "DP002", "DP003", "DP004", "DP005",
	"DP006", "DP007", "DP008", "DP009", "DP010", "DP011",
	"DP012", "DP013", "DP014", "DP015", "DP016", "DP017",
	"DP018", "DP019", "DP020", "DP021", "DP022", "DP023",
	"DP024", "DP025", "DP026", "DP027", "DP028", "DP029",
	"DP030", "DP031", "DP032", "DP033", "DP034", "DP035",
	"DP036", "DP037", "DP038", "DP039", "DP040", "DP041",
	"DP042", "DP043", "DP044", "DP045", "DP046", "DP047",
	"DP048", "DP049", "DP050", "DP051", "DP052", "DP053",
	"DP054", "DP055", "DP056", "DP057", "DP058", "DP059",
	"DP060", "DP061", "DP062", "DP063", "DP064", "DP065",
	"DP066", "DP067", "DP068", "DP069", "DP070", "DP071",
	"DP072", "DP073", "DP074", "DP075", "DP076", "DP077",
	"DP078", "DP079", "DP080", "DP081", "DP082", "DP083",
	"DP084", "DP085", "DP086", "DP087", "DP088", "DP089",
	"DP090", "DP091", "DP092", "DP093", "DP094", "DP095",
	"DP096", "DP097", "DP098", "DP099", "DP100", "DP101",
	"DP102", "DP103", "DP104", "DP105", "DP106", "DP107",
	"DP108", "DP109", "DP110", "DP111", "DP112", "DP113",
	"DP114", "DP115", "DP116", "DP117", "DP118", "DP119",
	"DP120", "DP121", "DP122", "DP123", "DP124", "DP125",
	"DP126", "DP127", "DP128", "DP129", "DP130", "DP131",
	"DP132", "DP133", "DP134", "DP135", "DP136", "DP137",
	"DP138", "DP139", "DP140", "DP141", "DP142", "DP143",
	"DP144", "DP145", "DP146", "DP147", "DP148", "DP149",
	"DP150", "DP151", "DP152", "DP153", "DP154", "DP155",
	"DP156", "DP157", "DP158", "DP159", "DP160", "DP161",
	"DP162", "DP163", "DP164", "DP165", "DP166", "DP167",
	"DP168", "DP169", "DP170", "DP171", "DP172", "DP173",
	"DP174", "DP175", "DP176", "DP177", "DP178", "DP179"};

// lexicographic order ---------------------------------------------------------
bool dominates(int[] e1, int[] e2){
	return e2[0]>e1[0] || (e2[0]==e1[0] && e2[1]>e1[1]);
}

// vertices --------------------------------------------------------------------
int A = 8;
int B = 4;

struct quadruple {
	real x;
	real y;
	real z;
	real t;
}

real[][] poly1 = new real[A][2];
for(int i = 0; i < A; ++i){
	poly1[i][0] = cos(i/A*2pi);
	poly1[i][1] = sin(i/A*2pi);
}
real[][] poly2 = new real[B][2];
for(int i = 0; i < B; ++i){
	poly2[i][0] = cos(pi/B+i/B*2pi);
	poly2[i][1] = sin(pi/B+i/B*2pi);
}
quadruple[][] vertices = new quadruple[A][B];
for(int i = 0; i < A; ++i){
	for(int j = 0; j < B; ++j){
		quadruple v;
		v.x = poly1[i][0]; v.y = poly1[i][1]; 
		v.z = poly2[j][0]; v.t = poly2[j][1];
		vertices[i][j] = v;
	}
}

// edges -----------------------------------------------------------------------
int[][][] edges;
for(int i = 0; i < A; ++i){
	for(int j = 0; j < B; ++j){
		int[] e = {i,j};
		int[] candidate = {i,(j-1)%B};
		if(dominates(e,candidate)){
			int[][] edge = {e,candidate}; 
			edges.push(edge);
		}
		int[] candidate = {i,(j+1)%B};
		if(dominates(e,candidate)){
			int[][] edge = {e,candidate}; 
			edges.push(edge);
		}
		int[] candidate = {(i-1)%A,j};
		if(dominates(e,candidate)){
			int[][] edge = {e,candidate}; 
			edges.push(edge);
		}
		int[] candidate = {(i+1)%A,j};
		if(dominates(e,candidate)){
			int[][] edge = {e,candidate}; 
			edges.push(edge);
		}
	}
}

// rotation in 4D space (right-isoclinic) --------------------------------------
quadruple rotate4d(real alpha, real beta, real xi, quadruple vec){
	real a = cos(xi);
	real b = sin(alpha)*cos(beta)*sin(xi);
	real c = sin(alpha)*sin(beta)*sin(xi);
	real d = cos(alpha)*sin(xi);
	real p = vec.x;
	real q = vec.y;
	real r = vec.z;
	real s = vec.t;
	quadruple out;
	out.x = a*p - b*q - c*r - d*s;
	out.y = a*q + b*p + c*s - d*r;
	out.z = a*r - b*s + c*p + d*q;
	out.t = a*s + b*r - c*q + d*p;
	return out;
}


// stereographic projection ----------------------------------------------------
triple stereog(quadruple A, real r){
	return (A.x, A.y, A.z) / (r - A.t);
}

// stereographic path ----------------------------------------------------------
path3 stereoPath(quadruple A, quadruple B, real r, int n){
    path3 out;
    for(int i = 0; i <= n; ++i){
    	real t = i/n;
    	quadruple M;
    	real x = (1-t)*A.x + t*B.x;
    	real y = (1-t)*A.y + t*B.y;
    	real z = (1-t)*A.z + t*B.z;
    	real t = (1-t)*A.t + t*B.t;
    	real lg = sqrt(x*x + y*y + z*z + t*t) / r;
    	M.x = x / lg; M.y = y / lg; M.z = z / lg; M.t = t / lg;
        out = out .. stereog(M, r);
    }
    return out;
}

// section transformation ------------------------------------------------------
transform T(path3 p3, real t, int n){
	triple M = relpoint(p3, t/(n/4));
	return scale(length(M)/15);	
}

// bounding box ----------------------------------------------------------------
real f=3, h = 4.5, g = 1.5;
path3 boundingbox = (-h,0,-f)--(-h,0,g)--(h,0,f)--(h,0,-g)--cycle;

// draw the duoprism -----------------------------------------------------------
int n = 100;
real r = sqrt(2);
real alpha = pi/2, beta = 0;

for(int file = 0; file < 180; ++file){
	real xi = 2*file*pi/180;
	picture pic;
	// draw bounding box
	draw(pic, boundingbox, rgb("363940ff")+opacity(0));
	// draw edges
	for(int k = 0; k < 2*A*B; ++k){
		quadruple A = vertices[edges[k][0][0]][edges[k][0][1]];
		quadruple B = vertices[edges[k][1][0]][edges[k][1][1]];
		path3 p3 = 
			stereoPath(rotate4d(alpha, beta, xi, A), rotate4d(alpha, beta, xi, B), r, n);
		transform S(real t){
			return T(p3, t, n);
		}
		draw(pic, tube(p3, unitcircle, S), rgb(139,0,139), render(compression=Low,merge=true));
	}
	// draw vertices
	for(int i = 0; i < A; ++i){
		for(int j = 0; j < B; ++j){
			triple Asg = stereog(rotate4d(alpha, beta, xi, vertices[i][j]), r);
			draw(pic, shift(Asg)*scale3(length(Asg)/10)*unitsphere, purple);
		}
	}
	// add and save picture
	add(pic);
	shipout(files[file], bbox(rgb("363940ff"), FillDraw(rgb("363940ff"))));
	erase();
}