Random walk in a 3d hexagonal close-packed (cannonball) lattice.

Let each integer position on the [Graphics:Images/hewxalk_gr_1.gif] axis of a cartesian [Graphics:Images/hewxalk_gr_2.gif] coordinate system be a lattice point.  In the [Graphics:Images/hewxalk_gr_3.gif] plane pack additional rows so that the new lattice points are eactly one unit from all their neighbors.  An angle of 60 degrees and radius one, or [Graphics:Images/hewxalk_gr_4.gif] in rectangular notation.

In the third dimension we need to find the point in the [Graphics:Images/hewxalk_gr_5.gif] half of the [Graphics:Images/hewxalk_gr_6.gif] plane that is exactly one unit from the origin and the point [Graphics:Images/hewxalk_gr_7.gif].

[Graphics:Images/hewxalk_gr_8.gif]
[Graphics:Images/hewxalk_gr_9.gif]
[Graphics:Images/hewxalk_gr_10.gif]

Store the vectors explicitly ...

[Graphics:Images/hewxalk_gr_11.gif]
[Graphics:Images/hewxalk_gr_12.gif]
[Graphics:Images/hewxalk_gr_13.gif]
[Graphics:Images/hewxalk_gr_14.gif]

These vectors are not orthogonal, so we cannot simply treat it as a coordinate system.  It will be sufficient to have an algorithm for converting a path description in terms of integer multiples of these vectors to Cartesian coordinates.  Vector addition will do, but we will develop an integer arithemtic-oriented algorithm.

Let us construct a synthetic orthogonal grid coordinate system.  In the [Graphics:Images/hewxalk_gr_15.gif] plane, the row of lattice points along the [Graphics:Images/hewxalk_gr_16.gif] axis is a one-dimensional unit.  Along the synthetic [Graphics:Images/hewxalk_gr_17.gif] line is a copy of this lattice row with each node displaced by one [Graphics:Images/hewxalk_gr_18.gif] from the corresponding node in the base row.  In the [Graphics:Images/hewxalk_gr_19.gif] row the nodes are again aligned with those in the base row: the appropriate offset is:

[Graphics:Images/hewxalk_gr_20.gif]
[Graphics:Images/hewxalk_gr_21.gif]

A similar treatment applies in the [Graphics:Images/hewxalk_gr_22.gif] direction.  The base [Graphics:Images/hewxalk_gr_23.gif] lattice plane is replicated at mutliple of [Graphics:Images/hewxalk_gr_24.gif], with 3 steps along [Graphics:Images/hewxalk_gr_25.gif] bringing us to a plane aligned exactly with the base plane:

[Graphics:Images/hewxalk_gr_26.gif]
[Graphics:Images/hewxalk_gr_27.gif]

We can exploit the above to efficiently compute the location given a path description.  It is sufficient to record the net motion in each of the lattice bond directions.  The Cartesian coordinates for a lattice point at [Graphics:Images/hewxalk_gr_28.gif] can be efficiently computed as:

[Graphics:Images/hewxalk_gr_29.gif]
[Graphics:Images/hewxalk_gr_30.gif]
[Graphics:Images/hewxalk_gr_31.gif]

Each lattice point has twelve neighbors at a distance of exactly one unit.  These can be generated by taking the differences of all ordered pairs selected without replacement from the set consisting of the three bond direction vectors and the zero vector:

[Graphics:Images/hewxalk_gr_32.gif]
[Graphics:Images/hewxalk_gr_33.gif]

We can denote one of these twelve neighbor vectors by a string of length three from the alphabet [Graphics:Images/hewxalk_gr_34.gif] where no two members of the string may have a positive product.  The vector corresponding to the string [Graphics:Images/hewxalk_gr_35.gif] is then [Graphics:Images/hewxalk_gr_36.gif], which is necessarily then one of the vectors constructed above.

A random walk in this latice can be generated by a Markov chain over these strings.  Suppose we wish to impose a stochastic limit on the curvature of the walk in 3-space.  Given the step vector that brought us to the current node, consider the projection (dot product) of each of the neighbor vectors onto this arrival vector.  The vector pointing to the previous node will have a projection of [Graphics:Images/hewxalk_gr_37.gif] while the arrival vector will of course have projection one.  There are two vectors orthogonal to the arrival vector.  Of the remaining eight vectors, four will have projection [Graphics:Images/hewxalk_gr_38.gif] and the other four will have projection [Graphics:Images/hewxalk_gr_39.gif].

[Graphics:Images/hewxalk_gr_40.gif]
[Graphics:Images/hewxalk_gr_41.gif]

It is a relatively simple matter to tabulate the projection relation in terms of the vector description string and to assign transition probabilities based on projection equivalence classes.  If angles more acute than [Graphics:Images/hewxalk_gr_42.gif] are to be avoided one may assign a probability to going straight and divide the complement equally among the remaining four forward vectors.

It is noteworthy that the four forward curve directions are not uniformly distributed about the  arrival vector.  If the observer is situated at the current node and looking in the direction of the arrival vector and oriented so that one of the forward curve vectors points toward the hoizon [Graphics:Images/hewxalk_gr_43.gif] left of center, then another points to the horizon [Graphics:Images/hewxalk_gr_44.gif] right of center.  Of the other two, one is above the horizon on one side while the last is below the horizon on the opposite side.  This picture is not changed when stepping along the arrival vector and can be thought of as a kind of bias in the lattice; a bias against the underrepresented quadrants of the local sky and ground.  This should have no macroscopic significance.

[Graphics:Images/hewxalk_gr_45.gif]
[Graphics:Images/hewxalk_gr_46.gif]

[Graphics:Images/hewxalk_gr_47.gif]

[Graphics:Images/hewxalk_gr_48.gif]

[Graphics:Images/hewxalk_gr_49.gif]

Here is the view of the forward vector targets described above.

[Graphics:Images/hewxalk_gr_50.gif]

[Graphics:Images/hewxalk_gr_51.gif]

In the simulation of vorticity filaments we need to know the midpoints of the vectors, so we will have our walk make two steps along each vector.  We will parameterize our walk by the probability of going directly forward, with the complementary events equally divided among the four off-forward directions.

We begin by tabulating the off-forward vectors.  If our arrival vector is desctibed by the string [Graphics:Images/hewxalk_gr_52.gif] then the simplest way to find the foward vectors is to generate a list of all possible descritions and select those that have the desired projection.

[Graphics:Images/hewxalk_gr_53.gif]
[Graphics:Images/hewxalk_gr_54.gif]
[Graphics:Images/hewxalk_gr_55.gif]
[Graphics:Images/hewxalk_gr_56.gif]

Given an arrival vector in string form we generate the departure vector description as follows:

[Graphics:Images/hewxalk_gr_57.gif]
[Graphics:Images/hewxalk_gr_58.gif]
[Graphics:Images/hewxalk_gr_59.gif]
[Graphics:Images/hewxalk_gr_60.gif]

Given a starting position and direction, both in path description form, we generate a walk of length [Graphics:Images/hewxalk_gr_61.gif] as follows:

[Graphics:Images/hewxalk_gr_62.gif]
[Graphics:Images/hewxalk_gr_63.gif]

[Graphics:Images/hewxalk_gr_64.gif]

We will generate our walks by starting at the origin and going in a randomly selected initial direction, then offsetting the entire walk to a randomly selected inital point near the origin.

[Graphics:Images/hewxalk_gr_65.gif]
[Graphics:Images/hewxalk_gr_66.gif]

To visualize a pair of walks embedded in the lattice it is sufficient to concatenate their list representations and feed the result to vwalk:

[Graphics:Images/hewxalk_gr_67.gif]

[Graphics:Images/hewxalk_gr_68.gif]

With all vectors having length two the dot product of any pair will be one of [Graphics:Images/hewxalk_gr_69.gif].  We now tabulate these dot products by descriptions string:

[Graphics:Images/hewxalk_gr_70.gif]
[Graphics:Images/hewxalk_gr_71.gif]

Given a pair of filaments (walks) in the lattice, we compute a statistic correspondiing to the interaction energy of the filaments.  For each pair of vectors with one member chosen from each filament, the energy is the dot product of the vectors divided by the distance between their midpoints.  We therefore begin our computation by partitioning the vectors of each filament into subsets by orientation.  In the resulting representation the orientation information will be implicit.

[Graphics:Images/hewxalk_gr_72.gif]
[Graphics:Images/hewxalk_gr_73.gif]

Given two walks in set representation, we can compute the total energy by constructing all pairs of vectors grouped by orientation.  There are 144 ways to choose one orientation from each filament, and we need a list of the distances between all pairs of vectors having the corresponding orientations.  We will avoid constructing this list for those orientation pairs that have are orthogonal, since they contribute nothing to the energy.

The contribution of each vector pair within a set having identical orientations is inversely proportional to the distance between center points.  To keep everything in integer form we will store the norm squared rather than the vector or its norm in for each pair.

[Graphics:Images/hewxalk_gr_74.gif]
[Graphics:Images/hewxalk_gr_75.gif]

Let's have a look at the spectrum for the pair of small filaments in the figure above.

[Graphics:Images/hewxalk_gr_76.gif]
[Graphics:Images/hewxalk_gr_77.gif]
[Graphics:Images/hewxalk_gr_78.gif]
[Graphics:Images/hewxalk_gr_79.gif]
[Graphics:Images/hewxalk_gr_80.gif]

[Graphics:Images/hewxalk_gr_81.gif]

To compute the total energy we need to divide the bin counts in this spectrum by the reciprocal of the corresponding norm.  Since the table is indexed by the square of the norm, let's build a table of the reciprocal square roots:

[Graphics:Images/hewxalk_gr_82.gif]

The energy then is given by the sum of the elements of the term-by-term product of these two lists.  Unfortunately this list is sorted with the largest elements first, statistically speaking, so it would be more numerically sensible to add it up in reverse order.

[Graphics:Images/hewxalk_gr_83.gif]
[Graphics:Images/hewxalk_gr_84.gif]

Let's see how long it takes to compute the energy for a larger scale example:

[Graphics:Images/hewxalk_gr_85.gif]
[Graphics:Images/hewxalk_gr_86.gif]

Let's see where the time is going:

[Graphics:Images/hewxalk_gr_87.gif]
[Graphics:Images/hewxalk_gr_88.gif]
[Graphics:Images/hewxalk_gr_89.gif]
[Graphics:Images/hewxalk_gr_90.gif]

Let's try to speed this up by building an optimized version of the difnorm2 function.  Here is the core of pathtoloc and the norm squared of the result:

[Graphics:Images/hewxalk_gr_91.gif]
[Graphics:Images/hewxalk_gr_92.gif]
[Graphics:Images/hewxalk_gr_93.gif]

As we observed, this is an integer expression.

[Graphics:Images/hewxalk_gr_94.gif]
[Graphics:Images/hewxalk_gr_95.gif]
[Graphics:Images/hewxalk_gr_96.gif]
[Graphics:Images/hewxalk_gr_97.gif]

To speed this up further will require a different approach to accumulating the spectrum or directly computing the sum without the intermediate spectral representation.  I have not yet found a better algorithm in Mathematica.  Note that the construction of the array of squared distances is only half the problem:

[Graphics:Images/hewxalk_gr_98.gif]
[Graphics:Images/hewxalk_gr_99.gif]
[Graphics:Images/hewxalk_gr_100.gif]
[Graphics:Images/hewxalk_gr_101.gif]
[Graphics:Images/hewxalk_gr_102.gif]

This version is the speed champion, but it does not create a spectrum.  Let's see if we can get a spectrum with reasonable additional overhead:

[Graphics:Images/hewxalk_gr_103.gif]
[Graphics:Images/hewxalk_gr_104.gif]
[Graphics:Images/hewxalk_gr_105.gif]
[Graphics:Images/hewxalk_gr_106.gif]

Well, at least it's correct.  It isn't appreciably faster than the original, though.


Converted by Mathematica      January 6, 2001