! whirl

! X(t) = x(t)i +y(t)j + z(t)k
  def x(t)=t*cos(t)
  def y(t)=t*sin(t)
  def z(t)= t

! V(t) = vx(t)i + vy(t)j + vz(t)k
  def vx(t)=-t*sin(t)+cos(t)
  def vy(t)=t*cos(t)+sin(t)
  def vz(t)=1

! A(t) = ax(t)i + ay(t)j + az(t)k
  def ax(t)= -2*sin(t)-t*cos(t)
  def ay(t)= 2*cos(t)-t*sin(t)
  def az(t)= 0




dim xscrn(3),yscrn(3)  ! unit orthogonal vectors which determine screen plane
dim p(3),q(3),r(3)     ! Scratch vectors.
dim u(3),v(3),o(3)
mat  o=zer             ! Zero vector.

def xs(a,b,c)          ! transforms 3 dim into screen coordinates:  xs,ys
   let ys = a*yscrn[1]+b*yscrn[2]+c*yscrn[3]    ! exploits side effect
   let xs = a*xscrn[1]+b*xscrn[2]+c*xscrn[3]
end def

def norm(a()) = sqr(dot(a,a))

sub eye(a,b,c)     ! Finds xscrn,yscrn orthogonal to (a,b,c) and such
  let xscrn[1]=-b  ! that the z-axis is straight up and down.
  let xscrn[2]= a
  let xscrn[3]= 0  ! This makes z-axis straight up and down.
  let yscrn[1]= -a
  let yscrn[2]= -b
  let yscrn[3]= (a*a + b*b)/c   ! needs c not equal zero
  mat xscrn=(1/norm(xscrn))*xscrn  ! normalize
  mat yscrn=(1/norm(yscrn))*yscrn
end sub

sub screendim(width,a,b,c)    ! (a,b,c) is the center of screen
  let height=(2/3)*width
  let xmin=xs(a,b,c)-(width/2)
  let xmax=xs(a,b,c)+(width/2)
  let ymin=ys-(height/2)
  let ymax=ys+(height/2)
end sub

sub axis(a,b,c)
    plot xs(0,0,0),ys ; xs(0,0,c),ys
    plot xs(0,0,0),ys ; xs(0,b,0),ys
    plot xs(0,0,0),ys ; xs(a,0,0),ys
    plot text, at xs(0,0,c),ys : "z"
    plot text, at xs(0,b,0),ys : "y"
    plot text, at xs(a,0,0),ys : "x"
end sub

sub maxaxis
  let inset= (1/20)*(xmax-xmin + ymax-ymin)
  let c3=ymax/yscrn[3]  - inset
  let b1=xmax/xscrn[2]
  let b2=ymin/yscrn[2]
  if b1<b2 then let b3=b1-inset else let b3=b2-inset
  let a1=xmin/xscrn[1]
  let a2=ymin/yscrn[1]
  if a1<a2 then let a3=a1-inset else let a3=a2-inset
  call axis(a3,b3,c3)
end sub

sub vector(p(),q())
    let s1=6/7              ! These two determines shape
    let s2=1/14             ! of arrow head.
    mat u=q-p
    let epsi=(s2)*norm(u)
    mat u=s1*u
    mat u=p+u               ! u = p + s1( q - p )
    plot xs(p(1),p(2),p(3)),ys       ; xs(q(1),q(2),q(3)),ys;
    plot xs(u(1),u(2),u(3)+epsi),ys  ; xs(u(1),u(2),u(3)-epsi),ys;
    plot xs(q(1),q(2),q(3)),ys
end sub

def random(a,b)=(b-a)*rnd+a
def randomint(a,b)=int(random(a,b+1))


sub init
  close #1
  call eye(5,2,3)
  call screendim(100,0,0,20)
  set mode "egahires"
  open #1 : screen 0,1,0,1
  window #1
  set window xmin,xmax,ymin,ymax
  clear
  box lines  xmin,xmax,ymin,ymax
  call maxaxis
end sub

!!     Main program
do
  call init
  let a = 0
  let b = 45
  let n= 2056
  let dt=(b-a)/n
  for t=a to b step dt
    plot xs(x(t),y(t),z(t)),ys;
  next t
  call init
  let t=a
  do until t>b
    set color "green"
    let maxi=randomint(30,70)
    let maxi=13
    for i=1 to maxi
      plot xs(x(t),y(t),z(t)),ys;
      let t=t+dt
    next i
      plot xs(x(t),y(t),z(t)),ys;
    let p(1)=vx(t)
    let p(2)=vy(t)
    let p(3)=vz(t)
    let q(1)=x(t)
    let q(2)=y(t)
    let q(3)=z(t)
    mat p=q+p
    set color "yellow"
    call vector(q,p)
    let p(1)=ax(t)
    let p(2)=ay(t)
    let p(3)=az(t)
    mat p=q+p
    set color "red"
    call vector(q,p)
    if key input then
       get key key
       get key key
    end if
loop
loop 
end