**==aa0001.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**   Program UNSCRAM     **
C***************************
C
C Purpose: To unscramble intensities for harmonics
C
C Input: .ge files and a .vc file
C        ( .ge files, i.e. .ge1 and .ge2 are produced by GENLAUE
C          and .vc (Victoreen Coefficient) files produced by AFSCALE)
C
      INCLUDE 'implicit_none.finc'
C
      INCLUDE 'overlaps.finc'
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C Var to select box or profile intensities
      COMMON /BP/ BOXPROF
      INTEGER   BOXPROF
C
C====== FUNCTIONS
C
      INTEGER LUNSTI, LUNSTO
C
C====== LOCALS:
C
      DOUBLE PRECISION     PARAMS(2,3,5)
 
      INTEGER  END_REC, I, END_POS
      INTEGER    IFAIL
 
      REAL     ACELL, BCELL, CCELL, CELLANG(3)
      REAL     ITOP, IBOT
      REAL     N_SIGS
      REAL     WAVMIN, WAVMAX
 
      INTEGER  NUM_OVERLAPS
      REAL     OVERLAP_INTENS(6,MAX_OVERLAPS)
      REAL     OVERLAP_SIGMAS(6,MAX_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(MAX_OVERLAPS)
      REAL     OVERLAP_DHKLS(MAX_OVERLAPS)
      REAL     OVERLAP_XY(2,MAX_OVERLAPS)
      INTEGER  OVERLAP_HKLS(MAX_OVERLAPS)
 
      CHARACTER*30  FNAME
 
      LOGICAL    SUPPRESS_OUTPUT, QUIT
C
C======  Set sensible values
C
      DATA  CELLANG /90., 90., 90./
      DATA  ITOP, IBOT /9000.0,100.0/
      DATA  N_SIGS /5/
      DATA  WAVMIN, WAVMAX /0.3,2.5/
C
C------------------------------------------------------------------------------
C
      CALL TRWSET(LUNSTI(0),LUNSTO(0),8,8,0)
      IN_TERM = LUNSTI(0)
      OUT_TERM = LUNSTO(0)
      IN_DATA = 2
      OUT_DATA = 1
      IFAIL = 0
      CALL CCPDPN(-8,'UNSCLOG','NEW','F',0,IFAIL)
      CALL TWRITE(' ')
      CALL TWRITE('Program:     UNSCRAM')
      CALL TWRITE('Last update: 06/01/92')
      CALL TWRITE('Authors:     S.Zurek, M.Z.Papiz, I.J.Clifton')
      CALL TWRITE(' ')
C
C====== Get generate file name
C
 20   CALL TPROMP('Enter root name for .ge1/.ge2 files: ')
      CALL TREAD(FNAME,QUIT)
      IF (FNAME.EQ.' ') GO TO 20
C
C====== Read in the Victoreen parameters
C
      CALL READ_VC(PARAMS)
C
C====== Open .ge files and find last record
C
      CALL OPEN_GE(IN_DATA,3,FNAME,END_POS,END_REC)
C
C====== Get rejection criteria and unit cell parameters
C
      CALL GET_REJECT_CRIT(ACELL,BCELL,CCELL,CELLANG,ITOP,IBOT,N_SIGS,
     +                     WAVMIN,WAVMAX)
C
C====== Use box or profile intensities?
C
      BOXPROF = 2
 30   CALL TASKI('Use box (1) or profile fitted (2) intensities',
     +           BOXPROF,QUIT)
      IF (BOXPROF.LE.0 .OR. BOXPROF.GT.2) GO TO 30
C
C====== Does user want to suppress bulk of output pouring out on the screen?
C
      SUPPRESS_OUTPUT = .TRUE.
      CALL TASKQ('Do you want to suppress output to terminal',
     +           SUPPRESS_OUTPUT,QUIT)
C
C====== Now unscramble all the spots
C
      CALL UNSCRAM_SPOTS(MAX_OVERLAPS,END_REC,ACELL,BCELL,CCELL,CELLANG,
     +                   ITOP,IBOT,N_SIGS,WAVMIN,WAVMAX,PARAMS,
     +                   NUM_OVERLAPS,OVERLAP_INTENS,OVERLAP_SIGMAS,
     +                   OVERLAP_LAMBDAS,OVERLAP_HKLS,OVERLAP_DHKLS,
     +                   OVERLAP_XY,SUPPRESS_OUTPUT)
C
      IF (NUM_OVERLAPS.GT.0) THEN
         CALL ANALYSE_UNSCRAM(NUM_OVERLAPS,OVERLAP_INTENS,
     +                        OVERLAP_LAMBDAS,OVERLAP_HKLS,
     +                        OVERLAP_DHKLS,SUPPRESS_OUTPUT)
C
         CALL OUTPUT_UNSCRAM(NUM_OVERLAPS,OVERLAP_INTENS,
     +                       OVERLAP_LAMBDAS,OVERLAP_HKLS,OVERLAP_DHKLS,
     +                       OVERLAP_SIGMAS,OVERLAP_XY,PARAMS,ACELL,
     +                       BCELL,CCELL,CELLANG)
      ELSE
         STOP ' No overlaps found'
      END IF
C
      STOP
      END
**==analyse_unscram.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************************
C**       ANALYSE_UNSCRAM       **
C*********************************
C
C
      SUBROUTINE  ANALYSE_UNSCRAM (NUM_OVERLAPS,
     +                             OVERLAP_INTENS, OVERLAP_LAMBDAS,
     +                             OVERLAP_HKLS, OVERLAP_DHKLS,
     +                             SUPPRESS_OUTPUT)
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  NUM_OVERLAPS
      REAL     OVERLAP_INTENS(6,NUM_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(NUM_OVERLAPS)
      INTEGER  OVERLAP_HKLS(NUM_OVERLAPS)
      REAL     OVERLAP_DHKLS(NUM_OVERLAPS)
      LOGICAL    SUPPRESS_OUTPUT
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      INTEGER  MAX_PAIRS, NUM_PAIRS
      PARAMETER (MAX_PAIRS=1000)
      INTEGER  SPOT1(MAX_PAIRS), SPOT2(MAX_PAIRS)
      INTEGER  S1, S2, I
 
      REAL     R_REJECT_THRESHOLD
 
      LOGICAL    REPLY, QUIT, OUTPUT_PAIRS
 
      CHARACTER*80 STR
C
      INTEGER       H,K,L, H2,K2,L2
      INTEGER  HKL(3), HKL2(3)
C
      EQUIVALENCE  (HKL(1),H)
      EQUIVALENCE  (HKL(2),K)
      EQUIVALENCE  (HKL(3),L)
C
      EQUIVALENCE  (HKL2(1),H2)
      EQUIVALENCE  (HKL2(2),K2)
      EQUIVALENCE  (HKL2(3),L2)
C
      DATA  R_REJECT_THRESHOLD /20.0/
C
C------------------------------------------------------------------------------
C
 20   CALL PAIR_UP(NUM_OVERLAPS,OVERLAP_INTENS,OVERLAP_HKLS,
     +             OVERLAP_LAMBDAS,MAX_PAIRS,NUM_PAIRS,SPOT1,SPOT2)
C
      IF (NUM_PAIRS.EQ.0) THEN
         CALL TWRITE(' ')
         CALL TWRITE('*** No Pairs ***')
         CALL TWRITE(' ')
      ELSE
C
         CALL TASKF('Enter R-factor rejection threshold',
     +              R_REJECT_THRESHOLD,1,QUIT)
C
         OUTPUT_PAIRS = .FALSE.
         CALL TASKQ('Output the pairs for the analysis to the log file',
     +              OUTPUT_PAIRS,QUIT)
C
         IF (.NOT.SUPPRESS_OUTPUT) WRITE (OUT_TERM,*) ' Pairs are:'
         IF (OUTPUT_PAIRS) WRITE (8,*) ' Pairs are:'
         DO 30 I = 1 , NUM_PAIRS
            S1 = SPOT1(I)
            S2 = SPOT2(I)
            CALL CCPBYI(HKL,OVERLAP_HKLS(S1),-3)
            CALL CCPBYI(HKL2,OVERLAP_HKLS(S2),-3)
            IF (.NOT.SUPPRESS_OUTPUT) WRITE (OUT_TERM,6010) S1 , H , K , 
     +          L , OVERLAP_LAMBDAS(S1) , S2 , H2 , K2 , L2 , 
     +          OVERLAP_LAMBDAS(S2)
            IF (OUTPUT_PAIRS) WRITE (8,6010) S1 , H , K , L , 
     +                               OVERLAP_LAMBDAS(S1) , S2 , H2 , 
     +                               K2 , L2 , OVERLAP_LAMBDAS(S2)
 30      CONTINUE
         CALL TWRITE(' ')
         WRITE (STR,6020) NUM_PAIRS
         CALL TWRITE(STR)
         CALL TWRITE(' ')
C
         CALL CALC_RSYMM(NUM_OVERLAPS,OVERLAP_INTENS,OVERLAP_LAMBDAS,
     +                   OVERLAP_HKLS,OVERLAP_DHKLS,NUM_PAIRS,SPOT1,
     +                   SPOT2,R_REJECT_THRESHOLD)
C
      END IF
C
      REPLY = .FALSE.
      CALL TASKQ('Do you want to change some parameters',REPLY,QUIT)
C
      IF (REPLY) GO TO 20
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT (2('Reflection:',I4,5X,3I5,F9.3,10X))
 6020 FORMAT ('A total of',I6,' pairs.')
      END
**==apply_reject_crit.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C************************************
C**       APPLY_REJECT_CRIT        **
C************************************
C
C
      INTEGER FUNCTION  APPLY_REJECT_CRIT ( MULTIPLICITY,
     +      INTENSITIES, SIGMAS, LAMBDAS, SPATIAL_OVERLAP,
     +      ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX,
     +      NUM_GOOD_FILMS )
C
C Purpose: Returns 1 if spot is rejected, 0 if OK
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  MULTIPLICITY
      REAL     INTENSITIES(6), SIGMAS(6), LAMBDAS(MULTIPLICITY)
      INTEGER  SPATIAL_OVERLAP
      REAL     ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX
      INTEGER  NUM_GOOD_FILMS
C
C ITOP : maximum intensity allowed
C IBOT : minimum intensity allowed
C
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      INTEGER  I
      REAL     WAV
C
C------------------------------------------------------------------------------
C
      APPLY_REJECT_CRIT = 0
C
      IF (SPATIAL_OVERLAP.EQ.0) THEN
C
C====== Check intensities
C
         NUM_GOOD_FILMS = 6
         DO 20 I = 1 , 6
            IF (INTENSITIES(I).GT.ITOP .OR. INTENSITIES(I).LT.IBOT .OR. 
     +          INTENSITIES(I).LT.N_SIGS*ABS(SIGMAS(I))) THEN
               INTENSITIES(I) = -1.0*ABS(INTENSITIES(I))
               NUM_GOOD_FILMS = NUM_GOOD_FILMS - 1
            END IF
 20      CONTINUE
C
C====== Check wavelengths of harmonics
C
         DO 30 I = 1 , MULTIPLICITY
C
            WAV = LAMBDAS(I)
C
            IF (WAV.LT.WAVMIN .OR. WAV.GT.WAVMAX) GO TO 40
C
            IF ((WAV.GE.0.486 .AND. WAV.LE.0.490) .OR. 
     +          (WAV.GE.0.900 .AND. WAV.LE.0.920)) GO TO 40
C
 30      CONTINUE
C
         RETURN
      END IF
C
C====== Failed test
C
 40   APPLY_REJECT_CRIT = 1
      RETURN
      END
**==calc_rfac.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**       CALC_RFAC        **
C****************************
C
C
      SUBROUTINE  CALC_RFAC (MEAN, DIFF, N, R)
C
C Purpose: To calculate R factor from given number of means and diffs
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  N
      REAL     MEAN(12,3,N), DIFF(12,3,N), R
C
C====== LOCALS:
C
      REAL     TOT_MEAN, TOT_DIFF
      INTEGER  I, J
C
C------------------------------------------------------------------------------
C
      TOT_MEAN = 0.0
      TOT_DIFF = 0.0
C
      DO 30 J = 1 , N
         DO 20 I = 1 , 12
            TOT_MEAN = TOT_MEAN + MEAN(I,1,J)
            TOT_DIFF = TOT_DIFF + DIFF(I,1,J)
 20      CONTINUE
 30   CONTINUE
C
      IF (TOT_MEAN.NE.0.0) THEN
         R = TOT_DIFF/TOT_MEAN*100.0
      ELSE
         R = -1.0
      END IF
C
      RETURN
      END
**==calc_rsymm.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       CALC_RSYMM        **
C*****************************
C
C
      SUBROUTINE  CALC_RSYMM (NUM_OVERLAPS,
     +                  OVERLAP_INTENS, OVERLAP_LAMBDAS,
     +                  OVERLAP_HKLS, OVERLAP_DHKLS,
     +                  NUM_PAIRS, SPOT1, SPOT2,
     +                  R_REJECT_THRESHOLD)
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  NUM_OVERLAPS
      REAL     OVERLAP_INTENS(6,NUM_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(NUM_OVERLAPS)
      INTEGER  OVERLAP_HKLS(NUM_OVERLAPS)
      REAL     OVERLAP_DHKLS(NUM_OVERLAPS)
      INTEGER  NUM_PAIRS
      INTEGER  SPOT1(NUM_PAIRS), SPOT2(NUM_PAIRS)
      REAL     R_REJECT_THRESHOLD
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      INTEGER  I, J, FILM
      INTEGER  IBIN, LB, DB, BIN(3), BIN_TYPE
      REAL     MEAN(12,3,6,5)
      REAL     DIFF(12,3,6,5)
      INTEGER  N(12,3,6,5)
      CHARACTER*80 STR
      REAL     R
      INTEGER  RB, R_BIN(9,6,5)
C
      INTEGER  S1, S2
      REAL     I1, I2, LAMBDA, DHKL, IMEAN, IDIFF
C
      INTEGER  MB, M_BIN1, M_BIN2
      INTEGER  H, K, L, M_BIN, HKL(4)
      INTEGER  H1, K1, L1, H2, K2, L2
C
      EQUIVALENCE  (HKL(1),H),(HKL(2),K),
     + (HKL(3),L),(HKL(4),M_BIN)
C
      INTEGER  NUM_DHKL_RANGES
      PARAMETER (NUM_DHKL_RANGES=9)
      REAL     DHKL_RANGES(NUM_DHKL_RANGES)
C
      DATA  DHKL_RANGES /20.0,10.0,7.0,5.0,4.0,3.33,2.9,2.6,0.0/
C
C------------------------------------------------------------------------------
C
      DO 100 MB = 1 , 5
C
         DO 90 FILM = 1 , 6
C
            DO 30 J = 1 , 3
               DO 20 I = 1 , 12
                  MEAN(I,J,FILM,MB) = 0.0
                  DIFF(I,J,FILM,MB) = 0.0
                  N(I,J,FILM,MB) = 0
 20            CONTINUE
 30         CONTINUE
C
            DO 40 I = 1 , 9
               R_BIN(I,FILM,MB) = 0
 40         CONTINUE
C
            DO 80 I = 1 , NUM_PAIRS
C
C====== Do basics
C
               S1 = SPOT1(I)
               S2 = SPOT2(I)
               I1 = OVERLAP_INTENS(FILM,S1)
               I2 = OVERLAP_INTENS(FILM,S2)
               LAMBDA = (OVERLAP_LAMBDAS(S1)+OVERLAP_LAMBDAS(S2))/2.0
               DHKL = (OVERLAP_DHKLS(S1)+OVERLAP_DHKLS(S2))/2.0
C
C====== Reject overlap if not in current multiplicity bin
C
               CALL CCPBYI(HKL,OVERLAP_HKLS(S1),-4)
               M_BIN1 = M_BIN
               M_BIN1 = IABS(M_BIN1)
               CALL CCPBYI(HKL,OVERLAP_HKLS(S2),-4)
               M_BIN2 = M_BIN
               M_BIN2 = IABS(M_BIN2)
C
C====== Check for pairs that have different multiplicities.
C
               IF (M_BIN1.NE.M_BIN2) THEN
                  CALL CCPBYI(HKL,OVERLAP_HKLS(S1),-4)
                  H1 = H
                  K1 = K
                  L1 = L
                  CALL CCPBYI(HKL,OVERLAP_HKLS(S2),-4)
                  H2 = H
                  K2 = K
                  L2 = L
                  WRITE (STR,6040)
                  CALL TWRITE(STR)
                  WRITE (STR,6050) S1 , S2
                  CALL TWRITE(STR)
                  WRITE (STR,6060) M_BIN1 , M_BIN2 , H1 , K1 , L1 , H2 , 
     +                             K2 , L2 , DHKL , LAMBDA
                  CALL TWRITE(STR)
                  GO TO 80
               END IF
               IF (M_BIN1.EQ.MB+1) THEN
C
C====== Reject negative intensities
C
                  IF (I1.GE.0.0 .AND. I2.GE.0.0) THEN
C
                     IMEAN = (I1+I2)/2.0
                     IDIFF = ABS(I1-I2)/2.0
C
C====== If Imean is zero then something funny is going on...
C
                     IF (IMEAN.EQ.0.0) THEN
                        WRITE (STR,6020) S1 , S2
                        CALL TWRITE(STR)
                        WRITE (STR,6030) I1 , I2
                        CALL TWRITE(STR)
                        CALL TWRITE(' ')
                        GO TO 80
                     END IF
C
C====== Work out individual R-factor
C
                     R = IDIFF/IMEAN*100.0
C
C====== Is it above rejection threshold ? - if so then discard from analysis
C
                     IF (R_REJECT_THRESHOLD.EQ.0.0 .OR. 
     +                   R.LE.R_REJECT_THRESHOLD) THEN
C
C====== Reflection is O.K.
C
                        RB = R/2.5 + 1
                        RB = MIN(MAX(RB,0),9)
                        R_BIN(RB,FILM,MB) = R_BIN(RB,FILM,MB) + 1
C
C====== Which bins?
C
                        IBIN = IMEAN/200.0 + 1
                        IF (IBIN.GT.11) IBIN = 11
                        BIN(1) = IBIN
C
                        LB = (LAMBDA-0.1)/0.2 + 1
                        IF (LB.LT.1) LB = 1
                        IF (LB.GT.12) LB = 12
                        BIN(2) = LB
C
                        DO 50 J = 1 , NUM_DHKL_RANGES
                           IF (DHKL.GT.DHKL_RANGES(J)) GO TO 60
 50                     CONTINUE
                        WRITE (STR,6010) DHKL
                        CALL TWRITE(STR)
                     END IF
                  END IF
               END IF
               GO TO 80
 60            DB = J
               BIN(3) = DB
C
C====== Now add to sub-totals
C
               DO 70 BIN_TYPE = 1 , 3
                  MEAN(BIN(BIN_TYPE),BIN_TYPE,FILM,MB)
     +               = MEAN(BIN(BIN_TYPE),BIN_TYPE,FILM,MB) + IMEAN
                  DIFF(BIN(BIN_TYPE),BIN_TYPE,FILM,MB)
     +               = DIFF(BIN(BIN_TYPE),BIN_TYPE,FILM,MB) + IDIFF
                  N(BIN(BIN_TYPE),BIN_TYPE,FILM,MB)
     +               = N(BIN(BIN_TYPE),BIN_TYPE,FILM,MB) + 1
 70            CONTINUE
C
 80         CONTINUE
C
 90      CONTINUE
C
 100  CONTINUE
C
C====== Now display results
C
      CALL DISP_RSYMM(MEAN,DIFF,N,R_BIN)
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Dhkl value zero or negative. Value =',G14.6)
 6020 FORMAT ('***  Imean zero at reflections',I6,' and',I6,'  ***')
 6030 FORMAT ('     Intensities are ',2G14.6)
 6040 FORMAT ('Found pair with different multiplicity')
 6050 FORMAT ('Spots: ',2I6,' Multiplicities: ',2I2)
 6060 FORMAT ('Indices: ',6I4,' D-hkl: ',F8.4,' Lambda: ',F8.4)
      END
**==disp_rsymm.f  
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE  DISP_RSYMM (MEAN, DIFF, N, R_BIN)
C
C Purpose: To calculate and display R-factors
C
      INCLUDE 'implicit_none.finc'
C
      REAL     MEAN(12,3,6,5)
      REAL     DIFF(12,3,6,5)
      INTEGER  N(12,3,6,5)
C
      INTEGER  R_BIN(9,6,5)
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      REAL     SUM, DENOM
      INTEGER  ISUM
C
      INTEGER  I, J, MB, FILM, TOTAL_N
      INTEGER  N_R_BIN(9), NR(12)
      REAL     R(12), RF
      CHARACTER*80 STR
      EXTERNAL SUM
C
C------------------------------------------------------------------------------
C
      DO 100 MB = 1 , 5
C
         CALL TWRITE(' ')
         CALL TWRITE(' ')
         CALL TWRITE(' ')
         WRITE (STR,6010) MB + 1
         CALL TWRITE(STR)
         WRITE (STR,6020)
         CALL TWRITE(STR)
         CALL TWRITE(' ')
         CALL TWRITE(' ')
C
         DO 50 FILM = 1 , 6
C
C====== Check to see if there are any pairs for this film and multiplicity
C====== by calculating overall R-factor for film
C
            CALL CALC_RFAC(MEAN(1,1,FILM,MB),DIFF(1,1,FILM,MB),1,RF)
            IF (RF.EQ.-1.0) THEN
               CALL TWRITE(' ')
               WRITE (STR,6260) FILM
               CALL TWRITE(STR)
               CALL TWRITE(' ')
               GO TO 50
            END IF
C
            DO 20 I = 1 , 11
               IF (MEAN(I,1,FILM,MB).NE.0.0) THEN
                  R(I) = DIFF(I,1,FILM,MB)/MEAN(I,1,FILM,MB)*100.0
               ELSE
                  R(I) = -1.0
               END IF
 20         CONTINUE
C
            CALL TWRITE(' ')
            WRITE (STR,6030) FILM
            CALL TWRITE(STR)
            WRITE (STR,6040)
            CALL TWRITE(STR)
            WRITE (STR,6050)
            CALL TWRITE(STR)
            WRITE (STR,6060) (N(I,1,FILM,MB),I=1,11)
            CALL TWRITE(STR)
            WRITE (STR,6070) (R(I),I=1,11)
            CALL TWRITE(STR)
C
            DO 30 I = 1 , 12
               IF (MEAN(I,2,FILM,MB).NE.0.0) THEN
                  R(I) = DIFF(I,2,FILM,MB)/MEAN(I,2,FILM,MB)*100.0
               ELSE
                  R(I) = -1.0
               END IF
 30         CONTINUE
C
            CALL TWRITE(' ')
            WRITE (STR,6120)
            CALL TWRITE(STR)
            WRITE (STR,6130)
            CALL TWRITE(STR)
            WRITE (STR,6140) (N(I,2,FILM,MB),I=1,12)
            CALL TWRITE(STR)
            WRITE (STR,6150) (R(I),I=1,12)
            CALL TWRITE(STR)
C
            DO 40 I = 1 , 9
               IF (MEAN(I,3,FILM,MB).NE.0.0) THEN
                  R(I) = DIFF(I,3,FILM,MB)/MEAN(I,3,FILM,MB)*100.0
               ELSE
                  R(I) = -1.0
               END IF
 40         CONTINUE
C
            CALL TWRITE(' ')
            WRITE (STR,6160)
            CALL TWRITE(STR)
            WRITE (STR,6170)
            CALL TWRITE(STR)
            WRITE (STR,6180) (N(I,3,FILM,MB),I=1,9)
            CALL TWRITE(STR)
            WRITE (STR,6190) (R(I),I=1,9)
            CALL TWRITE(STR)
C
C====== Display overall R-factor for film
C
            CALL TWRITE(' ')
            WRITE (STR,6200) RF
            CALL TWRITE(STR)
C
 50      CONTINUE
C
C====== Calculate overall R-factor for pack
C
         CALL CALC_RFAC(MEAN(1,1,1,MB),DIFF(1,1,1,MB),6,RF)
         IF (RF.LT.0.0) THEN
            CALL TWRITE(' ')
            WRITE (STR,6270)
            CALL TWRITE(STR)
            CALL TWRITE(' ')
            GO TO 100
         END IF
         CALL TWRITE(' ')
         WRITE (STR,6210) RF
         CALL TWRITE(STR)
C
C====== Display R-factors over a pack
C
         DO 60 I = 1 , 11
            DENOM = SUM(MEAN(I,1,1,MB),33,6)
            IF (DENOM.NE.0.0) THEN
               R(I) = SUM(DIFF(I,1,1,MB),33,6)/DENOM*100.0
            ELSE
               R(I) = -1.0
            END IF
            NR(I) = ISUM(N(I,1,1,MB),33,6)
 60      CONTINUE
         CALL TWRITE(' ')
         CALL TWRITE(' ')
         WRITE (STR,6080)
         CALL TWRITE(STR)
         CALL TWRITE(' ')
         WRITE (STR,6090)
         CALL TWRITE(STR)
         WRITE (STR,6100) (NR(I),I=1,11)
         CALL TWRITE(STR)
         WRITE (STR,6110) (R(I),I=1,11)
         CALL TWRITE(STR)
C
         DO 70 I = 1 , 12
            DENOM = SUM(MEAN(I,2,1,MB),33,6)
            IF (DENOM.NE.0.0) THEN
               R(I) = SUM(DIFF(I,2,1,MB),33,6)/DENOM*100.0
            ELSE
               R(I) = -1.0
            END IF
            NR(I) = ISUM(N(I,2,1,MB),33,6)
 70      CONTINUE
         CALL TWRITE(' ')
         WRITE (STR,6120)
         CALL TWRITE(STR)
         WRITE (STR,6130)
         CALL TWRITE(STR)
         WRITE (STR,6140) (NR(I),I=1,12)
         CALL TWRITE(STR)
         WRITE (STR,6150) (R(I),I=1,12)
         CALL TWRITE(STR)
C
         DO 80 I = 1 , 9
            DENOM = SUM(MEAN(I,3,1,MB),33,6)
            IF (DENOM.NE.0.0) THEN
               R(I) = SUM(DIFF(I,3,1,MB),33,6)/DENOM*100.0
            ELSE
               R(I) = -1.0
            END IF
            NR(I) = ISUM(N(I,3,1,MB),33,6)
 80      CONTINUE
         CALL TWRITE(' ')
         WRITE (STR,6160)
         CALL TWRITE(STR)
         WRITE (STR,6170)
         CALL TWRITE(STR)
         WRITE (STR,6180) (NR(I),I=1,9)
         CALL TWRITE(STR)
         WRITE (STR,6190) (R(I),I=1,9)
         CALL TWRITE(STR)
C
C====== Display R-factor distribution over pack
C
         DO 90 I = 1 , 9
            N_R_BIN(I) = ISUM(R_BIN(I,1,MB),9,6)
 90      CONTINUE
         CALL TWRITE(' ')
         WRITE (STR,6230)
         CALL TWRITE(STR)
         WRITE (STR,6240)
         CALL TWRITE(STR)
         WRITE (STR,6250) (N_R_BIN(I),I=1,9)
         CALL TWRITE(STR)
         CALL TWRITE(' ')
C
 100  CONTINUE
C
C====== Work out mega R-factor across pack and multiplicity
C
      CALL CALC_RFAC(MEAN(1,1,1,1),DIFF(1,1,1,1),6*5,RF)
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      WRITE (STR,6220) RF
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('  Multiplicity = ',I2)
 6020 FORMAT ('  -----------------')
 6030 FORMAT ('Film  ',I2)
 6040 FORMAT ('R-symms on intensity')
 6050 FORMAT ('  200   400   600   800  1000  1200  1400  1600',
     +        ' 1800  2000 >2000')
 6060 FORMAT (11I6)
 6070 FORMAT (11F6.1)
 6080 FORMAT ('For the PACK')
 6090 FORMAT ('  200   400   600   800  1000  1200  1400  1600',
     +        ' 1800  2000 >2000')
 6100 FORMAT (11I6)
 6110 FORMAT (11F6.1)
 6120 FORMAT ('R-symms on wavelength')
 6130 FORMAT (
     +     ' <0.3   0.5   0.7   0.9   1.1   1.3   1.5   1.7   1.9   2.1'
     +     ,'   2.3   >2.3')
 6140 FORMAT (12I6)
 6150 FORMAT (12F6.1)
 6160 FORMAT ('R-symms on dhkl')
 6170 FORMAT ('   >20    10     7     5     4  3.33   2.9   2.6  <2.6')
 6180 FORMAT (9I6)
 6190 FORMAT (9F6.1)
 6200 FORMAT ('Overall R-symm =',F6.1)
 6210 FORMAT ('Overall R-symm across pack =',F6.1)
 6220 FORMAT ('Overall R-symm across pack and all multiplicities = ',
     +        F6.1)
 6230 FORMAT ('R-symm distribution across pack')
 6240 FORMAT (' 0.0-2.5  2.5-5.0  5.0-7.5  7.5-10   10-12.5  ',
     +        '12.5-15  15-17.5  17.5-20  >20')
 6250 FORMAT (I6,8I9)
 6260 FORMAT ('   ***  No suitable reflections for film ',I2,'  ***')
 6270 FORMAT ('   ***  No suitable reflections at ',
     +        'this multiplicity  ***')
      END
**==get_ge_data.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       GET_GE_DATA        **
C******************************
C
C
      INTEGER FUNCTION  GET_GE_DATA (IN_DATA1, IN_DATA2, RECORD_NUM,
     +     X, Y, INTEN_BOX_S, SIG_BOX_S, INTEN_PROF_S, SIG_PROF_S,
     +     MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP,
     +     H_S, K_S, L_S, LAMBDAS, DMIN_THRESHOLD )
C
C Purpose: To read in data from a .ge format file
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  IN_DATA1, IN_DATA2, RECORD_NUM
      REAL     X, Y, INTEN_BOX_S(6), SIG_BOX_S(6),
     + INTEN_PROF_S(6), SIG_PROF_S(6)
      INTEGER  MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP
      INTEGER  H_S(MULTIPLICITY), K_S(MULTIPLICITY), L_S(MULTIPLICITY)
      REAL     LAMBDAS(MULTIPLICITY)
      REAL     DMIN_THRESHOLD
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      INTEGER NBITST
C
C====== LOCALS:
C
      INTEGER  I, REC_NUM, LINK_REC, MAX_MULTI
      REAL     BX, BY, BDMIN
      INTEGER  BINTENBS(6), BSIGBS(6), BINTENPS(6), BSIGPS(6)
      REAL     LAMBDA, LAMBDA1, LAMBDA2
      INTEGER       M_FLAGS, H, K, L, NIDX, H2, K2, L2
      INTEGER  IREAL
      REAL     RINT
      CHARACTER*80 STR
C
      INTEGER  BUFFER(19), BUFFER2(2)
      INTEGER       MHKL(4), HKL2(4)
C
      EQUIVALENCE (IREAL, RINT)
C
      EQUIVALENCE (BUFFER(1),BX)
      EQUIVALENCE (BUFFER(2),BY)
      EQUIVALENCE (BUFFER(16),LAMBDA1)
      EQUIVALENCE (BUFFER(18),LAMBDA2)
      EQUIVALENCE (BUFFER(19),BDMIN)
C
      EQUIVALENCE (MHKL(1),M_FLAGS)
      EQUIVALENCE (MHKL(2),H)
      EQUIVALENCE (MHKL(3),K)
      EQUIVALENCE (MHKL(4),L)
C
      EQUIVALENCE (HKL2(1),NIDX)
      EQUIVALENCE (HKL2(2),H2)
      EQUIVALENCE (HKL2(3),K2)
      EQUIVALENCE (HKL2(4),L2)
C
      EQUIVALENCE (BUFFER2(2),LAMBDA)
C
C------------------------------------------------------------------------------
C
      GET_GE_DATA = 0
C
      NODAL = 0
      SPATIAL_OVERLAP = 0
C
      MAX_MULTI = MULTIPLICITY
C
      READ (IN_DATA1,REC=RECORD_NUM) BUFFER
C
      X = BX
      Y = BY
C
      CALL CCPI2I(BINTENBS,BUFFER(3),6,.TRUE.,.FALSE.)
      CALL CCPI2I(BSIGBS,BUFFER(6),6,.TRUE.,.FALSE.)
      CALL CCPI2I(BINTENPS,BUFFER(9),6,.TRUE.,.FALSE.)
      CALL CCPI2I(BSIGPS,BUFFER(12),6,.TRUE.,.FALSE.)
      DO 20 I = 1 , 6
         INTEN_BOX_S(I) = BINTENBS(I)
         SIG_BOX_S(I) = BSIGBS(I)
         INTEN_PROF_S(I) = BINTENPS(I)
         SIG_PROF_S(I) = BSIGPS(I)
 20   CONTINUE
C
      CALL CCPBYI(MHKL,BUFFER(15),-4)
      CALL CCPBYI(HKL2,BUFFER(17),-4)
      H_S(1) = H
      K_S(1) = K
      L_S(1) = L
      H_S(2) = H2
      K_S(2) = K2
      L_S(2) = L2
C
      LAMBDAS(1) = LAMBDA1
      LAMBDAS(2) = LAMBDA2
C
      DMIN_THRESHOLD = BDMIN
C
      MULTIPLICITY = NBITST(M_FLAGS,0,5)
      IF (MULTIPLICITY.GT.MAX_MULTI) THEN
         WRITE (STR,6010) MULTIPLICITY - 1 , MAX_MULTI - 1
         CALL TWRITE(STR)
         MULTIPLICITY = MAX_MULTI
      END IF
C
      IF (NBITST(M_FLAGS,5,1).NE.0) THEN
         NODAL = 1
         NODIDX = NIDX
      ELSE
         NODIDX = 0
      END IF
      IF (NBITST(M_FLAGS,6,1).NE.0) SPATIAL_OVERLAP = 1
C
      IF (MULTIPLICITY.GT.2) THEN
         RINT = LAMBDAS(2)
         LINK_REC = IREAL
         DO 30 I = 2 , MULTIPLICITY
C
            READ (IN_DATA2,REC=LINK_REC) BUFFER2
C
            LINK_REC = LINK_REC + 1
C
            CALL CCPBYI(MHKL,BUFFER2(1),-4)
            H_S(I) = H
            K_S(I) = K
            L_S(I) = L
            LAMBDAS(I) = LAMBDA
C
 30      CONTINUE
      END IF
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Spot read that has',I3,' harmonics. Truncating ',
     +        'to just',I3,' harmonics.')
      END
**==get_overlaps.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**       GET_OVERLAPS        **
C*******************************
C
C
      INTEGER FUNCTION  GET_OVERLAPS ( IN_DATA1, IN_DATA2, RECORD_NUM,
     +          END_REC, MAX_MULTI,
     +          ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX,
     +          X, Y, INTENSITIES, SIGMAS, NUM_GOOD_FILMS,
     +          MULTIPLICITY, NODAL, SPATIAL_OVERLAP,
     +          H_S, K_S, L_S, LAMBDAS, C2TH )
C
C Purpose: Returns overlap data after applying rejection criteria
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  IN_DATA1, IN_DATA2, RECORD_NUM, END_REC, MAX_MULTI
C
      REAL     ITOP, IBOT, N_SIGS
      REAL     WAVMIN, WAVMAX
C
      REAL     X, Y, INTENSITIES(6), SIGMAS(6)
      INTEGER  NUM_GOOD_FILMS
      INTEGER  MULTIPLICITY, NODAL, SPATIAL_OVERLAP
      INTEGER  H_S(MAX_MULTI), K_S(MAX_MULTI), L_S(MAX_MULTI)
      REAL     LAMBDAS(MAX_MULTI), C2TH(MAX_MULTI)
C
C====== record_num : number of record to read from .ge1 file
C====== max_multi : maximum multiplicity of overlap allowed
C====== ITOP : maximum intensity allowed
C====== IBOT : minimum intensity allowed
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      COMMON /BP/ BOXPROF
      INTEGER   BOXPROF
C
C====== LOCALS:
C
      REAL     STHLSQ
      INTEGER  GET_GE_DATA, APPLY_REJECT_CRIT
      INTEGER  I, NODIDX
      REAL     WAV, INTENB(6), SIGB(6), INTENP(6), SIGP(6)
      REAL     H2, K2, L2, DSTAR, DHKL, THETA, SIN_THETA
      REAL     DMIN_THRESHOLD
C
C------------------------------------------------------------------------------
C
      GET_OVERLAPS = 0
C
C====== We want to know when mult>6 so need to use 7 here.
C
 20   MULTIPLICITY = 7
      RECORD_NUM = RECORD_NUM + 1
      IF (RECORD_NUM.GT.END_REC) THEN
         CALL TWRITE('Reached the end of file')
         GET_OVERLAPS = 1
         RETURN
      END IF
      IF (GET_GE_DATA(IN_DATA1,IN_DATA2,RECORD_NUM,X,Y,INTENB,SIGB,
     +    INTENP,SIGP,MULTIPLICITY,NODAL,NODIDX,SPATIAL_OVERLAP,H_S,K_S,
     +    L_S,LAMBDAS,DMIN_THRESHOLD).NE.0) GO TO 20
      IF (BOXPROF.EQ.1) THEN
         DO 30 I = 1 , 6
            INTENSITIES(I) = INTENB(I)
            SIGMAS(I) = SIGB(I)
 30      CONTINUE
      ELSE
         DO 40 I = 1 , 6
            INTENSITIES(I) = INTENP(I)
            SIGMAS(I) = SIGP(I)
 40      CONTINUE
      END IF
C
      IF (MULTIPLICITY.EQ.1) GO TO 20
C
C====== Now apply rejection criteria
C
      IF (APPLY_REJECT_CRIT(MULTIPLICITY,INTENSITIES,SIGMAS,LAMBDAS,
     +    SPATIAL_OVERLAP,ITOP,IBOT,N_SIGS,WAVMIN,WAVMAX,NUM_GOOD_FILMS)
     +    .EQ.1) GO TO 20
C
C====== If not enough good films then say so and reject
C
      IF (NUM_GOOD_FILMS.LT.MULTIPLICITY) GO TO 20
C
C====== Calculate cos(2*theta) for the harmonics
C
      DO 50 I = 1 , MULTIPLICITY
C
         WAV = LAMBDAS(I)
C
C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2
C====== NB This function uses implicit input from SETRSL
C
         DSTAR = 2*SQRT(STHLSQ(H_S(I),K_S(I),L_S(I)))
         DHKL = 1.0/DSTAR
         SIN_THETA = WAV*DSTAR/2.0
         THETA = ASIN(SIN_THETA)
         C2TH(I) = COS(2.0*THETA)
C
 50   CONTINUE
C
      RETURN
      END
**==get_reject_crit.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************************
C**       GET_REJECT_CRIT        **
C**********************************
C
C
      SUBROUTINE  GET_REJECT_CRIT (ACELL, BCELL, CCELL, CELLANG,
     +                             ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX)
C
C Purpose: To ge rejection criteria from user, as well as unit cell parameters
C
      INCLUDE 'implicit_none.finc'
C
      REAL     ACELL, BCELL, CCELL, CELLANG(3)
      REAL     ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      REAL  MILLER_INDEX, CP(6)
      LOGICAL QUIT
      INTEGER I
      CHARACTER*80 STR
C
C------------------------------------------------------------------------------
C
C====== Get unit cell parameters and wavelength range from .GE1 file.
C
      CALL READ_GE_PARAMS(CP,WAVMIN,WAVMAX)
      ACELL = CP(1)
      BCELL = CP(2)
      CCELL = CP(3)
      DO 20 I = 1 , 3
         CELLANG(I) = CP(I+3)
 20   CONTINUE
      CALL TWRITE(' ')
      WRITE (STR,6010)
      CALL TWRITE(STR)
      WRITE (STR,6020) ACELL , BCELL , CCELL , CELLANG
      CALL TWRITE(STR)
      WRITE (STR,6030) WAVMIN , WAVMAX
      CALL TWRITE(STR)
C
C====== Get intensity cut-offs
C
      CALL TASKF('Enter Itop for intensity range',ITOP,1,QUIT)
      CALL TASKF('Enter Ibot for intensity range',IBOT,1,QUIT)
C
C====== Get number of SD's before cut-off
C
      CALL TASKF('Enter number of SD`s before cut-off',N_SIGS,1,QUIT)
C
C====== This concludes the interactive input of parameters
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Cell Parameters:     a     b     c alpha  beta gamma')
 6020 FORMAT ('=============== ',6F6.1)
 6030 FORMAT ('Wavelength range   ',2F5.2)
      END
**==get_trans_factor.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***********************************
C**       GET_TRANS_FACTOR        **
C***********************************
C
C
      SUBROUTINE  GET_TRANS_FACTOR (PARAMS, DWAV, DC2TH, FILM, TF)
C
C Purpose: To calculate TRANSMISSION film factor for a given film and lambda
C
      INCLUDE 'implicit_none.finc'
C
      DOUBLE PRECISION     PARAMS(2,3,5), DWAV, DC2TH
      INTEGER  FILM
      DOUBLE PRECISION     TF
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      INTEGER  I, RANGE
      DOUBLE PRECISION     DWAV3, K, A, ARG, FF
      CHARACTER*80 STR
C
C FF : film factor (absorption)
C
C------------------------------------------------------------------------------
C
      TF = 1.0
      IF (FILM.EQ.1) RETURN
C
      IF (DWAV.LT.0.49) THEN
         RANGE = 1
      ELSE IF (DWAV.LT.0.92) THEN
         RANGE = 2
      ELSE IF (DWAV.GE.0.92) THEN
         RANGE = 3
      END IF
C
      DWAV3 = DWAV*DWAV*DWAV
C
      DO 20 I = 1 , FILM - 1
C
         K = PARAMS(1,RANGE,I)
         A = PARAMS(2,RANGE,I)
C
         ARG = A*DWAV3/DC2TH
         IF (ARG.GT.80) THEN
            WRITE (STR,6010)
            CALL TWRITE(STR)
            WRITE (STR,6020) ARG
            CALL TWRITE(STR)
            CALL TWRITE(' ')
            ARG = 80.0
         END IF
         FF = K*EXP(ARG)
C
         TF = TF/FF
C
 20   CONTINUE
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Possible overflow in calculating film factor, ',
     +        'K*exp(arg)')
 6020 FORMAT (' with arg =',G11.4,' Truncating arg to 80.00')
      END
**==lorentz.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       LORENTZ       **
C*************************
C
C
      REAL FUNCTION  LORENTZ (THETA)
C
C Purpose: Apply Lorentz correction
C
      INCLUDE 'implicit_none.finc'
      REAL THETA
C
C====== LOCALS:
C
      REAL  STH2
C
C------------------------------------------------------------------------------
C
      STH2 = SIN(THETA)**2
C
C====== lorentz = sth2       : changed by Pella 18.10.85
C
      LORENTZ = 1./STH2
      RETURN
      END
**==open_ge.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**       OPEN_GE        **
C**************************
C
C
      SUBROUTINE  OPEN_GE (IN_DATA1, IN_DATA2, FILE_NAME, END_POS,
     +                     END_REC)
C
C Purpose: To open .ge files using the given units and filename
C          Note that if the filename 'interactive' is supplied
C          then the user is prompted for a filename
C          Also note that filenames can have extensions, but these are removed
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER      IN_DATA1, IN_DATA2, END_REC
      CHARACTER*(*)  FILE_NAME
      INTEGER        END_POS
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      LOGICAL CCPEXS
C
C====== LOCALS:
C
      CHARACTER*80 STR
      INTEGER     I,IFAIL
      CHARACTER*60  FNAM, FNAME, FNAME1, FNAME2
      LOGICAL       EXISTS1, EXISTS2, QUIT
C
C FNAME : filename of .ge input file
C FNAME1 : filename with .ge1 appended
C FNAME2 : filename with .ge2 appended
C
C--------------------------------------------------------------------------
C
 20   IF (FILE_NAME.EQ.'interactive') THEN
         CALL TPROMP('Enter root name for .ge1/.ge2 files: ')
         CALL TREAD(FNAM,QUIT)
         CALL FILEXT(FNAM,'ge1',FNAME,I)
      ELSE
         CALL FILEXT(FILE_NAME,'ge1',FNAME,I)
      END IF
      FILE_NAME = FNAME
      END_POS = I
C
      FNAME1 = FNAME(1:I)//'.ge1'
      FNAME2 = FNAME(1:I)//'.ge2'
C
      EXISTS1 = CCPEXS(FNAME1)
      EXISTS2 = CCPEXS(FNAME2)
      IF (.NOT.(EXISTS1 .AND. EXISTS2)) THEN
         FILE_NAME = 'interactive'
         GO TO 20
      END IF
C
      IFAIL = 0
      CALL CCPDPN(-IN_DATA1,FNAME1,'READONLY','DU',19,IFAIL)
C
      IFAIL = 0
      CALL CCPDPN(-IN_DATA2,FNAME2,'READONLY','DU',2,IFAIL)
C
      WRITE (STR,6010) FNAME(1:I)
      CALL TWRITE(STR)
C
      READ (UNIT=IN_DATA1,REC=10) I , END_REC
      IF (END_REC.LT.1) THEN
         READ (UNIT=IN_DATA1,REC=9) I , I , END_REC
      END IF
      END_REC = END_REC + 10
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,'input ge file is ',A)
      END
