#spin6.plt
#パラメータなど
set parametric
set hidden3d
set xrange [-7 to 7]
set yrange [-7 to 7]
set zrange [-7 to 7]
set urange [-pi/2 to pi/2] #u,v are prepared for dipicting a ball
set vrange [0 to 2*pi]
m = 20 #end
t = 0 #start
r = 2 #radius
#a ball is created in this data table
#the former three dipict the ball, the latter dipicts a small ball stacking to the ball
set table "object.table"
set isosample 15,15
splot r*cos(u)*cos(v),r*cos(u)*sin(v),r*sin(u)#,r+0.5*cos(u)*cos(v),0.5*cos(u)*sin(v),0.5*sin(u)
unset isosample
unset table
unset parametric
x(t) = cos(t)*6 #center position X
y(t) = sin(t)*4 #center position Y
z(t) = r #center position z
dt = 0.2 #ticks
dx(t) = x(t+dt)-x(t)
dy(t) = y(t+dt)-y(t)
load "spin6sub3.plt" #definition of rotation
load "spin6sub.plt" #loop & plot
unset hidden3d
#spin6sub.plt
#グラフのプロット
#dx = x(t+dt)-x(t)
#dy = y(t+dt)-y(t)
#dl = sqrt(dx**2 + dy**2)
#the axis and angle of rotation
a= -dy(t)
b= dx(t)
c= 0
dl = sqrt(a**2+b**2+c**2)
psi = dl/r
#dipict the graph with the center positioned at (x(t),y(t),z(t))
splot "object.table" using ($1+x(t)):($2+y(t)):($3+z(t)) w l,0#,"object.table" using ($1+x(t)):($2+y(t)):(-$3-z(t)) w l
#create v1,v2,v3
load "spin6sub2.plt"
#save rotated positions
set table "newobject.table"
splot "object.table" using (q1x($1,$2,$3)):(q1y($1,$2,$3)):(q1z($1,$2,$3))
unset table
#next
set table "object.table"
splot "newobject.table"
unset table
t = t + dt
if (t<m) reread#spin6sub2.plt
#正規直交基底を３本選ぶ
nor(aa,bb,cc)=sqrt(aa**2 + bb**2 + cc**2)
#a= 1
#b= 0
#c= 0
if (a**2 + b**2 == 0 ) ph = 0
if (a**2 + b**2 >0 ) ph = -atan2(a,b)
#print ph
norm1 = nor(a,b,c)
v1x=a/norm1
v1y=b/norm1
v1z=c/norm1
#print v1x,v1y,v1z
#print nor(v1x,v1y,v1z)
v2x = -cos(ph)
v2y = -sin(ph)
v2z = 0
#print v2x,v2y,v2z
#print nor(v2x,v2y,v2z)
epx= v1y*v2z - v1z*v2y
epy= v1z*v2x - v1x*v2z
epz= v1x*v2y - v1y*v2x
norm2 = nor(epx,epy,epz)
v3x= epx/norm2
v3y= epy/norm2
v3z= epz/norm2
#print v3x,v3y,v3z
#print nor(v3x,v3y,v3z)
#ip12 = v1x*v2x + v1y*v2y + v1z*v2z
#ip23 = v2x*v3x + v2y*v3y + v2z*v3z
#ip31 = v3x*v1x + v3y*v1y + v3z*v1z
#print ip12,ip23,ip31
#spin6sub3.plt
#回転 点Pから点Qへ
p1x(aa,bb,cc)=aa
p1y(aa,bb,cc)=bb
p1z(aa,bb,cc)=cc
#position in v1,v2,v3 space
p2x(aa,bb,cc) = p1x(aa,bb,cc)*v1x+p1y(aa,bb,cc)*v1y+p1z(aa,bb,cc)*v1z
p2y(aa,bb,cc) = p1x(aa,bb,cc)*v2x+p1y(aa,bb,cc)*v2y+p1z(aa,bb,cc)*v2z
p2z(aa,bb,cc) = p1x(aa,bb,cc)*v3x+p1y(aa,bb,cc)*v3y+p1z(aa,bb,cc)*v3z
#psi-rotation of theta in v1,v2,v3 space
q2x(aa,bb,cc) = p2x(aa,bb,cc)
q2y(aa,bb,cc) = p2y(aa,bb,cc)*cos(psi)-p2z(aa,bb,cc)*sin(psi)
q2z(aa,bb,cc) = p2y(aa,bb,cc)*sin(psi)+p2z(aa,bb,cc)*cos(psi)
#position after roatation in ordinal xyz space
q1x(aa,bb,cc)=q2x(aa,bb,cc)*v1x+q2y(aa,bb,cc)*v2x+q2z(aa,bb,cc)*v3x
q1y(aa,bb,cc)=q2x(aa,bb,cc)*v1y+q2y(aa,bb,cc)*v2y+q2z(aa,bb,cc)*v3y
q1z(aa,bb,cc)=q2x(aa,bb,cc)*v1z+q2y(aa,bb,cc)*v2z+q2z(aa,bb,cc)*v3z
