**==aa0001.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**   Program PREPARE     **
C***************************
C
C Purpose: To prepare mtz file for harmonics
C
C Input: .ge files
C        ( .ge files, i.e. .ge1 and .ge2 are produced by INTLAUE)
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)
 
      LOGICAL    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,'PRELOG','NEW','F',0,IFAIL)
      CALL TWRITE(' ')
      CALL TWRITE('Program:     PREPARE')
      CALL TWRITE('Last update: 01/07/93')
      CALL TWRITE('Authors:     Q.Hao')
      CALL TWRITE(' ')
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,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====== Now prepare the spots
C
      CALL PREPARE_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)
C
      IF (NUM_OVERLAPS.LE.0) THEN
           STOP ' No overlaps found'                                                                                                  
      END IF                 
      STOP
      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
**==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
ccc      IF (NUM_GOOD_FILMS.LT.MULTIPLICITY) GO TO 20
      IF (NUM_GOOD_FILMS.LT.1) GO TO 20
      IF (MULTIPLICITY.GT.6) 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),WAV1,WAV2,R
      LOGICAL QUIT
      INTEGER I,J,NRANGE,ORDER
      CHARACTER*80 STR
      COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10)
C
C------------------------------------------------------------------------------
C
C====== Get unit cell parameters and wavelength range from .GE1 file.
C
      CALL READ_GE_PARAMS(CP,WAVMIN,WAVMAX)
      NRANGE = 2
      CALL TASKI('Enter Number of Wavelength ranges',NRANGE,QUIT)
      DO 10 I = 1,NRANGE
      READ(IN_TERM,*) WAV1(I),WAV2(I),ORDER(I),(R(I,J),J=1,ORDER(I)+1)
10    CONTINUE
      WAVMIN=WAV1(1)
      WAVMAX=WAV2(NRANGE)
      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
**==open_ge.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**       OPEN_GE        **
C**************************
C
C
      SUBROUTINE  OPEN_GE (IN_DATA1, IN_DATA2, END_POS,
     +                     END_REC)
C
C Purpose: To open .ge files using the given units and filename
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER      IN_DATA1, IN_DATA2, END_REC
      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
      LOGICAL       EXISTS1, EXISTS2, QUIT
C
C GE1 : filename with .ge1 appended
C GE2 : filename with .ge2 appended
C
C--------------------------------------------------------------------------
C
      IFAIL = 0
      CALL CCPDPN(-IN_DATA1,'GE1','READONLY','DU',19,IFAIL)
C
      IFAIL = 0
      CALL CCPDPN(-IN_DATA2,'GE2','READONLY','DU',2,IFAIL)
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
      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*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
      LOGICAL EXISTS, QUIT
      INTEGER COMM
      INTEGER I,J,K,TYPE,IFAIL
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('0 - single image plate')
      CALL TWRITE('1 - AL150')
      CALL TWRITE('2 - AL200')
      CALL TWRITE('3 - enter any other .VC file')
      CALL TWRITE(' ')
      COMM = 1
 20   CALL TASKI('Enter option',COMM,QUIT)
      IF (COMM.LT.0 .OR. COMM.GT.3) GO TO 20
C
C====== Now determine which option is to be used
C
      IF (COMM.EQ.0) THEN
         DO 25 K = 1 , 5
            DO 24 J = 1 , 3
                  PARAMS(1,J,K) = 1.
                  PARAMS(2,J,K) = 0.
 24         CONTINUE
 25      CONTINUE
         RETURN

      END IF

      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.3) THEN
         IFAIL = 1
         CALL CCPDPN(-10,'VC','READONLY','F',0,IFAIL)
         IF (IFAIL.EQ.-1) GO TO 220
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
220   STOP ' .vc file not found'
      END
**==sthlsq.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       STHLSQ        **
C*************************
C
C
      FUNCTION STHLSQ(IH, IK, IL )
C     ============================
C
      INTEGER IH,IK,IL
      COMMON /RECPLT/ COEFHH, COEFHK, COEFHL, COEFKK, COEFKL, COEFLL
