.PS
# graysurf.m4
threeD_init
NeedDpicTools

Overall: [

Surface:[
#                                  This just draws the facets back to front

  setview(20,30)
  dimen = 1.5
  
  xmax = 2*dimen
  ymax = 1.5*dimen
  zmax = 1*dimen
  
  define(`fn',
  `expe(-(`$1')/tt)*cos((`$2')/xmax*1.75*twopi_)*cos((`$1')/ymax*twopi_)')
  
  nx = 24; dx = xmax/nx
  ny = 30; dy = ymax/ny
  tt = 1.5
  fnmax = fn(0,0,0)
  
  Orig: (0,0)
  linethick_(1.2)
  arrow from Orig to 1.2 <Orig,(Project(xmax,0,0))> ; "$x$" below
  arrow from Orig to 1.2 <Orig,(Project(0,ymax,0))> ; "$y$" ljust
  arrow from Orig to 1.2 <Orig,(Project(0,0,zmax))> ; "$z$" rjust
  
  linethick_(0)
  for i = 0 to nx-1 do {
    x = xmax*i/nx
    for j = 0 to ny-1 do {
      y = ymax*j/ny
      fnxy = fn(x,y)
      shade((fnmax+fnxy)/fnmax/2,
        line from (Project(x,y,fnxy)) to (Project(x,(y+dy),`fn(x,(y+dy))'))\
          then to (Project((x+dx),(y+dy),`fn((x+dx),(y+dy))'))\
          then to (Project((x+dx),y,`fn((x+dx),y)'))\
          then to (Project(x,y,fnxy))
        ) } }
] #with .w at last [].e+(-0.2,0)

Torus: [
#                                  Calculate all the facet centres but draw
#                                  only a subset of them after sorting

  viewazimuth = 20                 # Set view angles in degrees
  viewelevation = 30

  setview(viewazimuth,viewelevation,0)

  sfactor = 0.7
  tradius = 2*sfactor
  sradius = 0.75*sfactor

define(`torus',`rot3Dz((`$1')*dtor_,
  sum3D(tradius,0,0,rot3Dy((`$2')*dtor_,sradius,0,0)))')

define(`tvisible',
 `dot3D(View3D, rot3Dz((`$1')*dtor_,rot3Dy((`$2')*dtor_,1,0,0)))')

Orig: (0,0)
X: arrow from Orig to (Project(tradius+sradius*4,0,0)); "x" rjust
Y: arrow from Orig to (Project(0,tradius+sradius*2,0)); "y" ljust

  dt = 10                          # major angle degrees per facet
  ds = 10                          # minor angle degrees per facet

                                   # create the arrays and sort
  n = 0 
  for t = dt to 360+dt/2 by dt do {
    for s=ds/2 to 360.1 by ds do { 
      if tvisible(t,s) >= 0 then {
      n +=1
      t[n] = t
      s[n] = s
      d[n] = dot3D(torus(t,s),view3D1,view3D2,view3D3) # view distance
      ix[n] = n
      } } }
#prval(n)
  dpquicksort(d,1,n,ix)

                                   # draw the facets
thinlines_
  for i = 1 to n do {
    tc = t[ix[i]]; sc = s[ix[i]]
    SE: (Project(torus(tc+dt/2,sc-ds/2)))
    SW: (Project(torus(tc-dt/2,sc-ds/2)))
    NW: (Project(torus(tc-dt/2,sc+ds/2)))
    NE: (Project(torus(tc+dt/2,sc+ds/2)))
    f = ((dcosine3D(3,torus(tc,sc))/sradius+1)/2)^2
    g = min(f,0.9)
    setrgb(g,g,g)
    line fill f ifpdf( invis ) \
      from SE to SW then to NW then to NE then to SE
#   line invis shaded rgbstring(f,f^3,f^3) \
#     f from SE to SW then to NW then to NE then to SE
    resetrgb
  }

thicklines_
  arrow from Orig to (Project(0,0,tradius*1.4)); "z" above
  dashline(from Orig to (Project(torus(0,0))),,,,G)
  arrow to X.end
  dashline(from Orig to (Project(torus(90,0))),,,,G)
  arrow to Y.end

] with .sw at last [].se+(0.2,0)

Hat: [
#                                  Sort as for the torus but draw edges too

  azimuth = 25                     # View angles in degrees
  elevation = 20
  setview(azimuth,elevation)
                                   # the 3D object: `hat(t,y)'
  define(`hat',`cosd($1)*($2),sind($1)*($2),(($2)^2-1)^2') 
  define(`gradient',`(-4*($2)*(($2)^2-1)*cosd(`$1')),
                     (-4*($2)*(($2)^2-1)*sind(`$1')),1')
  define edge {$2=dot3D(View3D,gradient($1,y))}
  define(`hatvis',
   `(($2 > 0.6) || (dot3D(View3D,gradient($1,$2)) > 0))')

                                   # grid parameters
  dang = 09
  maxy = 1.25
  dy = 0.05

Orig: (0,0)
X: arrow from Orig to (Project(maxy*1.2,0,0)); "x" rjust
Y: arrow from Orig to (Project(0,maxy*1.2,0)); "y" ljust
Xv: (Project(1,0,0))
Yv: (Project(0,1,0))
Zv: (Project(0,0,1))
                                   # create the arrays and sort
  n = 0 
  for t = dang/2 to 360.1 by dang do {
    for y=dy/2 to maxy by dy do {
      if hatvis(t,y) then {
      n +=1
      t[n] = t
      y[n] = y
      d[n] = dot3D(hat(t,y),View3D) #distance toward the front
      ix[n] = n
      } } } 
#prval(n)
  dpquicksort(d,1,n,ix)
                                   # draw the facets
  ellipse wid maxy*2 ht maxy*2*sind(elevation) \
    with .c at (Project(0,0,(maxy^2-1)^2)+(0,linethick bp__/2))
thinlines_
  for i = 1 to n do {
    tc = t[ix[i]]; yc = y[ix[i]]
    SE: (Project(hat(tc+dang/2,yc-dy/2)))
    SW: (Project(hat(tc-dang/2,yc-dy/2)))
    NW: (Project(hat(tc-dang/2,yc+dy/2)))
    NE: (Project(hat(tc+dang/2,yc+dy/2)))
    f = abs(yc^2-1)
    line invis fill f from SE to SW then to NW then to NE then to SE
  }

                                   # Outline where the surface goes invisible
mm = int(maxy/dy)
nr = -1
for i=1 to mm do { y = i/(mm+3)
  findroot(edge,(azimuth+2),(azimuth+180-2),1e-8,t)
  nr +=1
  R[mm+nr]: (Project(hat(t,y)))
  R[mm-1-nr]: (-R[mm+nr].x,R[mm+nr].y)
  }
fitcurve(R,mm*2-1)

                                   # Z axis
thicklines_
Z: arrow from Zv to (Project(0,0,1.2)); "z" ljust
  line dashed from Orig to Zv chop 0 chop dashwid/2
  line dashed from Orig to Xv
  line dashed from Orig to Yv

                                   # Partial rim
ellipsearc(maxy*2,maxy*2*sind(elevation),-pi_*1.4,pi_/3) \
  with .C at (Project(0,0,(maxy^2-1)^2)+(0,linethick bp__/2))
] scaled 1.25 with .s at 1st [].ne+(0,-0.75)
] scaled 0.85

.PE