**==output_unscram.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************************
C**       OUTPUT_UNSCRAM        **
C*********************************
C
C
      SUBROUTINE  OUTPUT_UNSCRAM (NUM_OVERLAPS, OVERLAP_INTENS,
     +      OVERLAP_LAMBDAS, OVERLAP_HKLS, OVERLAP_DHKLS,
     +      OVERLAP_SIGMAS,
     +      OVERLAP_XY,
     +      PARAMS, ACELL, BCELL, CCELL, CELLANG)
C
C Purpose: To output unscrambled reflections suitable as input for LAUESCAL
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  NUM_OVERLAPS
      REAL     OVERLAP_INTENS(6,NUM_OVERLAPS)
      REAL     OVERLAP_SIGMAS(6,NUM_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(NUM_OVERLAPS)
      REAL     OVERLAP_DHKLS(NUM_OVERLAPS)
      REAL     OVERLAP_XY(2,NUM_OVERLAPS)
      INTEGER  OVERLAP_HKLS(NUM_OVERLAPS)
C
      DOUBLE PRECISION     PARAMS(2,3,5)
      REAL     ACELL, BCELL, CCELL, CELLANG(3)
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS
C
      CHARACTER*80  STR
      CHARACTER*60  FNAM, OUT_FILE
C
      INTEGER  NWORDS, NVALS
      REAL     VALUES(10)
      CHARACTER*20  WORDS(10)
      INTEGER  CLASS, MULTI
      INTEGER  SCALE_SPOT
      LOGICAL    QUIT, MORE
      REAL     SCALE_FACTOR
      INTEGER  I, J, NS, IFAIL
      REAL     X, Y, INTENSITIES(6), SIGMAS(6)
      INTEGER  H, K, L
      INTEGER  I_INTENS, I_SIGMA
      REAL     LAMBDA, THETA, LOR, P, SCALED_INTENS, SCALED_SIGMA
C
      INTEGER  MAX_BAD_SPOTS, BAD_SPOT_COUNT
      PARAMETER (MAX_BAD_SPOTS=10)
C
      LOGICAL    APPLY_SIGMA, APPLY_LAMBDA, APPLY_MININT
      INTEGER  SIGMA_LIM, MININT
      REAL     LAMBDA_LLIM, LAMBDA_HLIM
C
C bad_spot_count: counts the number of times an integer bad_spot might occur
C                 and if this exceeds max_bad_spots then the routine aborts
C
      INTEGER       HKL(4)
C
C------------------------------------------------------------------------------
C
      BAD_SPOT_COUNT = 0
      APPLY_SIGMA = .false.
      APPLY_LAMBDA = .false.
      APPLY_MININT = .false.
      SIGMA_LIM = 0
      MININT = 0
      LAMBDA_LLIM = 0.0
      LAMBDA_HLIM = 999.0
C
C====== Open output file for LAUESCALE/LAUENORM
C
 20   CALL TPROMP('Enter name for output file (default ext=.afout): ')
      CALL TREAD(FNAM,QUIT)
      CALL FILEXT(FNAM,'afout',OUT_FILE,I)
      IFAIL = 1
      CALL CCPDPN(-4,OUT_FILE,'NEW','F',0,IFAIL)
      IF (IFAIL.EQ.-1) GO TO 20
      CALL TWRITE(' ')
      WRITE (STR,6020) OUT_FILE
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      CALL TWRITE(' ')
C
C====== Get class of spot to output
C
 30   CALL TPROMP('Enter class of spot to output (all,2,3 etc.): ')
      CALL TREAD(STR,QUIT)
      CALL STRTKN(STR,' ,',10,WORDS,NWORDS,VALUES,NVALS)
C
      IF (NVALS.GT.1) GO TO 30
      IF (NWORDS.GT.1) GO TO 30
      IF (NWORDS.EQ.NVALS) GO TO 30
C
      CLASS = -1
      IF (NWORDS.EQ.1) THEN
         IF (WORDS(1).EQ.'all' .OR. WORDS(1).EQ.'ALL') CLASS = 0
         IF (WORDS(1).EQ.'doublets' .OR. WORDS(1).EQ.'DOUBLETS')
     +       CLASS = 2
         IF (WORDS(1).EQ.'triplets' .OR. WORDS(1).EQ.'TRIPLETS')
     +       CLASS = 3
      ELSE
         CLASS = NINT(VALUES(1))
         IF (CLASS.LT.2 .OR. CLASS.GT.6) GO TO 30
      END IF
      IF (CLASS.EQ.-1) GO TO 30
      WRITE (STR,'(''Class ='',I2)') CLASS
      CALL TWRITE(STR)
C
C====== Get overall scaling factor
C
      SCALE_FACTOR = 100.
      CALL TASKF('Input scaling factor',SCALE_FACTOR,2,QUIT)
C
C====== Get rejection criteria
C
      CALL TASKQ('Do you want to apply a sigma rejection criterion',
     +           APPLY_SIGMA,QUIT)
      IF (APPLY_SIGMA) THEN
         CALL TASKI('Enter sigma limit (reject I<n*sigma)',SIGMA_LIM,
     +              QUIT)
      END IF
      CALL TASKQ(
     +    'Do you want to apply a minimum intensity rejection criterion'
     +    ,APPLY_MININT,QUIT)
      IF (APPLY_MININT) THEN
         CALL TASKI('Enter minimum intensity',MININT,QUIT)
      END IF
      CALL TASKQ('Do you want to apply a lambda rejection criterion',
     +           APPLY_LAMBDA,QUIT)
      IF (APPLY_LAMBDA) THEN
         CALL TASKF('Enter lower lambda limit',LAMBDA_LLIM,2,QUIT)
         CALL TASKF('Enter higher lambda limit',LAMBDA_HLIM,2,QUIT)
      END IF
C
      NS = 0
C
C====== Use CCP4 routine SETRSL to initialise storage for STHLSQ
C====== (STHLSQ is called inside scale_spot)
C
      CALL SETRSL(ACELL,BCELL,CCELL,CELLANG(1),CELLANG(2),CELLANG(3))
C
      DO 50 I = 1 , NUM_OVERLAPS
C
C====== Is this spot in the right output class ?
C
         CALL CCPBYI(HKL,OVERLAP_HKLS(I),-4)
         MULTI = HKL(4)
         MULTI = IABS(MULTI)
         IF (CLASS.EQ.0 .OR. CLASS.EQ.MULTI) THEN
C
            H = HKL(1)
            K = HKL(2)
            L = HKL(3)
C
            LAMBDA = OVERLAP_LAMBDAS(I)
            X = OVERLAP_XY(1,I)
            Y = OVERLAP_XY(2,I)
C
            DO 40 J = 1 , 6
               INTENSITIES(J) = OVERLAP_INTENS(J,I)
               SIGMAS(J) = OVERLAP_SIGMAS(J,I)
 40         CONTINUE
C
            IFAIL = SCALE_SPOT(PARAMS,SCALE_FACTOR,X,Y,INTENSITIES,
     +              SIGMAS,H,K,L,LAMBDA,THETA,SCALED_INTENS,
     +              SCALED_SIGMA,I_INTENS,I_SIGMA,LOR,P)
            IF (IFAIL.NE.0) THEN
C
C====== Bad spot - reject and maybe abort
C
               IF (IFAIL.NE.2) BAD_SPOT_COUNT = BAD_SPOT_COUNT + 1
               IF (BAD_SPOT_COUNT.GT.MAX_BAD_SPOTS) THEN
                  CALL TWRITE(' ')
                  WRITE (STR,6040) MAX_BAD_SPOTS
                  CALL TWRITE(STR)
                  WRITE (STR,6050)
                  CALL TWRITE(STR)
                  CALL TWRITE(' ')
                  CALL TWRITE(' ')
                  GO TO 60
               END IF
C
C====== Check the criteria
C
            ELSE IF (.NOT.(APPLY_SIGMA .AND. I_INTENS.LT.(SIGMA_LIM*
     +               I_SIGMA))) THEN
               IF (.NOT.(APPLY_MININT .AND. I_INTENS.LT.MININT)) THEN
                  IF (.NOT.(APPLY_LAMBDA .AND. (LAMBDA.LT.LAMBDA_LLIM
     +                .OR. LAMBDA.GT.LAMBDA_HLIM))) THEN
C
C====== Write out to output file
C
                     NS = NS + 1
                     WRITE (4,6010) H , K , L , LAMBDA , THETA , 
     +                              I_INTENS , I_SIGMA , X , Y
                  END IF
               END IF
            END IF
         END IF
C
 50   CONTINUE
C
 60   CALL TWRITE(' ')
      WRITE (STR,6030) NS
      CALL TWRITE(STR)
C
      CLOSE (4)
C
      MORE = .FALSE.
      CALL TASKQ('Do you want to output some more reflections',MORE,
     +            QUIT)       
      IF (MORE) THEN
         BAD_SPOT_COUNT = 0
         GO TO 20
      ELSE
C
         RETURN
      END IF
C
C====== Format statements
C
 6010 FORMAT (3I5,2F10.5,2I10,2F8.3)
 6020 FORMAT ('Writing output to ',A)
 6030 FORMAT ('Total number of output reflections',I6)
 6040 FORMAT ('More than',I4,' bad spots detected.')
 6050 FORMAT ('*****  ABORTING  *****')
      END
**==pair_up.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       PAIR_UP       **
C*************************
C
C
      SUBROUTINE  PAIR_UP (NUM_OVERLAPS,
     +            OVERLAP_INTENS, OVERLAP_HKLS, OVERLAP_LAMBDAS,
     +            MAX_PAIRS, NUM_PAIRS, SPOT1, SPOT2)
C
C Purpose: To pair up equivalent reflections
C          but reject
C                    1. mismatched wavelengths
C                    2. negative unscrambled intensities
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  NUM_OVERLAPS, MAX_PAIRS
      REAL     OVERLAP_INTENS(6,NUM_OVERLAPS)
      INTEGER  OVERLAP_HKLS(NUM_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(NUM_OVERLAPS)
      INTEGER  NUM_PAIRS
      INTEGER  SPOT1(MAX_PAIRS), SPOT2(MAX_PAIRS)
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      INTEGER    COMPARE_INDEX, I
      CHARACTER*1  REPLY
      CHARACTER*80 STR
      LOGICAL      QUIT
      REAL       LAMBDA_MIS_LIMIT, SEARCH_LAMBDA, LAMBDA_MIS
C
      INTEGER       H, K, L, MARKED, HKL(4)
      INTEGER       SEARCH_H, SEARCH_K, SEARCH_L
      INTEGER  OVERLAP_SPOT
C
      INTEGER  M2
C
      EQUIVALENCE (HKL(1),H)
      EQUIVALENCE (HKL(2),K)
      EQUIVALENCE (HKL(3),L)
      EQUIVALENCE (HKL(4),MARKED)
C
C------------------------------------------------------------------------------
C
C====== Which index to compare when pairing?
C
      COMPARE_INDEX = 1
      CALL TPROMP('Enter index for pairing up, h, k or l  [h]: ')
      CALL TREAD(REPLY,QUIT)
      IF (REPLY.EQ.'h' .OR. REPLY.EQ.'H') COMPARE_INDEX = 1
      IF (REPLY.EQ.'k' .OR. REPLY.EQ.'K') COMPARE_INDEX = 2
      IF (REPLY.EQ.'l' .OR. REPLY.EQ.'L') COMPARE_INDEX = 3
C
C====== Wavelength mismatch rejection limit
C
      LAMBDA_MIS_LIMIT = 0.05
      CALL TASKF('Enter wavelength mismatch limit',LAMBDA_MIS_LIMIT,2,
     +           QUIT)
C
      NUM_PAIRS = 0
      OVERLAP_SPOT = 0
C
C====== First unmark all spots
C
      DO 20 I = 1 , NUM_OVERLAPS
         CALL CCPBYI(HKL,OVERLAP_HKLS(I),-4)
         M2 = MARKED
         MARKED = IABS(M2)
         CALL CCPIBY(OVERLAP_HKLS(I),HKL,-4)
 20   CONTINUE
C
C====== Find first unmarked spot
C
 30   OVERLAP_SPOT = OVERLAP_SPOT + 1
      IF (OVERLAP_SPOT.GT.NUM_OVERLAPS) GO TO 70
      CALL CCPBYI(HKL,OVERLAP_HKLS(OVERLAP_SPOT),-4)
 40   IF (MARKED.LT.0) THEN
         OVERLAP_SPOT = OVERLAP_SPOT + 1
         IF (OVERLAP_SPOT.GT.NUM_OVERLAPS) GO TO 70
         CALL CCPBYI(HKL,OVERLAP_HKLS(OVERLAP_SPOT),-4)
         GO TO 40
      END IF
C
C====== Reject if intensity negative
C
      IF (OVERLAP_INTENS(1,OVERLAP_SPOT).GE.0.0) THEN
C
C====== Now look for partner
C
         CALL CCPBYI(HKL,OVERLAP_HKLS(OVERLAP_SPOT),-4)
         SEARCH_H = H
         SEARCH_K = K
         SEARCH_L = L
         IF (COMPARE_INDEX.EQ.1) SEARCH_H = -H
         IF (COMPARE_INDEX.EQ.2) SEARCH_K = -K
         IF (COMPARE_INDEX.EQ.3) SEARCH_L = -L
         SEARCH_LAMBDA = OVERLAP_LAMBDAS(OVERLAP_SPOT)
C
         DO 50 I = OVERLAP_SPOT + 1 , NUM_OVERLAPS
            CALL CCPBYI(HKL,OVERLAP_HKLS(I),-4)
C
C====== Reject pair if wavelength difference > limit
C
            LAMBDA_MIS = ABS(OVERLAP_LAMBDAS(I)-SEARCH_LAMBDA)
            IF (LAMBDA_MIS.LE.LAMBDA_MIS_LIMIT) THEN
C
C====== Reject negative intensities
C
               IF (OVERLAP_INTENS(1,I).GE.0.0) THEN
C
C====== Are the hkl's equivalent ?
C
                  IF ((H.EQ.SEARCH_H) .AND. (K.EQ.SEARCH_K) .AND. 
     +                (L.EQ.SEARCH_L) .AND. (MARKED.GE.0)) GO TO 60
               END IF
            END IF
C
C====== Fallen through loop so failed to find partner
C
 50      CONTINUE
      END IF
      GO TO 30
C
C====== Found partner, so add to list and mark spot
C
 60   NUM_PAIRS = NUM_PAIRS + 1
      IF (NUM_PAIRS.GT.MAX_PAIRS) THEN
         WRITE (STR,6010) MAX_PAIRS
         CALL TWRITE(STR)
         GO TO 70
      END IF
C
      SPOT1(NUM_PAIRS) = OVERLAP_SPOT
      SPOT2(NUM_PAIRS) = I
      CALL CCPBYI(HKL,OVERLAP_HKLS(I),-4)
      M2 = MARKED
      MARKED = -IABS(M2)
      CALL CCPIBY(OVERLAP_HKLS(I),HKL,-4)
      GO TO 30
C
C====== All pairs found
C
 70   RETURN
C
C====== Format statements
C
 6010 FORMAT ('Limit of',I6,' pairs reached. ',
     +        'Rsymms will be incomplete.')
      END
**==polar.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       POLAR         **
C*************************
C
C
      REAL FUNCTION  POLAR (TAU,AX,AY,THETA)
C
C Purpose: Apply polarization correction
C
      INCLUDE 'implicit_none.finc'
C
      REAL  TAU, AX, AY, THETA
C
C====== LOCALS:
C
      REAL  PSI, C2TH2, ALPHA, P0
C
C------------------------------------------------------------------------------
C
      IF (AY.EQ.0.) THEN
C
         POLAR = 1.0
C
      ELSE
C
         PSI = ATAN(AX/AY)
         C2TH2 = COS(2.0*THETA)**2
         P0 = (1.0+C2TH2)/2.0
         ALPHA = COS(2.0*PSI)*SIN(2*THETA)**2
         ALPHA = ALPHA/(1.0+C2TH2)
         POLAR = P0*(1.0-TAU*ALPHA)
C
      END IF
C
      RETURN
      END
**==read_ge_params.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************************
C**       READ_GE_PARAMS        **
C*********************************
C
C
      SUBROUTINE READ_GE_PARAMS (CP, WAVMIN, WAVMAX)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose:
C Read lambda limits and cell parameters from .GE1 file header. Expand
C cell parameters up to complete set (depends on xtal system) and convert
C to direct cell.
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      REAL WAVE, WAVMIN, WAVMAX
      REAL CP(6), BUFF (9), W_FS, BUFF2(3)
      REAL PARAMS(6)
      INTEGER XALNO, I
      INTEGER ISYLAT(2), ISL
      CHARACTER*8 FNAME
      CHARACTER*40 BUF
C
      READ (IN_DATA,REC=3) W_FS , BUFF , CP , ISL
      CALL CCPI2I(ISYLAT,ISL,2,.TRUE.,.FALSE.)
      XALNO = ISYLAT(1)
      READ (IN_DATA,REC=4) BUF , BUFF2 , WAVE , WAVMIN , WAVMAX
C
      IF (XALNO.EQ.1) THEN
C
C====== Triclinic
C
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.2) THEN
C
C====== Monoclinic (unique axis b)
C
         CP(5) = CP(4)
         CP(4) = 90.
         CP(6) = 90.
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.3) THEN
C
C====== Orthorhombic
C
         CP(4) = 90.
         CP(5) = 90.
         CP(6) = 90.
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.4) THEN
C
C====== Tetragonal
C
         CP(3) = CP(2)
         CP(2) = CP(1)
         CP(4) = 90.
         CP(5) = 90.
         CP(6) = 90.
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.5) THEN
C
C====== Hexagonal axes
C
         CP(3) = CP(2)
         CP(2) = CP(1)
         CP(4) = 90.
         CP(5) = 90.
         CP(6) = 60.
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.6) THEN
C
C====== Rhombohedral axes
C
         CP(4) = CP(2)
         CP(5) = CP(2)
         CP(6) = CP(2)
         CP(2) = CP(1)
         CP(3) = CP(1)
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE IF (XALNO.EQ.7) THEN
C
C====== Cubic
C
         CP(2) = CP(1)
         CP(3) = CP(1)
         CP(4) = 90.
         CP(5) = 90.
         CP(6) = 90.
         CALL RECCEL(PARAMS,CP,WAVE)
      ELSE
         STOP ' ERROR IN HEADER DATA '
      END IF
      DO 20 I = 1 , 6
         CP(I) = PARAMS(I)
 20   CONTINUE
      RETURN
      END
