typedef triple Transform3D(triple t); triple Identity(triple t) { return t; } triple Spherical(triple t) { real x = t.x*cos(t.y)*sin(t.z); real y = t.x*sin(t.y)*sin(t.z); real z = t.x*cos(t.z); return (x,y,z); } triple MaxwellCylindrical(triple t) { real x = 1/pi*(t.x+1+exp(t.x)*cos(t.y)); real y = 1/pi*(t.y +exp(t.x)*sin(t.y)); real z = t.z; return (x,y,z); } triple InvCassiniCylindrical(triple t) { real a = 1; real x = a*2^(1/2)/2*((exp(2*t.x) +2*exp(t.x)*cos(t.y)+1)^(1/2) + exp(t.x)*cos(t.y)+1)^(1/2)/(exp(2*t.x) +2*exp(t.x)*cos(t.y)+1)^(1/2); real y = a*2^(1/2)/2*((exp(2*t.x)+2*exp(t.x) *cos(t.y)+1)^(1/2) - exp(t.x)*cos(t.y)-1)^(1/2)/ (exp(2*t.x)+2*exp(t.x)*cos(t.y)+1)^(1/2); return (x,y,t.z); } triple BipolarCylindrical(triple t) { real h = cosh(t.y)-cos(t.x); if( h == 0) h = 1; real x = sinh(t.y)/ h; real y = sin(t.x) / h; return (x,y,t.z); } triple Ellipsoidal(triple t) { real a = 1, b = 1/2; real x = t.x*t.y*t.z/a/b; real y = abs((t.x^2-b^2)*(t.y^2-b^2)*(b^2-t.z^2)/(a^2-b^2))^(1/2)/b; real z = abs((t.x^2-a^2)*(a^2-t.y^2)*(a^2-t.z^2)/(a^2-b^2))^(1/2)/a; return (x,y,z); } //////////////////////////////////////////////// import palette; import graph3; typedef triple Function3D(pair t); currentlight = adobe; pair H = (-1, 0), L = ( 1, 7); int Size = 250, Grid = 24; projection SetSpecialView() { triple camera = (38.51, 11.61, -5.20); triple up = (-0.0135, -0.005, -0.0197); triple target = (0.0, 0.0, 0.0); return perspective(camera, up, target, autoadjust = false); } projection SetView(string view) { triple v = (0,1,0); // default if(view == "Top") v = ( 0, 0, 1); // +x+y else if(view == "Bottom") v = ( 0, 0,-1); // -x+y else if(view == "Front") v = ( 0, -1,0); // +x+z else if(view == "Back") v = ( 0, 1, 0); // -x+z else if(view == "Right") v = ( 1, 0, 0); // +y+z else if(view == "Left") v = (-1, 0, 0); // -y+z else if(view == "Special") return SetSpecialView(); return orthographic(v.x, v.y, v.z, showtarget = true); } void Plot3D(Function3D f, Transform3D t, string proj, string funname) { write("Plotting: " + funname + proj); size(Size, Size, IgnoreAspect); currentprojection = SetView(proj); triple tf(pair p) { return t(f(p)); } surface s = surface(tf, L, H, Grid); s.colors(palette(s.map(zpart), Rainbow())); draw(s, meshpen = black); // draw(Label("\$x\$",3/4,red),O--X,red,arrow=Arrow3); // draw(Label("\$y\$",3/4,red),O--Y,red,arrow=Arrow3); // draw(Label("\$z\$",3/4,red),O--Z,red,arrow=Arrow3); // dot((0,0,0),red); shipout(funname + proj); erase(); } void Plot3D(Function3D f, string proj, string funname) { Plot3D(f, Identity, proj, funname); } string[] VIEW = { "Top", "Front", "Right", "Bottom", "Back", "Left", "Special" }; void PlotAllViews(Function3D f, string funname) { for(int i = 0; i < VIEW.length; ++i) { Plot3D(f, VIEW[i], funname); } } triple Condor(pair t) { real y = t.y, x = t.x*y; real e = gamma(y+1), ymx = y - x, ypx = y + x; real a = gamma((ymx+1)/2), b = gamma((ymx+2)/2); real c = gamma((ypx+1)/2), d = gamma((ypx+2)/2); real A = cos(pi*ymx), B = cos(pi*ypx); real z = log(e)+log(a)*((A-1)/2)+log(b)*((-A-1)/2) +log(c)*((B-1)/2)+log(d)*((-B-1)/2); return (x,y,z); } Plot3D(Condor, InvCassiniCylindrical, "Front", "CondorInvCass"); Plot3D(Condor, BipolarCylindrical, "Right", "CondorBipolCyl"); Plot3D(Condor, Spherical, "Front", "CondorSph"); Plot3D(Condor, Ellipsoidal, "Front", "CondorEllip"); PlotAllViews(Condor, "Condor"); /********************* MAPLE *********************************** condor := proc(x,y) local a, b, c, d, e, alpha, beta; e := GAMMA(y+1); a := GAMMA(1/2*y-1/2*x+1/2); b := GAMMA(-1/2*x+1/2*y+1); c := GAMMA(1/2*x+1/2*y+1/2); d := GAMMA(1/2*x+1/2*y+1); alpha := cos(Pi*(y-x)); beta := cos(Pi*(y+x)); log(e)+log(a)*((alpha-1)/2)+log(b)* ((-alpha-1)/2)+log(c)*((beta-1)/2)+ log(d)*((-beta-1)/2) end; plot3dMaple := proc(view) local theta, phi; if(view = "top") then theta := -90; phi := 0 # +x+y elif(view = "bottom") then theta := 90; phi := 180 # -x+y elif(view = "front") then theta := -90; phi := 90 # +x+z elif(view = "back") then theta := 90; phi := 90 # -x+z elif(view = "right") then theta := 0; phi := 90 # +y+z elif(view = "left") then theta := 180; phi := 90 # -y+z elif(view = "special") then theta := -145; phi := -105 fi; plot3d(condor(x,y), x = -y..y, y = 0..8, #coords=ellipsoidal, axes=BOXED,orientation = [theta, phi], grid = [124,124]); end; ********************* MAPLE ***********************************/