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

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.