**==read_vc.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**       READ_VC        **
C**************************
C
C
      SUBROUTINE READ_VC(PARAMS)
C
C Purpose: To read in Victoreen parameters from unit 10
C
      INCLUDE 'implicit_none.finc'
C
      DOUBLE PRECISION PARAMS(2,3,5),AL150(2,3,5),AL200(2,3,5)
C
      COMMON /IO_STREAMS/IN_TERM,OUT_TERM,IN_DATA,OUT_DATA
      INTEGER IN_TERM,OUT_TERM,IN_DATA,OUT_DATA
C
C====== FUNCTION TYPES:
C
      LOGICAL CCPEXS
      INTEGER LENSTR
C
C====== LOCALS:
C
      CHARACTER*80 STR
      CHARACTER *60 FNAM, FILENAME
      LOGICAL EXISTS, QUIT
      INTEGER COMM
      INTEGER I,J,K,TYPE,IFAIL
      REAL   RTHICK
      DOUBLE PRECISION FILM(3),METAL(3),THICK
      DATA FILM/1.975,.682,.2856/
C FILM holds absorbancy coeffs. for film at 3 wavelength ranges
      DATA METAL/.0034415,.1198245,0/
C METAL holds absorb. coefs. for al + cu (absorb. per micron)
C
      DATA AL150/1,1.975,1,.682,1,.2856,1,2.4912,1,1.1982,1,.8018,1,
     +     2.4912,1,1.1982,1,.8018,1,2.4912,1,1.1982,1,.8018,1,7.1454,1,
     +     5.276,1,4.1812/
      DATA AL200/1,1.975,1,.682,1,.2856,1,2.6633,1,1.3703,1,.9739,1,
     +     2.6633,1,1.3703,1,.9739,1,2.6633,1,1.3703,1,.9739,1,7.1454,1,
     +     5.276,1,4.1812/
