The equations needed to construct the image above are, in fact, from a detailed analysis of Fuller's projection by Robert W. Gray . It is not simple to develop this projection using the mapcoords functions given elsewhere. Here we take a different approach in which we have developed the necessary Maple functions by translating Gray's C code from the web site listed above. We should also note that it will be very difficult to follow all the deails without referring to the information on Gray's web site.

First we identify the coordinates of the vertices of the triangles

> v[x] := [.420152426708710003, .995009439436241649, .518836730327364437, -.414682225320335218, -.515455959944041808, .355781402532944713, .414682225320335218, .515455959944041808, -.355781402532944713, -.995009439436241649, -.518836730327364437, -.420152426708710003];

> v[y] := [.78145249402782959e-1, -.91347795276427931e-1, .835420380378235850, .655962405434800777, -.381716898287133011, -.843580002466178147, -.655962405434800777, .381716898287133011, .843580002466178147, .91347795276427931e-1, -.835420380378235850, -.78145249402782959e-1];

> v[z] := [.904082550615019298, .40147175877166645e-1, .181331837557262454, .630675807891475371, .767200992517747538, .402234226602925571, -.630675807891475371, -.767200992517747538, -.402234226602925571, -.40147175877166645e-1, -.181331837557262454, -.904082550615019298];

Specify the angles through which each triangle is rotated.

> angles:=[240,300,0,60,180,300,300,0,300,60,60,120,60,0,0,60,0,120,120,300];

The x and y coordinates determined below are shifted so that each triangle is correctly positioned. The shift amounts are given below (except for the LCD triangles that are cut and moved separately).

> xb:=[2.0, 2.0,2.5,3.0, 2.5, 1.5, 1.0, 1.5, 1.5, 2.5, 3.5, 3.5, 4.0, 4.0, 5.0, 0.5, 1.0, 4.0, 4.5, 5.0];

> yb:=[7.0 / (2.0 * sqrt(3.0)),5.0 / (2.0 * sqrt(3.0)),2.0 / sqrt(3.0),5.0 / (2.0 * sqrt(3.0)),4.0 * sqrt(3.0) / 3.0,4.0 * sqrt(3.0) / 3.0,5.0 / (2.0 * sqrt(3.0)),2.0 / sqrt(3.0),1.0 / sqrt(3.0),1.0 / sqrt(3.0),1.0 / sqrt(3.0),2.0 / sqrt(3.0),5.0 / (2.0 * sqrt(3.0)),7.0 / (2.0 * sqrt(3.0)),7.0 / (2.0 * sqrt(3.0)),1.0 / sqrt(3.0),1.0 / (2.0 * sqrt(3.0)),1.0 / (2.0 * sqrt(3.0)),2.0 / sqrt(3.0),5.0 / (2.0 * sqrt(3.0))];

Some preliminary calculations and functions:

> hold := proc(v,i,j,k)
local letter;
eval(convert([seq(letter=(v[letter][i]+v[letter][j]+v[letter][k])/3,letter=[x,y,z])],table));
end;

> magn := (v,i,j,k)->sqrt(add(hold(v,i,j,k)[letter]^2,letter=[x,y,z]));

> centerpoint := proc(v,i,j,k)
local magdiv;
magdiv := (x,y)-> x/y;
map(magdiv,hold(v,i,j,k),magn(v,i,j,k))
end;

Specify which vertices are at the corners of each triangle.

> corners:=[[1, 3, 2], [1, 4, 3], [1, 5, 4], [1, 6, 5], [1, 2, 6], [2, 3, 8], [3, 9, 8], [3, 4, 9], [4, 10, 9], [4, 5, 10], [5, 11, 10], [5, 6, 11], [6, 7, 11], [2, 7, 6], [2, 8, 7], [8, 9, 12], [9, 10, 12], [10, 11, 12], [11, 7, 12], [8, 12, 7]];

