------------------

to main ; Chandelier by George Mills
3d_frame
perspective
draw
setpc 6 chain
end

to draw
setturtle -1
setxyz 500 500 1000
setturtle 0
repeat 30[rr 12
  repeat 180[
    setpc(se 255 255 255/180*repcount)
    arc2 -1 60]
      bulb
      pu arc2 180 60 pd]
end

to bulb
setpc 4
pu bk 20 pd
repeat 6[ellipse 8 18 rr 60]
pu fd 20 pd
end

to chain
repeat 9[ellipse 5 10
    pu fd 17 pd rr 180]
end

to procs
ed[main draw bulb chain procs]
end
------------------

to main ; 3D_chess by George Mills
3d_frame
perspective
; Get in position so we can see under the cube
setturtle -1 setxyz -500 -500 500
setturtle 0 setpc 6             ; Yellow lines
chess 128 4 heading             ;  Draw face 1
rr 90 chess 128 4 heading       ;  Draw face 2
rt 90 lr 90 chess 128 4 heading ;  Draw face 3
end

to chess :size :level :startheading
if :level<1[stop]
repeat 4[draw :startheading]
setfc ifelse 0=remainder round(heading-:startheading) 180 [4][1]
fill_it
end

to draw :startheading
fd :size rt 90
chess :size/2 :level-1 :startheading
end

to fill_it
pu rt 45 fd 10 pd
fill            ; Red or Blue Checks
pu bk 10 lt 45 pd
end

to procs
ed[main chess draw fill_it procs]
end
------------------

to main ; Hyperboloid by Olga
3d_frame ht
turnon3d 500 500 500
axis
setpc 6 draw 30 90 70
end

to turnon3d :xview :yview :zview
perspective
setturtle -1
setxyz :xview :yview :zview
setturtle 0
end

to axis
setpc 1
pd bk 110 fd 10 seth 90 label "-Y seth 0 fd 100
pd fd 105 pu fd 5 seth 90 label "+Y seth 0
home setpc 2 rt 90
pd bk 110 pu bk 40 label "-X fd 150
pd fd 90 label "+X
home setpc 4 downpitch 90
pd bk 100 pu bk 20 seth 90 label "+Z seth 0 fd 120
pd fd 100 pu fd 15 seth 90 label "-Z seth 0
home pd
end

to draw :a :c :h
for[alf 0 360 9][
   pu setposxyz point1
   pd setposxyz point2
   pu setposxyz point3
   pd setposxyz point4]
end

to point1
op(se  :h/:c*hyp.coor2+hyp.coor1
      -:h
      -:h/:c*hyp.coor1+hyp.coor2)
end

to point2
op(se -:h/:c*hyp.coor2+hyp.coor1
       :h
       :h/:c*hyp.coor1+hyp.coor2)
end

to point3
op(se -:h/:c*hyp.coor2+hyp.coor1
      -:h
       :h/:c*hyp.coor1+hyp.coor2)
end

to point4
op(se  :h/:c*hyp.coor2+hyp.coor1
       :h
      -:h/:c*hyp.coor1+hyp.coor2)
end

to hyp.coor1
op :a*cos :alf
end

to hyp.coor2
op :a*sin :alf
end

to procs
ed[main turnon3d axis draw point1 point2
  point3 point4 hyp.coor1 hyp.coor2 procs]
end
------------------

to main ; Sphere by George Mills
cs ht pu
perspective
clearpalette
ask -3 [setxyz 207 243 97]
setpc[0 0 128]
sphere 80 10
polyview
end

to sphere :rad :step
; Cover the surface of the sphere with polygons
repeat 360/:step[draw 0 rr :step]
end

to draw :i
if :i>180[stop]
down :i
localmake "PointA GetPoint :rad
down :step
localmake "PointB GetPoint :rad
up :step
up :i
rr :step
down :i
localmake "PointD GetPoint :rad
down :step
localmake "PointC GetPoint :rad
up :step
up :i
lr :step
localmake "PointE posxyz
setposxyz :PointA
pd poly pu
setposxyz :PointE
draw :i+:step
end

to poly
polystart
setposxyz :PointB
setposxyz :PointC
setposxyz :PointD
setposxyz :PointA
polyend
end

to GetPoint :rad
fd :rad
localmake "pos posxyz
bk :rad
output :pos
end

to procs
ed[main sphere draw poly getpoint procs]
end
------------------

to main ; Torus by George Mills
perspective
cs ht pu
setsc[192 192 192] setpc 4
; This is the color of the OBJECT not the color you'll see
torus 54 27 6
polyview
end

to torus :rad1 :rad2 :step
; Cover the surface of the torus with polygons
repeat 360/:step[slice :rad1 :rad2 :step rt :step]
end

to slice :rad1 :rad2 :step
; Draw an open ended cylinder
localmake "i 0
repeat 360/:step[draw]
end

to draw
fd :rad1
down :i
localmake "PointA GetPoint :rad2
down :step
localmake "PointB GetPoint :rad2
up :step
up :i
bk :rad1
rt :step
fd :rad1
down :i
localmake "PointD GetPoint :rad2
down :step
localmake "PointC GetPoint :rad2
up :step
up :i
bk :rad1
lt :step
localmake "PointE posxyz
setposxyz :PointA
pd poly pu
setposxyz :PointE
make "i :i+:step
end

to poly
polystart
setposxyz :PointB
setposxyz :PointC
setposxyz :PointD
setposxyz :PointA
polyend
end

to GetPoint :rad
fd :rad
localmake "pos posxyz
bk :rad
output :pos
end

to procs
ed[main torus slice draw poly getpoint procs]
end
------------------

to main ;z=sin_xy by George Mills
cs ht
perspective
setsc [192 192 192] setpc 4
draw
polyview
end

to draw
for[x 0 3 0.1] ~
  [for[y 0 3 0.1][
     localmake "x1 :x+0.0
     localmake "y1 :y+0.0
     localmake "x2 :x+0.1
     localmake "y2 :y+0.1
     localmake "z1 radsin :x1*:y1
     localmake "z2 radsin :x2*:y1
     localmake "z3 radsin :x2*:y2
     localmake "z4 radsin :x1*:y2
     pu
     setxyz :x1*45 :y1*45 :z1*45
     pd
     poly]]
end

to poly
polystart
setxyz :x2*45 :y1*45 :z2*45
setxyz :x2*45 :y2*45 :z3*45
setxyz :x1*45 :y2*45 :z4*45
setxyz :x1*45 :y1*45 :z1*45
polyend
end

to procs
ed[main draw poly procs]
end
------------------

