BeginPackage["UVW`RandWalk`"]
(* Version 2.0, 08/15/97 *)

RandomWalk::usage= "
RandomWalk[ListofVelocities,Deltat] represents the trajectory of a point 
in a square. ListofVelocities is a list of two-dimensional vectors, 
interpreted as consecutive speeds for the point. The point starts from 
the center with the first speed vector of the list. It changes its speed 
vector for the next one in the list at each integer multiple of deltat."



VectorField::usage= "
VectorField[arrayof2Dvelocities] represents graphically by segments on a 
grid the values of a discrete vector field. Arrayof2Dvelocities is a list 
with three levels. The first two correspond to the coordinates (i,j) of a 
point on the grid. The last level corresponds to the two coordinates of the 
vector attached to point (i,j) : Vx(i,j) , Vy(i,j) . The function represents 
the grid and the vector attached to each point."



VectorFieldTrajectory::usage= "
VectorFieldTrajectory[arrayof2Dvelocities, deltat, tmax] represents first a 
vector field on a grid by calling VectorField[arrayof2Dvelocities]. Then 
it draws the trajectory of a point starting at the center of the grid. At 
each integer multiple of deltat, the velocity vector of the point is changed 
for that of the vector field at the closest point on the grid. The trajectory 
is followed up to time tmax. The boundary conditions are periodic."

	   


Begin["`Private`"]

 

RandomWalk[listofvelocities_List,deltat_,opts___]:=

      Block[{points},

	   points = listofvelocities deltat;
	   points = FoldList[Plus,{1.,1.},points];
	   points = Mod[points,2.] - Table[{1.,1.},{Length[points]}];
	   
	   ListPlot[points,
		  opts,
		  Prolog     -> PointSize[0.001],
		  PlotRange  -> {{-1,1},{-1,1}},
		  Axes       -> False,
		  Frame      -> True,
		  FrameTicks -> None]
		       
	   ]
				       


VectorField[arrayofvelocities_List,opts___]:=

      Block[{nx,ny,vect,graph,i,j},
	   
	   nx   = Length[arrayofvelocities];
	   ny   = Length[First[arrayofvelocities]];
	   vect = N[arrayofvelocities]; 
	   vect = vect/ (2 Max[Abs[vect]]);
	   
	   graph = Table[Line[{{i,j},{i,j}+vect[[i,j]]}],{i,nx},{j,ny}];

	   Show[Graphics[graph],
			      opts,
			      Prolog  -> Thickness[0.001],
			      PlotRange  -> {{0,nx+1},{0,ny+1}},
			      Axes       -> False,
			      Frame      -> True,
			      FrameTicks -> None]
	   ]


VectorFieldTrajectory[arrayofvelocities_List,deltat_,tmax_]:=

      Block[{nx,ny,maxx,maxy,nstep,vect,x,y,point,traj,i,j,gvect,gtraj},

	   nx    = Length[arrayofvelocities];
	   maxx  = nx - 1.;
	   ny    = Length[First[arrayofvelocities]];
	   maxy  = ny - 1.;
	   nstep = Floor[tmax/deltat];
	   x = nx/2.;
	   y = ny/2.;
	   point={x,y};
	   traj = {point};
	   vect = N[arrayofvelocities];

	   Do [
	      (
	       point = point + deltat*vect[[ Round[x]+1,Round[y]+1 ]];  
	       x = Mod[First[point],maxx];                                              
	       y = Mod[Last[point],maxy];
	       point = {x,y};
	       traj = Append[traj,point]
	      ),
	      {nstep}
	      ];

     
	   traj=traj+Table[{0.5,0.5},{Length[traj]}];
	   
	   vect  = vect /(2 Max[Abs[vect]]);

	   gvect = Table[Line[{{i,j},{i,j}+vect[[i,j]]}],{i,nx},{j,ny}];


	   gtraj = Table[Point[traj[[i]]],{i,Length[traj]}]; 

	   Show[Graphics[gvect],Graphics[gtraj],
			  
			      Prolog  -> Thickness[0.001],
			      PlotRange  -> {{0,nx+1},{0,ny+1}},
			      Axes       -> False,
			      Frame      -> True,
			      FrameTicks -> None]        
	   
	   ];


End[ ]

EndPackage[ ]