' h1alib.bas - jackord@kw.igs.net - revised 10 Sep 02 - Liberty Basic v3.01
' path for fastest route between fixed points in a uniform gravitational field

' Initialize Window
  nomainwin
  button#1, "Vary Path", [vary], UL, 108, 10, 80, 20
  WindowWidth=410     ' pixel scale 0-400
  WindowHeight=229    ' pixel scale 0-200
  UpperLeftX=100: UpperLeftY=100
  open "Path Variation" for graphics_nsb as #1
  #1 "trapclose [quit]"
  dim y(257)
  y(0)=0: y(256)=5: gr=9.8
  #1 "cls ; down ; color red"                    ' Endpoints
  #1 "place 20 40 ; circle 5 ; place 276 90 ; circle 5"
  #1 "color blue; goto 20 40"                    ' Straight line path
  j=0: k=256: gosub [dt]
  #1 "place 300 35": #1 "\T(0) = "; left$(str$(d), 5)
  pi=4*atn(1)

[waitHere]
  wait

[vary]                                                 ' Vary path
  a=pi: fa=(1-cos(a))/(a-sin(a))-50/256                ' Find cycloid angle
  b=2*pi: fb=(1-cos(b))/(b-sin(b))-50/256
  while (b-a)>.00001                                   ' Bisection algorithm
    c=(b+a)/2: fc=(1-cos(c))/(c-sin(c))-50/256
    if fc*fa>0 then
      a=c: fa=fc
    else
      b=c: fb=fc
    end if
  wend
  c=(a+b)/2: r=y(256)/(1-cos(c)): w=(gr/r)^.5: t=c/w   ' Cycloid angle, r, t
  n=200: #1 "place 20 40 ; color green"          ' Draw cycloid
  for i=1 to n
    aa=i*c/n
    x1=20+10*(r*aa-r*sin(aa)): y1=40+10*r*(1-cos(aa))
    #1 "goto "; x1; " "; y1
  next i
  #1 "place 300 20": #1 "\T(A) = "; left$(str$(t), 5)
  #1 "color blue"
  dn=256
  for n=1 to 8                                         ' 2^n line segments
    dn=dn/2
    for ii=1 to 16                                     ' Iteration loop
      for i=dn to 256-dn step 2*dn                     ' Midpoints of
        gosub [yv]                                     '   previous segments
      next i
      if n>1 then
        for i=2*dn to 256-2*dn step 2*dn               ' Boundaries of
          gosub [yv]                                   '   previous segments
        next i
      end if
      scan
    next ii
    if n=8 then #1 "color red"                   ' Final path in red
    #1 "place 20 40"
    tt=0
    for i=dn to 256 step dn                            ' Path and travel time
      #1 "goto "; 20+i; " "; 40+int(10*y(i))
      j=i-dn: k=i: gosub [dt]: tt=tt+d
    next i
    #1 "place 300 "; 35+15*n
    #1 "\T("; n; ") = "; left$(str$(tt), 5)
  next n
goto [waitHere]

[yv]                                                   ' Vary mid-point y
  y(i)=(y(i-dn)+y(i+dn))/2                             '   to minimize t
  j=i-dn: k=i: gosub [dt]: t1=d
  j=i+dn: gosub [dt]: t1=t1+d
  dd=1
  while abs(dd)>.01                                    ' Step back and forth
    y(i)=y(i)+dd                                       '   across minimum
    j=i-dn: k=i: gosub [dt]: t=d
    j=i+dn: gosub [dt]: t=t+d
    if t>t1 then dd=0-dd/2
    t1=t
  wend
  y(i)=y(i)+2*dd
return

[dt]                                                   ' Time to traverse
  if j=0 then u=0 else u=(2*gr*y(j))^.5                '   single segment
  v=(2*gr*y(k))^.5
  dx=(k-j)/10: dy=y(k)-y(j)
  ds=(dx*dx+dy*dy)^.5
  d=2*ds/(u+v)
return

[quit]
  close #1
end