Work out locations of all centers:

> allcenters := [seq(centerpoint(v,op(c)),c=corners)]:

Some more general calculations:

> garc := 2.0 * arcsin( sqrt( 5. - sqrt(5.)) / sqrt(10.0) );
gt := garc / 2.0;
gdve := sqrt( 3 + sqrt(5.) ) / sqrt( 5. + sqrt(5.0) );
gel := sqrt(8.0) / sqrt(5.0 + sqrt(5.0));

Function to identify which triangle the current point is in.

> whichtriangle := proc(p)
local ca, letter, h_tri, h_dist, h_dist1, h_dist2, tri, i, lcd, h_lcd, h1;
global allcenters, corners;
h_tri := 0;
h_dist1 := 9999.0;

Which triangle face center is the closest to the given point is the triangle in which the given point is in
for i from 1 to 20 do
ca := allcenters[i];
for letter in [x,y,z] do h1[letter] := ca[letter] - p[letter]; od;
h_dist2 := sqrt(add(h1[letter]^2,letter=[x,y,z]));
if (h_dist2 < h_dist1) then
h_tri := i;
h_dist1 := h_dist2;
fi;
od;

Now we know the triangle.
tri := h_tri;
for i from 1 to 3 do
for letter in [x,y,z] do h1[letter] := p[letter]-v[letter][corners[tri][i]]; od;
h_dist[i] := sqrt(add(h1[letter]^2,letter=[x,y,z]));
od;

if ( (h_dist[1] <= h_dist[2]) and (h_dist[2] <= h_dist[3]) ) then h_lcd := 1 fi;
if ( (h_dist[1] <= h_dist[3]) and (h_dist[3] <= h_dist[2]) ) then h_lcd := 6 fi;
if ( (h_dist[2] <= h_dist[1]) and (h_dist[1] <= h_dist[3]) ) then h_lcd := 2 fi;
if ( (h_dist[2] <= h_dist[3]) and (h_dist[3] <= h_dist[1]) ) then h_lcd := 3 fi;
if ( (h_dist[3] <= h_dist[1]) and (h_dist[1] <= h_dist[2]) ) then h_lcd := 5 fi;
if ( (h_dist[3] <= h_dist[2]) and (h_dist[2] <= h_dist[1]) ) then h_lcd := 4 fi;
RETURN([tri,h_lcd]);
end:

> vc:=[1,1,1,1,1,2,3,3,4,4,5,5,6,2,2,8,9,10,11,8];

Some conversion routines.

> s_to_c := proc(theta, phi)
# Covert spherical polar coordinates to cartesian coordinates
# The angles are given in radians
local temp;
temp := [x = sin(theta) * cos(phi), y = sin(theta) * sin(phi), z = cos(theta)];
RETURN(eval(convert(evalf(temp),table)));
end;

> c_to_s := proc(x,y,z)
# convert cartesian coordinates into spherical polar coordinates.
# The angles are given in radians.
local a, lng, lat;
if (x>0.0 and y>0.0) then a := 0 fi;
if (x<0.0 and y>0.0) then a := Pi fi;
if (x<0.0 and y<0.0) then a := Pi fi;
if (x>0.0 and y<0.0) then a := 2*Pi fi;
lat := arccos(z);
if (x=0.0 and y>0.0)then lng := Pi/2 fi;
if (x=0.0 and y<0.0)then lng := 3*Pi/2 fi;
if (x>0.0 and y=0.0) then lng := 0 fi;
if (x<0.0 and y=0.0) then lng := Pi fi;
if (x <> 0.0 and y <> 0.0) then lng := arctan(y/x) + a fi;
RETURN(evalf([lat,lng]));
end:

Conversion to spherical coords. It appears to me that this routine (translated from Gray) uses the angle "up" from the equatorial plane as one coordinate.

