C
C     This code is distributed under the terms and conditions of the
C     CCP4 licence agreement as `Part ii)' software.  See the conditions
C     in the CCP4 manual for a copyright statement.
C
C************************************************************************
C
C      POLARRFN
C
C
C  Fast polar rotation function program, by W. KABSCH
C
C  15/9/93 PRE: output weighted mean & RMS to map file
C
C   1/1/90 PRE replaced s/r DCOSFD with more robust alternative
C
C  26/6/89 PRE test maximum & minimum amplitude after applying B-factor
C          Apply temperature factor correctly ie multiply F by
C          exp ( - B sintheta/lambda)
C
C   2/6/89 PRE allow peak threshold to be given as multiple of RMS
C        Interpolate peaks
C
C  Fancy input and output by Phil Evans, Feb 1986
C
C  25/5/89 PRE edited variable names ALFA, BETA, GAMA to PHI, OMEG, KAPA
C        Changed plotting to more portable interface (see routines PLT...
C          at end)
C
C  3/3/87  PRE Correction in s/r MAPRED, allowed title change in READ
C        Added JOIN option
C
C  30/7/86 PRE increased map array size, and check
C
C VAX  compile /NOPT (Fortran v4.3 anyway)
C      link  with F77LCFLIB, MAPLIB, PLOT84, CCPLIB
C
C************************************************************************
      PROGRAM POLARRFN
C     ================
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*80 TITLE,FNAME
      COMMON /TEXTC/ TITLE,FNAME
C
      PARAMETER (MAXNPK=100)
      COMMON /PEAKS/ NPEAK,PEAKS(MAXNPK),PKPOS(3,MAXNPK),LPEAK,
     .  LPKORD(MAXNPK),LPKOUT
C
C--     Open main read and write files.
C
      CALL MTZINI
      CALL CCPFYP
      IFAIL=0
C
      CALL CCPDPN(5,'DATA','READONLY','F',1,IFAIL)
      CALL CCPDPN(6,'PRINTER','PRINTER','F',1,IFAIL)
      CALL CCPRCS(6,'POLARRFN', '$Date: 2002/08/06 12:29:31 $')
C
C  Read control data
C
      CALL DATAIN
C
      IF(LPEAK.GT.0) THEN
C Option to print peaks only
            CALL PPEAKS
            CALL CCPERR(0,' Normal Termination from ppeaks')
      ENDIF
C
C  Option to read map calculated on a previous run
C
      IF(LREAD.GT.0) THEN
            CALL MAPRED
            CALL CCPERR(0,' Normal Termination from mapred')
      ENDIF
C
C Calculate ALMN coefficients
C
      IF(LSAVE.GE.0) THEN
            CALL CALMN
      ENDIF
C
C Sum rotation function
      CALL FRFSUM
      CALL CCPERR(0,' Normal termination from frfsum')
C
      END
C
C
C
      SUBROUTINE DATAIN
C     =================
C
C Read in all data
      INTEGER MCOLS
      PARAMETER (MCOLS=500)
      INTEGER NPARM
      PARAMETER (NPARM=200)
      INTEGER MAXNPK
      PARAMETER (MAXNPK=100)
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*200 LABIN_SAVE(2)
      CHARACTER*200 ERRLIN
      CHARACTER*80 TITLE,FNAME
      CHARACTER*80 FILIN(2)
      COMMON /TEXTC/ TITLE,FNAME
C
      COMMON /PEAKS/ NPEAK,PEAKS(MAXNPK),PKPOS(3,MAXNPK),LPEAK,
     .  LPKORD(MAXNPK),LPKOUT
C
      COMMON /PKFILE/ PKFNAM
      CHARACTER*80 PKFNAM
C
      LOGICAL NOTCAL
      CHARACTER*10 NAMSGP,NAMSGS(2),PGNAME,VERSNX
      REAL COLRNG(2,MCOLS)
C
C Maximum section size
      PARAMETER (NSMAP=33759)
C
C JOIN stuff
      PARAMETER (MXJOIN=30)
      COMMON /VJOIN/ NJOIN,PEAKJ(3,2,MXJOIN),DASH(2,MXJOIN)
C
      INTEGER LOOKUP(MCOLS)
      INTEGER LSPGRP,NSM,IPRINT
      REAL RR(3,3,6),RO(4,4),RF(4,4)
      CHARACTER CTPRGI(MCOLS)*1,LSPRGI(MCOLS,2)*30
C
      SAVE /PKFILE/, /VJOIN/
C
C     .. External Functions ..
      INTEGER LENSTR
      EXTERNAL LENSTR
C     ..
C     .. External Subroutines ..
      EXTERNAL CCPUPC,GTPINT,GTPREA,
     +         PARSER,RBRINV,LKYIN,
     +         LRASSN,LRINFO,LROPEN,LRSYMM
      DIMENSION CELL(6,2),CELCHK(6)
C     ..
C
C Things for PARSER ----
      CHARACTER KEY*4,CVALUE(NPARM)*4,LINE*600,LTYPE*1
      INTEGER IBEG(NPARM),IEND(NPARM),ITYP(NPARM),IDEC(NPARM),NTOK
      REAL FVALUE(NPARM)
      LOGICAL LEND
C
      PARAMETER (NKEY=20)
      CHARACTER KEYS(NKEY)*4
      CHARACTER*15  CODES(5)
      DATA KEYS/'TITL','SELF','CROS','RESO','CRYS','CELL','LABI',
     .   'LIMI','PRIN','NOPR','MAP ','PLOT','FIND','READ','SAVE',
     .   'SUM ','PEAK','JOIN','END',' '/
C
      DATA CODES/'a,c*xa,c*','b,a*xb,a*','c,b*xc,b*',
     . 'a+b,c*x(a+b),c*','a*,cxa*,c'/
C
C---- NLPRGI  =  number of input labels
C
      DATA NLPRGI,(LSPRGI(JJ,1),JJ=1,MCOLS)
     +               /5,'H','K','L','F','SIGF',495*' '/
      DATA (LSPRGI(JJ,2),JJ=1,MCOLS)
     +               /'H','K','L','F','SIGF',495*' '/
C
C
C---- This code signs which input columns are essential (LOOKUP)
C
      DATA CTPRGI/'H','H','H','F','Q',495*' '/
      DATA LOOKUP/-1,-1,-1,-1,-1,495*0/
      DATA FILIN/'HKLIN','HKLIN2'/
C
C Set defaults
C
      MAXFILES = 2
      NFILES_IN = 0
      NFILE = 0
      TITLE=' '
      FNAME='COEFFS'
      NCRYST=1
      IPHI1=0
      IPHI2=180
      NPHI=5
      IOMEG1=0
      IOMEG2=180
      NOMEG=5
      IKAPA1=0
      IKAPA2=180
      NKAPA=5
      KRKAP1=-1000
      KRKAP2=-1000
      ARAD=20.
      EPS=0.0
      RESMIN=10000.
      RESMAX=2.
      CELL(1,1)=0.0
      CELL(1,2)=0.0
      DO 2,I=1,2
      NCODES(I)=1
      LSPGPS(I)=1
      TFACS(I)=0.
      FRMINS(I)=0.
      FRMAXS(I)=1.0E20
C
CEJD      CALL ZEROI(SYM(1,1,1,I),16)
      DO 33,J=1,4
      DO 3,K=1,4
3     SYM(K,J,1,I)=0.0
33     SYM(J,J,1,I)=1.0
      NSYM(I)=0
2     CONTINUE
      LPRINT=-100
      LMAP=0
      LPLOT=0
      NCONTR=0
      NPEAK=0
      NJOIN=0
      PEAK=10.
      MAXPEK=10
      LPKRMS=0
      LREAD=0
      LSAVE=0
      LPEAK=0
      NOTCAL = .TRUE.
C
      NERR=0
      IXTL=1
C
C
1     NTOK=NPARM
      LINE=' '
      CALL PARSER(
     .  KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,LEND,.TRUE.)
C
      IF(LEND) GO TO 1000
C Ignore blank line
      IF(NTOK.EQ.0) GO TO 1
C
C Make uppercase
      CALL CCPUPC(KEY)
      DO 6,K=1,NKEY
        IF(KEY.EQ.KEYS(K)) GO TO 9
6     CONTINUE
C
      WRITE(6,7)
7     FORMAT(/' Unrecognized keyword'/)
      NERR=NERR+1
      GO TO 1
C      DATA KEYS/'TITL','SELF','CROS','RESO','CRYS','CELL','LABI',
C     .   'LIMI','PRIN','NOPR','MAP ','PLOT','FIND','READ','SAVE',
C     .   'SUM ','PEAK','JOIN','END',' '/
C
9     GO TO
     . (           10,    20,    30,    40,    50,    55,    60,
     .      70,    80,    90,   100,   110,   120,   130,   140,
     .     150,   160,   170,  1000,  1000),K
C
C            ********************************************************
C TITLe  read title
C
10    TITLE=LINE(IBEG(2):)
      GO TO 1
C
C            ********************************************************
C SELF, set flag for self rotation, read integration radius ARAD, smoothing
C       radius EPS
C
20    NCRYST=1
C
21    CALL GTPREA(2,ARAD,NTOK,ITYP,FVALUE)
      CALL GTPREA(3,EPS,NTOK,ITYP,FVALUE)
      GO TO 1
C
C            ********************************************************
C CROSS, set flag for self rotation, read integration radius ARAD, smoothing
C       radius EPS
30    NCRYST=2
      GO TO 21
C
C            ********************************************************
C RESOlution    read resolution limits in A, if only one treat as high
C              resolution limit
40    CONTINUE
C     read resolution limits in A, if only one treat as high
C     resolution limit
        ITOK = 2
C            ********************************************************
        CALL RDRESO(ITOK,ITYP,FVALUE,NTOK,RESMIN,RESMAX,SMIN,SMAX)
C            ********************************************************
        GO TO 1
C
C
C
C
C---- CRYSTAL   - subsiduary keywords; FILE  SYMM  ORTH  BFAC FLIMITS
C
C---- Check command line for
C---- subsidiary keyword FILE   or NUMBER - crystal number
C---- subsidiary keyword ORTH  - orthogonalisation code NCODE,
C---- subsidiary keyword BFAC  - temperature factor B
C                            Fused = Finput * exp(-Bssq/lsq)
C---- subsidiary keyword FLIMITS  - minimum F and maximum F(order unimportant)
C---- subsidiary keyword CELL  - a b c ( alpha beta  gamma - default = 90.0)
C                      rarely needed - just for peak analysis
C
C
50    CONTINUE
C
      NOTCAL = .FALSE.
C Set defaults
      F1 = 0.0
      F2 = 1.0E20
      NCODE = 1
      BB = 0.0
C
      ITOK = 1
51          ITOK = ITOK + 1
            IF (ITOK.GT.NTOK) THEN
              NCODES(IXTL) = NCODE
              TFACS(IXTL) = BB
              FRMINS(IXTL) = F1
              FRMAXS(IXTL) = F2
              GO TO 1
            END IF
              KEY = LINE(IBEG(ITOK) :IEND(ITOK))
              CALL CCPUPC(KEY)
              IF (KEY.EQ.'FILE'.OR.KEY.EQ.'NUMB') THEN
                ITOK = ITOK + 1
C            ********************************
                CALL GTPINT(ITOK,IXTL,NTOK,ITYP,FVALUE)
C            ********************************
C
C---- Check if allowed range of files ... 1 to 2 (is this enough ???)
C
                IF (IXTL.LE.0 .OR. IXTL.GT.MAXFILES) THEN
                  WRITE (6,FMT='(a,/,a,i4,a,i2,a,/,a)')
     +      ' Key_Word Error, For CRYSTAL FILE_NUMBER line:-',
     +      '  File_Number Given = ',IXTL,
     +      ' Is invalid only files 1 to ',MAXFILES,
     +      ' Allowed, Ignoring this line -',LINE(1:LENSTR(LINE))
                  NERR = NERR+1
                END IF
                GO TO 51
              END IF
C
C
              IF (KEY.EQ.'ORTH') THEN
                ITOK = ITOK + 1
                CALL GTPINT(ITOK,NCODE,NTOK,ITYP,FVALUE)
                GO TO 51
              END IF
C
              IF (KEY.EQ.'BFAC') THEN
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,BB,NTOK,ITYP,FVALUE)
                GO TO 51
              END IF
C
              IF (KEY.EQ.'FLIM') THEN
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,F1,NTOK,ITYP,FVALUE)
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,F2,NTOK,ITYP,FVALUE)
                IF(F1.LT.F2) GO TO 51
                F3 = F1
                F1 = F2
                F2 = F3
                GO TO 51
              END IF
C
              IF (KEY.EQ.'SYMM') THEN
C SYMMetry    read crystal number or  space group
                ITOK = ITOK + 1
C
C          ************************************************************
              CALL RDSYMM(ITOK,LINE,IBEG,IEND,ITYP,FVALUE,ITOK,NAMSGP,
     +            LSPGRP,PGNAME,NSYM(IXTL),NSYMP,SYM(1,1,1,IXTL))
C          ************************************************************
C
C   Only need primitive symmetry operators
                NSYM(IXTL) = NSYMP
                LSPGPS(IXTL) = LSPGRP
                NAMSGS(IXTL) = NAMSGP
                GO TO 51
              END IF
C
C CELL    read cell dimensions - not usually needed.
              IF (KEY.EQ.'CELL') THEN
                DO 52 I=4,6
52              CELL(I,IXTL)=90.
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,CELL(1,IXTL),NTOK,ITYP,FVALUE)
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,CELL(2,IXTL),NTOK,ITYP,FVALUE)
                ITOK = ITOK + 1
                CALL GTPREA(ITOK,CELL(3,IXTL),NTOK,ITYP,FVALUE)
                   IF(ITOK.LT.NTOK)THEN
                   ITOK = ITOK + 1
                   CALL GTPREA(ITOK,CELL(4,IXTL),NTOK,ITYP,FVALUE)
                   ITOK = ITOK + 1
                   CALL GTPREA(ITOK,CELL(5,IXTL),NTOK,ITYP,FVALUE)
                   ITOK = ITOK + 1
                   CALL GTPREA(ITOK,CELL(6,IXTL),NTOK,ITYP,FVALUE)
                   END IF
                GO TO 51
              END IF
C Unrecognized subkeyword
              WRITE (6, FMT='(/A,A,A/)')
     $             ' Key_word Error, unrecognized sub-keyword',
     $             ' for CRYSTAL command:  ',KEY
              NERR = NERR+1
              GO TO 51
C
C
C---- LABIN   (File_Number)  "a=b c=d ...."
C     ==========
C
  60  CONTINUE
              ITOK = 2
              KEY = LINE(IBEG(ITOK) :IEND(ITOK))
C
C          *************
      CALL CCPUPC(KEY)
C          *************
C
      IF (KEY(1:1).EQ.'F') THEN
C
C---- second arguement = 'FILE_NUMBER'
C     test if third argument is integer
C
        IF (ITYP(3).NE.2) THEN
          WRITE (6,FMT='(A,A,/,A,A,A,/,A)')
     +      ' Key_Word Error, Wrong type of input for  LABIN',' line:-',
     +      ' If Sub_KeyWord = FILE_NUMBER  then ',
     +      'third  Arguement MUST be an Integer  Ignoring this',
     +      ' line -',LINE(1:LENSTR(LINE))
          NERR = NERR+1
          GO TO 10
        END IF
C
C---- For an input file being given
C     Get File_Number and check if valid
C
C            ********************************
        CALL GTPINT(3,NFILE,NTOK,ITYP,FVALUE)
C            ********************************
C
        NFILES_IN = NFILES_IN + 1
Cc             max(nfiles_in,nfile)
C
C---- Check if allowed range of files ... 1 to 5 (is this enough ???)
C
        IF (NFILE.LE.0 .OR. NFILE.GT.MAXFILES) THEN
          WRITE (6,FMT='(a,/,a,i4,a,i2,a,/,a)')
     +      ' Key_Word Error, For LABIN FILE_NUMBER line:-',
     +      '  File_Number Given = ',NFILE,
     +      ' Is invalid only files 1 to ',MAXFILES,
     +      ' Allowed, Ignoring this line -',LINE(1:LENSTR(LINE))
          GO TO 10
        END IF
      END IF
C
C---- Now save labels for this input file
C
      LABIN_SAVE(NFILE) = 'LABIN  '//LINE(IBEG(4) :IEND(NTOK))
Cc      GOT_LABIN(NFILE)  = .TRUE.
C
C
      GO TO 1
C
C  CELL   read cell dimensions
C
55    CONTINUE
      GO TO 1
C
C
C
C  LIMIts  read limits and steps on phi,omega,kappa
C
70    CALL GTPINT(2,IPHI1,NTOK,ITYP,FVALUE)
      CALL GTPINT(3,IPHI2,NTOK,ITYP,FVALUE)
      CALL GTPINT(4,NPHI,NTOK,ITYP,FVALUE)
      CALL GTPINT(5,IOMEG1,NTOK,ITYP,FVALUE)
      CALL GTPINT(6,IOMEG2,NTOK,ITYP,FVALUE)
      CALL GTPINT(7,NOMEG,NTOK,ITYP,FVALUE)
      CALL GTPINT(8,IKAPA1,NTOK,ITYP,FVALUE)
      CALL GTPINT(9,IKAPA2,NTOK,ITYP,FVALUE)
      CALL GTPINT(10,NKAPA,NTOK,ITYP,FVALUE)
      GO TO 1
C
C PRINt    read print threshold
C
80    LPRINT=1
      CALL GTPINT(2,LPRINT,NTOK,ITYP,FVALUE)
      GO TO 1
C
C NOPRint  switch off map print
C
90    LPRINT=0
      GO TO 1
C
C MAP   switch to write map to map file
C
100   LMAP=1
      GO TO 1
C
C PLOT  read contour start, interval , NOTItle flag
C
110   LPLOT=1
      C1=10.
      C2=100.
      C=10.
      CALL GTPREA(2,C1,NTOK,ITYP,FVALUE)
      CALL GTPREA(3,C,NTOK,ITYP,FVALUE)
      IF(C1.GT.C2.OR.C.LE.0) THEN
            WRITE(6,9110) C1,C2,C
9110      FORMAT(/' Unreasonable contour limits:',3F10.2/)
            NERR=NERR+1
            GO TO 1
      ENDIF
C read NOTItle flag to suppress title plotting
      IF(NTOK.GE.4) THEN
            CALL CCPUPC(CVALUE(4))
            IF(CVALUE(4).EQ.'NOTI') LPLOT=-1
      ENDIF
      NCONTR=1
      CONTRS(1)=C1
111   CONTRS(NCONTR+1)=CONTRS(NCONTR)+C
      IF(CONTRS(NCONTR+1).GT.C2) GO TO 1
      NCONTR=NCONTR+1
      GO TO 111
C
C  FIND read peak threshold and maximum number of peaks
C    optional keyword "RMS" to indicate that threshold is in
C       multiples of RMS density rather than absolute
C    optional keyword OUTPUT filename to output unique peak list to file
C    NB RMS must precede OUTPUT if both are present
C
120   CALL GTPREA(2,PEAK,NTOK,ITYP,FVALUE)
      CALL GTPINT(3,MAXPEK,NTOK,ITYP,FVALUE)
      MAXPEK=MIN(MAXPEK,MAXNPK)
C  Process rest of line
      LPKRMS=0
      LPKOUT=0
      IF (NTOK .LT.  4) GO TO 1
      DO 121, J = 4, NTOK
         IF (CVALUE(J) .EQ. 'RMS') THEN
            LPKRMS=1
         ELSEIF(CVALUE(J) .EQ. 'OUTP') THEN
            LPKOUT=1
C Store rest of line as filename if given
            IF(NTOK .GT. J) THEN
               PKFNAM=LINE(IBEG(J+1):)
               GO TO 1
            ELSE
               PKFNAM='PEAKS'
            ENDIF
         ENDIF
121   CONTINUE
      GO TO 1
C
C READ   set flag to read map in, instead of calculating it
C        read 1st and last sections to process (kappa values)
C
130   LREAD=1
      CALL GTPINT(2,KRKAP1,NTOK,ITYP,FVALUE)
      CALL GTPINT(3,KRKAP2,NTOK,ITYP,FVALUE)
      IF (KRKAP2.LT.-999) THEN
         KRKAP2 = KRKAP1
      ENDIF
      IF (KRKAP1 .LT. -999 .OR. KRKAP2 .LT. KRKAP1) THEN
         WRITE (6, '(A)') ' *** Illegal read limits ***'
         NERR = NERR+1
      ENDIF
      GO TO 1
C
C SAVE  read file name for saving ALMN output
140   LSAVE=1
C
141   IF(NTOK.GT.1) THEN
            FNAME=LINE(IBEG(2):IEND(2))
C&&&VAX !!!
C Strip off any existing extension
            I=INDEX(FNAME,']')
            J=INDEX(FNAME(I+1:),'.')
            IF(J.GT.0) FNAME=FNAME(1:J-1)
      ENDIF
      GO TO 1
C
C SUM  read file name for reading previous ALMN results
C
150   LSAVE=-1
      GO TO 141
C
C PEAK read list of peaks omega, phi, kappa from following cards, terminated
C      by blank line or 'END'
C
160   LPEAK=1
      NPEAK=0
      IF (NTOK .GT. 1) THEN
C  Test second token to see if it the keyword MATRices
         CALL CCPUPC(CVALUE(2))
         IF (CVALUE(2) .EQ. 'MATR') THEN
            LPEAK=2
            WRITE(6,*) ' Matrices will be output for each peak'
         ENDIF
      ENDIF
161   NTOK=NPARM
      LINE=' '
      CALL PARSER(
     .  KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,LEND,.TRUE.)
C
      IF(LEND) GO TO 1000
C Make uppercase
      CALL CCPUPC(KEY)
      IF(NTOK.LE.0.OR.KEY.EQ.'END') GO TO 1
      IF(NPEAK.GE.MAXNPK) THEN
            WRITE(6,162) MAXNPK
162       FORMAT(/' Too many peaks given, maximum =',I4/)
            NERR=NERR+1
            GO TO 1
      ENDIF
      IF(NTOK.LT.3) THEN
            WRITE(6,163)
163       FORMAT(/' !!! Three angles must be given !!!')
            GO TO 161
      ENDIF
      NPEAK=NPEAK+1
      CALL GTPREA(1,PKPOS(1,NPEAK),NTOK,ITYP,FVALUE)
      CALL GTPREA(2,PKPOS(2,NPEAK),NTOK,ITYP,FVALUE)
      CALL GTPREA(3,PKPOS(3,NPEAK),NTOK,ITYP,FVALUE)
      LPKORD(NPEAK)=NPEAK
      GO TO 161
C
C
C JOIN  read peak positions to join up
C  kappa, omega1, phi1, omega2, phi2, [DASHED  repeat  dash]
C                                         (in fractions of radius)
C
170   IF(NJOIN.GE.MXJOIN) THEN
            WRITE(6,1170) MXJOIN
1170      FORMAT(/' !!!! No more than ',I4,' joins can be stored'//)
            GO TO 1
      ENDIF
      IF(NTOK.LT.6) THEN
            WRITE(6,1171)
1171      FORMAT(' Too few numbers for JOIN !!!')
            NERR=NERR+1
            GO TO 1
      ENDIF
C
      NJOIN=NJOIN+1
C read kappa
      CALL GTPREA(2,PEAKJ(3,1,NJOIN),NTOK,ITYP,FVALUE)
C Read 2 two sets of omega, phi
      K=2
      DO 171,J=1,2
      DO 171,I=1,2
      K=K+1
      CALL GTPREA(K,PEAKJ(I,J,NJOIN),NTOK,ITYP,FVALUE)
171   CONTINUE
C
C Check field 7 for keyword DASHed
      IF(NTOK.GE.7.AND.ITYP(7).EQ.1) THEN
            CALL CCPUPC(CVALUE(7))
            IF(CVALUE(7).EQ.'DASH') THEN
C    Set default repat and dash length
                  DASH(1,NJOIN)=0.05
                  DASH(2,NJOIN)=0.025
                  CALL GTPREA(8,DASH(1,NJOIN),NTOK,ITYP,FVALUE)
                  CALL GTPREA(9,DASH(2,NJOIN),NTOK,ITYP,FVALUE)
            ELSE
C     Solid
                  DASH(1,NJOIN)=0.0
            ENDIF
      ENDIF
      GO TO 1
C
C---------------------------------------------------------------------
C
C End of input
C
C Check for errors
1000  IF(NERR.GT.0) THEN
            WRITE(6,9901) NERR
9901      FORMAT(/1X,I6,' errors in input, see above !!!!'/)
            CALL CCPERR(1,' error in input')
      ENDIF
C   Open MTZ files
C---- Open input file for Read; print header...
C
      IF(NFILES_IN .GE.0) THEN
             WRITE(6,*) ' Using ',NFILES_IN,' files'
            DO 1005 IFIL=1,NFILES_IN
            IERR = -1
            IPRINT = 2
            LOOKUP(1) = -1
            LOOKUP(2) = -1
            LOOKUP(3) = -1
            LOOKUP(4) = -1
            LOOKUP(5) = -1
C
C           ==================================
            CALL LROPEN(IFIL,FILIN(IFIL),IPRINT,IERR)
C           ==================================
C
            ERRLIN = 'ERROR:HEADER NOT THERE, ILLEGAL FILE'
C
C                          =========================================
            IF (IERR.EQ.3) CALL LERROR(2,1,ERRLIN(1:LENSTR(ERRLIN)))
C                         =========================================
C
C
C
Cc         line= labin_save(nfile)
         LINE= LABIN_SAVE(IFIL)
C
C             ======================================================
         CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,FVALUE,CVALUE,IDEC,NTOK,
     +               LEND,.TRUE.)
         CALL LKYIN(IFIL,LSPRGI(1,IFIL),NLPRGI,NTOK,LINE,IBEG,IEND)
         CALL LRASSN(IFIL,LSPRGI(1,IFIL),NLPRGI,LOOKUP,CTPRGI)
C             ======================================================
C
C           ======================
            CALL LRCELL(IFIL,CELCHK)
C           ======================
C      check cell
            IF(CELL(1,IFIL).NE.0.0) THEN
              ACHK = 0.0
              DO 351 II = 1,6
              ACHK = ACHK +
     +               ABS(CELCHK(II) - CELL(II,IFIL)) /
     +                  (CELCHK(II) + CELL(II,IFIL))
              IF (ACHK.GT.0.02)CALL CCPERR(1,
     +        ' Large Difference in CELL DIMENSIONS')
              IF (ACHK.GT.0.01) WRITE(6,*)
     +         ' Small Difference in CELL DIMENSION'
351           CONTINUE
            ELSE
              DO 352 II = 1,6
              CELL(II,IFIL) =  CELCHK(II)
352           CONTINUE
            END IF
C
C---- Find out how many columns in input file?
C
C           ========================================
            CALL LRINFO(IFIL,VERSNX,NCOL,NREFTOT,COLRNG)
C           ========================================
C
C    Read symmetry and check
        CALL LRSYMI(IFIL,NSYMP,LTYPE,LSPGRP,NAMSGP,PGNAME)
        CALL LRSYMM(IFIL,NSM,SYM(1,1,1,IFIL) )
C
        IF(NSYM(IFIL).EQ.0 )THEN
C   Only need primitive symmetry operators
        NSYM(IFIL) = NSYMP
        LSPGPS(IFIL) = LSPGRP
        NAMSGS(IFIL) = NAMSGP
        END IF
C
        IF(NSYM(IFIL).GT.0 .AND. NSYM(IFIL).NE.NSYMP)THEN
C   Stop
        CALL CCPERR(1,' Inconsistent symmetry information')
        END IF
C
1005  CONTINUE
      END IF
C
       DO 1006 IFIL = 1,NCRYST
C  Call pgdefn to make sure you have a primitive set of symops.
C
C---- Check NCODE defined
C
         IF (NCODES(IFIL) .LT. 1 .OR. NCODES(IFIL) .GT. 6) THEN
           WRITE (6, 6101) J
6101       FORMAT(' !!!! Illegal or undefined NCODE for crystal ',I4/)
           CALL CCPERR(1,' error in ncode')
         ENDIF
C
C
C  Set up fractional to orthogonal matrix RR with appropriate NCODE.
      CALL RBFRO1(CELL(1,IFIL),VOLL,RR)
C
      DO 1007 I=1,3
      DO 1007 J=1,3
      RO(I,4)=0.0
      RO(4,I)=0.0
      RO(I,J)=RR(I,J,NCODES(IFIL))
      ORTHOG(I,J,IFIL)=RR(I,J,NCODES(IFIL))
1007     CONTINUE
      RO(4,4)=1.0
C
C  and its inverse
         CALL RBRINV(RO,RF)
      DO 1008 I=1,3
      DO 1008 J=1,3
      FRAC(I,J,IFIL)=RF(I,J)
1008     CONTINUE
1006  CONTINUE
C
C Check that crystal cards were given and LCF files opened
C unless READ or SUM options used
      IF(NOTCAL.AND.LREAD.EQ.0.AND.
     .              LSAVE.GE.0.AND.LPEAK.EQ.0) THEN
            WRITE(6,9902)
9902      FORMAT(/' No CRYSTAL card given !!!'/)
            CALL CCPERR(1,' error in crystal')
      ENDIF
C
C If 1 crystal, copy symmetry etc to crystal 2
      IF(NCRYST.EQ.1) THEN
            DO 9903 IS =1,NSYM(1)
            DO 9903 I =1,4
            DO 9903 J =1,4
9903        SYM(I,J,IS,2) = SYM(I,J,IS,1)
            NSYM(2)=NSYM(1)
            NCODES(2)=NCODES(1)
            LSPGPS(2)=LSPGPS(1)
CEJD            CALL MOVEI(ORTHOG(1,1,2),ORTHOG(1,1,1),9)
CEJD            CALL MOVEI(FRAC(1,1,2),FRAC(1,1,1),9)
            DO 9900 I =1,3
            DO 9900 J =1,3
            ORTHOG(J,I,2) = ORTHOG(J,I,1)
9900            FRAC(J,I,2) = FRAC(J,I,1)
      ENDIF
C
C  Default title
C
      IF(LENSTR(TITLE).EQ.0) THEN
            IF(NCRYST.EQ.1) THEN
                  TITLE='Self-rotation function'
            ELSE
                  TITLE='Cross-rotation function'
            ENDIF
      ENDIF
C
C Print control data
C
C  PEAK option
C
      IF(LPEAK.GT.0) THEN
            WRITE(6,9101) NPEAK
9101        FORMAT('1    Rotation function peak calculation'///
     .       ' All symmetry related peaks and matrices will be',
     .       ' calculated for ',I4,' input peaks'//)
            DO 1050,I=1,NCRYST
            J = LENSTR(NAMSGS(I))
            WRITE(6,9102) I,NCODES(I),CODES(NCODES(I)),LSPGPS(I),
     .                NAMSGS(I)(1:J),NSYM(I)
9102      FORMAT(/' Crystal',I3,': orthogonalisation code (NCODE) =',
     .         I2,'   x y z axes along ',A//' Space group ',I4,
     .         ' (',A,'),',I4,' symmetry operations'//)
            WRITE(6,9103) (CELL(J,I),J=1,6)
9103        FORMAT(' Cell dimensions: ',6F8.2/)
1050        CONTINUE
            RETURN
      ENDIF
C
      IF(LREAD.EQ.0) THEN
C Set print default if not set explicitly
            IF(LPRINT.LT.-99) LPRINT=1
            IF(NCRYST.EQ.1) THEN
                  WRITE(6,9001) 'Self',TITLE
9001      FORMAT('1',20X,
     .     'Polar angle rotation function:  ',A,'-rotation function'/
     . 21X,'======================================================='//
     .  ' Title:  ',A//)
            ELSE
                  WRITE(6,9001) 'Cross',TITLE
            ENDIF
C
            IF(ARAD.GT.RESMAX*5.83) THEN
C  Reset maximum resolution to integration radius/5.83
                  RESMAX=ARAD/5.83
                  WRITE(6,9002) RESMAX,ARAD
9002              FORMAT(/
     .             '  !!! Warning !!! maximum resolution reset to',
     .             F10.3,' = integration radius ',F10.3,'/ 5.83'/)
            ENDIF
C
            WRITE(6,9003) ARAD,RESMIN,RESMAX
9003        FORMAT(/' Integration radius = ',F9.2,'A',20X,
     .         'Resolution limits ',2F9.2)
C
            IF(EPS.GT.0.) WRITE(6,9004) EPS
9004      FORMAT(/' Radius for averaging (smoothing) ',F10.3)
C
            DO 1010,I=1,NCRYST
            J=LENSTR(NAMSGS(I))
            WRITE(6,9010) I,NCODES(I),CODES(NCODES(I)),LSPGPS(I),
     .         NAMSGS(I)(1:J),NSYM(I),TFACS(I),FRMINS(I),FRMAXS(I)
9010      FORMAT(/' Crystal',I3,': orthogonalisation code (NCODE) =',
     .         I2,'   x y z axes along ',A//' Space group ',I4,
     .         ' (',A,'),',I4,
     .         ' symmetry operations'//' Temperature factor = ',F10.3/
     .         5X,'Minimum F =',F10.0,'      Maximum F = ',G14.4/)
            WRITE(6,9015) (CELL(J,I),J=1,6)
9015        FORMAT(' Cell dimensions: ',6F8.2/)
1010        CONTINUE
C
            IF(LSAVE.GT.0) THEN
                  L=LENSTR(FNAME)
                  WRITE(6,9011) FNAME(1:L)
9011      FORMAT(/
     .       ' Output from ALMN stage will be saved in file(s) : ',A)
            ELSEIF(LSAVE.LT.0) THEN
                  L=LENSTR(FNAME)
                  WRITE(6,9012) FNAME(1:L)
9012      FORMAT(/
     .       ' Output from ALMN stage will be read from file(s) : ',A)
            ENDIF
C
C Check for maximum section size
            NSECSZ=((IPHI2-IPHI1)/NPHI+1)*((IOMEG2-IOMEG1)/NOMEG+1)
            IF(NSECSZ.GT.NSMAP) THEN
                  WRITE(6,9050) NSECSZ,NSMAP
9050      FORMAT(//' Section size too large, =',I8,
     .             '    maximum=',I8///)
                CALL CCPERR(1,' error in section')
            ENDIF
      ELSE
C  READ option
            WRITE(6,9020) TITLE
9020      FORMAT(/
     .          ' Previously calculated map will be read from MAPIN'//
     .          ' Title: ',A/)
C no save, no map output
            LSAVE=0
            LMAP=0
      ENDIF
C
C
      END
C************************************************************************
C******************* ND500-ALMN ********************** W.KABSCH 10-1982
C************************************************************************
C************************* FILES ****************************************
C
      SUBROUTINE CALMN
C     ===============
C
C
C OUTPUT FILE
C -----------
C THE OUTPUT DATA FILE CONTAINS BINARY DATA AS FOLLOWS:
C
C         L,M,N,ALMNR,ALMNI
C
C L,M,N     ARE THE INDICES OF THE SPHERICAL HARMONIC COEFFICIENT.
C            THE LAST RECORD HAS L=10000 TO INDICATE THE END OF VALID
C            DATA ON FILE.
C ALMNR,    ARE REAL AND IMAGINARY PART OF SPHERICAL HARMONIC
C ALMNI     EXPANSION COEFFICIENT
C
C ----------
C
C PARAMETERS
C    ARAD                 RADIUS (ANGSTROM) OF CUT-OFF SPHERE IN
C                         PATTERSON FUNCTION. THIS PROGRAM ALLOWS A
C                         LARGEST VALUE FOR L OF 30. THIS LEADS TO THE
C                         CONDITION ARAD<RESMAX*5.83  .
C    RESMIN               LOW  RESOLUTION LIMIT (ANGSTROM),
C    RESMAX               HIGH RESOLUTION LIMIT FOR ACCEPTING DATA.
C    EPS                  RADIUS (ANGSTROM) OF SPHERE WITHIN WHICH THE
C                         PATTERSON IS AVERAGED. EPS SHOULD BE LESS
C                         THAN 0.7*RESMAX.
C    FRMIN                REFLECTIONS WITH INTENSITY EQUAL OR LESS THAN
C                         FRMIN ARE NOT INCLUDED IN THE CALCULATION.
C    NEQ                  NUMBER OF EQUIVALENT POSITIONS IN SPACE GROUP
C                         OF PATTERSON. THIS INCLUDES THE UNIT MATRIX,
C                         BUT NOT POSITIONS RELATED BY A CENTRE OF
C                         SYMMETRY AT THE ORIGIN.
C    LIST                 0 - NO LISTING
C                         1 - LIST SPHERICAL HARMONIC COEFFICIENTS
C
C 3.)UNIT CELL ORIENTATION
C    BASIS(1),(2),(3)  X,Y,Z COMPONENTS (ANGSTROM) OF A-AXIS,
C    BASIS(4),(5),(6)  X,Y,Z COMPONENTS (ANGSTROM) OF B-AXIS,
C    BASIS(7),(8),(9)  X,Y,Z COMPONENTS (ANGSTROM) OF C-AXIS
C                         OF UNIT CELL RELATIVE TO REFERENCE FRAME.
C
C 4.)EQUIVALENT POSITIONS (NEQ OF THEM)
C    MEQ(1),...,(9)    H'=MEQ(1)*H+MEQ(2)*K+MEQ(3)*L,
C                         K'=MEQ(4)*H+MEQ(5)*K+MEQ(6)*L,
C                         L'=MEQ(7)*H+MEQ(8)*K+MEQ(9)*L,
C                         WHERE H',K',L' IS A SYMMETRY RELATED
C                         REFLECTION TO H,K,L.
C
C***********************************************************************
C
C
      INTEGER MCOLS
      PARAMETER (MCOLS=500)
C
      REAL ADATA(MCOLS),SS
      LOGICAL LOGMSS(MCOLS)
      LOGICAL EOF
      DIMENSION RSYM(4,4)
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*80 TITLE,FNAME
      COMMON /TEXTC/ TITLE,FNAME
C
      PARAMETER (NSRANG=10)
      DIMENSION NRRANG(NSRANG)
C
      INTEGER TAPOT,H
      DIMENSION BZ02(15),BZ04(15),BZ06(15),BZ08(15),BZ10(15),BZ12(15),
     1BZ14(15),BZ16(15),BZ18(15),BZ20(15),BZ22(15),BZ24(15),BZ26(15),
     2BZ28(15),BZ30(15),BZ(15,15)
      DIMENSION NMAX(15),IH(48),IK(48),IL(48),MEQ(216)
      COMMON /POLDAT/ FACTOR(15),FAC1(15,30),FAC2(15,30),BASIS(9),
     1   RBASIS(9),SPHBES(30),RLN(15),PLM(31),CMPHI(31),SMPHI(31),
     2   ARRAY(930),ALMN(10000)
      EQUIVALENCE (BZ02(1),BZ(1,1)),(BZ04(1),BZ(1,2)),(BZ06(1),BZ(1,3)),
     1(BZ08(1),BZ(1,4)),(BZ10(1),BZ(1,5)),(BZ12(1),BZ(1,6)),
     2(BZ14(1),BZ(1,7)),(BZ16(1),BZ(1,8)),(BZ18(1),BZ(1,9)),
     3(BZ20(1),BZ(1,10)),(BZ22(1),BZ(1,11)),(BZ24(1),BZ(1,12)),
     4(BZ26(1),BZ(1,13)),(BZ28(1),BZ(1,14)),(BZ30(1),BZ(1,15))
C
      CHARACTER NUMS(2)*1,FILEN*80,STATUS*8
      DATA NUMS/'1','2'/
C
C*****INITIALIZE TABLE OF BESSEL  ZEROS
      DATA BZ02/5.7635,9.0950,12.3229,15.5146,18.6890,21.8539,25.0128,
     @28.1678,31.3201,34.4705,37.6194,40.7671,43.9140,47.0601,50.2057/
      DATA BZ04/8.1826,11.7049,15.0395,18.3012,21.5254,24.7275,27.9156,
     @31.0939,34.2654,37.4317,40.5942,43.7536,46.9106,50.0656,53.2191/
      DATA BZ06/10.5128,14.2074,17.6479,20.9834,24.2627,27.5078,30.7304,
     @33.9371,37.1323,40.3189,43.4987,46.6733,49.8436,53.0105,56.1745/
      DATA BZ08/12.7908,16.6410,20.1824,23.5912,26.9270,30.2172,33.4768,
     @36.7145,39.9361,43.1454,46.3451,49.5371,52.7229,55.9035,59.0799/
      DATA BZ10/15.0335,19.0257,22.6626,26.1427,29.5346,32.8705,36.1681,
     @39.4382,42.6876,45.9212,49.1422,52.3531,55.5559,58.7517,61.9419/
      DATA BZ12/17.2504,21.3738,25.1010,28.6498,32.0966,35.4780,38.8140,
     @42.1169,45.3950,48.6537,51.8970,55.1279,58.3485,61.5606,64.7655/
      DATA BZ14/19.4477,23.6931,27.5057,31.1206,34.6211,38.0472,41.4214,
     @44.7574,48.0644,51.3488,54.6149,57.8663,61.1056,64.3345,67.5549/
      DATA BZ16/21.6292,25.9889,29.8828,33.5612,37.1139,40.5839,43.9958,
     @47.3648,50.7008,54.0110,57.3004,60.5728,63.8310,67.0774,70.3136/
      DATA BZ18/23.7979,28.2653,32.2365,35.9762,39.5795,43.0925,46.5416,
     @49.9432,53.3082,56.6442,59.9570,63.2505,66.5281,69.7921,73.0447/
      DATA BZ20/25.9557,30.5251,34.5704,38.3692,42.0215,45.5765,49.0622,
     @52.4960,55.8897,59.2516,62.5876,65.9025,69.1995,72.4814,75.7505/
      DATA BZ22/28.1043,32.7707,36.8871,40.7428,44.4428,48.0388,51.5603,
     @55.0259,58.4480,61.8356,65.1950,68.5311,71.8416,75.1476,78.4333/
      DATA BZ24/30.2450,35.0040,39.1886,43.0993,46.8457,50.4818,54.0384,
     @57.5353,60.9855,64.3985,67.7810,71.1384,74.4745,77.7925,81.0950/
      DATA BZ26/32.3787,37.2263,41.4766,45.4407,49.2322,52.9073,56.4985,
     @60.0260,63.5040,66.9421,70.3477,73.7261,77.0817,80.4179,83.7372/
      DATA BZ28/34.5063,39.4390,43.7527,47.7685,51.6039,55.3172,58.9421,
     @62.4999,66.0052,69.4682,72.8965,76.2960,79.6710,83.0251,86.3613/
      DATA BZ30/36.6284,41.6430,46.0180,50.0839,53.9621,57.7127,61.3707,
     @64.9582,68.4904,71.9780,75.4289,78.8492,82.2435,85.6156,88.9685/
C*****ASSIGN FORTRAN DEVICE UNIT NUMBERS
      DATA TAPOT,TWOPI/8,6.283185/
C
C
1000  FORMAT(A80)
1010  FORMAT(5F10.2,2I5)
1020  FORMAT(9F8.2)
1030  FORMAT(9I3)
1040  FORMAT(//' LMAX=',I2,10X,
     1 'NUMBER OF SPHERICAL HARMONIC COEFFICIENTS =',I6/1X,
     2 '         ,NMAX(L) ',15I4/)
1050  FORMAT(3I4,2E11.4)
C
      LIST=0
C      LIST=1
C
C
C---------------------------------------------------------------------------
C
C Loop on 1 or two crystals
C
      NCRYST=MIN(MAX(1,NCRYST),2)
      RESINT=1./((RESMAX**2)*FLOAT(NSRANG))
C
      DO 10,IXTL=1,NCRYST
C
C Open output file for ALMN coefficients
      L=LENSTR(FNAME)
      FILEN=FNAME(1:L)//NUMS(IXTL)
C Status normally SCRATCH, but NEW if SAVE option used
      IF(LSAVE.GT.0) THEN
            STATUS='NEW'
      ELSE
            STATUS='SCRATCH'
      ENDIF
      IFAIL=0
      CALL CCPDPN(TAPOT,FILEN,STATUS,'U',1,IFAIL)
C
C Set up orthogonalisation matrix  (basis set)
         DO 133 I=1,3
          BASIS(3*I -2) = ORTHOG(1,I,IXTL)
          BASIS(3*I -1) = ORTHOG(2,I,IXTL)
133       BASIS(3*I   ) = ORTHOG(3,I,IXTL)
CEJD      CALL MOVEI(BASIS(1),ORTHOG(1,1,IXTL),9)
      WRITE(6,1002) IXTL,(BASIS(I),BASIS(I+3),BASIS(I+6),I=1,3)
1002  FORMAT(/' Crystal number : ',I5/
     .  ' Basis set (orthogonalisation matrix) :',3F10.5,2(/39X,3F10.5))
C
C
C Invert and transpose symmetry to get reciprocal space matrices
      DO 14,L=1,NSYM(IXTL)
            CALL RBRINV(SYM(1,1,L,IXTL),RSYM)
            L0=(L-1)*9
            DO 13,I=1,3
            DO 13,J=1,3
                  K=(I-1)*3+J+L0
                  MEQ(K)=NINT(RSYM(J,I))
13          CONTINUE
            WRITE(6,1003) L,(MEQ(I+L0),I=1,9)
1003      FORMAT(/' Symmetry operation ',I3,': ',3I4,2(/25X,3I4)/)
14    CONTINUE
      NEQ=NSYM(IXTL)
      FRMIN=FRMINS(IXTL)
      FRMAX=FRMAXS(IXTL)
      TFAC2=TFACS(IXTL)/4.
C Counts of reflections in resolution ranges
CEJD      CALL ZEROI(NRRANG,NSRANG)
      DO 15 ISRANG = 1,NSRANG
15    NRRANG(ISRANG) = 0
C
C
C*****DETERMINE MAXIMUM VALUE OF L AND TOTAL NUMBER OF SPHERICAL
C*****HARMONIC COEFFICIENTS.
      A2PI=TWOPI*ARAD
      RESMIN=AMIN1(A2PI,RESMIN)
      RMAX=A2PI/RESMAX
      DO 50 LBY2=1,15
      IF (BZ(1,LBY2).GT.RMAX)GO TO 50
      LBY2MX=LBY2
      DO 40 N=1,15
      IF (BZ(N,LBY2).GT.RMAX)GO TO 50
      NMAX(LBY2)=N
40    CONTINUE
50    CONTINUE
      LMAX=LBY2MX+LBY2MX
      LMAX1=LMAX+1
      X=1.
      LGTH=0
      DO 70 LBY2=1,LBY2MX
      L=LBY2+LBY2
      L2=L+L
      X=X*FLOAT((L2-3)*(L2-1))/FLOAT(L2*(L2-2))
      FACTOR(LBY2)=SQRT(X*FLOAT(L2+1))
      LGTH=LGTH+(L+1)*NMAX(LBY2)
      DO 70 M1=1,L
      Y=(L+1-M1)*(L+M1)
      FAC1(LBY2,M1)=-FLOAT(M1+M1)/SQRT(Y)
70    FAC2(LBY2,M1)=-SQRT(FLOAT((L-M1)*(L+M1+1))/Y)
      LGTH2=LGTH+LGTH
      WRITE(6,1040)LMAX,LGTH,(NMAX(LBY2),LBY2=1,LBY2MX)
      IF (LGTH2.GT.10000)
     + CALL CCPERR(1, ' ERROR IN ARRAY SIZE')
C*****DETERMINE COMPONENTS OF RECIPROCAL CELL VECTORS
      RBASIS(1)=BASIS(5)*BASIS(9)-BASIS(8)*BASIS(6)
      RBASIS(2)=BASIS(6)*BASIS(7)-BASIS(9)*BASIS(4)
      RBASIS(3)=BASIS(4)*BASIS(8)-BASIS(7)*BASIS(5)
      RBASIS(4)=BASIS(8)*BASIS(3)-BASIS(2)*BASIS(9)
      RBASIS(5)=BASIS(9)*BASIS(1)-BASIS(3)*BASIS(7)
      RBASIS(6)=BASIS(7)*BASIS(2)-BASIS(1)*BASIS(8)
      RBASIS(7)=BASIS(2)*BASIS(6)-BASIS(5)*BASIS(3)
      RBASIS(8)=BASIS(3)*BASIS(4)-BASIS(6)*BASIS(1)
      RBASIS(9)=BASIS(1)*BASIS(5)-BASIS(4)*BASIS(2)
      X=BASIS(1)*RBASIS(1)+BASIS(4)*RBASIS(4)+BASIS(7)*RBASIS(7)
      DO 80 J1=1,9
80    RBASIS(J1)=RBASIS(J1)/X
      CONST=ARAD*SQRT(A2PI*4.)/X
C*****CLEAR ARRAY ALMN WHICH WILL BE USED FOR STORAGE OF SPHERICAL
C*****HARMONIC COEFFICIENTS ALMN.
      DO 90 J1=1,LGTH2
90    ALMN(J1)=0.
C
      NREFL=0
      NREJMX=0
C*****READ ALL REFLECTION INTENSITIES AND CALCULATE SPHERICAL
C*****HARMONIC COEFFICIENTS ALMN.
C
100   CONTINUE
             CALL LRREFF(IXTL,SS,ADATA,EOF)
             IF(EOF) GO TO 400
             CALL LRREFM(IXTL,LOGMSS)
      IF(LOGMSS(5)) GO TO 100
      IF(ABS(ADATA(5)).LE.0.0001) GO TO 100
      H=NINT(ADATA(1))
      K=NINT(ADATA(2))
      L=NINT(ADATA(3))
C*****CHECK IF REFLECTION IS WITHIN SPECIFIED RESOLUTION RANGE
C      X=H*RBASIS(1)+K*RBASIS(4)+L*RBASIS(7)
C      Y=H*RBASIS(2)+K*RBASIS(5)+L*RBASIS(8)
C      Z=H*RBASIS(3)+K*RBASIS(6)+L*RBASIS(9)
C      R2=X*X+Y*Y+Z*Z
      R2=SS
      R=SQRT(R2)
      IF ((R*RESMAX).GT.1.0)GO TO 100
      IF ((R*RESMIN).LT.1.0)GO TO 100
C***** Apply temperature factor before squaring F
C      FR=IA(5)
      FR=0.0
      IF(.NOT. LOGMSS(4)) FR=ADATA(4)
      IF(TFAC2.NE.0.) FR=FR*EXP(-TFAC2*R2)
C*****CHECK IF AMPLITUDE ABOVE THRESHOLD
      IF (FR.LE.FRMIN)GO TO 100
      IF(FR.GT.FRMAX) THEN
         NREJMX=NREJMX+1
         GO TO 100
      ENDIF
      FR=(FR**2)/10000.
C
      NREFL=NREFL+1
C Count in resolution ranges
      NR=INT(R2/RESINT)+1
      NR=MIN(NR,NSRANG)
      NRRANG(NR)=NRRANG(NR)+1
C***** ATTACH SIGMA-FACTOR TO INTENSITY
      T0=TWOPI*R*EPS
      IF (T0.GT.0.01)FR=FR*3.*(SIN(T0)-T0*COS(T0))/(T0*T0*T0)
C***** CALCULATE BESSELFUNCTIONS BY BACKWARDS RECURSION
      A2PIR=A2PI*R
      J0=A2PIR+12.
      J2=J0
      T1=0.
      T2=1.
      T0=J0+J0+1
210   T3=(J0+J0+1)*T2/A2PIR-T1
      ARRAY(J0)=T3
      T1=T2
      T2=T3
      T0=T0+(J0+J0-1)*T3*T3
      J0=J0-1
      IF (J0.GT.0)GO TO 210
      T1=1./SQRT(T0)
      DO 220 J1=2,LMAX1
      IF (J1.GT.J2) ARRAY(J1)=0.
220   SPHBES(J1-1)=T1*ARRAY(J1)
C*****GENERATE ALL SYMMETRY RELATED REFLECTIONS
      J0=0
      NRFL=1
      IH(1)=H
      IK(1)=K
      IL(1)=L
      DO 250 ISYM=1,NEQ
      IH0=H*MEQ(J0+1)+K*MEQ(J0+2)+L*MEQ(J0+3)
      IK0=H*MEQ(J0+4)+K*MEQ(J0+5)+L*MEQ(J0+6)
      IL0=H*MEQ(J0+7)+K*MEQ(J0+8)+L*MEQ(J0+9)
      DO 240 J2=1,NRFL
      IF (IH0.NE.IH(J2))GO TO 230
      IF (IK0.NE.IK(J2))GO TO 230
      IF (IL0.NE.IL(J2))GO TO 230
      GO TO 250
230   IF (-IH0.NE.IH(J2))GO TO 240
      IF (-IK0.NE.IK(J2))GO TO 240
      IF (-IL0.NE.IL(J2))GO TO 240
      GO TO 250
240   CONTINUE
      NRFL=NRFL+1
      IH(NRFL)=IH0
      IK(NRFL)=IK0
      IL(NRFL)=IL0
250   J0=J0+9
      DO 260 J0=1,930
260   ARRAY(J0)=0.
      DO 330 J0=1,NRFL
      X=IH(J0)*RBASIS(1)+IK(J0)*RBASIS(4)+IL(J0)*RBASIS(7)
      Y=IH(J0)*RBASIS(2)+IK(J0)*RBASIS(5)+IL(J0)*RBASIS(8)
      Z=IH(J0)*RBASIS(3)+IK(J0)*RBASIS(6)+IL(J0)*RBASIS(9)
      SINTH=SQRT(X*X+Y*Y)/R
      IF (SINTH.LT.0.01)GO TO 310
      COTTH=Z/(R*SINTH)
      PHI=ATAN2(Y,X)
      DO 280 M=1,LMAX1
      CMPHI(M)=COS((M-1)*PHI)
280   SMPHI(M)=SIN((M-1)*PHI)
      J1=0
      DO 300 LBY2=1,LBY2MX
      MM=LBY2+LBY2+1
      M=MM-1
      T1=FACTOR(LBY2)*(SINTH**M)
      T2=0.
      PLM(MM)=T1
290   T0=T1*COTTH*FAC1(LBY2,M)+T2*FAC2(LBY2,M)
      PLM(M)=T0
      T2=T1
      T1=T0
      M=M-1
      IF (M.GT.0)GO TO 290
      DO 300 M=1,MM
      J1=J1+1
      ARRAY(J1)=ARRAY(J1)+PLM(M)*CMPHI(M)
300   ARRAY(J1+465)=ARRAY(J1+465)-PLM(M)*SMPHI(M)
      GO TO 330
310   J1=1
      DO 320 LBY2=1,LBY2MX
      MM=LBY2+LBY2+1
      ARRAY(J1)=ARRAY(J1)+SQRT(FLOAT(MM+MM-1))
320   J1=J1+MM
330   CONTINUE
      J0=0
      J1=0
      T1=CONST*FR
      DO 360 LBY2=1,LBY2MX
      MM=LBY2+LBY2+1
      NN=NMAX(LBY2)
      T1=-T1
      DO 350 N=1,NN
      Z=BZ(N,LBY2)
      X=A2PIR-Z
      Y=A2PIR+Z
      IF (ABS(X).GT.0.1)GO TO 340
      Z=Z*A2PIR*SPHBES(MM-2)/(Y*(A2PIR+MM*X))
      GO TO 350
340   Z=Z*SPHBES(MM-1)/(X*Y)
350   RLN(N)=T1*Z
      DO 360 M=1,MM
      J1=J1+1
      DO 360 N=1,NN
      J0=J0+1
      ALMN(J0)=ALMN(J0)+ARRAY(J1)*RLN(N)
360   ALMN(J0+LGTH)=ALMN(J0+LGTH)+ARRAY(J1+465)*RLN(N)
      GO TO 100
C*****END OF MAIN LOOP************************************************
C*****OUTPUT SPHERICAL HARMONIC COEFFICIENTS ALMN
400   WRITE(6,5001) IXTL,NREFL,NREJMX,
     .  (I,NRRANG(I),I=1,NSRANG)
5001  FORMAT(/1X,' Crystal ',I3,8X,I8,' reflections read',
     . 6X,I8,' reflections too large'///
     .'  Numbers of reflections in ranges of sin**2 theta / lambda**2'
     .//30(1X,I6,5X,I6/))
      J0=0
      T0=0.
      DO 410 L=1,LGTH2
410   T0=AMAX1(T0,ABS(ALMN(L)))
      T0=T0*1.0E-8
      DO 420 LBY2=1,LBY2MX
      L=LBY2+LBY2
      MM=L+1
      NN=NMAX(LBY2)
      DO 420 M1=1,MM
      M=M1-1
      DO 420 N=1,NN
      J0=J0+1
      T1=ALMN(J0)
      T2=ALMN(J0+LGTH)
      IF (ABS(T1).LT.T0)T1=0.
      IF (ABS(T2).LT.T0)T2=0.
      IF (LIST.EQ.1)WRITE(6,1050)L,M,N,T1,T2
420   WRITE(TAPOT)L,M,N,T1,T2
C*****WRITE TERMINATOR RECORD
      L=10000
      WRITE(TAPOT)L,M,N,X,Y
      TAPOT=TAPOT+1
C
C---------------------------------------------------------------------------
C End of crystal loop
10    CONTINUE
C
      RETURN
      END
C
C
C
C
C
C
      SUBROUTINE FRFSUM
C     =================
C
C***********************************************************************
C***********************************************************************
C***********************************************************************
C
C           ND500-FRFSUM                  W.KABSCH  10-1982
C
C***********************************************************************
C***********************************************************************
C***********************************************************************
C***** FAST ROTATION FUNCTION SUMMATION PROGRAM
C       COMPARES THE ROTATED PATTERSON DENSITY OF THE SECOND CRYSTAL
C       WITH THE UNROTATED PATTERSON DENSITY OF THE FIRST CRYSTAL. THE
C       ROTATION-ANGLE IS KAPPA ABOUT A ROTATION AXIS U=(UX,UY,UZ),
C       WITH POLAR COORDINATES PHI,OMEGA, SUCH THAT
C       UX=SIN(omega)*COS(phi)
C       UY=SIN(omega)*SIN(phi)
C       UZ=COS(omega).
C       THE PRINTED VALUES OF THE OVERLAP-INTEGRALS APPEAR IN SECTIONS
C       AT CONSTANT KAPPA, EACH LINE AT A CONSTANT VALUE OF PHI.
C       OMEGA INCREASES FROM LEFT TO RIGHT, PHI INCREASES DOWNWARDS.
C
C************************* FILES **************************************
C
C
C 2.)CONTROL PARAMETER CARD                  FORMAT(10I5,E15.5)
C    IPHI1                  STARTING VALUE OF PHI (DEGREES)
C    IPHI2                  MAXIMUM VALUE OF PHI (DEGREES)
C    NPHI                   STEP SIZE
C    IOMEG1                  STARTING VALUE OF OMEG
C    IOMEG2                  MAXIMUM VALUE OF OMEG
C    NOMEG                   STEP SIZE IN OMEG
C    IKAPA1                  STARTING VALUE OF KAPA
C    IKAPA2                  MAXIMUM VALUE OF KAPA
C    NKAPA                   STEP SIZE IN KAPA
C    LIST                    0 : NO LISTING,
C                            1 : LIST EXPANSION COEFFICIENTS ALMN,BLMN.
C    RMAX                    MAXIMUM VALUE OF OVERLAP INTEGRAL.
C                            RMAX MUST BE ZERO , OTHERWISE IT IS
C                            ASSUMED THAT A PREVIOUSLY COMPUTED
C                            MAP NEEDS TO BE PRINTED AGAIN USING THE
C                            KNOWN VALUE OF RMAX.
C
C***********************************************************************
C
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*80 TITLE,FNAME
      COMMON /TEXTC/ TITLE,FNAME
C
      COMMON /RHMOUT/ WRMEAN, WRMS
      SAVE /RHMOUT/
      REAL WRMEAN, WRMS
C
C Maximum section size
      PARAMETER (NSMAP=33759)
C
      CHARACTER FILEN*80
      DOUBLE PRECISION D1,D2,D3,DSINB,DCOTB,DLNSIN,DLNCOT,DLM1M2(31,61),
     @ DBINOM,DFACT(15,31)
      DIMENSION F(4),FF(61,15,4),FRQ(61),FIQ(61),R(360),DBINOM(60,31),
     @ RROT(NSMAP),CCR(61),CCI(61)
      REAL COST(360),SINT(360),FRPQ(31,121),FIPQ(31,121)
C ejd      INTEGER*2 IREC,IBUF
CEJD       DIMENSION REC(4),IREC(11)
C  Hell  - previous unnamed common in CALMN has  all arrays  real.
C  CALMN common:
C      COMMON /POLDAT/ FACTOR(15),FAC1(15,30),FAC2(15,30),BASIS(9),
C     1   RBASIS(9),SPHBES(30),RLN(15),PLM(31),CMPHI(31),SMPHI(31),
C     2   ARRAY(930),ALMN(10000)
C    Now it contains this.....
C  ejd - Changed FACT(..) to DFACT(..); BINOM(..) to DBINOM(..)
C ejd - removed ibuf(..) from COMMON - space not so important now!
C  FRFSUM common:
C      COMMON /POLDT2/ DFACT(15,31),COST(360),SINT(360),DLM1M2(31,61),
C     @ FRPQ(31,121),FIPQ(31,121)
C     @ FRPQ(31,121),FIPQ(31,121),IBUF(2200)
C        EQUIVALENCE (IREC(4),REC(1))
C  ejd Change output routine.. IBUF IREC now integer*4
             DIMENSION AOBUF(1400),REC(4),IREC(3)
C             DIMENSION IOBUF(1400),REC(4),IREC(3)
             COMMON /IOREC/ IBUF(600),RBUF(800)
C             EQUIVALENCE (IOBUF(1),IBUF(1)),(IOBUF(601),RBUF(1))
C
C Map file stuff
      DIMENSION IUVW(3),MXYZ(3),DCELL(6)
C
      DATA ITAPE1,ITAPE2,ITAPE3,IOUT/8,9,64,10/
      DATA TWOPI/6.283185/
C
C
C
      LIST=0
C      LIST=1
      RMAX=0.
      RMIN=0.
      RMEAN=0.
      RRMS=0.
      EPS=1.0E-10
C
C***********************************************************************
C
      IF(LSAVE.GE.0) THEN
            REWIND ITAPE1
            IF(NCRYST.EQ.2) REWIND ITAPE2
      ELSE
C  SUM option, read ALMN files from previous run
            L=LENSTR(FNAME)
            FILEN=FNAME(1:L)//'1'
            IFAIL=0
            CALL CCPDPN(ITAPE1,FILEN,'READONLY','U',1,IFAIL)
            IF(NCRYST.EQ.2) THEN
            FILEN=FNAME(1:L)//'2'
              IFAIL=0
              CALL CCPDPN(ITAPE2,FILEN,'READONLY','U',1,IFAIL)
            ENDIF
      ENDIF
C
C If writing out map file, open it
C Note on header for map file
C   Since this map is not on a crystallographic cell, and the sampling
C   intervals do not have be a factor of 360, the 'sampling intervals'
C   are set to the interval in degrees, the 'cell' is 360,360,360,90,90,90
C   and the start points on each axis are in degrees. However the numbers
C   of points in each direction are correct (as they must be), so the
C   'stop points' are actually the start point + (number of points - 1)
      IF(LMAP.GT.0) THEN
         DO 1,I=1,3
            IUVW(I)=I
            DCELL(I)=360.
            DCELL(I+3)=90.
 1       CONTINUE
         MXYZ(1)=360/NPHI
         MXYZ(2)=360/NOMEG
         MXYZ(3)=360/NKAPA
         NSEC=(IKAPA2-IKAPA1)/NKAPA+1
         NU2=(IPHI2-IPHI1)/NPHI+IPHI1
         NV2=(IOMEG2-IOMEG1)/NOMEG+IOMEG1
         CALL MWRHDR(4,TITLE,NSEC,IUVW,MXYZ,IKAPA1,
     .        IPHI1,NU2,IOMEG1,NV2,DCELL,1,2)
      ENDIF
C
C
C
      LALF=(IPHI2-IPHI1)/NPHI+1
      LBET=1+((IOMEG2-IOMEG1)/NOMEG)
      LGAM=(IKAPA2-IKAPA1)/NKAPA+1
      NPOINT=LALF*LBET*LGAM
C
C          ******************************
      CALL QOPEN(IOUT,'ROTFUN','SCRATCH')
C          ******************************
C
      IF (IOUT .LE. 0) THEN
         CALL CCPERR(1,
     $        '*** POLARRFN: failed to open file ROTFUN ***')
      ENDIF
C
C          ********************************
      CALL QMODE(IOUT,2,NCH)
      CALL QOPEN(ITAPE3,'ITAPE3','SCRATCH')
C          ********************************
C
      IF (ITAPE3 .LE. 0) THEN
         CALL CCPERR(1,
     $        '*** POLARRFN: failed to open file ITAPE3 ***')
      ENDIF
C
C          *******************
      CALL QMODE(ITAPE3,2,NCH)
C          *******************
C
C  Ejd Use IBSCR to count BLOCKS written out;
C          IPSCR to count within each BLOCK.
C
       IPSCR=0
      IBSCR=0
C
C---- Calculate binomial coefficients
C
      DBINOM(1,1)=1.0D+0
C
C
      DO 10 N=1,59
      M1=(N+1)/2
      DBINOM(N+1,1)=1.0D+0
      DBINOM(N,M1+1)=DBINOM(N,N+1-M1)
      DO 10 M=1,M1
10    DBINOM(N+1,M+1)=DBINOM(N,M+1)+DBINOM(N,M)
      DO 20 LBY2=1,15
      L1=LBY2+LBY2+1
      DO 20 M=1,L1
20    DFACT(LBY2,M)=5.0D-1*DLOG(DBINOM(4*LBY2,L1-M+1))
C
C---- set up sin/cos - tables
C
      DO 30 L=1,360
      COST(L)=COS(L/57.29578)
30    SINT(L)=SIN(L/57.29578)
C
C---- form clm1m2 coefficients
C
      L0=0
      LEND=0
40    READ(ITAPE1)L,M,N,(F(J1),J1=1,2)
      IF (L.EQ.10000)GO TO 120
      F(3)=F(1)
      F(4)=F(2)
      IF (NCRYST.EQ.2)READ(ITAPE2)L,M,N,(F(J1),J1=3,4)
      IF (LIST.EQ.1)WRITE(6,1050)L,M,N,F
1050   FORMAT(3I8,4E15.4)
      IF (L.NE.L0)GO TO 70
      IF (N.GT.NMAX)NMAX=N
50    S1=1-2*MOD(M,2)
C
C
      DO 60 J1=1,4
60    FF(M+31,N,J1)=F(J1)
C
C
      FF(31-M,N,3)=S1*F(3)
      FF(31-M,N,4)=-S1*F(4)
      GO TO 40
70    IF (L0.EQ.0)GO TO 100
      L1=L0+1
      L21=L0+L1
C
C
      DO 90 M1=L1,L21
      MM1=M1-L1
      J1=31+MM1
      DO 80 M2=1,L21
      CCR(M2)=0.
      CCI(M2)=0.
      J2=31+M2-L1
C
C
      DO 80 N1=1,NMAX
      CCR(M2)=CCR(M2)+FF(J1,N1,1)*FF(J2,N1,3)+FF(J1,N1,2)*FF(J2,N1,4)
80    CCI(M2)=CCI(M2)+FF(J1,N1,1)*FF(J2,N1,4)-FF(J1,N1,2)*FF(J2,N1,3)
C
C
      DO 90 M2=1,L21
      MM2=M2-L1
      IREC(1)=L0
      IREC(2)=MM1
      IREC(3)=MM2
      REC(1)=CCR(M2)
      REC(2)=CCI(M2)
      REC(3)=CCR(L1-MM2)
      REC(4)=CCI(L1-MM2)
C
C
        IBUF(3*IPSCR + 1) = IREC(1)
        IBUF(3*IPSCR + 2) = IREC(2)
        IBUF(3*IPSCR + 3) = IREC(3)
C
C
       DO 140 J3=1,4
140    RBUF(4*IPSCR+J3)=REC(J3)
C
C
       IPSCR=IPSCR+1
       IF (IPSCR.NE.200)GO TO 90
C
C
        DO 141 J3=1,600
141     AOBUF(J3)=IBUF(J3)
        DO 142 J3=1,800
142     AOBUF(600+J3)=RBUF(J3)
      CALL QWRITE(ITAPE3,AOBUF,1400)
C
C
      IPSCR=0
      IBSCR=IBSCR + 1
90    CONTINUE
100   IF (LEND.EQ.1)GO TO 130
      L0=L
      NMAX=N
      DO 110 M1=1,61
      DO 110 N1=1,15
      DO 110 J1=1,4
110   FF(M1,N1,J1)=0.
      GO TO 50
120   CONTINUE
      LEND=1
      GO TO 70
130   IBUF(3*IPSCR+1)=10000
C
C
        DO 131 J3=1,600
131     AOBUF(J3)=IBUF(J3)
        DO 132 J3=1,800
132     AOBUF(600+J3)=RBUF(J3)
      CALL QWRITE(ITAPE3,AOBUF,1400)
C
C----  Ejd  Keep count of number of records.
C
        NBSCR = IBSCR + 1
      CLOSE(ITAPE1)
      IF (NCRYST.EQ.2)CLOSE(ITAPE2)
C
      LMAX1=L0+1
      L2MAX=L0+L0
      L2MAX1=LMAX1+L0
      L4MAX1=4*L0+1
C
C---- Calculate rotation function at given eulerian angles
C
      IB=0
      IOMEG=IOMEG1
170   CONTINUE
C
C----  Ejd  - Now to unpack things..
C                IPSCR=2200-11
C
      IBSCR=-1
C
C
      CALL QSEEK(ITAPE3,1,1,1)
C
C
      IF ((IOMEG.EQ.0).OR.(IOMEG.EQ.180))GO TO 180
      D1=DBLE(IOMEG)/1.145915590261646D+2
      DSINB=DSIN(D1)
        IF (ABS(DSINB).LT.EPS) WRITE(6,3406) DSINB
 3406   FORMAT(' DSINB ',E13.5)
      DCOTB=DCOS(D1)/DSINB
      DLNSIN=DLOG(DSINB)
      DLNCOT=DLOG(DCOTB)
C
C
180   DO 190 IP=1,LMAX1
      DO 190 IQ=1,L4MAX1
      FRPQ(IP,IQ)=0.
190   FIPQ(IP,IQ)=0.
C
C---- sum over l and m2 to form fpq
C
      L0=-1
C
C
200    CONTINUE
C
C
       DO 160 IBSCR = 1,NBSCR
       CALL QREAD(ITAPE3,AOBUF,1400,IEND)
        DO 161 J3=1,600
161     IBUF(J3) = NINT(AOBUF(J3))
        DO 162 J3=1,800
162     RBUF(J3) = AOBUF(600+J3)
       IF(IEND.NE.0) CALL CCPERR(1,  ' stop ITAPE3')
C
C150   CONTINUE
      DO 155 IRC = 1,200
C
C      Unpack array IOBUF
C
       IREC(1) = IBUF(3*IRC -2)
       IREC(2) = IBUF(3*IRC -1)
       IREC(3) = IBUF(3*IRC )
       REC(1) = RBUF(4*IRC -3)
       REC(2) = RBUF(4*IRC -2)
       REC(3) = RBUF(4*IRC -1)
       REC(4) = RBUF(4*IRC )
       L=IREC(1)
      IF (L.EQ.10000)GO TO 400
      M1=IREC(2)
      M2=IREC(3)
      S1=REC(1)
      S2=REC(2)
      S3=REC(3)
      S4=REC(4)
      IF (L.EQ.L0)GO TO 210
C
C  All this done once only.
C
230   L0=L
      L1=L0+1
      L21=L0+L1
      IF (IOMEG.EQ.0)GO TO 290
      IF (IOMEG.EQ.180)GO TO 320
C
C----  calculate starting values
C
      LBY2=L0/2
      D1=DBLE(L21-1)*DLNSIN
      DO 240 MM2=1,L21
      J1=1+IABS(MM2-L1)
      D2=D1+DBLE(MM2-1)*DLNCOT+DFACT(LBY2,J1)
      D3=0.0D+0
      IF (D2.GT.-160.0D+0)D3=DEXP(D2)
240   DLM1M2(L1,MM2)=D3
C
C---- apply recurrence relation
C
      DO 250 MM1=1,L0
      MMM1=L1-MM1
      M21=MMM1+MMM1-1
      D1=DSQRT(DBLE((L0+MMM1)*MM1))
        IF (ABS(D1).LT.EPS) WRITE(6,3414) D1
 3414   FORMAT(' D1 ',E13.5)
      DO 250 MM2=1,M21
      MMM2=MM2-MMM1
250   DLM1M2(MMM1,MMM2+L1)=(DBLE(MMM1-MMM2)*DCOTB*
     1 DLM1M2(MMM1+1,MMM2+L1)-DSQRT(DBLE((L0+MMM2)*(L1-MMM2)))*
     2 DLM1M2(MMM1+1,MMM2+L0))/D1
C
C---- extend by symmetry
C
      J1=-1
      DO 270 MM1=1,L0
      J2=J1
      DO 260 MM2=1,MM1
      DLM1M2(MM2,L1-MM1)=DLM1M2(MM1+1,L1-MM2+1)
      DLM1M2(MM2,L1+MM1)=DLM1M2(MM1+1,L0+MM2)*DBLE(J2)
260   J2=-J2
270   J1=-J1
      DO 280 MM1=1,L1
      DO 280 MM2=1,L21
      IF (DABS(DLM1M2(MM1,MM2)).LT.1.0D-06)DLM1M2(MM1,MM2)=0.0D+00
280   CONTINUE
      GO TO 210
C
C---- special case omega=0
C
290   DO 310 MM1=1,L1
      DO 300 MM2=1,L21
300   DLM1M2(MM1,MM2)=0.0D+0
310   DLM1M2(MM1,MM1+L0)=1.0D+0
      GO TO 210
C
C----special case omega=180
C
320   D1=1.0D+0
      DO 340 MM1=1,L1
      DO 330 MM2=1,L21
330   DLM1M2(MM1,MM2)=0.0D+0
      DLM1M2(MM1,L1+1-MM1)=D1
340   D1=-D1
      GO TO 210
C
210   Y=1.
      IF (MOD(M1+M2,2).NE.0)Y=-1.
      J1=M1-M2+L2MAX1
      J2=L2MAX1-M1-M2
      DO 220 IP=1,L1
      X=SNGL(DLM1M2(IP,L1+M1)*DLM1M2(IP,L1+M2))
      FRPQ(IP,J1)=FRPQ(IP,J1)+S1*X
      FIPQ(IP,J1)=FIPQ(IP,J1)+S2*X
      IF (M1.EQ.0)GO TO 220
      X=SNGL(DLM1M2(IP,L1-M1)*DLM1M2(IP,L1+M2))*Y
      FRPQ(IP,J2)=FRPQ(IP,J2)+S3*X
      FIPQ(IP,J2)=FIPQ(IP,J2)-S4*X
220   CONTINUE
C SUM never used???
      SUM = FRPQ(1,J2)+FIPQ(1,J2)
155    CONTINUE
160    CONTINUE
C
C---- sum over p to form fq at given omega and kappa
C
400   CONTINUE
      IKAPA=IKAPA1
410   DO 430 IQ=1,L2MAX1
      X=FRPQ(1,L2MAX+IQ)
      Y=FIPQ(1,L2MAX+IQ)
      J1=IKAPA
      DO 420 IP=2,LMAX1
        J2=MOD(J1,360)
        IF (J2.LE.0)J2=J2+360
        S1=FRPQ(IP,IQ+L2MAX)
        S2=FRPQ(IP,L2MAX1-IQ+1)
        S3=FIPQ(IP,IQ+L2MAX)
        S4=FIPQ(IP,L2MAX1-IQ+1)
        X=X+(S1+S2)*COST(J2)+(S3+S4)*SINT(J2)
        Y=Y+(S3-S4)*COST(J2)-(S1-S2)*SINT(J2)
        J1=J1+IKAPA
420   CONTINUE
      FRQ(IQ)=X
430   FIQ(IQ)=Y
C
C---- Final sum over q
C
      L=0
500   IPHI=IPHI1
510   X=FRQ(1)*0.5
      J1=IPHI
      DO 520 IQ=2,L2MAX1
        J2=MOD(J1,360)
        IF (J2.LE.0) J2=J2+360
        X=X+FRQ(IQ)*COST(J2)+FIQ(IQ)*SINT(J2)
        J1=J1+IPHI
520   CONTINUE
      X=2.*X
      IF (RMAX.LT.X) RMAX=X
      RMEAN = RMEAN+X
      RRMS = RRMS+X*X
      L=L+1
      R(L)=X
      IPHI=IPHI+NPHI
      IF (IPHI.LE.IPHI2)GO TO 510
C
      CALL QWRITE(IOUT,R,LALF)
C
      IB=IB+1
      IKAPA=IKAPA+NKAPA
      IF (IKAPA.LE.IKAPA2)GO TO 410
      IOMEG=IOMEG+NOMEG
      IF (IOMEG.LE.IOMEG2)GO TO 170
C
      WRITE(6,1010) RMAX
1010  FORMAT(/' Maximum R = ',E15.4/)
      RSCALE=100./RMAX
      RMEAN=RMEAN*RSCALE
      RRMS=SQRT(RRMS*RSCALE*RSCALE/NPOINT)
C
      RWT = 0.0
      WRMS = 0.0
      WRMEAN = 0.0
C
C
C---- Map output
C
600   JKAPA=IKAPA1
      NSKAP=(IKAPA2-IKAPA1)/NKAPA+1
      DO 610,JKAPA=IKAPA1,IKAPA2,NKAPA
        WW = 1.0 - COS(JKAPA* TWOPI/360.0)
        IF(JKAPA.EQ.180.0) WW = 0.5*WW
        IB=(JKAPA-IKAPA1)/NKAPA+1
        L=0
        DO 620,IOMEG=IOMEG1,IOMEG2,NOMEG
        WV = WW*SIN(IOMEG* TWOPI/360.0)
        IF(IOMEG.EQ.180) WV = 0.5*WV
          CALL QSEEK(IOUT,IB,1,LALF)
          CALL QREAD(IOUT,R,LALF,IEND)
          IF(IEND.NE.0) CALL CCPERR(1, 'stop IOUT')
          IB=IB+NSKAP
          DO 640 J1=1,LALF
            WU=WV
            IF (J1.EQ.1 .OR. (J1-1)*NPHI.EQ.180) WU=.5*WU
            AJ2=R(J1)*RSCALE
            WRMEAN = WRMEAN + WU*AJ2
            WRMS = WRMS + AJ2*AJ2*WU
            RWT = RWT + WU
            IF (AJ2.LT.-99.) AJ2=-99
            RMIN=MIN(RMIN,AJ2)
            RROT(J1+L)=AJ2
640       CONTINUE
          L=L+LALF
620     CONTINUE
        IF(LPRINT.GT.0) CALL PRNSEC(RROT,LALF,LBET,JKAPA)
C
C Find peaks in section, plot etc - NB MAPSEC calls s/r PEAKPK
C which uses RRMS to select peaks
C
        CALL MAPSEC(RROT,LALF,LBET,JKAPA)
C end of section output loop
610   CONTINUE
C
C
      WRMEAN = WRMEAN/RWT
      WRMS = WRMS/RWT - WRMEAN**2
      IF (WRMS .GE. 0.0) THEN
         WRMS = SQRT(WRMS)
      ELSE
         WRMS = 0.0
      ENDIF
C Search last section for peaks
      CALL MAPSEC(RROT,LALF,LBET,-1)
C Print out peaks
      CALL PPEAKS
C
C Print weighted rms and mean: these also go into map header
      WRITE(6,FMT='(A,F7.3/A,F7.3)') ' Weighted mean of map',WRMEAN,
     1                               ' Weighted rms',WRMS
C
      CALL QCLOSE(IOUT)
      CALL QCLOSE(ITAPE3)
      CALL CCPERR(0,' Normal Termination from Polarrfn.for')
      END
C
      SUBROUTINE PRNSEC(RROT,LALF,LBET,JKAPA)
C     =========================================
C
C Print out section
C
      DIMENSION RROT(*)
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
C
      INTEGER OMEGMK(100),PHIMK(100)
      CHARACTER*3 IT,LPRNT(60)
C
      SAVE IFIRST,NBMK
      DATA IFIRST/0/
C
C
C
1020  FORMAT(//,' Rotation angle kappa =',I3/
     @ '   polar coordinates of rotation '
     @,'Axes are phi (vertical) and omega (horizontal)'/)
 1028   FORMAT(4X,'P'/4X,'h'/4X,'i',10X,'O m e g a'/)
1030  FORMAT(7X,42I3)
1031  FORMAT(6X,42A)
 1057  FORMAT(1X,I4,'> ',60A)
C
C
      IF(IFIRST.EQ.0) THEN
C--Setting arrays PHIMK() and OMEGMK() for marker lines.
        NBMK=1+((IOMEG2-IOMEG1)/NOMEG)
        DO 2146 IMK=1,NBMK
          OMEGMK(IMK)=IOMEG1+NOMEG*(IMK-1)
 2146     CONTINUE
        NAMK=1+((IPHI2-IPHI1)/NPHI)
        DO 2148 IMK=1,NAMK
          PHIMK(IMK)=IPHI1+NPHI*(IMK-1)
 2148     CONTINUE
        IFIRST=1
      ENDIF
C
C
C  Write out section
C
      L=LALF*(LBET-1)+1
      SMAX=-1000.
C
      WRITE(6,1020)JKAPA
        WRITE(6,1028)
        WRITE(6,1030) (OMEGMK(IMK),IMK=1,NBMK)
        IT='  V'
        WRITE(6,1031) (IT,IMK=1,NBMK)
      DO 660 J1=1,LALF
           KPR=0
           DO 2341 J2=1,L,LALF
                R=RROT(J1+J2-1)
                IF(R.GT.SMAX) THEN
                        SMAX=R
                        J1M=J1
                        J2M=(J2-1)/LALF
                ENDIF
                N9=NINT(R)
                  IF(N9.LT.0) THEN
                        IT='  -'
                  ELSEIF(N9.LT.LPRINT) THEN
                        IT='  .'
                  ELSE
                        WRITE(IT,2009) N9
2009            FORMAT(I3)
                  ENDIF
                  KPR=KPR+1
                  LPRNT(KPR)=IT
 2341      CONTINUE
 660   WRITE(6,1057) PHIMK(J1),(LPRNT(J2),J2=1,KPR)
        IT='  ^'
        WRITE(6,1031) (IT,IMK=1,NBMK)
        WRITE(6,1030) (OMEGMK(IMK),IMK=1,NBMK)
C
      J1=(J1M-1)*NPHI+IPHI1
      J2=J2M*NOMEG+IOMEG1
      WRITE(6,1035) SMAX,J2,J1
1035  FORMAT(/'    Maximum on this section = ',F8.1,'   at omega =',I4,
     . '  phi =',I4)
      RETURN
      END
C
C
      SUBROUTINE MAPRED
C     =================
C
C Read back file of rotation function for plotting etc
C Use sections KRKAP1 to KRKAP2
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*80 TITLE,FNAME,QTITLE
      COMMON /TEXTC/ TITLE,FNAME
C
C Maximum section size
      PARAMETER (NSMAP=33759)
C
      DIMENSION IUVW(3),MXYZ(3),DCELL(6),RROT(NSMAP)
C
C Open map file
      CALL MRDHDR(3,'MAPIN',QTITLE,NSEC,IUVW,MXYZ,IKAPA1,
     .  IPHI1,NU2,IOMEG1,NV2,DCELL,LS,LMODE,RMIN,RMAX,RMEAN,RRMS)
C
C Keep input title if specified
      IF(TITLE.EQ.' ') TITLE=QTITLE
C Set parameters from file. See FRFSUM for comments on what is in map header
      NPHI=360/MXYZ(1)
      NOMEG=360/MXYZ(2)
      NKAPA=360/MXYZ(3)
      NU=NU2-IPHI1+1
      NV=NV2-IOMEG1+1
      IKAPA2=IKAPA1+(NSEC-1)*NKAPA
      IPHI2=IPHI1+(NU-1)*NPHI
      IOMEG2=IOMEG1+(NV-1)*NOMEG
      IF(KRKAP1.LT.-999) KRKAP1=IKAPA1
      IF(KRKAP2.LT.-999) KRKAP2=IKAPA2
C Check for maximum section size
      NSECSZ=NU*NV
      IF(NSECSZ.GT.NSMAP) THEN
                  WRITE(6,9050) NSECSZ,NSMAP
9050      FORMAT(//' Section size too large, =',I8,
     .             '    maximum=',I8///)
          CALL CCPERR(1,' error on section 9050')
      ENDIF
C
C Read sections KRkap1 to KRkap2
C
C   Corrected next line 3/3/87 PRE
      JSEC=(KRKAP1-IKAPA1)/NKAPA+IKAPA1
      CALL MPOSN(3,JSEC)
      DO 10,JKAPA=KRKAP1,KRKAP2,NKAPA
      CALL MGULP(3,RROT,IER)
      IF(IER.NE.0) GO TO 11
C Print section
      IF(LPRINT.GT.0) CALL PRNSEC(RROT,NU,NV,JKAPA)
C Plot and find peaks as required
      CALL MAPSEC(RROT,NU,NV,JKAPA)
10    CONTINUE
C
11    CALL MAPSEC(RROT,NU,NV,-1)
      CALL PPEAKS
      RETURN
      END
C
C
      SUBROUTINE MAPSEC(RROT,LA,LB,IG)
C     ================================
C
C  Deal with section kappa = IG in RROT.
C
C   1. write out map section to map file if LMAP .gt. 0
C
C   2. plot if LPLOT .ne. 0
C
C   3. search for peaks in PREVIOUS section: if IG=-1 search last section
CEJD LA - number of points along Phi
CEJD LB - number of points along Omega - iomeg1,ipmeg2,nomeg
CEJD IG = Kappa
C
C
      DIMENSION RROT(LA,LB)
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
      COMMON /RHMOUT/ WRMEAN, WRMS
      SAVE /RHMOUT/, RMAP
      REAL WRMEAN, WRMS
C
C Maximum section size
      PARAMETER (NSMAP=33759)
C
      DIMENSION PLTMAP(NSMAP),RMAP(3*NSMAP)
C
C
      IF(IG.LT.0) THEN
         IF (LMAP .GT. 0) THEN
C Close map file if required
            CALL MCLOSC(4,RMIN,100., WRMEAN, WRMS)
         ENDIF
C Close plot file if required
         IF(LPLOT.NE.0) CALL PLTSTP
         GO TO 10
      ENDIF
C
      IF(LMAP.GT.0) THEN
C Write map section to map file
            CALL MSPEW(4,RROT)
      ENDIF
C
      IF(LPLOT.NE.0) THEN
C  Plot section
            NLA=360/NPHI+1
            NLB=90/NOMEG+1
            IF(NLA*NLB.GT.NSMAP) CALL CCPERR(1,  'stop NSMAP!!!')
            CALL PLOTMP(RROT,LA,LB,IG,PLTMAP,NLA,NLB)
      ENDIF
C
C Pick peaks on PREVIOUS section
10    IF(LA*LB.GT.NSMAP) CALL CCPERR(1, 'stop-10 NSMAP!!!')
      CALL PEAKPK(RROT,LA,LB,IG,RMAP)
C
      RETURN
      END
C
C
C
      SUBROUTINE PLOTMP(RROT,LA,LB,IG,PLTMAP,NLA,NLB)
C     ============================================
C
C Plot section as stereogram
C
C First generate array PLTMAP, using symmetry as necessary,
C then plot (this part of program taken from RFPLOT)
C  If LPLOT .lt. 0, don't plot any labels
C
C In the absence of symmetry, each kappa section is plotted as
C two hemispheres omega = 0 to 90, and 90 to 180. Subroutine
C CHKSYM works out the symmeytry of the rotation function, and sets
C flag LSYMFG to indicate whether one or two hemispheres need to
C be plotted. LSYMFG = 2  no symmetry, plot both hemispheres
C                    = 1  symmetry sufficient to reduce asymmetric unit
C                         to one hemisphere
C                    =-2  asymmetric unit omega = 0 to 180, phi = 0 to 180
C                         plot two quarter spheres in same circle
C
C
      DIMENSION RROT(LA,LB),PLTMAP(NLB,NLA)
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
      CHARACTER*80 TITLE,FNAME
      COMMON /TEXTC/ TITLE,FNAME
C
      PARAMETER (IPTMAX=200)
      DIMENSION XXX(IPTMAX),YYY(IPTMAX)
C Maximum section size
      PARAMETER (NSMAP=33759)
C
      LOGICAL BITS(2*NSMAP)
C
      CHARACTER*80 LINE
      SAVE IFIRST
C
C Set up plotting viewport and window
C
C  Relative character sizes for tick marks and titles
      REAL  CSZ1, CSZ2
      REAL VWPORT(2,2), WINDOW(2,2)
C    WINDOW  are the corresponding coordinates in user units, which
C           set radius of stereogram = 1
C                   x1    y1    x2   y2
      DATA WINDOW/ -1.3, -1.3, 1.3, 1.3/
C
      DATA  CSZ1, CSZ2/0.8,1.1/
C
C    VWPORT  are corners of plotting space in range 0-1 on x & y
C                   x1    y1    x2    y2
      DATA VWPORT/ 0.05, 0.05, 0.95, 0.95/
      DATA IFIRST/0/
C
      CHSZ = 0
C
      IF(IOMEG1.LT.0.OR.IOMEG2.GT.180) THEN
            WRITE(6,1001) IOMEG1,IOMEG2
1001      FORMAT(/' Plotting routine cannot cope with omega outside',
     .        ' range 0 to 180, omega limits =',2I5/)
            RETURN
      ENDIF
C
C On first entry, start plotting
C
      IF(IFIRST.EQ.0) THEN
            CALL PLTINI
C Set plotting transform
            CALL PLTWIN(VWPORT, WINDOW)
C Find rotation function symmetry from crystal symmetries
            CALL CHKSYM
      ENDIF
C
C
C If necessary, loop for two halves of kappa section
C   LSYMFG = 1,2 or -2 see above
C Kappa = 180 has symmetry type 1 (180-omega, 180+phi)
      NHALF=IABS(LSYMFG)
      IF(IG.EQ.180) THEN
            NHALF=1
C Add symmetry type 1 if not already present
            IF(NCRYST.EQ.2) THEN
                  MRSYM=MRSYM+1
                  CALL SETSYM(JRSYM(1,1,MRSYM),1)
            ENDIF
      ENDIF
C
      DO 200,LHALF=1,NHALF
C
C Copy appropriate part of RROT to PLTMAP using symmetry as necessary
C First define limits required on omega(omeg) and phi(phi) in PLTMAP
      IF(LHALF.EQ.1) THEN
            IB1=0
            IB2=90
      ELSE
            IB1=90
            IB2=180
      ENDIF
      IF(LSYMFG.LT.0.AND.IG.NE.180) THEN
C LSYMFG = -2, plot 2 half hemispheres in same circle
            MLA=(NLA+1)/2
            IA1=0
            IA2=180
            LSPLIT=1
      ELSE
            MLA=NLA
            IA1=0
            IA2=360
            LSPLIT=0
      ENDIF
C
Cejd      CALL ZEROI(PLTMAP,MLA*NLB)
      DO 201 JPLT = 1,NLB
      DO 201 IPLT = 1,MLA
201      PLTMAP(JPLT,IPLT) = 0
C
C Loop points in output array PLTMAP
      KB=1
      DO 11,IB=IB1,IB2,NOMEG
      KA=1
      DO 12,IA=IA1,IA2,NPHI
C Put point IA,IB into calculated area if possible: if not JA=0
      CALL SYMPUT(JA,JB,IA,IB)
      IF(JA.GT.0) PLTMAP(KB,KA)=RROT(JA,JB)
      KA=KA+1
12    CONTINUE
      KB=KB+1
11    CONTINUE
C
C Start new picture except for LSYMFG=-2
      IF(LHALF.EQ.1.OR.LSYMFG.GT.1) THEN
            IF (IFIRST .NE. 0) CALL PLTPIC
            IFIRST=1
C
C
C   Now draw the ring around the edge, using a do loop of drawto's
C
            CONV=360./6.28318
            CALL PLTMVU(1.,0.)
            DO 110 K=1,101
            FI=(3.60*FLOAT(K-1))/CONV
            CALL PLTDWU(COS(FI),SIN(FI))
  110       CONTINUE
C
C  Plot tick marks
            TICK=0.01
            CALL PLTMVU(1.,0.)
            CALL PLTDBU(TICK,0.)
            IF(LPLOT.GT.0) THEN
                  CALL PLTTNF(1.+2.*TICK,0.,X,Y)
                  CALL PLTCTX('0.0',X,Y,CSZ1,1,0.0)
            ENDIF
C
            CALL PLTMVU(0.,1.)
            CALL PLTDBU(0.,TICK)
            IF(LPLOT.GT.0) THEN
                  CALL PLTTNF(0.,1.+4.*TICK,X,Y)
                  IF(LSPLIT.GT.0) THEN
                     CALL PLTCTX(
     .                  '90.0, omega 0 to 90',X,Y,CSZ1,1,0.0)
                  ELSE
                     CALL PLTCTX('90.0',X,Y,CSZ1,2,0.0)
                  ENDIF
            ENDIF
C
            CALL PLTMVU(-1.,0.)
            CALL PLTDBU(-TICK,0.)
            IF(LPLOT.GT.0) THEN
                  CALL PLTTNF(-1.-2.*TICK,0.,X,Y)
                  CALL PLTCTX('180.0',X,Y,CSZ1,3,0.0)
            ENDIF
C
            CALL PLTMVU(0.,-1.)
            CALL PLTDBU(0.,-TICK)
            IF(LPLOT.GT.0) THEN
                  CALL PLTTNF(0.,-1.-4.*TICK-CHSZ,X,Y)
                  IF(LSPLIT.GT.0) THEN
                        CALL PLTCTX(
     .                   '90.0, omega 90 to 180',X,Y,CSZ1,1,0.0)
                  ELSE
                        CALL PLTCTX('-90.0',X,Y,CSZ1,2,0.0)
                  ENDIF
            ENDIF
C
C Draw line across section
            IF(LSPLIT.GT.0) THEN
                  CALL PLTMVU(1.,0.)
                  CALL PLTDWU(-1.,0.)
            ENDIF
C  Plot title centred on 0,-1.1
            IF(LPLOT.GT.0) THEN
                  CALL PLTTNF(0.,-1.14,X,Y)
                  CALL PLTCTX(TITLE,X,Y,CSZ2,2,0.0)
                  XX=0.
                  IF(LSYMFG.EQ.2) XX=-0.8
                  CALL PLTTNF(XX,-1.24,X,Y)
                  WRITE (LINE, 6001) IG
6001              FORMAT('Section kappa =  ',I4)
                  L = LENSTR(LINE) + 1
                  IF(LSYMFG.EQ.2) THEN
                        LINE(L:) =
     .                   ', omega = 0 to 90, omega = 90 to 180'
                  ENDIF
                  L = LENSTR(LINE)
                  CALL PLTCTX(LINE(1:L),X,Y,CSZ2,2,0.0)
            ENDIF
      ENDIF
C
C Plot join lines (if any)
      CALL PLJOIN(IG,NPHI,NOMEG)
C
C Set solid lines for contour
      CALL PLTDSH(0,0.,0.,0.)
C
      CALL HPCNTR(PLTMAP,NLB,MLA,BITS,CONTRS,NCONTR,XXX,YYY,IPTMAX)
C
C
200   CONTINUE
C
C
C- - - - -- - - - - - - -    - - - - - - - - - - - - - - - -
C Remove special symmetry for kappa=180
      IF(NCRYST.EQ.2.AND.IG.EQ.180) MRSYM=MRSYM-1
C
      RETURN
      END
C
C
      SUBROUTINE PLJOIN(IG,NPHI,NOMEG)
C     =================================
C
C Plot JOIN lines for section kappa=IG
C  NPHI  sampling interval in degrees along phi
C  NOMEG sampling interval in degrees along omega
C
C Lines are drawn between specified points along the great circle
C JOIN stuff
      PARAMETER (MXJOIN=30)
      COMMON /VJOIN/ NJOIN,PEAKJ(3,2,MXJOIN),DASH(2,MXJOIN)
      SAVE /VJOIN/
C PEAKJ(I,J,K) contains peaks
C  I=1,2,3  omega, phi, kappa
C  J=1,2    point 1, point 2  (kappa2=kappa1)
C
      DIMENSION P(3),Q(3),R(3),T(3,3),TINV(3,3),V1(3),V2(3)
C
C Line buffer
      PARAMETER (MAXLIN=200)
      DIMENSION X(MAXLIN),Y(MAXLIN)
C
C
      CONV=ATAN(1.)/45.
C
      DO 1,K=1,NJOIN
C Check all joins, is this one in current kappa section?
      IF(NINT(PEAKJ(3,1,K)).NE.IG) GO TO 1
C
C We have a peak pair K
C  Direction cosines of points 1 & 2
      CALL DRCCOS(V1,PEAKJ(1,1,K))
      CALL DRCCOS(V2,PEAKJ(1,2,K))
C Define coordinate frame pqr with great circle perpendicular to z
C Put p axis along x
      CALL VSET(P,V1)
C r axis orthogonal to V1 and V2
      CALL CROSS(R,V1,V2)
      CALL UNIT(R)
C q axis completes set
      CALL CROSS(Q,R,P)
      CALL UNIT(Q)
C Transformation matrix T takes point in pqr frame to xyz
      CALL VSET(T(1,1),P)
      CALL VSET(T(1,2),Q)
      CALL VSET(T(1,3),R)
C Invert (transpose) T
      CALL TRANSP(TINV,T)
C
C Theta is angle in great circle from p (V1) towards q (v2)
C  Point 1 (V1) is at theta=0
C  Get theta for point 2
      CALL MATVEC(V1,TINV,V2)
      THETA2=ATAN2(V1(2),V1(1))
C Stepsize in radians
      VSTEP=SIGN(0.02,THETA2)
      N=NINT(ABS(THETA2)/VSTEP)
      VSTEP=THETA2/FLOAT(N)
      NSTEP=N+1
      IF(NSTEP.GT.MAXLIN) THEN
            WRITE(6,1001) MAXLIN
1001      FORMAT(' More than',I6,'lines in PLJOIN')
            NSTEP=MAXLIN
      ENDIF
C
C Loop to draw line between points
      DO 10,N=1,NSTEP
      THETA=(N-1)*VSTEP
      P(1)=COS(THETA)
      P(2)=SIN(THETA)
      P(3)=0.0
      CALL MATVEC(V1,T,P)
      SOMEGA=SQRT(V1(1)**2+V1(2)**2)
C Store lines, x = omega, y= phi, stored as grid units to match what comes
C out of the contour routine
      X(N)=ATAN2(SOMEGA,V1(3))/(CONV*NOMEG)
      IF(ABS(V1(1)).LT.0.000001.AND.ABS(V1(1)).LT.0.000001) THEN
            Y(N)=0.0
      ELSE
            Y(N)=ATAN2(V1(2),V1(1))/(CONV*NPHI)
      ENDIF
10    CONTINUE
C
C Set dash parameters
      IF(DASH(1,K).EQ.0.) THEN
C Solid
            M=0
      ELSE
C Dashed
            M=1
      ENDIF
      CALL PLTDSH(M,DASH(1,K),DASH(2,K),0.)
C Plot
      CALL LINES(X,Y,NSTEP)
C
1     CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE DRCCOS(V,P)
C     ======================
C
C Set dircetion cosines V of axis along omega = P(1), phi = P(2)
C
      DIMENSION V(3),P(2)
C
      CONV=ATAN(1.)/45.
C
      OMEGA=P(1)*CONV
      PHI=P(2)*CONV
C
      SOMEGA=SIN(OMEGA)
      V(1)=SOMEGA*COS(PHI)
      V(2)=SOMEGA*SIN(PHI)
      V(3)=COS(OMEGA)
      RETURN
      END
C
C
C
      SUBROUTINE CHKSYM
C     =================
C
C Check symmetry of the two crystals and set up symmetry of rotation function
C
C Sets MRSYM  number of symmetry operations
C      JRSYM  rotation function symmetry operators, matrices to multiply
C             (omega   phi  1) in degrees.
C
C Possible operations arise as follows:
C
C 1.  180-omega, 180+phi      self-rotation
C
C 2.  180-omega,  -phi        x dyads in both crystals
C
C 3.  180-omega, 180-phi      y dyads in both crystals
C
C 4.  omega,     180+phi      z dyads in both crystals
C
C Other operations may arise as combinations of these
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
      DIMENSION JRMAT(3,3),LDYAD(24,2),LPERM(4)
      DATA LPERM/2,3,1,2/
C
C Starting flags for no symmetry
      LSYMFG=2
      MRSYM=0
C
C Is this a self rotation function?
      IF(NCRYST.EQ.1) THEN
            MRSYM=1
            LSYMFG=1
            CALL SETSYM(JRSYM(1,1,1),1)
      ENDIF
C
C Find out which crystal symmetry operations are dyads
      DO 10,K=1,NCRYST
      DO 11,L=1,NSYM(K)
      LDYAD(L,K)=0
C A dyad matrix is diagonal with  Trace = -1
      TRACE=0.
      DO 12,I=1,3
12    TRACE=TRACE+SYM(I,I,L,K)
      IF(ABS(TRACE+1.).GT.0.0001) GO TO 11
C  Trace=-1, check off diagonals
      M=0
      DO 13,I=1,3
      IF(SYM(I,I,L,K).GT.0.) M=I
      DO 13,J=1,3
      IF(I.NE.J) THEN
            IF(ABS(SYM(I,J,L,K)).GT.0.0001) GO TO 11
      ENDIF
13    CONTINUE
C Yes we have a dyad along axis M
C  Allow for axis permutation from NCODE = 2 or 3
      IF(NCODES(K).EQ.2) THEN
            M=LPERM(M+1)
      ELSEIF(NCODES(K).EQ.3) THEN
            M=LPERM(M)
      ENDIF
      LDYAD(L,K)=M
C If self-rotation, dyads must be paralle in the two crystals
      IF(NCRYST.EQ.1) THEN
C  so add symmetry to list
            MRSYM=MRSYM+1
            CALL SETSYM(JRSYM(1,1,MRSYM),M+1)
      ENDIF
11    CONTINUE
10    CONTINUE
C
C Now check for parallel dyads in the 2 crystals for cross-rotation
      IF(NCRYST.EQ.2) THEN
C  Loop 1st crystal operations
            DO 21,L=1,NSYM(1)
            IF(LDYAD(L,1).GT.0) THEN
C   Check 2nd crystal
                  DO 22,K=1,NSYM(2)
                  IF(LDYAD(K,2).EQ.LDYAD(L,1)) THEN
C    We have parallel dyads, so add symmetry
                        M=LDYAD(L,1)+1
                        MRSYM=MRSYM+1
                        CALL SETSYM(JRSYM(1,1,MRSYM),M)
C    Set asymmetric unit flag
                        IF(M.LE.3) THEN
                              LSYMFG=1
                        ELSEIF(LSYMFG.EQ.2) THEN
                              LSYMFG=-2
                        ENDIF
                  ENDIF
22                CONTINUE
            ENDIF
21          CONTINUE
      ENDIF
C
C
C We now have MRSYM symmetry operations: see if combinations generate
C any new ones
      IF(MRSYM.LE.1) RETURN
C Loop pairs of operations: any new ones generated need not be tested
      MRS=MRSYM
      DO 100,K=1,MRS-1
      DO 110,L=K+1,MRS
      CALL MATMLI(JRMAT,JRSYM(1,1,L),JRSYM(1,1,K))
C 'Translations' can have 360 added or subtracted to bring into range 0 to 180
      DO 111,I=1,2
      J=JRMAT(I,3)
      IF(J.GT.359) J=J-360
      IF(J.LT.-179) J=J+360
      JRMAT(I,3)=J
111   CONTINUE
C Check if we already have the same matrix
      DO 121,M=1,MRSYM
      DO 120,I=1,3
      DO 120,J=1,3
      IF(JRMAT(I,J).NE.JRSYM(I,J,M)) GO TO 121
120   CONTINUE
C Same as this matrix M so ignore
      GO TO 110
121   CONTINUE
C
C New one, so add to list
      MRSYM=MRSYM+1
      DO 125,I=1,3
      DO 125,J=1,3
125   JRSYM(I,J,MRSYM)=JRMAT(I,J)
C
110   CONTINUE
100   CONTINUE
      RETURN
      END
C
C
C
      SUBROUTINE SETSYM(JRSYM,M)
C     ==========================
C
C Set symmetry matrix JRSYM to symmetry type M
C
C 1.  180-omega, 180+phi      self-rotation
C
C 2.  180-omega,  -phi        x dyads in both crystals
C
C 3.  180-omega, 180-phi      y dyads in both crystals
C
C 4.  omega,     180+phi      z dyads in both crystals
C
      DIMENSION JRSYM(3,3)
C
      IF(M.LE.0) RETURN
C
      DO 1,I=1,3
      DO 1,J=1,3
1     JRSYM(I,J)=0
      JRSYM(3,3)=1
C
      IF(M.EQ.1) THEN
C 180-omega, 180+phi
            JRSYM(1,1)=-1
            JRSYM(1,3)=180
            JRSYM(2,2)=+1
            JRSYM(2,3)=180
      ELSEIF(M.EQ.2) THEN
C 180-omega, -phi
            JRSYM(1,1)=-1
            JRSYM(1,3)=180
            JRSYM(2,2)=-1
            JRSYM(2,3)=0
      ELSEIF(M.EQ.3) THEN
C 180-omega, 180-phi
            JRSYM(1,1)=-1
            JRSYM(1,3)=180
            JRSYM(2,2)=-1
            JRSYM(2,3)=180
      ELSEIF(M.EQ.4) THEN
C omega, 180+phi
            JRSYM(1,1)=1
            JRSYM(1,3)=0
            JRSYM(2,2)=1
            JRSYM(2,3)=180
      ENDIF
C
      RETURN
      END
C
C
C
      SUBROUTINE SYMPUT(JA,JB,IA,IB)
C     ==============================
C
C Put point IA, IB (phi, omega values in degrees) into calculated
C area, using symmetry if required, return array index JA, JB.
C If no corresponding point, return JA=0
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
C INSIDE is a statement function to decide if a point is inside the
C area calculated
      LOGICAL INSIDE
      INSIDE(LA,LB)=(LA.GE.IPHI1.AND.LA.LE.IPHI2).AND.
     .              (LB.GE.IOMEG1.AND.LB.LE.IOMEG2)
C
C
      DO 10,L=0,MRSYM
C
      IF(L.EQ.0) THEN
C Identity
            IAA=IA
            IBB=IB
      ELSE
            IAA=JRSYM(2,2,L)*IA+JRSYM(2,3,L)
            IBB=JRSYM(1,1,L)*IB+JRSYM(1,3,L)
      ENDIF
C
      IF(INSIDE(IAA,IBB)) GO TO 100
C Check 'translations' on phi
      IF(IAA.GT.0) THEN
            IAA=IAA-360
      ELSE
            IAA=IAA+360
      ENDIF
      IF(INSIDE(IAA,IBB)) GO TO 100
10    CONTINUE
C
C Not found
      JA=0
      RETURN
C
C Valid point found, calculate array indices
100   JA=(IAA-IPHI1)/NPHI+1
      JB=(IBB-IOMEG1)/NOMEG+1
      RETURN
      END
C
C
C
      SUBROUTINE HPCNTR (A,M,N,BITS,CONT,NCONT,X,Y,IPTMAX)
C
C***************************************
C
C                  CONTUR ROUTINE
C
C***************************************
C
C  A IS FUNCTION TO BE CONTOURED, DIMENSIONED A(M,N) OR A(M*N)
C  M = NUMBER OF ROWS OF A
C  N = NUMBER OF COLUMN
C  BITS IS LOGICAL ARRAY DIMENSIONED AT LEAST 2*M*N IN CALLING ROUTINE.
C  CONT  ARRAY OF NCONT CONTOUR LEVELS
C  X,Y  ARRAYS OF DIMENSION IPTMAX USED IN SUBROUTINE TO STORE CONTOUR LINES
C
C  IDIV = NUMBER OF LINE SEGMENTS BETWEEN CONTOUR LEVEL CROSSINGS IN UNIT GRID
C     CELL.  IF >= 1, THEN ROUTINE FINDS IDIV-1 COORDINATES BETWEEN CROSSINGS,
C     BY LINEAR INTERPOLATION ON THE UNIT CELL.
C  CONTUR OUTPUT VIA CALLS TO SUBROUTINE APLOT (X,Y,IPT)
C     X,Y ARE COORDINATES WITH 0<=X<M AND 0<=Y<N.
C  SUBROUTINE CONSTL(CL) IS CALLED TO SET UP PLOTTING STYLE FOR CONTOUR LEVEL CL
C
      DIMENSION IDIR(4), ISIDE(5), IVCT(4), JVCT(4), KVCT(4), LVCT(4),
     *   SIDE(4), NVCT(4)
      LOGICAL BITS
      DIMENSION A(*),BITS(*),C9(4),CONT(NCONT)
      DIMENSION X(*),Y(*)
      DATA  IDIR/3,4,1,2/, ISIDE/1,2,3,4,1/, IVCT/0,-1,0,1/,
     * JVCT/-1,0,1,0/, SIDE/0.,0.,1.0,1.0/, KVCT/4,1,2,3/, LVCT/2,3,4,1/
C
      IDIV=1
      MN=M*N
      MN2=2*MN
C      WRITE(6,1008)M,N,MN,MN2
C 1008 FORMAT (/' M,N,MN,MN2 :',4I5)
      DIV=IDIV
      ICONT=1
C
C LOOP CONTOUR LEVELS CL
1     CL=CONT(ICONT)
      ICONT=ICONT+1
C
C
C
        DO 10 I=1,MN2
10    BITS(I)=.FALSE.
      NVCT(1)=0
      NVCT(2)=MN
      NVCT(3)=M
      NVCT(4)=MN+1
      IPT=1
      MM=M-1
      NN=N-1
C     SEARCH FOR CONTOUR CROSSING BETWEEN ADJACENT COLUMN OF ARRAY A(I,J)
      I=0
      J=1
      ISUB=0
      JSUB=0
      IRTN=1
  100 IF (J .GT. N) GO TO 140
  110 IF (I .GE. MM) GO TO 130
      I=I+1
      ISUB=ISUB+1
      JSUB=JSUB+1
      IF (A(ISUB)-CL) 115,600,120
115   IF (A(ISUB+1)-CL) 110,110,125
120   IF (A(ISUB+1)-CL) 125,110,110
125   IF (BITS(JSUB+NVCT(1)))  GO TO 110
      XSTART=(CL-A(ISUB))/(A(ISUB+1)-A(ISUB))
      YSTART=0
      GO TO 200
  130 I=0
      ISUB=ISUB+1
      JSUB=JSUB+1
      J=J+1
      GO TO 100
C     SEARCH FOR CONTOUR CROSSING BETWEEN ADJACENT ROWS OF ARRAY A(I,J)
  140 I=0
      J=1
      JSUB=0
      ISUB=0
      IRTN=2
  150 IF (J .GT. NN) GO TO 190
  160 IF (I .GE. M) GO TO 180
      I=I+1
      ISUB=ISUB+1
      JSUB=JSUB+1
      IF (A(ISUB)-CL) 165,160,170
165   IF (A(ISUB+M)-CL) 160,160,175
170   IF (A(ISUB+M)-CL) 175,160,160
175   IF (BITS(JSUB+NVCT(2)))  GO TO 160
      YSTART=(CL-A(ISUB))/(A(ISUB+M)-A(ISUB))
      XSTART=0
      GO TO 200
  180 I=0
      J=J+1
      GO TO 150
190   IF(ICONT.GT.NCONT)RETURN
      GO TO 1
C
C
C     BEGIN FOLLOWING CONTOUR LINE... SAVE INDICIES FOR RETURN TO SEARCH
  200 ISAVE=I
      JSAVE=J
      ISUBSV=ISUB
      JSUBSV=JSUB
      XSAVE=XSTART
      YSAVE=YSTART
      X(1)=XSTART+FLOAT(I-1)
      Y(1)=YSTART+FLOAT(J-1)
      IENT=IRTN
      IRS=0
      GO TO 250
C     DUMP LINE AND FOLLOW CONTOUR LINE ON OPPOSITE SIDE OF STARTING PIONT
C     WHEN USED A SECOND TIME THIS ENTRY RETURNS TO SEARCH
  205 IRS=1
  210 IF (IPT .GT. 1) CALL  APLOT (X,Y,IPT)
      IPT=1
      I=ISAVE
      J=JSAVE
      ISUB=ISUBSV
      JSUB=JSUBSV
      XSTART=XSAVE
      YSTART=YSAVE
      X(1)=XSTART+FLOAT(I-1)
      Y(1)=YSTART+FLOAT(J-1)
      IF (IRS.NE.0) GO TO (110,160), IRTN
      IEXIT=IRTN
      IRS=1
      GO TO 240
C     RETURN FROM FOLLOWING CONTOUR LINE THROUGH A CELL
230   IF (BITS(JSUB+NVCT(IEXIT)))  GO TO 205
  240 I=I+IVCT(IEXIT)
      J=J+JVCT(IEXIT)
      JSUB=I+(J-1)*M
      ISUB=JSUB
      IENT=IDIR(IEXIT)
250   BITS(JSUB+NVCT(IENT))=.TRUE.
      IF (I.LT.1 .OR. I.GT.MM .OR. J.LT.1 .OR. J.GT.NN)  GO TO 210
C     FIND CONTOUR CROSSING IN NEW CELL
      IF (ISUB+1.GT.MN .OR. ISUB+M.GT.MN
     1     .OR. ISUB+1+M.GT.MN)  GO TO 210
      C9(1)=A(ISUB+1)
      C9(2)=A(ISUB)
      C9(3)=A(ISUB+M)
      C9(4)=A(ISUB+1+M)
      JRTN=1
      ICNT=1
      JCNT=1
      DO 290 IROUND=1,4
      IF (IROUND .EQ. IENT) GO TO 290
      I1=ISIDE(IROUND)
      I2=ISIDE(IROUND+1)
      IF (C9(I1)-CL) 270,285,275
  270 IF (C9(I2)-CL) 290,290,280
  275 IF (C9(I2)-CL) 280,290,290
  280 IEXIT=IROUND
      ICNT=ICNT+1
      GO TO 290
  285 JEXIT=IROUND
      JCNT=JCNT+1
  290 CONTINUE
      GO TO (300,310,700,210), JCNT
  300 GO TO (210,320,210,800), ICNT
  310 GO TO (710,320,210,210), ICNT
  320 GO TO (330,340,350,360), IENT
  330 GO TO (210,410,500,410), IEXIT
  340 GO TO (510,210,510,400), IEXIT
  350 GO TO (500,410,210,410), IEXIT
  360 GO TO (510,400,510,210), IEXIT
C     FOLLOW CONTOUR LINE ACROSS A CELL TO A SIDE
  400 XSTART=SIDE(IENT)
  410 XFIN=SIDE(IEXIT)
      XINC=(XFIN-XSTART)/DIV
      XBASE=FLOAT(I-1)
      YBASE=FLOAT(J-1)
      A1=CL-C9(2)
      A2=C9(1)-C9(2)
      A3=C9(3)-C9(2)
      A4=C9(2)-C9(1)+C9(4)-C9(3)
      DO 440 INTERP=1,IDIV
      XSTART=XSTART+XINC
      YSTART=(A1-A2*XSTART)/(A3+A4*XSTART)
      IF (IPT.LT.IPTMAX)  GO TO 430
      CALL APLOT(X, Y, IPT)
      X(1)=X(IPT)
      Y(1)=Y(IPT)
      IPT=1
  430 IPT=IPT+1
      X(IPT)=XBASE+XSTART
      Y(IPT)=YBASE+YSTART
  440 CONTINUE
      GO TO (230,210,615,635), JRTN
  500 YSTART=SIDE(IENT)
C     FOLLOW CONTOUR LINE ACROSS A CELL TO A TOP OR BOTTOM
  510 YFIN=SIDE(IEXIT)
      XBASE=FLOAT(I-1)
      YINC=(YFIN-YSTART)/DIV
      YBASE=FLOAT(J-1)
      A1=CL-C9(2)
      A2=C9(3)-C9(2)
      A3=C9(1)-C9(2)
      A4=C9(2)-C9(1)+C9(4)-C9(3)
      DO 540 INTERP=1,IDIV
      YSTART=YSTART+YINC
      XSTART=(A1-A2*YSTART)/(A3+A4*YSTART)
      IF (IPT.LT.IPTMAX) GO TO 530
      CALL APLOT(X, Y, IPT)
      X(1)=X(IPT)
      Y(1)=Y(IPT)
      IPT=1
  530 IPT=IPT+1
      X(IPT)=XBASE+XSTART
      Y(IPT)=YBASE+YSTART
  540 CONTINUE
      GO TO (230,210,615,635), JRTN
C     FOLLOW CONTOUR LINE FROM CORNER TO CORNER
600   K1=ISUB-M
      K2=ISUB+1-M
      K3=ISUB+1
      K4=ISUB+1+M
      K5=ISUB+M
      K6=ISUB-1+M
      K7=ISUB-1
      C1=A(K1)
      C2=A(K2)
      C3=A(K3)
      C4=A(K4)
      C5=A(K5)
      C6=A(K6)
      C7=A(K7)
      C8=A(ISUB)
      IF (ISUB.LT.1 .OR. ISUB.GT.MN) GO TO 640
      X(1)=FLOAT(I-1)
      Y(1)=FLOAT(J-1)
      IF (J .EQ. 1 .OR. J .EQ. M)  GO TO 610
      IF (K1.LT.1 .OR. K1.GT.MN)  GO TO 610
      IF (K2.LT.1 .OR. K2.GT.MN)  GO TO 610
      IF (K3.LT.1 .OR. K3.GT.MN)  GO TO 610
      IF (K4.LT.1 .OR. K4.GT.MN)  GO TO 610
      IF (K5.LT.1 .OR. K5.GT.MN)  GO TO 610
      IF (C3.NE.CL)  GO TO 610
      IF (C1 .EQ. CL .AND. C2 .EQ. CL .AND.
     *    C4 .EQ. CL .AND. C5 .EQ. CL)  GO TO 610
      X(2)=X(1)+1.
      Y(2)=Y(1)
      CALL APLOT (X, Y, 2)
      GO TO 620
  610 IF (J .EQ. 1)  GO TO 620
      IF (K1.LT.1 .OR. K1.GT.MN)  GO TO 620
      IF (K2.LT.1 .OR. K2.GT.MN)  GO TO 620
      IF (K3.LT.1 .OR. K3.GT.MN)  GO TO 620
      IF (C2 .NE. CL)  GO TO 620
      IF (C1 .EQ. CL .OR. C3 .EQ. CL)  GO TO 620
      IF (C1 .GT. CL .AND. C3 .GT. CL .OR.
     *    C1 .LT. CL .AND. C3 .LT. CL)  GO TO 620
      C9(1)=C2
      C9(2)=C1
      C9(3)=C8
      C9(4)=C3
      J=J-1
      JRTN=3
      IENT=3
      IEXIT=1
      GO TO 500
  615 IF (IPT .GT. 1)  CALL APLOT (X, Y, IPT)
      IPT=1
      J=J+1
      X(1)=FLOAT(I-1)
      Y(1)=FLOAT(J-1)
  620  IF (J .EQ. M .OR. I .EQ. 1)  GO TO 630
      IF (K3.LT.1 .OR. K3.GT.MN)  GO TO 630
      IF (K4.LT.1 .OR. K4.GT.MN)  GO TO 630
      IF (K5.LT.1 .OR. K5.GT.MN)  GO TO 630
      IF (K6.LT.1 .OR. K6.GT.MN)  GO TO 630
      IF (K7.LT.1 .OR. K7.GT.MN)  GO TO 630
      IF (C5 .NE. CL)  GO TO 630
      IF (C3 .EQ. CL .AND. C4 .EQ. CL .AND.
     *    C6 .EQ. CL .AND. C7 .EQ. CL)   GO TO 630
      X(2)=X(1)
      Y(2)=Y(1)+1.
      CALL APLOT (X, Y, 2)
      GO TO 640
  630 IF (J .EQ. M)  GO TO 640
      IF (K3.LT.1 .OR. K3.GT.MN)  GO TO 640
      IF (K4.LT.1 .OR. K4.GT.MN)  GO TO 640
      IF (K5.LT.1 .OR. K5.GT.MN)  GO TO 640
      IF (C4 .NE. CL)  GO TO 640
      IF (C3 .EQ. CL .OR. C5 .EQ. CL)  GO TO 640
      IF (C3 .GT. CL .AND. C5 .GT. CL .OR.
     *    C3 .LT. CL .AND. C5 .LT. CL)  GO TO 640
      C9(1)=C3
      C9(2)=C8
      C9(3)=C5
      C9(4)=C4
      JRTN=4
      IENT=1
      IEXIT=3
      GO TO 500
  635 IF (IPT .GT. 1) CALL APLOT (X, Y, IPT)
      IPT=1
      X(1)=FLOAT(I-1)
      Y(1)=FLOAT(J-1)
  640 GO TO (110,160), IRTN
C    FOLLOW CONTOUR LINE FROM SIDE TO CORNER OR CORNERS
  700 JRTN=2
      IOPP=IDIR(IENT)
      I1=ISIDE(IOPP)
      I2=ISIDE(IOPP+1)
      IEXIT=IOPP
      C9(I1)=C9(KVCT(I1))
      C9(I2)=C9(LVCT(I2))
      GO TO 320
  710 JRTN=2
      IEXIT=JEXIT
      GO TO 320
C     FOLLOW CONTOUR LINE THROUGH SADDLE POINT
  800 IOPP=IDIR(IENT)
      I1=ISIDE(IENT)
      C1=C9(I1)
      I2=ISIDE(IENT+1)
      C2=C9(I2)
      I3=ISIDE(IOPP)
      C3=C9(I3)
      I4=ISIDE(IOPP+1)
      C4=C9(I4)
      IF ((C1-CL)/(C1-C2) .EQ. (C4-CL)/(C4-C3))  GO TO 820
      IF ((C1-CL)/(C1-C4) .GT. (C2-CL)/(C2-C3))  GO TO 810
      IEXIT=I4
      GO TO 320
  810 IEXIT=I2
      GO TO 320
  820 C9(I3)=C9(I2)
      C9(I4)=C9(I1)
      IEXIT=I3
      GO TO 320
      END
C
      SUBROUTINE APLOT(X,Y,IPT)
C     =========================
C
C JOIN IPT COORDINATE PAIRS X,Y OF CONTOUR CROSSINGS, IN GRID UNITS
C FOR M*N GRID , COORDINATES ARE IN RANGE 0,0 TO M,N
C
      DIMENSION X(IPT),Y(IPT)
C
C SAVE 1ST AND LAST POINTS BEFORE TRANSFORMATION
C - - - - - -- - - - - - - - - - - - - - - - - - - - - -
C
C RESET ALL SORTS OF TRANSFORMATIONS AND SCALE ETC TO USE THIS
C SUBROUTINE BLINDLY IN RFPLOT
C
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
      X1=X(1)
      Y1=Y(1)
      XL=X(IPT)
      YL=Y(IPT)
C
C DRAW LINES
C
      ILPPP = 10
      IF (ILPPP.GT.IPT) ILPPP = IPT
C     WRITE(6,1000)(X(INNN),Y(INNN),INNN=1,ILPPP)
C 1000 FORMAT(// ' MADE IT INTO APLOT',15F5.1)
      CALL LINES(X,Y,IPT)
C RESTORE 1ST & LAST POINTS
      X(1)=X1
      Y(1)=Y1
      X(IPT)=XL
      Y(IPT)=YL
      RETURN
C
      END
C
      SUBROUTINE LINES(X,Y,NPT)
C     =========================
C
C DRAW LINES TO JOIN NPT POINTS X,Y
C
C Note y = phi, x = omega, in grid units
C
      DIMENSION X(NPT),Y(NPT)
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
      CONV=360.0/6.28318
      IF(NPT.LE.1) RETURN
C
C  STEREOGRAPH
C
      DO 13 ILLL=1,NPT
C
      OMEGA=X(ILLL)*FLOAT(NOMEG)
      IF(LHALF.EQ.2) OMEGA=90.-OMEGA
      OMEGA=OMEGA/CONV
      R=SIN(OMEGA)/(COS(OMEGA)+1.0)
C
      PHI=Y(ILLL)*FLOAT(NPHI)
      IF(LHALF.EQ.2.AND.LSYMFG.LT.0) PHI=180.+PHI
      PHI=PHI/CONV
C
      X(ILLL)=R*COS(PHI)
13    Y(ILLL)=R*SIN(PHI)
C
      CALL PLTLIN(X,Y,NPT)
      RETURN
      END
C
C
C
      SUBROUTINE PEAKPK(RROT,LA,LB,IG,RMAP)
C     =====================================
C
C Find peaks in previous section
C   Threshold is PEAK if LPKRMS=0, PEAK*RMS if LPKRMS=1
C
C IG is kappa (kappa) for current section, IG=-1 for search of final section
C
C RMAP is used as a rotating buffer, RMAP(,,I1), RMAP(,,I2), RMAP(,,I3),
C      contain 3 sections, I2 is the one to be searched
C
C
      DIMENSION RROT(LA,LB),RMAP(0:LA+1,0:LB+1,3)
C
C
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
      PARAMETER (MAXNPK=100)
      COMMON /PEAKS/ NPEAK,PEAKS(MAXNPK),PKPOS(3,MAXNPK),LPEAK,
     .  LPKORD(MAXNPK),LPKOUT
C
      REAL OFFSET(3),OFFMAX,PEAKMX(-1:+1,-1:+1,-1:+1)
      REAL RMAX
      SAVE I1,I2,I3,KKAPA,OFFMAX
C
      DATA I2/0/
C Largest acceptable offset from central grid point
      DATA OFFMAX / 0.6 /
C
      IF (MAXPEK .EQ. 0) RETURN
C Set threshold
      THRHLD=PEAK
      IF (LPKRMS.GT.0) THRHLD=THRHLD*RRMS
C
      IF(IG.LT.0) THEN
C  If I1 = -1, this is first & last section (1 section only), so
C   generate dummy section all = -100.
            IF (I1 .LT. 0) THEN
                I1 = 1
                DO 1, J = 0, LB+1
                DO 1, I = 0, LA+1
1       RMAP(I,J,I1) = -100.
            ENDIF
C For final section, copy slot 1 to slot 3
                DO 2, J = 0, LB+1
                DO 2, I = 0, LA+1
2       RMAP(I,J,I3) =  RMAP(I,J,I1)
CEJD            CALL MOVEI(RMAP(0,0,I3), RMAP(0,0,I1), (LA+2)*(LB+2))
            GO TO 10
      ENDIF
C
      IF(I2.EQ.0) THEN
C First section, just store section
            I1=-1
            I2=2
            I3=3
C KKAPA is kappa (kappa) for section next to be searched (slot I2)
            KKAPA=IG
            DO 5, J = 1, LB
            DO 5, I = 1, LA
             RMAP(I,J,I2)=RROT(I,J)
CEJD               CALL MOVEI(RMAP(1,J,I2),RROT(1,J),LA)
5           CONTINUE
C  To handle edge effects, each section is extended by copying the
C  2nd & penultimate rows and columns to give a extra row & column
C  all round. This is equivalent to assuming each edge is a mirror
C  plane, which is unlikely to be true, however it will force all
C  edge peaks to lie on the edge plane. Similarly, sandwich the 1st
C  section, between 2 copies of the 2nd section, and the last section
C  between 2 copies of the penultimate
C
            CALL EXSECT(RMAP(0,0,I2), LA, LB)
C
C  Set up interpolation matrix
            CALL PKSETM
C
            RETURN
      ENDIF
C
C Store new section in slot I3
      DO 6, J = 1, LB
            DO 6, I = 1, LA
             RMAP(I,J,I3)=RROT(I,J)
CEJD         CALL MOVEI(RMAP(1,J,I3),RROT(1,J),LA)
6     CONTINUE
      CALL EXSECT(RMAP(0,0,I3), LA, LB)
C
      IF (I1 .LT. 0) THEN
C  2nd section, copy into slot I1
            I1 = 1
                DO 12, J = 0, LB+1
                DO 12, I = 0, LA+1
12      RMAP(I,J,I1) =  RMAP(I,J,I3)
CEJD            CALL MOVEI(RMAP(0,0,I1), RMAP(0,0,I3), (LA+2)*(LB+2))
      ENDIF
C
C Now we have a section KKAPA in slot I2, and maybe in slots I1 and I3 if slot
C numbers are .gt. 0, so search section I2
C First check for special cases kappa=0
10    IF(KKAPA.EQ.0.OR.KKAPA.EQ.360) THEN
            NA1=1
            NA2=1
            NB1=1
            NB2=1
      ELSE
            NA1=1
            NA2=LA
            NB1=1
            NB2=LB
      ENDIF
C
C Search for peaks above threshold THRHLD
      DO 20,I=NA1,NA2
      DO 21,J=NB1,NB2
      R=RMAP(I,J,I2)
C Is this point above THRHLD?
      IF(R.GT.THRHLD) THEN
C Gather 3 x 3 x 3 grid centred on this point
            RMAX = -100.
            DO 25, L2 = -1, +1
            DO 25, L1 = -1, +1
            PEAKMX(L1,L2,-1) = RMAP(I+L1,J+L2,I1)
            PEAKMX(L1,L2, 0) = RMAP(I+L1,J+L2,I2)
            PEAKMX(L1,L2,+1) = RMAP(I+L1,J+L2,I3)
            DO 26, L3 = -1, +1
            RMAX = MAX (RMAX, PEAKMX(L1,L2,L3))
26        CONTINUE
25        CONTINUE
C
C     Is the central point the highest?
            IF ( RMAX .GT. PEAKMX(0,0,0)) GO TO 21
C  No neighbours higher, so we have a peak
            N=NPEAK+1
            IF(N.GT.MAXPEK) THEN
C  We've already got enough peaks, so find smallest if smaller than the
C  new one
                  A=R
                  N=0
                  DO 22,K=1,NPEAK
                  IF(PEAKS(K).LT.A) THEN
                        A=PEAKS(K)
                        N=K
                  ENDIF
22                CONTINUE
                  IF(N.EQ.0) GO TO 21
            ELSE
C  Append to list
                  NPEAK=N
            ENDIF
C
C Store peak in slot N
            PEAKS(N)=R
C Interpolate position
            CALL PKFITP(PEAKMX, OFFMAX, OFFSET, ISTAT)
C
            PKPOS(1,N)=(J-1+OFFSET(2))*NOMEG+IOMEG1
            PKPOS(2,N)=(I-1+OFFSET(1))*NPHI+IPHI1
            PKPOS(3,N)=KKAPA+OFFSET(3)*NKAPA
      ENDIF
C
21    CONTINUE
C
C If omeg limits are 0 or 180, reset omeg limits after first row
C This removes some redundancy, but is by no means foolproof.
C More complicated tests are needed
      IF(I.EQ.1) THEN
            IF(IOMEG1.EQ.0) NB1=NB1+1
            IF(IOMEG2.EQ.180) NB2=NB2-1
      ENDIF
C
20    CONTINUE
C
C Set pointers for next section
      I1=I2
      I2=I3
      I3=I2+1
      IF(I3.GT.3) I3=1
      KKAPA=IG
C
C
      RETURN
C
      END
C
C
C
C
C
      SUBROUTINE EXSECT(RMAP, LA, LB)
C     ===============================
C
C Extend section by 1 row & column all round, assuming mirror plane
C on edge. This is likely not to be true, but at least will force
C edge peaks to lie on the edge
C
      INTEGER LA, LB
      REAL RMAP(0:LA+1,0:LB+1)
C
      INTEGER I
C
      DO 10, I = 1, LA
      RMAP(I, 0) = RMAP(I, 2)
10    RMAP(I, LB+1) = RMAP(I, LB-1)
C
      DO 20, I = 0, LB+1
      RMAP(0, I) = RMAP(2, I)
20    RMAP(LA+1, I) = RMAP(LA-1, I)
C
      RETURN
      END
C
C
C
      SUBROUTINE PKSETM
C     =================
C
C Set up normal matrix etc for 3D paraboloid fit
C
C Observational equation is
C  r(calc) = a x**2 + b y**2 + c z**2 + d xy + e yz + f zx +
C            g x + h y + i z + j
C
C Observation vector
C  xT = ( x**2  y**2  z**2  xy  yz  zx   x   y   z   1 )
C
C Observation matrix
C
C  A  =  (  x(1)T  )          for each point on 3 x 3 x 3 grid
C        (  x(2)T  )
C        ( . . . . )
C        (  x(n)T  )
C
      INTEGER NP, NE, NG
      PARAMETER ( NP = 10, NE = 3, NG = NE*NE*NE)
      REAL  AT(NP, NG), ATA(NP, NP)
C
      COMMON /PKFIT/ SMAT(NP, NG), A(NG, NP)
      REAL SMAT, A
C
      INTEGER IX, IY, IZ, I, J
      REAL X, Y, Z, DET
      INTEGER W1(NP), W2(NP)
C
      SAVE /PKFIT/
C
C  Loop points on 3 x 3 x 3 grid
      I = 0
      DO 10, IZ = -1, +1
      Z = IZ
      DO 10, IY = -1, +1
      Y = IY
      DO 10, IX = -1, +1
      X = IX
      I = I + 1
C Matrix elements
      A(I, 1) = X**2
      A(I, 2) = Y**2
      A(I, 3) = Z**2
      A(I, 4) = X * Y
      A(I, 5) = Y * Z
      A(I, 6) = Z * X
      A(I, 7) = X
      A(I, 8) = Y
      A(I, 9) = Z
      A(I, 10) = 1.0
10    CONTINUE
C
C Form A(transpose)
      DO 20, I = 1, NG
      DO 20, J = 1, NP
20    AT(J, I) = A(I, J)
C
C A(tranpose) * A
      CALL GMPRD(AT, A, ATA, NP, NG, NP)
C
C
C Invert in place
      CALL MINVN(ATA, NP, DET, W1, W2)
      IF (DET .EQ. 0.0) THEN
        CALL CCPERR(1, 'PKSETM: singular matrix')
      ENDIF
C
C  SMAT = (AT A)**1 AT
      CALL GMPRD(ATA, AT, SMAT, NP, NP, NG)
C
      RETURN
      END
C
C
C
      SUBROUTINE PKFITP (PEAKMX, OFFMAX, OFFSET, ISTAT)
C     =================================================
C
C Fit paraboloid to 3 x 3 x 3 grid in PEAKMX (central point assumed to
C be maximum), return position of maximum in OFFSET (grid units)
C Subroutine PKSETM must have been called first to calculate solution
C matrix SMAT
C
C On entry:
C  PEAKMX       values of peak  (3 x 3 x 3 grid)
C  OFFMAX       maximum value of offset
C
C On exit:
C  OFFSET(3)    offsets in grid units
C  ISTAT        = 1 if offset reset to maximum, else = 0
C
C
      REAL PEAKMX(-1:+1,-1:+1,-1:+1), OFFMAX, OFFSET(3)
      INTEGER ISTAT
C
      INTEGER NP, NE, NG
      PARAMETER ( NP = 10, NE = 3, NG = NE*NE*NE)
C
      COMMON /PKFIT/ SMAT(NP, NG), A(NG, NP)
      REAL SMAT, A
C
      INTEGER I
      REAL P(NP), B(3,3), G(3), DET
      INTEGER W1(3), W2(3)
C
      SAVE /PKFIT/
C
C
C Solve for parameters P(1)->P(10), = a,b,c,d,e,f,g,h,i,j
C
C    p  = [ (AT A)**-1 AT ]  r
C where  (AT A)**-1 AT = SMAT is calculated in PKSETM
C        r = PEAKMX, observation vector
C Note that in PKSETM, the order of grid points corresponding to the rows
C of matrix A must correspond to the order of values in PEAKMX, ie x fast,
C y medium, z slow
C
      CALL GMPRD (SMAT, PEAKMX, P, NP, NG, 1)
C
C Set up matrix B and vector g for solving position of maximum from the
C paraboloid parameters
C
C           (  2a   d   f  )                    (  -g  )
C     B  =  (   d  2b   e  )              g  =  (  -h  )
C           (   f   e  2c  )                    (  -i  )
C
C   Then position of maximum  OFFSET  =  x  =  B**-1 g
C
      B(1,1) = 2. * P(1)
      B(1,2) = P(4)
      B(2,1) = P(4)
      B(1,3) = P(6)
      B(3,1) = P(6)
      B(2,2) = 2. * P(2)
      B(2,3) = P(5)
      B(3,2) = P(5)
      B(3,3) = 2. * P(3)
C
      G(1) = -P(7)
      G(2) = -P(8)
      G(3) = -P(9)
C
C Invert B
      CALL MINVN(B, 3, DET, W1, W2)
      IF (DET .EQ. 0.0) THEN
        CALL CCPERR(1, 'PKFITP: singular matrix')
      ENDIF
C
      CALL GMPRD(B, G, OFFSET, 3, 3, 1)
C
      ISTAT = 0
      DO 50, I = 1, 3
      IF (ABS(OFFSET(I)) .GT. OFFMAX) THEN
C Reset to maximum, set flag
            ISTAT = 1
            OFFSET(I) = SIGN(OFFMAX, OFFSET(I))
      ENDIF
50    CONTINUE
C
      RETURN
      END
C
      SUBROUTINE SORT1(N,A0,A1)
C     =========================
C
C****
C**** N IS THE LENGTH OF THE ARRAY TO BE SORTED I.E. THE NUMBER OF
C**** ITEMS TO BE SORTED
C**** THE ELEMENTS OF ARRAY A0(N) ARE SORTED INTO INCREASING ORDER OF
C**** MAGNITUDE.  ARRAYS A1 MAY BE CARRIED WITH A0 AS REQUIRED
C****
      DIMENSION A0(N),A1(N)
      INTEGER A1
C
      NO2 = N/2
      DO 100 I=1,NO2
      C = A0(I)
      D = C
      IH = I
      K = I
      IP1 = I+1
      NP1 = N+1-I
      DO 20 J=IP1,NP1
      A = A0(J)
      IF ( A.GE.C ) GO TO 15
      IH = J
      C = A
      GO TO 20
   15 IF ( D.GE.A ) GO TO 20
      D = A
      K = J
   20 CONTINUE
      IF ( C.EQ.D ) GO TO 200
      IF ( IH.EQ.I ) GO TO 30
      C = A1(IH)
      A1(IH) = A1(I)
      A1(I) = C
      C = A0(IH)
      A0(IH) = A0(I)
      A0(I) = C
   30 IF ( K.EQ.I ) K=IH
      IF ( K.EQ.NP1 ) GO TO 100
      D = A1(K)
      A1(K) = A1(NP1)
      A1(NP1) = D
      D = A0(K)
      A0(K) = A0(NP1)
      A0(NP1) = D
  100 CONTINUE
  200 RETURN
      END
C
C
C
      SUBROUTINE PPEAKS
C     =================
C
C Print all peaks found
C
      COMMON /MATRIX/ ORTHOG(3,3,2),FRAC(3,3,2),SYM(4,4,192,2),NSYM(2),
     . MRSYM,JRSYM(3,3,12),LSYMFG,LHALF
      COMMON /ALMFRF/ NCRYST,IPHI1,IPHI2,NPHI,IOMEG1,IOMEG2,NOMEG,
     . IKAPA1,IKAPA2,NKAPA,ARAD,EPS,RESMIN,RESMAX,NCODES(2),LSPGPS(2),
     . TFACS(2),FRMINS(2),FRMAXS(2),LPRINT,LMAP,LPLOT,CONTRS(50),
     . NCONTR,LREAD,LSAVE,PEAK,MAXPEK,LPKRMS,KRKAP1,KRKAP2,
     . RMIN,RMEAN,RRMS,RWT,NPOINT
C
      COMMON /PKFILE/ PKFNAM
      CHARACTER*80 PKFNAM
      SAVE /PKFILE/
C
      CHARACTER*80 TITLE,FNAME
      COMMON /TEXTC/ TITLE,FNAME
C
      PARAMETER (MAXNPK=100)
      COMMON /PEAKS/ NPEAK,PEAKS(MAXNPK),PKPOS(3,MAXNPK),LPEAK,
     .  LPKORD(MAXNPK),LPKOUT
C
      DIMENSION AMINV(4,4)
      DIMENSION ARSYM(3,3,24,2),AJUNK4(3,3),TRSYM(3,3)
     1   ,TRMAT(3,3),DC(5),AJUNK2(3,3),AJUNK3(3,3),AJUNK1(3,3)
      REAL EPKPOS(3,MAXNPK),PKHGHT(MAXNPK)
C
      CHARACTER*15  CODES(5)
      EXTERNAL MATMUL
      DATA CODES/'a,c*xa,c*','b,a*xb,a*','c,b*xc,b*',
     . 'a+b,c*x(a+b),c*','a*,cxa*,c'/
C
      IF (NPEAK .LE. 0) RETURN
C
      IF (LPEAK .LE. 0) THEN
C Sort reflection pointers into ascending order of magnitude
C   Initialize pointers and copy peak height. Note that the sort
C   does not reorder storage arrays PEAKS & PKPOS, but just creates
C   a list of pointers
         DO 200, I = 1, NPEAK
            PKHGHT(I) = PEAKS(I)
200      LPKORD(I) = I
         CALL SORT1(NPEAK, PKHGHT, LPKORD)
C
      ELSE
C  PEAK option, set peak order LPKORD to reverse, so that peaks are
C    printed out in same order as input
         DO 210, I = 1, NPEAK
            LPKORD(I) = NPEAK - I + 1
210      CONTINUE
      ENDIF
C
      LI=6
C
C Print symmetry
      CALL PRTSYM(LI,1,NSYM(1),SYM(1,1,1,1))
      IF (NCRYST.EQ.2) CALL PRTSYM(LI,2,NSYM(2),SYM(1,1,1,2))
C
      IF(LPEAK.GT.0) THEN
         WRITE(LI,9011) NPEAK
9011     FORMAT('1     Positions of ',I4,' peaks'//)
      ELSEIF(LPKRMS.EQ.0) THEN
         WRITE(LI,9012) NPEAK,PEAK
9012     FORMAT('1     Positions of ',I4,' peaks above',F8.0//)
      ELSE
         WRITE(LI,9013) NPEAK,PEAK,RRMS
9013     FORMAT('1     Positions of ',I4,' peaks above',F8.0,
     .    ' * RMS = ',F9.3//)
      ENDIF
      WRITE(LI,9014)
9014  FORMAT(
     . 24X,' Eulerian angles                        Polar angles'//
     . 20X,'  Alpha   Beta   Gamma      Peak      ',
     . ' Omega    Phi  Kappa          Direction cosines'/
     . ' Symmetry: 1   2'//)
C
      JSYM1=NSYM(1)
      JSYM2=NSYM(2)
C  Set up matrices for crystal 1
      DO 1 ISYM1=1,JSYM1
      DO 3 I=1,3
      DO 3 J=1,3
3     AJUNK3(I,J)=SYM(I,J,ISYM1,1)
C
      CALL MATMUL(AJUNK4,AJUNK3,FRAC(1,1,1))
      CALL MATMUL(AJUNK3,ORTHOG(1,1,1),AJUNK4)
      DO 4 I=1,3
      DO 4 J=1,3
4     ARSYM(I,J,ISYM1,1)=AJUNK3(I,J)
C  ARSYM(,,1 to JSYM1) contains  Orthog . Sym . Frac
1     CONTINUE
C
C  Set up matrices for crystal 2
      DO 5 ISYM2=1,JSYM2
C    inverse symmetry matrix AMINV  (4x4)
      CALL RBRINV(SYM(1,1,ISYM2,2),AMINV)
C    copy to 3x3
      DO 66 I=1,3
      DO 66 J=1,3
      AJUNK3(J,I)=AMINV(J,I)
66    CONTINUE
C
      CALL MATMUL(AJUNK4,AJUNK3,FRAC(1,1,2))
      CALL MATMUL(AJUNK3,ORTHOG(1,1,2),AJUNK4)
      DO 7 I=1,3
      DO 7 J=1,3
7     ARSYM(I,J,ISYM2,2)=AJUNK3(I,J)
C  ARSYM(,,ISYM2,2) contains  Orthog . Sym**-1 . Frac
5     CONTINUE
C
C  Loop all peaks in list
C
      DO 100 JP=1,NPEAK
C  Get peak number in descending order of size
      J=LPKORD(NPEAK-JP+1)
C
      WRITE(LI,15) JP
15    FORMAT(/1X,' Peak ',I4)
C
C Get rotation matrix corresponding to angles in list
C
      CALL POLMAT(TRMAT,PKPOS(1,J))
C
C  Convert to Eulerian for output file of unique peaks
      CALL EULERA(TRMAT,ALPHA,BETA,GAMMA)
      CALL RNGANG(ALPHA)
      EPKPOS(1,JP)=ALPHA
      CALL RNGANG(BETA)
      EPKPOS(2,JP)=BETA
      CALL RNGANG(GAMMA)
      EPKPOS(3,JP)=GAMMA
C
C
      DO 99 ISYM1=1,JSYM1
      DO 103 I=1,3
      DO 103 JJ=1,3
103   AJUNK3(I,JJ)=ARSYM(I,JJ,ISYM1,1)
C
      DO 98 ISYM2=1,JSYM2
       DO 102 I=1,3
      DO 102 JJ=1,3
102   AJUNK2(I,JJ)=ARSYM(I,JJ,ISYM2,2)
      CALL MATMUL(AJUNK1,TRMAT,AJUNK2)
      CALL MATMUL(TRSYM,AJUNK3,AJUNK1)
      TRACE=TRSYM(1,1)+TRSYM(2,2)+TRSYM(3,3)
      IF (TRACE.GE.2.9999) THEN
         WRITE(LI,116) PEAKS(J),(PKPOS(I,J),I=1,3)
116      FORMAT('  Origin peak   ',28X,F8.1,5X,3F7.1)
         GO TO 100
      ENDIF
C  Get Euler angles
      CALL EULERA(TRSYM,ALPHA,BETA,GAMMA)
C  Get Kappa & direction cosines
       CALL DCOSFD(TRSYM, DC, AKAPPA)
C  Get polar angles
      CALL POLANG(DC,OMEGA,PHI)
C Put all angles in range 0-360 degrees
      CALL RNGANG(ALPHA)
      CALL RNGANG(BETA)
      CALL RNGANG(GAMMA)
      CALL RNGANG(OMEGA)
      CALL RNGANG(PHI)
      CALL RNGANG(AKAPPA)
      WRITE(LI,14) ISYM1,ISYM2,ALPHA,BETA,GAMMA,PEAKS(J),OMEGA,
     1  PHI,AKAPPA,DC(1),DC(2),DC(3)
14    FORMAT(8X,2I4,4X,3F7.1,3X,F8.1,5X,3F7.1,6X,3F8.4)
C
C
      IF(LPEAK.GT.1) THEN
C Output matrices if peaks input and MATRICES keyword given
        CALL MINV3(AJUNK2,TRSYM,D)
        WRITE(LI,114)
     1  ((TRSYM(I,JI),JI=1,3),(AJUNK2(I,K),K=1,3),I=1,3)
114     FORMAT(/'  Rotation matrix and its inverse',
     +    /(3(3F10.5,10X,3F10.5/))/)
      ENDIF
C
98    CONTINUE
99    CONTINUE
 100    CONTINUE
C
      WRITE(LI,9022)
9022  FORMAT(///' The rotation given by the angles of a peak rotates',
     . ' coordinates in the orthogonal frame'/
     . ' of crystal 2 to the orthogonal frame of crystal 1.'/
     . ' Beware of axis permutations introduced by NCODE = 2, 3 or 4'//
     . /' Rotation matrices are defined as follows:'///
     . ' ( l**2+(m**2+n**2)cos k     lm(1-cos k)-nsin k     ',
     .  'nl(1-cos k)+msin k )'/
     . ' ( lm(1-cos k)+nsin k      m**2+(l**2+n**2)cos k    ',
     .  'mn(1-cos k)-lsin k )',5X,'Polar angles'/
     . ' ( nl(1-cos k)-msin k        mn(1-cos k)+lsin k   ',
     .  'n*2+(l**2+m**2)cos k )'//
     . '     where l m n are the direction cosines of the axis about',
     . ' which the rotation k = kappa takes place.'//
     .   31X,'( l )             ( sin omega   cos phi  )'/
     .   31X,'( m )     =       ( sin omega   sin phi  )'/
     .   31X,'( n )             (       cos omega      )'/)
C
      WRITE(LI,9023)
9023  FORMAT(//
     . ' ( cosa cosb cosg - sina sing      -cosa cosb sing - sina cosg',
     .  '      cosa sinb )'/
     . ' ( sina cosb cosg + cosa sing      -sina cosb sing + cosa cosg',
     .  '      sina sinb )',5X,'Eulerian angles Alpha, Beta, Gamma'/
     . ' (         -sinb cosg                         sinb sing       ',
     .  '         cosb   )'/)
C
C
C------------------------------------------------------------------------
C Write unique peaks out to file
C
      IF (LPKOUT .GT. 0) THEN
         L=LENSTR(PKFNAM)
         WRITE(LI, 61001) PKFNAM(1:L)
61001    FORMAT(//' Unique peaks written out to file ',A)
         LUN=1
         CALL CCPDPN(LUN, PKFNAM(1:L),'NEW','F',0,0)
C
         L=LENSTR(TITLE)
         WRITE (LUN, 61000) TITLE(1:L)
61000    FORMAT('!   Peaks from POLARRFN'/'TITLE ',A)
         IF (NCRYST .EQ. 1) THEN
            WRITE (LUN, 61002) 'SELF'
         ELSEIF (NCRYST .EQ. 1) THEN
            WRITE (LUN, 61002) 'CROSS'
61002       FORMAT(A)
         ENDIF
         DO 1001, L=1,NCRYST
            WRITE (LUN, 61003) NCODES(L),CODES(NCODES(L))
61003       FORMAT(
     .        'CRYSTAL',I5,'       ! Ncode, orthogonalization ',A)
            WRITE (LUN, 61004) L, NSYM(L)
61004       FORMAT('SYMORTHOGONAL',2I5)
C
            DO 1002, K=1,NSYM(L)
               WRITE (LUN, 61005) K,((ARSYM(I,J,K,L), J=1,3), I=1,3)
61005          FORMAT('!   Symmetry operation ',I5/
     .                     2(3F12.5/),3F12.5)
1002        CONTINUE
1001      CONTINUE
C
C Now the peaks
C
          WRITE (LUN, 61006) NPEAK
61006     FORMAT('PEAKS',I8/
     .     '!         Alpha      Beta     Gamma     Height')
          DO 1010, L = 1, NPEAK
             K = LPKORD(NPEAK-L+1)
             WRITE (LUN, 61007) L, (EPKPOS(I, L), I = 1, 3), PEAKS(K)
61007        FORMAT(I6,3F10.1,F10.2)
1010      CONTINUE
          CLOSE(UNIT=LUN,STATUS='KEEP')
      ENDIF
C
      RETURN
      END
C
C
      SUBROUTINE PRTSYM(IOUT,ICRYST,NSYM,ROT)
C     =======================================
C
C Print symmetry matrices
C
      INTEGER IOUT,ICRYST,NSYM
      REAL ROT(4,4,NSYM)
C
      INTEGER J,M,L
C
      WRITE(IOUT,2) ICRYST
2     FORMAT(///' Unpermuted symmetry operations for crystal ',I3,': ')
C
      DO 3 J=1,NSYM
         WRITE(IOUT,4) J,((ROT(L,M,J),M=1,4),L=1,3)
4        FORMAT(/21X,'Symmetry matrix',I5,5X,3F10.5,5X,F10.5/
     .    2(46X,3F10.5,5X,F10.5/))
C
3     CONTINUE
      RETURN
      END
C
      SUBROUTINE EULERA(TRSYM,ALPHA,BETA,GAMMA)
C     =========================================
C
C Get Euler angles ALPHA,BETA,GAMMA from matrix TRSYM
C
      REAL TRSYM(3,3),ALPHA,BETA,GAMMA
C
      CONV=45./ATAN(1.)
      IF(ABS(TRSYM(3,3)).GE.0.99999) GO TO 13
C
      BETA=CONV*ACOS(TRSYM(3,3))
C
      R13=TRSYM(1,3)
      R23=TRSYM(2,3)
      R31=TRSYM(3,1)
      R32=TRSYM(3,2)
      GAMMA=CONV*ATAN2(R32,-R31)
      ALPHA=CONV*ATAN2(R23,R13)
      RETURN
C
C
C  BETA = 0 OR 180  CAN ONLY FIND ALPHA + OR - GAMMA.
13    BETA=0.
      IF(TRSYM(3,3).LT.0.) BETA=180.
      GAMMA=CONV*ATAN2(TRSYM(2,1),TRSYM(2,2))
      ALPHA=0.
      RETURN
      END
C
C
C
      SUBROUTINE DCOSFD(ROT,DIRCOS,KAPPA)
C     ===================================
C******************************************************
C       Given a rotation matrix ROTN, expressing a rotation
C       through an angle KAPPA right-handedly about an axis
C       with direction cosines DIRCOS(I),
C       DCOSFD determines DIRCOS() and KAPPA
C       KAPPA is returned in degrees in the range 0 to 180
C
C       Expression for rotation matrix is (see J & J, "Mathl.
C               Phys.", p. 122):-
C
C  cw + n(1)n(1)(1-cw)    n(1)n(2)(1-cw)-n(3)sw  n(1)n(3)(1-cw)+n(3)sw
C  n(1)n(2)(1-cw)+n(3)sw  cw + n(2)n(2)(1-cw)    n(2)n(3)(1-cw)-n(1)sw
C  n(3)n(1)(1-cw)-n(2)sw  n(3)n(2)(1-cw)+n(1)sw  cw + n(3)n(3)(1-cw)
C
C       where cw = cos(KAPPA), sw = sin(KAPPA),
C       & n(1), n(2), n(3) are the direction cosines.
C
C*******************************************************}
C
        REAL ROT(3,3),DIRCOS(3),S(3),P(3),D(3)
        REAL KAPPA
        REAL DIFF,PI,TRACE,R2,R3,SINKAP,COSKAP,TERM
        INTEGER NS(3),NP(3),NCS(3),ND(3),KNTRL,J,JMX
C
C
          DIFF=1.0E-04
          PI=4.0*ATAN(1.)
          KNTRL=0
C       Trace = 1 + 2.cos w
          TRACE=ROT(1,1)+ROT(2,2)+ROT(3,3)
C       S(1)=2.n(1).sin w
C       S(2)=2.n(2).sin w
C       S(3)=2.n(3).sin w
          S(1)=ROT(3,2)-ROT(2,3)
          S(2)=ROT(1,3)-ROT(3,1)
          S(3)=ROT(2,1)-ROT(1,2)
          CALL ORDR3(KNTRL,S,NS)
C--Is biggest S() zero (i.e. is sin w = 0)?
          IF (ABS(S(NS(1))).GT.DIFF) THEN
            R2=S(NS(2))/S(NS(1))
            R3=S(NS(3))/S(NS(1))
            DIRCOS(NS(1))=1.0/SQRT(1 + R2**2 + R3**2)
            DIRCOS(NS(2))=R2*DIRCOS(NS(1))
            DIRCOS(NS(3))=R3*DIRCOS(NS(1))
            GO TO 123
          END IF
C--Calculation when sin w = 0 (esp. when w = 180').
C       P(1)=n(2).n(3).{0(W=0)or2(w=180)}
C       P(2)=n(3).n(1).{0(W=0)or2(w=180)}
C       P(3)=n(1).n(2).{0(W=0)or2(w=180)}
          P(1)=ROT(3,2)
          P(2)=ROT(1,3)
          P(3)=ROT(2,1)
          CALL ORDR3(KNTRL,P,NP)
C--Is biggest P() zero (i.e. are all off-diag. terms zero)?
          IF (ABS(P(NP(1))).LT.DIFF) THEN
            IF (TRACE.GT.0.0) THEN
C--Matrix is unit matrix.
                KAPPA=0.0
                DIRCOS(1)=0.0
                DIRCOS(2)=0.0
                DIRCOS(3)=1.0
                RETURN
            END IF
C--Trace -ve, so dyad about x,y or z.
            KAPPA=PI
            DO 40 J=1,3
 40           D(J)=ROT(J,J)+1.0
            CALL ORDR3(KNTRL,D,ND)
            DIRCOS(ND(1))=1.0
            DIRCOS(ND(2))=0.0
            DIRCOS(ND(3))=0.0
            RETURN
          END IF
        IF (ABS(P(NP(2))).LT.DIFF) THEN
C--One d.c. is zero.
          DIRCOS(NP(1))=0.0
          DIRCOS(NP(2))=SQRT(0.5*(1.0+ROT(NP(2),NP(2))))
          DIRCOS(NP(3))=SQRT(0.5*(1.0+ROT(NP(3),NP(3))))
          IF (P(NP(1)).LT.0.0) DIRCOS(NP(2))=-DIRCOS(NP(2))
          GO TO 123
        END IF
C--All d.c's are nonzero.
C       R2=n(NP(1))/n(NP(2))
C       R3=n(NP(1))/n(NP(3))
          R2=P(NP(2))/P(NP(1))
          R3=P(NP(3))/P(NP(1))
          DIRCOS(NP(2))=1.0/SQRT(1 + R2**2 + (R2/R3)**2)
          DIRCOS(NP(3))=1.0/SQRT(1 + R3**2 + (R3/R2)**2)
C--Check: don't take sqrt of negative values
          TERM=1.0-DIRCOS(NP(2))**2-DIRCOS(NP(3))**2
          IF (TERM.LT.0.0) THEN
            IF (ABS(TERM).GT.DIFF)
     +          CALL CCPERR (1,'Negative value in DCOSFD')
            TERM=ABS(TERM)
          END IF
          DIRCOS(NP(1))=SQRT(TERM)
C--Adjust signs of DIRCOS().
          JMX=0
          DO 20 J=1,3
            IF (P(J).GT.0.0) THEN
              IF (JMX.GT.0) JMX=-1
              IF (JMX.EQ.0) JMX=J
            END IF
 20       CONTINUE
          IF (JMX.GT.0) DIRCOS(JMX)=-DIRCOS(JMX)
          IF (JMX.EQ.0) WRITE(*,1358)
 1358     FORMAT(' All P()s. are negative or zero')
C--Given d.c's, calculate KAPPA.
123       KNTRL=1
          CALL ORDR3(KNTRL,DIRCOS,NCS)
C       Find KAPPA.
          SINKAP=ROT(NCS(3),NCS(2))-ROT(NCS(2),NCS(3))
          SINKAP=0.5*SINKAP/DIRCOS(NCS(1))
          COSKAP=0.5*(TRACE-1.0)
          KAPPA=ATAN2(SINKAP,COSKAP)
          KAPPA=KAPPA*180.0/PI
C  If kappa negative, negate kappa and invert DIRCOS
          IF (KAPPA .LT. 0.0) THEN
             KAPPA = -KAPPA
             DIRCOS(1) = -DIRCOS(1)
             DIRCOS(2) = -DIRCOS(2)
             DIRCOS(3) = -DIRCOS(3)
          ENDIF
        RETURN
        END
C
C
C
C
C
      SUBROUTINE ORDR3(KNTRL,Q,NQ)
C     ============================
C
C*************************************************************
C       Find NQ(1) so that |Q(NQ(1)| is biggest Q.
C       If KNTRL=1, make NQ(1),NQ(2),NQ(3) a positive permutn.
C       If KNTRL= anything else,
C       find nos. NQ() so that|Q(NQ(1)|>|Q(NQ(2))|>|Q(NQ(3))|
C
C*************************************************************
        REAL Q(3)
        INTEGER NQ(3),KNTRL,J
C
C       Find no. NQ(1) of Q() with largest |Q()|.
          NQ(1)=1
          IF (ABS(Q(2)).GT.ABS(Q(1))) NQ(1)=2
          IF (ABS(Q(3)).GT.ABS(Q(NQ(1)))) NQ(1)=3
C       Find nos. for other two Q()'s & order so |Q(NQ(2))|>|Q(NQ(3))|
          NQ(2)=MOD(NQ(1),3)+1
          NQ(3)=MOD(NQ(2),3)+1
          IF (KNTRL.EQ.1) GO TO 100
            IF (ABS(Q(NQ(3))).GT.ABS(Q(NQ(2)))) THEN
                J=NQ(3)
                NQ(3)=NQ(2)
                NQ(2)=J
            END IF
 100    RETURN
        END
C
C
C
      SUBROUTINE POLANG(DC,OMEGA,PHI)
C     ===============================
C
C Find polar angles OMEGA & PHI (degrees) from direction cosines DC
C  OMEGA in range 0 to 180, PHI in range 0 to 360
C
      REAL DC(3), OMEGA, PHI
C
      REAL CONV, SOMEGA
C
      CONV=45./ATAN(1.)
C
       SOMEGA=SQRT(ABS(1.-DC(3)*DC(3)))
       OMEGA=CONV*ATAN2(SOMEGA,DC(3))
       IF(DC(1).EQ.0..AND.DC(2).EQ.0.) THEN
             PHI=0.0
       ELSE
             PHI=CONV*ATAN2(DC(2),DC(1))
       ENDIF
C
      RETURN
      END
C
C
C
      SUBROUTINE RNGANG(A)
C     ====================
C
C Put angle A into range 0-360 degrees
C
      REAL A
C
1     IF (A .LT. 0.0) THEN
            A = A + 360.0
            GO TO 1
      ENDIF
2     IF(A .GE. 360.0) THEN
            A = A - 360.0
            GO TO 2
      ENDIF
      RETURN
      END
C
C
C
      SUBROUTINE POLMAT(ROTN,ANGLES)
C     ==============================
C
C**************************************************************
C       Sets rotation matrix ROTN expressing a rotation through an
C       angle AKAPPA right-handedly about an axis
C       with polar angles OMEGA, PHI (OMEGA with Z-axis, projection
C       giving angle PHI with X-axis).
C       These angles give direction cosines DIRCOS(I).
C***************************************************************
C
      REAL ROTN(3,3),DIRCOS(3),ANGLES(3)
      REAL OMEGA,PHI,AKAPPA
C
C
      CONV=ATAN(1.)/45.
      OMEGA=ANGLES(1)
      PHI=ANGLES(2)
      AKAPPA=ANGLES(3)
      SNOM=SIN(OMEGA*CONV)
      CSOM=COS(OMEGA*CONV)
      SNPH=SIN(PHI*CONV)
      CSPH=COS(PHI*CONV)
      SNKA=SIN(AKAPPA*CONV)
      CSKA=COS(AKAPPA*CONV)
      DIRCOS(1)=SNOM*CSPH
      DIRCOS(2)=SNOM*SNPH
      DIRCOS(3)=CSOM
      DO 40 I=1,3
            DO 30 J=1,3
                K1=6-I-J
                K=K1
                IF ((K1.LT.1).OR.(K1.GT.3)) K=3
                EPSIJK=((I-J)*(J-K)*(K-I))/2
                ROTN(I,J)=DIRCOS(I)*DIRCOS(J)*
     $             (1.0-CSKA)-EPSIJK*DIRCOS(K)*SNKA
                IF (I.EQ.J) ROTN(I,J)=ROTN(I,J)+CSKA
 30         CONTINUE
 40       CONTINUE
      RETURN
      END