C
C------------------------------------------------------------------------------
C
      TYPE = 1
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      CALL TWRITE('Options available on loading film factors:-')
      CALL TWRITE(' ')
      CALL TWRITE('1 - AL150')
      CALL TWRITE('2 - AL200')
      CALL TWRITE('3 - enter any other .VC file')
      CALL TWRITE('4 - define your own parameters')
      CALL TWRITE(' ')
      COMM = 1
 20   CALL TASKI('Enter option',COMM,QUIT)
      IF (COMM.LT.1 .OR. COMM.GT.4) GO TO 20
C
C====== Now determine which option is to be used
C
      IF (COMM.EQ.1) THEN
         DO 50 K = 1 , 5
            DO 40 J = 1 , 3
               DO 30 I = 1 , 2
                  PARAMS(I,J,K) = AL150(I,J,K)
 30            CONTINUE
 40         CONTINUE
 50      CONTINUE
         RETURN
 
      END IF
 
      IF (COMM.EQ.2) THEN
         DO 80 K = 1 , 5
            DO 70 J = 1 , 3
               DO 60 I = 1 , 2
                  PARAMS(I,J,K) = AL200(I,J,K)
 60            CONTINUE
 70         CONTINUE
 80      CONTINUE
         RETURN
 
      END IF
 
      IF (COMM.EQ.4) THEN
