11

I would like to produce an Asymptote version of the scene portrayed here: https://www.youtube.com/watch?v=kE6aC51LMv0 The output should preferably be vector graphics. Nevertheless, this is how far I got:

//import settings;
settings.outformat = "png";
settings.render = 32;
//settings.antialias = 32;
settings.prc = false;

import graph3;

size (6cm, 0);
//triple cam = (50, 70, 10);
//triple lookAt = O;
//currentprojection = perspective(cam, target = lookAt);

//currentprojection=perspective(
//camera=(80,13,30),
//up=(Z),
//target=(1,1,1),
//zoom=0.01,
//angle=0,
//autoadjust=false;

currentprojection=perspective(
camera=(72.0192787736514,-39.5035415631105,24.2748528122902),
up=(-0.0105388868685802,0.00596444269179909,0.0425371172389632),
target=(0.999999999999988,0.999999999999999,1),
zoom=0.00872308215245425,
angle=6.49698942740468,
autoadjust=false);

//real x_centre = 3;
//real y_centre = 3;
real x_centre = 4;
real y_centre = 0;
real rad = 1.5;
triple vertex=(0,0,1);
path p2 = (x_centre,y_centre-rad)..(x_centre+rad,y_centre)..(x_centre,y_centre+rad)..(x_centre-rad,y_centre)..cycle;
path3 p = path3 (p2);
triple f(pair t)
{
    real s=1-t.x; // scale factor for cross section
    return(t.x*vertex+xscale3(s)*yscale3(s)*point(p,t.y*length(p)));
}

surface sph=shift(0,0,.5)*scale3(.5)*unitsphere;

draw(surface(f,(0, 0),(1,1),30,30,Spline),gray+opacity(.3));
draw(sph,paleblue);
draw(p,gray);

triple[] inter;
for(int i=0; i<size(p2);++i)
{
    //draw(point(p,i)--vertex,gray(.3));
    //path3 intray=point(p,i)--vertex;
    //inter.push(intersectionpoints(intray,sph);
}

triple map_to_sphere(triple P1)
{
    real xs=(2*P1.x)/(1+((P1.x)^2)+((P1.y)^2));
    real ys=(2*P1.y)/(1+((P1.x)^2)+((P1.y)^2));
    real zs=(-1+((P1.x)^2)+((P1.y)^2))/(1+((P1.x)^2)+((P1.y)^2));
    return(xs,ys,zs);
}
triple tst=(5.5,0,0);
dot((map_to_sphere(tst)),2bp+green);

xaxis3("$x$",Arrow3);
yaxis3("$y$",Arrow3);
//zaxis3("$z$",Arrow3);

And the output: enter image description here Explanation: Since I failed to use the intersectionpoints() function (getting a persistent syntax error), I resorted to the following tactics: Create a circle in the plane and draw the conic surface. For 4 of the nodes of the unitcircle, find their projections on the unitsphere. Construct a circle based on these 4 new points. Unfortunately, I got stuck, as can be seen in the rendered image. The problem is that the green point should be on the top of the circle on the surface of the sphere, and right now it's floating in space. I guess I implemented the math in a wrong way (but when I try to solve it explicitly I only get a bigger mess). Since I'm a novice I used code from these two sources:

http://asy.gmaths.net/forum/geometrie-espace-f9/intersection-sphere-cone-oblique-t48.html?sid=64f141c0adc9c39e5cb6abe9c44e7de6

http://sourceforge.net/p/asymptote/discussion/409349/thread/8d7f6aa2/

