  'GEODESICS.INC

  uses glo2\geometry


  function Gauss2D(single x,single z) as single
  =============================================
  static single e=2.718281828 'Euler's number
  single g=-.5*(x*x+z*z)
  return -1+2*pow(e,g)
  end function


  procedure DrawNoisyGauss()
  ==========================

  static as single e,g,x,y,z,nx,ny,nz,sr,st,tt,r,texx,texz,texs,htn,ght
  static as sys    en,i,j,k,v,noiseinit,p1,p2
  '
  '--------------------------------------------
  en=30      'number of quads along axis
  sr=5       'size of surface for x and z axes
  htn=0.25   'height noise
  texs=0.15  'texture map step
  ght=1      'gauss height factor
  '
  st=sr/en   'x and z step
  tt=st*.25  'x and z jitter
  '
  'TEXTURE DATA
  '
  % texsize 0x10000
  '
  static sys vv
  static long tx1[texsize]
  '
  if not vv
    for i=1 to texsize
      v=noiseb i
      mov eax,v
      shl eax,8
      'mov al,v
      'shl eax,8
      mov al,v
      mov v,eax
      'v=&hffff00ff 'set lilac color
      tx1[i]=v
    next
    vv=1
    MakeTexture @tx1,0x100,0x100,texn[2]
  end if
  '
  glBindTexture GL_TEXTURE_2D,texn[2]
  '
  z=sr*.5
  texz=.1
  '
  'lilac coloring
  'single c[4]=>(1,0,1,1)
  'glmaterialfv GL_FRONT,GL_DIFFUSE, c
  
  '
  for i = 1 to en
    '
    x=-sr*.5
    texx=.1
    '
    for j = 1 to en
      '
      '
      'PLOT
      '----
      '
      glBegin GL_QUAD_STRIP
      '
      for k=0 to 1
        '
        p2=i*en+j+k
        p1=p2-en
        '
        y=noise(p1)*htn'+htn*.3
        '
        nx = 0 : ny = y : nz = 0
        '
        y+=Gauss2D(x,z)*ght
        '
        glNormal3f nx, ny, nz
        glTexCoord2f texx,texz
        glVertex3f x+noise(p1)*tt, y, z+noise(p1+1)*tt
        '
        y=noise(p2)*htn'+htn*.3
        '
        nx = 0 : ny = y : nz = 0
        '
        y+=Gauss2D(x,z-st)*ght
        '
        glNormal3f nx, ny, nz
        glTexCoord2f texx,texz+texs
        glVertex3f x+noise(p2)*tt, y, z-st+noise(p2+1)*tt
        'glVertex3f x, y, z
        '
        x+=st
        texx+=texs
        '
      next
      '
      glEnd()
      '
      x-=st
      texx-=texs
      if texx>=1 then texx-=1
      '
    next
    '
    z-=st
    texz+=texs
    if texz>=1 then texz-=1
    '
  next
  '
  end procedure


  procedure DrawGauss()
  =====================
  '
  static as single e,x,y,z,nx,ny,nz
  static as sys    i,j,k 
  '
  e=2.718281828 'Euler's number
  z=2.5
  for i = 1 to 20
    '
    x=-2.5
    glBegin GL_QUAD_STRIP
    '
    for j = 1 to 20
      '
      '
      'PLOT
      '----
      '
      y=Gauss2D(x,z)
      nx = x
      ny = y+1
      nz = z
      '
      glNormal3f nx, ny, nz
      glVertex3f x, y, z
      '
      z-=.25
      '
      y=Gauss2D(x,z)
      nx = x
      ny = y+1
      nz = z
      '
      glNormal3f nx, ny, nz
      glVertex3f x, y, z
      '
      z+=.25
      x+=.25
      '
    next
    '
    glEnd()
    z-=.25
    '
  next
  '
  end procedure



  procedure Spirodome2(double*ev, sys n)
  ======================================
  '
  'n sides eg: 10
  'ev fraction of pi .15 .. .25
  '
  double a,b,r,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
  double a1,a2,a3,a4,r1,r2,r3,r4
  double na,nx,ny,nz
  double ha,cr,rd,w1,w2,l1,l2,s1,s2,h1,h2
  '
  sys e,i,j,k,alt
  '
  'PRELIM CALCULATIONS
  '-------------------
  '
  ha=pi()/n           'half angle
  a=ha*2              'panel step angle
  '
  rd=.5               'set initial half radius (centre of initial panel)
  cr=rd               'radius at mid panel
  r=cr/cos ha         'radius for mid panel points
  w1=r*sin ha         'intial panel half width
  ev*=pi()          'set initial slope
  l1= rd / cos ev     'initial panel half length
  s1=sqr(w1*w1+l1*l1) 'fix side length
  h1=rd*tan ev        'first panel height contribution
  ny=cos ev           'normal y component
                      '
  r1=0                'apex radius
  r2=r                'mid panel points radius
  r3=cr+cr            'panel tail radius
                      '
  y1=2                'initial dome height
  y2=y1-h1            'initial mid panel height
  y3=y2-h1            'initial panel tail height
                      '
  rd=0                'radius step difference
  a1=-ha              'initial face-on displacement
  e=6                 'number of tiers
  '
  'EACH TIER
  '---------
  '
  for j=1 to e
    '
    'EACH PANEL
    '----------
    '
    x1=r2*sin a1
    z1=r2*cos a1
    '
    for i=1 to n
      a2=a1+a
      x2=r2*sin a2
      z2=r2*cos a2
      '
      'face normal
      '
      nx=sin na
      nz=cos na
      '
      'TIP TRIANGLE
      '------------
      '
      a3=a1+ha
      x3=r1*sin a3
      z3=r1*cos a3
      if j=e then
        glbegin GL_TRIANGLES
        glNormal3f nx,ny,nz
        glVertex3f x1,y2,z1
        glVertex3f x2,y2,z2
        glVertex3f x3,y1,z3
        glend
      end if
      '
      'TAIL TRIANGLE
      '-------------
      '
      x4=r3*sin a3
      z4=r3*cos a3
      '
      if j<e then
        glbegin GL_QUADS
        glNormal3f nx,ny,nz
        glVertex3f x3,y1,z3
        glVertex3f x1,y2,z1
        glVertex3f x4,y3,z4
        glVertex3f x2,y2,z2
        glend
      end if
      '
      '
      'ANGLES FOR NEXT PANEL
      '
      a1=a2 : x1=x2 : z1=z2
      na=a1+ha 'angle for next normal
      '
    next
    '
    'next tier
    '---------
    '
    'CALCULATIONS FOR NEXT PANEL SET
    '-------------------------------
    '
    w1=r2*sin ha        'new panel half width
    l1=sqr(s1*s1-w1*w1) 'new panel half length
    r1=r2               'new r1
    r2=r3               'new r2     
    rd=r2*cos(ha)-r1    'rdiff
    r3=r1+rd+rd         'new r3
    ev=acos(rd/l1)      'elevation     
    ny=rd/l1            'normal y
    h1=rd*tan ev        'panel height contribution
    y1=y2               'new y1
    y2=y3               'new y2
    y3=y1-h1-h1         'new y3
    '
    alt=1-alt 'alternate starting angle for next tier
    if alt then a1=0 else a1=-ha
    na=a1+ha
    '
  next
  '
  end procedure


  procedure Inter(vector *v, int n1,n2,n3,n4,float f)
  ===================================================
  'n1 n2 end points
  'n3 n4 interpolation point
  float g=1.0-f
  v[n3].x=v[n1].x+f*(v[n2].x-v[n1].x)
  v[n3].y=v[n1].y+f*(v[n2].y-v[n1].y)
  v[n3].z=v[n1].z+f*(v[n2].z-v[n1].z)
  v[n4].x=v[n1].x+g*(v[n2].x-v[n1].x)
  v[n4].y=v[n1].y+g*(v[n2].y-v[n1].y)
  v[n4].z=v[n1].z+g*(v[n2].z-v[n1].z)
  end procedure


  procedure TetraVecs(vector*v)
  =============================
  '
  float sa=0.5 'sin 30
  float ca=sqrt(0.75) 'cos 30
  v[1]={  0,   1,  0} 'top
  v[2]={ ca, -sa,  0} '0
  v[3]={-sa, -sa,-ca} '120
  v[4]={-sa, -sa, ca} '240
  '
  end procedure


  procedure OctaVecs(vector*v)
  ============================
  '
  v[1]={ 0, 1, 0}
  v[2]={ 1, 0, 0}
  v[3]={ 0, 0,-1}
  v[4]={-1, 0, 0}
  v[5]={ 0, 0, 1}
  v[6]={ 0,-1, 0}
  '
  end procedure



  procedure IcosaVecs(vector*v)
  =============================
  '
  'vector    v[12]
  float x1,x2,x3,y1,y2,y3,z1,z2,z3
  float a,b,r,a1,a3,r1,ha
  int i,j
  '
  'PRELIM CALCULATIONS
  '
  b=atn 2     'angle subtended by 1 tirangle side (1st lateral angle)
  r1=sin b    'radius for all coords except top and bottom
  ha=0.2*pi() 'half step angle (36 degrees)
  a1=-ha      'start longitudinal angle
  a=ha+ha     'longitudinal step angle (72 degrees)
  y1=cos b    'lateral level for coords 2..6
  b=pi()-b    '2nd lateral angle
  y2=cos b    'lateral angle for coords 7..11
  '
  x1=r1*sin a1
  z1=r1*cos a1
  for i=1 to 5
    '
    j=i+5
    '   
    x1=r1*sin a1
    z1=r1*cos a1
    a3=a1+ha
    x3=r1*sin a3
    z3=r1*cos a3
    '
    'icosahedral triangle vectors
    v[i]<=x1,y1,z1
    v[j]<=x3,y2,z3
    x1=x3 : z1=z3
    a1+=a
    '
  next
  v[11]={0.0, 1.0, 0.0} 'top
  v[12]={0.0,-1.0, 0.0} 'bottom
  '
  end procedure

  float rn=1.0
  '
  procedure NormalizeVec(vector*w)
  ======================================
  float f=sqrt( w.x*w.x + w.y*w.y + w.z*w.z )
  if f
    f=rn/f
  endif
  w.x*=f
  w.y*=f
  w.z*=f
  end procedure
  '
  '
  procedure CentreVec(vector*w,*v, int n1,n2,n3)
  ==============================================
  w.x=(v[n1].x+v[n2].x+v[n3].x)/3
  w.y=(v[n1].y+v[n2].y+v[n3].y)/3
  w.z=(v[n1].z+v[n2].z+v[n3].z)/3
  NormalizeVec w
  end procedure


  procedure CubeVecs(vector*w)
  ============================
  vector v[6]
  octavecs v
  CentreVec(w[1],v,1,2,3)
  CentreVec(w[2],v,1,3,4)
  CentreVec(w[3],v,1,4,5)
  CentreVec(w[4],v,1,5,2)
  '
  CentreVec(w[5],v,6,3,2)
  CentreVec(w[6],v,6,4,3)
  CentreVec(w[7],v,6,5,4)
  CentreVec(w[8],v,6,2,5)
  end procedure


  procedure DodecVecs(vector*w)
  =============================
  vector v[12]
  IcosaVecs v
  'upper
  CentreVec(w[1],v,11,1,2)
  CentreVec(w[2],v,11,2,3)
  CentreVec(w[3],v,11,3,4)
  CentreVec(w[4],v,11,4,5)
  CentreVec(w[5],v,11,5,1)
  'upper middle
  CentreVec(w[6],v,2,1,6)
  CentreVec(w[7],v,3,2,7)
  CentreVec(w[8],v,4,3,8)
  CentreVec(w[9],v,5,4,9)
  CentreVec(w[10],v,1,5,10)
  'lower middle
  CentreVec(w[11],v,6,7,2)
  CentreVec(w[12],v,7,8,3)
  CentreVec(w[13],v,8,9,4)
  CentreVec(w[14],v,9,10,5)
  CentreVec(w[15],v,10,6,1)
  'lower
  CentreVec(w[16],v,12,7,6)
  CentreVec(w[17],v,12,8,7)
  CentreVec(w[18],v,12,9,8)
  CentreVec(w[19],v,12,10,9)
  CentreVec(w[20],v,12,6,10)
  end procedure


  procedure RenderTriangle(vector *v1,*v2,*v3)
  ============================================
  single nx,ny,nz
  nx=(v1.x+v2.x+v3.x)*.3333
  ny=(v1.y+v2.y+v3.y)*.3333
  nz=(v1.z+v2.z+v3.z)*.3333
  glbegin GL_TRIANGLES
  glNormal3f nx, ny, nz
  glVertex3f v1.x, v1.y, v1.z
  glVertex3f v2.x, v2.y, v2.z
  glVertex3f v3.x, v3.y, v3.z
  glend
  end procedure


  procedure RenderTriangle(vector *v,int n1,n2,n3)
  ================================================
  RenderTriangle v[n1],v[n2],v[n3]
  end procedure


  procedure RenderSquare(vector *v1,*v2,*v3,*v4)
  ==============================================
  single nx,ny,nz
  nx=(v1.x+v2.x+v3.x+v4.x)*.25
  ny=(v1.y+v2.y+v3.y+v4.y)*.25
  nz=(v1.z+v2.z+v3.z+v4.z)*.25
  glbegin GL_QUADS
  glNormal3f nx, ny, nz
  glVertex3fv v1
  glVertex3fv v2
  glVertex3fv v3
  glVertex3fv v4
  glend
  end procedure


  procedure RenderSquare(vector *v,int n1,n2,n3,n4)
  =================================================
  RenderSquare v[n1],v[n2],v[n3],v[n4]
  end procedure


  procedure RenderPent(vector *v1,*v2,*v3,*v4,v5)
  ===============================================
  single nx,ny,nz
  nx=(v1.x+v2.x+v3.x+v4.x+v5.x)*.2
  ny=(v1.y+v2.y+v3.y+v4.y+v5.y)*.2
  nz=(v1.z+v2.z+v3.z+v4.z+v5.z)*.2
  glbegin GL_POLYGON
  glNormal3f nx, ny, nz
  glVertex3fv v1
  glVertex3fv v2
  glVertex3fv v3
  glVertex3fv v4
  glVertex3fv v5
  glend
  end procedure


  procedure RenderPent(vector *v,int n1,n2,n3,n4,n5)
  ==================================================
  RenderPent v[n1],v[n2],v[n3],v[n4],v[n5]
  end procedure


  procedure RenderCube()
  ======================
  vector w[8]
  CubeVecs w
  RenderSquare w,1,2,3,4 'top
  RenderSquare w,8,7,6,5 'bottom
  RenderSquare w,2,1,5,6 'back
  RenderSquare w,3,2,6,7 'left
  RenderSquare w,4,3,7,8 'front
  RenderSquare w,1,4,8,5 'right
  end procedure


  procedure RenderIcos(int n=0)
  =============================
  'indexbase 1
  vector v[12]
  Icosavecs v
  '
  'TOP
  RenderTriangle v,11,1,2
  RenderTriangle v,11,2,3
  RenderTriangle v,11,3,4
  RenderTriangle v,11,4,5
  RenderTriangle v,11,5,1
  'UPPER MIDDLE
  RenderTriangle v,1,6,2
  RenderTriangle v,2,7,3
  RenderTriangle v,3,8,4
  RenderTriangle v,4,9,5
  RenderTriangle v,5,10,1
  'LOWER MIDDLE
  RenderTriangle v,6,7,2
  RenderTriangle v,7,8,3
  RenderTriangle v,8,9,4
  RenderTriangle v,9,10,5
  RenderTriangle v,10,6,1
  'BOTTOM
  RenderTriangle v,12,6,10
  RenderTriangle v,12,10,9
  RenderTriangle v,12,9,8
  RenderTriangle v,12,8,7
  RenderTriangle v,12,7,6
  end procedure


  procedure RenderDoDec()
  =======================
  vector w[20]
  DoDecVecs w
  RenderPent w,1,2,3,4,5 'top
  '
  RenderPent w,2,1,6,11,7 'upper 1
  RenderPent w,3,2,7,12,8 'upper 2
  RenderPent w,4,3,8,13,9 'upper 3
  RenderPent w,5,4,9,14,10 'upper 4
  RenderPent w,1,5,10,15,6 'upper 5
  '
  RenderPent w,20,16,11,6,15 'lower 1
  RenderPent w,16,17,12,7,11 'lower 2
  RenderPent w,17,18,13,8,12 'lower 3
  RenderPent w,18,19,14,9,13 'lower 4
  RenderPent w,19,20,15,10,14 'lower 5
  '
  RenderPent w,16,17,18,19,20 'bottom
  end procedure


  procedure Pentagonn(vector *v1,*v2,*v3,*v4,*v5)
  ===================================================
  float nx,ny,nz
  nx=(v1.x+v2.x+v3.x+v4.x+v5.x)*.2
  ny=(v1.y+v2.y+v3.y+v4.y+v5.y)*.2
  nz=(v1.z+v2.z+v3.z+v4.z+v5.z)*.2
  glBegin GL_POLYGON
  glnormal3f nx,ny,nz
  glvertex3fv v1
  glvertex3fv v2
  glvertex3fv v3
  glvertex3fv v4
  glvertex3fv v5
  glEnd
  end procedure


  procedure Hexagonn(vector *v1,*v2,*v3,*v4,*v5,*v6)
  ==================================================
  float nx,ny,nz
  nx=(v1.x+v2.x+v3.x+v4.x+v5.x+v6.x)*.1667
  ny=(v1.y+v2.y+v3.y+v4.y+v5.y+v6.y)*.1667
  nz=(v1.z+v2.z+v3.z+v4.z+v5.z+v6.z)*.1667
  glBegin GL_POLYGON
  glnormal3f nx,ny,nz
  glvertex3fv v1
  glvertex3fv v2
  glvertex3fv v3
  glvertex3fv v4
  glvertex3fv v5
  glvertex3fv v6
  glEnd
  end procedure


  procedure RenderBuckyBall(float f=0.33)
  =======================================
  single nx,ny,nz
  vector v[80] '20 60
  Icosavecs v '1 to 12
  '
  'TOP
  Inter v,11,1,21,26,f
  Inter v,11,2,22,27,f
  Inter v,11,3,23,28,f
  Inter v,11,4,24,29,f
  Inter v,11,5,25,30,f
  '
  Inter v,1,2,31,36,f
  Inter v,2,3,32,37,f
  Inter v,3,4,33,38,f
  Inter v,4,5,34,39,f
  Inter v,5,1,35,40,f
  '
  'MIDDLE
  '
  Inter v,1,6,41,46,f
  Inter v,2,7,42,47,f
  Inter v,3,8,43,48,f
  Inter v,4,9,44,49,f
  Inter v,5,10,45,50,f
  '
  Inter v,2,6,51,56,f
  Inter v,3,7,52,57,f
  Inter v,4,8,53,58,f
  Inter v,5,9,54,59,f
  Inter v,1,10,55,60,f
  '
  Inter v,6,7,61,66,f
  Inter v,7,8,62,67,f
  Inter v,8,9,63,68,f
  Inter v,9,10,64,69,f
  Inter v,10,6,65,70,f
  '
  'BOTTOM
  Inter v,12,6,76,71,f
  Inter v,12,7,77,72,f
  Inter v,12,8,78,73,f
  Inter v,12,9,79,74,f
  Inter v,12,10,80,75,f

  Pentagonn v[21],v[22],v[23],v[24],v[25] 'P11 TOP
  Pentagonn v[76],v[77],v[78],v[79],v[80] 'P12 BOTTOM
  '
  Pentagonn v[26],v[40],v[55],v[41],v[31] 'P1
  Pentagonn v[27],v[36],v[56],v[42],v[32] 'P2
  Pentagonn v[28],v[37],v[57],v[43],v[33] 'P3
  Pentagonn v[29],v[38],v[58],v[44],v[34] 'P4
  Pentagonn v[30],v[39],v[59],v[45],v[35] 'P5
  '
  Pentagonn v[46],v[70],v[71],v[61],v[56] 'P6
  Pentagonn v[47],v[66],v[72],v[62],v[57] 'P7
  Pentagonn v[48],v[67],v[73],v[63],v[58] 'P8
  Pentagonn v[49],v[68],v[74],v[64],v[59] 'P9
  Pentagonn v[50],v[69],v[75],v[65],v[60] 'P10
  '
  Hexagonn v[21],v[26],v[31],v[36],v[27],v[22] 'H1
  Hexagonn v[22],v[27],v[32],v[37],v[28],v[23] 'H2
  Hexagonn v[23],v[28],v[33],v[38],v[29],v[24] 'H3
  Hexagonn v[24],v[29],v[34],v[39],v[30],v[25] 'H4
  Hexagonn v[25],v[30],v[35],v[40],v[26],v[21] 'H5
  '
  Hexagonn v[31],v[41],v[46],v[56],v[51],v[36] 'H6
  Hexagonn v[32],v[42],v[47],v[57],v[52],v[37] 'H7
  Hexagonn v[33],v[43],v[48],v[58],v[53],v[38] 'H8
  Hexagonn v[34],v[44],v[49],v[59],v[54],v[39] 'H9
  Hexagonn v[35],v[45],v[50],v[60],v[55],v[40] '10
  '
  Hexagonn v[47],v[42],v[51],v[56],v[61],v[66] 'H11
  Hexagonn v[48],v[43],v[52],v[57],v[62],v[67] 'H12
  Hexagonn v[49],v[44],v[53],v[58],v[63],v[68] 'H13
  Hexagonn v[50],v[45],v[54],v[59],v[64],v[69] 'H14
  Hexagonn v[46],v[41],v[55],v[60],v[65],v[70] 'H15
  '
  Hexagonn v[66],v[61],v[71],v[76],v[77],v[72] 'H16
  Hexagonn v[67],v[62],v[72],v[77],v[78],v[73] 'H17
  Hexagonn v[68],v[63],v[73],v[78],v[79],v[74] 'H18
  Hexagonn v[69],v[64],v[74],v[79],v[80],v[75] 'H19
  Hexagonn v[70],v[65],v[75],v[80],v[76],v[71] 'H20

  end procedure



  procedure Render60(int n=0)
  ===========================
  'indexbase 1
  vector v[40]
  Icosavecs v      '1 to 12
  'rn=0.886838
  DoDecVecs v[21]  '21 to 40
  'rn=1.0
  '
  'TOP
  RenderTriangle v,11,21,22
  RenderTriangle v,11,22,23
  RenderTriangle v,11,23,24
  RenderTriangle v,11,24,25
  RenderTriangle v,11,25,21
  '
  RenderTriangle v,2,22,21
  RenderTriangle v,3,23,22
  RenderTriangle v,4,24,23
  RenderTriangle v,5,25,24
  RenderTriangle v,1,21,25
  '
  RenderTriangle v,1,26,21
  RenderTriangle v,2,27,22
  RenderTriangle v,3,28,23
  RenderTriangle v,4,29,24
  RenderTriangle v,5,30,25
  '
  RenderTriangle v,2,21,26
  RenderTriangle v,3,22,27
  RenderTriangle v,4,23,28
  RenderTriangle v,5,24,29
  RenderTriangle v,1,25,30
  '
  RenderTriangle v,2,26,31
  RenderTriangle v,3,27,32
  RenderTriangle v,4,28,33
  RenderTriangle v,5,29,34
  RenderTriangle v,1,30,35
  '
  RenderTriangle v,6,31,26
  RenderTriangle v,7,32,27
  RenderTriangle v,8,33,28
  RenderTriangle v,9,34,29
  RenderTriangle v,10,35,30
  '
  RenderTriangle v,1,35,26
  RenderTriangle v,2,31,27
  RenderTriangle v,3,32,28
  RenderTriangle v,4,33,29
  RenderTriangle v,5,34,30
  '
  RenderTriangle v,6,26,35
  RenderTriangle v,7,27,31
  RenderTriangle v,8,28,32
  RenderTriangle v,9,29,33
  RenderTriangle v,10,30,34
  '
  RenderTriangle v,6,36,31
  RenderTriangle v,7,37,32
  RenderTriangle v,8,38,33
  RenderTriangle v,9,39,34
  RenderTriangle v,10,40,35
  '
  RenderTriangle v,7,31,36
  RenderTriangle v,8,32,37
  RenderTriangle v,9,33,38
  RenderTriangle v,10,34,39
  RenderTriangle v,6,35,40
  '
  RenderTriangle v,7,36,37
  RenderTriangle v,8,37,38
  RenderTriangle v,9,38,39
  RenderTriangle v,10,39,40
  RenderTriangle v,6,40,36
  '
  RenderTriangle v,12,37,36
  RenderTriangle v,12,38,37
  RenderTriangle v,12,39,38
  RenderTriangle v,12,40,39
  RenderTriangle v,12,36,40

  end procedure








  procedure VectorSnub(vector *v1,*v2,*v3, single*sa)
  ===================================================
  double f,h,ds,ha,s1,w
  '
  'DISTANCE BETWEEN POINTS
  '
  v3.x=v2.x-v1.x
  v3.y=v2.y-v1.y
  v3.z=v2.z-v1.z
  ds=sqr( v3.x*v3.x + v3.y*v3.y + v3.z*v3.z )
  '
  ha=asin(ds*.5)   'HALF ANGLE BETWEEN POINTS
  h=1-cos(sa)      'DEPTH OF SNUB
  'w=h/tan(ha)     'WIDTH OF SNUB SEGMENT
  s1=h/sin(ha)     'DISTANCE OF SNUB POINT FROM V1
  f=s1/ds          'PROPORTION OF DS
  '
  'LOCATION OF SNUB POINT
  '
  v3.x=v1.x+v3.x*f
  v3.y=v1.y+v3.y*f
  v3.z=v1.z+v3.z*f
  '
  end procedure


  procedure RenderSnubTri(vector *v1,*v2,*v3, single sa)
  ======================================================
  '
  vector va,vb,vc,vn
  '
  def plot
  glVertex3f va.x, va.y, va.z
  end def

  '
  vc.x=0.33333*(v1.x+v2.x+v3.x)
  vc.y=0.33333*(v1.y+v2.y+v3.y)
  vc.z=0.33333*(v1.z+v2.z+v3.z)
  '
  glBegin GL_TRIANGLE_FAN
  '
  glNormal3f vc.x, vc.y, vc.z
  glVertex3f vc.x, vc.y, vc.z
  '
  VectorSnub v1,v2,va,sa
  plot
  VectorSnub v2,v1,va, sa
  plot
  VectorSnub v2,v3,va, sa
  plot
  VectorSnub v3,v2,va, sa
  plot
  VectorSnub v3,v1,va, sa
  plot
  VectorSnub v1,v3,va, sa
  plot
  VectorSnub v1,v2,va, sa
  plot
  glEnd
  end procedure


  procedure RenderSnubTriP(vector *v1,*v2,*v3, single sa)
  =======================================================
  '
  vector va,vb,vc,vn
  '
  def plot
  SurfaceNormal vn,va,vb,vc
  glNormal3f vn.x, vn.y, vn.z
  glVertex3f va.x, va.y, va.z
  glVertex3f vb.x, vb.y, vb.z
  glVertex3f vc.x, vc.y, vc.z
  va<=vb.x,vb.y,vb.z
  end def

  '
  vc.x=0.333333*(v1.x+v2.x+v3.x)
  vc.y=0.333333*(v1.y+v2.y+v3.y)
  vc.z=0.333333*(v1.z+v2.z+v3.z)
  '
  glBegin GL_TRIANGLES
  '
  VectorSnub v1,v2,va,sa
  VectorSnub v2,v1,vb,sa
  plot
  VectorSnub v2,v3,vb, sa
  plot
  VectorSnub v3,v2,vb, sa
  plot
  VectorSnub v3,v1,vb, sa
  plot
  VectorSnub v1,v3,vb, sa
  plot
  VectorSnub v1,v2,vb, sa
  plot
  glEnd
  end procedure


  type GeodesicPanelSpec
  =====================
  int p    'recursive pattern
  int q    'panel form
  int n    'recursion level
  int i    'serial index
  sys d    'pointer to vector array
  int c    'vector count
  end type



  procedure RecursTr4(vector *v1,*v2,*v3, sys n, GeodesicPanelSpec*g)
  ===================================================================
  '
  vector v4,v5,v6
  single sa
  '
  if n<1 then
    '
    'generate triangle
    '
    sa=.15 'radians
    select case g.q 'panel form
    case 1 : RenderTriangle v1,v2,v3
    case 2 : RenderSnubTri v1,v2,v3,sa
    case 3 : RenderSnubTriP v1,v2,v3,sa
    end select
    g.i++
    exit procedure
  end if
  n--
  '
  VectorBetween v1,v2,v4  'mid left
  VectorBetween v2,v3,v5  'mid bottom
  VectorBetween v3,v1,v6  'mid right
  '
  recursTr4 v1,v4,v6,n,g 'triangle top
  recursTr4 v4,v2,v5,n,g 'triangle lower left
  recursTr4 v6,v5,v3,n,g 'triangle lower right
  recursTr4 v5,v6,v4,n,g 'triangle centre inverted
  '
  end procedure



  procedure GeoRender(int a,b,c,GeodesicPanelSpec*g)
  ==================================================
  indexbase 1
  vector*v : @v=g.d
  sys n=g.n
  select case g.p
  case 4 : RecursTr4 v[a],v[b],v[c],n,g
  end select
  end procedure


  procedure central(sys pv,ia,ib,ic,id)
  =====================================
  static double ft=atn(2)                       'arc of isos triangle side
  static double fc=atn(tan(pi()*.2)*sin(ft*.5)) '20.9...
  static double fs=cos(fc)/cos(ft*.25)
  static double th=1/3
  vector *a,*b,*c,*d
  %sv sizeof vector 'internal itr metaval sizeof
  '
  @a=pv+ia*sv-sv
  @b=pv+ib*sv-sv
  @c=pv+ic*sv-sv
  @d=pv+id*sv-sv
  '
  d.x=th*(a.x+b.x+c.x)
  d.y=th*(a.y+b.y+c.y)
  d.z=th*(a.z+b.z+c.z)
  'radius adjustment
  double r=sqr(d.x*d.x + d.y*d.y + d.z*d.z)
  double f=fs/r
  d.x *= f
  d.y *= f
  d.z *= f
  end procedure
  '
  procedure FiboDiamondGeo(vector *v)
  ===================================
  'requires min 32 vector array
  sys p=@v
  '
  IcosaVecs v 'returns the 12 icosahedral points (vectors)
  '
  'APPEND (ICOS FACE CENTRE) POINTS FOR DIAMOND
  '
  'TOP
  Central p,11,1,2,13
  Central p,11,2,3,14
  Central p,11,3,4,15
  Central p,11,4,5,16
  Central p,11,5,1,17
  'UPPER MIDDLE
  Central p,1,6,2,18
  Central p,2,7,3,19
  Central p,3,8,4,20
  Central p,4,9,5,21
  Central p,5,10,1,22
  'LOWER MIDDLE
  Central p,6,7,2,23
  Central p,7,8,3,24
  Central p,8,9,4,25
  Central p,9,10,5,26
  Central p,10,6,1,27
  'BOTTOM
  Central p,12,6,10,28
  Central p,12,10,9,29
  Central p,12,9,8,30
  Central p,12,8,7,31
  Central p,12,7,6,32
  end procedure


  procedure RenderIcosGeo(int n)
  ==============================
  'requires min 12 vector array
  indexbase 1
  vector v[12]
  Icosavecs v
  GeodesicPanelSpec g
  '
  g.p=4    'recursive pattern
  g.q=3    'panel form '3
  g.n=n    'recursion level
  g.i=0    'serial index
  g.d=@v   'pointer to vector array
  g.c=20   'vector count
  '
  'TOP
  GeoRender 11,1,2,g
  GeoRender 11,2,3,g
  GeoRender 11,3,4,g
  GeoRender 11,4,5,g
  GeoRender 11,5,1,g
  'UPPER MIDDLE
  GeoRender 1,6,2,g
  GeoRender 2,7,3,g
  GeoRender 3,8,4,g
  GeoRender 4,9,5,g
  GeoRender 5,10,1,g
  'LOWER MIDDLE
  GeoRender 6,7,2,g
  GeoRender 7,8,3,g
  GeoRender 8,9,4,g
  GeoRender 9,10,5,g
  GeoRender 10,6,1,g
  'BOTTOM
  GeoRender 12,6,10,g
  GeoRender 12,10,9,g
  GeoRender 12,9,8,g
  GeoRender 12,8,7,g
  GeoRender 12,7,6,g
  end procedure


  procedure RenderOctaGeo(int n)
  ==============================
  'requires min 6 vector array
  indexbase 1
  vector v[6]
  Octavecs v
  GeodesicPanelSpec g
  '
  g.p=4    'recursive pattern
  g.q=1    'panel form '3
  g.n=n    'recursion level
  g.i=0    'serial index
  g.d=@v   'pointer to vector array
  g.c=20   'vector count
  '
  'TOP
  GeoRender 1,2,3,g
  GeoRender 1,3,4,g
  GeoRender 1,4,5,g
  GeoRender 1,5,2,g
  'BOTTOM
  GeoRender 6,3,2,g
  GeoRender 6,4,3,g
  GeoRender 6,5,4,g
  GeoRender 6,2,5,g
  end procedure


  procedure RenderTetraGeo(int n)
  ===============================
  'requires min 4 vector array
  indexbase 1
  vector v[4]
  Tetravecs v
  GeodesicPanelSpec g
  '
  g.p=4    'recursive pattern
  g.q=1    'panel form '3
  g.n=n    'recursion level
  g.i=0    'serial index
  g.d=@v   'pointer to vector array
  g.c=20   'vector count
  '
  'TOP
  GeoRender 1,2,3,g
  GeoRender 1,3,4,g
  GeoRender 1,4,2,g
  'BOTTOM
  GeoRender 4,3,2,g
  end procedure


