  'GEOMETRY.INC  

  type vector3f single x,y,z
  type vector3d  double x,y,z

  'adopt single precision
  '======================

  typedef vector3f     vector


  type triangle vector p1,p2,p3


  sub VectorScale(vector *v1,*v2)
  ===============================
  single r
  r=1/sqr(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z)
  v2.x*=r
  v2.y*=r
  v2.z*=r
  end sub

  sub VectorBetween(vector *v1,*v2,*v3)
  =====================================
  v3.x=0.5*(v1.x+v2.x)
  v3.y=0.5*(v1.y+v2.y)
  v3.z=0.5*(v1.z+v2.z)
  vectorscale v3,v3
  end sub


  sub VectorAlong(vector *v1,*v2,*v3, single*f)
  =============================================
  single vd
  vd=v2.x-v1.x : v3.x=v1.x+vd*f
  vd=v2.y-v1.y : v3.y=v1.y+vd*f
  vd=v2.z-v1.z : v3.z=v1.z+vd*f
  vectorscale v3,v3
  end sub


  sub TriCentre(vector *v1,*v2,*v3,*v4)
  =====================================
  v4.x=0.5*(v1.x+v2.x+v3.x)
  v4.y=0.5*(v1.y+v2.y+v3.y)
  v4.z=0.5*(v1.z+v2.z+v3.z)
  vectorscale v4,v4
  end sub



  'http://www.mvps.org/directx/articles/math/dot/index.htm
  'http://en.wikipedia.org/wiki/Dot_product


  function DotProduct(vector *u,*v) as single
  ===========================================
  return u.x*v.x + u.y*v.y + u.z*v.z
  end function



/*
  Surface Normals
  ===============

  These are vectors perpendicular to a surface.
  They are mostly used to calculate light reflected from 3d surfaces.

  http://www.opengl.org/wiki/Calculating_a_Surface_Normal
  discussion on use of normals:
  http://www.opengl.org/resources/features/KilgardTechniques/oglpitfall/


  Pseudocode
  ----------

  Begin Function CalculateSurfaceNormal (Input Triangle) Returns Vector

  Set Vector U to (Triangle.p2 minus Triangle.p1)
  Set Vector V to (Triangle.p3 minus Triangle.p1)

  Set Normal.x to (multiply U.y by V.z) minus (multiply U.z by V.y)
  Set Normal.y to (multiply U.z by V.x) minus (multiply U.x by V.z)
  Set Normal.z to (multiply U.x by V.y) minus (multiply U.y by V.x)

  Returning Normal

  End Function
*/


  sub normal(vector *p1,*p2,*p3,*n)
  =================================

  'assumes points are defined anticlockwise

  vector3f u,v

  'set Vector U to (Triangle.p2 minus Triangle.p1)
  'set Vector V to (Triangle.p3 minus Triangle.p1)

  u.x=p2.x-p1.x
  u.y=p2.y-p1.y
  u.z=p2.z-p1.z
  v.x=p3.x-p1.x
  v.y=p3.y-p1.y
  v.z=p3.z-p1.z


  'Set Normal.x to (multiply U.y by V.z) minus (multiply U.z by V.y)
  'Set Normal.y to (multiply U.z by V.x) minus (multiply U.x by V.z)
  'Set Normal.z to (multiply U.x by V.y) minus (multiply U.y by V.x)

  n.x=(u.y*v.z)-(u.z*v.y)
  n.y=(u.z*v.x)-(u.x*v.z)
  n.z=(u.x*v.y)-(u.y*v.x)

  'then scale the normal to unity

  'VectorScale n,n

  end sub