For the record I'm using Win 7 64-bit. I am also aware of the ghostscript bug if that's somehow important (that is, it's fixed).

Side note: I would really appreciate an explained solution.

Update

I added the lines O.G. kindly suggested, and this time Asymptote starts, but never finishes the process. I figured my previous code sample was too long so this time I created a MWE (or a minimum bugged example in this case). Here it is:

settings.render = 4;
settings.prc = false;

import graph3;
size (6cm, 0);

surface sph=shift(0,0,.5)*scale3(.5)*unitsphere;
path3 line=(0,0,1)--(1,1,0);
triple[] tmp;
//tmp[0]=(1,1,1); // Just to make sure I understand how arrays work :)
tmp=intersectionpoints(line,sph,1e-6); // Asymptote freezes because of this line
draw(sph,lightblue+opacity(.8));
draw(line);
dot(tmp[0],1bp+deepgreen);

What happens is that the Asymptote process (asy.exe *32 in my case) consumes only around 30Mb of memory, but seems to freeze (something like an infinite loop, except here is no loop).

So, I'm wondering is this an internal Asymptote bug or am I still incorrectly using the intersectionpoints() function?

P.S. Any solutions to the original problem are welcome (whether or not they use the mentioned problematic function).

Update 2

When the fuzz value is changed from 1e-6 to 1e-4 the code works properly. For the record I run a 64bit Win7 OS on a Core i7 machine with 16GB of RAM. It was reported by O.G. that the code works correctly with the fuzz of 1e-6 on a 64bit Linux OS with an Core i5.

  • Haha i did the same with pgfplots. In fact i started to write a complex number library for this but dumped it after TikZ switched to Lua :) – percusse Apr 12 '15 at 19:36
  • Could you provide the entire asymptote code ? – O.G. Apr 12 '15 at 20:22
  • This is all the code needed to generate the image shown. I've included the links to where I took the code. – blaze Apr 12 '15 at 20:52
  • x_centre, y_centre, etc are not defined. The code cannot be compiled. – O.G. Apr 12 '15 at 21:02
  • I think it is possible to calculate coordinates in pure math... – Symbol 1 Apr 13 '15 at 2:06
1

The intersection is not difficult to compute. Just write the parametrisation of the line, put inside the condition "belong to the sphere". You find two solutions (one is (0,0,1)) and the other one is what you want. No math here :) I prefer beer (instead of reputation).

I modify you code.

settings.render = 0;
settings.prc = false;

import graph3;
import solids;
size (6cm, 0);
//currentprojection=perspective(10*2,0.5*2,6*2);
currentprojection=perspective(10,0.5,6);
dotfactor=2.5;
//grid3(gridroutine=XYXgrid(),Step=1,step=.5);

pen def_front=cmyk(black)+linewidth(.5);
pen def_front_red=cmyk(red)+linewidth(.3);
pen def_back=gray(.5)+dotted+linewidth(.3)+linetype(new real[] {4,2});
pen def_back_solid=gray(.5)+linewidth(.3);
pen def_back_red=cmyk(red)+dotted+linewidth(.3)+linetype(new real[] {4,2});

real x_centre = 1.5;
real y_centre = 1.5;
real rad = 1;
triple vertex=(0,0,1);
path p2 = (x_centre,y_centre-rad)..(x_centre+rad,y_centre)..(x_centre,y_centre+rad)..(x_centre-rad,y_centre)..cycle;
path3 p = path3(p2);

surface sph=shift(0,0,.5)*scale3(.5)*unitsphere;
revolution b=sphere((0,0,0.5),0.5); // Used because of silhouette()
//path3 line=(0,0,1)--(1,1,0);
//triple[] tmp;
//tmp[0]=(1,1,1); // Just to make sure I understand how arrays work :)
//tmp=intersectionpoints(line,sph,1e-4); // Asymptote freezes on this line

draw(b,1,frontpen=def_front,backpen=def_back,longitudinalpen=nullpen); // Draw equator
draw(p,def_front_red); // Draw the circle in the plane

triple[] inter;
int cntr = 0;
for (int i = 0; i < size(p2); ++i)
  {
    path3 intray=point(p,i)--vertex;
    triple pM=point(p,i);
    triple pint=(0,0,1)+scale3(1/(1+(pM.x)^2+(pM.y)^2))*(pM.x,pM.y,-1);
    inter.push(pint);
    draw(point(p,i)--inter[cntr],def_back_solid);
    draw(inter[cntr]--vertex,def_back);
    cntr = cntr + 1;
  }
draw(inter[0]..inter[1]..inter[2]..inter[3]..cycle,def_front_red);
draw(b.silhouette(),def_front);
draw(inter[1]--point(p,1),def_back_solid); // Correct the incorrect overlap of the red circle and the lines
draw(inter[2]--point(p,2),def_back_solid);
draw(inter[3]--point(p,3),def_back_solid);
draw((3,-1.5,0)--(3,3,0)--(-2.5,3,0)--(-2.5,-1.5,0)--cycle,def_back_solid);

