% Torus.m4
.PS
gen_init
threeD_init
copy "HOMELIB_/dpictools.pic"

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

  setview(viewazimuth,viewelevation,0)

  tradius = 2
  sradius = 0.5
# sradius = 0.75

define(`tsection',`rot3Dz(`$1',rot3Dy(`$2',sradius,0,0))')
define(`torus',`rot3Dz(`$1',sum3D(tradius,0,0,rot3Dy(`$2',sradius,0,0)))')
#efine(`tvisible',`dot3D(tsection(`$1',`$2'),View3D)')
define tvisible {$3 = 0
 [ d3D = dot3D(tsection(`$1',`$2'),View3D)
   $3 := d3D
#  X:  project(PtoBase3D(0,project(torus(`$1',`$2'))))
   ] }

define xvis { tvisible($1,yang,$2) }
define yvis { tvisible(zang,$1,$2) }

define(`pr3',`print sprintf("%g, %g, %g",`$1',`$2',`$3')')
#pr3(0,project(0,1,0))
#pr3(0,project(PtoBase3D(0,project(0,1,0))))

  arrow from (0,0) to project(tradius,0,0); "x" rjust
  arrow from (0,0) to project(0,tradius,0); "y" ljust
  arrow from (0,0) to project(0,0,tradius); "z" above

  nz = 36
  ny = 44
   linethick = 2.5; ny = 128

define innerpart {
  tvisible($1,yang,vv)
  if v*vv < 0 then {
    findroot(xvis,zz,$1, 1e-6, z)
    $1 = z
    if v >= 0 then { n += 1; P[n]: (project(torus(z,yang)))
      if n > 2 then {fitcurve(P,n,)}
      n = -1 }
    v = -v }
  zz = $1
  if v >= 0 then {
    n +=1
    P[n]: (project(torus($1,yang)))
#   if n==0 then { sprintf("%g",y) at P[0] }
    }
  }

  zz0 = -pi_+viewazimuth*dtor_
  for y=0 to ny/2 do { yang = (90-y/ny*360)*dtor_
#   uu = dot3D(View3D,tsection(zz0,yang))
    tvisible(0,yang,uu)
    if uu > 0 then { u = uu/sradius
      setrgb(u,u,u) # print u
      }
    n = -1
    zz = zz0
    tvisible(zz,yang,v)
    loop( zang = zz0-pi_/2-pi_/nz, zang < zz0+pi_/2+pi_/nz, zang+=twopi_/nz,
      innerpart(zang))
    if n > 2 then { fitcurve(P,n,); n = -1 } 
    loop( zang = zz0+pi_/2+pi_/nz, zang < zz0+3*pi_/2+pi_/nz, zang+=twopi_/nz,
      innerpart(zang))
    if n > 2 then { fitcurve(P,n,); n = -1 } 
    }

if 1 < 0 then {
  nq = -1
  ns = -1
  for z=0 to nz do { zang = z/nz*twopi_
    n = -1
    yy = 0
    tvisible(zang,yy,v)
    loop( yang = 0, yang < (1+2/ny)*twopi_, yang += twopi_/ny,
      tvisible(zang,yang,vv)
      if v*vv < 0 then {
        findroot(yvis,yy,yang, 1e-6, y)
        yang = y
        if v >= 0 then { n += 1; P[n]: (project(torus(zang,yang)))
          ns +=1; S[ns]: P[n]
          if n > 2 then { fitcurve(P,n)}
          n = -1 } \
        else { nq +=1; Q[nq]: (project(torus(zang,yang))) }
        v = -v }
      yy = yang
      if v >= 0 then {
        n +=1
        P[n]: (project(torus(zang,yang)))
        }
      )
    if n > 2 then { fitcurve(P,n) } 
    }
    fitcurve(Q,nq)
    fitcurve(S,ns)
  }

.PE