'Program to calculate crystal orientation resulting from rotation about
'    lab axes, expressed as total rotation from std position in standard
'    rotx, roty, rotz.
f$="#####.###"
5 pi=4*ATN(1)
6 delta=.001:'closeness to Ewald's sphere within which spot will be observed
7 'f$="###. ###. ###.    #.#####  #.#####  #.#####  ####.##  ####.##"
'F2$="###.######"
9 distance=100 :'mm xtal to film
10 lambda=1.08:'angstroms
11  lamdainv=1/lambda
12 a(1)=221.4:a(2)=a(1):a(3)=386.95:'angstroms
12 'a(1)=212:a(2)=a(1):a(3)=352:'angstroms
13 gamma(1)=pi/3:gamma(2)=pi/2:gamma(3)=pi/2
14 'for i=1 to 3:b(i)=1/a(i):next i
16 b(1,1)=COS(pi/6)/(a(1)*COS(pi/6)):b(2,1)=SIN(pi/6)/(a(1)*COS(pi/6)):b(3,1)=0
     ' 1 = a . a* = |a| |a*|cos 30		(a*, b* sep by 60, a orthog b*, -> a,a* sep by 30)
     ' |a|=212
     ' |a*|=1/(212 cos 30)
     'take x axis halfway between a*, b*, i.e. both 30 from x axis. Component along x is |a*|cos 30.
     'Component along y is |a*|sin(30)
     'gives 4.72, 2.72, 0  E-3 reciprocal angstroms
17 b(1,2)=b(1,1):b(2,2)=-b(2,1):b(3,2)=0
18 b(1,3)=0:b(2,3)=0:b(3,3)=1/a(3):'0,0, 2.84E-3 rec. angstrom
'19 FOR I=1 TO 3: FOR j=1 TO 3:PRINT b(I,j);:NEXT j:PRINT:NEXT I:PRINT




50 'read rot angles:
'52 READ rotx,roty,rotz
'PRINT USING f$;rotx, roty,rotz
DATA -245.49, 61.32,0
INPUT "standard rotation angles for current position";rotx, roty,rotz
56 rotx=rotx*pi/180:roty=roty*pi/180:rotz=rotz*pi/180
58 rmx(1,1)=1:rmx(2,2)=COS(rotx):rmx(3,2)=SIN(rotx):rmx(2,3)=-rmx(3,2):rmx(3,3)=rmx(2,2)
60 rmy(2,2)=1:rmy(1,1)=COS(roty):rmy(3,1)=SIN(roty):rmy(1,3)=-rmy(3,1):rmy(3,3)=rmy(1,1)
62 rmz(3,3)=1:rmz(1,1)=COS(rotz):rmz(2,1)=SIN(rotz):rmz(1,2)=-rmz(2,1):rmz(2,2)=rmz(1,1)
'multiply the matrices together, put result in prod().
 FOR i=1 TO 3:FOR K=1 TO 3:temp(i,K)=0
 FOR j=1 TO 3:temp(i,K)=temp(i,K)+rmy(i,j)*rmz(j,K)
 NEXT j:NEXT K:NEXT i
FOR i=1 TO 3:FOR K=1 TO 3:prod(i,K)=0
FOR j=1 TO 3:prod(i,K)=prod(i,K)+rmx(i,j)*temp(j,K)
NEXT j:NEXT K:NEXT i

postrot:
INPUT"further rotate about x,y, or z";a$
INPUT"how many degrees";x
x=x*pi/180
IF UCASE$(a$)="X" THEN GOSUB rotx
IF UCASE$(a$)="Y" THEN GOSUB roty
IF UCASE$(a$)="Z" THEN GOSUB rotz
mult:
FOR i=1 TO 3:FOR K=1 TO 3:temp(i,K)=0
 FOR j=1 TO 3:temp(i,K)=temp(i,K)+rmx(i,j)*prod(j,K)
 NEXT j:NEXT K:NEXT i

'932 FOR i=1 TO 3:FOR j=1 TO 3:PRINT prod(i,j),:NEXT j:PRINT:NEXT i:PRINT

 FOR i=1 TO 3: FOR j=1 TO 3: prod(i,j)=temp(i,j):NEXT j:NEXT i

'926 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmz(i,j),:NEXT j:PRINT:NEXT i:PRINT

'928 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmy(i,j),:NEXT j:PRINT:NEXT i:PRINT

'930 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmx(i,j),:NEXT j:PRINT:NEXT i:PRINT

'932 FOR i=1 TO 3:FOR j=1 TO 3:PRINT prod(i,j),:NEXT j:PRINT:NEXT i:PRINT


Factor:' R into 3 rotation matrices
 x=-prod(1,3):GOSUB arcsin:b=y:'angle of rotation about y
  IF b>pi THEN b=b-2*pi
  a=ATN(-prod(2,3)/prod(3,3))
  IF prod(3,3)/COS(b)<0 THEN a=a+pi
  IF a>pi THEN a=a-2*pi
  c=ATN(-prod(1,2)/prod(1,1))
  IF prod(1,1)/COS(b)<0 THEN c=c+pi
  IF c>pi THEN c=c-2*pi
  f2$="rotation about x,y,z:####.### ####.### ####.###"
