!
! motion of the "hyperion" (actually a dumbbell
! theta() = dumbbell angle   omega() = dumbbell angular velocity t() = time
!
! made using kepler.tru
!
program lyapunov
   option nolet
   library "sgfunc.trc"
   library "sglib.trc"
   dim tha(100000),thb(100000), oma(100000),omb(100000), t(100000)
   dim delt(100000),delo(100000)
   call initialize(xc,vxc,yc,vyc,tha,thb,oma,omb,t,dt,mr,tot,st1$,st2$,st3$)
   call calculate(xc,vxc,yc,vyc,tha,thb,oma,omb,t,dt,mr,tot,st1$,st2$,st3$,delt,delo)
   call display(delt,delo,t,dt,mr)
end
! initialize variables
! tha(),thb() = dumbbell angles  oma(),omb() = dumbbell angular velocities
sub initialize(xc,vxc,yc,vyc,tha(),thb(),oma(),omb(),t(),dt,mr,tot,st1$,st2$,st3$)
   input prompt "mass ratio of saturn to sun -> ": mr 
   input prompt "initial x position (in AU) -> ": xc
   input prompt "initial y position (in AU) -> ": yc
   input prompt "initial x velocity (in AU/yr) -> ": vxc
   input prompt "initial y velocity (in AU/yr) -> ": vyc
   input prompt "initial angle for series #1 (in radians) -> ": tha(1)
   input prompt "initial omega for #1 (in radians/s) -> ": oma(1)
   input prompt "initial angle for series #2 (in radians) -> ": thb(1)
   input prompt "initial omega for #2 (in radians/s) -> ": omb(1)
   input prompt "time step (in yr) -> ": dt
   input prompt "total time (in yrs) -> ": tot
   input prompt "output file (all t) -> ": st1$
   input prompt "output file (t for max in del theta) -> ": st2$
   input prompt "output file (t for max in del omega) -> ": st3$
   t(1) = 0
end sub
! use the RK2 method
sub calculate(xc,vxc,yc,vyc,tha(),thb(),oma(),omb(),t(),dt,mr,tot,st1$,st2$,st3$,delt(),delo())
   open #1: name st1$, organization text, create newold
   erase #1
   open #2: name st2$, organization text, create newold
   erase #2
   open #3: name st3$, organization text, create newold
   erase #3
   delt(1) = tha(1)-thb(1)
   x1=delt(1)-int(delt(1)/6.283185307)*6.283185307
   if x1 > 3.141592654 then x1=x1-6.283185307
   if x1 < -3.141592654 then x1=x1+6.283185307
   delt(1) = abs(x1)
   delo(1) = abs(oma(1)-omb(1))
   print #1: t(1),delt(1),delo(1)
   i = 0
   do
      i = i + 1
      t(i+1) = t(i) + dt
      call dv(xc,yc,tha(i),thb(i),dvxc,dvyc,doma,domb,mr)
      oma1 = oma(i)+0.5*dt*doma 
      tha1 = tha(i)+0.5*dt*oma(i)
      omb1 = omb(i)+0.5*dt*domb 
      thb1 = thb(i)+0.5*dt*omb(i)
      xc1 = xc+0.5*dt*vxc
      vxc1 = vxc+0.5*dt*dvxc
      yc1 = yc+0.5*dt*vyc
      vyc1 = vyc+0.5*dt*dvyc
      call dv(xc1,yc1,tha1,thb1,dvxc2,dvyc2,doma2,domb2,mr)
      oma(i+1) = oma(i)+dt*doma2
      tha(i+1) = tha(i)+dt*oma1
      omb(i+1) = omb(i)+dt*domb2
      thb(i+1) = thb(i)+dt*omb1
      xc = xc+dt*vxc1
      vxc = vxc+dt*dvxc2
      yc = yc+dt*vyc1
      vyc = vyc+dt*dvyc2
      delt(i+1) = tha(i+1)-thb(i+1)
      x1=delt(i+1)-int(delt(i+1)/6.283185307)*6.283185307
      if x1 > 3.141592654 then x1=x1-6.283185307
      if x1 < -3.141592654 then x1=x1+6.283185307
      delt(i+1) = abs(x1)
      delo(i+1) = abs(oma(i+1)-omb(i+1))
      print #1: t(i+1),delt(i+1),delo(i+1)
      if i+1 >= 3 then
	if (delt(i) > delt(i-1)) and (delt(i+1) < delt(i)) then
	    print #2: t(i),delt(i)
	end if
	if (delo(i) > delo(i-1)) and (delo(i+1) < delo(i)) then
	    print #3: t(i),delo(i)
	end if
      end if
   loop until (t(i+1) >= tot or i+1 >= 100000)
   close #1
   close #2
   close #3
   mat redim oma(i+1),omb(i+1),tha(i+1),thb(i+1),t(i+1),delt(i+1),delo(i+1)
end sub

sub dv(xc0,yc0,tha0,thb0,dvxc,dvyc,doma,domb,mr)
   r = sqr(xc0^2+yc0^2)
   dvxc = -(4*pi^2*mr*xc0)/r^3
   dvyc = -(4*pi^2*mr*yc0)/r^3
   doma = -(12*pi^2*mr)*(xc0*sin(tha0)-yc0*cos(tha0))*(xc0*cos(tha0)+yc0*sin(tha0))/r^5
   domb = -(12*pi^2*mr)*(xc0*sin(thb0)-yc0*cos(thb0))*(xc0*cos(thb0)+yc0*sin(thb0))/r^5
end sub

! display delt and delo in separate windows
sub display(delt(),delo(),t(),dt,mr)
   set background color "white"
   call setcanvas("white")
   set color "black"
   clear
   open #1: screen 0,.9,.94,1
   set cursor 1,1
   print "mass ratio = "; mr;", time step = ";dt
   print "initial angle difference = "; delt(1)
   close #1
   open #1: screen 0,.9,0,.45    ! bottom window
   call settitle("Angular velocity difference vs. time")
   call sethlabel("Time (yr)")
   call setvlabel("Delta Omega")
   call datagraph(t,delo,1,1,"black")
   close #1
   open #1: screen 0,.9,.49,.94    ! top window
   call settitle("Angle difference vs. time")
   call sethlabel("Time (yr)")
   call setvlabel("Delta Theta")
   call datagraph(t,delt,1,1,"black")
   close #1
   get key z
end sub