C     CALCULATE (SIN(THETA)/LAMBDA)**2 FROM H,K,L; COEF'S SET BY CALL TO
C        SETRSL : GOOD FOR ANY KIND OF AXES
C
C-------------------------------------------------------------------------------
C
      STHLSQ = IH*IH*COEFHH + IH*IK*COEFHK + IH*IL*COEFHL + 
     +         IK*IK*COEFKK + IK*IL*COEFKL + IL*IL*COEFLL
      RETURN
      END
**==prepare_spots.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C********************************
C**       PREPARE_SPOTS         **
C********************************
C
C
      SUBROUTINE  PREPARE_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)
C
C Purpose: To prepare 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
C
      COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10)
      INTEGER NRANGE,ORDER
      REAL WAV1,WAV2,R
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER           IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      COMMON /SYMM1/PGNAME
      CHARACTER*10 PGNAME
      CHARACTER*10 SPGN
      CHARACTER*6 PGN 
      CHARACTER*1 LTP
      INTEGER NSY,NSG,ISYM
C
C====== LOCALS:
C
      REAL     STHLSQ
      INTEGER  GET_OVERLAPS
      LOGICAL    QUIT,MTZ
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,RESOL
C
      REAL X, Y, INTENSITIES(6), SIGMAS(6),rat(6),ratall
      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)
      REAL     sigtemp,lamcur,rattt
      INTEGER    IFAIL
C
CCC      INTEGER IPOINT(5) 
CCC      CHARACTER*30 LABEL(5)
CCC      CHARACTER*1 CTYP(5) 
CCC      LOGICAL EOF
      REAL ADATA(5)
C
      EQUIVALENCE (HKL(1),H)
      EQUIVALENCE (HKL(2),K)
      EQUIVALENCE (HKL(3),L)
      EQUIVALENCE (HKL(4),MULTI)
C
CCC      DATA LABEL/'H','K','L','F','SIGF'/
CCC      DATA CTYP/'H','H','H','F','Q'/  
C
C------------------------------------------------------------------------------
C
      NUM_LOGGED = 0
      NUM_SPOT_LOG = 0
      CALL TASKI('How many spots do you want to output on log file',
     +           NUM_SPOT_LOG,QUIT)
CCC      DO 75  I =1,5
CCC      IPOINT(I)=-1 
CCC75    CONTINUE
	CALL LROPEN (3,'HKLIN',1,IFAIL)
	IF (IFAIL.EQ.1) STOP 
        CALL LRSYMI(3,NSY,LTP,NSG,SPGN,PGN)
        PGNAME = PGN(3:)
         CALL SYMMAT 
CCC	CALL LRASSN (3,LABEL,5,IPOINT,CTYP)
        CALL LWOPEN (3,'HKLOUT')
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
C
      NUM_OVERLAPS = 0
      RECORD_NUM = 10
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 (NUM_LOGGED.LT.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 (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN
               WRITE (8,6040) H_S(I) , K_S(I) , L_S(I) , LAMBDAS(I)
            END IF
            NUM_OVERLAPS = NUM_OVERLAPS + 1
         IF (PGNAME(1:1).EQ.'3' .OR. PGNAME(1:1).EQ.'6'
     +       .OR. PGNAME(1:1).EQ.'4') THEN
            CALL TURN3(H_S(I),K_S(I),L_S(I),ISYM)
         ELSE
            CALL TURNIP(H_S(I),K_S(I),L_S(I),ISYM)
         END IF
         ADATA(1) = H_S(I)  
         ADATA(2) = K_S(I)  
         ADATA(3) = L_S(I)  
         ADATA(4) = NUM_OVERLAPS*10.
         ADATA(5) = SQRT(ADATA(4))
         CALL LWREFL (3,ADATA)
 30      CONTINUE
C
C
        NUM_LOGGED = NUM_LOGGED + 1
        GOTO 20
       END IF
      CALL LWCLOS(3,1) 
      WRITE (8,6120) NUM_OVERLAPS
C
      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)
 6120 FORMAT (/' Number of reflections read  =',I6/)
      END
