PROGRAM MORPHO C************************************************************************* C* Revision July 1991 * C* * C* Program MORPHO * C* * C* by Mario Nardelli * C* Dipartimento di Chimica Generale ed Inorganica, Chimica Analitica, * C* Chimica Fisica della Universita' degli Studi di Parma, * C* Centro di Studio per la Strutturistica Diffrattometrica del CNR, * C* Viale delle Scienze, 78, I-43100 PARMA (Italy) * C* * C* Literature: M. Nardelli, Proceedings of the Sixth European Crystallo- * C* graphic Meeting, Barcelona, Spain, 28 July - 1 August 1980, p.184. * C* __________________ * C* * C* This program calculates the linear axial ratio (a:b=1:c), the linear * C* interaxial angles (alpha, beta, gamma), the polar axial ratio (p0:q0: * C* r0=1), the polar interaxial angles (lambda, mu, nu), the gnomonic * C* projection constants (p0',q0', x0',y0') and coordinates (x',y') from * C* the two-circle goniometer (phi, ro) measurements and from X-ray * C* crystal data. In addition, if the distances of the faces from a point * C* inside the crystal are given, it calculates the crystal corners. * C* Stereographic and gnomonic projections are also plotted. * C* Also X-ray data are considered. * C* Input * C* * C* 1. Control card (free format): M,KM,KS,KX,KI,KV,KST,KFI * C* M= Number of cards for the title (max 5) * C* KM=0 goniometer measurements are not given * C* KM=1 goniometer measurements are given * C* KM=2 goniometer measurements and distances from the 'center' of * C* the crystal are given * C* KS=0 the e.s.d.'s of goniometer measurements are not * C* given * C* KS=1 the e.s.d.'s of goniometer measurements are given * C* KX=0 no X-ray crystal data are given * C* KX=1 X-ray crystal data are given * C* KI=0 the projection data calculated from X-ray data refer to all * C* the possible faces up to a maximum index, IM * C* KI=1 the projection data calculated from X-ray data refer only * C* to a given set of faces * C* KV=0 corners are not requested * C* KV=1 corners are requested * C* KST=0 the face indices are not given * C* KST=1 the face indices are given * C* KFI=0 the vertical-circle readings have already been transformed * C* into phi's (convention adopted: phi(010)=0) * C* KFI=1 the vertical-circle readings have not been transformed * C* into phi's * C* 2. Title cards (max.5)(not given if M=0): TITLE (80A1) * C* 3. Control card (free format): N,KDG,KV0 * C* N=number of faces (max 100) * C* KDG=0 if the angle values are introduced as degrees,minutes, * C* seconds (DDD.MMSSSS), * C* KDG=1 if the angle values are introduced as decimal degrees * C* KV0=0 if the V0 value must be deduced only from the 010 face * C* KV0=1 if the V0 value must be deduced from all the measurements * C* V0 is the vertical circle reading for the 010 face * C* Hexagonal crystals are treated as triclinic, but, if the Bravais * C* indices are given and the vertical-circle readings are such that * C* VC(100)=30 deg. and VC(010) is not given, the program rotates * C* the phi values by 30 deg., so: phi(010)=0, phi(100)=60, and * C* nu=60 deg. * C* 4. Data cards (free format)(not given if KM=0): h,k,l, V(I), RO(I), * C* SV(I),SR(I),DST(I) * C* h,k,l= MILLER indices (not given if KST=0) * C* V(I)= Vertical-circle readings (they must be such as to have * C* V(010) less than 180 deg.) * C* RO(I)= Horizontal-circle readings * C* SV(I)= Standard deviation of V(I) (not given if KS=0) * C* SR(I)= Standard deviation of RO(I) (not given if KS=0) * C* DST(I)= Distance of the face from a point inside the crystal * C* (given only if KM=2) * C* 5. X-ray data card (free format)(not given if KX=0): a,b,c,sigma(a), * C* sigma(b),sigma(c),alpha,beta,gamma,sigma(alpha),sigma(beta), * C* sigma(gamma) * C* 6. Radius card (free format)(not given if KX=0): R,N (or IM) * C* R= Radius of the stereographic circle for the table of coordinates* C* N= Number of faces (not given if KX=1 or KM=1 or KM=2) * C* IM= Maximum absolute value of the indices (not given if KI=1). * C* Care must be taken in choosing the value of IM, as the number * C* of faces greatly increases with an increase in IM, e.g. it is * C* 26 for IM=1, 98 for IM=2, 290 for IM=3, 650 for IM=4, etc. * * C* 7. Line-printer card (free format):PO,PV * C* PO= width in mm of a line-printer character * C* PV= interline distance in mm of the line-printer * C* 8. Plot card (free format)(not given if KV=0):N,RZ,RY,KDG * C* N= Number of faces (given only if KM=0) * C* RZ= Rotation angle around Z for drawing (if RZ=0 the angle is * C* set as 18.43 deg.)(not given if KV=0) * C* RY= Rotation angle around Y for drawing (if RY=0 the angle is * C* set as 9.47 deg.) (not given if KV=0) * C* KDG, as previously defined * C* 9. Data card (free format)(given only if KM=0 and KV=1): h,k,l, * C* FI(I),RO(I),DST(I) * C* * C* Subroutines and functions required: * C* (1) DEGRD * C* (2) DNSYST REQUIRES: (2 BIS) STNDV * C* (3) PRTRST * C* (4) COSTER REQUIRES: (4 BIS) PRGJT * C* (5) CORNER * C* (6) CNTR REQUIRES: (6 BIS) PLANG * C* (7) PLTCRY * C* (8) PLTSTE * C* (9) PLTGNO * C************************************************************************* CHARACTER*1 TIT(400) CHARACTER*14 FILIN,FILOUT COMMON/LBL1/C(3),A(3),SC(3),SA(3) 1 /LBL2/IN(3,100) 2 /LBL3/FI(100),RO(100) 3 /LBL4/SV(100),SR(100) 4 /LBL5/XV(3,200),XVR(3,200) 5 /LBL6/IV(3,200),IE(2,300) 6 /LBL7/V(100) 7 /LBL8/DST(100) 8 /LBL9/MSD(300),CCF(2,100) 9 /GRP1/N,NV,NE A /GRP2/X1(100),Y1(100),SX(100),SY(100) B /GRP3/AR(3),SAR(3),SNN,CSN,FI0,PSN,PCN C /GRP4/X01,SX01,Y01,SY01,P0,SP0,Q0,SQ0, 1 P01,SP01,Q01,SQ01,SRO0,CRO0 D /GRP5/IND(3,17) DATA IP,IO,RD/13,5,0.0174532925/ WRITE(*,500) 500 FORMAT(/' Key in the name of the input file'/) READ(*,'(A14)')FILIN WRITE(*,501) 501 FORMAT(/' Key in the name of the output file'/) READ(*,'(A14)')FILOUT OPEN(UNIT=13,FILE=FILIN,FORM='FORMATTED') REWIND(UNIT=13) OPEN(UNIT=5,FILE=FILOUT,FORM='FORMATTED') REWIND(UNIT=5) READ(IP,*)M,KM,KS,KX,KI,KV,KST,KFI IF(M.EQ.0)GO TO 1 L=M*80 READ(IP,101)(TIT(J),J=1,L) 101 FORMAT(80A1) WRITE(IO,102)(TIT(J),J=1,L) 102 FORMAT(1X,80A1,/(1X,80A1)) 1 DO 2 I=1,100 SV(I)=5E-4 2 SR(I)=5E-4 IF(KM.EQ.0)GO TO 16 WRITE(IO,119) 119 FORMAT(/9X,'Two-circle optical goniometer measurements') IF(KST.EQ.0)GO TO 44 IF(KM.NE.2)WRITE(IO,111) 111 FORMAT(/16X,'Faces',6X,'Goniometer readings') IF(KM.EQ.2)WRITE(IO,112) 112 FORMAT(/12X,'Faces',6X,'Goniometer readings',6X,'Distances') 44 CONTINUE READ(IP,*)N,KDG,KV0 KDG=KDG+1 KV0=KV0+1 IF(KST.EQ.0)GO TO 45 IF(KM.EQ.1)GO TO 3 IF(KM.EQ.2.AND.KV.EQ.1)GO TO 7 3 IF(KS.EQ.0)GO TO 5 DO 4 I=1,N 4 READ(IP,*)(IN(J,I),J=1,3),V(I),RO(I),SV(I),SR(I) GO TO 37 5 DO 6 I=1,N 6 READ(IP,*)(IN(J,I),J=1,3),V(I),RO(I) GO TO 11 7 IF(KS.EQ.0)GO TO 9 DO 8 I=1,N 8 READ(IP,*)(IN(J,I),J=1,3),V(I),RO(I),SV(I),SR(I),DST(I) GO TO 39 9 DO 10 I=1,N 10 READ(IP,*)(IN(J,I),J=1,3),V(I),RO(I),DST(I) GO TO 35 11 WRITE(IO,113) 113 FORMAT(10X,'N',4X,'h k l',7X,'Vc',12X,'Ro') DO 34 I=1,N 34 WRITE(IO,114)I,(IN(J,I),J=1,3),V(I),RO(I) 114 FORMAT(8X,I3,2X,3I3,3X,F7.2,7X,F7.2) GO TO 41 35 WRITE(IO,115) 115 FORMAT(6X,'N',4X,'h k l',7X,'Vc',12X,'Ro',10X,'d') DO 36 I=1,N 36 WRITE(IO,116)I,(IN(J,I),J=1,3),V(I),RO(I),DST(I) 116 FORMAT(4X,I3,2X,3I3,3X,F7.2,7X,F7.2,7X,F5.2) GO TO 41 37 WRITE(IO,113) DO 38 I=1,N 38 WRITE(IO,117)I,(IN(J,I),J=1,3),V(I),SV(I),RO(I),SR(I) 117 FORMAT(8X,I3,2X,3I3,3X,F7.2,'(',F4.2,')',1X,F7.2,'(', 1 F4.2,')') GO TO 41 39 WRITE(IO,115) DO 40 I=1,N 40 WRITE(IO,118)I,(IN(J,I),J=1,3),V(I),SV(I),RO(I),SR(I),DST(I) 118 FORMAT(4X,I3,2X,3I3,3X,F7.2,'(',F4.2,')',1X,F7.2,'(', 1 F4.2,')',F6.2) GO TO 41 45 WRITE(IO,120) 120 FORMAT(/17X,'N',12X,'Vc',8X,'Ro') DO 46 I=1,N READ(IP,*)V(I),RO(I) 46 WRITE(IO,121)I,V(I),RO(I) 121 FORMAT(15X,I3,5X,2(3X,F7.3)) 41 CONTINUE GO TO(12,14),KDG 12 DO 13 I=1,N B=V(I) V(I)=DEGRD(B) B=RO(I) RO(I)=DEGRD(B) IF(KS.EQ.0)GO TO 13 B=SV(I) SV(I)=DEGRD(B) B=SR(I) SR(I)=DEGRD(B) 13 CONTINUE 14 DO 15 I=1,N V(I)=V(I)*RD RO(I)=RO(I)*RD SV(I)=SV(I)*RD SR(I)=SR(I)*RD IF(KS.EQ.0)SV(I)=1E-2 IF(KS.EQ.0)SR(I)=1E-2 15 CONTINUE IF(KST.EQ.0)GO TO 57 C-----Calculation of morphological parameters from two-circle goniometer C-----data CALL DNSYST(KV0,N,KFI,KPR) IF(KPR.EQ.1)GO TO 16 C-----Calculation of gnomonic and stereographic coordinates and C-----printing of morphological parameters CALL PRTRST(RD,KS,N) GO TO 16 57 DO 58 I=1,N 58 FI(I)=V(I) 16 IF(KX.EQ.0)GO TO 21 C-----Calculation of morphological parameters and gnomonic and C-----stereographic coordinates from X-ray crystal data READ(IP,*)C,SC,A,SA IF(KI.EQ.0)GO TO 19 IF(KM.EQ.0)GO TO 17 READ(IP,*)R GO TO 18 17 READ(IP,*)R,N READ(IP,*)((IN(J,I),J=1,3),I=1,N) 18 KDC=1 GO TO 20 19 READ(IP,*)R,IM KDC=0 IF(IM.EQ.0)IM=1 20 CALL COSTER(R,KDC,IM,N,X01,Y01,PSN,PCN,Q01) 21 READ(IP,*)PO,PV IF(KV.EQ.0)GO TO 33 IF(KM.EQ.2)GO TO 42 WRITE(IO,103) 103 FORMAT(/14X,'Faces',6X,'Goniometer readings',5X, 1 'Distances'/8X,'N',4X,'h k l',7X,'Phi',11X,'Ro',9X,'d') 42 KM=KM+1 GO TO(22,26,27),KM C-----Calculation of corner coordinates and edges 22 READ(IP,*)N,RZ,RY,KDG KDG=KDG+1 DO 25 I=1,N READ(IP,*)(IN(J,I),J=1,3),FI(I),RO(I),DST(I) GO TO(23,24),KDG 23 B=FI(I) FI(I)=DEGRD(B) B=RO(I) RO(I)=DEGRD(B) 24 FI(I)=FI(I)*RD 25 RO(I)=RO(I)*RD GO TO 28 26 READ(IP,*)RZ,RY READ(IP,*)(DST(I),I=1,N) GO TO 28 27 READ(IP,*)RZ,RY 28 IF(N.LT.4)GO TO 32 IF(KM.EQ.3)GO TO 43 DO 29 I=1,N B=FI(I)/RD D=RO(I)/RD 29 WRITE(IO,104)I,(IN(J,I),J=1,3),B,D,DST(I) 104 FORMAT(6X,I3,2X,3I3,3X,F7.2,7X,F7.2,6X,F5.2) 43 CONTINUE CALL CORNER(RZ,RY) CALL CNTR WRITE(IO,108) 108 FORMAT(/3X,' Coordinates of the centres of the faces in the ', 1 'drawing'//12X,'N',5X,'h k l',8X,'yr',10X,'zr') DO 31 I=1,N 31 WRITE(IO,109)I,(IN(J,I),J=1,3),(CCF(K,I),K=1,2) 109 FORMAT(10X,I3,3X,3I3,2(5X,F7.3)) CALL PLTCRY(PO,PV,RZ,RY) GO TO 33 32 WRITE(IO,110) 110 FORMAT(//' The number of the faces is not enough' / 1 1X,' to define a polyhedron') 33 IF(KM.EQ.0.AND.KX.EQ.0)GO TO 49 KFL=0 IF(KM.EQ.0)KFL=1 54 WRITE(IO,122) 122 FORMAT(///' Stereographic projection (radius = 10 cm)') IF(KFL.EQ.0)GO TO 50 WRITE(IO,124) 124 FORMAT(' from X-ray data'/) GO TO 51 50 WRITE(IO,123) 123 FORMAT(' from optical goniometer data'/) 51 IF(KFL.EQ.0)GO TO 55 IF(KX.EQ.0)GO TO 49 KST=1 N=17 DO 48 I=1,N FI(I)=SV(I) RO(I)=SR(I) DO 48 J=1,3 48 IN(J,I)=IND(J,I) 55 CALL PLTSTE(N,PO,PV,KST) IF(KFL.EQ.0.AND.KPR.EQ.1)GO TO 56 WRITE(IO,125) 125 FORMAT(///' Gnomonic projection') IF(KFL.EQ.0)GO TO 52 WRITE(IO,124) GO TO 53 52 WRITE(IO,123) 53 CALL PLTGNO(N,PO,PV,KST,X01,Y01,PSN,PCN,Q01) 56 IF(KFL.EQ.1)GO TO 49 KFL=KFL+1 IF(KX.EQ.0)GO TO 49 GO TO 54 49 CONTINUE STOP END FUNCTION DEGRD(A) IA=A D1=IA DM=(A-D1)*100. IDM=DM D2=IDM DS=DM-D2 DEGRD=D1+D2/60.+DS/36. RETURN END SUBROUTINE STNDV(SX2,SW,XM,N,SGX,SX) C-----This routine calculates the weighted mean of a set of data C-----with its e.s.d. and puts some parameters equal to zero XM=SX/SW A=0. IF(N.GT.1)A=(SX2/SW-XM**2)/(N-1) SGX=SQRT(A+1E-4/SW) N=0 SX=0 SX2=0 SW=0 RETURN END SUBROUTINE PLANG(Y1,Z1,Y2,Z2,Y3,Z3,AG) C-----This routine calculates an angle in a plane D1=Y1-Y2 D2=Z1-Z2 D3=Y1-Y3 D4=Z1-Z3 S1=SQRT(D1**2+D2**2) S2=SQRT(D3**2+D4**2) D1=D1/S1 D2=D2/S1 D3=D3/S2 D4=D4/S2 CSA=D1*D3+D2*D4 IF(ABS(CSA).LT.1E-6)AA=1. IF(ABS(CSA).GT.1.)AA=0. IF(ABS(CSA).GT.1E-6.OR.ABS(CSA).LT.1.)AA=1-CSA**2 IF(AA.LT.0.)AA=0. AG=ATAN2(SQRT(AA),CSA) RETURN END SUBROUTINE PRGJT(NH,X01,Y01,P01,Q01,RO0,FI0,SNNU,CSNU,R, 1 FN,RN) C-----This routine calculates the gnomonic and stereographic C-----coordinates of a face DIMENSION NH(3),AIN(3) RD=0.0174532925 IT=5 C-----Calculation of gnomonic x' & y' coordinates IF(NH(1).EQ.0.AND.NH(2).EQ.0)GO TO 5 IF(NH(3).LT.0)GO TO 15 DO 1 I=1,3 1 AIN(I)=NH(I) GO TO 16 15 DO 17 I=1,3 17 AIN(I)=-NH(I) 16 CONTINUE IF(NH(3).EQ.0)GO TO 6 D1=AIN(1)/AIN(3) D2=AIN(2)/AIN(3) IF(ABS(D1).LT.1E-7)GO TO 10 X1=X01+D1*P01*SNNU GO TO 11 10 X1=X01 11 IF(ABS(D1).LT.1E-7.AND.ABS(D2).LT.1E-7)GO TO 12 Y1=Y01+D1*P01*CSNU+D2*Q01 GO TO 13 12 Y1=Y01 13 CONTINUE C-----Calculation of stereographic Phi & Ro coordinates IF(ABS(Y1).LT.1E-7)GO TO 2 FI=ATAN2(X1,Y1) IF(FI.LT.0.)FI=360*RD+FI GO TO 3 2 IF(X1.GE.0.)FI=90.*RD IF(X1.LT.0.)FI=270.*RD 3 SFI=SIN(FI) IF(ABS(SFI).LT.1E-7)GO TO 4 RO=ATAN2(X1,SFI) IF(RO.LT.0.)RO=RO+180*RD IF(X1.EQ.0..AND.Y1.LT.0.)RO=RO-180.*RD GO TO 8 4 RO=ATAN(Y1) GO TO 8 14 FI=ATAN2(X01,Y01) RO=ATAN2(X01,SIN(FI)) GO TO 8 6 IF(NH(1).EQ.0)GO TO 7 GO TO 9 7 IF(NH(2).GT.0)FI=0. IF(NH(2).LT.0)FI=180.*RD RO=90.*RD GO TO 8 9 DEN=(AIN(2)/AIN(1))*(Q01/P01)+CSNU FI=ATAN2(SNNU,DEN) IF(NH(1).LT.0)FI=FI+180.*RD RO=90.*RD GO TO 19 5 FI=FI0 X1=X01 Y1=Y01 RO=RO0 8 CONTINUE IF(NH(3).LT.0)GO TO 18 GO TO 19 18 FI=FI+180.*RD RO=180.*RD-RO IF(RO.GT.(180.*RD))RO=360.*RD-RO 19 CONTINUE IF(RO.LT.0.)GO TO 22 GO TO 23 22 RO=-RO IF(ABS(FI-180.*RD).GT.1E-4)FI=FI+180.*RD 23 CONTINUE C-----Calculation of stereographic XS & YS coordinates RS=RO/2. IF(RO.GT.(90.*RD))RS=(180.*RD-RO)/2. R1=R*SIN(RS)/COS(RS) XS=R1*SIN(FI) YS=R1*COS(FI) IF(FI.GE.(360.*RD))FI=FI-360.*RD FN=FI RN=RO FI=FI/RD RO=RO/RD IF(NH(3).LE.0)GO TO 20 WRITE(IT,100)NH,X1,Y1,FI,RO,XS,YS 100 FORMAT(1X,3I3,3X,2(1X,F7.3),2X,2(1X,F7.2),2X,2(F7.3)) GO TO 21 20 WRITE(IT,101)NH,FI,RO,XS,YS 101 FORMAT(1X,3I3,8X,'-',7X,'-',5X,2(F7.2,1X),1X,2F7.3) 21 CONTINUE RETURN END SUBROUTINE DNSYST(KV0,N,KFI,KPR) C-----This routine calculates the morphological parameters from C-----optical goniometer measurements IMPLICIT INTEGER(H) COMMON/LBL1/C(3),A(3),SC(3),SA(3) 1 /LBL2/IN(3,100) 2 /LBL3/FI(100),RO(100) 3 /LBL4/SV(100),SR(100) 4 /LBL7/V(100) 3 /GRP2/X1(100),Y1(100),SX(100),SY(100) 3 /GRP3/AR(3),SAR(3),SNN,CSN,FI0,PSN,PCN 7 /GRP4/X01,SX01,Y01,SY01,P0,SP0,Q0,SQ0, 8 P01,SP01,Q01,SQ01,SRO0,CRO0 DATA IT,P2/5,1.57079632679/ IF(N.LT.3)GO TO 52 NM1=N-1 NM2=N-2 V2=0. SV2=0. SW1=0. NT=0 KVT=0 IF(KFI.EQ.0)GO TO 6 C-----Calculation of V0 DO 1 I=1,N IF(IN(1,I).EQ.0.AND.IN(2,I).EQ.1.AND.IN(3,I).EQ.0)GO TO 2 1 CONTINUE KV0=2 GO TO 3 2 V0=V(I) VPR=V(I) SV0=SV(I) KVT=1 W=1./(SV(I)*1E2)**2 IF(ABS(V(I)).LE.1E-6)SV2=0. IF(ABS(V(I)).GT.1E-6)SV2=V(I)**2*W V2=V(I)*W SW1=W NT=NT+1 3 IF(KV0.EQ.1)GO TO 6 DO 56 I=1,NM1 IP1=I+1 IF(IN(1,I).EQ.0.OR.IN(3,I).LE.0.)GO TO 56 DO 4 J=IP1,N IF(IN(1,J).EQ.0.OR.IN(3,J).LE.0)GO TO 4 IF((FLOAT(IN(1,I))/FLOAT(IN(3,I))).NE.(FLOAT(IN(1,J))/ 1 FLOAT(IN(3,J))))GO TO 4 SN1=SIN(V(I)) SN2=SIN(V(J)) CS1=COS(V(I)) CS2=COS(V(J)) SNR1=SIN(RO(I)) SNR2=SIN(RO(J)) CSR1=COS(RO(I)) CSR2=COS(RO(J)) TG1=SNR1/CSR1 TG2=SNR2/CSR2 AN=SN2*TG2-SN1*TG1 DN=CS2*TG2-CS1*TG1 IF(ABS(DN).LT.1E-7)GO TO 4 V3=ATAN2(AN,DN) IF(V3.GE.4.*P2)V3=V3-4.*P2 IF(V3.GE.2.*P2.AND.AN.LT.0..AND.DN.LT.0.)V3=V3-2.*P2 IF(V3.GT.3.*P2)V3=V3-4.*P2 IF(ABS(V3-2.*P2).LT.1E-5)V3=0. IF(KVT.EQ.1.AND.ABS(V0-V3).GT..01)GO TO 4 CSV0=COS(V3) TGV0=SIN(V3)/CSV0 W=((CSV0**2)/DN)**2 W=W*((TG1*(CS1+SN1*TGV0)*SV(I))**2+ 1 (TG2*(CS2+SN2*TGV0)*SV(J))**2+ 2 ((SN1-CS1*TGV0)*SR(I)/CSR1**2)**2+ 3 ((SN2-CS2*TGV0)*SR(J)/CSR2**2)**2) W=1./(W*1E4) NT=NT+1 SW1=SW1+W IF(ABS(V3).LT.1E-7)GO TO 4 V2=V2+V3*W SV2=SV2+V3**2*W 4 CONTINUE 56 CONTINUE DO 58 I=1,NM2 IF(IN(3,I).NE.0)GO TO 58 IP1=I+1 DO 57 J=IP1,NM1 IF(IN(3,J).NE.0)GO TO 57 JP1=J+1 DO 5 K=JP1,N IF(IN(3,K).NE.0)GO TO 5 IA=(IN(1,J)*IN(2,K)-IN(1,K)*IN(2,J))*IN(1,I) IDM=(IN(1,I)*IN(2,K)-IN(1,K)*IN(2,I))*IN(1,J) IF(IA.EQ.0.OR.IDM.EQ.0)GO TO 5 Q=FLOAT(IA)/FLOAT(IDM) SN1=SIN(V(J)-V(I)) SN2=SIN(V(K)-V(J)) CS1=COS(V(J)-V(I)) CS2=COS(V(K)-V(J)) IF(ABS(SN1).LT.1E-7.OR.ABS(SN2).LT.1E-7)GO TO 5 TG1=CS1/SN1 TG2=CS2/SN2 TGV0=Q*TG2-(1.-Q)*TG1 IF(ABS(TGV0).LT.1E-7)GO TO 5 AA=ATAN(1./TGV0) B=SIN(AA)**4 W=B*(((Q/SN2**2+(1.-Q)/SN1**2)*SV(J))**2+(Q/SN2**2*SV(K)) 1 **2+((1.-Q)/SN1**2*SV(I))**2)-SV(J)**2 W=ABS(1./(W*1E4)) V3=AA+V(J) IF(V3.GE.4.*P2)V3=V3-4.*P2 IF(V3.GE.2.*P2.AND.V3.LT.3.*P2)V3=V3-2.*P2 IF(V3.GT.3.*P2)V3=V3-4.*P2 IF(ABS(V3-2.*P2).LT.1E-5)V3=0. IF(KVT.EQ.1.AND.ABS(V0-V3).GT..01)GO TO 5 NT=NT+1 SW1=SW1+W IF(ABS(V3).LT.1E-7)GO TO 5 V2=V2+V3*W SV2=SV2+V3**2*W 5 CONTINUE 57 CONTINUE 58 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,V0,NT,SV0,V2) 55 CONTINUE IF(KVT.EQ.1.AND.ABS(V0-VPR).GT..001)V0=VPR C-----Calculation of Phi's, and x' & y' 6 DO 8 I=1,N FI(I)=V(I) IF(KFI.EQ.1)FI(I)=V(I)-V0 SV(I)=SQRT(SV(I)**2+SV0**2) IF(ABS(FI(I)).LT.1E-3)FI(I)=0. IF(FI(I).LT.0.)FI(I)=4.*P2+FI(I) X1(I)=0. Y1(I)=0. SX(I)=0. SY(I)=0. SN1=SIN(FI(I)) CS1=COS(FI(I)) SN2=SIN(RO(I)) CS2=COS(RO(I)) IF(ABS(CS2).LT.1E-7)GO TO 8 TG1=SN2/CS2 X1(I)=SN1*TG1 Y1(I)=CS1*TG1 SX(I)=SQRT((Y1(I)*SV(I))**2+(SN1/CS2**2*SR(I))**2) SY(I)=SQRT((X1(I)*SV(I))**2+(CS1/CS2**2*SR(I))**2) 8 CONTINUE NT=0 V2=0. SV2=0. SW1=0. C-----Calculation of p0'*sin(nu) DO 59 I=1,NM1 IF(IN(3,I).LE.0)GO TO 59 IP1=I+1 DO 9 J=IP1,N IF(IN(3,J).LE.0)GO TO 9 LL=IN(3,I)*IN(3,J) HL=IN(1,I)*IN(3,J)-IN(1,J)*IN(3,I) IF(HL.EQ.0)GO TO 9 AA=X1(I)-X1(J) IF(ABS(AA).LT.1E-7)GO TO 9 PSN=LL*AA/HL NT=NT+1 LL=LL*LL HL=HL*HL AA=FLOAT(LL)/FLOAT(HL) B=((SX(I)*1E2)**2+(SX(J)*1E2)**2)*AA**2 IF(B.LT.1E-7)B=1E-6 W=1./B V2=V2+PSN*W SV2=SV2+(PSN**2)*W SW1=SW1+W 9 CONTINUE 59 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,PSN,NT,SPSN,V2) C-----Calculation of p0'*cos(nu) DO 61 I=1,NM2 IF(IN(3,I).LE.0)GO TO 61 IP1=I+1 DO 60 J=IP1,NM1 IF(IN(3,J).LE.0)GO TO 60 JP1=J+1 DO 10 K=JP1,N IF(IN(3,K).EQ.0)GO TO 10 LL1=IN(3,I)*IN(3,J) LL2=IN(3,J)*IN(3,K) KL1=IN(2,J)*IN(3,K) KL2=IN(2,K)*IN(3,J) KL3=IN(2,I)*IN(3,J) KL4=IN(2,J)*IN(3,I) HL1=IN(1,I)*IN(3,J) HL2=IN(1,J)*IN(3,I) HL3=IN(1,J)*IN(3,K) HL4=IN(1,K)*IN(3,J) AN=LL1*(Y1(I)-Y1(J))*(KL1-KL2)-LL2*(Y1(J)-Y1(K))*(KL3-KL4) HL=(HL1-HL2)*(KL1-KL2)-(HL3-HL4)*(KL3-KL4) IF(HL.EQ.0.)GO TO 10 V4=AN/HL L1=LL1*(KL1-KL2) L2=LL2*(KL3-KL4) AA=SY(I)*1E2 B=SY(J)*1E2 CC=SY(K)*1E2 W=(L1**2*(AA**2+B**2)+L2**2*(B**2+CC**2))/HL**2 W=1./W NT=NT+1 V2=V2+V4*W SV2=SV2+(V4**2)*W SW1=SW1+W 10 CONTINUE 60 CONTINUE 61 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,PCN,NT,SPCN,V2) C-----Calculation of q0' DO 62 I=1,NM1 IF(IN(3,I).LE.0)GO TO 62 IP1=I+1 DO 11 J=IP1,N IF(IN(3,J).LE.0)GO TO 11 HL1=IN(1,J)*IN(3,I) HL2=IN(1,I)*IN(3,J) LL1=IN(3,I)*IN(3,J) KL1=IN(2,I)*IN(3,J) KL2=IN(2,J)*IN(3,I) AN=PCN*(HL1-HL2)+LL1*(Y1(I)-Y1(J)) HL=KL1-KL2 IF(HL.EQ.0)GO TO 11 V4=ABS(AN/HL) IF(SPCN.LE.1E-6)AA=0. IF(SPCN.GT.1E-6)AA=((HL1-HL2)*SPCN*1E2)**2 B=LL1**2*((SY(I)*1E2)**2+(SY(J)*1E2)**2) IF((AA+B).LE.1E-7)W=0. IF((AA+B).GT.1E-7)W=HL**2/(AA+B) NT=NT+1 V2=V2+V4*W SV2=SV2+V4**2*W SW1=SW1+W 11 CONTINUE 62 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,Q01,NT,SQ01,V2) C-----Calculation of x0' DO 12 I=1,N IF(IN(3,I).LE.0)GO TO 12 X01=X1(I)-IN(1,I)*PSN/IN(3,I) AA=FLOAT(IN(1,I))/FLOAT(IN(3,I)) IF(SX(I).LE.1E-7)B=0. IF(SX(I).GT.1E-7)B=(SX(I)*1E2)**2 IF(IN(1,I).EQ.0.OR.SPSN.LE.1E-7)CC=0. IF(IN(1,I).NE.0.AND.SPSN.GT.1E-7)CC=(AA*SPSN*1E2)**2 IF((B+CC).LE.1E-7)W=0. IF((B+CC).GT.1E-7)W=1./(B+CC) NT=NT+1 V2=V2+X01*W SV2=SV2+X01**2*W SW1=SW1+W 12 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,X01,NT,SX01,V2) C-----Calculation of y0' DO 13 I=1,N IF(IN(3,I).LE.0)GO TO 13 Y01=Y1(I)-(IN(2,I)*Q01)/IN(3,I)-(IN(1,I)*PCN)/IN(3,I) AA=FLOAT(IN(2,I))/FLOAT(IN(3,I)) B=FLOAT(IN(1,I))/FLOAT(IN(3,I)) IF(SY(I).LE.1E-7)CC=0. IF(SY(I).GT.1E-7)CC=(SY(I)*1E2)**2 IF(IN(2,I).EQ.0.OR.SQ01.LE.1E-7)D=0. IF(IN(2,I).NE.0.AND.SQ01.GT.1E-7)D=(AA*SQ01*1E2)**2 IF(IN(1,I).EQ.0.OR.SPCN.LE.1E-7)E=0. IF(IN(1,I).NE.0.AND.SPCN.GT.1E-7)E=(B*SPCN*1E2)**2 IF((CC+D+E).LE.1E-7)W=0. IF((CC+D+E).GT.1E-7)W=1./(CC+D+E) NT=NT+1 V2=V2+Y01*W SV2=SV2+Y01**2*W SW1=SW1+W 13 CONTINUE IF(NT.EQ.0)GO TO 52 CALL STNDV(SV2,SW1,Y01,NT,SY01,V2) C-----Calculation of nu IF(ABS(PCN).LT.1E-7)GO TO 14 AR(3)=ATAN2(PSN,PCN) IF(ABS(AR(3)-P2).LT.2E-4)AR(3)=P2 SNN=SIN(AR(3)) CSN=COS(AR(3)) IF(ABS(CSN).LT.1E-5)GO TO 14 SAR(3)=ABS(SNN*CSN*SQRT((SPSN/PSN)**2+(SPCN/PCN)**2)) GO TO 15 14 AR(3)=P2 SAR(3)=1E-6 CSN=0. SNN=1. 15 CONTINUE C-----Calculation of p0' P01=ABS(PSN/SNN) SP01=ABS(P01*SQRT((SPSN/PSN)**2+(CSN/SNN*SAR(3))**2)) C-----Calculation of Phi0 & Ro0 IF(ABS(X01).LT.1E-7.AND.ABS(Y01).LT.1E-7)GO TO 18 IF(ABS(X01).LT.1E-7)GO TO 19 IF(ABS(Y01).LT.1E-7)GO TO 20 GO TO 21 18 FI0=0. RO0=0. SGR0=1E-6 SGF0=1E-6 GO TO 22 19 FI0=0. RO0=ATAN(Y01) SGF0=1E-6 CRO0=COS(RO0) SGR0=CRO0**2*SY01 GO TO 22 20 FI0=P2 RO0=ATAN(X01) SGF0=1E-6 CRO0=COS(RO0) SGR0=CRO0**2*SX01 GO TO 22 21 FI0=ATAN2(X01,Y01) SFI0=SIN(FI0) CFI0=COS(FI0) IF(ABS(X01).GE.ABS(Y01))RO0=ATAN2(X01,SFI0) IF(ABS(Y01).GT.ABS(X01))RO0=ATAN2(Y01,CFI0) SRO0=SIN(RO0) CRO0=COS(RO0) SGF0=CFI0**2/Y01*SQRT(SX01**2+(X01/Y01*SY01)**2) IF(ABS(X01).GE.ABS(Y01))SGR0=CRO0**2/SFI0*SQRT(SX01**2+ 1 (X01*CFI0/SFI0*SGF0)**2) IF(ABS(Y01).GT.ABS(X01))SGR0=CRO0**2/CFI0*SQRT(SY01**2+ 1 (Y01*SFI0/CFI0*SGF0)**2) GO TO 33 22 CONTINUE SFI0=SIN(FI0) CFI0=COS(FI0) SRO0=SIN(RO0) CRO0=COS(RO0) 33 CONTINUE C-----Calculation of x0 & y0 X0=SRO0*SFI0 Y0=SRO0*CFI0 SX0=SQRT((CRO0*SFI0*SGR0)**2+(SRO0*CFI0*SGF0)**2) SY0=SQRT((CRO0*CFI0*SGR0)**2+(SRO0*SFI0*SGF0)**2) C-----Calculation of p0 & q0 P0=ABS(P01*CRO0) Q0=ABS(Q01*CRO0) SP0=SQRT((CRO0*SP01)**2+(P01*SRO0*SGR0)**2) SQ0=SQRT((CRO0*SQ01)**2*(Q01*SRO0*SGR0)**2) C-----Calculation of lambda IF(ABS(Y0).LT.1E-6)GO TO 23 SNL=SQRT(1.-Y0**2) CSL=Y0 AR(1)=ATAN2(SNL,CSL) SAR(1)=SY0/SNL GO TO 24 23 AR(1)=P2 SAR(1)=1E-6 SNL=1. CSL=0. 24 CONTINUE C-----Calculation of mu CSM=X0*SNN IF(ABS(Y0).GT.1E-7.AND.ABS(CSM).GT.1E-7)CSM=Y0*CSN+X0*SNN IF(ABS(CSM).LE.1E-7)GO TO 25 SNM=SQRT(1.-CSM**2) AR(2)=ATAN2(SNM,CSM) SAR(2)=SQRT((CSN*SY0)**2+(SNN*SX0)**2+((X0*CSN-Y0*SNN) 1 *SAR(3))**2)/SNM GO TO 26 25 AR(2)=P2 SAR(2)=1E-6 CSM=0. SNM=1. 26 CONTINUE C-----Calculation of the axial ratios a/b & c/b C(1)=ABS(Q01*SNL/(P01*SNM)) C(3)=ABS(Q01*CRO0*SNN/SNM) B=(SQ01/Q01)**2 D=(CSM/SNM*SAR(2))**2 SC(1)=C(1)*SQRT(B+(CSL/SNL*SAR(1))**2+(SP01/P01)**2+D) SC(3)=C(3)*SQRT(B+(SRO0/CRO0*SGR0)**2+(CSN/SNN*SAR(3))**2 1 +D) C-----Calculation of alpha, beta & gamma SG=(AR(1)+AR(2)+AR(3))/2. SNSG=SIN(SG) SS1=SIN(SG-AR(1)) SS2=SIN(SG-AR(2)) SS3=SIN(SG-AR(3)) AA=SG+SG SS4=SIN(AA-AR(1)) SS5=SIN(AA-AR(2)) SS6=SIN(AA-AR(3)) SNA=SQRT((SNSG*SS1)/(SNM*SNN)) SNB=SQRT((SNSG*SS2)/(SNN*SNL)) SNG=SQRT((SNSG*SS3)/(SNL*SNM)) IF(ABS(SNA-1.).LT.1E-6)GO TO 27 A(1)=2.*ATAN2(SNA,SQRT(1.-SNA**2)) GO TO 28 27 A(1)=P2 28 CONTINUE IF(ABS(SNB-1.).LT.1E-6)GO TO 29 A(2)=2.*ATAN2(SNB,SQRT(1.-SNB**2)) GO TO 30 29 A(2)=P2 30 CONTINUE IF(ABS(SNG-1.).LT.1E-6)GO TO 31 A(3)=2.*ATAN2(SNG,SQRT(1.-SNG**2)) GO TO 32 31 A(3)=P2 32 CONTINUE SA(1)=ABS(SQRT((SNL/2.*SAR(1))**2+((SS4/2.+SNSG*SS1*CSM/ 1 SNM)*SAR(2))**2+((SS4/2.+SNSG*SS1*CSN/SNN)*SAR(3))**2)/ 2 (SNA*SNM*SNN)) SA(2)=ABS(SQRT((SNM/2.*SAR(2))**2+((SS5/2.+SNSG*SS2*CSL/ 1 SNL)*SAR(1))**2+((SS5/2.+SNSG*SS2*CSN/SNN)*SAR(3))**2)/ 2 (SNB*SNL*SNN)) SA(3)=ABS(SQRT((SNN/2.*SAR(3))**2+((SS6/2.+SNSG*SS3*CSL/ 1 SNL)*SAR(1))**2+((SS6/2.+SNSG*SS3*CSM/SNM)*SAR(2))**2)/ 2 (SNG*SNL*SNM)) KPR=0 GO TO 53 52 WRITE(IT,113) 113 FORMAT(//' The goniometer data are not enough to', 1 ' calculate the morphological'/' parameters') DO 64 I=1,N X1(I)=0. Y1(I)=0. IF(ABS(RO(I)-P2).LT.1E-4)GO TO 64 TG1=SIN(RO(I))/COS(RO(I)) X1(I)=SIN(FI(I))*TG1 Y1(I)=COS(FI(I))*TG1 KPR=1 64 CONTINUE 53 RETURN END SUBROUTINE PRTRST(RD,KS,N) C-----This routine calculates Phi(calc.)'s & Ro(calc.)'s C-----and the angles formed by the faces with (100),(010),(001) C-----and prints these results and those given by routine DNSYST DIMENSION FC(100),RC(100) COMMON/LBL1/C(3),A(3),SC(3),SA(3) 1 /LBL2/IN(3,100) 2 /LBL3/FI(100),RO(100) 3 /GRP2/X1(100),Y1(100),SX(100),SY(100) 4 /GRP3/AR(3),SAR(3),SNN,CSN,FI0,PSN,PCN 5 /GRP4/X01,SX01,Y01,SY01,P0,SP0,Q0,SQ0, 1 P01,SP01,Q01,SQ01,SRO0,CRO0 DATA IT,P2/5,1.57079632679/ C-----Calculation of Phi(calc.) & Ro(calc.) DO 43 I=1,N IF(IN(3,I).EQ.0)GO TO 37 ANM=X01+IN(1,I)*PSN/IN(3,I) DEN=Y01+IN(2,I)*Q01/IN(3,I)+IN(1,I)*PCN/IN(3,I) IF(ABS(ANM).LT.1E-7.AND.ABS(DEN).LT.1E-7)GO TO 33 IF(ABS(ANM).LT.1E-7)GO TO 34 IF(ABS(DEN).LT.1E-7)GO TO 35 FC(I)=ATAN2(ANM,DEN) D1=SIN(FC(I)) IF(ABS(D1).LT.1E-7)GO TO 36 RC(I)=ATAN2(ANM,D1) GO TO 42 33 FC(I)=0. RC(I)=0. GO TO 43 34 FC(I)=0. RC(I)=ATAN(ABS(DEN)) GO TO 43 35 IF(ANM.GT.0.)FC(I)=P2 IF(ANM.LT.0.)FC(I)=3.*P2 RC(I)=ATAN(ABS(ANM)) GO TO 42 36 D1=COS(FC(I)) RC(I)=ATAN2(DEN,D1) GO TO 42 37 IF(IN(1,I).EQ.0)GO TO 39 IF(IN(2,I).EQ.0)GO TO 40 DEN=IN(2,I)*Q01/(IN(1,I)*P01)+CSN IF(ABS(DEN).LT.1E-7)GO TO 38 FC(I)=ATAN2(SNN,DEN) IF(IN(1,I).LT.0)FC(I)=FC(I)+2.*P2 IF(FC(I).GE.(4.*P2))FC(I)=FC(I)-(4.*P2) GO TO 41 38 FC(I)=P2 GO TO 41 39 IF(IN(2,I).GT.0)FC(I)=0. IF(IN(2,I).LT.0)FC(I)=2.*P2 GO TO 41 40 FC(I)=AR(3) IF(IN(1,I).LT.0)FC(I)=AR(3)+2.*P2 41 RC(I)=P2 42 IF(FC(I).LT.0)FC(I)=4.*P2+FC(I) 43 CONTINUE FI1=AR(3) DO 54 I=1,3 A(I)=A(I)/RD SA(I)=SA(I)/RD AR(I)=AR(I)/RD 54 SAR(I)=SAR(I)/RD WRITE(IT,113) 113 FORMAT(//10X,'Morphological parameters from optical data') IF(KS.EQ.0)GO TO 44 WRITE(IT,100)C(1),SC(1),(A(I),SA(I),I=1,2),C(3),SC(3),A(3), 1 SA(3) 100 FORMAT(/6X,'a/b =',F8.5,'(',F7.5,')',7X,'alpha =',F8.3,'(', 1 F5.3,')'/6X,'b/b = 1',22X,'beta =',F8.3,'(',F5.3,')'/6X, 2 'c/b =',F8.5,'(',F7.5,')',7X,'gamma =',F8.3,'(',F5.3,')'/) WRITE(IT,101)P0,SP0,AR(1),SAR(1),Q0,SQ0,(AR(I),SAR(I),I=2,3) 101 FORMAT(7X,'p0 =',F8.5,'(',F7.5,')',7X,'lambda=',F8.3,'(', 1 F5.3,')'/7X,'q0 =',F8.5,'(',F7.5,')',7X,'mu',4X,'=',F8.3, 2 '(',F5.3,')'/7X,'r0 = 1',22X,'nu',4X,'=',F8.3,'(',F5.3,')'/) WRITE(IT,102)P01,SP01,X01,SX01,Q01,SQ01,Y01,SY01 102 FORMAT(7X,'p0''=',F8.5,'(',F7.5,')',10X,'x0''=',F9.5,'(',F7.5, 1 ')'/7X,'q0''=',F8.5,'(',F7.5,')',10X,'y0''=',F9.5,'(',F7.5,')'/) GO TO 45 44 WRITE(IT,103)C(1),(A(I),I=1,2),C(3),A(3) 103 FORMAT(/6X,'a/b =',F8.5,16X,'alpha =',F8.3/6X,'b/b = 1',22X, 1 'beta =',F8.3/6X,'c/b =',F8.5,16X,'gamma =',F8.3/) WRITE(IT,104)P0,AR(1),Q0,(AR(I),I=2,3) 104 FORMAT(7X,'p0 =',F8.5,16X,'lambda=',F8.3/7X,'q0 =',F8.5,16X, 1 'mu',4X,'=',F8.3/7X,'r0 = 1',22X,'nu',4X,'=',F8.3/) WRITE(IT,105)P01,X01,Q01,Y01 105 FORMAT(7X,'p0''=',F8.5,19X,'x0''=',F9.5/7X,'q0''=',F8.5,19X, 1 'y0''=',F9.5/) 45 WRITE(IT,106) 106 FORMAT(/11X,'Coordinates of gnomonic projection poles'// 1 8X,'h',2X,'k',2X,'l',9X,'x''',18X,'y''') DO 48 I=1,N IF(ABS(RO(I)-P2).LE.1E-4.OR.IN(3,I).LT.0)GO TO 47 IF(KS.EQ.0)GO TO 46 WRITE(IT,107)(IN(J,I),J=1,3),X1(I),SX(I),Y1(I),SY(I) 107 FORMAT(6X,3I3,1X,2(2X,F9.5,'(',F7.5,')')) GO TO 48 46 WRITE(IT,108)(IN(J,I),J=1,3),X1(I),Y1(I) 108 FORMAT(6X,3I3,3X,F9.5,11X,F9.5) GO TO 48 47 WRITE(IT,109)(IN(J,I),J=1,3) 109 FORMAT(6X,3I3,9X,'-',19X,'-') 48 CONTINUE WRITE(IT,110) 110 FORMAT(//2X,'Stereographic projection coordinates and ', 1 'interfacial angles'//4X,'Faces',5X,'Observed',7X,'Calculated', 2 7X,'angles with'/3X,'h',2X,'k',2X,'l',2(4X,'Phi', 3 4X,'Ro',2X),3X,'(100)',2X,'(010)',2X,'(001)') DO 51 I=1,N IF(RC(I).GT.(2.*P2))RC(I)=RC(I)-2.*P2 FIO=FI(I)/RD ROO=RO(I)/RD FIC=FC(I)/RD ROC=RC(I)/RD IF(ROC.LT.0) ROC=ROC+180. IF(ABS(FIO-360.).LT..01)FIO=0. IF(ABS(FIC-360.).LT..01)FIC=0. SNR=SIN(RO(I)) CSR=COS(RO(I)) CSA=SNR*COS(FI1-FI(I)) CSB=SNR*COS(FI(I)) CSC=CSR*CRO0+SNR*SRO0*COS(FI0-FI(I)) AA=1.-CSA**2 BB=1.-CSB**2 CC=1.-CSC**2 IF(AA.LT.0.)AA=0. IF(BB.LT.0.)BB=0. IF(CC.LT.0.)CC=0. SNA=SQRT(AA) SNB=SQRT(BB) SNC=SQRT(CC) D1=ATAN2(SNA,CSA)/RD D2=ATAN2(SNB,CSB)/RD D3=ATAN2(SNC,CSC)/RD DIF=ROO-ROC IF(ABS(DIF).LT.10) GO TO 2 IF(DIF.GT.0) ROC=180-ROC IF(DIF.LT.0) ROC=180+ROC 2 IF(IN(3,I).GT.0)GO TO 1 DIF=FIO-FIC IF(ABS(DIF).LT.10) GO TO 1 IF(DIF.GT.0) FIC=FIC+180 IF(DIF.LT.0) FIC=FIC-180 1 IF(ABS(ROO).LT.1E-3)GO TO 49 GO TO 50 49 WRITE(IT,111)(IN(J,I),J=1,3),ROO,ROC,D1,D2,D3 111 FORMAT(1X,3I3,5X,'-',3X,F6.2,5X,'-',3X,F6.2,1X,3(1X,F6.2)) GO TO 51 50 WRITE(IT,112)(IN(J,I),J=1,3),FIO,ROO,FIC,ROC,D1,D2,D3 112 FORMAT(1X,3I3,2(1X,F7.2,1X,F6.2),1X,3(1X,F6.2)) 51 CONTINUE RETURN END SUBROUTINE COSTER(R,KDC,IM,NI,X01,Y01,PSN,PCN,Q01) C-----This routine calculates the coordinates for the C-----stereographic projection of the crystal faces from C-----the X-ray crystal data DIMENSION SN(3),CS(3),SNR(3),CSR(3),NH(3),AG(3),ARG(3) DIMENSION SAR(3),AB(3) COMMON/LBL1/C(3),A(3),SC(3),SA(3) 1 /LBL2/ID(3,100) 2 /LBL4/FC(100),RC(100) 3 /GRP5/IND(3,17) DATA IN,IOT/13,5/ RD=0.0174532925 P2=90.*RD DO 1 I=1,3 AG(I)=A(I) A(I)=A(I)*RD SA(I)=SA(I)*RD ARG(I)=(SC(I)/C(I))**2 SN(I)=SIN(A(I)) CS(I)=COS(A(I)) 1 CONTINUE C-----Calculation of reciprocal lattice angles IF(ABS(CS(2)).LT.1E-7.OR.ABS(CS(3)).LT.1E-7) GO TO 20 AA=CS(2)*CS(3) GO TO 21 20 AA=0. 21 CSR(1)=(AA-CS(1))/(SN(2)*SN(3)) IF(ABS(CS(1)).LT.1E-7.OR.ABS(CS(3)).LT.1E-7) GO TO 22 AA=CS(1)*CS(3) GO TO 23 22 AA=0. 23 CSR(2)=(AA-CS(2))/(SN(1)*SN(3)) IF(ABS(CS(1)).LT.1E-7.OR.ABS(CS(2)).LT.1E-7) GO TO 24 AA=CS(1)*CS(2) GO TO 25 24 AA=0. 25 CSR(3)=(AA-CS(3))/(SN(1)*SN(2)) DO 2 I=1,3 IF(ABS(CSR(I)).LT.1E-7)GO TO 14 SNR(I)=SQRT(1.-CSR(I)**2) GO TO 2 14 CSR(I)=0. SNR(I)=1. 2 CONTINUE IF(SA(1).LT.1E-7)GO TO 62 IF(ABS(CSR(1)).LT.1E-7.OR.(ABS(CS(2)).LT.1E-7.AND. 1 ABS(CS(3)).LT.1E-7)) GO TO 26 IF(ABS(CS(2)).LT.1E-7.AND.ABS(CS(3)).GE.1E-7) GO TO 27 IF(ABS(CS(2)).GE.1E-7.AND.ABS(CS(3)).LT.1E-7) GO TO 28 AA=CSR(1)*CS(2) BB=CSR(1)*CS(3) GO TO 29 26 AA=0. BB=0. GO TO 29 27 AA=0. BB=CSR(1)*CS(3) GO TO 29 28 AA=CSR(1)*CS(2) BB=0. 29 SAR(1)=SQRT(((CS(3)/SN(3)+AA/SN(2))*SA(2))**2+ 1 ((CS(2)/SN(2)+BB/SN(3))*SA(3))**2+ 2 (SN(1)/(SN(2)*SN(3))*SA(1))**2)/SNR(1) GO TO 63 62 SAR(1)=0. 63 IF(SA(2).LT.1E-7) GO TO 64 IF(ABS(CSR(2)).LT.1E-7.OR.(ABS(CS(3)).LT.1E-7.AND. 1 ABS(CS(1)).LT.1E-7)) GO TO 30 IF(ABS(CS(3)).LT.1E-7.AND.ABS(CS(1)).GE.1E-7) GO TO 31 IF(ABS(CS(3)).GE.1E-7.AND.ABS(CS(1)).LT.1E-7) GO TO 32 AA=CSR(2)*CS(3) BB=CSR(2)*CS(1) GO TO 33 30 AA=0. BB=0. GO TO 33 31 AA=0. BB=CSR(2)*CS(1) GO TO 33 32 AA=CSR(2)*CS(3) BB=0. 33 SAR(2)=SQRT(((CS(1)/SN(1)+AA/SN(3))*SA(3))**2+ 1 ((CS(3)/SN(3)+BB/SN(1))*SA(1))**2+ 2 (SN(2)/(SN(1)*SN(3))*SA(2))**2)/SNR(2) GO TO 65 64 SAR(2)=0. 65 IF(SA(3).LT.1E-7) GO TO 66 IF(ABS(CSR(3)).LT.1E-7.OR.(ABS(CS(1)).LT.1E-7.AND. 1 ABS(CS(2)).LT.1E-7)) GO TO 34 IF(ABS(CS(1)).LT.1E-7.AND.ABS(CS(2)).GE.1E-7) GO TO 35 IF(ABS(CS(1)).GE.1E-7.AND.ABS(CS(2)).LT.1E-7) GO TO 36 AA=CSR(3)*CS(1) BB=CSR(3)*CS(2) GO TO 37 34 AA=0. BB=0. GO TO 37 35 AA=0. BB=CSR(3)*CS(2) GO TO 37 36 AA=CSR(3)*CS(1) BB=0. 37 SAR(3)=SQRT(((CS(2)/SN(2)+AA/SN(1))*SA(1))**2+ 1 ((CS(1)/SN(1)+BB/SN(2))*SA(2))**2+ 2 (SN(3)/(SN(1)*SN(2))*SA(3))**2)/SNR(3) GO TO 67 66 SAR(3)=0. 67 CONTINUE C-----Calculation of x0 & y0 Y0=CSR(1) IF(ABS(CSR(2)).LT.1E-7.AND.ABS(CSR(3)).LT.1E-7. 1 AND.ABS(Y0).LT.1E-7)GO TO 15 IF(ABS(Y0).LT.1E-7.OR.ABS(CSR(3)).LT.1E-7)GO TO 3 X0=(CSR(2)-Y0*CSR(3))/SNR(3) GO TO 68 3 X0=CSR(2)/SNR(3) 68 CONTINUE GO TO 16 15 X0=0. 16 CONTINUE SY0=SNR(1)*SAR(1) AA=0. IF(SAR(2).GT.1E-7)AA=(SNR(2)/SNR(3)*SAR(2))**2 BB=0. IF(SY0.GT.1E-7.AND.ABS(CSR(3)).GT.1E-7)BB=(CSR(3)/SNR(3)*SY0)**2 IF(SAR(3).LE.1E-7.OR.(X0.LE.1E-7.AND. 1 (Y0.LE.1E-7.OR.ABS(CSR(3)).LE.1E-7)))CC=0. IF(SAR(3).GT.1E-7.AND.(ABS(Y0).GT.1E-7.OR.(ABS(X0).GT.1E-7 1 .AND.ABS(CSR(3)).GT.1E-7)))CC=((Y0-X0*CSR(3)/SNR(3)) 2 *SAR(3))**2 SX0=SQRT(AA+BB+CC) C-----Calculation of Phi0 & Ro0 IF(ABS(X0).LT.1E-7.AND.ABS(Y0).LT.1E-7)GO TO 4 IF(ABS(X0).LT.1E-7)GO TO 5 IF(ABS(Y0).LT.1E-7)GO TO 6 FI0=ATAN2(X0,Y0) SF0=SIN(FI0) CF0=COS(FI0) IF(ABS(SF0).LT.1E-7)GO TO 5 SR0=X0/SF0 CR0=SQRT(1.-SR0**2) RO0=ATAN2(SR0,CR0) IF(SX0.LE.1E-7.OR.ABS(X0).LE.1E-7)AA=0. IF(SX0.GT.1E-7.OR.ABS(X0).GT.1E-7)AA=SX0/X0 IF(SY0.LE.1E-7.OR.ABS(Y0).LE.1E-7)BB=0. IF(SY0.GT.1E-7.OR.ABS(Y0).GT.1E-7)BB=SY0/Y0 IF(ABS(AA).LT.1E-7)GO TO 44 AA=AA*AA GO TO 45 44 AA=0. 45 IF(ABS(BB).LT.1E-7)GO TO 46 BB=BB*BB GO TO 47 46 BB=0. 47 IF(ABS(SF0).LE.1E-7.OR.ABS(CF0).LE.1E-7.OR. 1 SQRT(AA+BB).LE.1E-7)GO TO 38 SGF0=ABS(SF0*CF0*SQRT(AA+BB)) GO TO 49 38 SGF0=0. IF(ABS(SF0).LT.1E-7)GO TO 69 49 IF(ABS(CF0).LE.1E-7.OR.SGF0.LE.1E-7)BB=0. IF(ABS(CF0).GT.1E-7.AND.SGF0.GT.1E-7)BB=(CF0/SF0*SGF0)**2 51 SGR0=ABS(SR0/CR0*SQRT(AA+BB)) GO TO 70 69 SGR0=0. 70 CONTINUE GO TO 7 4 X0=0. Y0=0. SGR0=0. SGF0=0. RO0=0. FI0=0. GO TO 7 5 FI0=0. SF0=SIN(FI0) CF0=COS(FI0) RO0=ASIN(Y0) CR0=COS(RO0) SGF0=0. SGR0=0. IF(ABS(CR0).GT.1E-7.AND.SY0.GT.1E-7)SGR0=SY0/CR0 GO TO 7 6 FI0=90.*RD SF0=SIN(FI0) CF0=COS(FI0) RO0=ASIN(X0) CR0=COS(RO0) SGF0=0. SGR0=0. IF(ABS(CR0).GT.1E-7.AND.SX0.GT.1E-7)SGR0=SX0/CR0 7 CONTINUE KDR=0 SR0=SIN(RO0) CR0=COS(RO0) C-----Calculation of p0, q0, p0', q0', x0', y0' IF(ABS(A(1)-90.*RD).LT.1E-6.AND.ABS(A(2)-90.*RD) 1 .LT.1E-6.AND.ABS(A(3)-90.*RD).LT.1E-6)GO TO 17 GO TO 18 17 X01=0. Y01=0. P0=C(3)/C(1) Q0=C(3)/C(2) P01=P0 Q01=Q0 SX01=0. SY01=0. SP0=P0*SQRT(ARG(1)+ARG(3)) SQ0=Q0*SQRT(ARG(2)+ARG(3)) SP01=SP0 SQ01=SQ0 KDR=1 GO TO 71 18 CONTINUE TGR=SIN(RO0)/COS(RO0) X01=SIN(FI0)*TGR Y01=COS(FI0)*TGR D1=CR0*SNR(3) Q01=(C(3)*SNR(2))/(C(2)*D1) P01=(C(3)*SNR(1))/(C(1)*D1) P0=(C(3)*SN(1))/(C(1)*SN(3)) Q0=(C(3)*SN(2))/(SN(3)*C(2)) 71 CONTINUE DO 52 I=1,3 IF(ABS(CS(I)).LT.1E-7.OR.SA(I).LT.1E-7)GO TO 53 AB(I)=(CS(I)/SN(I)*SA(I))**2 GO TO 52 53 AB(I)=0. 52 CONTINUE SRC1=(C(1)/C(2))*SQRT(ARG(1)+ARG(2)) SRC2=(C(3)/C(2))*SQRT(ARG(3)+ARG(2)) IF(KDR.EQ.1)GO TO 72 SP0=P0*SQRT(ARG(1)+ARG(3)+AB(1)+AB(3)) SQ0=Q0*SQRT(ARG(2)+ARG(3)+AB(2)+AB(3)) DO 54 I=1,3 IF(ABS(CSR(I)).LT.1E-7.OR.SAR(I).LT.1E-7)GO TO 55 AB(I)=(CSR(I)/SNR(I)*SAR(I))**2 GO TO 54 55 AB(I)=0. 54 CONTINUE IF(ABS(CR0).LT.1E-7.OR.ABS(SR0).LT.1E-7.OR. 1 SGR0.LT.1E-7)GO TO 56 AA=(SR0/CR0*SGR0)**2 GO TO 57 56 AA=0. 57 SP01=P01*SQRT(ARG(1)+ARG(3)+AB(1)+AA+AB(3)) SQ01=Q01*SQRT(ARG(2)+ARG(3)+AB(2)+AA+AB(3)) AA=0. IF(ABS(CF0).GT.1E-7.AND.ABS(SR0).GT.1E-7.AND.SGF0.GT.1E-7. 1 AND.CR0.GT.1E-7)AA=(CF0*SR0/CR0*SGF0)**2 BB=0. IF(ABS(SF0).GT.1E-7.AND.ABS(CR0).GT.1E-7.AND.SGR0.GT.1E-7) 1 BB=(SF0/CR0**2*SGR0)**2 SX01=SQRT(AA+BB) AA=0. IF(ABS(FI0).GT.1E-7.AND.ABS(RO0).GT.1E-7.AND.SGF0.GT.1E-7) 1 AA=(SF0*SR0/CR0*SGF0)**2 BB=0. IF(ABS(FI0-P2).GT.1E-7.AND.SGR0.GT.1E-7)BB=((CF0+Y01*SR0)/CR0* 1 SGR0)**2 SY01=SQRT(AA+BB) 72 CONTINUE SNNU=SNR(3) CSNU=CSR(3) RS=C(2) DO 19 I=1,3 ARG(I)=ATAN2(SNR(I),CSR(I))/RD C(I)=C(I)/RS SA(I)=SA(I)/RD SAR(I)=SAR(I)/RD 19 CONTINUE WRITE(IOT,106) 106 FORMAT(//8X,'Morphological parameters from X-ray crystal data'/) WRITE(IOT,101)C(1),SRC1,AG(1),SA(1),AG(2),SA(2),C(3),SRC2,AG(3), 1 SA(3) 101 FORMAT(7X,'a/b =',F8.5,'(',F7.5,')',7X,'alpha =',F8.3,'(', 1 F5.3,')'/7X,'b/b = 1',22X,'beta =',F8.3,'(',F5.3,')'/ 2 7X,'c/b =',F8.5,'(',F7.5,')',7X,'gamma =',F8.3,'(',F5.3, 3 ')'/) WRITE(IOT,102)P0,SP0,ARG(1),SAR(1),Q0,SQ0,ARG(2),SAR(2), 1 ARG(3),SAR(3) 102 FORMAT(8X,'p0 =',F8.5,'(',F7.5,')',7X,'lambda=',F8.3,'(', 1 F5.3,')'/8X,'q0 =',F8.5,'(',F7.5,')',7X,'mu',4X,'=',F8.3, 2 '(',F5.3,')'/8X,'r0 = 1',22X,'nu',4X,'=',F8.3,'(',F5.3,')'/) WRITE(IOT,103)P01,SP01,X01,SX01,Q01,SQ01,Y01,SY01 103 FORMAT(8X,'p0''=',F8.5,'(',F7.5,')',10X,'x0''=',F9.5,'(', 1 F7.5,')'/8X,'q0''=',F8.5,'(',F7.5,')',10X,'y0''=',F9.5,'(', 2 F7.5,')'//) WRITE(IOT,100) 100 FORMAT(1X,'Parameters & coordinates for gnomonic & ', 1 'stereographic projections'/) WRITE(IOT,104) 104 FORMAT(4X,'Faces',3X,'Gnomonic coordinates',3X,'Stereographic ', 1 'coordinates'/3X,'h',2X,'k',2X,'l',8X,'x''',6X,'y''',7X,'Phi' 2 ,5X,'Ro',7X,'xs',5X,'ys') C-----Set up faces NFC=0 NSF=0 IF(KDC.EQ.0) GO TO 8 IF(KDC.EQ.1) GO TO 10 8 IMP=IM+1 KM=2*IM+1 DO 9 I=1,KM IH=I-IM-1 DO 9 J=1,KM IK=J-IM-1 DO 9 K=1,KM IL=K-IM-1 C-----Reject faces IF(IH.EQ.0.AND.IK.EQ.0.AND.IL.EQ.0)GO TO 9 IF(ABS(IH).NE.1.AND.IK.EQ.0.AND.IL.EQ.0)GO TO 9 IF(IH.EQ.0.AND.ABS(IK).NE.1.AND.IL.EQ.0)GO TO 9 IF(IH.EQ.0.AND.IK.EQ.0.AND.ABS(IL).NE.1)GO TO 9 IF(IL.EQ.0.AND.ABS(IH).EQ.ABS(IK).AND.ABS(IH).NE.1)GO TO 9 IF(IK.EQ.0.AND.ABS(IH).EQ.ABS(IL).AND.ABS(IH).NE.1)GO TO 9 IF(IH.EQ.0.AND.ABS(IK).EQ.ABS(IL).AND.ABS(IK).NE.1)GO TO 9 IF(ABS(IH).EQ.ABS(IK).AND.ABS(IH).EQ.ABS(IL).AND.ABS(IH) 1 .NE.1)GO TO 9 NH(1)=IH NH(2)=IK NH(3)=IL C-----Calculation of gnomonic and stereographic coordinates CALL PRGJT(NH,X01,Y01,P01,Q01,RO0,FI0,SNNU,CSNU,R,FN,RN) NFC=NFC+1 IF(NH(3).LT.0)GO TO 9 DO 39 L=1,3 IF(ABS(NH(L)).GT.1)GO TO 9 39 CONTINUE NSF=NSF+1 FC(NSF)=FN RC(NSF)=RN DO 40 L=1,3 40 IND(L,NSF)=NH(L) 9 CONTINUE GO TO 13 10 DO 11 I=1,NI DO 12 J=1,3 12 NH(J)=ID(J,I) CALL PRGJT(NH,X01,Y01,P01,Q01,RO0,FI0,SNNU,CSNU,R,FN,RN) NFC=NFC+1 FC(I)=FN RC(I)=RN 11 CONTINUE 13 CONTINUE WRITE(IOT,105)NFC 105 FORMAT(/1X,'Number of faces',I6) PSN=P01*SNNU PCN=P01*CSNU RETURN END SUBROUTINE CORNER(RZ,RY) C-----This routine calculates the coordinates of the vertices of C-----a crystal, referred to a cartesian system of reference with C-----the origin at a point inside the crystal, using the two-circle C-----goniometer Phi & Ro readings DIMENSION CD(3,3),R(3,3),DT(3),X(3) COMMON/LBL2/IN(3,100) 1 /LBL3/FI(100),RO(100) 2 /LBL5/XV(3,200),XVR(3,200) 3 /LBL6/IV(3,200),IE(2,300) 4 /LBL8/DST(100) 5 /GRP1/N,NV,NE DATA INN,IT/13,5/ RD=0.0174532925 IF(ABS(RZ).LT.1E-7.AND.ABS(RY).LT.1E-7)GO TO 1 RZ=RZ*RD RY=RY*RD GO TO 2 1 RZ=0.3217165 RY=0.1652303 2 CONTINUE C-----Set up the rotation matrix R(1,2)=SIN(RZ) R(2,2)=COS(RZ) R(3,1)=-SIN(RY) R(3,3)=COS(RY) R(1,1)=R(2,2)*R(3,3) R(1,3)=-R(2,2)*R(3,1) R(2,1)=-R(1,2)*R(3,3) R(2,3)=R(1,2)*R(3,1) R(3,2)=0. WRITE(IT,100) 100 FORMAT(//' Coordinates of the vertices ', 1 '(x,y,z are referred to a cartesian'/' system of ', 2 'reference with the origin inside the crystal)'//4X,'N ', 3 'h1 k1 l1 h2 k2 l2 h3 k3 l3',6X,'x',7X,'y',7X,'z') C-----Calculation of corner coordinates NV=0 NM2=N-2 NM1=N-1 DO 10 I=1,NM2 IP1=I+1 DO 10 J=IP1,NM1 JP1=J+1 DO 10 K=JP1,N SNR=SIN(RO(I)) CD(1,1)=SNR*SIN(FI(I)) CD(2,1)=SNR*COS(FI(I)) CD(3,1)=COS(RO(I)) DT(1)=DST(I) SNR=SIN(RO(J)) CD(1,2)=SNR*SIN(FI(J)) CD(2,2)=SNR*COS(FI(J)) CD(3,2)=COS(RO(J)) DT(2)=DST(J) SNR=SIN(RO(K)) CD(1,3)=SNR*SIN(FI(K)) CD(2,3)=SNR*COS(FI(K)) CD(3,3)=COS(RO(K)) DT(3)=DST(K) IF(ABS(CD(1,1)).LE.1E-5.AND.ABS(CD(1,2)).LE.1E-5.AND. 1 ABS(CD(1,3)).LE.1E-5)GO TO 10 IF(ABS(CD(2,1)).LE.1E-5.AND.ABS(CD(2,2)).LE.1E-5.AND. 1 ABS(CD(2,3)).LE.1E-5)GO TO 10 IF(ABS(CD(3,1)).LE.1E-5.AND.ABS(CD(3,2)).LE.1E-5.AND. 1 ABS(CD(3,3)).LE.1E-5)GO TO 10 IF(ABS(CD(1,1)).LE.1E-5.AND.ABS(CD(1,2)).GT.1E-5)GO TO 3 IF(ABS(CD(1,1)).LE.1E-5.AND.ABS(CD(1,3)).GT.1E-5)GO TO 5 GO TO 7 3 DO 4 L=1,3 RS=CD(L,1) CD(L,1)=CD(L,2) 4 CD(L,2)=RS RS=DT(1) DT(1)=DT(2) DT(2)=RS GO TO 7 5 DO 6 L=1,3 RS=CD(L,1) CD(L,1)=CD(L,3) 6 CD(L,3)=RS RS=DT(1) DT(1)=DT(3) DT(3)=RS 7 R1=CD(1,2)/CD(1,1) R2=CD(1,3)/CD(1,1) A=R1*CD(2,1)-CD(2,2) B=R2*CD(3,1)-CD(3,3) C=R1*CD(3,1)-CD(3,2) D=R2*CD(2,1)-CD(2,3) E=R1*DT(1)-DT(2) F=R2*DT(1)-DT(3) DN=A*B-C*D IF(ABS(DN).LT.1E-7)GO TO 10 X(2)=(E*B-C*F)/DN X(3)=(A*F-D*E)/DN X(1)=(DT(1)-CD(2,1)*X(2)-CD(3,1)*X(3))/CD(1,1) C-----Reject corners falling outside the crystal DO 8 L=1,N SNR=SIN(RO(L)) A=SNR*SIN(FI(L)) B=SNR*COS(FI(L)) C=COS(RO(L)) D=DST(L)-(X(1)*A+X(2)*B+X(3)*C) IF(D.LT.-0.0001)GO TO 10 8 CONTINUE NV=NV+1 C-----Rotate the crystal DO 9 L=1,3 XV(L,NV)=X(L) XVR(L,NV)=0 DO 9 LL=1,3 9 XVR(L,NV)=XVR(L,NV)+R(L,LL)*X(LL) IV(1,NV)=I IV(2,NV)=J IV(3,NV)=K WRITE(IT,101)NV,(IN(L,I),L=1,3),(IN(L,J),L=1,3), 1 (IN(L,K),L=1,3),(XV(L,NV),L=1,3) 101 FORMAT(2X,I3,3(1X,3I3),2X,3F8.4) 10 CONTINUE RZ=RZ/RD RY=RY/RD WRITE(IT,106)RZ,RY 106 FORMAT(//' Coordinates of the vertices after rotation of ', 1 'the crystal'/' by',F6.2,' degrees around z and by',F6.2, 2 ' degrees around y'//14X,'N',8X,'xr',8X,'yr',8X,'zr') DO 18 I=1,NV 18 WRITE(IT,107)I,(XVR(L,I),L=1,3) 107 FORMAT(12X,I3,2X,3F10.4) C-----Define crystal edges NE=0 NVM=NV-1 DO 15 I=1,NVM IP1=I+1 DO 15 J=IP1,NV KDL=0 DO 11 L=1,3 DO 11 LL=1,3 11 IF(IV(L,I).EQ.IV(LL,J))KDL=KDL+1 IF(KDL.EQ.2)GO TO 12 GO TO 15 12 IF(NE.EQ.0)GO TO 14 DO 13 K=1,NE IF((IE(1,K).EQ.I.AND.IE(2,K).EQ.J).OR. 1 (IE(1,K).EQ.J.AND.IE(2,K).EQ.I))GO TO 15 13 CONTINUE 14 NE=NE+1 IE(1,NE)=I IE(2,NE)=J 15 CONTINUE WRITE(IT,103)N,NV,NE 103 FORMAT(/7X,' The Euler''s criterion for convex polyhedra:'/ 1 8X,I3,' faces + ',I3,' vertices = ',I3,' edges + 2') IF((N+NV).NE.(NE+2))GO TO 16 WRITE(IT,104) 104 FORMAT(7X,' is fulfilled') GO TO 17 16 WRITE(IT,105) 105 FORMAT(7X,' is not fulfilled') STOP 17 RETURN END SUBROUTINE CNTR C-----This routine indicates, for a given set of corners of a C-----crystal, if the edges must be drawn with solid or dotted C-----lines, and gives the coordinates of the centres of the C-----faces in the drawings DIMENSION IC(100),LE(3),ANG(6),LV(2,6) COMMON/LBL5/XV(3,200),XVR(3,200) 1 /LBL6/IV(3,200),IE(2,300) 2 /LBL9/MSD(300),CCF(2,100) 3 /GRP1/N,NV,NE IO=5 C-----Choose the corner most distant from the center DM=0. DO 1 I=1,NV D=SQRT(XVR(2,I)**2+XVR(3,I)**2) IF(D.LT.DM)GO TO 1 DM=D IMX=I 1 CONTINUE NC=1 IC(NC)=IMX C-----Choose the edges making the largest polygon around center N1=0 DO 2 I=1,NE IF(IC(1).NE.IE(1,I).AND.IC(1).NE.IE(2,I))GO TO 2 N1=N1+1 LE(N1)=IE(1,I) IF(IC(1).EQ.IE(1,I))LE(N1)=IE(2,I) 2 CONTINUE N1M=N1-1 N2=0 DO 3 I=1,N1M DO 3 J=1,N1 Y2=XVR(2,LE(I)) Z2=XVR(3,LE(I)) Y1=XVR(2,IC(1)) Z1=XVR(3,IC(1)) Y3=XVR(2,LE(J)) Z3=XVR(3,LE(J)) CALL PLANG(Y1,Z1,Y2,Z2,Y3,Z3,AG) N2=N2+1 ANG(N2)=AG LV(1,N2)=LE(I) LV(2,N2)=LE(J) 3 CONTINUE DM=0. DO 4 I=1,N2 IF(ANG(I).LT.DM)GO TO 4 DM=ANG(I) MX1=LV(1,I) MX2=LV(2,I) 4 CONTINUE D1=SQRT(XVR(2,MX1)**2+XVR(3,MX1)**2) D2=SQRT(XVR(2,MX2)**2+XVR(3,MX2)**2) IF(D1.GE.D2)IC(2)=MX1 IF(D1.LT.D2)IC(2)=MX2 NC=2 6 CONTINUE N1=0 DO 10 I=1,NE IF(IC(NC).NE.IE(1,I).AND.IC(NC).NE.IE(2,I))GO TO 10 IF(IE(1,I).EQ.IC(NC-1).OR.IE(2,I).EQ.IC(NC-1))GO TO 10 N1=N1+1 IF(IC(NC).EQ.IE(1,I))GO TO 8 Y3=XVR(2,IE(1,I)) Z3=XVR(3,IE(1,I)) LE(N1)=IE(1,I) GO TO 9 8 Y3=XVR(2,IE(2,I)) Z3=XVR(3,IE(2,I)) LE(N1)=IE(2,I) 9 Y1=XVR(2,IC(NC)) Z1=XVR(3,IC(NC)) Y2=XVR(2,IC(NC-1)) Z2=XVR(3,IC(NC-1)) CALL PLANG(Y1,Z1,Y2,Z2,Y3,Z3,AG) ANG(N1)=AG 10 CONTINUE DM=0 DO 11 I=1,N1 IF(ANG(I).LT.DM)GO TO 11 DM=ANG(I) IMX=LE(I) 11 CONTINUE IF(IMX.EQ.IC(1))GO TO 31 IF(NC.GT.NV)GO TO 7 NC=NC+1 IC(NC)=IMX GO TO 6 7 WRITE(IO,100) 100 FORMAT(//39H THE EDGES DO NOT FORM A CLOSED CONTOUR) STOP 31 CONTINUE C-----Assign codes to edges (for the drawing) DO 30 I=1,NE DO 24 J=1,NC IF(IE(1,I).EQ.IC(J))GO TO 25 24 CONTINUE ID1=0 GO TO 26 25 ID1=1 26 DO 27 J=1,NC IF(IE(2,I).EQ.IC(J))GO TO 28 27 CONTINUE ID2=0 GO TO 29 28 ID2=1 29 MSD(I)=1H+ IF(ID1.EQ.0.AND.ID2.EQ.1.AND.XVR(1,IE(1,I)).LT.0.)MSD(I)=1H- IF(ID1.EQ.1.AND.ID2.EQ.0.AND.XVR(1,IE(2,I)).LT.0.)MSD(I)=1H- IF(ID1.EQ.0.AND.ID2.EQ.0.AND.XVR(1,IE(1,I)).LT.0..AND. 1 XVR(1,IE(2,I)).LT.0.)MSD(I)=1H- 30 CONTINUE C-----Find the centers of the faces in the drawing DO 18 I=1,N NVF=0 DO 12 L=1,2 12 CCF(L,I)=0. DO 16 J=1,NV DO 13 K=1,3 IF(IV(K,J).EQ.I)GO TO 14 13 CONTINUE GO TO 16 14 NVF=NVF+1 DO 15 L=1,2 15 CCF(L,I)=CCF(L,I)+XVR(L+1,J) 16 CONTINUE DO 17 L=1,2 17 CCF(L,I)=CCF(L,I)/NVF 18 CONTINUE C-----Print edges WRITE(IO,101) 101 FORMAT(//' Edges (a (+) means that the edge should be ', 1 'drawn by a solid'/' line, a (-) by a dotted line)'/) NR=MOD(NE,4) NM=NE-NR DO 19 I=1,NM,4 19 WRITE(IO,102)(IE(J,I),J=1,2),MSD(I),(IE(J,I+1),J=1,2), 1 MSD(I+1),(IE(J,I+2),J=1,2),MSD(I+2),(IE(J,I+3),J=1,2), 2 MSD(I+3) 102 FORMAT(1X,4(5X,I2,'-',I2,' (',A1,')')) NR=NR+1 GO TO(23,20,21,22),NR 20 WRITE(IO,103)(IE(J,NM+1),J=1,2),MSD(NM+1) 103 FORMAT(1X,5X,I2,'-',I2,' (',A1,')') GO TO 23 21 WRITE(IO,104)(IE(J,NM+1),J=1,2),MSD(NM+1),(IE(J,NM+2), 1 J=1,2),MSD(NM+2) 104 FORMAT(1X,2(5X,I2,'-',I2,' (',A1,')')) GO TO 23 22 WRITE(IO,105)(IE(J,NM+1),J=1,2),MSD(NM+1),(IE(J,NM+2), 1 J=1,2),MSD(NM+2),(IE(J,NM+3),J=1,2),MSD(NM+3) 105 FORMAT(1X,3(5X,I2,'-',I2,' (',A1,')')) 23 CONTINUE RETURN END SUBROUTINE PLTCRY(PO,PV,RZ,RY) C-----This routine plots the vertices of a crystal CHARACTER*1 LS(10) DIMENSION IDV(3,100) CHARACTER*1 LN(90) COMMON/LBL5/XV(3,200),XVR(3,200) 1 /GRP1/N,NV,NE DATA LS/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ IO=5 WRITE(IO,100)RZ,RY 100 FORMAT(/' View of the crystal after rotation by',F6.2, 1 ' degrees around z'/' and by',F6.2,' degrees around', 2 ' y (orthogonal axes)'//) YMN=XVR(2,1) ZMN=XVR(3,1) YMX=XVR(2,1) ZMX=XVR(3,1) DO 1 I=1,NV IF(XVR(2,I).LT.YMN)YMN=XVR(2,I) IF(XVR(2,I).GT.YMX)YMX=XVR(2,I) IF(XVR(3,I).LT.ZMN)ZMN=XVR(3,I) IF(XVR(3,I).GT.ZMX)ZMX=XVR(3,I) 1 CONTINUE IF(((YMX-YMN)-(ZMX-ZMN)).LE.0.)GO TO 3 CO=78.*PO/(YMX-YMN) GO TO 4 3 CO=48.*PV/(ZMX-ZMN) 4 CONTINUE DO 5 I=1,NV IDV(1,I)=I IDV(2,I)=((XVR(2,I)-YMN)*CO/PO+.5)+1 5 IDV(3,I)=((XVR(3,I)-ZMN)*CO/PV+.5)+1 NM1=NV-1 DO 7 I=1,NM1 IP1=I+1 DO 6 J=IP1,NV IF(IDV(3,J).LE.IDV(3,I))GO TO 6 DO 18 K=1,3 IR=IDV(K,I) IDV(K,I)=IDV(K,J) IDV(K,J)=IR 18 CONTINUE 6 CONTINUE 7 CONTINUE DO 10 I=1,NM1 IP1=I+1 DO 10 J=IP1,NV IF(IDV(3,I).EQ.IDV(3,J).AND.IDV(2,I).GT.IDV(2,J))GO TO 8 GO TO 10 8 DO 9 K=1,3 IR=IDV(K,I) IDV(K,I)=IDV(K,J) 9 IDV(K,J)=IR 10 CONTINUE IVM=IDV(3,1)+1 DO 17 K=1,IVM DO 11 KL=1,90 11 LN(KL)=1H DO 16 I=1,NV IF((IVM-IDV(3,I)).NE.K)GO TO 16 LN(IDV(2,I))=1H+ IF(LN(IDV(2,I)+1).NE.' '.AND.LN(IDV(2,I)+2).NE.' ')GO TO 16 IDC=IDV(1,I)/10 IIN=IDV(1,I)-IDC*10 DO 12 KK=1,10 IF(IDC.EQ.(KK-1))GO TO 13 12 CONTINUE 13 LN(IDV(2,I)+1)=LS(KK) DO 14 KK=1,10 IF(IIN.EQ.(KK-1))GO TO 15 14 CONTINUE 15 LN(IDV(2,I)+2)=LS(KK) 16 CONTINUE WRITE(IO,101)LN 101 FORMAT(2X,90A1) 17 CONTINUE RETURN END SUBROUTINE PLTSTE(N,PO,PV,KST) C-----This routine plots a stereographic projection DIMENSION IS(2,100),ID(3,100),NO(100) CHARACTER*1 LS(10) CHARACTER*1 LN(90) COMMON/LBL2/IN(3,100) 1 /LBL3/FI(100),RO(100) DATA LS/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ IO=5 NC=(200./PO+0.5)+1 NL=(200./PV+0.5)+1 NS=0 DO 2 I=1,N NO(I)=I IF(KST.EQ.0)GO TO 16 IF(IN(3,I).LT.0)GO TO 2 16 NS=NS+1 DO 1 J=1,3 1 ID(J,NS)=IN(J,I) A=RO(I)/2. R=100.*SIN(A)/COS(A) X=R*SIN(FI(I))+100. Y=R*COS(FI(I))+100. IS(1,NS)=(X/PV+0.5)+1 IS(2,NS)=(Y/PO+0.5)+1 2 CONTINUE NSM=NS-1 DO 4 I=1,NSM IP1=I+1 DO 4 J=IP1,NS IF(IS(1,I).LE.IS(1,J))GO TO 4 IR=IS(1,I) IS(1,I)=IS(1,J) IS(1,J)=IR IR=IS(2,I) IS(2,I)=IS(2,J) IS(2,J)=IR IF(KST.EQ.0)GO TO 17 DO 3 K=1,3 IR=ID(K,I) ID(K,I)=ID(K,J) ID(K,J)=IR 3 CONTINUE GO TO 4 17 IR=NO(I) NO(I)=NO(J) NO(J)=IR 4 CONTINUE DO 6 I=1,NSM IP1=I+1 DO 6 J=IP1,NS IF(IS(1,I).NE.IS(1,J))GO TO 6 IF(IS(2,I).LE.IS(2,J))GO TO 6 IR=IS(2,I) IS(2,I)=IS(2,J) IS(2,J)=IR IF(KST.EQ.0)GO TO 18 DO 5 K=1,3 IR=ID(K,I) ID(K,I)=ID(K,J) ID(K,J)=IR 5 CONTINUE GO TO 6 18 IR=NO(I) NO(I)=NO(J) NO(J)=IR 6 CONTINUE IC=(100./PO+0.5)+1 IL=(100./PV+0.5)+1 DO 15 I=1,NL DO 7 J=1,90 7 LN(J)=1H IF(I.EQ.1)LN(IC)=1H+ IF(I.NE.IL)GO TO 8 LN(1)=1H+ LN(IC)=1H+ LN(NC)=1H+ 8 CONTINUE IF(I.EQ.NL)LN(IC)=1H+ DO 9 K=1,NS IF(IS(1,K).NE.I)GO TO 9 LN(IS(2,K))=1H* 9 CONTINUE DO 14 K=1,NS IF(IS(1,K).NE.I)GO TO 14 IF(KST.EQ.0)GO TO 19 DO 10 J=1,6 IF(LN(IS(2,K)+J).NE.' ')GO TO 14 10 CONTINUE DO 13 J=1,3 IF(ID(J,K).LT.0)LN(IS(2,K)+2*J-1)=1H- DO 11 L=1,10 IF(ABS(ID(J,K)).EQ.(L-1))GO TO 12 11 CONTINUE 12 LN(IS(2,K)+2*J)=LS(L) 13 CONTINUE GO TO 14 19 IF(LN(IS(2,K)+1).NE.' '.OR.LN(IS(2,K)+2).NE.' ')GO TO 14 IDC=NO(K)/10 IIN=NO(K)-IDC*10 DO 20 L=1,10 IF(IDC.EQ.(L-1))GO TO 21 20 CONTINUE 21 LN(IS(2,K)+1)=LS(L) DO 22 L=1,10 IF(IIN.EQ.(L-1))GO TO 23 22 CONTINUE 23 LN(IS(2,K)+2)=LS(L) 14 CONTINUE WRITE(IO,101)LN 101 FORMAT(1X,90A1) 15 CONTINUE RETURN END SUBROUTINE PLTGNO(N,PO,PV,KST,X01,Y01,PSN,PCN,Q01) C-----This routine plots a gnomonic projection DIMENSION NO(100),ID(3,100),X(100),Y(100), 1 IG(2,100),IK(100) CHARACTER*1 LS(10) CHARACTER*1 LN(110) COMMON/LBL2/IN(3,100) 1 /LBL3/FI(100),RO(100) DATA LS/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ IO=5 P2=1.570796327 RMX=0. DO 12 I=1,N IF(RO(I).GT.P2)GO TO 12 IF(ABS(RO(I)-P2).LT.1E-4)GO TO 12 IF(RO(I).GT.RMX)RMX=RO(I) 12 CONTINUE IF(1.1*RMX.GE.P2)RMX=RMX/1.1-.05 NG=0 DO 3 I=1,N NO(I)=I IF(KST.EQ.0)GO TO 1 IF(IN(3,I).LT.0)GO TO 3 1 IF(RO(I)-P2.GT.1E-3)GO TO 3 NG=NG+1 IK(NG)=0 IF(ABS(RO(I)-P2).LT.1E-4)IK(NG)=1 DO 2 J=1,3 2 ID(J,NG)=IN(J,I) X(NG)=0. Y(NG)=0. CSR=COS(RO(I)) IF(ABS(CSR).LT.1E-4)CSR=COS(1.1*RMX) TNG=SIN(RO(I))/CSR X(NG)=SIN(FI(I))*TNG Y(NG)=COS(FI(I))*TNG IF(KST.EQ.0)GO TO 3 IF(ABS(RO(I)-P2).LT.1E-4)X(NG)=X(NG)+X01 IF(ABS(RO(I)-P2).LT.1E-4)Y(NG)=Y(NG)+Y01 3 CONTINUE XMN=X(1) XMX=X(1) DO 4 I=1,NG IF(X(I).LT.XMN)XMN=X(I) IF(X(I).GT.XMX)XMX=X(I) 4 CONTINUE YMN=Y(1) YMX=Y(1) DO 5 I=1,NG IF(Y(I).LT.YMN)YMN=Y(I) IF(Y(I).GT.YMX)YMX=Y(I) 5 CONTINUE IF(XMN.LT.0)GO TO 6 X0=0. GO TO 8 6 DO 7 I=1,NG 7 X(I)=X(I)-XMN X0=-XMN 8 IF(YMN.LT.0)GO TO 9 Y0=0. GO TO 11 9 DO 10 I=1,NG 10 Y(I)=Y(I)-YMN Y0=-YMN 11 DX=XMX-XMN DY=YMX-YMN IF(DX.GE.DY)FT=200./DX IF(DY.GT.DX)FT=245./DY IF(FT*DX.GT.200.)FT=(2.-FT*DX/200.)*FT IF(KST.EQ.0)GO TO 36 IHM=ID(1,1) IKM=ID(2,1) DO 13 I=1,NG IF(ABS(ID(1,I)).GT.IHM)IHM=ID(1,I) IF(ABS(ID(2,I)).GT.IKM)IKM=ID(2,I) 13 CONTINUE IHM=IHM+1 IKM=IKM+1 IHP=2*IHM+1 IKP=2*IKM+1 IPSN=(PSN*FT/PV+.5) IPCN=(PCN*FT/PO+.5) IQ01=(Q01*FT/PO+.5) LX01=((X01+X0)*FT/PV+.5)+1 LY01=((Y01+Y0)*FT/PO+.5)+1 36 CONTINUE LX0=(X0*FT/PV+0.5)+1 LY0=(Y0*FT/PO+0.5)+1 DO 15 I=1,NG IG(1,I)=(X(I)*FT/PV+0.5)+1 15 IG(2,I)=(Y(I)*FT/PO+0.5)+1 NGM=NG-1 DO 18 I=1,NGM IP1=I+1 DO 18 J=IP1,NG IF(IG(1,I).LE.IG(1,J))GO TO 18 IR=IG(1,I) IG(1,I)=IG(1,J) IG(1,J)=IR IR=IG(2,I) IG(2,I)=IG(2,J) IG(2,J)=IR IR=IK(I) IK(I)=IK(J) IK(J)=IR IF(KST.EQ.0)GO TO 17 DO 16 K=1,3 IR=ID(K,I) ID(K,I)=ID(K,J) 16 ID(K,J)=IR GO TO 18 17 IR=NO(I) NO(I)=NO(J) NO(J)=IR 18 CONTINUE DO 21 I=1,NGM IP1=I+1 DO 21 J=IP1,NG IF(IG(1,I).NE.IG(1,J))GO TO 21 IF(IG(2,I).LE.IG(2,J))GO TO 21 IR=IG(2,I) IG(2,I)=IG(2,J) IG(2,J)=IR IR=IK(I) IK(I)=IK(J) IK(J)=IR IF(KST.EQ.0)GO TO 20 DO 19 K=1,3 IR=ID(K,I) ID(K,I)=ID(K,J) 19 ID(K,J)=IR GO TO 21 20 IR=NO(I) NO(I)=NO(J) NO(J)=IR 21 CONTINUE IC=(245./PO+0.5)+1 IL=(210./PV+0.5)+1 DO 34 I=1,IL DO 22 J=1,110 22 LN(J)=1H IF(KST.EQ.0)GO TO 37 IF(LX01.EQ.I)LN(LY01)=1H+ 37 CONTINUE IF(LX0.EQ.I)LN(LY0)=1HO IF(IL.EQ.I)LN(LY0)=1HX IF(IL.EQ.I)LN(LY0+1)=1HR IF(LX0.EQ.I)LN(IC)=1HY IF(LX0.EQ.I)LN(IC+1)=1HR DO 23 K=1,NG IF(IG(1,K).NE.I)GO TO 23 IF(IK(K).EQ.0)LN(IG(2,K))=1H* IF(IK(K).EQ.1)LN(IG(2,K))=1H. 23 CONTINUE DO 33 K=1,NG IF(IG(1,K).NE.I)GO TO 33 IF(KST.EQ.0)GO TO 28 DO 24 J=1,6 IF(LN(IG(2,K)+J).NE.' ')GO TO 33 24 CONTINUE DO 27 J=1,3 IF(ID(J,K).LT.0)LN(IG(2,K)+2*J-1)=1H- DO 25 L=1,10 IF(ABS(ID(J,K)).EQ.(L-1))GO TO 26 25 CONTINUE 26 LN(IG(2,K)+2*J)=LS(L) 27 CONTINUE GO TO 33 28 IF(LN(IG(2,K)+1).NE.' '.OR.LN(IG(2,K)+2).NE.' ')GO TO 33 IDC=NO(K)/10 IIN=NO(K)-IDC*10 DO 29 L=1,10 IF(IDC.EQ.(L-1))GO TO 30 29 CONTINUE 30 LN(IG(2,K)+1)=LS(L) DO 31 L=1,10 IF(IIN.EQ.(L-1))GO TO 32 31 CONTINUE 32 LN(IG(2,K)+2)=LS(L) 33 CONTINUE IF(KST.EQ.0)GO TO 38 DO 35 K=1,IHP IX1=LX01+(K-IHM-1)*IPSN IF(IX1.NE.I)GO TO 35 DO 14 L=1,IKP IY1=LY01+(K-IHM-1)*IPCN+(L-IKM-1)*IQ01 IF(IY1.LT.1)GO TO 14 IF(IY1.GT.110)GO TO 35 IF(LN(IY1).EQ.' ')LN(IY1)=1H+ 14 CONTINUE 35 CONTINUE 38 CONTINUE WRITE(IO,100)LN 100 FORMAT(1X,110A1) 34 CONTINUE RETURN END