Page 384 - Handbook Of Multiphase Flow Assurance

Computer code (Makogon, 1994, 1997) 383

RRT3 = 1.D0/RT3
RT32 = RT3/2.D0
RRT32 = 1.D0/RT32
DXms = rt2*DNINT(Xn(Norgn)*rrt2 -RRT3*Yn(Norgn)*rrt2)
: + rt2*DNINT(RRT32*Yn(Norgn)*rrt2)*0.5D0
DYms = rt2*DNINT(RRT32*Yn(Norgn)*rrt2)*RT32
if(DXms.ne.0.D0 .OR. DYms.ne.0.D0) then
do 498 i=1, natom
Xn(I) = Xn(I) - DXms
498 Yn(I) = Yn(I) - DYms
endif
endif
RETURN
END
C
C-------------------------------------------------------------------
C
SUBROUTINE NEWORI
IMPLICIT REAL*8(A-H,O-Z)
COMMON /types / iatomtype(1000,2),nbondfound,natom,ngroups,ncall
common /matrix/ r11,r12,r13,r21,r22,r23,r31,r32,r33,r41,r42,r43
COMMON /NATOMS/ Norgn, Nwatrs, M1, M2, M3, M4
COMMON /NWPOS / xn(1000), yn(1000), zn(1000)
COMMON /STEPS / Deltrn, Delrot, Delpiv, Acctrn, Accrot, Accpiv
COMMON /SPACE / xms(1000), yms(1000), zms(1000)
common /rotate/ A1,A2,A3,v1,v2,v3,alpha
C
C identify call type (rotation move)
ncall=2
C ... get random numbers to determine axis of rotation
delx = ROULET()
dely = ROULET()
delz = ROULET()

c choose angle of rotation; and two coordinates on rotation axis
alpha=DELROT*ROULET()
a1=xms(Norgn)
a2=yms(Norgn)
a3=zms(Norgn)
v1=a1+delx
v2=a2+dely
v3=a3+delz

c prepare the rotation matrix
c call Initrotate(A1,A2,A3,v1,v2,v3,alpha)