#
# Maple routines to create animated orbits
#
#   Written by Richard Bernatz     September 28, 2000
#
orbits2d1:=proc(x1,y1,dt,n_steps,hmin,hmax,vmin,vmax)
local rkpts1,rklist1,plist1,ptsplot1,trajplt1,i;
rkpts1:=rungekuttahf([eq1,eq2],[0,x1,y1],dt,n_steps):
rklist1:=makelist(rkpts1,2,3):
plist1:=[]:
  for i from 1 to n_steps+1 do
  plist1:=[op(plist1),[rklist1[i][1],rklist1[i][2]]]:
  ptsplot1:=plot({[plist1[i][1],plist1[i][2]]},style=point,symbol=circle,
                   color=blue,view=[hmin..hmax,vmin..vmax]):
    trajplt1:=plot(plist1,color=blue);
    ptstraj||i:=display([ptsplot1,trajplt1]):
    od: 
 display([ptstraj||(1..n_steps+1)],view=[hmin..hmax,vmin..vmax],insequence=true);
 end:
# 
# Animate the phase portrait showing the path. Two paths.
# 
orbits2d2:=proc(x1,y1,x2,y2,dt,n_steps,hmin,hmax,vmin,vmax)
local rkpts1, rkpts2,rklist1,rklist2,plist1,plist2,
          ptsplot1,ptsplot2,trajplt1,trajplt2,i;
 rkpts1:=rungekuttahf([eq1,eq2],[0,x1,y1],dt,n_steps):
 rklist1:=makelist(rkpts1,2,3):
 rkpts2:=rungekuttahf([eq1,eq2],[0,x2,y2],dt,n_steps):
 rklist2:=makelist(rkpts2,2,3):
 plist1:=[]: plist2:=[]:
 for i from 1 to n_steps+1 do
    plist1:=[op(plist1),[rklist1[i][1],rklist1[i][2]]]:
   
    ptsplot1:=plot({[plist1[i][1],plist1[i][2]]},style=point,symbol=circle,
                   color=blue,view=[hmin..hmax,vmin..vmax]):
    trajplt1:=plot(plist1,color=blue);
    plist2:=[op(plist2),[rklist2[i][1],rklist2[i][2]]]:
    ptsplot2:=plot({[plist2[i][1],plist2[i][2]]},style=point,symbol=circle,
                   color=red):
    trajplt2:=plot(plist2,color=red);
    ptstraj||i:=display([ptsplot1,trajplt1,ptsplot2,trajplt2],view=[hmin..hmax,vmin..vmax]):
    od: 
 display([ptstraj||(1..n_steps+1)],view=[hmin..hmax,vmin..vmax],insequence=true);
 end:
#
# Animate a single path of a three equation system. The call to the routine must specify the plane
#  of the plot through the choice of variables
# 
orbits3x1:=proc(x1,y1,z1,dt,n_steps,var1,var2)
  local rkpts1,rklist1,plist1,ptsplot1,trajplt1,i;
  rkpts1:=rungekuttahf([eq1,eq2,eq3],[0,x1,y1,z1],dt,n_steps):
  rklist1:=makelist(rkpts1,var1,var2):
  plist1:=[]:
  for i from 1 to n_steps+1 do
     plist1:=[op(plist1),[rklist1[i][1],rklist1[i][2]]]:
     ptsplot1:=plot({[plist1[i][1],plist1[i][2]]},style=point,symbol=circle,
                    color=blue):
     trajplt1:=plot(plist1,color=blue);
     ptstraj||i:=display([ptsplot1,trajplt1]):
     od: 
  display([ptstraj||(1..n_steps+1)],insequence=true);
  end:
