  '---------------+++-----------------------------
  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


  '-------------------
  sub 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 sub


  '--------------
  sub 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 sub



  '-------------------------------
  sub 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 sub


  '-----------------
  sub Icos(vector*v)
  '=================
  '
  'vector    v[12]
  single x1,x2,x3,y1,y2,y3,z1,z2,z3
  single a,b,r,a1,a3,r1,ha
  sys 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
  '
  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, 1,0
  v[12]<=0,-1,0
  '
  end sub



  '-------------------------------------
  sub 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 sub


  '--------------------------------------------
  sub 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 sub


  '-----------------------------------------------
  sub 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 sub


  '------------------------------------------------
  sub RenderSnubTriP(vector *v1,*v2,*v3, single sa)
  '================================================
  '
  vector va,vb,vc,vn
  '
  def plot
  normal va,vb,vc,vn
  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 sub


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



  '------------------------------------------------------------
  sub 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 sub
  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 sub



  '-------------------------------------------
  Sub GeoRender(sys a,b,c,GeodesicPanelSpec*g)
  '===========================================
  indexbase 1
  vector*v : @v=g.d
  'if a>g.c or b>g.c or c>g.c then FlagError : exit sub
  sys n=g.n
  select case g.p
  case 4 : RecursTr4 v[a],v[b],v[c],n,g
  end select
  end sub

  sub 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 sub
  '
  '----------------------------
  sub FiboDiamondGeo(vector *v)
  '============================
  'requires min 32 vector array
  sys p=@v
  '
  Icos 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 sub


  '-----------------------
  sub RenderIcosGeo(sys n)
  '=======================
  'requires min 12 vector array
  indexbase 1
  vector v[12]
  Icos v
  GeodesicPanelSpec g
  '
  g.p=4    'recursive pattern
  g.q=3    'panel form
  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 sub