C
C====== Material and thickness must be I/P for
C====== each gap. Values are calc'd + written
C
         DO 110 K = 1 , 5
            CALL TWRITE(' ')
            WRITE (STR,6020) K
            CALL TWRITE(STR)
 90         CALL TASKI('Enter foil type (1-Al 2-Cu 3-NONE)',TYPE,QUIT)
            IF (TYPE.LT.1 .OR. TYPE.GT.3) GO TO 90
            IF (TYPE.NE.3) THEN
               RTHICK = THICK
               CALL TASKF('Enter foil thickness in microns',RTHICK,2,
     +                    QUIT)
               THICK = RTHICK
 
            ELSE
               THICK = 0.
            END IF
 
            DO 100 J = 1 , 3
C
C====== For each wavelength range
C
               PARAMS(1,J,K) = 1
               PARAMS(2,J,K) = FILM(J) + METAL(TYPE)*THICK
 100        CONTINUE
 110     CONTINUE
      END IF
 
      IF (COMM.EQ.3) THEN
 120     CALL TPROMP('Name of Victoreen coeffs file (default ext=.vc): '
     +               )
         CALL TREAD(FNAM,QUIT)
         CALL FILEXT(FNAM,'vc',FILENAME,I)
         EXISTS = CCPEXS(FILENAME)
         IF (.NOT.EXISTS) GO TO 120
         IFAIL = 1
         CALL CCPDPN(-10,FILENAME,'READONLY','F',0,IFAIL)
         IF (IFAIL.EQ.-1) GO TO 120
         WRITE (STR,6010) FILENAME(1:LENSTR(FILENAME))
         CALL TWRITE(STR)