/*

  sub recurst4(vector *v1,*v2,*v3, sys n, triangle*t,sys*i)
  '========================================================
    '
    vector v4,v5,v6
    '
    if n<1 then
      '
      'generate triangle
      '
      t[i]<=v1.x,v1.y,v1.z,v2.x,v2.y,v2.z,v3.x,v3.y,v3.z
      i+=1
      exit sub
    end if
    n-=1
    '
    VectorBetween v1,v2,v4  'mid left
    VectorBetween v2,v3,v5  'mid bottom
    VectorBetween v3,v1,v6  'mid right
    '
    recurst4 v1,v4,v6,n,t,i 'triangle top
    recurst4 v4,v2,v5,n,t,i 'triangle lower left
    recurst4 v6,v5,v3,n,t,i 'triangle lower right
    recurst4 v5,v6,v4,n,t,i 'triangle centre inverted
    '
  end sub


  sub recurst9(vector *v1,*v2,*v3, sys n, triangle*t,sys*i)
  =========================================================
  '
  vector v4,v5,v6,v7,v8,v9,v10
  double p=1/3
  '
  if n<1 then
    '
    'generate triangle
    '
    t[i]<=v1.x,v1.y,v1.z,v2.x,v2.y,v2.z,v3.x,v3.y,v3.z
    i++
    exit sub
  end if
  n--
  '
  VectorAlong v1,v2,v4,p 'upper left
  VectorAlong v2,v1,v5,p 'mid left
  VectorAlong v2,v3,v6,p 'left bottom
  VectorAlong v3,v2,v7,p 'right bottom
  VectorAlong v3,v1,v8,p 'lower right
  VectorAlong v1,v3,v9,p 'upper right
  TriCentre v1,v2,v3,v10 'centre
  '
  'recurse for 9 triangles
  '
  recurst9 v1, v4,v9,n,t,i '1 outer top
  recurst9 v2, v6,v5,n,t,i '2 outer lower left
  recurst9 v3, v8,v7,n,t,i '3 outer lower right
                           '
  recurst9 v10,v9,v4,n,t,i '4 inner top
  recurst9 v10,v4,v5,n,t,i '5 inner upper left
  recurst9 v10,v5,v6,n,t,i '6 inner lower left
  recurst9 v10,v6,v7,n,t,i '7 inner bottom
  recurst9 v10,v7,v8,n,t,i '8 inner lower right
  recurst9 v10,v8,v9,n,t,i '9 inner upper right
    '
  end sub

*/

  sub CylEnd(single r, y, s)
  ==========================
  '
  single x,z
  sys i,j,st=s
  '
  glBegin GL_TRIANGLE_FAN 
  glNormal3f  0, y,  0
  glVertex3f  0, y,  0
  for i=0 to 360 step st
    if y<0 then
      j=360-i
    else
      j=i
    end if
    x=r*cos(rad(j))
    z=r*sin(rad(j))
    glVertex3f  x,  y,  z
  next
  glEnd
  '
  end sub


  sub Cylinder(single r1, r2, y, m)
  =================================
  '
  single s,c,r,x1,z1,x2,z2,yn
  sys i,st=m
  '
  'SIDE SURFACE
  '
  glBegin GL_QUAD_STRIP
  for i=0 to 360 step st
    r=rad(i)
    c=cos(r)
    s=sin(r)
    x1=r1*c
    z1=r1*s
    x2=r2*c
    z2=r2*s
    yn=sin()
    glNormal3f  x2,  r1-r2,  z2
    glVertex3f  x2,  y,  z2
    glVertex3f  x1, -y,  z1
  next
  glEnd
  '
  'TOP
  '
  Cylend r2,y,m
  '
  'BOTTOM
  '
  Cylend r1,-y,m
  '
  end sub


  sub HemiSphere(single rw, rh, m)
  ================================
  '
  single s,c,r,r1,r2,x1,x2,y1,y2,z1,z2,yn
  sys e,i,j,st=m
  '
  e=90-st
  for j=0 to e step st
    r=rad(j)
    r1=rw*cos(r)
    y1=rh*sin(r)
    r=rad(j+st)
    r2=rw*cos(r)
    y2=rh*sin(r)
    glBegin GL_QUAD_STRIP
    for i=0 to 360 step st
      r=rad(i)
      c=cos(r)
      s=sin(r)
      x1=r1*c
      z1=r1*s
      x2=r2*c
      z2=r2*s
      yn=sin()
      glNormal3f  x2, y2,  z2
      glVertex3f  x2, y2,  z2
      glNormal3f  x2, y1,  z1
      glVertex3f  x1, y1,  z1
    next
    glEnd
  next
  '
  end sub


  sub Cuboid(single x, y, z)
  ==========================
  '
  glBegin GL_QUADS
  '

  'FRONT
  '
  glNormal3f  0,  0,  z
  glVertex3f -x, -y,  z
  glVertex3f -x,  y,  z
  glVertex3f  x,  y,  z
  glVertex3f  x, -y,  z
  '
  'RIGHT
  '
  glNormal3f  x,  0,  0
  glVertex3f  x, -y,  z
  glVertex3f  x,  y,  z
  glVertex3f  x,  y, -z
  glVertex3f  x, -y, -z
  '
  'BACK
  '
  glNormal3f  0,  0, -z
  glVertex3f  x, -y, -z
  glVertex3f  x,  y, -z
  glVertex3f -x,  y, -z
  glVertex3f -x, -y, -z
  '
  'LEFT
  '
  glNormal3f -x,  0,  0
  glVertex3f -x, -y, -z
  glVertex3f -x,  y, -z
  glVertex3f -x,  y,  z
  glVertex3f -x, -y,  z
  '
  'TOP
  '
  glNormal3f  0,  y,  0
  glVertex3f -x,  y,  z
  glVertex3f -x,  y, -z
  glVertex3f  x,  y, -z
  glVertex3f  x,  y,  z
  '
  'BOTTOM
  '
  glNormal3f  0, -y,  0
  glVertex3f  x, -y,  z
  glVertex3f  x, -y, -z
  glVertex3f -x, -y, -z
  glVertex3f -x, -y,  z
  '
  glend
      
  end sub



  sub Toroid(single major1 , minor1, major2 , minor2, ro, yd)
  ===========================================================

    static as double
    '
    twopi          = pi()*2,
    torus_major_res= 32,
    torus_minor_res= 32
    '
    static sys i,j,k,ie,je
    '
    static single a, b, s, t, x, y, z, nx, ny, nz, gscale,tmc,tmd,tme,major,minor,ja
    '
    'Draw the torus
    '--------------
    '
    ie=torus_major_res-1
    je=ro*torus_major_res
    '
    for i = 0 to ie
      '
      glBegin GL_QUAD_STRIP
      '
      for j = 0 to je
        '
        for k = 1 to 0 step -1
          ja=j/je
          major=major1*(1-ja)+major2*ja
          minor=minor1*(1-ja)+minor2*ja
          '
          s = mod( i+k,torus_major_res + 0.5)
          t = mod(j,torus_minor_res)
          '
          'CALCULATE POINT ON SURFACE
          '--------------------------
          '
          tmd=s*twopi/torus_minor_res
          tme=t*twopi/torus_major_res
          tmc=major+minor * cos tmd
          '
          x = tmc * cos tme
          y = minor * sin tmd
          z = tmc * sin tme
          '
          'CALCULATE SURFACE NORMAL
          '------------------------
          '
          a=major * cos tme
          nx = x - a
          ny = y
          a=major * sin tme
          nz = z - a
          '
          '
          glNormal3f nx, ny, nz
          glVertex3f x, y+yd*ja, z
          '
        next
        '
      next
      '
      glend
      '
    next
    '
  end sub

  sub Torus(single major, minor)
  ==============================
  toroid major,minor,major,minor,1,0
  end sub


  sub polycone(single h, sys n)
  =============================
  single a,a1,a2,x1,x2,y1,y2,z1,z2,na,nx,ny,nz
  '
  a=2*pi()/n
  a1=-a*.5
  x1=sin a1
  z1=cos a1
  if h<=0 then h=.001
  ny=1/h
  '
  int i
  for i=1 to n
    a2=a1+a
    x2=sin a2
    z2=cos a2
    '
    'face normal
    '
    nx=sin na
    nz=cos na
    
    '
    glbegin GL_TRIANGLES
    glNormal3f nx,ny,nz
    glVertex3f x1,0,z1
    glVertex3f x2,0,z2
    glVertex3f 0,h,0
    glend
    '
    a1=a2 : x1=x2 : y1=y2 : z1=z2
    na+=a
  next
  '
  end sub

  sub polysides(single h, sys n)
  ==============================
  single a,a1,a2,x1,x2,y1,y2,z1,z2,na,nx,ny,nz
  '
  a=2*pi()/n
  a1=-a*.5
  x1=sin a1
  z1=cos a1
  if h<=0 then h=.001
  ny=0
  int i
  for i=1 to n
    a2=a1+a
    x2=sin a2
    z2=cos a2
    '
    'face normal
    '
    nx=sin na
    nz=cos na
    '
    glbegin GL_QUADS
    glNormal3f nx,ny,nz
    glVertex3f x1,0,z1
    glVertex3f x2,0,z2
    glVertex3f x2,h,z2
    glVertex3f x1,h,z1
    glend
    '
    a1=a2 : x1=x2 : y1=y2 : z1=z2
    na+=a
  next
  '
  end sub


