- Email: [email protected]
Software for modeling the motion of rigid triaxial ellipsoidal particles in viscous flow
Computers & Geosciences Vol. 20, No. 3, pp. 409~,24, 1994
Pergamon
0098-3004(93)E0017-A
Copyright © 1994ElsevierScienceLtd Printed in Great Britain.All rights reserved 0098-3004/94$7.00+ 0.00
S O F T W A R E FOR M O D E L I N G THE M O T I O N OF R I G I D T R I A X I A L ELLIPSOIDAL PARTICLES IN VISCOUS FLOW JOSEF JE~,EK Department of Applied Mathematics and Computer Science, Faculty of Science, Charles University, Albertov 6, 128 43 Praha 2, Czech Republic (Received 31 October 1993)
Abstract--Three MS-FORTRAN programs are presented which enable study of the behavior of rigid triaxial ellipsoidal particles in slowly moving viscous flow. The modeling is based on numerical solution of Jeffery's equations. These programs allow the investigation of individual particle motion as well as the evolution of a multiparticle system composed of noninteracting particles. In the combined pure shear and simple shear flow, it is possible to estimate the character of multiparticle system evolution by diagrams of time derivatives of Euler angles. The programs are interactive with graphical output on the screen and run on a standard IBM-PC with VGA. Some examples of graphical output of the modeling are presented. Key Words: Ellipsoidal particle motion, Viscous flow, Multiparticle system motion, Pure shear, Simple
shear
INTRODUCTION In the interpretation of preferred orientations of fabric elements in tectonites the model of Jeffery (1923) is used. Jeffery's model describes the behavior of a rigid ellipsoidal particle in slow viscous flow. Given the flow geometry (velocity gradient tensor), the particle axial ratio, and its starting position, the trajectory of the particle during the flow can be determined as the solution of Jeffery's equations. Analytical solutions of these equations of the particle motion have been determined only for biaxial particles (i.e. ellipsoids of revolution) and restricted types of flow geometry: simple shear flow (Jeffery, 1923), pure shear flow (Gay, 1968) axial flattening (Debat and others, 1975; Tulis, 1976), and general orthorhombic flow (Fernandez, 1988). These analytical solutions were used in modeling the motion of an individual biaxial particle and also in modeling multiparticle systems composed of noninteracting particles (Reed and Tryggvason, 1974; Tullis, 1976). In the situation of a triaxial ellipsoidal particle, numerical techniques can be used for solving Jeffery's equations as in Hinch and Leal (1979). Freeman (1985), and Fernandez and Laporte (1991). The results of investigating the motion of an individual triaxial particle with different initial orientations in different types of flow show a significant deviation of the particle behavior from that of a biaxial particle and the need of further work in this field. Modeling of a multiparticle system composed of noninteracting rind ellipsoidal triaxial particles has been also suggested by Freeman (1985), but such a work has not appeared up to this time.
The aim of this paper is the presentation of software enabling the user to investigate the motion of an individual triaxial ellipsoidal particle as well as the evolution of a multiparticle system composed of noninteracting triaxial particles, in different flow geometries. The programs presented could serve also as a tool for estimating the possible development of preferred orientation in real conditions. They were written in MS-FORTRAN, v. 5.0, and run on a standard IBM PC. The Intel 80486 or mathematical coprocessor is recommended for satisfactory speed. The programs are interactive, demand necessary input parameters from console, and yield graphical output of results on the screen. Program JEFF listed in the Appendix enables the modeling of an individual particle motion. The motion of the particle is drawn on the screen as the path of one of the particle axes in the selected type of projection. For solving the differential equations of the particle motion, subroutine DVERK from the IMSL Library, performing the fourth- and fifth-order Runge-Kutta method, was used. This subroutine is not listed here. The user also should provide his/her own subroutine performing print-screen on his type of printer. Program MUPLES allows the user to model and visualize a multiparticle system evolution. Because of its affinity to the program JEFF and the finite extent of the paper it is described only briefly and not listed, and the program TIDER, drawmg patterns of Euler angles derivatives, as well. These diagrams help in estimating the character of multiparticle system evolution, especially in the situation of combined pure shear and simple shear model,
409
410
J. JE~EK
which can be used as an approximation in many real situations. THE MODEL OF ,IEFFERY
The equations describing our model were given by Jeffery (1923). In his work, he considered a slowly moving Newtonian viscous fluid in which sufficiently small rigid and noninteracting particles of an ellipsoidal shape are immersed. Starting with the Navier-Stokes equations of motion of incompressible viscous fluid and neglecting the terms of the second order, Jeffery determined the equations of motion of the individual particle. Consider a right-handed orthogonal coordinate system (x ~, x ~, x ~) positioned at the particle mass center with axes in the directions fixed in space (further termed "fixed coordinate system") and let another coordinate system (x~, x2, x3) be coincident with the axes a, b, c of the ellipsoidal particle (further termed "rotating coordinate system"). Denote the velocity gradient tensor characterizing undisturbed flow at a large distance from the particle by V ~ = Ov~lt3x~, and suppose that V~ is constant in time and obeys the equation of continuity. Denote further by E~ and f ~ the symmetric and antisymmetric part of the velocity gradient tensor, respectively. Then, in the rotating coordinate system, the Jeffery equations can be written as ~01 = Bt E32 - ~'~32, 092 = B2E13
-
-
~'~13'
(1)
These equations can be solved numerically as an initial value problem. Given the velocity gradient tensor, axial ratio of the particle, and initial values of the Euler angles 4~, 0, ~ (the initial position of the particle), we compute the time evolution of these angles, that is the particle rotation.
INDIVIDUAL PARTICLE MOTION (Program JEFF)
For the numerical solution of the Equations (2) and display of the resulting particle motion, program J E F F was written. The program demands the following input parameters from console: --velocity gradient tensor V ~, --semiaxes of the ellipsoidal particle a, b, c, --initial Euler angles ~b, 0, ~, of the particle, - - t i m e increment for the Runge-Kutta method, - - n u m b e r of time steps performed at once (before the program again demands for continuation), --which of the particle axes will be traced (its path will be drawn on the screen), - - t y p e of projection for the display of this particle path (stereographic, Lambert, or perspective p l o ~ s e e Figs. 1 and 2), --optionally, if the c-axis is traced, the position of the a-axis can be sketched along the c-axis path, which allows to imagine the triaxial particle rotation around its c-axis.
093 = B3 E21 - ~21, where eg~, i = 1,2,3, are the angular velocities of the particle around its axes and the vector B l = (b 2 _ c 2)/(b 2 + c 2),
B2 = (c 2 - a 2)/(c 2 + a 2), B 3 = (a 2 _ b 2)/(a 2 + b 2)
corresponds to the axial ratio of the ellipsoidal particle. An analytical solution of Jeffery's equations has been determined only for biaxial particle in some flow geometry as was mentioned in the Introduction. These equations, however, can be reformulated to the form convenient for numerical solution. Describing the rotation of the particle with respect to the fixed coordinate system by Euler angles q% 0, ~b and using the known relations between the angular velocity ~0 and the time derivatives of the Euler angles we obtain three differential equations corresponding to (1) 4; = (~ol sin @ + ~o2cos ~, )/sin 0, = ol cos ~ - o~2 sin $, = ~o3 -
4; c o s 0.
(2)
The program is interactive. The path of a traced particle is drawn on the screen for the given number of time increments of the Runge-Kutta method and after that the program asks if continuation or printscreen is desired. If no continuation is selected, the program offers to rerun the model or exit. If rerun is desired, the user can keep some parameters from the preceding run. For solving the differential Equations (2), subroutine DVERK of the IMSL Library, performing the fourth- and fifth-order Runge-Kutta method, was precise by a comparison of its results with some solutions known analytically. The recommended value of the time step for the Runge-Kutta method is 0.1 Subroutine DVERK is not included in the listing of program J E F F (see Appendix) and neither is subroutine TOPRN for print-screen. The user should provide his/her own subroutine performing print-screen on his/her type of printer. Standard VGA is required in subroutines for graphical output on the screen. Examples of use of the program J E F F are in Figures 1 and 2. In Figure 1 one of classical Jeffery's orbits of a biaxial particle is drawn in perspective plot, representing the fact that biaxial particles in simple shear flow follow closed orbits and their
Rigid triaxial ellipsoidal particle motion
411
Ellipsoidal Particle Motion
Uelocity Gradient .000 .050 .000 Axial
Tensor
.000 .000 .OOO .OOO .000 .000 Ratio
1.00 : 1.00 : 3 . 0 0 Euler Rngles
Initial
O.
30.
O.
Time Increment :
.100
Time S t e p
SO00
:
Traced R x i s
:
C
/~//' ,
"~-'~-~--'--~-~
///// /
Figure 1. Jeffery's orbit of biaxial particle in perspective plot. See description of program JEFF for explanation of parameters of model.
motion is periodic. In Figure 2 a triaxial particle was traced. In the Lambert equal-area projection used, the a-axis has been sketched along the c-axis path, which allows to imagine the rotation of the particle around its c-axis. The same simple shear flow as in Figure 1 was used, but the triaxial particle motion is no longer periodic.
Ellipsoidal
Particle
MULTIPARTICLE SYSTEM EVOLUTION
(Program MUPLES) Supposing that there are no interactions among the individual particles rotating in viscous fluid, the multiparticle system can be investigated by Equations (2). For this purpose program M U P L E S was written.
X
Motion
I
Uelocity .000 .050 .000 Axial
Gradient
Tensor
~
.000 .OOO .OOO . ~ ,090 .000 Ratio
.67
Initial G.
: 1.50
: 3.00
¥
Failer nmjles 30.
0.
Time Increment :
.IO0
Time Step :
SO00
Traced
7
/
Axis : C
Figure 2. Behavior of triaxial particle in same flow as in Figure 1. In Lambert equal-area projection used, a-axis has been sketched along c-axis path.
412
J. JE~EK
Four - Groups System Motion
x
U e l o c i t y Gradient Tensor .010 .000 .050 - . 0 1 0 .000 .000
.GO0 .OO0 .000
o
[] \
o 0
1
:
2
:
3
o
1
:
Z
:
4
,~
1
:
4
:
6
[]
1
:
Z
:
6
STEP :
/
X
o
[]
\\o ~
A
o ~o ~>o o A
<6 no ~ ~
[]
\
0
~
o
o
6
1
/ /,'
100
Figure 3. Modeling of multiparticle system evolution by program MUPLES. Multiparticle system is composed of four groups of particles, each of them having about 100 particles of same axial ratio. Current positions of c-axes of particles are marked by symbols in Lambert equal-area projection. Different symbols correspond to different axial ratios. Because of the finite extent of this paper, the program MUPLES is not listed and only a brief description will be given here. The program can be created using some subroutines of program JEFF. The program input parameters are the velocity gradient tensor, time increment for the Runge-Kutta method, number of time steps performed at once, particle axis to be traced and type of projection. Initial Euler angles of the multiparticle system can be read from an ASCII file or randomly generated. A convenient selection is an initial is•tropic distribution of angles tk, 0 (similar to Sanderson and Meneilly, 1981) and a random distribution of the angle ~. During the time evolution of the multiparticle system studied, after the number of Runge-Kutta steps computed at once, the selected axis positions for all the particles are marked in the desired type of projection. Optionally, the orientation tensor (for definition, see e.g. Sanderson and Meneilly, 1981) and its eigenvalues and eigenvectors or intensity of the fabric (Lisle, 1985) can be computed and presented. In Figure 3 an example of output from program MUPLES is shown. A multiparticle system was composed of four groups of particles, each of them having about 100 particles of the same axial ratio. The system started with an initial is•tropic distribution of angles ~, 0 and a random distribution of the angle ~, which determine the position of the particle. EVALUATING
THE EULER ANGLES DERIVATIVES
TIME
(Program TIDER) The mentioned programs JEFF and MUPLES
compute and visualize the time evolution of a oneor multiparticle system. Program TIDER allows another type of view on the system development based on the fact that this model is stationary. The velocity gradient tensor V~. is constant in the fixed coordinate system and thus for particles of a given axial ratio the right-hand sides of the Equations (2) depend only on the Euler angles, that is on the instantaneous position of particles. The analysis of this dependence can help estimate the future or final state of the one- or multiparticle system. This is valid, for example in the combined pure shear and simple shear model given by the velocity gradient tensor
v~=
~
-~
0
0
1
(3)
where ~ and 7 are constant in time. Such a model may approximate many real geological situations. Therefore program TIDER was written.to p!ot the distribution patterns of the derivatives ~, 0, ~ given by the Equations (2) in an equal-area stere•graphic projection. Input parameters are the velocity gradient tensor V ~and the lengths of the particle semiaxes a, b, c. For the possibility of efficient print-screen the dithering technique was used in plotting these diagrams. Because of the finite extent of this paper, program TIDER is not listed here. The program uses some subroutines of program JEFF. As an example of the usefulness of program TIDER, Figure 4 shows a ~b-derivative pattern corresponding to a combined pure shear and simple shear flow. The
Rigid triaxial ellipsoidal particle motion
dPH/dt
Field
Oelocitg 6 r a d i e n t .010 .000 .050 -.el0 •0 0 0 .000 Axial
413
Tensor
.0~ .000 .000
Batio .33
: 1.50
: 6.00
dPH/dt
< -Z.O00 < -1.75o < -1.5oo <
< < <
< <
- 1 .zso -1.000
-no -.500 -.250 . O~
•
I M
.....
Figure 4. ~b-Derivative pattern corresponding to combined pure shear and simple shear flow. In gray-scale used white corresponds to nonnegative values of ~b and all gray grades correspond to negative values of
9.
values of ~b-derivative were computed for an isotropic distribution of angles 4), 0 and a random distribution of the angle ~k, and then plotted in the stereogram. In the gray-scale used the white corresponds to nonnegative values of ~b and gray grades correspond to negative values of qk The ~b-derivative pattern can be considered as representation of a multiparticle system composed of triaxial particles of the same axial ratio, with the position of particles given by the angles 4), 0, and ~b. The value of ~b has a strong influence on the next motion of the particles. In the next evolution of the system, the particles having nonnegative values of ~b will remain in their position or will rotate in such a way that their c-axis projection will move clockwise. On the other side, the projection of the c-axis of the particles having negative values of ~b will move counterclockwise. Ana!yzing simultaneously the patterns of the derivatives q~, 0, we can deduce the type of evolution of a given system. In the given example a stable preferred orientation will arise during the multiparticle system evolution. CONCLUSION The presented programs J E F F , MUPLES, and T I D E R can be used for a detailed discussion o f the one- and multiparticle system evolution in different flow geometries. The possibility of the preferred orientation development in the combined pure shear and simple shear flow was investigated by
Je2ek and others (1994). As a natural example to which the Jeffery's model can be applied, ultramylonites were mentioned. In the ultramylonites the matrix/clast proportion may be high. Thus, it is possible to investigate the rotation of relict clasts in the time of the development of fabric using the Jeffery's model.
Acknowledgments--I would like to thank Karel Schulmann who initiated this work. He and Radek Melka introduced me into this problem and had clear requirements on the capabilities of programs. Also, I thank Karel Segeth who kindly reviewed the paper.
REFERENCES Debat, P., Sirieys, P., Deramond, J., and Soula, J. C., 1975, Paleodrformations d'un massif orthogenessique (massif des Cammazes. Montagne Noire Occidentale. France): Tectonophysics, v. 28, no. 3, p. 159-183. Fernandez, A., 1988, Strain analysis from shape preferred orientation in magmatic flows: Bull. Geol. Inst. Univ. Uppsala, v. 14, p. 61~57. Fernandez, A., and Laporte, D., 1991, Significance of low symmetry fabrics in magrnatic flows: Jour. Struct. Geology, v. 13, no. 3, p. 337 347. Freeman, B., 1985, The motion of rigid ellipsoidal particles in slow flows: Tectonophysics, v. 113, no. 2, p. 163-183. Gay, N. C., 1968, The motion of rigid particles embedded in a viscous fluid during pure shear deformation of the fluid: Tectonophysics, v. 5, no. 2, p. 81-88. Hinch, E. J., and Leal, L. G., 1979, Rotation of small non-axisymetric particles in a simple shear flow: Jour. Fluid Mech., v. 92, pt. 3, p. 591~08.
414
J. JE~EK
Jeffery, G. B., 1923, The motion of ellipsoidal particles immersed in a viscous fluid: Proc. Roy. Soc. London, Ser. A, v. 102, p. 161 179. Je~ek, J., Melka, R., Schulmann, K., and Venera, Z., 1994, The behaviour of rigid triaxial ellipsoidal particles in viscous flows-modeling of fabric evolution in multiparticle system: Tectonophysics, v. 229, p. 165-180. Lisle, R. J., 1985, The use of the orientation tensor for the description and statistical testing of fabrics: Jour. Struct. Geology, v. 7, no. 1, p. 115-117.
Reed, L. J., and Tryggvason, E., 1974, Preferred orientations of rigid particles in a viscous matrix deformed by pure shear and simple shear: Tectonophysics, v. 24, no. 2, p. 85-98. Sanderson, D. J., and Meneilly, A. W., 1981, Analysis of three-dimensional strain modified uniform distributions: andalusite fabrics from a granite aureole: Jour. Struct. Geology, v. 3, no. 2, p. 109-116. Tullis, T. E., 1976, Experiments on the origin of slatey cleavage and schistosity: Geol. Soc. America Bull., v. 87, no. 5, p. 745-753.
APPENDIX Program Listing C C C C C C C C C C C C C C
P R O G R A M
JEFF
N U M E R I C A L S O L U T I O N OF J E F F E R Y ' S E Q U A T I O N S F O R T H E M O T I O N OF A R I G I D E L L I P S O I D A L P A R T I C L E IN S L O W F L O W S R U N S O N I B M - T Y P E PC W I T H VGA. M A T H E M A T I C A L C O P R O C E S S O R OR I N T E L 80486 IS R E C O M M E N D E D F O R S A T I S F A C T O R Y S P E E D INPUT
PARAMETERS:
SEE
SUBROUTINE
INPPAR
P R O G R A M U S E S S U B R O U T I N E D V E R K F R O M I M S L (NOT L I S T E D HERE) U S E R S H O U L D I N C L U D E A S U B R O U T I N E P E R F O R M I N G P R I N T - S C R E E N O N HIS O W N T Y P E OF P R I N T E R
***************************************************************************
INCLUDE 'FGRAPH.FI' INCLUDE 'FGRAPH.FD' INTEGER* 2 DUMMY C H A R A C T E R ANS, S T R 6 * 6 DOUBLE PRECISION X,Y(3),XEND,TOL,YPR(3),C(24),W(3,9),0(3) RECORD /RCCOORD/ CURPOS EXTERNAL FCN D I M E N S I O N D V ( 3 , 3 ) ,EF (3,3) ,OF(3,3) ,AX(3), B(3), YS (3) ,YPER(3) C O M M O N /O/ O D A T A LUP/15/ DEFAULT DATA DATA DATA
VALUES
OF
INPUT
((DV(I,J),J=I,3),I=I,3)/0.,0.,0.,0.05,0.,0.,0.,0.,0./
AX/I.,I.,3./,YS/O.,30.,O./,TI/O.I/,NINC/IO00/ YPER/-45.,45.,90/
S T A R T
S E S S I O N
i0 C O N T I N U E WRITE(*,11) ii F O R M A T * (///' M O V E M E N T .
e
INPUT
OF
RIGID
ELIPSOIDAL
PARTICLE
IN S L O W
*******************************************************
PARAMETERS
CALL
PARAMETERS
OF
THE
MODEL
I N P P A R ( D V , YS,AX, Y P E R , T I , N I N C , ITP,IAX, ISA)
C O M P U T E S Y M M E T R I C A N D A N T I S Y M M E T R I C P A R T OF THE V E L O C I T Y G R A D I E N T T E N S O R IN THE C O O R D I N A T E S Y S T E M F I X E D IN S P A C E DO 20 I = i , 3 DO 20 J = l , 3 EF(I,J)=(DV(I,J)+DV(J,I))*.5 OF(I,J)=(DV(I,J)-DV(J,I))*.5 20 C O N T I N U E PREPARE
VECTOR
B CHARACTERIZING
THE
PARTICLE
SHAPE
FLOWS'/
Rigid triaxial ellipsoidal particle motion
415
AI2=AX(1)*AX(1) A22=AX(2)*AX(2) A32=AX(3)*AX(3) B(1)=(A22-A32)/(A22+A32) B(2)=(A32-AI2)/(A32+AI2) B(3)=(AI2-A22)/(AI2+A22) SET PARAMETERS AND INITIAL VALUES FOR RUNGE KUTTA (SUBROUTINE X ..... I N I T I A L V A L U E OF I N D E P E N D E N T V A R I A B L E (TIME) Y(.) .. I N I T I A L V A L U E S OF E U L E R A N G L E S ( S T A R T I N G P O S I T I O N OF THE P A R T I C L E )
DVERK)
N=3 TOL=. 0001 IND=I NW=3 X=0. DO 50 I=l, 3 Y(I) =YS (I) *3. 1 4 1 5 9 2 6 / 1 8 0 . Y P R (I) = Y P E R (I) *3. 1 4 1 5 9 2 6 / 1 8 0 . 50 C O N T I N U E S W I T C H TO G R A P H I C S , P L O T THE S T E R E O G R A M B O U N D A R Y W R I T E P A R A M E T E R S OF THE M O D E L O N THE S C R E E N CALL PLOT
INTPLOT(DV,YS,AX,YPR,TI,X0,Y0,DX,DY,
INITIAL
POSITION
OF THE
PARTICLE
ISTEP=0 CALL POPLOT(Y,YPR,X0,Y0,DX,
S O L V I N G
IAX, ITP)
DY, ITP, IAX, 0,ISTEP)
J E F F E R Y'S
E Q U A T I O N S
I00 C O N T I N U E DUMMY=SETTEXTCOLOR(15) CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT('... running
RUNGE
KUTTA
DO 200
LOOP
FOR NINC
INCREMENTS
,) OF T I M E
J=I,NINC
C C
XEND
... V A L U E
OF T I M E
AT W H I C H
SOLUTION
IS D E S I R E D
C XEND=X+TI C C C
EXCLUDE
THE
SINGULAR
CASE
THETA.EQ.0.
I F ( Y ( 2 ) . N E . 0 . ) GO TO i01 DUMMY = SETTEXTCOLOR(15) CALL SETTEXTPOSITION(27,I,CURPOS) CALL OUTTEXT('Euler Angle THETA = 0 GO TO 2000 C C C C
C O M P U T E THE V E C T O R A R O U N D ITS A X E S i01
C C C C C C C
O CONTAINING
ANGULAR
')
VELOCITIES
PARTICLE
CONTINUE CALL PREP(Y,B,EF,OF)
ONE S T E P OF R U N G E K U T T A THE C A L L OF S U B R O U T I N E D V E R K F R O M I M S L L I B R A R Y . F O R THE D E T A I L E D D E S C R I P T I O N OF ITS P A R A M E T E R S , SEE M A N U A L OF I M S L (THE N A M E S OF P A R A M E T E R S H E L D F O R S I M P L I C I T Y ) . THE S U B R O U T I N E D V E R K U S E S THE F U N C T I O N F C N (LISTED HERE) W H I C H E V A L U A T E S F I R S T D E R I V A T I V E S OF E U L E R A N G L E S C A L L D V E R K ( N , FCN, X , Y , X E N D , T O L , I N D , C , N W , W , IER) IF(IND.GE.0.AND.IER.LE.0) GO TO 192 ERRONEOUS
CAGEO 20/3
OF THE
K
RETURN
FROM
RUNGE
KUTTA
416
J. JE~EK DUMMY = SETTEXTCOLOR(15) CALL SETTEXTPOSITION(27,I,CURPOS) C A L L O U T T E X T ( ' E R R O R IN R U N G E K U T T A CALL SETTEXTPOSITION(28,1,CURPOS) C A L L O U T T E X T ( ' I N D = ') W R I T E (STR6, '(I6)') I N D CALL OUTTEXT(STR6) CALL SETTEXTPOSITION(28,18,CURPOS) C A L L O U T T E X T ( ' I E R = ') W R I T E (STR6, '(I6)') IER CALL OUTTEXT(STR6) GO T O 2 0 0 0 PROJECT 102
AND
PLOT THE
RESULTING
!
,)
POINT
CONTINUE C A L L P © P L O T ( Y , Y P R , X 0 , Y 0 , D X , D Y , ITP, IAX, ISA, ISTEP)
WRITE
TIME
STEP
(STEP OF
RUNGE
KUTTA)
ON THE
SCREEN
ISTEP=ISTEP+I DUMMY=SETTEXTCOLOR(6) CALL SETTEXTPOSITION(24,18,CURPOS) W R I T E (STR6, ' (I6)') I S T E P CALL OUTTEXT(STR6) 200 C O N T I N U E C
C
DECISION
TO C O N T I N U E
OR
PRINT-SCREEN
C i000 C O N T I N U E DUMMY=SETTEXTCOLOR(15) CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT( ' CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT('CONTINUE ? (Yes/Print/No) R E A D ( * , '(A) ') A N S PRINT-SCREEN (USER S H O U L D
ON PRINTER OR TO A FILE INSERT HIS OWN SUBROUTINE
IF(ANS.NE.'p'.AND.ANS.NE.'P') CALL TOPRN(LUP) G O T O i000
GO TO
'
)
')
TOPRN) i001
C i001 C O N T I N U E IF(ANS.EQ.'n'.OR.ANS.EQ.'N') G O T O i00 C C C
EXIT FROM
THE
PROGRAM
GO TO
2000
OR RERUN MODEL WITH NEW PARAMETERS
2000 C O N T I N U E CALL SETTEXTPOSITION(29,I,CURPOS) C A L L O U T T E X T ( ' C h o o s e E to e x i t or R to r e r u n m o d e l R E A D ( * , ' (A)') A N S DUMMY=SETVIDEOMODE($DEFAULTMODE) IF(ANS.EQ.'r'.OR.ANS.EQ.'R')
GO TO
:
')
i0
STOP END C
........................................................................
SUBROUTINE
INPPAR
READS PARAMETERS LAST CALL OF THE
FROM CONSOLE (DEFAULT VALUES S U B R O U T I N E C A N BE HELD)
OR THE VALUES
OUTPUT PARAMETERS: D V ( i j) .. V E L O C I T Y G R A D I E N T T E N S O R d v ( j ) / d x ( i ) YS(I .... I N I T I A L V A L U E OF E U L E R A N G L E PHI YS(2 .... I N I T I A L V A L U E OF E U L E R A N G L E T H E T A YS(3 .... I N I T I A L V A L U E OF E U L E R A N G L E PSI AX(I .... S E M I A X I S a OF A N E L L I P S O I D A L P A R T I C L E AX(2 .... S E M I A X I S b AX(3 .... S E M I A X I S c Y P E R ..... A N G L E S C O N T R O L L I N G T H E O R I E N T A T I O N OF C O O R D I N A T E IN T H E C A S E OF P E R S P E C T I V E P L O T
FROM THE
AXES
Rigid triaxial ellipsoidal particle motion TI ....... N I N C ..... ITP IAX ISA
C
...... ...... ......
TIME INCREMENT FOR RUNGE KUTTA T H E M O T I O N O F T H E P A R T I C L E IS T R A C E D F O R N I N C I N C R E M E N T S O F TIME, A F T E R T H A T T H E P R O G R A M D E M A N D S F O R C O N T I N U A T I O N T Y P E OF P R O J E C T I O N TRACED AXIS OF THE ELLIPSOIDAL PARTICLE P A R A M E T E R T E L L I N G IF T H E a - A X I S IS T O BE S K E T C H E D A L O N G T H E P A T H OF T H E c - A X I S
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
S U B R O U T I N E I N P P A R ( D V , Y S , A X , Y P E R , T I , N I N C , I T P , IAX, ISA) DIMENSION DV(3,3),AX(3),YS(3),YPER(3) WRITE
i0
CURRENT
PARAMETERS
AND ENTER
NEW VALUES
WRITE(*,10) ((DV(I,J),J=I,3),I=I,3),AX,YS,TI,NINC FORMAT(' CURRENT PARAMETERS ARE '/ * ' Velocity gradient tensor : ',2(3F6.3/37X),3F6.3// * ' S e m i a x e s of the p a r t i c l e a , b , c : ' , 3 F 6 . 2 / * Initial Euler angles PHI,TH,PSI: ',3F6.0// * ' Time increment for Runge Kutta : ',F6.3/ * ' N u m b e r of i n c r e m e n t s : ',I6//) WRITE(*,*) * ' E N T E R N E W P A R A M E T E R S (or "," f o r k e e p i n g c u r r e n t v a l u e s ) ' WRITE(*, (A\)) V e l o c i t y g r a d i e n t t e n s o r (1-st row) : READ(*,* ( D V ( I , J ) , J = I 3) WRITE(*, (A\)) (2-nd row) : READ(*,* ( D V ( 2 , J ) , J = I 3) W R I T E ( * , (A\)') (3-rd row) : READ(*,* ( D V ( 3 , J ) , J = I 3) Particle semiaxes a,b,c : W R I T E ( * , (A\)') READ(*,*) AX Initial Euler angles PHI,TH,PSI : W R I T E ( * , ' (A\)') R E A D ( * , * ) YS Time increment for Runge Kutta : WRITE(*,'(A\) ) R E A D ( * , * ) TI N u m b e r of i n c r e m e n t s : WRITE(*,'(A\) ) READ(*,*) NINC WRITE(*,*)
READ THE TYPE OF
PROJECTION
ITP=0 20 W R I T E ( * , '(A\)') ¢ &' T Y P E OF P R O J E C T I O N (S=Stereographic,L=Lambert,P=Perspective R E A D ( * , ' ( A ) ') Q IF(Q.EQ.'S'.OR.Q.EQ.'s') ITP=2 IF(Q.EQ.'L'.OR.Q.EQ.'I') ITP=3 IF(Q.EQ.'P'.OR.Q.EQ.'p') ITP=4 I F ( I T P . E Q . 0 ) G O T O 20 IF(ITP.EQ.4) THEN WRITE(*,30) YPER 30 F O R M A T ( / ' A n g l e s c o n t r o l i n g a x e s o r i e n t a t i o n are : ',3F6.0/) WRITE(*,'(A\)') * ' E N T E R N E W V A L U E S (or ",,," ) : READ(*,*) YPER E N D IF WRITE(*,*) CHOOSE
AN AXIS
T O BE T R A C E D
(a/b/c)
IAX=0 40 W R I T E ( * ' (A\) ') ' T R A C E D A X I S (a/b/c) R E A D ( * , '(A)') Q IF(Q.EQ.'A'.OR.Q.EQ.'a') IAX=I IF(Q.EQ.'B'.OR.Q.EQ.'b') IAX=2 IF(Q.EQ.'C'.OR.Q.EQ.'c') IAX=3 I F ( I A X . E Q . 0 ) G O T O 40 SKETCH
THE
a-AXIS
ALONG
T H E P A T H OF T H E
ISA=0 IF(ITP.NE.4.AND.IAX.EQ.3) THEN W R I T E ( * , '(A\)') ' S K E T C H a - A X I S R E A D ( * , '(A)') Q IF(Q.EQ.'Y'.OR.Q.EQ.'y') ISA=I E N D IF RETURN END C
?
:
c-AXIS
(Y/N)
?
:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
417
418
J. JE.~EK
SUBROUTINE
INTPLOT
SWITCHES TO GRAPHICS, PLOTS B O U N D A R Y AND P A R A M E T E R S ON THE SCREEN
OF THE STEREOGRAM,
O U T P U T PARAMETERS: X0,Y0,DX, DY ... P O S I T I O N AND WIDTH OF THE S T E R E O G R A M (USED IN OTHER DRAWING SUBROUTINES) INPUT
PARAMETERS
(SEE SUBROUTINE
INPPAR
FOR THEIR
ON THE SCREEN
DESCRIPTION)
........................................................................
S U B R O U T I N E INTPLOT(DV, Y S , A X , Y P R , T I , X 0 , Y 0 , D X , D Y , INCLUDE 'FGRAPH.FD' INTEGER*2 D U M M Y , X M DIMENSION DV(3,3),AX(3),YS(3) DOUBLE P R E C I S I O N YPR(3) CHARACTER*30 LIST/"T'TMS RMN'H20WI6B"/ CHARACTER STR*20,STR6*6 C H A R A C T E R A X I S ( 3 ) * 2 /' A',' B',' C'/ RECORD /RCCOORD/ CURPOS RECORD /XYCOORD/ XY C C C
SET THE V I D E O M O D E
(VGA,
GRAPHICS,
640x480
PIXELS,
IAX,ITP)
16 COLORS)
DUMMY=SETVIDEOMODE($VRESI6COLOR) COLORED LOCATE
TITLE
ON THE S C R E E N
AND R E G I S T E R
FONTS
IF AVAILABLE
IN THE ACTUAL
DIRECTORY
IF(REGISTERFONTS('TMSRB.FON').LT.0) GO TO i00 IF(SETFONT(LIST).LE.0) GO TO 100 DUMMY=SETCOLOR(3) CALL MOVETO(0,0,XY) CALL O U T G T E X T ( ' E I I i p s o i d a l Particle Motion') CALL M O V E T O ( 0 , 1 8 , X Y ) CALL OUTGTEXT(' ') GO TO 200 THE CASE
TITLE
PARAMETERS:
ALL O T H E R C
WRITES
THE FONTS
WERE NOT FOUND
i00 CONTINUE D U M M Y = SETTEXTCOLOR(3) CALL S E T T E X T P O S I T I O N ( I , I , CALL O U T T E X T ( ' E l l i p s o i d a l CALL SETTEXTPOSITION(2,1, CALL OUTTEXT(' SET P O S I T I O N
AND WIDTH
CURPOS) Particle CURPOS)
Motion') ')
OF THE S T E R E O G R A M
200 C O N T I N U E X0=410. Y0=230. DX=200. DY=200. DRAW B O U N D A R Y
OF THE S T E R E O G R A M
(CIRCLE
DUMMY=SETCOLOR(2) DUMMY=ELLIPSE($GBORDER, * INT(X0-DX),INT(Y0-DY)
INT(X0+DX),INT(Y0+DY))
IF(ITP.NE.4) THEN DUMMY=SETCOLOR(15) DO XM=I,10 DUMMY=SETPIXEL(INT(X0),INT(Y0-DY)-XM) D U M M Y = S E T P I X E L ( I N T ( X 0 + D X ) + X M , INT(Y0)) END DO CALL S E T T E X T P O S I T I O N ( I , 5 2 , C U R P O S ) CALL OUTTEXT('N') CALL OUTTEXT('X') CALL S E T T E X T P O S I T I O N ( 1 5 , 8 0 , C U R P O S ) CALL OUTTEXT('E') CALL OUTTEXT('Y') END IF
Rigid triaxial ellipsoidal particle motion
DRAW THE AXES
FOR P E R S P E C T I V E
IF(ITP.EQ.4) C C C
SET THE C O L O R
419
PLOT
CALL AXPLOT(YPR, X0,Y0,DX, DY) FOR T R A C I N G
PARTICLE
DUMMY=SETCOLOR(15) C C C
WRITE
PARAMETERS
OF THE MODEL
ON THE SCREEN
DUMMY=SETTEXTCOLOR(6) CALL S E T T E X T P O S I T I O N ( 5 , I , C U R P O S ) CALL O U T T E X T ( ' V e l o c i t y Gradient Tensor') DO J = 1,3 CALL S E T T E X T P O S I T I O N ( J + 6 , I , C U R P O S ) WRITE (STR,' (3(IX,F5.3))') (DV(J,K),K=I,3) CALL OUTTEXT(STR) END DO CALL S E T T E X T P O S I T I O N ( 1 2 , 1 , C U R P O S CALL O U T T E X T ( ' A x i a l Ratio ') WRITE(STR,'(2(F4.2,3H : ),F4.2)' AX CALL S E T T E X T P O S I T I O N ( 1 4 , 3 , C U R P O S CALL OUTTEXT(STR) CALL S E T T E X T P O S I T I O N ( 1 7 , I , C U R P O S ) CALL O U T T E X T ( ' I n i t i a l Euler Angles') CALL S E T T E X T P O S I T I O N ( 1 9 , I , C U R P O S ) WRITE (STR,' (3(IX,F5.0)) ') (YS(J),J=I,3) CALL OUTTEXT(STR) C ISTEP=0 CALL S E T T E X T P O S I T I O N ( 2 2 , I , C U R P O S ) CALL O U T T E X T ( ' T i m e I n c r e m e n t :') CALL S E T T E X T P O S I T I O N ( 2 2 , 1 9 , C U R P O S ) WRITE (STR6,'(F5.3,1H)') TI CALL OUTTEXT(STR6) CALL S E T T E X T P O S I T I O N ( 2 4 , I , C U R P O S ) CALL O U T T E X T ( ' T i m e Step :') CALL S E T T E X T P O S I T I O N ( 2 4 , 1 8 , C U R P O S ) WRITE (STR6, '(I6)') ISTEP CALL OUTTEXT(STR6) C CALL S E T T E X T P O S I T I O N ( 2 6 , I , C U R P O S ) CALL O U T T E X T ( ' T r a c e d Axis : ') CALL S E T T E X T P O S I T I O N ( 2 6 , 1 5 , C U R P O S ) WRITE (STR6,'(A2)') AXIS(IAX) CALL OUTTEXT(STR6) RETURN END C
.........................................................................
C C
SUBROUTINE
AXPLOT
C C
AUXILIARY
SUBROUTINE
FOR DRAWING
AXES
IN THE P E R S P E C T I V E
C C
.........................................................................
SUBROUTINE AXPLOT(YPR,X0,Y0,DX,DY) INCLUDE 'FGRAPH.FD' INTEGER*2 DUMMY,XM, YM, XX,YY,D RECORD /XYCOORD/ XY DOUBLE P R E C I S I O N Y P R ( 3 ) , R ( 3 , 3 ) , X A ( 3 ) D I M E N S I O N XF(3) DATA D / 5 / , D D / I . 0 5 / XM=INT(X0) YM=INT(Y0) CALL RMAT(YPR,R) DO i0 I=i,3 XF(I)=R(I,I) i0 CONTINUE XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) DUMMY=SETCOLOR(10) CALL MOVETO(XM,YM, XY) D U M M Y = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL M O V E T O ( X X + D , Y Y + D , XY) D U M M Y = LINETO(XX-D, YY-D)
PLOT
420
J. JE~EK CALL M O V E T O ( X X + D , Y Y - D , X Y ) DUMMY LINETO(XX-D,YY+D) DO 20 I=i,3 XF(I)=R(I, 2) 20 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) C A L L MOVETO(XM, YM,XY) DUMMY = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL MOVETO(XX+D,YY-D,XY) DUMMY = LINETO(XX,YY) DUMMY = LINETO(XX-D,YY-D) CALL MOVETO(XX, YY, XY) D U M M Y = LINETO(XX, YY+D) DO 30 I=i,3 XF(I)=R(I,3) 30 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) CALL M O V E T O ( X M , Y M , X Y ) DUMMY = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL MOVETO(XX-D,YY-D,XY) DUMMY LINETO(XX+D,YY-D) DUMMY LINETO(XX-D,YY+D) D U M M Y = LINETO(XX+D, YY+D) =
= =
C DUMMY=SETCOLOR(II) XA(3)=0. DO 70 K = I , 3 6 0 " 4
XA(1)=COS(3.1415926/180.*FLOAT(K-I)/4.) XA(2)=SIN(3.1415926/180.*FLOAT(K-I)/4.) DO 50 I=i,3 XF(I)=0. DO 50 J=l,3 XE(I)=XF(I)+XA(J)*R(I,J) 50 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) IF(XF(1).GE.0.) DUMMY=SETPIXEL(XX,YY) IF(XF(1).LT.0.AND.(INT(FLOAT(K)/8.)*8).EQ.K) 70 C O N T I N U E RETURN END C
C C C C C C C C C C C C C
D U M M Y = S E T P I X E L ( X X , YY)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
SUBROUTINE
POPLOT
P R O J E C T S A N D PLOTS ONE POINT OF THE T R A C E D AXIS PATH TO ONE STEP OF RUNGE KUTTA)
(CORRESPONDING
INPUT PARAMETERS: Y ... V E C T O R C O N T A I N I N G E U L E R ANGLES PHI, T H E T A AND PSI D E S C R I B I N G THE I N S T A N T A N E O U S P O S I T I O N OF THE P A R T I C L E ALL O T H E R P A R A M E T E R S (FOR THEIR D E S C R I P T I O N SEE S U B R O U T I N E I N P P A R A N D INTBLOT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
S U B R O U T I N E POPLOT(Y, YPR, X0,Y0,DX,DY, ITP,IAX,ISA, ISTEP) I N C L U D E 'FGRAPH.FD' D O U B L E P R E C I S I O N Y(3),YPR(3) I N T E G E R * 2 DUMMY,XM, YM, XX,YY R E C O R D / X Y C O O R D / XY D I M E N S I O N XF(3) D A T A IDS/200/, ISK/100/ C C C
COMPUTE
COORDINATES
OF THE P L O T T E D POINT O N THE S C R E E N
CALL P R O J C T ( Y , X F , Y P R , IAX, ITP) XX=INT(XF(2)*DX+X0) YY=INT(-XF(1)*DY+Y0) GO TO (2,2,2,4), ITP PLOT POINT
IN THE S T E R E O G R A P H I C
OR L A M B E R T
PROJECTION
Rigid triaxial ellipsoidal particle motion
421
2 CONTINUE DUMMY=SETPIXEL(XX,YY) IF(ISA.EQ.0) RETURN C C
C C
SKETCH THE a-AXIS EACH ISK STEPS)
ALONG
THE
c-AXIS
PATH
(THE a - A X I S
IF(INT(FLOAT(ISTEP)/FLOAT(ISK))*ISK.EQ.ISTEP) C A L L P R O J C T ( Y , X F , Y P R , I,I) XM=INT(XF(2)*DX/10.)+XX YM=INT(-XF(1)*DY/10.)+YY C A L L M O V E T O ( X X , Y Y , XY) DUMMY = LINETO(XM,YM) CALL MOVETO(XX, YY,XY) E N D IF RETURN
IS S K E T C H E D
AFTER
THEN
C C C
THE CASE OF OF T H E P A T H
PERSPECTIVE IS P I C T U R E D
PLOT - ONLY EACH IDS-th POINT O N T H E I N V I S I B L E S I D E OF T H E S P H E R E
C 4 CONTINUE IF(XF(3).GE.0.) DUMMY=SETPIXEL(XX,YY) IF(XF(3).LT.0.AND.(INT(FLOAT(ISTEP)/FLOAT(IDS))*IDS.EQ.ISTEP)) * D U M M Y = S E T P I X E L ( X X , YY) RETURN END C
C C C C C C C C C
C C C
C C C C C
.........................................................................
SUBROUTINE
PROJCT
PROJECTS THE POINT CORRESPONDING TO THE TRACED TYPE OF PROJECTION OR IN THE PERSPECTIVE PLOT
INPUT PARAMETERS: .... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, T H E T A , PSI Y P R .. A N G L E S C O N T R O L L I N G T H E P O S I T I O N OF C O O R D I N A T E A X E S PERSPECTIVE PLOT I A X .. P A R A M E T E R T E L L I N G W H I C H P A R T I C L E A X I S IS T R A C E D I T P .. T Y P E OF P R O J E C T I O N
.........................................................................
C DIRECTIONAL
COSINES
OF T H E a - A X I S
C 1 CONTINUE XF(1)=CPH*CPS-SPH*CTH*SPS
XF(2)=SPH*CPS+CPH*CTH*SPS XF(3)=STH*SPS GO TO 5 C C
DIRECTIONAL
COSINES
OF THE
b-AXIS
C 2 CONTINUE
XF(1)=-CPH*SPS-SPH*CTH*CPS XF(2)=-SPH*SPS+CPH*CTH*CPS XF(3) = S T H * C P S G O TO 5 C C C
IN
OUTPUT PARAMETERS: XF ... C O O R D I N A T E S O F T H E P R O J E C T E D P O I N T C O R R E S P O N D I N G TO T H E T Y P E OF P R O J E C T I O N G I V E N (XF(1) A N D XF(2) A R E T H E Y - C O O R D I N A T E A N D T H E X - C O O R D I N A T E IN T H E P L O T T E D C I R C L E )
S U B R O U T I N E P R O J C T ( Y , X F , Y P R , I A X , ITP) DOUBLE PRECISION Y(3),R(3,3),YPR(3) DIMENSION XF(3),XA(3) SPH=SIN(Y(1)) CPH=COS(Y(1)) STH=SIN(Y(2)) CTH=COS(Y(2)) SPS=SIN(Y(3)) CPS=COS(Y(3)) G O TO (1,2,3), I A X C
IN T H E D E S I R E D
Y
C C
AXIS
DIRECTIONAL
COSINES
3 CONTINUE XF(1)=SPH*STH
OF T H E
c-AXIS
422
J. JE~EK XF(2)=-CPH*STH XF(3)=CTH 5 CONTINUE IF(ITP.NE.2.AND.ITP.NE.3)
GO TO
i0
L O W E R H E M I S P H E R E P R O J E C T I O N W I L L BE C O N S I D E R E D IN THE S T E R E O G R A P H I C OR L A M B E R T P R O J E C T I O N ( X F ( 3 ) . L E . 0 ) IF(XF(3).LE.0.) XF(1)=-XF(1) XF(2)=-XF(2) XF(3)=-XF(3) RECOMPUTE i0 GO TO
GO TO
COORDINATES
ORTHOGONAL
OF
I0
XF F O R THE
(11,12,13,14),
CASE
PARTIAL
TYPE
OF P R O J E C T I O N
ITP
PROJECTION
ii C O N T I N U E RETURN STEREOGRAPHIC
PROJECTION
(WULF)
12 C O N T I N U E XF(1)=XF(1)/(I.-XF(3)) XF(2)=XF(2)/(I.-XF(3)) RETURN EQUAL-AREA
PROJECTION
(LAMBERT)
13 C O N T I N U E QQ=XF(1)*XF(1)+XF(2)*XF(2) IF(QQoEQ.0.) RETURN Q=SQRT(QQ+(XF(3)+I)*(XF(3)+I)) QQ=SQRT(QQ) Q=Q/(QQ*SQRT(2.)) XF(1)=XF(1)*Q XF(2)=XF(2)*Q RETURN PERSPECTIVE
PLOT
14 C O N T I N U E DO 17 I=1,3 XA(I)=XF(I) 17 C O N T I N U E C A L L R M A T ( Y P R , R) DO 20 I=i,3 XF(I)=0. DO 20 J = l , 3 XF(I)=XF(I)+XA(J)*R(I,J) 20 C O N T I N U E XA(1)=XF(1) XF(1)=XF(3) XF (3)=XA(1) RETURN END C
......................................................................
SUBROUTINE
PREP
COMPUTES VECTOR A R O U N D ITS A X E S
O C O N T A I N I N G A N G U L A R V E L O C I T I E S OF THE P A R T I C L E IN THE C O O R D I N A T E S Y S T E M F I X E D IN THE P A R T I C L E
INPUT PARAMETERS: Y ....... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, THETA, PSI B ....... V E C T O R C H A R A C T E R I Z I N G THE S H A P E OF THE P A R T I C L E E F , O F ... S Y M M E T R I C A N D A N T I S Y M M E T R I C PART, R E S P E C T I V E L Y , OF THE V E L O C I T Y G R A D I E N T T E N S O R IN THE C O O R D I N A T E S Y S T E M F I X E D IN T H E S P A C E OUTPUT PARAMETERS: O ... V E C T O R C O N T A I N I N G A N G U L A R V E L O C I T I E S A R O U N D ITS A X E S ( O U T P U T V I A COMMON)
OF THE
PARTICLE
C ......................................................................
423
Rigid triaxial ellipsoidal particle motion SUBROUTINE PREP(Y,B,EF,OF) DOUBLE PRECISION Y(3),O(3),EB(3,3),OB(3,3),R(3,3) DIMENSION EF(3,3),OF(3,3),B(3) C O M M O N /0/ 0 COMPUTE CALL
THE
ROTATION
MATRIX
RMAT(Y,R)
C O M P U T E S Y M M E T R I C A N D A N T I S Y M M E T R I C P A R T OF THE V E L O C I T Y T E N S O R IN T H E R O T A T I N G C O O R D I N A T E S Y S T E M
GRADIENT
DO DO
20 I = i , 3 20 J = l , 3 OB(I,J)=.0 EB(I,J)=.0 DO 15 K = I , 3 DO 15 L = I , 3 RR=R(I,K)*R(J,L) EB(I,J)=EB(I,J)+RR*EF(K,L) OB(I,J)=OB(I,J)+RR*OF(K,L) 15 CONTINUE 20 C O N T I N U E EVALUATE
VECTOR
0
O(1)=B(1)*EB(3,2)-OB(3,2) O(2)=B(2)*EB(I,3)-OB(I,3) O(3)=B(3)*EB(2,1)-OB(2,1) RETURN END C
.........................................................................
C C
SUBROUTINE
RMAT
C C
COMPUTES
THE
ROTATION
MATRIX
C C C
C C C C C
INPUT PARAMETERS: ... V E C T O R C O N T A I N I N G
Y
EULER
ANGLES
PHI,
THETA,
PSI
OUTPUT PARAMETERS: R ... R O T A T I O N M A T R I X .........................................................................
SUBROUTINE RMAT(Y,R) DOUBLE PRECISION Y(3),R(3,3) S P H = S I N ( Y i) C P H = C O S ( Y i) S T H = S I N ( Y 2) C T H = C O S ( Y 2) S P S = S I N ( Y 3) C P S = C O S ( Y 3) R i,i = C P H * C P S - S P H * C T H * S P S R 1,2 = S P H * C P S + C P H * C T H * S P S R 1,3 = S T H * S P S R 2,1 = - C P H * S P S - S P H * C T H * C P S R 2,2 = - S P H * S P S + C P H * C T H * C P S R 2,3 = S T H * C P S R 3,1 = SPH*STH R 3,2 = - C P H * S T H R 3 , 3 = CTH RETURN END C
.........................................................................
SUBROUTINE COMPUTES
FCN
VALUES
OF F I R S T
DERIVATIVES
OF E U L E R
ANGLES
INPUT PARAMETERS: N ... L E N G T H OF THE V E C T O R Y Y ... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, THETA, PSI X ... N O T U S E D (HELD F O R C O M P A T I B I L I T Y W I T H S U B R O U T I N E DVERK) O ... V E C T O R C O N T A I N I N G A N G U L A R V E L O C I T I E S OF THE P A R T I C L E A R O U N D ITS A X E S (INPUT V I A C O M M O N ) OUTPUT YP ...
PARAMETERS: FIRST DERIVATIVES
OF E U L E R
ANGLES
424
J. JE~EK
SUBROUTINE FCN(N,X,Y, YP) DOUBLE PRECISION X,Y(N),YP(N),O(3) COMMON /O/ O YP(1)=(O(1)*SIN(Y(3))+O(2)*COS(Y(3)))/SIN(Y(2)) YP(2)=O(1)*COS(Y(3))-O(2)*SIN(Y(3)) YP(3)=O(3)-YP(1)*COS(Y(2)) RETURN END
Pergamon
0098-3004(93)E0017-A
Copyright © 1994ElsevierScienceLtd Printed in Great Britain.All rights reserved 0098-3004/94$7.00+ 0.00
S O F T W A R E FOR M O D E L I N G THE M O T I O N OF R I G I D T R I A X I A L ELLIPSOIDAL PARTICLES IN VISCOUS FLOW JOSEF JE~,EK Department of Applied Mathematics and Computer Science, Faculty of Science, Charles University, Albertov 6, 128 43 Praha 2, Czech Republic (Received 31 October 1993)
Abstract--Three MS-FORTRAN programs are presented which enable study of the behavior of rigid triaxial ellipsoidal particles in slowly moving viscous flow. The modeling is based on numerical solution of Jeffery's equations. These programs allow the investigation of individual particle motion as well as the evolution of a multiparticle system composed of noninteracting particles. In the combined pure shear and simple shear flow, it is possible to estimate the character of multiparticle system evolution by diagrams of time derivatives of Euler angles. The programs are interactive with graphical output on the screen and run on a standard IBM-PC with VGA. Some examples of graphical output of the modeling are presented. Key Words: Ellipsoidal particle motion, Viscous flow, Multiparticle system motion, Pure shear, Simple
shear
INTRODUCTION In the interpretation of preferred orientations of fabric elements in tectonites the model of Jeffery (1923) is used. Jeffery's model describes the behavior of a rigid ellipsoidal particle in slow viscous flow. Given the flow geometry (velocity gradient tensor), the particle axial ratio, and its starting position, the trajectory of the particle during the flow can be determined as the solution of Jeffery's equations. Analytical solutions of these equations of the particle motion have been determined only for biaxial particles (i.e. ellipsoids of revolution) and restricted types of flow geometry: simple shear flow (Jeffery, 1923), pure shear flow (Gay, 1968) axial flattening (Debat and others, 1975; Tulis, 1976), and general orthorhombic flow (Fernandez, 1988). These analytical solutions were used in modeling the motion of an individual biaxial particle and also in modeling multiparticle systems composed of noninteracting particles (Reed and Tryggvason, 1974; Tullis, 1976). In the situation of a triaxial ellipsoidal particle, numerical techniques can be used for solving Jeffery's equations as in Hinch and Leal (1979). Freeman (1985), and Fernandez and Laporte (1991). The results of investigating the motion of an individual triaxial particle with different initial orientations in different types of flow show a significant deviation of the particle behavior from that of a biaxial particle and the need of further work in this field. Modeling of a multiparticle system composed of noninteracting rind ellipsoidal triaxial particles has been also suggested by Freeman (1985), but such a work has not appeared up to this time.
The aim of this paper is the presentation of software enabling the user to investigate the motion of an individual triaxial ellipsoidal particle as well as the evolution of a multiparticle system composed of noninteracting triaxial particles, in different flow geometries. The programs presented could serve also as a tool for estimating the possible development of preferred orientation in real conditions. They were written in MS-FORTRAN, v. 5.0, and run on a standard IBM PC. The Intel 80486 or mathematical coprocessor is recommended for satisfactory speed. The programs are interactive, demand necessary input parameters from console, and yield graphical output of results on the screen. Program JEFF listed in the Appendix enables the modeling of an individual particle motion. The motion of the particle is drawn on the screen as the path of one of the particle axes in the selected type of projection. For solving the differential equations of the particle motion, subroutine DVERK from the IMSL Library, performing the fourth- and fifth-order Runge-Kutta method, was used. This subroutine is not listed here. The user also should provide his/her own subroutine performing print-screen on his type of printer. Program MUPLES allows the user to model and visualize a multiparticle system evolution. Because of its affinity to the program JEFF and the finite extent of the paper it is described only briefly and not listed, and the program TIDER, drawmg patterns of Euler angles derivatives, as well. These diagrams help in estimating the character of multiparticle system evolution, especially in the situation of combined pure shear and simple shear model,
409
410
J. JE~EK
which can be used as an approximation in many real situations. THE MODEL OF ,IEFFERY
The equations describing our model were given by Jeffery (1923). In his work, he considered a slowly moving Newtonian viscous fluid in which sufficiently small rigid and noninteracting particles of an ellipsoidal shape are immersed. Starting with the Navier-Stokes equations of motion of incompressible viscous fluid and neglecting the terms of the second order, Jeffery determined the equations of motion of the individual particle. Consider a right-handed orthogonal coordinate system (x ~, x ~, x ~) positioned at the particle mass center with axes in the directions fixed in space (further termed "fixed coordinate system") and let another coordinate system (x~, x2, x3) be coincident with the axes a, b, c of the ellipsoidal particle (further termed "rotating coordinate system"). Denote the velocity gradient tensor characterizing undisturbed flow at a large distance from the particle by V ~ = Ov~lt3x~, and suppose that V~ is constant in time and obeys the equation of continuity. Denote further by E~ and f ~ the symmetric and antisymmetric part of the velocity gradient tensor, respectively. Then, in the rotating coordinate system, the Jeffery equations can be written as ~01 = Bt E32 - ~'~32, 092 = B2E13
-
-
~'~13'
(1)
These equations can be solved numerically as an initial value problem. Given the velocity gradient tensor, axial ratio of the particle, and initial values of the Euler angles 4~, 0, ~ (the initial position of the particle), we compute the time evolution of these angles, that is the particle rotation.
INDIVIDUAL PARTICLE MOTION (Program JEFF)
For the numerical solution of the Equations (2) and display of the resulting particle motion, program J E F F was written. The program demands the following input parameters from console: --velocity gradient tensor V ~, --semiaxes of the ellipsoidal particle a, b, c, --initial Euler angles ~b, 0, ~, of the particle, - - t i m e increment for the Runge-Kutta method, - - n u m b e r of time steps performed at once (before the program again demands for continuation), --which of the particle axes will be traced (its path will be drawn on the screen), - - t y p e of projection for the display of this particle path (stereographic, Lambert, or perspective p l o ~ s e e Figs. 1 and 2), --optionally, if the c-axis is traced, the position of the a-axis can be sketched along the c-axis path, which allows to imagine the triaxial particle rotation around its c-axis.
093 = B3 E21 - ~21, where eg~, i = 1,2,3, are the angular velocities of the particle around its axes and the vector B l = (b 2 _ c 2)/(b 2 + c 2),
B2 = (c 2 - a 2)/(c 2 + a 2), B 3 = (a 2 _ b 2)/(a 2 + b 2)
corresponds to the axial ratio of the ellipsoidal particle. An analytical solution of Jeffery's equations has been determined only for biaxial particle in some flow geometry as was mentioned in the Introduction. These equations, however, can be reformulated to the form convenient for numerical solution. Describing the rotation of the particle with respect to the fixed coordinate system by Euler angles q% 0, ~b and using the known relations between the angular velocity ~0 and the time derivatives of the Euler angles we obtain three differential equations corresponding to (1) 4; = (~ol sin @ + ~o2cos ~, )/sin 0, = ol cos ~ - o~2 sin $, = ~o3 -
4; c o s 0.
(2)
The program is interactive. The path of a traced particle is drawn on the screen for the given number of time increments of the Runge-Kutta method and after that the program asks if continuation or printscreen is desired. If no continuation is selected, the program offers to rerun the model or exit. If rerun is desired, the user can keep some parameters from the preceding run. For solving the differential Equations (2), subroutine DVERK of the IMSL Library, performing the fourth- and fifth-order Runge-Kutta method, was precise by a comparison of its results with some solutions known analytically. The recommended value of the time step for the Runge-Kutta method is 0.1 Subroutine DVERK is not included in the listing of program J E F F (see Appendix) and neither is subroutine TOPRN for print-screen. The user should provide his/her own subroutine performing print-screen on his/her type of printer. Standard VGA is required in subroutines for graphical output on the screen. Examples of use of the program J E F F are in Figures 1 and 2. In Figure 1 one of classical Jeffery's orbits of a biaxial particle is drawn in perspective plot, representing the fact that biaxial particles in simple shear flow follow closed orbits and their
Rigid triaxial ellipsoidal particle motion
411
Ellipsoidal Particle Motion
Uelocity Gradient .000 .050 .000 Axial
Tensor
.000 .000 .OOO .OOO .000 .000 Ratio
1.00 : 1.00 : 3 . 0 0 Euler Rngles
Initial
O.
30.
O.
Time Increment :
.100
Time S t e p
SO00
:
Traced R x i s
:
C
/~//' ,
"~-'~-~--'--~-~
///// /
Figure 1. Jeffery's orbit of biaxial particle in perspective plot. See description of program JEFF for explanation of parameters of model.
motion is periodic. In Figure 2 a triaxial particle was traced. In the Lambert equal-area projection used, the a-axis has been sketched along the c-axis path, which allows to imagine the rotation of the particle around its c-axis. The same simple shear flow as in Figure 1 was used, but the triaxial particle motion is no longer periodic.
Ellipsoidal
Particle
MULTIPARTICLE SYSTEM EVOLUTION
(Program MUPLES) Supposing that there are no interactions among the individual particles rotating in viscous fluid, the multiparticle system can be investigated by Equations (2). For this purpose program M U P L E S was written.
X
Motion
I
Uelocity .000 .050 .000 Axial
Gradient
Tensor
~
.000 .OOO .OOO . ~ ,090 .000 Ratio
.67
Initial G.
: 1.50
: 3.00
¥
Failer nmjles 30.
0.
Time Increment :
.IO0
Time Step :
SO00
Traced
7
/
Axis : C
Figure 2. Behavior of triaxial particle in same flow as in Figure 1. In Lambert equal-area projection used, a-axis has been sketched along c-axis path.
412
J. JE~EK
Four - Groups System Motion
x
U e l o c i t y Gradient Tensor .010 .000 .050 - . 0 1 0 .000 .000
.GO0 .OO0 .000
o
[] \
o 0
1
:
2
:
3
o
1
:
Z
:
4
,~
1
:
4
:
6
[]
1
:
Z
:
6
STEP :
/
X
o
[]
\\o ~
A
o ~o ~>o o A
<6 no ~ ~
[]
\
0
~
o
o
6
1
/ /,'
100
Figure 3. Modeling of multiparticle system evolution by program MUPLES. Multiparticle system is composed of four groups of particles, each of them having about 100 particles of same axial ratio. Current positions of c-axes of particles are marked by symbols in Lambert equal-area projection. Different symbols correspond to different axial ratios. Because of the finite extent of this paper, the program MUPLES is not listed and only a brief description will be given here. The program can be created using some subroutines of program JEFF. The program input parameters are the velocity gradient tensor, time increment for the Runge-Kutta method, number of time steps performed at once, particle axis to be traced and type of projection. Initial Euler angles of the multiparticle system can be read from an ASCII file or randomly generated. A convenient selection is an initial is•tropic distribution of angles tk, 0 (similar to Sanderson and Meneilly, 1981) and a random distribution of the angle ~. During the time evolution of the multiparticle system studied, after the number of Runge-Kutta steps computed at once, the selected axis positions for all the particles are marked in the desired type of projection. Optionally, the orientation tensor (for definition, see e.g. Sanderson and Meneilly, 1981) and its eigenvalues and eigenvectors or intensity of the fabric (Lisle, 1985) can be computed and presented. In Figure 3 an example of output from program MUPLES is shown. A multiparticle system was composed of four groups of particles, each of them having about 100 particles of the same axial ratio. The system started with an initial is•tropic distribution of angles ~, 0 and a random distribution of the angle ~, which determine the position of the particle. EVALUATING
THE EULER ANGLES DERIVATIVES
TIME
(Program TIDER) The mentioned programs JEFF and MUPLES
compute and visualize the time evolution of a oneor multiparticle system. Program TIDER allows another type of view on the system development based on the fact that this model is stationary. The velocity gradient tensor V~. is constant in the fixed coordinate system and thus for particles of a given axial ratio the right-hand sides of the Equations (2) depend only on the Euler angles, that is on the instantaneous position of particles. The analysis of this dependence can help estimate the future or final state of the one- or multiparticle system. This is valid, for example in the combined pure shear and simple shear model given by the velocity gradient tensor
v~=
~
-~
0
0
1
(3)
where ~ and 7 are constant in time. Such a model may approximate many real geological situations. Therefore program TIDER was written.to p!ot the distribution patterns of the derivatives ~, 0, ~ given by the Equations (2) in an equal-area stere•graphic projection. Input parameters are the velocity gradient tensor V ~and the lengths of the particle semiaxes a, b, c. For the possibility of efficient print-screen the dithering technique was used in plotting these diagrams. Because of the finite extent of this paper, program TIDER is not listed here. The program uses some subroutines of program JEFF. As an example of the usefulness of program TIDER, Figure 4 shows a ~b-derivative pattern corresponding to a combined pure shear and simple shear flow. The
Rigid triaxial ellipsoidal particle motion
dPH/dt
Field
Oelocitg 6 r a d i e n t .010 .000 .050 -.el0 •0 0 0 .000 Axial
413
Tensor
.0~ .000 .000
Batio .33
: 1.50
: 6.00
dPH/dt
< -Z.O00 < -1.75o < -1.5oo <
< < <
< <
- 1 .zso -1.000
-no -.500 -.250 . O~
•
I M
.....
Figure 4. ~b-Derivative pattern corresponding to combined pure shear and simple shear flow. In gray-scale used white corresponds to nonnegative values of ~b and all gray grades correspond to negative values of
9.
values of ~b-derivative were computed for an isotropic distribution of angles 4), 0 and a random distribution of the angle ~k, and then plotted in the stereogram. In the gray-scale used the white corresponds to nonnegative values of ~b and gray grades correspond to negative values of qk The ~b-derivative pattern can be considered as representation of a multiparticle system composed of triaxial particles of the same axial ratio, with the position of particles given by the angles 4), 0, and ~b. The value of ~b has a strong influence on the next motion of the particles. In the next evolution of the system, the particles having nonnegative values of ~b will remain in their position or will rotate in such a way that their c-axis projection will move clockwise. On the other side, the projection of the c-axis of the particles having negative values of ~b will move counterclockwise. Ana!yzing simultaneously the patterns of the derivatives q~, 0, we can deduce the type of evolution of a given system. In the given example a stable preferred orientation will arise during the multiparticle system evolution. CONCLUSION The presented programs J E F F , MUPLES, and T I D E R can be used for a detailed discussion o f the one- and multiparticle system evolution in different flow geometries. The possibility of the preferred orientation development in the combined pure shear and simple shear flow was investigated by
Je2ek and others (1994). As a natural example to which the Jeffery's model can be applied, ultramylonites were mentioned. In the ultramylonites the matrix/clast proportion may be high. Thus, it is possible to investigate the rotation of relict clasts in the time of the development of fabric using the Jeffery's model.
Acknowledgments--I would like to thank Karel Schulmann who initiated this work. He and Radek Melka introduced me into this problem and had clear requirements on the capabilities of programs. Also, I thank Karel Segeth who kindly reviewed the paper.
REFERENCES Debat, P., Sirieys, P., Deramond, J., and Soula, J. C., 1975, Paleodrformations d'un massif orthogenessique (massif des Cammazes. Montagne Noire Occidentale. France): Tectonophysics, v. 28, no. 3, p. 159-183. Fernandez, A., 1988, Strain analysis from shape preferred orientation in magmatic flows: Bull. Geol. Inst. Univ. Uppsala, v. 14, p. 61~57. Fernandez, A., and Laporte, D., 1991, Significance of low symmetry fabrics in magrnatic flows: Jour. Struct. Geology, v. 13, no. 3, p. 337 347. Freeman, B., 1985, The motion of rigid ellipsoidal particles in slow flows: Tectonophysics, v. 113, no. 2, p. 163-183. Gay, N. C., 1968, The motion of rigid particles embedded in a viscous fluid during pure shear deformation of the fluid: Tectonophysics, v. 5, no. 2, p. 81-88. Hinch, E. J., and Leal, L. G., 1979, Rotation of small non-axisymetric particles in a simple shear flow: Jour. Fluid Mech., v. 92, pt. 3, p. 591~08.
414
J. JE~EK
Jeffery, G. B., 1923, The motion of ellipsoidal particles immersed in a viscous fluid: Proc. Roy. Soc. London, Ser. A, v. 102, p. 161 179. Je~ek, J., Melka, R., Schulmann, K., and Venera, Z., 1994, The behaviour of rigid triaxial ellipsoidal particles in viscous flows-modeling of fabric evolution in multiparticle system: Tectonophysics, v. 229, p. 165-180. Lisle, R. J., 1985, The use of the orientation tensor for the description and statistical testing of fabrics: Jour. Struct. Geology, v. 7, no. 1, p. 115-117.
Reed, L. J., and Tryggvason, E., 1974, Preferred orientations of rigid particles in a viscous matrix deformed by pure shear and simple shear: Tectonophysics, v. 24, no. 2, p. 85-98. Sanderson, D. J., and Meneilly, A. W., 1981, Analysis of three-dimensional strain modified uniform distributions: andalusite fabrics from a granite aureole: Jour. Struct. Geology, v. 3, no. 2, p. 109-116. Tullis, T. E., 1976, Experiments on the origin of slatey cleavage and schistosity: Geol. Soc. America Bull., v. 87, no. 5, p. 745-753.
APPENDIX Program Listing C C C C C C C C C C C C C C
P R O G R A M
JEFF
N U M E R I C A L S O L U T I O N OF J E F F E R Y ' S E Q U A T I O N S F O R T H E M O T I O N OF A R I G I D E L L I P S O I D A L P A R T I C L E IN S L O W F L O W S R U N S O N I B M - T Y P E PC W I T H VGA. M A T H E M A T I C A L C O P R O C E S S O R OR I N T E L 80486 IS R E C O M M E N D E D F O R S A T I S F A C T O R Y S P E E D INPUT
PARAMETERS:
SEE
SUBROUTINE
INPPAR
P R O G R A M U S E S S U B R O U T I N E D V E R K F R O M I M S L (NOT L I S T E D HERE) U S E R S H O U L D I N C L U D E A S U B R O U T I N E P E R F O R M I N G P R I N T - S C R E E N O N HIS O W N T Y P E OF P R I N T E R
***************************************************************************
INCLUDE 'FGRAPH.FI' INCLUDE 'FGRAPH.FD' INTEGER* 2 DUMMY C H A R A C T E R ANS, S T R 6 * 6 DOUBLE PRECISION X,Y(3),XEND,TOL,YPR(3),C(24),W(3,9),0(3) RECORD /RCCOORD/ CURPOS EXTERNAL FCN D I M E N S I O N D V ( 3 , 3 ) ,EF (3,3) ,OF(3,3) ,AX(3), B(3), YS (3) ,YPER(3) C O M M O N /O/ O D A T A LUP/15/ DEFAULT DATA DATA DATA
VALUES
OF
INPUT
((DV(I,J),J=I,3),I=I,3)/0.,0.,0.,0.05,0.,0.,0.,0.,0./
AX/I.,I.,3./,YS/O.,30.,O./,TI/O.I/,NINC/IO00/ YPER/-45.,45.,90/
S T A R T
S E S S I O N
i0 C O N T I N U E WRITE(*,11) ii F O R M A T * (///' M O V E M E N T .
e
INPUT
OF
RIGID
ELIPSOIDAL
PARTICLE
IN S L O W
*******************************************************
PARAMETERS
CALL
PARAMETERS
OF
THE
MODEL
I N P P A R ( D V , YS,AX, Y P E R , T I , N I N C , ITP,IAX, ISA)
C O M P U T E S Y M M E T R I C A N D A N T I S Y M M E T R I C P A R T OF THE V E L O C I T Y G R A D I E N T T E N S O R IN THE C O O R D I N A T E S Y S T E M F I X E D IN S P A C E DO 20 I = i , 3 DO 20 J = l , 3 EF(I,J)=(DV(I,J)+DV(J,I))*.5 OF(I,J)=(DV(I,J)-DV(J,I))*.5 20 C O N T I N U E PREPARE
VECTOR
B CHARACTERIZING
THE
PARTICLE
SHAPE
FLOWS'/
Rigid triaxial ellipsoidal particle motion
415
AI2=AX(1)*AX(1) A22=AX(2)*AX(2) A32=AX(3)*AX(3) B(1)=(A22-A32)/(A22+A32) B(2)=(A32-AI2)/(A32+AI2) B(3)=(AI2-A22)/(AI2+A22) SET PARAMETERS AND INITIAL VALUES FOR RUNGE KUTTA (SUBROUTINE X ..... I N I T I A L V A L U E OF I N D E P E N D E N T V A R I A B L E (TIME) Y(.) .. I N I T I A L V A L U E S OF E U L E R A N G L E S ( S T A R T I N G P O S I T I O N OF THE P A R T I C L E )
DVERK)
N=3 TOL=. 0001 IND=I NW=3 X=0. DO 50 I=l, 3 Y(I) =YS (I) *3. 1 4 1 5 9 2 6 / 1 8 0 . Y P R (I) = Y P E R (I) *3. 1 4 1 5 9 2 6 / 1 8 0 . 50 C O N T I N U E S W I T C H TO G R A P H I C S , P L O T THE S T E R E O G R A M B O U N D A R Y W R I T E P A R A M E T E R S OF THE M O D E L O N THE S C R E E N CALL PLOT
INTPLOT(DV,YS,AX,YPR,TI,X0,Y0,DX,DY,
INITIAL
POSITION
OF THE
PARTICLE
ISTEP=0 CALL POPLOT(Y,YPR,X0,Y0,DX,
S O L V I N G
IAX, ITP)
DY, ITP, IAX, 0,ISTEP)
J E F F E R Y'S
E Q U A T I O N S
I00 C O N T I N U E DUMMY=SETTEXTCOLOR(15) CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT('... running
RUNGE
KUTTA
DO 200
LOOP
FOR NINC
INCREMENTS
,) OF T I M E
J=I,NINC
C C
XEND
... V A L U E
OF T I M E
AT W H I C H
SOLUTION
IS D E S I R E D
C XEND=X+TI C C C
EXCLUDE
THE
SINGULAR
CASE
THETA.EQ.0.
I F ( Y ( 2 ) . N E . 0 . ) GO TO i01 DUMMY = SETTEXTCOLOR(15) CALL SETTEXTPOSITION(27,I,CURPOS) CALL OUTTEXT('Euler Angle THETA = 0 GO TO 2000 C C C C
C O M P U T E THE V E C T O R A R O U N D ITS A X E S i01
C C C C C C C
O CONTAINING
ANGULAR
')
VELOCITIES
PARTICLE
CONTINUE CALL PREP(Y,B,EF,OF)
ONE S T E P OF R U N G E K U T T A THE C A L L OF S U B R O U T I N E D V E R K F R O M I M S L L I B R A R Y . F O R THE D E T A I L E D D E S C R I P T I O N OF ITS P A R A M E T E R S , SEE M A N U A L OF I M S L (THE N A M E S OF P A R A M E T E R S H E L D F O R S I M P L I C I T Y ) . THE S U B R O U T I N E D V E R K U S E S THE F U N C T I O N F C N (LISTED HERE) W H I C H E V A L U A T E S F I R S T D E R I V A T I V E S OF E U L E R A N G L E S C A L L D V E R K ( N , FCN, X , Y , X E N D , T O L , I N D , C , N W , W , IER) IF(IND.GE.0.AND.IER.LE.0) GO TO 192 ERRONEOUS
CAGEO 20/3
OF THE
K
RETURN
FROM
RUNGE
KUTTA
416
J. JE~EK DUMMY = SETTEXTCOLOR(15) CALL SETTEXTPOSITION(27,I,CURPOS) C A L L O U T T E X T ( ' E R R O R IN R U N G E K U T T A CALL SETTEXTPOSITION(28,1,CURPOS) C A L L O U T T E X T ( ' I N D = ') W R I T E (STR6, '(I6)') I N D CALL OUTTEXT(STR6) CALL SETTEXTPOSITION(28,18,CURPOS) C A L L O U T T E X T ( ' I E R = ') W R I T E (STR6, '(I6)') IER CALL OUTTEXT(STR6) GO T O 2 0 0 0 PROJECT 102
AND
PLOT THE
RESULTING
!
,)
POINT
CONTINUE C A L L P © P L O T ( Y , Y P R , X 0 , Y 0 , D X , D Y , ITP, IAX, ISA, ISTEP)
WRITE
TIME
STEP
(STEP OF
RUNGE
KUTTA)
ON THE
SCREEN
ISTEP=ISTEP+I DUMMY=SETTEXTCOLOR(6) CALL SETTEXTPOSITION(24,18,CURPOS) W R I T E (STR6, ' (I6)') I S T E P CALL OUTTEXT(STR6) 200 C O N T I N U E C
C
DECISION
TO C O N T I N U E
OR
PRINT-SCREEN
C i000 C O N T I N U E DUMMY=SETTEXTCOLOR(15) CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT( ' CALL SETTEXTPOSITION(29,I,CURPOS) CALL OUTTEXT('CONTINUE ? (Yes/Print/No) R E A D ( * , '(A) ') A N S PRINT-SCREEN (USER S H O U L D
ON PRINTER OR TO A FILE INSERT HIS OWN SUBROUTINE
IF(ANS.NE.'p'.AND.ANS.NE.'P') CALL TOPRN(LUP) G O T O i000
GO TO
'
)
')
TOPRN) i001
C i001 C O N T I N U E IF(ANS.EQ.'n'.OR.ANS.EQ.'N') G O T O i00 C C C
EXIT FROM
THE
PROGRAM
GO TO
2000
OR RERUN MODEL WITH NEW PARAMETERS
2000 C O N T I N U E CALL SETTEXTPOSITION(29,I,CURPOS) C A L L O U T T E X T ( ' C h o o s e E to e x i t or R to r e r u n m o d e l R E A D ( * , ' (A)') A N S DUMMY=SETVIDEOMODE($DEFAULTMODE) IF(ANS.EQ.'r'.OR.ANS.EQ.'R')
GO TO
:
')
i0
STOP END C
........................................................................
SUBROUTINE
INPPAR
READS PARAMETERS LAST CALL OF THE
FROM CONSOLE (DEFAULT VALUES S U B R O U T I N E C A N BE HELD)
OR THE VALUES
OUTPUT PARAMETERS: D V ( i j) .. V E L O C I T Y G R A D I E N T T E N S O R d v ( j ) / d x ( i ) YS(I .... I N I T I A L V A L U E OF E U L E R A N G L E PHI YS(2 .... I N I T I A L V A L U E OF E U L E R A N G L E T H E T A YS(3 .... I N I T I A L V A L U E OF E U L E R A N G L E PSI AX(I .... S E M I A X I S a OF A N E L L I P S O I D A L P A R T I C L E AX(2 .... S E M I A X I S b AX(3 .... S E M I A X I S c Y P E R ..... A N G L E S C O N T R O L L I N G T H E O R I E N T A T I O N OF C O O R D I N A T E IN T H E C A S E OF P E R S P E C T I V E P L O T
FROM THE
AXES
Rigid triaxial ellipsoidal particle motion TI ....... N I N C ..... ITP IAX ISA
C
...... ...... ......
TIME INCREMENT FOR RUNGE KUTTA T H E M O T I O N O F T H E P A R T I C L E IS T R A C E D F O R N I N C I N C R E M E N T S O F TIME, A F T E R T H A T T H E P R O G R A M D E M A N D S F O R C O N T I N U A T I O N T Y P E OF P R O J E C T I O N TRACED AXIS OF THE ELLIPSOIDAL PARTICLE P A R A M E T E R T E L L I N G IF T H E a - A X I S IS T O BE S K E T C H E D A L O N G T H E P A T H OF T H E c - A X I S
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
S U B R O U T I N E I N P P A R ( D V , Y S , A X , Y P E R , T I , N I N C , I T P , IAX, ISA) DIMENSION DV(3,3),AX(3),YS(3),YPER(3) WRITE
i0
CURRENT
PARAMETERS
AND ENTER
NEW VALUES
WRITE(*,10) ((DV(I,J),J=I,3),I=I,3),AX,YS,TI,NINC FORMAT(' CURRENT PARAMETERS ARE '/ * ' Velocity gradient tensor : ',2(3F6.3/37X),3F6.3// * ' S e m i a x e s of the p a r t i c l e a , b , c : ' , 3 F 6 . 2 / * Initial Euler angles PHI,TH,PSI: ',3F6.0// * ' Time increment for Runge Kutta : ',F6.3/ * ' N u m b e r of i n c r e m e n t s : ',I6//) WRITE(*,*) * ' E N T E R N E W P A R A M E T E R S (or "," f o r k e e p i n g c u r r e n t v a l u e s ) ' WRITE(*, (A\)) V e l o c i t y g r a d i e n t t e n s o r (1-st row) : READ(*,* ( D V ( I , J ) , J = I 3) WRITE(*, (A\)) (2-nd row) : READ(*,* ( D V ( 2 , J ) , J = I 3) W R I T E ( * , (A\)') (3-rd row) : READ(*,* ( D V ( 3 , J ) , J = I 3) Particle semiaxes a,b,c : W R I T E ( * , (A\)') READ(*,*) AX Initial Euler angles PHI,TH,PSI : W R I T E ( * , ' (A\)') R E A D ( * , * ) YS Time increment for Runge Kutta : WRITE(*,'(A\) ) R E A D ( * , * ) TI N u m b e r of i n c r e m e n t s : WRITE(*,'(A\) ) READ(*,*) NINC WRITE(*,*)
READ THE TYPE OF
PROJECTION
ITP=0 20 W R I T E ( * , '(A\)') ¢ &' T Y P E OF P R O J E C T I O N (S=Stereographic,L=Lambert,P=Perspective R E A D ( * , ' ( A ) ') Q IF(Q.EQ.'S'.OR.Q.EQ.'s') ITP=2 IF(Q.EQ.'L'.OR.Q.EQ.'I') ITP=3 IF(Q.EQ.'P'.OR.Q.EQ.'p') ITP=4 I F ( I T P . E Q . 0 ) G O T O 20 IF(ITP.EQ.4) THEN WRITE(*,30) YPER 30 F O R M A T ( / ' A n g l e s c o n t r o l i n g a x e s o r i e n t a t i o n are : ',3F6.0/) WRITE(*,'(A\)') * ' E N T E R N E W V A L U E S (or ",,," ) : READ(*,*) YPER E N D IF WRITE(*,*) CHOOSE
AN AXIS
T O BE T R A C E D
(a/b/c)
IAX=0 40 W R I T E ( * ' (A\) ') ' T R A C E D A X I S (a/b/c) R E A D ( * , '(A)') Q IF(Q.EQ.'A'.OR.Q.EQ.'a') IAX=I IF(Q.EQ.'B'.OR.Q.EQ.'b') IAX=2 IF(Q.EQ.'C'.OR.Q.EQ.'c') IAX=3 I F ( I A X . E Q . 0 ) G O T O 40 SKETCH
THE
a-AXIS
ALONG
T H E P A T H OF T H E
ISA=0 IF(ITP.NE.4.AND.IAX.EQ.3) THEN W R I T E ( * , '(A\)') ' S K E T C H a - A X I S R E A D ( * , '(A)') Q IF(Q.EQ.'Y'.OR.Q.EQ.'y') ISA=I E N D IF RETURN END C
?
:
c-AXIS
(Y/N)
?
:
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
417
418
J. JE.~EK
SUBROUTINE
INTPLOT
SWITCHES TO GRAPHICS, PLOTS B O U N D A R Y AND P A R A M E T E R S ON THE SCREEN
OF THE STEREOGRAM,
O U T P U T PARAMETERS: X0,Y0,DX, DY ... P O S I T I O N AND WIDTH OF THE S T E R E O G R A M (USED IN OTHER DRAWING SUBROUTINES) INPUT
PARAMETERS
(SEE SUBROUTINE
INPPAR
FOR THEIR
ON THE SCREEN
DESCRIPTION)
........................................................................
S U B R O U T I N E INTPLOT(DV, Y S , A X , Y P R , T I , X 0 , Y 0 , D X , D Y , INCLUDE 'FGRAPH.FD' INTEGER*2 D U M M Y , X M DIMENSION DV(3,3),AX(3),YS(3) DOUBLE P R E C I S I O N YPR(3) CHARACTER*30 LIST/"T'TMS RMN'H20WI6B"/ CHARACTER STR*20,STR6*6 C H A R A C T E R A X I S ( 3 ) * 2 /' A',' B',' C'/ RECORD /RCCOORD/ CURPOS RECORD /XYCOORD/ XY C C C
SET THE V I D E O M O D E
(VGA,
GRAPHICS,
640x480
PIXELS,
IAX,ITP)
16 COLORS)
DUMMY=SETVIDEOMODE($VRESI6COLOR) COLORED LOCATE
TITLE
ON THE S C R E E N
AND R E G I S T E R
FONTS
IF AVAILABLE
IN THE ACTUAL
DIRECTORY
IF(REGISTERFONTS('TMSRB.FON').LT.0) GO TO i00 IF(SETFONT(LIST).LE.0) GO TO 100 DUMMY=SETCOLOR(3) CALL MOVETO(0,0,XY) CALL O U T G T E X T ( ' E I I i p s o i d a l Particle Motion') CALL M O V E T O ( 0 , 1 8 , X Y ) CALL OUTGTEXT(' ') GO TO 200 THE CASE
TITLE
PARAMETERS:
ALL O T H E R C
WRITES
THE FONTS
WERE NOT FOUND
i00 CONTINUE D U M M Y = SETTEXTCOLOR(3) CALL S E T T E X T P O S I T I O N ( I , I , CALL O U T T E X T ( ' E l l i p s o i d a l CALL SETTEXTPOSITION(2,1, CALL OUTTEXT(' SET P O S I T I O N
AND WIDTH
CURPOS) Particle CURPOS)
Motion') ')
OF THE S T E R E O G R A M
200 C O N T I N U E X0=410. Y0=230. DX=200. DY=200. DRAW B O U N D A R Y
OF THE S T E R E O G R A M
(CIRCLE
DUMMY=SETCOLOR(2) DUMMY=ELLIPSE($GBORDER, * INT(X0-DX),INT(Y0-DY)
INT(X0+DX),INT(Y0+DY))
IF(ITP.NE.4) THEN DUMMY=SETCOLOR(15) DO XM=I,10 DUMMY=SETPIXEL(INT(X0),INT(Y0-DY)-XM) D U M M Y = S E T P I X E L ( I N T ( X 0 + D X ) + X M , INT(Y0)) END DO CALL S E T T E X T P O S I T I O N ( I , 5 2 , C U R P O S ) CALL OUTTEXT('N') CALL OUTTEXT('X') CALL S E T T E X T P O S I T I O N ( 1 5 , 8 0 , C U R P O S ) CALL OUTTEXT('E') CALL OUTTEXT('Y') END IF
Rigid triaxial ellipsoidal particle motion
DRAW THE AXES
FOR P E R S P E C T I V E
IF(ITP.EQ.4) C C C
SET THE C O L O R
419
PLOT
CALL AXPLOT(YPR, X0,Y0,DX, DY) FOR T R A C I N G
PARTICLE
DUMMY=SETCOLOR(15) C C C
WRITE
PARAMETERS
OF THE MODEL
ON THE SCREEN
DUMMY=SETTEXTCOLOR(6) CALL S E T T E X T P O S I T I O N ( 5 , I , C U R P O S ) CALL O U T T E X T ( ' V e l o c i t y Gradient Tensor') DO J = 1,3 CALL S E T T E X T P O S I T I O N ( J + 6 , I , C U R P O S ) WRITE (STR,' (3(IX,F5.3))') (DV(J,K),K=I,3) CALL OUTTEXT(STR) END DO CALL S E T T E X T P O S I T I O N ( 1 2 , 1 , C U R P O S CALL O U T T E X T ( ' A x i a l Ratio ') WRITE(STR,'(2(F4.2,3H : ),F4.2)' AX CALL S E T T E X T P O S I T I O N ( 1 4 , 3 , C U R P O S CALL OUTTEXT(STR) CALL S E T T E X T P O S I T I O N ( 1 7 , I , C U R P O S ) CALL O U T T E X T ( ' I n i t i a l Euler Angles') CALL S E T T E X T P O S I T I O N ( 1 9 , I , C U R P O S ) WRITE (STR,' (3(IX,F5.0)) ') (YS(J),J=I,3) CALL OUTTEXT(STR) C ISTEP=0 CALL S E T T E X T P O S I T I O N ( 2 2 , I , C U R P O S ) CALL O U T T E X T ( ' T i m e I n c r e m e n t :') CALL S E T T E X T P O S I T I O N ( 2 2 , 1 9 , C U R P O S ) WRITE (STR6,'(F5.3,1H)') TI CALL OUTTEXT(STR6) CALL S E T T E X T P O S I T I O N ( 2 4 , I , C U R P O S ) CALL O U T T E X T ( ' T i m e Step :') CALL S E T T E X T P O S I T I O N ( 2 4 , 1 8 , C U R P O S ) WRITE (STR6, '(I6)') ISTEP CALL OUTTEXT(STR6) C CALL S E T T E X T P O S I T I O N ( 2 6 , I , C U R P O S ) CALL O U T T E X T ( ' T r a c e d Axis : ') CALL S E T T E X T P O S I T I O N ( 2 6 , 1 5 , C U R P O S ) WRITE (STR6,'(A2)') AXIS(IAX) CALL OUTTEXT(STR6) RETURN END C
.........................................................................
C C
SUBROUTINE
AXPLOT
C C
AUXILIARY
SUBROUTINE
FOR DRAWING
AXES
IN THE P E R S P E C T I V E
C C
.........................................................................
SUBROUTINE AXPLOT(YPR,X0,Y0,DX,DY) INCLUDE 'FGRAPH.FD' INTEGER*2 DUMMY,XM, YM, XX,YY,D RECORD /XYCOORD/ XY DOUBLE P R E C I S I O N Y P R ( 3 ) , R ( 3 , 3 ) , X A ( 3 ) D I M E N S I O N XF(3) DATA D / 5 / , D D / I . 0 5 / XM=INT(X0) YM=INT(Y0) CALL RMAT(YPR,R) DO i0 I=i,3 XF(I)=R(I,I) i0 CONTINUE XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) DUMMY=SETCOLOR(10) CALL MOVETO(XM,YM, XY) D U M M Y = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL M O V E T O ( X X + D , Y Y + D , XY) D U M M Y = LINETO(XX-D, YY-D)
PLOT
420
J. JE~EK CALL M O V E T O ( X X + D , Y Y - D , X Y ) DUMMY LINETO(XX-D,YY+D) DO 20 I=i,3 XF(I)=R(I, 2) 20 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) C A L L MOVETO(XM, YM,XY) DUMMY = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL MOVETO(XX+D,YY-D,XY) DUMMY = LINETO(XX,YY) DUMMY = LINETO(XX-D,YY-D) CALL MOVETO(XX, YY, XY) D U M M Y = LINETO(XX, YY+D) DO 30 I=i,3 XF(I)=R(I,3) 30 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) CALL M O V E T O ( X M , Y M , X Y ) DUMMY = LINETO(XX,YY) XX=INT(XF(2)*DX*DD+X0) YY=INT(-XF(3)*DY*DD+Y0) CALL MOVETO(XX-D,YY-D,XY) DUMMY LINETO(XX+D,YY-D) DUMMY LINETO(XX-D,YY+D) D U M M Y = LINETO(XX+D, YY+D) =
= =
C DUMMY=SETCOLOR(II) XA(3)=0. DO 70 K = I , 3 6 0 " 4
XA(1)=COS(3.1415926/180.*FLOAT(K-I)/4.) XA(2)=SIN(3.1415926/180.*FLOAT(K-I)/4.) DO 50 I=i,3 XF(I)=0. DO 50 J=l,3 XE(I)=XF(I)+XA(J)*R(I,J) 50 C O N T I N U E XX=INT(XF(2)*DX+X0) YY=INT(-XF(3)*DY+Y0) IF(XF(1).GE.0.) DUMMY=SETPIXEL(XX,YY) IF(XF(1).LT.0.AND.(INT(FLOAT(K)/8.)*8).EQ.K) 70 C O N T I N U E RETURN END C
C C C C C C C C C C C C C
D U M M Y = S E T P I X E L ( X X , YY)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
SUBROUTINE
POPLOT
P R O J E C T S A N D PLOTS ONE POINT OF THE T R A C E D AXIS PATH TO ONE STEP OF RUNGE KUTTA)
(CORRESPONDING
INPUT PARAMETERS: Y ... V E C T O R C O N T A I N I N G E U L E R ANGLES PHI, T H E T A AND PSI D E S C R I B I N G THE I N S T A N T A N E O U S P O S I T I O N OF THE P A R T I C L E ALL O T H E R P A R A M E T E R S (FOR THEIR D E S C R I P T I O N SEE S U B R O U T I N E I N P P A R A N D INTBLOT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
S U B R O U T I N E POPLOT(Y, YPR, X0,Y0,DX,DY, ITP,IAX,ISA, ISTEP) I N C L U D E 'FGRAPH.FD' D O U B L E P R E C I S I O N Y(3),YPR(3) I N T E G E R * 2 DUMMY,XM, YM, XX,YY R E C O R D / X Y C O O R D / XY D I M E N S I O N XF(3) D A T A IDS/200/, ISK/100/ C C C
COMPUTE
COORDINATES
OF THE P L O T T E D POINT O N THE S C R E E N
CALL P R O J C T ( Y , X F , Y P R , IAX, ITP) XX=INT(XF(2)*DX+X0) YY=INT(-XF(1)*DY+Y0) GO TO (2,2,2,4), ITP PLOT POINT
IN THE S T E R E O G R A P H I C
OR L A M B E R T
PROJECTION
Rigid triaxial ellipsoidal particle motion
421
2 CONTINUE DUMMY=SETPIXEL(XX,YY) IF(ISA.EQ.0) RETURN C C
C C
SKETCH THE a-AXIS EACH ISK STEPS)
ALONG
THE
c-AXIS
PATH
(THE a - A X I S
IF(INT(FLOAT(ISTEP)/FLOAT(ISK))*ISK.EQ.ISTEP) C A L L P R O J C T ( Y , X F , Y P R , I,I) XM=INT(XF(2)*DX/10.)+XX YM=INT(-XF(1)*DY/10.)+YY C A L L M O V E T O ( X X , Y Y , XY) DUMMY = LINETO(XM,YM) CALL MOVETO(XX, YY,XY) E N D IF RETURN
IS S K E T C H E D
AFTER
THEN
C C C
THE CASE OF OF T H E P A T H
PERSPECTIVE IS P I C T U R E D
PLOT - ONLY EACH IDS-th POINT O N T H E I N V I S I B L E S I D E OF T H E S P H E R E
C 4 CONTINUE IF(XF(3).GE.0.) DUMMY=SETPIXEL(XX,YY) IF(XF(3).LT.0.AND.(INT(FLOAT(ISTEP)/FLOAT(IDS))*IDS.EQ.ISTEP)) * D U M M Y = S E T P I X E L ( X X , YY) RETURN END C
C C C C C C C C C
C C C
C C C C C
.........................................................................
SUBROUTINE
PROJCT
PROJECTS THE POINT CORRESPONDING TO THE TRACED TYPE OF PROJECTION OR IN THE PERSPECTIVE PLOT
INPUT PARAMETERS: .... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, T H E T A , PSI Y P R .. A N G L E S C O N T R O L L I N G T H E P O S I T I O N OF C O O R D I N A T E A X E S PERSPECTIVE PLOT I A X .. P A R A M E T E R T E L L I N G W H I C H P A R T I C L E A X I S IS T R A C E D I T P .. T Y P E OF P R O J E C T I O N
.........................................................................
C DIRECTIONAL
COSINES
OF T H E a - A X I S
C 1 CONTINUE XF(1)=CPH*CPS-SPH*CTH*SPS
XF(2)=SPH*CPS+CPH*CTH*SPS XF(3)=STH*SPS GO TO 5 C C
DIRECTIONAL
COSINES
OF THE
b-AXIS
C 2 CONTINUE
XF(1)=-CPH*SPS-SPH*CTH*CPS XF(2)=-SPH*SPS+CPH*CTH*CPS XF(3) = S T H * C P S G O TO 5 C C C
IN
OUTPUT PARAMETERS: XF ... C O O R D I N A T E S O F T H E P R O J E C T E D P O I N T C O R R E S P O N D I N G TO T H E T Y P E OF P R O J E C T I O N G I V E N (XF(1) A N D XF(2) A R E T H E Y - C O O R D I N A T E A N D T H E X - C O O R D I N A T E IN T H E P L O T T E D C I R C L E )
S U B R O U T I N E P R O J C T ( Y , X F , Y P R , I A X , ITP) DOUBLE PRECISION Y(3),R(3,3),YPR(3) DIMENSION XF(3),XA(3) SPH=SIN(Y(1)) CPH=COS(Y(1)) STH=SIN(Y(2)) CTH=COS(Y(2)) SPS=SIN(Y(3)) CPS=COS(Y(3)) G O TO (1,2,3), I A X C
IN T H E D E S I R E D
Y
C C
AXIS
DIRECTIONAL
COSINES
3 CONTINUE XF(1)=SPH*STH
OF T H E
c-AXIS
422
J. JE~EK XF(2)=-CPH*STH XF(3)=CTH 5 CONTINUE IF(ITP.NE.2.AND.ITP.NE.3)
GO TO
i0
L O W E R H E M I S P H E R E P R O J E C T I O N W I L L BE C O N S I D E R E D IN THE S T E R E O G R A P H I C OR L A M B E R T P R O J E C T I O N ( X F ( 3 ) . L E . 0 ) IF(XF(3).LE.0.) XF(1)=-XF(1) XF(2)=-XF(2) XF(3)=-XF(3) RECOMPUTE i0 GO TO
GO TO
COORDINATES
ORTHOGONAL
OF
I0
XF F O R THE
(11,12,13,14),
CASE
PARTIAL
TYPE
OF P R O J E C T I O N
ITP
PROJECTION
ii C O N T I N U E RETURN STEREOGRAPHIC
PROJECTION
(WULF)
12 C O N T I N U E XF(1)=XF(1)/(I.-XF(3)) XF(2)=XF(2)/(I.-XF(3)) RETURN EQUAL-AREA
PROJECTION
(LAMBERT)
13 C O N T I N U E QQ=XF(1)*XF(1)+XF(2)*XF(2) IF(QQoEQ.0.) RETURN Q=SQRT(QQ+(XF(3)+I)*(XF(3)+I)) QQ=SQRT(QQ) Q=Q/(QQ*SQRT(2.)) XF(1)=XF(1)*Q XF(2)=XF(2)*Q RETURN PERSPECTIVE
PLOT
14 C O N T I N U E DO 17 I=1,3 XA(I)=XF(I) 17 C O N T I N U E C A L L R M A T ( Y P R , R) DO 20 I=i,3 XF(I)=0. DO 20 J = l , 3 XF(I)=XF(I)+XA(J)*R(I,J) 20 C O N T I N U E XA(1)=XF(1) XF(1)=XF(3) XF (3)=XA(1) RETURN END C
......................................................................
SUBROUTINE
PREP
COMPUTES VECTOR A R O U N D ITS A X E S
O C O N T A I N I N G A N G U L A R V E L O C I T I E S OF THE P A R T I C L E IN THE C O O R D I N A T E S Y S T E M F I X E D IN THE P A R T I C L E
INPUT PARAMETERS: Y ....... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, THETA, PSI B ....... V E C T O R C H A R A C T E R I Z I N G THE S H A P E OF THE P A R T I C L E E F , O F ... S Y M M E T R I C A N D A N T I S Y M M E T R I C PART, R E S P E C T I V E L Y , OF THE V E L O C I T Y G R A D I E N T T E N S O R IN THE C O O R D I N A T E S Y S T E M F I X E D IN T H E S P A C E OUTPUT PARAMETERS: O ... V E C T O R C O N T A I N I N G A N G U L A R V E L O C I T I E S A R O U N D ITS A X E S ( O U T P U T V I A COMMON)
OF THE
PARTICLE
C ......................................................................
423
Rigid triaxial ellipsoidal particle motion SUBROUTINE PREP(Y,B,EF,OF) DOUBLE PRECISION Y(3),O(3),EB(3,3),OB(3,3),R(3,3) DIMENSION EF(3,3),OF(3,3),B(3) C O M M O N /0/ 0 COMPUTE CALL
THE
ROTATION
MATRIX
RMAT(Y,R)
C O M P U T E S Y M M E T R I C A N D A N T I S Y M M E T R I C P A R T OF THE V E L O C I T Y T E N S O R IN T H E R O T A T I N G C O O R D I N A T E S Y S T E M
GRADIENT
DO DO
20 I = i , 3 20 J = l , 3 OB(I,J)=.0 EB(I,J)=.0 DO 15 K = I , 3 DO 15 L = I , 3 RR=R(I,K)*R(J,L) EB(I,J)=EB(I,J)+RR*EF(K,L) OB(I,J)=OB(I,J)+RR*OF(K,L) 15 CONTINUE 20 C O N T I N U E EVALUATE
VECTOR
0
O(1)=B(1)*EB(3,2)-OB(3,2) O(2)=B(2)*EB(I,3)-OB(I,3) O(3)=B(3)*EB(2,1)-OB(2,1) RETURN END C
.........................................................................
C C
SUBROUTINE
RMAT
C C
COMPUTES
THE
ROTATION
MATRIX
C C C
C C C C C
INPUT PARAMETERS: ... V E C T O R C O N T A I N I N G
Y
EULER
ANGLES
PHI,
THETA,
PSI
OUTPUT PARAMETERS: R ... R O T A T I O N M A T R I X .........................................................................
SUBROUTINE RMAT(Y,R) DOUBLE PRECISION Y(3),R(3,3) S P H = S I N ( Y i) C P H = C O S ( Y i) S T H = S I N ( Y 2) C T H = C O S ( Y 2) S P S = S I N ( Y 3) C P S = C O S ( Y 3) R i,i = C P H * C P S - S P H * C T H * S P S R 1,2 = S P H * C P S + C P H * C T H * S P S R 1,3 = S T H * S P S R 2,1 = - C P H * S P S - S P H * C T H * C P S R 2,2 = - S P H * S P S + C P H * C T H * C P S R 2,3 = S T H * C P S R 3,1 = SPH*STH R 3,2 = - C P H * S T H R 3 , 3 = CTH RETURN END C
.........................................................................
SUBROUTINE COMPUTES
FCN
VALUES
OF F I R S T
DERIVATIVES
OF E U L E R
ANGLES
INPUT PARAMETERS: N ... L E N G T H OF THE V E C T O R Y Y ... V E C T O R C O N T A I N I N G E U L E R A N G L E S PHI, THETA, PSI X ... N O T U S E D (HELD F O R C O M P A T I B I L I T Y W I T H S U B R O U T I N E DVERK) O ... V E C T O R C O N T A I N I N G A N G U L A R V E L O C I T I E S OF THE P A R T I C L E A R O U N D ITS A X E S (INPUT V I A C O M M O N ) OUTPUT YP ...
PARAMETERS: FIRST DERIVATIVES
OF E U L E R
ANGLES
424
J. JE~EK
SUBROUTINE FCN(N,X,Y, YP) DOUBLE PRECISION X,Y(N),YP(N),O(3) COMMON /O/ O YP(1)=(O(1)*SIN(Y(3))+O(2)*COS(Y(3)))/SIN(Y(2)) YP(2)=O(1)*COS(Y(3))-O(2)*SIN(Y(3)) YP(3)=O(3)-YP(1)*COS(Y(2)) RETURN END