cntr = 0;
for (int i = 0; i < size(p2); ++i)
  {
    dot(point(p,i),cmyk(red));
    dot(inter[cntr],cmyk(red));
    cntr = cntr + 1;
  }
  • Accepted because of beer :) – blaze Apr 15 '15 at 21:36
  • Thanks ! As far as mathematics/algorithms are concerned, there are a lot of informations about sphere/ray or line, sphere/plane, etc... intersections. – O.G. Apr 16 '15 at 7:27
9

Note: I edited the code a couple of times (mostly small cosmetic changes). So, here is the code I finally used (pardon all the commented code, part of it is dedicated to creating raster output). I didn't comment out my code since I was the asker :) but if there is need for some explanation post it as a comment and I'll be glad to answer.

settings.render = 0;
settings.prc = false;

import graph3;
import solids;
size (6cm, 0);
//currentprojection=perspective(10*2,0.5*2,6*2);
currentprojection=perspective(10,0.5,6);
dotfactor=2.5;
//grid3(gridroutine=XYXgrid(),Step=1,step=.5);

pen def_front=cmyk(black)+linewidth(.5);
pen def_front_red=cmyk(red)+linewidth(.3);
pen def_back=gray(.5)+dotted+linewidth(.3)+linetype(new real[] {4,2});
pen def_back_solid=gray(.5)+linewidth(.3);
pen def_back_red=cmyk(red)+dotted+linewidth(.3)+linetype(new real[] {4,2});

real x_centre = 1.5;
real y_centre = 1.5;
real rad = 1;
triple vertex=(0,0,1);
path p2 = (x_centre,y_centre-rad)..(x_centre+rad,y_centre)..(x_centre,y_centre+rad)..(x_centre-rad,y_centre)..cycle;
path3 p = path3(p2);

surface sph=shift(0,0,.5)*scale3(.5)*unitsphere;
revolution b=sphere((0,0,0.5),0.5); // Used because of silhouette()
//path3 line=(0,0,1)--(1,1,0);
//triple[] tmp;
//tmp[0]=(1,1,1); // Just to make sure I understand how arrays work :)
//tmp=intersectionpoints(line,sph,1e-4); // Asymptote freezes on this line

draw(b,1,frontpen=def_front,backpen=def_back,longitudinalpen=nullpen); // Draw equator
draw(p,def_front_red); // Draw the circle in the plane

triple[] inter;
int cntr = 0;
for (int i = 0; i < size(p2); ++i)
{
    path3 intray=point(p,i)--vertex;
    triple[] tmp=intersectionpoints(intray,sph,1e-4);
    for (int j=0; j<tmp.length; ++j)
    {
        inter.push(tmp[j]);
    }
    draw(point(p,i)--inter[cntr],def_back_solid);
    draw(inter[cntr]--vertex,def_back);
    cntr = cntr + 2;
}

draw(inter[0]..inter[2]..inter[4]..inter[6]..cycle,def_front_red);
//draw(sph,surfacepen=material(diffusepen=0.5*blue + opacity(0.5), emissivepen=0.5*white));
//draw(line);
//dot(inter[0],2bp+deepgreen);
//dot(inter[2],2bp+deepgreen);
//dot(inter[4],2bp+deepgreen);
//dot(inter[6],2bp+deepgreen);
draw(b.silhouette(),def_front);
draw(inter[2]--point(p,1),def_back_solid); // Correct the incorrect overlap of the red circle and the lines
draw(inter[4]--point(p,2),def_back_solid);
draw(inter[6]--point(p,3),def_back_solid);
draw((3,-1.5,0)--(3,3,0)--(-2.5,3,0)--(-2.5,-1.5,0)--cycle,def_back_solid);

cntr = 0;
for (int i = 0; i < size(p2); ++i)
{
    dot(point(p,i),cmyk(red));
    dot(inter[cntr],cmyk(red));
    cntr = cntr + 2;
}

Which produces the output: enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.