PRINT USING f2$; 180*a/pi,180*b/pi,180*c/pi

 GOTO postrot
rotx:
rmx(1,2)=0:rmx(1,3)=0:rmx(2,1)=0:rmx(3,1)=0: rmx(1,1)=1
rmx(2,2)=COS(x):rmx(3,2)=SIN(x):rmx(2,3)=-rmx(3,2):rmx(3,3)=rmx(2,2)
RETURN
roty:
rmx(1,2)=0:rmx(3,2)=0:rmx(2,1)=0:rmx(2,3)=0: rmx(2,2)=1
rmx(1,1)=COS(x):rmx(3,1)=SIN(x):rmx(1,3)=-rmx(3,1):rmx(3,3)=rmx(1,1)
RETURN
rotz:
rmx(1,3)=0:rmx(2,3)=0:rmx(3,1)=0:rmx(3,2)=0:rmx(3,3)=1
rmx(1,1)=COS(x):rmx(2,1)=SIN(x):rmx(1,2)=-rmx(2,1):rmx(2,2)=rmx(1,1)
RETURN

arcsin:
IF ABS(x)>1.01 THEN PRINT"invalid arg for arcsin!";x:STOP
IF ABS(x)>=1 THEN y=pi/2:IF x<0 THEN y=-y:RETURN
y=ATN(x/SQR(1-x*x))
RETURN
arccos:
IF ABS(x)>1.01 THEN PRINT"invalid arg for arccos!";x:STOP
IF ABS(x)>=1 THEN y=pi/2:IF x<0 THEN y=-y:RETURN
IF x=0 THEN y=pi/2:RETURN
y=ATN(SQR(1-x*x)/x)
RETURN

notused:
'INPUT"Enter further rotation about x, y, and z (z then y then x)";rotx,roty,rotz
x=x*pi/180
'rotx=rotx*pi/180:roty=roty*pi/180:rotz=rotz*pi/180
 rmx(1,1)=1:rmx(2,2)=COS(rotx):rmx(3,2)=SIN(rotx):rmx(2,3)=-rmx(3,2):rmx(3,3)=rmx(2,2)
 rmy(2,2)=1:rmy(1,1)=COS(roty):rmy(3,1)=SIN(roty):rmy(1,3)=-rmy(3,1):rmy(3,3)=rmy(1,1)
 rmz(3,3)=1:rmz(1,1)=COS(rotz):rmz(2,1)=SIN(rotz):rmz(1,2)=-rmz(2,1):rmz(2,2)=rmz(1,1)

 FOR i=1 TO 3:FOR K=1 TO 3:temp(i,K)=0
 FOR j=1 TO 3:temp(i,K)=temp(i,K)+rmz(i,j)*prod(j,K)
 NEXT j:NEXT K:NEXT i

 FOR i=1 TO 3:FOR K=1 TO 3:prod(i,K)=0
 FOR j=1 TO 3:prod(i,K)=prod(i,K)+rmy(i,j)*temp(j,K)
 NEXT j:NEXT K:NEXT i

 FOR i=1 TO 3:FOR K=1 TO 3:temp(i,K)=0
 FOR j=1 TO 3:temp(i,K)=temp(i,K)+rmx(i,j)*prod(j,K)
 NEXT j:NEXT K:NEXT i
 
 FOR i=1 TO 3: FOR j=1 TO 3 prod(i,j)=temp(i,j):NEXT j:NEXT i


 END



926 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmz(i,j),:NEXT j:PRINT:NEXT i:PRINT

928 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmy(i,j),:NEXT j:PRINT:NEXT i:PRINT

930 FOR i=1 TO 3:FOR j=1 TO 3:PRINT rmx(i,j),:NEXT j:PRINT:NEXT i:PRINT

932 FOR i=1 TO 3:FOR j=1 TO 3:PRINT prod(i,j),:NEXT j:PRINT:NEXT i:PRINT
OPEN"o",3,"rmat 2"
PRINT#3,"rot matrix derived"
FOR i=1 TO 3:FOR j=1 TO 3:PRINT#3, USING f2$;prod(i,j),:NEXT:PRINT#3,:NEXT
PRINT#3,"RXR:"
FOR i=1 TO 3
FOR j=1 TO 3
x=0:FOR K=1 TO 3
x=x+prod(i,K)*prod(j,K)
NEXT K
PRINT#3, USING f2$; x,
NEXT j:PRINT#3,
NEXT i
CLOSE#3


946 RETURN

rotation:
950 FOR i=1 TO 3:xlab(i)=0
952 FOR j=1 TO 3:xlab(i)=xlab(i)+prod(i,j)*x(j)
954 NEXT j:NEXT i
956 FOR i=1 TO 3:x(i)=xlab(i):NEXT i
960 RETURN