C
         DO 150 K = 1 , 5
            DO 140 J = 1 , 3
               DO 130 I = 1 , 2
                  READ (10,*) PARAMS(I,J,K)
 130           CONTINUE
 140        CONTINUE
 150     CONTINUE
C
         CLOSE (10)
      END IF
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Victoreen coeffs file is ',A)
 6020 FORMAT ('********** GAP ',I4,' ***************')
      END
**==scale_spot.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       SCALE_SPOT        **
C*****************************
C
C
      INTEGER  FUNCTION  SCALE_SPOT (PARAMS, SCALE_FACTOR,
     +                          X, Y, INTENSITIES, SIGMAS, H, K, L,
     +                          LAMBDA, THETA, SCALED_INTENS,
     +                          SCALED_SIGMA, I_INTENS,
     +                          I_SIGMA, LOR,P)
C
C Purpose: To scale one spot combining the intensities from 6 films
C
      INCLUDE 'implicit_none.finc'
C
      DOUBLE PRECISION     PARAMS(2,3,5)
      REAL     SCALE_FACTOR
C
      REAL     X, Y, INTENSITIES(6), SIGMAS(6)
      INTEGER  H, K, L
      REAL     LAMBDA
C
      REAL     THETA, SCALED_INTENS, SCALED_SIGMA, LOR, P
      INTEGER  I_INTENS, I_SIGMA
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      CHARACTER*80 STR
      REAL     STHLSQ
      REAL     LORENTZ, POLAR
      REAL     KP, T
      PARAMETER (KP=0.105, T=1.0)