#
#Animate trajectories. Three equation system and two paths.
#
orbits3x2:=proc(x1,y1,z1,x2,y2,z2,dt,n_steps,var1,var2)
  local rkpts1,rkpts2,rklist1,rklist2,plist1,plist2, 
          ptsplot1,ptsplot2,trajplt1,trajplt2,i;
  rkpts1:=rungekuttahf([eq1,eq2,eq3],[0,x1,y1,z1],dt,n_steps):
  rklist1:=makelist(rkpts1,var1,var2):
  rkpts2:=rungekuttahf([eq1,eq2,eq3],[0,x2,y2,z2],dt,n_steps):
  rklist2:=makelist(rkpts2,var1,var2):
  plist1:=[]: plist2:=[]:
  for i from 1 to n_steps+1 do
     plist1:=[op(plist1),[rklist1[i][1],rklist1[i][2]]]:
     ptsplot1:=plot({[plist1[i][1],plist1[i][2]]},style=point,symbol=circle,
                    color=blue):
     trajplt1:=plot(plist1,color=blue);
     plist2:=[op(plist2),[rklist2[i][1],rklist2[i][2]]]:
     ptsplot2:=plot({[plist2[i][1],plist2[i][2]]},style=point,symbol=circle,
                    color=red):
     trajplt2:=plot(plist2,color=red);
     ptstraj||i:=display([ptsplot1,trajplt1,ptsplot2,trajplt2]):
     od: 
  display([ptstraj||(1..n_steps+1)],insequence=true);
  end:
#
# A single path in 3D
#
orbits3D1:=proc(x1,y1,z1,dt,n_steps)
  local rkpts1,rklist1,whole_list1,plist3D1, 
          ptsplot1,trajplt1,i;
  rkpts1:=rungekuttahf([eq1,eq2,eq3],[0,x1,y1,z1],dt,n_steps):
  whole_list1:=[seq([rkpts1[i][2],rkpts1[i][3],rkpts1[i][4]],i=0..n_steps)]:
  plist3D1:=[]: 
  for i from 1 to n_steps+1 do
     plist3D1:=[op(plist3D1),[whole_list1[i][1],whole_list1[i][2],whole_list1[i][3]]]:
#     ptsplot1:=plot3D({[plist3D1[i][1],plist3D1[i][2],plist3D1[i][3]]},style=point,
#               symbol=circle,color=blue):
     trajplt1:=display(curve(plist3D1),color=blue,axes=normal,labels=[x,y,z]);
     ptstraj||i:=display([trajplt1]):
     od: 
  display([ptstraj||(1..n_steps+1)],insequence=true);
  end:
#
# Two paths in 3D
#
orbits3D2:=proc(x1,y1,z1,x2,y2,z2,dt,n_steps,rkpts1,rkpts2)
  local whole_list1,whole_list2,plist3D1,plist3D2, 
          ptsplot1,ptsplot2,trajplt1,trajplt2,i;
  rkpts1:=rungekuttahf([eq1,eq2,eq3],[0,x1,y1,z1],dt,n_steps):
  rkpts2:=rungekuttahf([eq1,eq2,eq3],[0,x2,y2,z2],dt,n_steps):
  whole_list1:=[seq([rkpts1[i][2],rkpts1[i][3],rkpts1[i][4]],i=0..n_steps)]:
  whole_list2:=[seq([rkpts2[i][2],rkpts2[i][3],rkpts2[i][4]],i=0..n_steps)]:
  plist3D1:=[]: plist3D2:=[]:
  for i from 1 to n_steps+1 do
     plist3D1:=[op(plist3D1),[whole_list1[i][1],whole_list1[i][2],whole_list1[i][3]]]:
     plist3D2:=[op(plist3D2),[whole_list2[i][1],whole_list2[i][2],whole_list2[i][3]]]:
#     ptsplot1:=plot3D({[plist3D1[i][1],plist3D1[i][2],plist3D1[i][3]]},style=point,
#               symbol=circle,color=blue):
     trajplt1:=display(curve(plist3D1),color=blue,axes=normal,labels=[x,y,z]);
     trajplt2:=display(curve(plist3D2),color=red,axes=normal,labels=[x,y,z]);
     ptstraj||i:=display([trajplt1,trajplt2]):
     od: 
  display([ptstraj||(1..n_steps+1)],insequence=true);
  end: