For rotation of porphyroblasts, we have adopted the developed by Jeffery (1922). Details of the mathematical description can also be found in [Jezek et al, 1999]. Figure 1a illustrates the velocity field and velocity gradient tensor (`V'`

) with the external coordinate system (x'y'z'). Jeffery derived the sets of equations for rigid objects rotating in a viscous fluid matrix;_{ij}

_{1}= B

_{1}E'

_{32}- Ω'

_{32}

_{2}= B

_{2}E'

_{13}- Ω'

_{13}

_{3}= B

_{3}E'

_{21}- Ω'

_{21}

where `w`

represents the angular velocities of an rigid objects around the internal axes (xyz, Figure 1b), and where:
_{i}

and:

The parameters `B`

, _{1}`B`

, and _{2}`B`

which are related to the shape of rigid objects, respectively, represents _{3}`(b`

, ^{2}-c^{2})/(b^{2}+c^{2})`(c`

, and ^{2}-a^{2})/(c^{2}+a^{2})`(a`

, and ^{2}-b^{2})/(a^{2}+b^{2})`a`

, `b`

, and `c`

are axes of an ellipsoid object. Description for the orientation of rigid objects can be made much simpler by adopting Euler angles that relate internal coordinate system for rigid objects (xyz) and external coordinate system (x'y'z'). Freeman (1985) provided the solution for the rotation of rigid objects with Euler angles;

_{1}·sin y + w

_{2}·cos y) / sin q

_{1}·cos y - w

_{2}·sin y

_{3}- cos q

and `f`

, `q`

and `y`

represents the Euler angles defined in Figure 1b and the quantities divided by `dt`

represent the rate of changes in three Euler angles.

In two dimensions, an area occupied by a crystal with orthorhombic symmetry (a crystal with two sets of mutually perpendicular faces) can be represented with two sets of two parallel lines (Figure 2). If we let (`h`

) to represent the Miller index of a face, the region between the two parallel faces can be described as,_{1}k_{1}0

_{1}·x + k

_{1}·y + 0·z < |A|

where A is a parameter related to the distance between the face and the origin. Similarly, for another face with the Miller index of (`h`

), the region between these two faces is,_{2}k_{2}0

_{2}·x + k

_{2}·y + 0·z < |B|

where B is also a parameter similar to A. If a point satisfies the above two inequalities simultaneously, then the point is within the region surrounded by four faces. If we assign larger values for A and B (e.g. new A = old A + growth amount), the region surrounded by the four faces can expand, thus the crystal grows.

For three dimensions, we need three inequalities to define a volume occupied by a crystal with the faces of the Miller indices (`h`

), (_{1}k_{1}l_{1}`h`

), and (_{2}k_{2}l_{2}`h`

)._{3}k_{3}l_{3}

_{1}·x + k

_{1}·y + l

_{1}·z < |A|

_{2}·x + k

_{2}·y + l

_{2}·z < |B|

_{3}·x + k

_{3}·y + l

_{3}·z < |C|

where A, B, and C are parameters related to the distance between the face and the origin. A rectangular-box-shaped porphyroblast can grow by assigning larger values for A, B and C. The shape of the porphyroblast can be determined by relative magnitudes of A, B, and C.

The data structure for the model (Figure 2) is based on pixels (or volumetric pixels). The pixels within the model system have their coordinate information (x'y'z') and additional material information. The material information for each pixel can be among the following four types;

porphyroblast,

internal foliation within porphyroblast,

matrix,

external foliation within matrix (or future inclusions when they are entrapped).

Development of internal foliation patterns is performed by repeated growth and rotation (Figure 3). During growth, parameters of A, B, and C (described in the earlier section) becomes larger in each step of growth. Depending on the material information, the program replaces matrix with porphyroblast or external foliation with internal foliation as the crystal faces expand. To apply rotation for a porphyroblast, either porphyroblast or matrix can rotate. Here, we have adopted matrix rotation because it is much simpler but still gives the same internal fabric. If pixels have material information of porphyroblast or internal foliation, they remain stationary. On the other hand, pixels occupied by external foliation rotate backward with the amount calculated using Euler angles.

One of the key assumptions made in this study is no dragging of external foliation by the rotation of a porphyroblast. This assumption is physically unrealistic, but we made this assumption because this makes algorithm much simpler and still produces complicated enough textures that may give some insights when interpreting natural textures. Therefore, the results made in this study would be simpler compared with the results made by studies that consider dragging of external foliation [Jezek et al., 1999]. However, these studies also have limitations because of the assumption of Newtonian rheology of the matrix. For successful modelling of natural rocks, a model has to be fully mechanical in order to model non-linear viscous flow in the matrix.