
  '===================================
  'SOLVING SPHERICAL GEOMETRY PROBLEMS
  'WITH SUCCESSIVE APPROXIMATION
  '===================================

  '
  'global variables for transient use
  '
  double a,b,c,d,e,h
  '
  '---------------------------------------------------------
  function SplitTriangle( ans as double ) as double callback
  '=========================================================
  '
  '
  'SPHERICAL TRIANGLE
  '
  '        ^
  '       /|\
  '      / | \
  '    a/  |  \b
  '    /   |   \
  '   /    |h?  \
  '  /     |     \
  '  - -----------
  '     d?    e?
  '
  ' c=d+e
  ' measurements in radians
  ' cos(h) = cos(a)/cos(d) = cos(b)/cos(e)
  '
  ' split into 2 spherical right angle triangles
  ' a and b should be less than pi/2
  ' c should be less than a+b
  '
  d=ans
  e=c-d
  'equations:  cos(h) = cos(a)/cos(d) = cos(b)/cos(e)
  '
  return cos(a)/cos(d) - cos(b)/cos(e) 'null expression returning error
  '
  end function

  declare function ptr feedback(ans as double) as double
  '
  '-----------------------------------------------------
  function approx(sys pfunction, double a1,a2) as double
  '=====================================================
  '
  if not pfunction then print "invalid feedback function" : return 0
  '
  @feedback=pfunction
  '
  double i1,i2,scale
  sys c
  '
  'a1             'first estimate of answer
  'a2             'second estimate of answer
  i1=feedback(a1) 'error from first estimate
  i2=feedback(a2) 'error from second estimate
  '
  do
    c++
    if (abs(a2-a1)<1e-16) then exit do 'successfully reached precision limit
    if c>1000 then print "Unable to resolve" : exit do 'iteration limit
    '
    if i1=i2 then exit do 'feedback scale limit reached
    scale=(a2-a1)/(i2-i1)
    a1=a2
    a2-=i2*scale
    i1=i2
    i2=feedback(a2) 
  end do
  '
  'print str(c) " iterations"  '
  return a2
  '
  end function

  '
  'MAIN
  '====
  '
  'set params
  '
  a=atn 2 'arc of icosahedral triangle (well known value to test)
  b=a
  c=a
  '
  'calc
  '
  d=approx (@SplitTriangle,1,2)
  h=acos(cos(a)/cos(d))
  '
  'display
  '-------
  '
  string cr=chr(13,10)
  print "Split Triangle Results:" cr +
  "d=" cr str(d) " radians" cr+
  str(deg(d)) " degrees" cr +
  "a/2=" cr str(a/2) " radians" cr+
  str(deg(a/2)) " degrees" cr +
  "h=" cr str(h) " radians" cr+
  str(deg(h)) " degrees" cr 

  