C
      REAL     MEAN_INTENS, MEAN_SIG, SIG1, SIG2
C
      REAL     SUM_OF_WEIGHTS, WEIGHTED_SUM, WEIGHTS(6)
      REAL     SUM_OF_WEIGHTS2, WEIGHTED_RESID2
C
      REAL     FACTOR
C
      INTEGER  I, J
C
      REAL     DSTAR, DHKL, SIN_THETA
C
C------------------------------------------------------------------------------
C
      SCALE_SPOT = 0
C
C====== Check for overload
C
      DO 20 J = 1 , 6
         IF (ABS(INTENSITIES(J)).EQ.9999.0 .OR. ABS(SIGMAS(J))
     +       .EQ.9999.0) THEN
            INTENSITIES(J) = 0
            SIGMAS(J) = 1E6
         END IF
 20   CONTINUE
C
C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2
C====== NB This function uses implicit input from SETRSL
C
      DSTAR = 2*SQRT(STHLSQ(H,K,L))
      DHKL = 1.0/DSTAR
      SIN_THETA = LAMBDA*DSTAR/2.0
C
C====== Check that sin theta makes sense
C
      IF (ABS(SIN_THETA).GT.1.0) THEN
         CALL TWRITE(' ')
         WRITE (STR,6030) SIN_THETA
         CALL TWRITE(STR)
         WRITE (STR,6040) H , K , L
         CALL TWRITE(STR)
         CALL TWRITE(' ')
         SIN_THETA = SIGN(1.0,SIN_THETA)
      END IF
      THETA = ASIN(SIN_THETA)
C
      CALL SCALE_UP(PARAMS,INTENSITIES,SIGMAS,LAMBDA,THETA)
C
      SUM_OF_WEIGHTS = 0.0
      WEIGHTED_SUM = 0.0
C
      DO 30 J = 1 , 6
         IF (SIGMAS(J).GE.1E6) THEN
            WEIGHTS(J) = 0.0
         ELSE
            WEIGHTS(J) = 1.0/(SIGMAS(J)*SIGMAS(J))
         END IF
         SUM_OF_WEIGHTS = SUM_OF_WEIGHTS + WEIGHTS(J)
         WEIGHTED_SUM = WEIGHTED_SUM + WEIGHTS(J)*INTENSITIES(J)
 30   CONTINUE
C
      IF (SUM_OF_WEIGHTS.LT.1E-6) THEN
         SCALE_SPOT = 2
         GO TO 50
      END IF
C
      MEAN_INTENS = WEIGHTED_SUM/SUM_OF_WEIGHTS
      SCALED_INTENS = MEAN_INTENS
      SIG1 = SQRT(1.0/SUM_OF_WEIGHTS)
C
      WEIGHTED_RESID2 = 0.0
      SUM_OF_WEIGHTS2 = 0.0
C
      DO 40 J = 1 , 6
         WEIGHTED_RESID2 = WEIGHTED_RESID2 + 
     +                     ((INTENSITIES(J)-MEAN_INTENS)*WEIGHTS(J))**2
         SUM_OF_WEIGHTS2 = SUM_OF_WEIGHTS2 + WEIGHTS(J)**2
 40   CONTINUE
C
      SIG2 = SQRT(WEIGHTED_RESID2/SUM_OF_WEIGHTS2)
      IF (SIG2.GT.SIG1) THEN
         MEAN_SIG = SIG2
      ELSE
         MEAN_SIG = SIG1
      END IF
      SCALED_SIGMA = MEAN_SIG
C
      LOR = LORENTZ(THETA)
      P = POLAR(1.0,X,Y,THETA)
      FACTOR = SCALE_FACTOR*EXP(KP*T/COS(2.0*THETA))/(LOR*P)
C
      MEAN_INTENS = MEAN_INTENS*FACTOR
      MEAN_SIG = MEAN_SIG*FACTOR
C
      IF (MEAN_INTENS.GE.1E9 .OR. MEAN_SIG.GE.1E9) THEN
         CALL TWRITE(' ')
         WRITE (STR,6010) MEAN_INTENS
         CALL TWRITE(STR)
         WRITE (STR,6020)
         CALL TWRITE(STR)
         SCALE_SPOT = 1
      ELSE
         I_INTENS = MEAN_INTENS
         I_SIGMA = MEAN_SIG
      END IF
C
 50   RETURN
C
C====== Format statements
C
 6010 FORMAT ('Intensity =',G11.4,', danger of integer overflow. ',
     +        'Change scaling factor.')
 6020 FORMAT ('** Ignoring this spot **'/)
 6030 FORMAT ('Warning - lambda/2d (=sin theta) =',F5.1,
     +        ' -setting theta to 90')
 6040 FORMAT ('(h,k,l) = (',I4,',',I4,',',I4,')')
      END
**==scale_up.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       SCALE_UP        **
C***************************
C
C
      SUBROUTINE  SCALE_UP (PARAMS, INTENSITIES, SIGMAS, LAMBDA, THETA)
C
C Purpose: Scale up film intensities and sigmas to A film
C
      INCLUDE 'implicit_none.finc'
C
      DOUBLE PRECISION     PARAMS(2,3,5)
      REAL     INTENSITIES(6), SIGMAS(6)
      REAL     LAMBDA, THETA
C
C====== LOCALS
C
      DOUBLE PRECISION     DWAV, DC2TH, TF
      REAL     FF
      INTEGER  FILM
C
C------------------------------------------------------------------------------
C
      DWAV = LAMBDA
      DC2TH = COS(2.0*THETA)
C
C====== Calculate total scales
C
      DO 20 FILM = 2 , 6
C
C====== Unscram has reset bad sigmas to 1e6
C
         IF (SIGMAS(FILM).LT.1E6) THEN
            CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF)
            IF (TF.LT.1E-5) THEN
               FF = 0.0
            ELSE
               FF = 1.0/TF
            END IF
            INTENSITIES(FILM) = INTENSITIES(FILM)*FF
            SIGMAS(FILM) = SIGMAS(FILM)*FF
C
C====== Normal bad (-9999) sigmas caught here
C
            IF (SIGMAS(FILM).LE.0.0) SIGMAS(FILM) = 1E6
         END IF
 20   CONTINUE
C
      RETURN
C
      END
**==unscram_spots.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C********************************
C**       UNSCRAM_SPOTS        **
C********************************
C
C
      SUBROUTINE  UNSCRAM_SPOTS (MAX_OVERLAPS, END_REC,
     +                           ACELL, BCELL, CCELL, CELLANG,
     +                           ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX,
     +                           PARAMS,
     +                           NUM_OVERLAPS, OVERLAP_INTENS,
     +                           OVERLAP_SIGMAS, OVERLAP_LAMBDAS,
     +                           OVERLAP_HKLS, OVERLAP_DHKLS,
     +                           OVERLAP_XY,
     +                           SUPPRESS_OUTPUT)
C
C Purpose: To unscramble intensities for harmonics
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER  MAX_OVERLAPS, END_REC
C
      REAL     ACELL, BCELL, CCELL, CELLANG(3)
      REAL     ITOP, IBOT
      REAL     N_SIGS
      REAL     WAVMIN, WAVMAX
C
      DOUBLE PRECISION     PARAMS(2,3,5)
C
      INTEGER  NUM_OVERLAPS
      REAL     OVERLAP_INTENS(6,MAX_OVERLAPS)
      REAL     OVERLAP_SIGMAS(6,MAX_OVERLAPS)
      REAL     OVERLAP_LAMBDAS(MAX_OVERLAPS)
      REAL     OVERLAP_DHKLS(MAX_OVERLAPS)
      REAL OVERLAP_XY(2,MAX_OVERLAPS)
      INTEGER  OVERLAP_HKLS(MAX_OVERLAPS)
C
      LOGICAL    SUPPRESS_OUTPUT
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
      REAL     STHLSQ
      INTEGER  GET_OVERLAPS
      LOGICAL    QUIT
C
      INTEGER  NUM_SPOT_LOG,NUM_LOGGED
C
      INTEGER  OLD_OUT_TERM, RECORD_NUM, SPOT_NUM
C
      INTEGER  MAX_MULTI
      PARAMETER (MAX_MULTI=6)
C
      REAL     IIN(6), SIGIN(6), DELTAS(6), W_DELTAS(6)
      REAL     TF_MAT(6,MAX_MULTI)
      REAL     IOUT(MAX_MULTI), SIGOUT(MAX_MULTI), IFILM
C
      DOUBLE PRECISION     TF, DWAV, DC2TH
      INTEGER  FILM, ROW
      INTEGER  I, J, M, NUM_NEG
      LOGICAL    IS_NEG
C
      INTEGER       H, K, L, MULTI, HKL(4)
      REAL     DSTAR, DHKL
C
      REAL     X, Y, INTENSITIES(6), SIGMAS(6)
      INTEGER  NUM_GOOD_FILMS
      INTEGER  MULTIPLICITY, NODAL, SPATIAL_OVERLAP
      INTEGER  H_S(MAX_MULTI+1), K_S(MAX_MULTI+1), L_S(MAX_MULTI+1)
      REAL     LAMBDAS(MAX_MULTI+1), C2TH(MAX_MULTI+1)
      INTEGER    IFAIL
C
      EQUIVALENCE (HKL(1),H)
      EQUIVALENCE (HKL(2),K)
      EQUIVALENCE (HKL(3),L)
      EQUIVALENCE (HKL(4),MULTI)
C
C------------------------------------------------------------------------------
C
      NUM_LOGGED = 0
      NUM_SPOT_LOG = 20
      CALL TASKI('How many spots do you want to output on log file',
     +           NUM_SPOT_LOG,QUIT)
C
C====== Use CCP4 routine SETRSL to initialise storage for STHLSQ
C
      CALL SETRSL(ACELL,BCELL,CCELL,CELLANG(1),CELLANG(2),CELLANG(3))
C
C====== Loop through entire file, unscrambling on the way
C
      NUM_OVERLAPS = 0
      RECORD_NUM = 10
      NUM_NEG = 0
C
 20   IF (GET_OVERLAPS(IN_DATA,3,RECORD_NUM,END_REC,MAX_MULTI,ITOP
     +    ,IBOT,N_SIGS,WAVMIN,WAVMAX,X,Y,INTENSITIES,SIGMAS,
     +    NUM_GOOD_FILMS,MULTIPLICITY,NODAL,SPATIAL_OVERLAP,H_S,K_S,L_S,
     +    LAMBDAS,C2TH).EQ.0) THEN
C
C======  Write out spot info
C
         SPOT_NUM = RECORD_NUM - 10
         IF (.NOT.SUPPRESS_OUTPUT) THEN
            WRITE (OUT_TERM,6010) SPOT_NUM , X , Y , MULTIPLICITY , 
     +                            NODAL
            WRITE (OUT_TERM,6020) (INTENSITIES(I),I=1,6) , 
     +                            (SIGMAS(I),I=1,6)
            WRITE (OUT_TERM,6030)
         END IF
         IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
            WRITE (8,6010) SPOT_NUM , X , Y , MULTIPLICITY , NODAL
            WRITE (8,6020) (INTENSITIES(I),I=1,6) , (SIGMAS(I),I=1,6)
            WRITE (8,6030)
         END IF
         DO 30 I = 1 , MULTIPLICITY
            IF (.NOT.SUPPRESS_OUTPUT) THEN
               WRITE (OUT_TERM,6040) H_S(I) , K_S(I) , L_S(I) , 
     +                               LAMBDAS(I)
            END IF
            IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
               WRITE (8,6040) H_S(I) , K_S(I) , L_S(I) , LAMBDAS(I)
            END IF
 30      CONTINUE