to main [:size "small] ;Orange Peel by Tom Lynn
		       ; Enhanced by George Mills
; Use '(main "large)' to draw a larger version.

local "small
ifelse :size="small [make "small "true][make "small "false]

cs setsc 0
perspective
if :small [ask -1 [setposxyz [-400 -400 -400]]]
setlight [0 0.2] ask -3 [pu setposxyz [100 100 100]]
rt 150 ; Change this to change orientation of peel

setpc [255 160 40]
ifelse :small [orangepeel 60 15][orangepeel 200 15]
; Uncomment the rt 180 in the next line to turn it into a ball.
setpc 4 ;rt 180
ifelse :small [orangepeel 60 15][orangepeel 200 15]
polyview
ifelse :small [
  repeat 36 [ask -1 [fd 35 down 10]]
  ;Create a GIF file
  ;  setactivearea [-90 -90 90 90]
  ;  gifsave "smallpeel.gif
  ][
  repeat 36 [ask -1 [fd 200 down 10]]
  ]
end

to orangepeel :radius :strips
; Tracks a spiralling path down around a sphere.

  localmake "v 0
  localmake "r 0
  localmake "hanglestep 7.5
  localmake "vanglestep :hanglestep/:strips
  localmake "stripheight 360*:vanglestep/:hanglestep
  localmake "oldlowerpos 0
  localmake "oldnextpos 0

  setturtle 0
    ht
    localmake "startpos posxyz
    localmake "startorientation orientation

  setturtle 1
    ht pu
    setposxyz :startpos
    setorientation :startorientation

  setturtle 2
    ht
  until [and not (:v < 180) not (:r < 350)] [
    setturtle 1
      pu
      rt :r
      ifelse :v < 180 [up (:v - 90)][up 90]
      fd :radius
      localmake "nextpos posxyz

      if (:v < 180 + :stripheight/2) [
        bk :radius
        ifelse :v > :stripheight/2
          [ifelse :v < 180 [down :stripheight/2][down :stripheight/2-:v+180]]
          [down :v]
        fd :radius
        localmake "lowerpos posxyz
        ask 2 [
          ; if :oldnextpos has been set then build polygons
          if listp :oldnextpos [
            ifelse (:v + :stripheight/2)< 180 [
              pu
              setposxyz :lowerpos
              pd
             polystart
                setposxyz :nextpos
                setposxyz :oldnextpos
                setposxyz :oldlowerpos
                setposxyz :lowerpos
              polyend
            ] [
              ; Split the quadrilateral into 2 triangles
              ; This guarantees all polygons are planar and convex
              pu
              setposxyz :lowerpos
              pd

              polystart
                setposxyz :nextpos
                setposxyz :oldnextpos
                setposxyz :lowerpos
              polyend

              polystart
                setposxyz :oldnextpos
                setposxyz :oldlowerpos
                setposxyz :lowerpos
              polyend
            ]
          ]

          make "oldnextpos :nextpos
          make "oldlowerpos :lowerpos
        ]
        bk :radius
        ifelse :v > :stripheight/2
          [ifelse :v < 180 [up :stripheight/2][up :stripheight/2-:v+180]]
          [up :v]
        fd :radius
      ]

      bk :radius
      ifelse :v < 180 [down (:v - 90)][down 90]
      lt :r

    setturtle 0
      pu ; ifelse posxyz = [0 0 0][pu][pd] ; Draw main track
      setposxyz :nextpos

    localmake "r :r + :hanglestep
    if :r > 360 [localmake "r :r - 360]
    localmake "v :v + :vanglestep
  ]

  ; Reset turtle 0 to original state
  setturtle 0
    pu
    setposxyz :startpos
    setorientation :startorientation
end
------------------

to main ; Punched_Cubes by Tom Lynn

  cs perspective setsc [0 0 0] setlight [0 0.3] ask -3 [setposxyz [-1000 -1000 400]]
  ht
  setcam [150 150 850]
  ask -2 [setposxyz [150 150 -150]]

  (holeycube 300 100 [95 127 95] [63 95 63] 4)
  setposxyz [100 100 -500] setorientation [0 0 0]
  (holeycube 100 50 [95 127 191] [95 127 191] 4)
  polyview
  flyaround [150 150 -150] 1000 10
end

to angleto :pxyz
; Returns the angle to turn right by in order to point towards :pxyz,
; (assuming that :pxyz is in the same plane).

  localmake "startpos posxyz
  if :startpos = :pxyz [output 0]
  localmake "hvec vecsub :pxyz :startpos
  localmake "penstate pendownp

  pu fd 1 localmake "fdvec vecsub posxyz :startpos bk 1
  rt 90 fd 1 localmake "rtvec vecsub posxyz :startpos bk 1 lt 90
  local "theta
  localmake "costheta (dotprod :hvec :fdvec) / (veclen :hvec)
  if 1 < (abs :costheta) [
    ifelse :costheta < 0 [make "costheta -1][make "costheta 1]
  ]
  make "theta arccos :costheta
  if (dotprod :hvec :rtvec) < 0 [make "theta 360-:theta]

  ifelse :penstate [pd][pu]
  output :theta
end

to dotprod :p :q
  output apply "sum (map "product :p :q)
end

to flyaround :centre :radius :spq
  localmake "a 90/:spq
  localmake "step :radius*(sin :a)

  ; spq = slices per quarter
  ;   a = slice angle

  setturtle -1

  ; Should really go down :a/2 to centre properly,
  ; but this causes problems with the view-flipping.

  ; The "ask -2"s are needed to cope with a view-flipping
  ; effect which I don't really understand.  Comment them
  ; out to see what I mean.

  repeat :spq [fd :step down :a]
  ask -2 [up 90 rt 180]
  repeat :spq [fd :step down :a]
  ask -2 [down 90 rt 180]
  repeat :spq [fd :step down :a]
  ask -2 [up 90 rt 180]
  repeat :spq [fd :step down :a]
end

to hole :size [:slices 8]
; Makes a square with a circular hole of same diameter as square
; side length, :size.  This can be used together with holeface and
; tube to create a 3d object with a cylindrical hole.  Circle is
; represented as a regular polygon with 4*:slices sides.

  localmake "sliceangle 90/:slices
  localmake "holeside :size*sin (:sliceangle/2)
  local [k orient corner lastpos]

  repeat 4 [
    make "corner posxyz
    pu fd :size/2 make "lastpos posxyz bk :size/2
    repeat :slices [
      pu make "orient orientation
      pd polystart
        make "k (angleto :lastpos)
        rt :k fd distancexyz :lastpos lt :k
        make "k 180-(repcount-0.5)*:sliceangle
        rt :k fd :holeside make "lastpos posxyz lt :k
        make "k (angleto :corner)
        rt :k fd distancexyz :corner lt :k
      polyend
    ]
    pu setposxyz :corner setorientation :orient
    fd :size rt 90
  ]
end

to holeface :side :diameter [:slices 8]
; Creates a square of sidelength :side with a central circular hole
; of diameter :diameter.  In fact, the circle is a :slices sided polygon.

  localmake "diag (0.5*sqrt 2)*(:side-:diameter)

  repeat 4 [
    pd polystart
      fd :side
      rt 135 fd :diag
      rt 45 fd :diameter
      rt 45 fd :diag
      rt 135
    polyend
    pu fd :side rt 90
  ]
  pu rt 45 fd :diag lt 45
  (hole :diameter :slices)
  pu rt 45 bk :diag lt 45
end

to holeycube :squareside :holediam :cubecol :tubecol [:slices 8]
  repeat 2 [
    setpc :cubecol (holeface :squareside :holediam :slices)
    pu rt 45 fd 0.5*:squareside*sqrt 2 lt 45
    setpc :tubecol (tube :holediam :squareside :slices)
    pu rt 45 bk 0.5*:squareside*sqrt 2 lt 45
    fd :squareside down 90
    setpc :cubecol pd polystart repeat 4 [fd :squareside rt 90] polyend
    pu fd :squareside down 90
  ]
  pu lt 90 down 90
  repeat 2 [
    pd polystart repeat 4 [fd :squareside rt 90] polyend
    pu repeat 2 [fd :squareside down 90]
  ]
end

to reorient :rpy
; Restores an orientation returned by orientation.

  setorientation [0 0 0]
  rr (first :rpy)
  up (first butfirst :rpy)
  rt (first butfirst butfirst :rpy)
end

to setcam :pxyz
  ask -1 [setposxyz :pxyz]
end

to tube :d :h [:slices 8]
; Draws an open-ended cylindrical tube, diameter :d, height :h
; with the centre of the top around the current turtle position
; and the tube perpendicular below it.  Circular cross-section
; is modelled by a regular polygon with 4*:slices sides.

  localmake "r :d/2
  localmake "halfsliceangle 45/:slices
  localmake "side :d*sin (:halfsliceangle)
  repeat 4*:slices [
    pu fd :r rt :halfsliceangle down 90
    pd polystart
      repeat 2 [fd :h rt 90 fd :side rt 90]
    polyend
    pu up 90 lt :halfsliceangle bk :r
    rt 2*:halfsliceangle
  ]
end

to vecadd :a :b
; Adds vectors :a and :b
  output (map [?1 + ?2] :a :b)
end

to veclen :v
  localmake "a first :v
  localmake "b first butfirst :v
  localmake "c first butfirst butfirst :v
  output sqrt (:a*:a + :b*:b +:c*:c)
end

to vecscale :s :v
; Multiplies vector :v by scalar :s
  output map (list "? "* :s) :v
end

to vecsub :p :q
  output (map [?1 - ?2] :p :q)
end
------------------

to main ; Sombrero by MHElhefni
cs ht setpc 4
perspective
grid 10 .5
ask -3[setxyz 300 200 150]
polyview
end

to grid :a :stp
for[x -:a :a :stp][
  for[z -:a :a :stp][
    polystart
    pd
                       setxyz :x*15 y*15 :z*15
    make "x :x+ :stp   setxyz :x*15 y*15 :z*15
    make "z :z+ :stp   setxyz :x*15 y*15 :z*15
    make "x :x- :stp   setxyz :x*15 y*15 :z*15
    make "z :z- :stp   setxyz :x*15 y*15 :z*15
    polyend]
       pu]
end

to y
localmake "r sqrt :x* :x+:z* :z+.001
op 8*(sin :r*180/pi)/ :r
end

to procs
ed[main grid y procs]
end
------------------

to main ; Sphere2 by MHElhefny
reset setpc 1
perspective
sphere2 45 sphere1 70
shade 80 160 30
end

to sphere2 :a
for[theta 90 180 10][for [fi 0 360 10][
polystart
pd
sppos :a :fi :theta
sppos :a :fi :theta+10
sppos :a :fi+10 :theta+10
sppos :a :fi+10 :theta
polyend
]pu]
end

to sphere1 :a
for[fi 70 180 10][for [theta 0 360 10][
polystart
pd
sppos :a :fi :theta
sppos :a :fi :theta+10
sppos :a :fi+10 :theta+10
sppos :a :fi+10 :theta
polyend
]pu]
end

to sppos :r :fi :theta
setxyz :r*(sin :theta)*cos :fi :r*cos :theta :r*(sin :theta)*sin :fi
end

to shade :c :d :e
ask -3 [setxyz :c :d :e]
polyview
end

to procs
ed[main sphere2 sphere1 sppos shade procs]
end
------------------

to main ; Grail by MHElhefni
reset setpc [0 0 128]
perspective
grid 30 60
shade 200 150 30
polyview
end

to grid :a :b
for[fi 0 360 10][for [r -:b :a 10][
polystart pd
cylpos :r :fi y
make "r :r+ 10    cylpos :r :fi y
make "fi :fi+ 10  cylpos :r :fi y
make "r :r- 10    cylpos :r :fi y
make "fi :fi- 10  cylpos :r :fi y
polyend
pd]pu]
end

to cylpos :r :fi :y
setxyz :r*cos :fi :y-20 :r*sin :fi
end

to shade :f :g :h
ask -3[setxyz :f :g :h]
polyview
end

to y
op ifelse :r >0 [-10*sqrt(:r)][:r* :r/40]
end

to procs
ed[main grid cylpos shade y procs]
end
------------------

to main ; Spiral by MHElhefni
reset setpc[128 0 128]
perspective
grid 15
grid1 15
shade 250 80 20
polyview
end

to grid :step
for[x 0 360 :step][
   for[z 0 620 :step][
polystart pd
                 setxyz x y+18 z
make "x :x+:step setxyz x y+18 z
make "z :z+:step setxyz x y+18 z
make "x :x-:step setxyz x y+18 z
polyend
]pu]
end

to grid1 :step
for[x 0 360 :step][
   for[z 0 700 :step][
polystart pd
                 setxyz x1 y+18 z1
make "x :x+:step setxyz x1 y+18 z1
make "z :z+:step setxyz x1 y+18 z1
make "x :x-:step setxyz x1 y+18 z1
polyend
]pu]
end

to shade :f :g :h
ask -3[setxyz :f :g :h]
polyview
end

to x
op :x*(sin :z)/6
end

to x1
op 130 *(sin :x)/6
end

to y
op -70+ :z/6
end

to z
op :x*(cos :z)/6
end

to z1
op 130*(cos :x)/6
end

to procs
ed[main grid grid1 shade x x1 y z z1 procs]
end
------------------

to main ; Pretzel by Mike Sandy
;NOTE DRAW USES 'THROW' TO TRAP A PARTICULAR ERROR
;reset
setsc [0 0 0]
perspective
setpc [255 0 0]
setlight [.35 .6]
setturtle -3 setposxyz [1000 900 0]
setturtle 0
cs
ht

;PRETZEL
draw [0 360*:rs 360/:n]~
     [:size1*(:dr*(cos :u)+:h*cos (:dr/:rs*:u))]~
     [:size1*(:dr*(sin :u)-:h*sin (:dr/:rs*:u))]~
     [:size1*sin (:fr*:u+180*:lor)]~
     [rs 2 rl 3 n 90 size 42
      lor 1  h .2*:rs dr :rl-:rs size1 :size/:dr
      fr :rl/:rs]~
     [offset [0 0 0] eyepos [0 0 700]]~
     [bradius .25*:size b.incr 15]
; rs AND rl SHOULD BE SMALL, RELATIVELY PRIME INTEGERS. DO NOT PUT rs = rl
; fh (0 - 1) SIZE OF h RELATIVE TO rs
; size CONTROLS SIZE OF PLOT
; n NUMBER OF SEGMENTS/CYCLE FOR WHOLE CURVE
; lor (values 0 or 1) DETERMINES WHETHER KNOT IS LEFT OR RIGHT HANDED

show readword ;PRESS ENTER
cs
;AMMONITE!
draw[-9 2*pi*:n .12]~
    [:size*:c*(exp :a*:u)*((exp :b)+1)*radsin :u]~
    [:size*:c*(exp :a*:u)*((exp :b)+1)*radcos :u]~
    [0]~
    [n 2.2 size 204 c .05 a .1 b :a*2*pi] ;n 2.2 size 250 ~
    [offset [-12 0 0] eyepos [0 -800 1000]]~
    [bradius :size*:c*(exp :a*:u)*((exp :b)-1) b.incr 10]
show readword
cs
;CONICAL SPIRAL
draw[0 2*pi*:n 0.07]~
    [:size*:a*:u/:k*radsin :u]~
    [:size*:a*:u/:k*radcos :u]~
    [:size*:u/:k]~
    [n 2.9 size 57 a .4 c 1 k 2*pi*:c] ;c .7~
    [offset [0 0 -95] eyepos [0 1000 100]]~
    [bradius 1.5*:u b.incr 15]
show readword
cs
;7-KNOT FROM LISSAJOUS CURVE
draw[-.2 2*pi*:k 0.1]~
    [:size*radsin (:a*:u+:b*pi)]~
    [:size*radsin :u]~
    [:size*0.2*radcos :u*7/3]~
    [size 84 k 3 a 2/:k b 0 ]~
    [offset [0 0 0] eyepos [0 0 -1000]]~
    [bradius .15*:size b.incr 15]
show readword
cs

end

to draw :urange.l :xexpr.l :yexpr.l :zexpr.l ~
        :fnpars.l :plotpars.l :bradius.l

; :urange.l - RANGE LIST FOR THE VARIABLE u - [START END INCR]
; :u CAN ONLY APPEAR IN THE PARAMETRIC EQUATIONS AND BRADIUS VALUE
; :xexpr.l etc. - PARAMETRIC EQUNS
; :fnpars - AS LIST [var1 val1 var2 val2.. ]
; A VARIABLE MAY BE DEFINED IN TERMS OF PREVIOUSLY DEFINED VARIABLES/PARAMETERS
; BUT u IS NOT INCLUDED ANYWHERE IN THIS LIST
; :plotpars.l - [offset [  ] eyepos [  ]]
; eyepos TURTLE -1 POSITION
; bradius - RADIUS OF BAND. :u CAN BE INCLUDED IN ITS VALUE
; bradius AND b.incr MUST BE DECLARED
; b.incr band incremental angle in degreeS

; EXCEPT IN THE CASE OF A SIMPLE CURVE PLOT
; IF :BRADIUS IS TOO SMALL COMPARED TO :U THE PLOT IS STOPPED

localmake "eyepos [0 0 1000]                            ;SETS DEFAULT VALUE
localmake "offset [0 0 0]                               ;SETS DEFAULT VALUE
localmake "b.incr run (list last :bradius.l)            ;ASSIGNS BAND INCREMENT
localmake "br.l butlast butlast bf :bradius.l           ;:BR.L bradius FORMULA
local "bradius
assign.val :fnpars.l assign.val :plotpars.l             ;ASSIGN VARIABLE VALUES
setturtle -1 setposxyz :eyepos
setturtle -2 setposxyz [0.1 0 0]
setturtle 0
localmake "u.start run (list first :urange.l)
localmake "u :u.start                                   ;SET :U TO START
make "urange.l bf :urange.l                             ;REMOVE 
localmake "u.incr  run (list last :urange.l)            ;SET UP U INCREMENT
localmake "u.end run butlast :urange.l                  ;SET UP U.END
localmake "iposn vecsum :offset (list x y z)            ;STORES FIRST CURVE POINT
localmake "iposn0 []
localmake "edge2.l [] localmake "edge1.l []             ;STORES FOR BAND EDGE COORD
localmake "fl 0                                         ;FLAG FOR START
pu
repeat (round (:u.end-:u.start)/:u.incr)+1 ~
        [make "u :u+:u.incr
         make "bradius run :br.l                        ;ASSIGN :BRADIUS FROM FORMULA
         localmake "tcoord vecsum :offset (list x  y z) ;FIND NEXT CURVE POINT
         ifelse :bradius=0~
            [setposxyz :iposn pd setposxyz :tcoord pu]   ;CURVE ONLY PLOT ~
            [if (and (abs :u)>8*:u.incr~
                 (sqrt sumsq vecdiff :iposn :tcoord)>1.5*:bradius)~
                 [(throw "ERROR [BRADIUS TOO SMALL, OR U INCR TOO LARGE])]
             edge2 :bradius                              ;GENERATES EDGE VALUES
             if :fl=1[band]
             make "edge1.l :edge2.l make "edge2.l []
             make "fl 1
             make "iposn0 :iposn]
             make "iposn :tcoord
         ]
polyview
erns
end

to assign.val :val.list                                ;SETS UP PARAMETERS
if empty? :val.list [ stop]
make (first :val.list)~
     run (list first bf :val.list)
assign.val bf bf :val.list
end

to edge2 :r                                             ;MAKES LIST OF BAND EDGE COORDS
localmake "norm vecdiff :iposn :tcoord                  ;SPHERICAL COORD OF NORMAL TO BAND
localmake "s (sqrt sumsq :norm)
localmake "theta (atan first :norm first bf :norm)
localmake "phi arccos (last :norm)/:s
localmake "ang 0
repeat (round (360/:b.incr))+1~
       [make "edge2.l
             fput (vecsum :iposn (list xc :ang yc :ang zc :ang))
                  :edge2.l
        make "ang :ang+:b.incr]
Make "edge2.l fput last :edge2.l :edge2.l
end

to edge1 :r                                             ;AS edge2 BUT START PI OUT OF PHASE
                                                        ;CORRECTS ANOMALY
localmake "norm vecdiff :iposn0 :iposn
localmake "s (sqrt sumsq :norm)
localmake "theta (atan first :norm first bf :norm)
localmake "phi arccos (last :norm)/:s                   ;CHANGE THE OVERLAP BETWEEN BANDS HERE
localmake "ang 180
make "edge1.l []
repeat (round (360/:b.incr))+1~
       [make "edge1.l
              fput (vecsum :iposn0 (list xc :ang yc :ang zc :ang))
                   :edge1.l
        make "ang :ang+:b.incr]
Make "edge1.l fput last :edge1.l :edge1.l
end

to band
;localmake "norm1 vecdiff :iposn0 :iposn
local [p2 q2]
localmake "p1 first :edge1.l localmake "q1 first :edge2.l
localmake "dist (sqrt sumsq vecdiff :p1 :q1)
if :dist>2*:bradius                                     ;ALLOWS FOR OUT OF PHASE CONDITION~
           [edge1 :bradius
            make "p1 first :edge1.l make "q1 first :edge2.l]
hband bf :edge1.l bf :edge2.l
end

to hband :l1 :l2
if empty? :l1[stop]
make "p2 first :l1
make "q2 first :l2
rect
make "p1 :p2 make "q1 :q2
hband bf :l1 bf :l2
end

to sumsq :vec
if empty? :vec [op 0 stop]
localmake "1stel first :vec
op :1stel*:1stel+sumsq bf :vec
end

to rect
polystart
pu
setposxyz :p1
pd
setposxyz :p2
setposxyz :q2
setposxyz :q1
setposxyz :p1
pu
polyend
end

to txz :r :ang
op (list :r*(cos :ang)*cos :phi :r*(sin :ang))         ;X,Z FOR BAND AFTER Z ROTN
end

to xc :ang                                             ;X BAND AFTER Y ROTN
op (first txz :r :ang)*(cos :theta)-(last txz :r :ang)*sin :theta
end

to yc :ang
op (first txz :r :ang)*(sin :theta)+(last txz :r :ang)*cos :theta
end

to zc :ang
op  -:r*(cos :ang)*sin :phi
end

to swap
make "p1 :p0
make "p3 :p2
end

to x                                                    ;X FOR CURVE
op run :xexpr.l
end

to y
op run :yexpr.l
end

to z
op run :zexpr.l
end

to vecsum :v0 :v1
if empty? :v0 [[]stop]
op fput (first :v1)+(first :v0) vecsum bf :v0 bf :v1
end

to vecdiff :v0 :v1
if empty? :v0 [[]stop]
op fput (first :v1)-(first :v0) vecdiff bf :v0 bf :v1
end

to atan :x :y
if (and :x=0 :y=0)[op 0 stop]
if (and :x=0 :y>0 )[op 90 stop]
if (and :x=0 :y<0)[op -90 stop]
op (arctan :x :y)
end

to procs
ed[main draw  assign.val edge2 edge1 band hband sumsq rect
   txz xc yc zc  swap x y z vecsum vecdiff atan procs]
end
------------------

to main ; Earth by Doug Merrill
;Interval, step and delay can be changed.  Interval is the angle of each turn as it rotates.
;Step is the number of degrees on each square on the globe.  Delay is the delay between frames
;in the movie measured in milliseconds.

   make "interval 10
   make "step 5
   make "delay 0

   perspective
   cs
   setsc [192 192 192]
   clearpalette
   ask -3 [setxyz 207 243 97] ;sets position of light source
   ;ask -1 [setxyz -600 0 0] ;sets viewpoint - default is 400 400 600

   ht
   pu

   ; Time drawing
   localmake "start timemilli
   ;make "append "False
   ;To make a move un-comment the next few lines and start with rot=0
   make "rot 180 ;centered on Africa
   ;repeat 360/:interval ~
   ;[
   Sphere 85 :step :rot ;earth4 changed radius to fit in frame
   ;make "rot :rot + :interval
   (gifsave "earth4.gif :delay :append 0)
   ;make "append "True
   ;]

   ; Display elapsed drawing
   (print "Drawing "time (timemilli - :start) / 1000)
end

to GetPoint :rad
   fd :rad
   localmake "pos posxyz
   bk :rad
   output :pos
end

to Land :long :lat
;These if statements indicate location of land
;When defining a rectange using inequalities add 1 to the lat farther to the north
;and subtract 1 from the long farther to the west Antarctica
     if (and (and :lat > -90 :lat < -69) (and :long > -1 :long < 60)) [setpc :color2 stop]
     if (and (and :lat > -90 :lat < -64) (and :long > 59 :long < 165)) [setpc :color2 stop]
     if (and (and :lat > -90 :lat < -74) (and :long > 164 :long < 300)) [setpc :color2 stop]
     if (and (and :lat > -90 :lat < -64) (and :long > 299 :long < 305)) [setpc :color2 stop]
     if (and :lat = 60 :long = 305) [setpc :color2 stop]
     if (and :lat = -85 (and :long > 304 :long < 330)) [setpc :color2  stop]
     if (and (and :lat > -90 :lat < -74) (and :long > 329 :long < 360)) [setpc :color2 stop]

     ;North America
     if (and (and :lat > 60 :lat < 71) (and :long > 204 :long < 235)) [setpc :color2 stop]
     if (and (and :lat > 55 :lat < 71) (and :long > 234 :long < 240)) [setpc :color2 stop]
     if (and (and :lat > 50 :lat < 76) (and :long > 239 :long < 245)) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 81) (and :long > 244 :long < 270)) [setpc :color2 stop]
     if (and (and :lat > 30 :lat < 41) (and :long > 249 :long < 290)) [setpc :color2 stop]
     if (and (and :lat > 65 :lat < 86) (and :long > 269 :long < 290)) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 66) :long = 270) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 56) (and :long > 274 :long < 285)) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 51) :long = 285) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 56) :long = 290) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 61) (and :long > 294 :long < 305)) [setpc :color2 stop]
     if (and :lat = 70 (and :long > 289 :long < 300)) [setpc :color2 stop]
     if (and :lat = 55 :long = 305) [setpc :color2 stop]
     if (and :lat = 30 (and :long > 254 :long < 270)) [setpc :color2 stop]
     if (and (and :lat > 15 :lat < 26) (and :long > 259 :long < 270)) [setpc :color2 stop]
     if (and :lat = 30 :long = 285) [setpc :color2 stop]
     if (and :lat = 20 (and :long > 269 :long < 280)) [setpc :color2 stop]
     if (and :lat = 15 :long = 280) [setpc :color2 stop]
     if (and :lat = 10 :long = 285) [setpc :color2 stop]

     ;Greenland
     if (and (and :lat > 75 :lat < 86) (and :long > 299 :long < 345)) [setpc :color2 stop]
     if (and :lat = 74 :long = 310) [setpc :color2 stop]
     if (and (and :lat > 65 :lat < 76) :long = 315) [setpc :color2 stop]
     if (and (and :lat > 65 :lat < 76) (and :long > 324 :long < 340)) [setpc :color2 stop]

     ;South America
     if (and :lat = 0 :long = 285) [setpc :color2 stop]
     if (and (and :lat > -15 :lat < 16) (and :long > 289 :long < 300)) [setpc :color2 stop]
     if (and (and :lat < -14 :lat > -55) :long = 295) [setpc :color2 stop]
     if (and (and :lat > -30 :lat < 1) (and :long > 299 :long < 320)) [setpc :color2 stop]
     if (and (and :lat < 11 :lat > 0) :long = 300) [setpc :color2 stop]
     if (and :lat = 5 (and :long > 304 :long < 315)) [setpc :color2 stop]
     if (and (and :lat > -25 :lat < 1) (and :long > 319 :long < 330)) [setpc :color2 stop]
     if (and (and :lat > -40 :lat < -29) (and :long > 299 :long < 315)) [setpc :color2 stop]

     ;Africa and Saudi Arabia
     if (and (and :lat > 5 :lat < 26) (and :long > -1 :long < 55)) [setpc :color2 stop]
     if (and (and :lat > 5 :lat < 31) :long = 355) [setpc :color2 stop]
     if (and (and :lat > 10 :lat < 26) :long = 350) [setpc :color2 stop]
     if (and (and :lat > 25 :lat < 36) (and :long > -1 :long < 25)) [setpc :color2 stop]
     if (and :lat = 30 (and :long > 24 :long < 55)) [setpc :color2 stop]
     if (and (and :lat > 15 :lat < 26) :long = 55) [setpc :color2 stop]
     if (and :lat = 25 :long = 60) [setpc :color2 stop]
     if (and :lat = 10 :long = 55) [setpc :color2 stop]
     if (and :lat = 5 (and :long > 4 :long < 55)) [setpc :color2 stop]
     if (and (and :lat > -15 :lat < 1) (and :long > 19 :long < 50)) [setpc :color2 stop]
     if (and (and :lat > -25 :lat < -14) (and :long > 19 :long < 45)) [setpc :color2 stop]
     if (and (and :lat > -35 :lat < -24) (and :long > 24 :long < 35)) [setpc :color2 stop]
     if (and :lat = -25 :long = 35) [setpc :color2 stop]
     if (and (and :lat > -30 :lat < -19) :long = 50) [setpc :color2 stop]

     ;Europe
     if (and :lat = 45 :long = 355) [setpc :color2 stop]
     if (and (and :lat > 35 :lat < 46) :long = 0) [setpc :color2 stop]
     if (and (and :lat > 40 :lat < 51) :long = 5) [setpc :color2 stop]
     if (and :lat = 60 :long = 0) [setpc :color2 stop]
     if (and (and :lat > 45 :lat < 56) (and :long > 9 :long < 45)) [setpc :color2 stop]
     if (and :lat = 45 (and :long > 19 :long < 35)) [setpc :color2 stop]
     if (and :lat = 40 (and :long > 34 :long < 60)) [setpc :color2 stop]
     if (and :lat = 35 (and :long > 44 :long < 60)) [setpc :color2 stop]
     if (and :lat = 45 (and :long > 49 :long < 60)) [setpc :color2 stop]
     if (and (and :lat > 45 :lat < 66) (and :long > 44 :long < 60)) [setpc :color2 stop]
     if (and (and :lat > 55 :lat < 71) (and :long > 29 :long < 40)) [setpc :color2 stop]
     if (and (and :lat > 55 :lat < 66) :long = 40) [setpc :color2 stop]
     if (and (and :lat > 55 :lat < 71) :long = 20) [setpc :color2 stop]
     if (and :lat = 65 :long = 15) [setpc :color2 stop]
     if (and :lat = 70 :long = 25) [setpc :color2 stop]

     ;Asia
     if (and (and :lat > 30 :lat < 71) :long = 60) [setpc :color2 stop]
     if (and (and :lat > 25 :lat < 71) (and :long > 64 :long < 125)) [setpc :color2 stop]
     if (and :lat = 75 :long = 75) [setpc :color2 stop]
     if (and :lat = 80 (and :long > 99 :long < 110)) [setpc :color2 stop]
     if (and :lat = 75 (and :long > 89 :long < 135)) [setpc :color2 stop]
     if (and (and :lat > 60 :lat < 71) (and :long > 124 :long < 185)) [setpc :color2 stop]
     if (and :lat = 70 (and :long > 184 :long < 195)) [setpc :color2 stop]
     if (and :lat = 60 :long = 165) [setpc :color2 stop]
     if (and :lat = 60 (and :long > 124 :long < 145)) [setpc :color2 stop]
     if (and :lat = 45 (and :long > 124 :long < 140)) [setpc :color2 stop]
     if (and (and :lat > 45 :lat < 56) (and :long > 124 :long < 150)) [setpc :color2 stop]
     if (and :lat = 40 :long = 145) [setpc :color2 stop]
     if (and :lat = 35 :long = 140) [setpc :color2 stop]
     if (and :lat = 25 (and :long > 74 :long < 120)) [setpc :color2 stop]
     if (and :lat = 20 (and :long > 79 :long < 90)) [setpc :color2 stop]
     if (and :lat = 15 :long = 80) [setpc :color2 stop]
     if (and :lat = 10 :long = 85) [setpc :color2 stop]
     if (and :lat = 15 :long = 110) [setpc :color2 stop]
     if (and :lat = 20 (and :long > 99 :long < 115)) [setpc :color2 stop]
     if (and (and :lat > 0 :lat < 16) :long = 100) [setpc :color2 stop]
     if (and :lat = 15 :long = 130) [setpc :color2 stop]
     if (and (and :lat > -5 :lat < 6) (and :long > 114 :long < 125)) [setpc :color2 stop]

     ;Australia
     if (and :lat = 0 :long = 145) [setpc :color2 stop]
     if (and :lat = -5 :long = 150) [setpc :color2 stop]
     if (and :lat = -40 :long = 175) [setpc :color2 stop]
     if (and :lat = -35 :long = 180) [setpc :color2 stop]
     if (and (and :lat > -30 :lat < -14) (and :long > 129 :long < 155)) [setpc :color2 stop]
     if (and (and :lat > -35 :lat < -19) (and :long > 119 :long < 130)) [setpc :color2 stop]
     if (and :lat = -30 :long = 130) [setpc :color2 stop]
     if (and :lat = -10 :long = 135) [setpc :color2 stop]
     if (and (and :lat > -40 :lat < -19) :long = 155) [setpc :color2 stop]
     if (and (and :lat > -40 :lat < -29) :long = 150) [setpc :color2 stop]
     if (and :lat = -30 :long = 145) [setpc :color2 stop]

     ; the following are junk
     ;if (and :lat = 0 :long = 285) [setpc :color2 stop]
     ;if (and (and :lat > 0 :lat < 31) :long = 200) [setpc :color2 stop]
     ;if (and :lat = 0 (and :long > 179 :long < 220)) [setpc :color2 stop]
     ;if (and (and :lat > 0 :lat < 51) (and :long > 259 :long < 320)) [setpc :color2 stop]
     ;if (and :lat = -25 :long = 35) [setpc :color2 stop]
