# Stereographic projection of a circle to a plane in Asymptote

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,

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,

//real x_centre = 3;
//real y_centre = 3;
real x_centre = 4;
real y_centre = 0;
triple vertex=(0,0,1);
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: 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

``````

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 :) Commented Apr 12, 2015 at 19:36
• Could you provide the entire asymptote code ?
– O.G.
Commented Apr 12, 2015 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. Commented Apr 12, 2015 at 20:52
• `x_centre`, `y_centre`, etc are not defined. The code cannot be compiled.
– O.G.
Commented Apr 12, 2015 at 21:02
• I think it is possible to calculate coordinates in pure math... Commented Apr 13, 2015 at 2:06

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;
triple vertex=(0,0,1);
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 :) Commented Apr 15, 2015 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.
Commented Apr 16, 2015 at 7:27

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;
triple vertex=(0,0,1);
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: