# Author Michael Monagan, 2001.
# modified by Mohammad ali Ebrahimi, 2005
macro( ARROWHULL1 = `GraphTheory/vectorplot3d/arrowhull1`,
ARROWHULL2 = `GraphTheory/vectorplot3d/arrowhull2`,
ARROWHULL3 = `GraphTheory/vectorplot3d/arrowstem3`,
NORM = `GraphTheory/vectorplot3d/norm`,
FAST = `GraphTheory/vectorplot3d/A.v+b`,
vectorplot3d = `GraphTheory/vectorplot3d` ):
# Temporaries needed for loading
macro( N=`GraphTheory/vectorplot3d/N`,
H=`GraphTheory/vectorplot3d/H`,
A=`GraphTheory/vectorplot3d/A` ):
#-----------------------------------------------------------------------
NORM := proc(v) option inline; (v[1]^2+v[2]^2+v[3]^2)^(1/2) end:
#-----------------------------------------------------------------------
vectorplot3d := proc()
local opts,sx,sy,sz,A,p,w,lw,uw,t,dis,as,u,v,b,M,U;
if nargs=1 or not type(args[2],{list,vector,Vector}) then return vectorplot3d([0,0,0],args) fi;
(p,w) := args[1],args[2];
if not type(p,list) then p := [p[1],p[2],p[3]] fi;
if not type(w,list) then w := [w[1],w[2],w[3]] fi;
(p,w) := evalf(p),evalf(w);
if not type(p,[numeric$3]) then error "vector of three real constants expected" fi;
if not type(w,[numeric$3]) then error "vector of three real constants expected" fi;
lw := NORM(w);
uw := w/lw;
t := p+1.05*w;
opts := [args[3..nargs]];
if hasoption( opts, 'arrowsize', 'dis', 'opts') then dis:=evalf(dis); fi;
if hasoption( opts, 'arrowstyle', 'as', 'opts' ) then else as := 'smooth' fi;
if not member(as,{'onedir','nodir','twodir'}) then error "invalid arrowstyle"; fi;
if w[1]=0.0 then
u := [1.0,0.0,0.0];
if w[2]=0.0 then
v := [0.0,1.0,0.0];
else
v := [0.0,-w[3]/w[2],1.0];
v := v/NORM(v);
fi
else
u := [-w[2]/w[1],1.0,0.0];
u := u/NORM(u);
v := [-w[3]/w[1],0.0,1.0];
v := v - u[1]*v[1]*u;
v := v/NORM(v);
fi;
sx := 0.3*dis;
sy := 0.3*dis;
sz := lw;
A := hfarray(1..3,1..3,[[u[1]*sx,v[1]*sy,uw[1]*sz],
[u[2]*sx,v[2]*sy,uw[2]*sz],
[u[3]*sx,v[3]*sy,uw[3]*sz]]);
b := hfarray(1..3,p);
# Compute v' = A.v+b for each point v in ARROWHULL
M := proc(P,A,T) local n; n := op([2,2,2],P); evalhf(FAST(A,P,n,T)) end;
if as='nodir' then
U := seq(M(u,A,b),u=ARROWHULL1);
elif as='twodir' then
U := seq(M(u,A,b),u=ARROWHULL2);
else
U := seq(M(u,A,b),u=ARROWHULL3);
fi;
return(U);
end:
#-----------------------------------------------------------------------
FAST := proc(A,P,n,T) local i,j,B;
B := hfarray(1..n,1..3);
# Compute B = Transpose (A.P + )
for i to 3 do
for j to n do
B[j,i] := T[i]+A[i,1]*P[1,j]+A[i,2]*P[2,j]+A[i,3]*P[3,j];
od;
od;
B;
end:
#----------------------------------------------------------------------
#no direction
A := evalf(Pi/4): N := 8: H := 1.0: ARROWHULL1 :=
POLYGONS(
# an octagon approximating a disk at the base of the arrow head
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ),
# two perpendicular rectangles approximating a cylinder for the body of the arrow
[<-0.04,0.0,0.0>,<-0.04,0.0,1.0>,<+0.04,0.0,1.0>,<+0.04,0.0,0.0>],
[<-0.02,-0.03464,0.0>,<-0.02,-0.03464,1.0>,<0.02,0.03464,1.0>,<0.02,0.03464,0.0>],
[<-0.02,0.03464,0.0>,<-0.02,0.03464,1.0>,<0.02,-0.03464,1.0>,<0.02,-0.03464,0.0>],
# an octagon approximating a disk at the bottom of the arrow
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] )
):
ARROWHULL1 := map(Matrix,ARROWHULL1):
ARROWHULL1 := map(convert,ARROWHULL1,listlist):
ARROWHULL1 := map(hfarray,ARROWHULL1):
#-----------------------------------------------------------------------
#two directiom
A := evalf(Pi/4): N := 8: H := 1.0:
ARROWHULL2:= POLYGONS(
# two perpendicular triangles approximating a cone for the head of the arrow
[<+0.04,0.0,0.8>,<-0.04,0.0,0.8>,<-0.1,0.0,0.65>,<0.1,0.0,0.65>],
[<0.02,0.03464,0.8>,<-0.02,-0.03464,0.8>,<-0.05,-0.0866,0.65>,<0.05,0.0866,0.65>],
[<0.02,-0.03464,0.8>,<-0.02,0.03464,0.8>,<-0.05,0.0866,0.65>,<0.05,-0.0866,0.65>],
[<+0.04,0.0,0.2>,<-0.04,0.0,0.2>,<-0.1,0.0,0.35>,<0.1,0.0,0.35>],
[<0.02,0.03464,0.2>,<-0.02,-0.03464,0.2>,<-0.05,-0.0866,0.35>,<0.05,0.0866,0.35>],
[<0.02,-0.03464,0.2>,<-0.02,0.03464,0.2>,<-0.05,0.0866,0.35>,<0.05,-0.0866,0.35>],
# an octagon approximating a disk at the base of the arrow head
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ),
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ),
# two perpendicular rectangles approximating a cylinder for the body of the arrow
[<-0.04,0.0,0.8>,<-0.04,0.0,1.0>,<+0.04,0.0,1.0>,<+0.04,0.0,0.8>],
[<-0.02,-0.03464,0.8>,<-0.02,-0.03464,1.0>,<0.02,0.03464,1.0>,<0.02,0.03464,0.8>],
[<-0.02,0.03464,0.8>,<-0.02,0.03464,1.0>,<0.02,-0.03464,1.0>,<0.02,-0.03464,0.8>],
[<-0.04,0.0,0.35>,<-0.04,0.0,0.65>,<+0.04,0.0,0.65>,<+0.04,0.0,0.35>],
[<-0.02,-0.03464,0.35>,<-0.02,-0.03464,0.65>,<0.02,0.03464,0.65>,<0.02,0.03464,0.35>],
[<-0.02,0.03464,0.35>,<-0.02,0.03464,0.65>,<0.02,-0.03464,0.65>,<0.02,-0.03464,0.35>],
[<-0.04,0.0,0.0>,<-0.04,0.0,0.2>,<+0.04,0.0,0.2>,<+0.04,0.0,0.0>],
[<-0.02,-0.03464,0.0>,<-0.02,-0.03464,0.2>,<0.02,0.03464,0.2>,<0.02,0.03464,0.0>],
[<-0.02,0.03464,0.0>,<-0.02,0.03464,0.2>,<0.02,-0.03464,0.2>,<0.02,-0.03464,0.0>],
# an octagon approximating a disk at the bottom of the arrow
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ) ):
ARROWHULL2 := map(Matrix,ARROWHULL2):
ARROWHULL2 := map(convert,ARROWHULL2,listlist):
ARROWHULL2 := map(hfarray,ARROWHULL2):
#-----------------------------------------------------------------------
#one direction
A := evalf(Pi/4): N := 8: H := 1.0: ARROWHULL3 :=
CURVES(
# two perpendicular triangles approximating a cone for the head of the arrow
[<0.04,0.0,0.55>,<-0.04,0.0,0.55>,<-0.1,0.0,0.45>,<0.1,0.0,0.45>],
[<0.02,0.03464,0.55>,<-0.02,-0.03464,0.55>,<-0.05,-0.0866,0.45>,<0.05,0.0866,0.45>],
[<0.02,-0.03464,0.55>,<-0.02,0.03464,0.55>,<-0.05,0.0866,0.45>,<0.05,-0.0866,0.45>],
# an octagon approximating a disk at the base of the arrow head
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ),
# two perpendicular rectangles approximating a cylinder for the body of the arrow
[<-0.04,0.0,0.0>,<-0.04,0.0,0.45>,<+0.04,0.0,0.45>,<+0.04,0.0,0.0>],
[<-0.02,-0.03464,0.0>,<-0.02,-0.03464,0.45>,<0.02,0.03464,0.45>,<0.02,0.03464,0.0>],
[<-0.02,0.03464,0.0>,<-0.02,0.03464,0.45>,<0.02,-0.03464,0.45>,<0.02,-0.03464,0.0>],
[<-0.04,0.0,0.55>,<-0.04,0.0,1.0>,<+0.04,0.0,1.0>,<+0.04,0.0,0.55>],
[<-0.02,-0.03464,0.55>,<-0.02,-0.03464,1.0>,<0.02,0.03464,1.0>,<0.02,0.03464,0.55>],
[<-0.02,0.03464,0.55>,<-0.02,0.03464,1.0>,<0.02,-0.03464,1.0>,<0.02,-0.03464,0.55>],
# an octagon approximating a disk at the bottom of the arrow
evalf( [seq( , theta=[seq( i*A,i=0..N-1 )] )] ) ):
ARROWHULL3 := map(Matrix,ARROWHULL3):
ARROWHULL3 := map(convert,ARROWHULL3,listlist):
ARROWHULL3 := map(hfarray,ARROWHULL3):
macro( N=N, H=H, A=A ):
#----------------------------------------------------------------------
#savelib('NORM');
#savelib('vectorplot3d');
#savelib('FAST');
#savelib('ARROWHULL1');
#savelib('ARROWHULL2');
#savelib('ARROWHULL3');
#----------------------------------------------------------------------