> `convert/spherical` := proc(P::list(numeric))
local h_theta, h_phi, theta, phi;
h_theta := 90.0 - P[2];
h_phi := P[1];
if (P[1] < 0.0) then h_phi := P[1] + 360.0 fi;
theta := h_theta * Pi/180;
phi := h_phi * Pi/180;
RETURN(evalf([theta,phi]));
end:

Some angular rotations (these will be replaced by Maple's own functions when all of this works).

> rotate:=proc(angle, p)
local a;
a := angle * Pi/180;
RETURN(eval(convert(evalf([x = p[x] * cos(a) -p[y] * sin(a), y = p[x] * sin(a) +p[y] * cos(a)]),table)));
end;

> r2 :=proc(axis, alpha, h)
local a, b, c, temp;
a := h[x];
b := h[y];
c := h[z];
if (axis = 1) then
temp:=[x=a,y= b * cos(alpha) + c * sin(alpha), z=c * cos(alpha) - b * sin(alpha)]
elif (axis = 2) then
temp:=[x=a * cos(alpha) - c * sin(alpha), y=b,z=a * sin(alpha) + c * cos(alpha)]
elif (axis = 3) then
temp:=[x=a * cos(alpha) + b * sin(alpha), y=b * cos(alpha) - a * sin(alpha), z=c]
fi;
RETURN(eval(convert(evalf(temp),table)));
end:

The main routine for the point calculation.

> dymax_point := proc(tri, lcd, pp, country)
local h0, h1, v1, temp, hlat,p,letter,
hlng, axis,gz, gs, gxp,gyp,ga1p,ga2p,ga3p,ga1,ga2,ga3,gx, gy, px, py, case;
global vc, allcenters,v;
p := eval(pp);
case := tri;
h0 := p;
v1 := vc[tri];
h1 := table([seq(letter=v[letter][v1],letter=[x,y,z])]);
temp:=c_to_s(allcenters[tri][x],allcenters[tri][y],allcenters[tri][z]);
hlat := temp[1]; hlng := temp[2];
axis := 3;
h0:=r2(axis,hlng,h0);
h1:=r2(axis,hlng,h1);
axis := 2;
h0:=r2(axis,hlat,h0);
h1:=r2(axis,hlat,h1);
temp:=c_to_s(h1[x],h1[y],h1[z]);
hlng :=evalf(temp[2] - Pi/2);
hlat:=temp[1];
axis := 3;
h0:=r2(axis,hlng,h0);
gz := sqrt(1 - h0[x] * h0[x] - h0[y] * h0[y]);
gs := sqrt( 5. + 2 * sqrt(5.) ) / ( gz * sqrt(15.) );
gxp := h0[x] * gs ;
gyp := h0[y] * gs ;
ga1p := 2.0 * gyp / sqrt(3.0) + (gel / 3.0) ;
ga2p := gxp - (gyp / sqrt(3.0)) + (gel / 3.0) ;
ga3p := (gel / 3.0) - gxp - (gyp / sqrt(3.0));
ga1 := gt + arctan( (ga1p - 0.5 * gel) / gdve);
ga2 := gt + arctan( (ga2p - 0.5 * gel) / gdve);
ga3 := gt + arctan( (ga3p - 0.5 * gel) / gdve);
gx := 0.5 * (ga2 - ga3) ;
gy := (1.0 / (2.0 * sqrt(3.0)) ) * (2 * ga1 - ga2 - ga3);
p[x] := gx / garc;
p[y] := gy / garc;
if (case <> 16 and case <> 9) or (case = 9 and lcd > 2) or (case = 16 and lcd < 4) then
p := rotate(angles[case],p);
px := p[x] + xb[case];
py := p[y] + yb[case];
elif case = 9 then
p := rotate(0,p);
px := p[x] + 2.0;
py := p[y] + 1.0 / (2.0 * sqrt(3.0));
elif case = 16 then
p := rotate(0,p);
px := p[x] + 5.5;
py := p[y] + 2.0 / sqrt(3.0);
fi;
RETURN(evalf([px,py]));
end:

>

> `convert/dymax` := proc(P)
local temp, temp2,tri,hcld;
# Convert the given (long.,lat.) coordinate into spherical
# polar coordinates (r, theta, phi) with radius=1. Angles are converted to radians.
temp := `convert/spherical`(P);
# convert the spherical polar coordinates into cartesian (x, y, z) coordinates.
temp:=s_to_c(temp[1],temp[2]);
# determine which of the 20 spherical icosahedron triangles the given point is in and the LCD triangle.
temp2:=whichtriangle(temp);
tri := temp2[1];
hcld:=temp2[2];
# Determine the corresponding Fuller map plane (x, y) point.
temp := dymax_point(tri, hcld, temp);
RETURN(eval(temp));
end:

Here is the world map.

>

> libname:=`c:/earth`,libname;

> with(Maps);

> dmap:=plot([seq(map(convert,shore[k],dymax),k=0..1866)],color=black):

> plots[display]({dmap},scaling=constrained);

To make a grid we proceed as follows:

> xvalues:=[seq(k*30,k=-6..6)]; yvalues:=[-90,seq(k*30,k=-2..2),90];

> grid:=graticule(xvalues,yvalues):

> dymaxgrid:= proc(areamap::PLOT)
local newmap, allcurves, numcurves, datalist, numpoints, slist, i, ii, k, newpoints, delta_grid;
newmap := areamap;
if nargs = 1 then delta_grid := 50 else delta_grid := args[2] fi;
allcurves := select(has, [op(areamap)], CURVES);
numcurves := nops(allcurves);
for i to numcurves do for datalist in [op(allcurves[i])] do
if type(datalist, list) then
newpoints := map(convert,datalist,dymax);
newmap := subs(datalist = newpoints, newmap)
fi
od
od;
RETURN(newmap)
end:

> dymaxgrid(grid):

To overlay the triangular coordinate grid we make a list of the latitude and longitude of the points on R.W. Gray's list.ys

> gridpts := [[10.53620,64.7],[-5.24539,2.300882],[58.15771,10.447378],[122.3,39.1],[-143.47849,50.103201],[-67.13233,+23.717925],[36.52151,-50.1032],[112.86767,-23.71793],[174.75461,-2.300882],[-121.84229,-10.44735],[-57.7,-39.1],[-169.4638,-64.7]];

The conversion of these points to the Fuller system gives

> dpts:=map(convert,gridpts,dymax);

The triangular grid is created by joining the necessary points with stright lines. Note how some of the triangles must be cut and moved elewhere.

The height of a triangle is

> yval:=sqrt(3.0)/2;

Here is a plot of the triangular grid.

> g3:=plot({[[0.5,0],[1.5,0],[1.5,yval*2/3],[2,yval],[2,yval/3],[2.5,0]],[[0.75,yval/2],[0,yval],[5.5,yval],[5.5,yval*2]],[[0.5,yval*2],[5.5,yval*2]],[[1,yval*3],[3,yval*3]],
[[1,yval*3],[2.5,0],[3,yval]],[[2,yval*3],[0.5,0]],[[3,yval*3],[2,yval]],[[2,yval*3],[3.5,0],[4,yval],[4.5,0],[3.5,0]],[[0.5,yval*2],[1.5,0]],
seq([[k+0.5,yval*2],[k,yval]],k=3..5),seq([[k+0.5,yval*2],[k+1,yval]],k=3..4),
[[3.5,yval*2],[4,yval*3],[4.5,yval*2]], [[4.5,yval*2],[5,yval*3],[5.5,yval*2]]},
color=blue):

> dgrid:=plots[display]({g3},scaling=constrained): dgrid;

To create the map at the start of these notes we simply combine the grid above with the map.

> plots[display]({dgrid,dmap},scaling=constrained,axes=none);

>