# Fill enclosed region on sphere

Given the longitude and latitude coordinates of multiple points on a sphere, I can connect them using circular arc segments by using the `Arc` command in Asymptote:

``````path3 arc(triple c,real r,real theta1,real phi1,real theta2,real phi2,triple normal=O);
``````

As an example, given `(-10,0)`, `(0,10)`, `(10,0)` and `(-10,0)` in degrees, I can use the following code to draw four black arcs on the sphere to make it look like I've connected the four points from start to finish (and cycled).

``````settings.render = 0;

import solids;
import three;
size(6cm);

real zenith  = pi/12.0;
real azimuth = pi/12.0;
currentprojection = perspective(cos(azimuth)*cos(zenith),
sin(azimuth)*cos(zenith),
sin(zenith));
real r = 1;
real ar = 1.2;

path3 myarc = Arc(c=O,normal=X, v1=-Z*r, v2=Z*r, n=12);
surface sphere = surface(myarc, angle1=0, angle2=360, c=O, axis=Z, n=12);

draw((0,0,  0)--(ar,0,0),   red+linewidth(1pt));
draw((0,0,  0)--(0,ar,0), green+linewidth(1pt));
draw((0,0,-ar)--(0,0,ar),  blue+linewidth(1pt));

label("\$S\$",(0,0,-ar),S);
draw(sphere, surfacepen=material(white+opacity(0.8),ambientpen=white));
label("\$N\$",(0,0,ar),N);