C
C====== Put the valid intensities into input vector
C
         ROW = 0
         DO 50 FILM = 1 , 6
            IF (INTENSITIES(FILM).GE.0.0) THEN
               ROW = ROW + 1
               IIN(ROW) = INTENSITIES(FILM)
               SIGIN(ROW) = 1.0/SIGMAS(FILM)**2
               DO 40 M = 1 , MULTIPLICITY
                  DWAV = LAMBDAS(M)
                  DC2TH = C2TH(M)
                  CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF)
                  TF_MAT(ROW,M) = TF
 40            CONTINUE
            END IF
 50      CONTINUE
C
C====== Display transmission factor matrix
C
         CALL MTX_PACK(TF_MAT,TF_MAT,6,MAX_MULTI,ROW,MULTIPLICITY)
         IF (.NOT.SUPPRESS_OUTPUT) THEN
            WRITE (OUT_TERM,6090)
            CALL MTX_WRITE(OUT_TERM,TF_MAT,ROW,MULTIPLICITY)
         END IF
         IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
            WRITE (8,6090)
            CALL MTX_WRITE(8,TF_MAT,ROW,MULTIPLICITY)
         END IF
C
C====== Unscramble!
C
         CALL UNSORD(IIN,SIGIN,TF_MAT,IOUT,SIGOUT,NUM_GOOD_FILMS,
     +               MULTIPLICITY)
C
C====== Resulting unscrambled spot in iout array
C====== set reject flag if any negative components
C
         IS_NEG = .FALSE.
         DO 60 I = 1 , MULTIPLICITY
            IF (IOUT(I).LT.0) IS_NEG = .TRUE.
 60      CONTINUE
C
C====== Display results
C
         IF (.NOT.SUPPRESS_OUTPUT) THEN
            WRITE (OUT_TERM,6050)
         END IF
         IF (NUM_LOGGED.LE.NUM_SPOT_LOG) WRITE (8,6050)
C
         DO 70 FILM = 1 , 6
            IF (INTENSITIES(FILM).GE.0.0) THEN
               DELTAS(FILM) = -1.0*INTENSITIES(FILM)
            ELSE
               DELTAS(FILM) = -9999.0
            END IF
 70      CONTINUE
C
         DO 90 I = 1 , MULTIPLICITY
C
C====== Record intensities, lambdas, hkls and dhkls for R-symm calculation
C
            NUM_OVERLAPS = NUM_OVERLAPS + 1
            IF (NUM_OVERLAPS.GT.MAX_OVERLAPS) THEN
               IF (.NOT.SUPPRESS_OUTPUT) THEN
                  WRITE (OUT_TERM,*) ' Reached overlap limit of' , 
     +                               MAX_OVERLAPS
               END IF
               IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
                  WRITE (8,*) ' Reached overlap limit of' , MAX_OVERLAPS
               END IF
               GO TO 130
            END IF
C
            OVERLAP_LAMBDAS(NUM_OVERLAPS) = LAMBDAS(I)
C
            H = H_S(I)
            K = K_S(I)
            L = L_S(I)
            MULTI = MULTIPLICITY
            CALL CCPIBY(OVERLAP_HKLS(NUM_OVERLAPS),HKL,-4)
C
C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2
C====== NB This function uses implicit input from SETRSL
C
            DSTAR = 2*SQRT(STHLSQ(H_S(I),K_S(I),L_S(I)))
C
            DHKL = 1.0/DSTAR
            IF (DHKL.LE.0.0) THEN
               IF (.NOT.SUPPRESS_OUTPUT) THEN
                  WRITE (OUT_TERM,*) ' Dhkl = ' , DHKL , ' OV=' , 
     +                               NUM_OVERLAPS
               END IF
               WRITE (8,*) ' Dhkl = ' , DHKL , ' OV=' , NUM_OVERLAPS
            END IF
            OVERLAP_DHKLS(NUM_OVERLAPS) = DHKL
            OVERLAP_XY(1,NUM_OVERLAPS) = X
            OVERLAP_XY(2,NUM_OVERLAPS) = Y
C
            DWAV = LAMBDAS(I)
            DC2TH = C2TH(I)
            IFILM = IOUT(I)
            DO 80 FILM = 1 , 6
C
C====== Reject spots which got bad sigmas by setting them to 0
C====== -following programs will trap this.
C
               IF (SIGOUT(I).EQ.-9999.0) THEN
                  OVERLAP_INTENS(FILM,NUM_OVERLAPS) = 0.0
                  OVERLAP_SIGMAS(FILM,NUM_OVERLAPS) = 0.0
               ELSE
                  CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF)
                  OVERLAP_INTENS(FILM,NUM_OVERLAPS) = TF*IFILM
                  OVERLAP_SIGMAS(FILM,NUM_OVERLAPS) = TF*SIGOUT(I)
                  IF (DELTAS(FILM).NE.-9999.0) THEN
                     DELTAS(FILM) = DELTAS(FILM)
     +                              + OVERLAP_INTENS(FILM,NUM_OVERLAPS)
                  END IF
               END IF
 80         CONTINUE
            IF (.NOT.SUPPRESS_OUTPUT) THEN
               WRITE (OUT_TERM,6060) LAMBDAS(I) , 
     +                               (OVERLAP_INTENS(J,NUM_OVERLAPS),
     +                               J=1,6)
               WRITE (OUT_TERM,6070) (OVERLAP_SIGMAS(J,NUM_OVERLAPS),J=1
     +                               ,6)
            END IF
            IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
               WRITE (8,6060) LAMBDAS(I) , 
     +                        (OVERLAP_INTENS(J,NUM_OVERLAPS),J=1,6)
               WRITE (8,6070) (OVERLAP_SIGMAS(J,NUM_OVERLAPS),J=1,6)
            END IF
 90      CONTINUE
         IF (.NOT.SUPPRESS_OUTPUT) THEN
            WRITE (OUT_TERM,6100) (DELTAS(I),I=1,6)
         END IF
         IF (NUM_LOGGED.LE.NUM_SPOT_LOG) WRITE (8,6100)
     +       (DELTAS(I),I=1,6)
         DO 100 I = 1 , 6
            IF (DELTAS(I).NE.-9999.0) THEN
               W_DELTAS(I) = DELTAS(I)/(DELTAS(I)+2*INTENSITIES(I))
            ELSE
               W_DELTAS(I) = -9999.0
            END IF
 100     CONTINUE
         IF (.NOT.SUPPRESS_OUTPUT) THEN
            WRITE (OUT_TERM,6110) (W_DELTAS(I),I=1,6)
            WRITE (OUT_TERM,6080)
         END IF
         IF (NUM_LOGGED.LE.NUM_SPOT_LOG) THEN
            WRITE (8,6110) (W_DELTAS(I),I=1,6)
            WRITE (8,6080)
         END IF
C
C====== That was the last log file output for this spot so inc num_logged.
C
         NUM_LOGGED = NUM_LOGGED + 1
C
C======  Loop for more spots
C
C====== Reset counter if need to reject this harmonic
C====== because of negative components
C
         IF (IS_NEG) THEN
            NUM_OVERLAPS = NUM_OVERLAPS - MULTIPLICITY
            NUM_NEG = NUM_NEG + 1
         ELSE
C
C====== Only save top film measurement
C
            DO 120 J = 1 , MULTIPLICITY
               DO 110 I = 2 , 6
                  OVERLAP_SIGMAS(I,NUM_OVERLAPS-MULTIPLICITY+J) = -9999.
 110           CONTINUE
 120        CONTINUE
C
         END IF
         GO TO 20
      END IF
C
C======  Now reached end of .ge file or too many overlaps
C
 130  IF (.NOT.SUPPRESS_OUTPUT) THEN
         WRITE (OUT_TERM,6120) NUM_OVERLAPS
      END IF
      IF (NUM_LOGGED.LE.NUM_SPOT_LOG) WRITE (8,6120) NUM_OVERLAPS
C
      IF (.NOT.SUPPRESS_OUTPUT) THEN
         WRITE (OUT_TERM,6130) NUM_NEG
      END IF
      IF (NUM_LOGGED.LE.NUM_SPOT_LOG) WRITE (8,6130) NUM_NEG
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,130('-'),/,1X,'Spot ',I5,'    X = ',G11.4,' Y = ',
     +        G11.4,' Multiplicity = ',I2,'   Nodal flag = ',I1)
 6020 FORMAT (' Observed intensities on A..F films, and sigmas'/,1X,6X,
     +        1X,6G14.6,/,8X,6G14.6)
 6030 FORMAT (1X,'    h    k    l       lambda')
 6040 FORMAT (1X,3I5,5X,G11.4)
 6050 FORMAT (1X,'Wavelength,  Unscrambled intensities on A..F films, ',
     +        'and sigma')
 6060 FORMAT (1X,F6.3,',',6G14.6)
 6070 FORMAT (1X,'sigma= ',6G14.6)
 6080 FORMAT (1X,130('.'))
 6090 FORMAT (9X,'Transmission Factor Matrix')
 6100 FORMAT (1X,'Deltas=',6G14.6)
 6110 FORMAT (1X,'W. Del=',6G14.6)
 6120 FORMAT (/' Number of reflections read and unscrambled =',I6/)
 6130 FORMAT (/' Number of rejected reflections (negative component',1X,
     +        ')=',I6/)
      END
**==unsord.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       UNSORD        **
C*************************
C
C
      SUBROUTINE  UNSORD (Y, RECIP_SIG2, TF_MAT, X, INTENS_SIGMA, M, N)
C
C Purpose: Unscramble harmonics from an intensity value
C          We have n harmonics and m good films
C
       INCLUDE 'implicit_none.finc'
C
       INTEGER  M, N
       REAL     Y(M), TF_MAT(M,N), RECIP_SIG2(M)
       REAL     X(N), INTENS_SIGMA(N)
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C====== LOCALS:
C
       REAL     ROWCOL
C
       REAL     CX(6,6), CY(6,6), GY(6,6), CT(6,6), A(6,6)
       REAL     SCR1(6,6), SCR2(6,6)
       REAL     R, INTSIG
       INTEGER  LL, IND
C
C------------------------------------------------------------------------------
C
C====== Make weight matrix GY
C
      CALL MTX_DIAG(RECIP_SIG2,GY,M)
C
C
C====== Make reciprocal scale matrix in A
C
      CALL MTX_NEGATE(TF_MAT,A,M,N)
C
C
C====== Need to solve  A*X=Y , by least squares.
C====== Y known, X unknown intensities
C
      CALL MTX_MULT(GY,A,SCR1,M,M,N)
C
C
      CALL MTX_MULT_TR(A,SCR1,SCR2,N,M,N)
C
C
      CALL MTX_MULT_TR(A,GY,SCR1,N,M,M)
C
C
      CALL MTX_INV(SCR2,CX,N)
C
C
      CALL MTX_MULT(CX,SCR1,CT,N,N,M)
C
C
      CALL MTX_MULT(CT,Y,X,N,M,1)
      CALL MTX_SCA_MULT(X,X,-1.,N,1)
C
C
      DO 20 LL = 1 , N
         INTSIG = ROWCOL(CX,N,N,LL,LL)
C
C======Flag spots with bad sigmas (this happens v. occasionally)
C
         IF (INTSIG.GE.0) THEN
            INTENS_SIGMA(LL) = SQRT(INTSIG)
         ELSE
            INTENS_SIGMA(LL) = -9999.0
         END IF
 20   CONTINUE
C
      RETURN
      END
