      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