pair[] region = new pair[] {(-10,0),(0,10),(10,0),(0,-10),(-10,0)};
for(int i=1; i<region.length; ++i){
pair one, two;
if(region[i-1].y < region[i].y){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else if(region[i-1].y > region[i].y){
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}else if(region[i-1].x > region[i].x){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else{
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}
path3 temp = arc(O,r,one.y,one.x,two.y,two.x);
draw(temp,black+linewidth(.5pt));
}
``````

Notes: the `if - else if - else if - else` block is used to draw the shorter one between the two possible arcs between each pair of point.

The code above produces this image without the 3D viewer being enabled.

Now, I'd like to fill that region on the sphere with any color (for example, red with 50% transparency). How can I achieve this using Asymptote?

• I'm aware that the sphere is completely covering the red axis, but I would like to think about it on my own for now. :) – Frenzy Li Jul 25 '16 at 11:13
• I have solved the issue in my first comment by drawing two hemispheres. I draw a hemisphere from -90 to 90 degrees, then I draw the red axis, after which I draw the second hemisphere from 90 to 270 degrees. – Frenzy Li Jul 25 '16 at 11:49
• You want to create a red surface covering the region you want to "fill." Since you're using `settings.render=0`, you need to draw the red surface after you've drawn the sphere. If you were using rasterized rendering, you'd instead make the red surface part of a sphere with a slightly larger radius (e.g., 1.001 instead of 1) so that the program knows which surface is closer to the viewer. – Charles Staats Jul 26 '16 at 17:05
• Additional note: since you're using `arc` and generating arcs of less than 90 degrees, I'm guessing each arc in the black patch consists of only one segment. You should be able to string these together into a cyclic path of length four using the `&` connector: `path outline = arc1 & arc2 & arc3 & arc4 & cycle;`. At this point, `surface(patch(outline))` might give you a reasonable surface to fill the patch, at least with `render=0`. – Charles Staats Jul 26 '16 at 17:08
• @CharlesStaats Is it that: by drawing a cyclic path and patching it won't necessarily produce the spherical region, but an approximation to it that will land inside the sphere? – Frenzy Li Jul 26 '16 at 17:12

Sorry I made a confusion in the first answer about the region (parallel circles and great circles), you can see the code at the end (the approximation is ok for small value of angle <30).

About the second solution, i.e. construct the region with a single patch, computing the derivative is easy (from the arc computation of Asymptote) but the obtained surface is far from the desired result (try to draw the `oo3` surface in the code). I think that it is possible by decomposing the region into 4 or 8 parts but in any case it is crucial to compute the internal control points of the Bézier patch. It needs some knowledge about Bézier patch and surface approximation.

So I tried to complete the first solution: a parametric definition of the region. Since I do not want to define myself, I used the `Arc` definition of Asymptote. The strange parametrization (see the definition of f) is needed to avoid some numerical artefacts in the interpolation and you can recognize the Tchebychev points. It seems to be ok, for angle less or equal than 80. Please consider the code

``````import solids;
import three;
size(6cm);

real zenith  = pi/12.0;
real azimuth = pi/12.0;
currentprojection = perspective(cos(azimuth)*cos(zenith),
sin(azimuth)*cos(zenith),
sin(zenith));
real r = 1;
real ar = 1.2;

path3 myarc = Arc(c=O,normal=X, v1=-Z*r, v2=Z*r, n=12);
surface sphere = surface(myarc, angle1=0, angle2=360, c=O, axis=Z, n=12);

draw((0,0,  0)--(ar,0,0),   red+linewidth(1pt));
draw((0,0,  0)--(0,ar,0), green+linewidth(1pt));
draw((0,0,-ar)--(0,0,ar),  blue+linewidth(1pt));

label("\$S\$",(0,0,-ar),S);
draw(sphere, surfacepen=material(white+opacity(0.8),ambientpen=white));
label("\$N\$",(0,0,ar),N);
real angle=80;
path3[] aar;
pair[] region = new pair[] {(-angle,0),(0,angle),(angle,0),(0,-angle),(-angle,0)};
for(int i=1; i<region.length; ++i){
pair one, two;
if(region[i-1].y < region[i].y){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else if(region[i-1].y > region[i].y){
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}else if(region[i-1].x > region[i].x){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else{
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}
path3 temp = Arc(O,r,one.y,one.x,two.y,two.x,32);
draw(temp,black+linewidth(.5pt));
aar.push(temp);
}
dot((cos(angle/180*pi),sin(angle/180*pi),0));
triple AA=(cos(angle/180*pi),sin(angle/180*pi),0);
dot((cos(angle/180*pi),sin(-angle/180*pi),0),blue);
triple BB=(cos(angle/180*pi),sin(-angle/180*pi),0);
dot((cos(angle/180*pi),0,-sin(angle/180*pi)),red);
triple CC=(cos(angle/180*pi),0,-sin(angle/180*pi));
dot((cos(angle/180*pi),0,sin(angle/180*pi)),green);
triple DD=(cos(angle/180*pi),0,sin(angle/180*pi));

real angr=angle/180*pi;
triple f(pair t)
{
triple x1=point(aar[1],(1+cos((2*t.x+1)*pi))/2*length(aar[1]));
triple x2=point(aar[3],(1+cos((2*t.x+1)*pi))/2*length(aar[3]));
path3 temp_arc=Arc(O,x1,x2,32);
return(point(temp_arc, (1+cos((2*t.y+1)*pi))/2*length(temp_arc)));
}
surface s=surface(f,(0,0),(1,1),32,32,Spline);
draw(shift(0.005*ar,0,0)*s,2bp+blue);

patch oo3=patch(AA{-dir(aar[1])}..{-dir(aar[1],0)}DD{dir(aar[0],0)}..{dir(aar[0])}BB{dir(aar[3],0)}..{dir(aar[3])}CC{-dir(aar[2])}..{-dir(aar[2],0)}cycle);
``````

and the result (angle=80) with OpenGL renderer.

Previous anwser : With a few computations, it is possible to describe the region through a parametric surface and then to draw it.

In the spirit of Charles Staats indication, you can have an approximation of the surface using the `patch` constructor. By specifying the tangent at the vertices the approximation is ok for low angle. The best should be to compute the last four control points of the surface (see the definition of the sphere in `three_surface.asy`). To avoid some artefacts the surface is shifted and it works with `render=0` and the OpenGL renderer.

``````import solids;
import three;
size(6cm);

real zenith  = pi/12.0;
real azimuth = pi/12.0;
currentprojection = perspective(cos(azimuth)*cos(zenith),
sin(azimuth)*cos(zenith),
sin(zenith));
real r = 1;
real ar = 1.2;

path3 myarc = Arc(c=O,normal=X, v1=-Z*r, v2=Z*r, n=12);
surface sphere = surface(myarc, angle1=0, angle2=360, c=O, axis=Z, n=12);

draw((0,0,  0)--(ar,0,0),   red+linewidth(1pt));
draw((0,0,  0)--(0,ar,0), green+linewidth(1pt));
draw((0,0,-ar)--(0,0,ar),  blue+linewidth(1pt));

label("\$S\$",(0,0,-ar),S);
draw(sphere, surfacepen=material(white+opacity(0.8),ambientpen=white));
label("\$N\$",(0,0,ar),N);

pair[] region = new pair[] {(-10,0),(0,10),(10,0),(0,-10),(-10,0)};
triple[] mypoints;
for(int i=1; i<region.length; ++i){
pair one, two;
if(region[i-1].y < region[i].y){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else if(region[i-1].y > region[i].y){
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}else if(region[i-1].x > region[i].x){
one = (region[i]  .x, 90 - region[i]  .y);
two = (region[i-1].x, 90 - region[i-1].y);
}else{
one = (region[i-1].x, 90 - region[i-1].y);
two = (region[i]  .x, 90 - region[i]  .y);
}
path3 temp = arc(O,r,one.y,one.x,two.y,two.x);
draw(temp,black+linewidth(.5pt));
}
real angle=10;
dot((cos(angle/180*pi),sin(angle/180*pi),0));
triple AA=(cos(angle/180*pi),sin(angle/180*pi),0);
dot((cos(angle/180*pi),sin(-angle/180*pi),0),blue);
triple BB=(cos(angle/180*pi),sin(-angle/180*pi),0);
dot((cos(angle/180*pi),0,-sin(angle/180*pi)),red);
triple CC=(cos(angle/180*pi),0,-sin(angle/180*pi));
dot((cos(angle/180*pi),0,sin(angle/180*pi)),green);
triple DD=(cos(angle/180*pi),0,sin(angle/180*pi));

patch oo=patch(AA{cross(AA,AA-CC)}..{cross(DD,DD-BB)}DD{cross(DD,DD-AA)}.. {cross(BB,BB-CC)}BB{cross(BB,BB-DD)}..{cross(CC,CC-AA)}CC{cross(CC,AA-DD)}..{cross(AA,AA-DD)}cycle);
draw(shift(.005*ar,0,0)*surface(oo),blue);
``````

And the result

For angle=45, the difference is apparent (with the OpenGL renderer)...

• This actually works pretty well for small angles, but as you say, angle=45 or angle=90, the patch doesn't have the correct tangents. I think if the tangent leaving `AA` would be in the same plane with the great circle determined by AA and BB, the difference will be gone. – Frenzy Li Aug 4 '16 at 13:25
• Not sure that the tangents are wrong. A Bezier patch has 16 control points, here 12 are defined and the last 4 are determined by Asymptote. You can observe a small difference between the two definitions `real a=4/3*(sqrt(2)-1);patch oo1=patch(X{Y}..{-X}Y{Z}..{-Y}Z..Z{X}..{-Z}cycle, new triple[] {(1,a,a),(a,1,a),(a^2,a,1), (a,a^2,1)}); patch oo2=patch(X{Y}..{-X}Y{Z}..{-Y}Z..Z{X}..{-Z}cycle);` `oo1` is `octant` defined in `three_surface.asy` to construct the sphere. – O.G. Aug 4 '16 at 15:28
• @Frenzy DT : You're right, the tangents should be modified. – O.G. Aug 4 '16 at 16:54