end

to Slice :rad :step
   make "color1 [0 0 128]
   make "color2 [185 139 97]
   ; Draw an "orange" slice (just the outside surface)
   make "i 0
   make "lat 90
   repeat 180/:step [sliceguts :rad :step ]
end

to Sliceguts :rad :step
     down :i
     localmake "PointA GetPoint :rad
     down :step
     localmake "PointB GetPoint :rad
     up :step
     up :i
     rr :step
     down :i
     localmake "PointD GetPoint :rad
     down :step
     localmake "PointC GetPoint :rad
     up :step
     up :i
     lr :step
     localmake "PointE posxyz
     setposxyz :PointA
     setpc :color1
     Land :long :lat
     pd
     polystart
     setposxyz :PointB
     setposxyz :PointC
     setposxyz :PointD
     setposxyz :PointA
     polyend
     polyview
     pu
     setposxyz :PointE
     make "i :i + :step
     make "lat :lat - :step
end

to Sphere :rad :step :rot
   cs
   make "long 0 - :rot
   ; Cover the surface of the sphere with polygons
   repeat 360/:step ~
      [if :long < 0 [make "long :long + 360]
      if (or :long > 360 :long = 360) [make "long :long - 360]
      Slice :rad :step
      rr :step
      make "long :long + :step
      ]
end
------------------

to main ; Horns by George Mills
perspective
cs setsc [0 0 0]
clearpalette
ask -3 [setxyz 207 243 97]
ht pu
repeat 3~
   [
    repeat 6 [horn 50 10 rt 360/6]
    rr 360/6
   ]
polyview
pd
end

to horn :rad :step
; Cover the surface of the horn with polygons
; (progressing from narrow end to wide end)
; The 12 here in the repeat is tied to the
; *12 after the loop and the 120 in setpc
repeat 12~
   [
    ; This formula for choosing color assumes
    ; 12 loops
    setpc (list 20 120-repcount*10 repcount*10)
    slice (repcount-1)*2 :step
    fd :step
   ]
bk :step*12
end

to slice :rad :step
; this generates one slice accross the horn
; (a tappered cyclinder)
repeat 360/:step ~
   [
    down 90
    localmake "PointA GetPoint func :rad
    rt :step
    localmake "PointB GetPoint func :rad
    lt :step
    up 90

    fd :step

    down 90
    localmake "PointD GetPoint func :rad+2
    rt :step
    localmake "PointC GetPoint func :rad+2
    lt :step
    up 90

    bk :step

    localmake "PointE posxyz
    setposxyz :PointA
    pd
    polystart
    setposxyz :PointB setposxyz :PointC
    setposxyz :PointD setposxyz :PointA
    polyend
    pu
    setposxyz :PointE
    rr :step
   ]
end

to func :x
; This is the main function that creates the horn
; shape (which is effectively rotated in 3D)
; Note that the 12 here and in the repeat argument
; in HORN are tied together and really
; should all be functions of one another.
output :x + power 1.3 :x-12
end

to GetPoint :rad
fd :rad
localmake "pos posxyz
bk :rad
output :pos
end

to procs
ed[main horn slice func GetPoint procs]
end
------------------

to main ; Icosahedron by George Mills
; An Icosahedron built by connecting vertices of 3 perpendicular golden rectangles
; Inspired by Stan Munson
cs
perspective
setsc 0
setpc [255 0 0]
ask -1 [setxyz 100 200 300]
ask -3 [setxyz 0 455 834]

; draw 3 golden rectangles perpendicular to one another (returning the vertices of each)

localmake "set1 goldenrectangle 100
rr 90
lt 90
localmake "set2 goldenrectangle 100
rr 90
lt 90
localmake "set3 goldenrectangle 100
rr 90
lt 90

; Now connect the vertcies
connectsets :set1 :set2 :set3

; Render it
polyview
end

to connectpair :p1 :p2
; This is a cool routine, it will connect any two points in space with a cylinder (a pipe)

localmake "rad 1    ; radius of pipe
localmake "sides 16 ; number of sides of pipe
local "q1
local "q2
local "q3
local "q4
pu

; Position our selves at one end pointing directly at the other
setposxyz :p1
setorientation towardsxyz :p2
localmake "dist distancexyz :p2

; draw panels (polygons) around the center line bewteen p1 and p2
repeat :sides~
   [
   ; Go to the edge of the cylinder and record the coordinate
   down 90
   fd :rad
   make "q1 posxyz
   bk :rad
   up 90

   ; Now go to the other end of teh cylinder and do the same
   fd :dist
   down 90
   fd :rad
   make "q2 posxyz
   bk :rad
   up 90
   bk :dist

   ; Now roll (one panel) and record two more coordinates the same way
   rr 360/:sides

   ; Go to the edge of the cylinder and record the coordinate
   down 90
   fd :rad
   make "q3 posxyz
   bk :rad
   up 90

   ; Now go to the other end of teh cylinder and do the same
   fd :dist
   down 90
   fd :rad
   make "q4 posxyz
   bk :rad
   up 90
   bk :dist

   ; Now connect the 4 vertices into a polygon
   setposxyz :q1
   pd
   polystart
   setposxyz :q2
   setposxyz :q4
   setposxyz :q3
   setposxyz :q1
   pu
   polyend

   ; Get back into the center for the next panel
   setposxyz :p1
   ]
pd
end

to connectpentagon :p1 :p2 :p3 :p4 :p5
; Make a polygon from the 5 points
   setposxyz :p1
   pd
   polystart
   setposxyz :p2
   setposxyz :p3
   setposxyz :p4
   setposxyz :p5
   setposxyz :p1
   pu
   polyend
end

to connectpentagons :set1 :set2 :set3
setpc [0 100 0]
; These connect all the vertices of each pentagon (little bit of trial and error)

connectpentagon item 1 :set1 item 2 :set1 item 2 :set3 item 1 :set2 item 3 :set3
connectpentagon item 1 :set1 item 2 :set1 item 1 :set3 item 4 :set2 item 4 :set3

connectpentagon item 3 :set1 item 4 :set1 item 4 :set3 item 3 :set2 item 1 :set3
connectpentagon item 3 :set1 item 4 :set1 item 3 :set3 item 2 :set2 item 2 :set3

connectpentagon item 1 :set2 item 2 :set2 item 2 :set1 item 1 :set3 item 3 :set1
connectpentagon item 1 :set2 item 2 :set2 item 1 :set1 item 4 :set3 item 4 :set1

connectpentagon item 3 :set2 item 4 :set2 item 4 :set1 item 3 :set3 item 1 :set1
connectpentagon item 3 :set2 item 4 :set2 item 3 :set1 item 2 :set3 item 2 :set1

connectpentagon item 1 :set3 item 2 :set3 item 2 :set2 item 1 :set1 item 3 :set2
connectpentagon item 1 :set3 item 2 :set3 item 1 :set2 item 4 :set1 item 4 :set2

connectpentagon item 3 :set3 item 4 :set3 item 4 :set2 item 3 :set1 item 1 :set2
connectpentagon item 3 :set3 item 4 :set3 item 3 :set2 item 2 :set1 item 2 :set2
end

to connectsets :set1 :set2 :set3
; These connect all the vertices of the icosahedron (little bit of trial and error)
setpc [0 0 255]

; Set 1

connectpair item 1 :set1 item 2 :set1
connectpair item 3 :set1 item 4 :set1

connectpair item 1 :set1 item 3 :set3
connectpair item 1 :set1 item 4 :set3

connectpair item 2 :set1 item 1 :set3
connectpair item 2 :set1 item 2 :set3

connectpair item 4 :set1 item 3 :set3
connectpair item 4 :set1 item 4 :set3

connectpair item 3 :set1 item 1 :set3
connectpair item 3 :set1 item 2 :set3

; Set 2

connectpair item 1 :set2 item 2 :set2
connectpair item 3 :set2 item 4 :set2

connectpair item 1 :set2 item 3 :set1
connectpair item 1 :set2 item 4 :set1

connectpair item 2 :set2 item 1 :set1
connectpair item 2 :set2 item 2 :set1

connectpair item 4 :set2 item 3 :set1
connectpair item 4 :set2 item 4 :set1

connectpair item 3 :set2 item 1 :set1
connectpair item 3 :set2 item 2 :set1

; Set 3

connectpair item 1 :set3 item 2 :set3
connectpair item 3 :set3 item 4 :set3

connectpair item 1 :set3 item 3 :set2
connectpair item 1 :set3 item 4 :set2

connectpair item 2 :set3 item 1 :set2
connectpair item 2 :set3 item 2 :set2

connectpair item 4 :set3 item 3 :set2
connectpair item 4 :set3 item 4 :set2

connectpair item 3 :set3 item 1 :set2
connectpair item 3 :set3 item 2 :set2

end

to goldenratio
output 1/2*(1 + sqrt 5)
end

to goldenrectangle :l
  ; Draw a golden rectangle :l by :l*goldenratio
  ; This routine looks more complex because it draws it with respect to
  ; its center and it also records the vertices to be output.
  localmake "vertices []
  pu
  fd :l*goldenratio/2
  lt 90
  fd :l/2
  rt 180
  pd
  polystart
  repeat 2~
    [
    make "vertices lput posxyz :vertices
    fd :l
    rt 90
    make "vertices lput posxyz :vertices
    fd :l*goldenratio
    rt 90
    ]
  polyend
  pu
  rt 180
  bk :l/2
  rt 90
  bk :l*goldenratio/2
  pd
  output :vertices
end
------------------

to main ; Pentagons by George Mills
; An Icosahedron built by connecting vertices of
; 3 perpendicular golden rectangles
; Inspired by Stan Munson
; This version shows the 12 pentagons inside
;  the Icosahedron instead of the rectangles
cs
perspective
setsc 0
setpc [100 0 0]
ask -1 [setposxyz [100 200 300]]
ask -3 [setposxyz [72 65 164]]

; draw 3 golden rectangles perpendicular to one another
; (returning the vertices of each)

localmake "set1 goldenrectangle 100
rr 90
lt 90
localmake "set2 goldenrectangle 100
rr 90
lt 90
localmake "set3 goldenrectangle 100
rr 90
lt 90

; Now connect the vertcies with pipes
connectsets :set1 :set2 :set3

; Now shade the internal pentagons
connectpentagons :set1 :set2 :set3

; Render it
polyview
;messagebox [] []
;repeat 72 [ask -1 [fd 50 down 5] polyview]
;repeat 100 [ask -3 [setxyz random 300 random 300 random 300]
;            ask -3 [show posxyz] polyview wait 60]
end

to connectpair :p1 :p2
; This is a cool routine, it will connect any
; two points in space with a cylinder (a pipe)

localmake "rad 1    ; radius of pipe
localmake "sides 16 ; number of sides of pipe
local "q1
local "q2
local "q3
local "q4
pu

; Position our selves at one end pointing directly at the other
setposxyz :p1
setorientation towardsxyz :p2
localmake "dist distancexyz :p2

; draw panels (polygons) around the center line bewteen p1 and p2
repeat :sides~
   [
   ; Go to the edge of the cylinder and record the coordinate
   down 90
   fd :rad
   make "q1 posxyz
   bk :rad
   up 90

   ; Now go to the other end of teh cylinder and do the same
   fd :dist
   down 90
   fd :rad
   make "q2 posxyz
   bk :rad
   up 90
   bk :dist

   ; Now roll (one panel) and record two more coordinates the same way
   rr 360/:sides

   ; Go to the edge of the cylinder and record the coordinate
   down 90
   fd :rad
   make "q3 posxyz
   bk :rad
   up 90

   ; Now go to the other end of teh cylinder and do the same
   fd :dist
   down 90
   fd :rad
   make "q4 posxyz
   bk :rad
   up 90
   bk :dist

   ; Now connect the 4 vertices into a polygon
   setposxyz :q1
   pd
   polystart
   setposxyz :q2
   setposxyz :q4
   setposxyz :q3
   setposxyz :q1
   pu
   polyend

   ; Get back into the center for the next panel
   setposxyz :p1
   ]
pd
end

to connectpentagon :p1 :p2 :p3 :p4 :p5
; Make a polygon from the 5 points
   setposxyz :p1
   pd
   polystart
   setposxyz :p2
   setposxyz :p3
   setposxyz :p4
   setposxyz :p5
   setposxyz :p1
   pu
   polyend
end

to connectpentagons :set1 :set2 :set3
setpc [0 100 0]
; These connect all the vertices of each pentagon (little bit of trial and error)

connectpentagon item 1 :set1 item 2 :set1 item 2 :set3 item 1 :set2 item 3 :set3
connectpentagon item 1 :set1 item 2 :set1 item 1 :set3 item 4 :set2 item 4 :set3

connectpentagon item 3 :set1 item 4 :set1 item 4 :set3 item 3 :set2 item 1 :set3
connectpentagon item 3 :set1 item 4 :set1 item 3 :set3 item 2 :set2 item 2 :set3

connectpentagon item 1 :set2 item 2 :set2 item 2 :set1 item 1 :set3 item 3 :set1
connectpentagon item 1 :set2 item 2 :set2 item 1 :set1 item 4 :set3 item 4 :set1

connectpentagon item 3 :set2 item 4 :set2 item 4 :set1 item 3 :set3 item 1 :set1
connectpentagon item 3 :set2 item 4 :set2 item 3 :set1 item 2 :set3 item 2 :set1

connectpentagon item 1 :set3 item 2 :set3 item 2 :set2 item 1 :set1 item 3 :set2
connectpentagon item 1 :set3 item 2 :set3 item 1 :set2 item 4 :set1 item 4 :set2

connectpentagon item 3 :set3 item 4 :set3 item 4 :set2 item 3 :set1 item 1 :set2
connectpentagon item 3 :set3 item 4 :set3 item 3 :set2 item 2 :set1 item 2 :set2
end

to connectsets :set1 :set2 :set3
; These connect all the vertices of the icosahedron (little bit of trial and error)
setpc [0 0 100]

; Set 1

connectpair item 1 :set1 item 2 :set1
connectpair item 3 :set1 item 4 :set1

connectpair item 1 :set1 item 3 :set3
connectpair item 1 :set1 item 4 :set3

connectpair item 2 :set1 item 1 :set3
connectpair item 2 :set1 item 2 :set3

connectpair item 4 :set1 item 3 :set3
connectpair item 4 :set1 item 4 :set3

connectpair item 3 :set1 item 1 :set3
connectpair item 3 :set1 item 2 :set3

; Set 2

connectpair item 1 :set2 item 2 :set2
connectpair item 3 :set2 item 4 :set2

connectpair item 1 :set2 item 3 :set1
connectpair item 1 :set2 item 4 :set1

connectpair item 2 :set2 item 1 :set1
connectpair item 2 :set2 item 2 :set1

connectpair item 4 :set2 item 3 :set1
connectpair item 4 :set2 item 4 :set1

connectpair item 3 :set2 item 1 :set1
connectpair item 3 :set2 item 2 :set1

; Set 3

connectpair item 1 :set3 item 2 :set3
connectpair item 3 :set3 item 4 :set3

connectpair item 1 :set3 item 3 :set2
connectpair item 1 :set3 item 4 :set2

connectpair item 2 :set3 item 1 :set2
connectpair item 2 :set3 item 2 :set2

connectpair item 4 :set3 item 3 :set2
connectpair item 4 :set3 item 4 :set2

connectpair item 3 :set3 item 1 :set2
connectpair item 3 :set3 item 2 :set2

end

to goldenratio
output 1/2*(1 + sqrt 5)
end

to goldenrectangle :l
  ; Draw a golden rectangle :l by :l*goldenratio
  ; This routine looks more complex because it draws it with respect to
  ; its center and it also records the vertices to be output.
  ; Note that this version does not make polygons
  ; of the rectangles (polystart/polyend
  ; is commented out).
  localmake "vertices []
  pu
  fd :l*goldenratio/2
  lt 90
  fd :l/2
  rt 180
  pd
; polystart
  repeat 2~
    [
    make "vertices lput posxyz :vertices
    fd :l
    rt 90
    make "vertices lput posxyz :vertices
    fd :l*goldenratio
    rt 90
    ]
; polyend
  pu
  rt 180
  bk :l/2
  rt 90
  bk :l*goldenratio/2
  pd
  output :vertices
end
------------------

to main :c :fsize :st.num ; Polyhedra by Mike Sandy
;IF YOU DO NOT HAVE "LET IN YOUR LOGOLIB
;THEN REMOVE "MY. FROM THE .MACRO MY.LET
; :c (0 or 1)
; :fsize (1 gives basic plot size)
; :st.num polyhedron number
perspective
ask -1 [setposxyz [950 -50 700]]
ask -3 [setposxyz [1000 800 200]]
setturtle 1
setlight [.15 .45]
setsc 0 ;[245 245 245]
cs
ht
 let[[tau (1+sqrt 5)/2][pr1 (list 1+:tau 1+:tau )]
     [pr2 (list  -1*(1+2*:tau) 1+:tau)][pr3 (list 1+:tau -1*(1+2*:tau))]
     [q1 (list -1*:tau 1+:tau)][q2 (list 0 -1*:tau)]
     [q3 (list 1+:tau 0)][qi1 (list 1+:tau (-1)*:tau )]
     [qi2 (list 0 1+:tau)][qi3 (list -1*:tau 0)]
     [ic1 1.17082039324994][ic2 0.723606797749979]
     [stellns [[[1][]][[2][]][[3 4][]][[5 6 7][5 6]][[8 9 10][9 10]]
             [[11 12][]][[13][]][[3 5][5]][[5 6 9 10][5 6 9 10]]
             [[10 12][10]][[3 6 9 10][6 9 10]][[3 6 9 12][6 9]][[5 6 9 12][5 6 9]]
             [[4 6 7][6]][[7 8][]][[8 9 11][9]][[4 6 8][6]]
             [[4 6 9 11][6 9]][[7 9 11][9]][[4 5][5]][[7 9 10][9 10]]
             [[8 9 12][9]][[4 6 9 10][6 9 10]][[4 6 9 12][6 9]]
             [[7 9 12][9]][[3 6 7][6]][[5 6 8][5 6]][[10 11][10]]
             [[3 6 8][6]][[3 6 9 11][6 9]][[5 9 11][5 9]][[9 10][5 6]]
             [[3 5 9 10][6]][[4 5 9 10][6]][[9 12][5 6 10]][[3 5 9 12][6 10]]
             [[4 5 9 12][6 10]][[8 10 11][5 6 9]][[3 5 8 10 11][6 9]]
             [[4 5 8 10 11][6 9]][[7 10 11][5 6 9]][[3 5 7 10 11][6 9]]
             [[4 5 7 10 11][6 9]][[4 6 7 9 10][5]][[3 6 7 9 10][5]]
             [[5 6 7 9 10][]][[4 6 7 9 12][5 10]][[3 6 7 9 12][5 10]]
             [[5 6 7 9 12][10]][[4 6 8 9 10][5]][[3 6 8 9 10][5]]
             [[5 6 8 9 10][]][[4 6 8 9 12][5 10]][[3 6 8 9 12][5 10]]
             [[5 6 8 9 12][10]][[4 6 11][5 9]][[3 6 10 11][5 9]][[5 6 10 11][9]]
            ]]

     [lines (list (list :pr1 :pr2)(list :pr1 :qi3)(list :pr1 :q2)
               (list :pr1 :pr3)(list :pr2 :pr3)(list :pr2 :qi1)
               (list :pr2 :q3)(list :pr3 :qi2)(list :pr3 :q1)
               (list :q1 :q2)(list :q1 :qi1)(list :q1 :q3)
               (list :q2 :q3)(list :q2 :qi2)(list :q3 :qi3)
               (list :qi1 :qi2)(list :qi1 :qi3)(list :qi2 :qi3))]
     [s.faces1 [[[2 14 9][15 6 3][8 11 7]];EACH SUB.LIST IS A POLYGON
                [[2 9 7][9 6 2][6 3 9][3 8 6][7 2 8][8 7 3]]
                [[3 6 13 17][2 10 18 9][8 12 16 7]]
                [[9 17 14][7 18 11][3 16 15][15 10 6][2 14 12][8 11 13]]
                [[14 17 3][18 9 11][16 7 15]]
                [[18 10 7][16 12 3][17 13 9]]
                [[15 14 17 10][14 11 18 12][11 15 16 13]]
                [[13 17 8 9][10 18 6 7][12 16 2 3]]
                [[16 15 13][17 14 10][18 11 12]]
                [[13 3 9][7 3 12][9 7 10]]
                [[2 16 1][12 3 4][8 17 4][13 9 5][6 18 5][10 7 1]]
                [[17 10 5][18 12 1][16 13 4]]
                [[4 16 5][5 4 10][5 17 1][1 5 12][1 18 4][4 1 13]
                 [4 12 17][5 13 18][1 10 16]]
              ]]
     [s.faces2 [[[2 14 9][15 6 3][8 11 7]]
                [[2 9 7][9 6 2][6 3 9][3 8 6][7 2 8][8 7 3]]
                [[3 6 13 17][2 10 18 9][8 12 16 7]]
                [[9 17 14][7 18 11][3 16 15][15 10 6][2 14 12][8 11 13]]
                [[10 15 2][12 14 8][13 11 6]]
                [[10 6 18][12 2 16][13 8 17]]
                [[15 14 17 10][14 11 18 12][11 15 16 13]]
                [[13 17 8 9][10 18 6 7][12 16 2 3]]
                [[13 16 11][15 10 17][14 12 18]]
                [[17 8 6][16 2 8][2 18 6]]
                [[2 16 1][12 3 4][8 17 4][13 9 5][6 18 5][10 7 1]]
                [[17 10 5][18 12 1][16 13 4]]
                [[4 16 5][5 4 10][5 17 1][1 5 12][1 18 4][4 1 13][4 12 17][5 13 18][1 10 16]]
              ]]
     [ref.pts (list(list (list :ic2 0 -1*:ic1) (list 0 -1*:ic1 -1*:ic2) (list -1*:ic2 0 -1*:ic1))
               (list (list 0 :ic1 -1*:ic2) (list :ic2 0 -1*:ic1) (list -1*:ic2 0 -1*:ic1))
               (list (list :ic1 :ic2 0) (list :ic2 0 -1*:ic1) (list 0 :ic1 -1*:ic2))
               (list (list :ic1 -1*:ic2 0) (list :ic2 0 -1*:ic1) (list :ic1 :ic2 0))
               (list (list :ic2 0 -1*:ic1) (list :ic1 -1*:ic2 0) (list 0 -1*:ic1 -1*:ic2))
               (list (list -1*:ic2 0 -1*:ic1) (list -1*:ic1 :ic2 0) (list 0 :ic1 -1*:ic2))
               (list (list -1*:ic1 -1*:ic2 0) (list -1*:ic2 0 -1*:ic1) (list 0 -1*:ic1 -1*:ic2))
               (list (list -1*:ic1 :ic2 0) (list -1*:ic2 0 -1*:ic1) (list -1*:ic1 -1*:ic2 0))
               (list (list 0 :ic1 -1*:ic2) (list 0 :ic1 :ic2) (list :ic1 :ic2 0))
               (list (list 0 :ic1 :ic2) (list 0 :ic1 -1*:ic2) (list -1*:ic1 :ic2 0))
               (list (list :ic1 :ic2 0) (list :ic2 0 :ic1) (list :ic1 -1*:ic2 0))
               (list (list :ic2 0 :ic1) (list :ic1 :ic2 0) (list 0 :ic1 :ic2))
               (list (list 0 -1*:ic1 :ic2) (list 0 -1*:ic1 -1*:ic2) (list :ic1 -1*:ic2 0))
               (list (list -1*:ic1 -1*:ic2 0) (list 0 -1*:ic1 -1*:ic2) (list 0 -1*:ic1 :ic2))
               (list (list :ic1 -1*:ic2 0) (list :ic2 0 :ic1) (list 0 -1*:ic1 :ic2))
               (list (list 0 -1*:ic1 :ic2) (list -1*:ic2 0 :ic1) (list -1*:ic1 -1*:ic2 0))
               (list (list 0 -1*:ic1 :ic2) (list :ic2 0 :ic1) (list -1*:ic2 0 :ic1))
               (list (list -1*:ic1 -1*:ic2 0) (list -1*:ic2 0 :ic1) (list -1*:ic1 :ic2 0))
               (list (list 0 :ic1 :ic2) (list -1*:ic1 :ic2 0) (list -1*:ic2 0 :ic1))
               (list (list -1*:ic2 0 :ic1) (list :ic2 0 :ic1) (list 0 :ic1 :ic2)))
       ]
     [colors [[255 0 0][91 238 168][0 255 0][0 159 0][0 0 255]
              [103 103 185][255 88 0][255 255 0][52 179 197][255 28 255]
              [158 26 230][0 51 255][251 224 132][106 251 4][255 182 108]
              [147 147 0][210 0 105][165 85 73][58 160 153][140 140 140]]]
     [size 50]
    ]
 plot.poly :c :fsize :st.num

end

to plot.poly :c :fsize :st.num
let[[st.num :st.num-1]]
(ifelse :st.num=0
   [ifelse :c=0[setpc item 12 :colors][setpc item :c :colors]
     set.size :fsize :st.num draw.faces 0 map "1st.last :ref.pts]
   [set.size :fsize :st.num
    let[[f.coeff1  map[s.face.coeff item ? :s.faces1] first item :st.num :stellns]
        [f.coeff2 map[s.face.coeff item ? :s.faces2] last item :st.num :stellns]]
      repeat 20 [ifelse :c=0[setpc item 12 :colors][setpc item repcount :colors]
                 let[[f.coord1 map [(map [poly.coord item repcount :ref.pts ?] ?)] :f.coeff1]
                     [f.coord2 map [(map [poly.coord item repcount :ref.pts ?] ?)] :f.coeff2]
                    ]
                 polyhed.face :f.coord1  polyhed.face :f.coord2]])
end

to polyhed.face :f.list
if emptyp :f.list [stop]
draw.faces 0 first :f.list
polyhed.face   bf :f.list
end

to set.size :f :st.num
let[[set1 [0 1 2 3]]
     [set2 [4 8 9 10 11 12 13
            14 20 21 23 24 25 26 32
            33 34 35 36 37 44 45 46
            47 48 49 ]]
     [set3 [5 6 15 16 17 18 19 22 27
            28 29 30 31 38 39 40 41
            42 43 50 51 52 53 54
            55 56 57 58]]
     ]
 if memberp :st.num :set1[make "size :f*120 stop]
 if memberp :st.num :set2[make "size :f*85 stop]
 (ifelse memberp :st.num :set3[make "size :f*50]
    [make "size :f*25])
end

to centroid :list
 let[[cn count :list]]
 op h.cent :list 0 0 0
end

to h.cent :list :sx :sy :sz
 if emptyp :list[op (list :sx/:cn :sy/:cn :sz/:cn)]
 let[[1st first :list]
     [x first :1st]
     [y first bf :1st]
     [z last :1st]]
 op h.cent bf :list :x+:sx :y+:sy :z+:sz
end

to draw.faces  :col.incr :list ;[f1 f2..]
 if emptyp :list[pu polyview stop]
 polystart
 pu draw.face  first :list
 polyend
 draw.faces    :col.incr bf :list
end

to draw.face :pts.list
 if emptyp :pts.list[pu stop]
 setposxyz scal.mult :size first :pts.list
 pd draw.face bf :pts.list
end

to intersect :list1 :list2
 if emptyp :list1[op[]]
 let[[1st first :list1]]
 ifelse memberp :1st :list2[
   op :1st][op intersect bf :list1 :list2]
end

to line.pr :list
 (op fput (list first :list last :list)
        h.line.pr :list)
end

to h.line.pr :list
 if emptyp bf :list[op[]]
 (op fput (list first :list first bf :list)
       h.line.pr bf :list)
end

to poly.coeff :nlist ;eg [1 2 3]allows for other than triangles
 op h.poly.coeff 1st.last line.pr :nlist;[[c1 c2 c3][c4..]..[c1 c2 c3]]
end

to h.poly.coeff :pr.list
 if emptyp :pr.list[op[]]
 (op fput line.inter first :pr.list
         h.poly.coeff bf :pr.list)
end

to s.face.coeff :face ;[[[c1..]..[c1..]] [[c2..]..[c2..]]...]
 op map "poly.coeff :face
end

to pt.coord :ref.pts :coeff;[c1 c2 c3],coord for single point
 if emptyp :ref.pts [op[0 0 0]]
 op v.sum scal.mult first :coeff first :ref.pts pt.coord bf :ref.pts bf :coeff
end

to poly.coord :ref.pts :poly;SUB.FACE [coeff1 coeff2 coeff3 coeff1]
 if emptyp :poly[op[]]
 (op fput pt.coord :ref.pts first :poly
        poly.coord :ref.pts bf :poly)
end

to s.face.coord :ref.pts :face.num
 (op map [poly.coord :ref.pts ?]
       s.face.coeff item :face.num :s.faces)
end

to det :v1 :v2
 let[[a1 first :v1]
     [b1 last :v1]
     [a2 first :v2]
     [b2 last :v2]]
 op :b1*:a2-:a1*:b2
end

to line.inter :l.pr
 let[[l1 item first :l.pr :lines]
     [l2 item last :l.pr :lines]
     [inter intersect :l1 :l2]]
 (ifelse emptyp :inter
    [let[[p1 first :l1];[m1 n1 r1]
         [p2 last :l1]
         [p3 first :l2]
         [p4 last :l2]
         [p12 v.diff :p1 :p2]
         [p34 v.diff :p3 :p4]
         [p24 v.diff :p2 :p4]
         [delta  det :p12 :p34]
         [a -1*(det :p24 :p34)/:delta]
         [tvec v.sum scal.mult :a v.diff :p1 :p2 :p2]
         [1st first :tvec]
         [2nd last :tvec]]
      op lput 1-:1st-:2nd :tvec]
    [let [[1st first :inter]
          [2nd last :inter]]
     op lput 1-:1st-:2nd :inter])
end

to v.sum :vec1 :vec2
 if emptyp :vec1[op[]]
 (op fput (first :vec1)+first :vec2
          v.sum bf :vec1 bf :vec2)
end

to v.diff :vec1 :vec2
 if emptyp :vec1[op[]]
 (op fput (first :vec1)-first :vec2
          v.diff bf :vec1 bf :vec2)
end

to scal.mult :k :list
 if emptyp :list[op []]
 op fput :k*first :list scal.mult :k bf :list
end

to 1st.last :list
 op lput first :list :list
end

.macro my.let :list
 if empty? :list [op []]
 op (list "local "first first :list
          "make "first first :list
                "run bf first :list
           "let bf :list)
end

to procs
ed[main plot.poly
  polyhed.face set.size
  centroid h.cent
  draw.faces draw.face
  intersect line.pr h.line.pr
  poly.coeff h.poly.coeff s.face.coeff
  pt.coord poly.coord s.face.coord
  det line.inter
  v.sum v.diff scal.mult
  1st.last my.let procs]
end
------------------