**==aa0001.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       AFSCALE       **
C*************************
C
C Program: AFSCALE
C
C Purpose: To fit film factor vs wavelength curve given data from INTLAUE
C
C
C     Fit curve to ratio of Ifilma to Ifilmb for laue films
C
C
      INCLUDE 'implicit_none.finc'
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION     PARAMS(2,3,5)
      COMMON /PARAMK/ K
      DOUBLE PRECISION K
      COMMON /RANGES/ WAVMIN,WAVMID1,WAVMID2,WAVMAX
      DOUBLE PRECISION WAVMIN,WAVMID1,WAVMID2,WAVMAX
C
      COMMON /BP/ BOXPROF
      INTEGER   BOXPROF
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C IN_TERM : main input stream (usually the terminal =5)
C OUT_TERM : main output (printer) stream (usually =6)
C IN_DATA : input stream for data (to be written to files)
C OUT_DATA : output stream for data (to be written to files)
C
C  I/O Unit Usage:
C  ---------------
C
C      1      -
C      2      .ge1 file
C      3      .ge2 file
C      4      output file for scaled data
C      5      terminal input
C      6      terminal output
C      7      -
C      8      log file
C      9      temp scratch file for singlets
C      10      film parameters file
C
C
C====== Set up arrays to replace unit 9 scratch files
C
      INTEGER   ARRAY_SIZE
      PARAMETER (ARRAY_SIZE = 40000)
      REAL   WAV(ARRAY_SIZE), FACTOR(ARRAY_SIZE), WEIGHT(ARRAY_SIZE)
      REAL   AVAL(ARRAY_SIZE), BVAL(ARRAY_SIZE), C2TH(ARRAY_SIZE)
      REAL   DHKL(ARRAY_SIZE), DMIN_THRESHOLD(ARRAY_SIZE)
      INTEGER   IH(ARRAY_SIZE), IK(ARRAY_SIZE), IL(ARRAY_SIZE)
C
C====== Array_size, and array names are to be given as args. to functions or
C====== routines that require them.
C
      INTEGER LUNSTI, LUNSTO
C
C====== LOCALS:
C
      LOGICAL      EOF, ERR
      CHARACTER*80 STR
      CHARACTER*20 TRM
      INTEGER      ITTYP
      INTEGER      FITTED, NSPOT, FILM, END_REC, END_POS
C
C FITTED : flag indicating that data has been fitted (1=fit made)
C NSPOT : total number of accepted singlets
C FILM : front film of current pair of films being used
C
      REAL         WVMIN, WVMAX, ACELL, BCELL, CCELL, CELLANG(3)
      REAL         IAMIN, IBMIN, ACON, BCON
      INTEGER      IREC
      CHARACTER*1  COMM, FILM_X, FILM_Y
      CHARACTER*6  FILMS
      CHARACTER*11 INTER
      LOGICAL      DONE_OP, DO_OP
      INTEGER      DATA, I, J, IFAIL
      DOUBLE PRECISION  ALL_PARAMS(2,3,5)
      DOUBLE PRECISION  WWMIN, WWMAX
C
C COMM : command letter
C DATA : data subroutine - sets data parameters and collects data
C
      DATA  FILMS /'ABCDEF'/
      DATA  DONE_OP/.FALSE./
      DATA  INTER/'interactive'/
C
C--------------------------------------------------------------------------
C
C====== Initialisations
C
      K = 1.0
      BOXPROF = 1
      IN_TERM = LUNSTI(0)
      OUT_TERM = LUNSTO(0)
      IN_DATA = 2
      OUT_DATA = 1
      WAVMIN = 0.0
      WAVMID1 = 0.49
      WAVMID2 = 0.92
      WAVMAX = 2.0
C
C====== Open log file
C
      IFAIL = 0
      CALL CCPDPN(-8,'AFLOG','NEW','F',0,IFAIL)
      CALL TRWSET(LUNSTI(0),LUNSTO(0),8,8,0)
C
      CALL TWRITE(' ')
      CALL TWRITE('Program:     AFSCALE')
      CALL TWRITE('Last update: 02/06/92')
      CALL TWRITE('Authors:     I.J.Clifton, S.Zurek, N. Gunn')
      CALL TWRITE(' ')
C
C====== Open .ge files
C
      CALL OPEN_GE(IN_DATA,3,INTER,END_POS,END_REC)
C
C====== Initialise parameters
C
      FITTED = 0
      FILM = 0
      WVMIN = 0.3
      WVMAX = 2.5
      WWMIN = WVMIN
      WWMAX = WVMAX
C
      CELLANG(1) = 90.
      CELLANG(2) = 90.
      CELLANG(3) = 90.
C
      IAMIN = 100.0
      IBMIN = 100.0
C
      ACON = 0.0
      BCON = 0.0
C
      NSPOT = 0
C
C====== Graphics initialise
C
      ITTYP = 0
      CALL UGTENV('TERM',TRM)
      IF (TRM(1:5).EQ.'xterm') ITTYP = 1
      CALL QGINIT(15,OUT_TERM,ITTYP,LUNSTI(0),LUNSTO(0))
C
C====== Menu
C
      COMM = 'D'
C
C====== Read in singlets from .ge file
C
 20   IF (COMM.EQ.'d' .OR. COMM.EQ.'D') THEN
         IF (DATA(WVMIN,WVMAX,IAMIN,IBMIN,ACON,BCON,ACELL,BCELL,CCELL,
     +       CELLANG,NSPOT,FILM,ARRAY_SIZE,WAV,FACTOR,WEIGHT,AVAL,BVAL,
     +       IH,IK,IL,C2TH,DHKL,DMIN_THRESHOLD,FITTED).EQ.1)
     +       THEN
             WRITE (STR,6070)
             CALL TWRITE(STR)
         END IF
         CALL TWRITE(' ')
         WWMIN = WVMIN
         WWMAX = WVMAX
      END IF
C
C====== Fit option - note if no parameters have been read in
C====== then we force their reading on the user.
C====== This means we never have some parameters determined
C====== and others zero.
C
      IF (COMM.EQ.'f' .OR. COMM.EQ.'F') THEN
         IF (NSPOT.GT.0) CALL FIT(NSPOT,FILM,ARRAY_SIZE,WAV,FACTOR,
     +                            WEIGHT,AVAL,BVAL,IH,IK,IL,C2TH,DHKL,
     +                            DMIN_THRESHOLD)
      END IF
C
      IF (COMM.EQ.'s' .OR. COMM.EQ.'S')
     +    CALL SCALE(NSPOT,FILM,ARRAY_SIZE,WAV,FACTOR,WEIGHT,AVAL,BVAL,
     +    IH,IK,IL,C2TH,DHKL,DMIN_THRESHOLD)
C
      IF (COMM.EQ.'o' .OR. COMM.EQ.'O') THEN
         CALL OUTPUT(ACELL,BCELL,CCELL,CELLANG)
         DONE_OP = .TRUE.
      END IF
C
      IF (COMM.EQ.'w' .OR. COMM.EQ.'W') THEN
         CALL WRITE_FILM_FACTORS
      END IF
C
C====== All the following commands require a film pair to be read in
C
      IF (FILM.EQ.0) THEN
         CALL TWRITE(' ')
         WRITE (STR,6030)
         CALL TWRITE(STR)
         CALL TWRITE(' ')
C
C====== Plot points and also curve if fitted
C
      END IF
C
 30   IF (FILM.NE.0) THEN
         FILM_X = FILMS(FILM:FILM)
         FILM_Y = FILMS(FILM+1:FILM+1)
         CALL TWRITE(' ')
         WRITE (STR,6010) FILM , FILM_X , FILM_Y
         CALL TWRITE(STR)
         CALL TLOG(.FALSE.)
         CALL TWRITE('OPTIONS AVAILABLE-')
         CALL TWRITE(' ')
         CALL TWRITE('    E - Exit')
         CALL TWRITE('    D - Define selection criteria')
         CALL TWRITE('    F - Fit Victoreen')
         CALL TWRITE('    P - Plot scaling v. lambda')
         CALL TWRITE('    S - Calculate Scaling R-factor')
         CALL TWRITE('    R - Calculate symmetry R-factors')
         CALL TWRITE('    O - Output scaled-together intensities')
         CALL TWRITE('    W - Write film factors')
         CALL TWRITE('    L - Load film factors')
         CALL TWRITE(' ')
         CALL TLOG(.TRUE.)
         CALL TPROMP('Enter option: ')
      ELSE
         CALL TWRITE(' ')
         WRITE (STR,6020)
         CALL TWRITE(STR)
         CALL TLOG(.FALSE.)
         CALL TWRITE('OPTIONS AVAILABLE-')
         CALL TWRITE(' ')
         CALL TWRITE('    E - Exit')
         CALL TWRITE('    D - Define selection criteria')
         CALL TWRITE('    F - Fit Victoreen')
         CALL TWRITE('    P - Plot scaling v. lambda')
         CALL TWRITE('    S - Calculate Scaling R-factor')
         CALL TWRITE('    R - Calculate symmetry R-factors')
         CALL TWRITE('    O - Output scaled-together intensities')
         CALL TWRITE('    W - Write film factors')
         CALL TWRITE('    L - Load film factors')
         CALL TWRITE(' ')
         CALL TLOG(.TRUE.)
         CALL TPROMP('Enter option: ')
      END IF
C
      CALL TREAD(COMM,ERR)
C
C====== Interpret command
C
      IF (COMM.EQ.'e' .OR. COMM.EQ.'E') THEN
         IF (.NOT.DONE_OP) THEN
            CALL TASKQ('Do you want to write an output file',DO_OP,EOF)
            IF (DO_OP) CALL OUTPUT(ACELL,BCELL,CCELL,CELLANG)
         END IF
         CALL QGCLOS
         STOP
      END IF
C
C====== Load film factors - note only change those of current film if fitted
C
      IF (COMM.EQ.'l' .OR. COMM.EQ.'L') THEN
         IF (FITTED.EQ.1) THEN
            WRITE (STR,6050) FILM_X , FILM_Y
            CALL TWRITE(STR)
            CALL READ_VC(ALL_PARAMS)
            DO 50 I = 1 , 2
               DO 40 J = 1 , 3
                  PARAMS(I,J,FILM) = ALL_PARAMS(I,J,FILM)
 40            CONTINUE
 50         CONTINUE
         ELSE
            WRITE (STR,6060)
            CALL TWRITE(STR)
            CALL READ_VC(PARAMS)
         END IF
         FITTED = 1
      END IF
C
      IF (COMM.EQ.'r' .OR. COMM.EQ.'R')
     +    CALL RSYMM(NSPOT,ARRAY_SIZE,WAV,FACTOR,WEIGHT,AVAL,BVAL,IH,IK,
     +    IL,C2TH,DHKL,DMIN_THRESHOLD)
C
      IF (COMM.EQ.'p' .OR. COMM.EQ.'P') THEN
         CALL PLOT(FITTED,NSPOT,FILM,WWMIN,WWMAX,ARRAY_SIZE,WAV,FACTOR)
         CALL QGALPH
      END IF
C
C====== Note - all the following commands require
C====== the parameters to be determined either by reading in
C====== or by fitting.
C
      IF (FITTED.EQ.0) THEN
         CALL TWRITE(' ')
         WRITE (STR,6040)
         CALL TWRITE(STR)
         CALL TWRITE(' ')
         GO TO 30
      END IF
      GO TO 20
C
C====== Format statements
C
 6010 FORMAT ('Current film:',I2,' i.e. scaling from ',A1,' to ',A1)
 6020 FORMAT ('Currently no selected film.')
 6030 FORMAT ('!!! No films have been read yet. Use D(ata) option.')
 6040 FORMAT ('!!! No starting parameters have been loaded yet. ',
     +        'Please load them.')
 6050 FORMAT ('Loading new parameters for scaling between films ',A1,
     +        ' to ',A1,' only.')
 6060 FORMAT ('Loading new parameters for ALL films.')
 6070 FORMAT ('!!! Error whilst collecting data.')
      END
**==calfun.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       CALFUN        **
C*************************
C
C
      SUBROUTINE  CALFUN (M, N, FVECC, A)
C
C Purpose: HARWELL routine uses this to fit victoreen
C
C Updates: Changed function minimised to Rscale-like function
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    M, N
      DOUBLE PRECISION  A, FVECC(M)
C
      INTEGER    NUMSPOTMAX
      PARAMETER (NUMSPOTMAX = 12000)
C
      COMMON /CALFUN_PARAMS/ FUNCT_TYPE, POWER, ARRAY
      INTEGER                FUNCT_TYPE
      REAL                   POWER
      DOUBLE PRECISION       ARRAY(5,NUMSPOTMAX)
C
      COMMON /PARAMK/ K
      DOUBLE PRECISION K
C
C FUNCT_TYPE : determines how to interpret array
C            : if funct_type = 1 then
C ARRAY(1,I) : wavelength
C ARRAY(2,I) : film to film ratio
C ARRAY(3,I) : 1/sigma2
C            : if funct_type = 2 then
C ARRAY(1,I) : wavelength
C ARRAY(2,I) : aval - intensity value on X film
C ARRAY(3,I) : bval - intensity value on Y film
C ARRAY(4,I) : c2th
C ARRAY(5,I) : weighting -- could be 1/sigma2 or 1/aval
C
C====== LOCALS:
C
      DOUBLE PRECISION     VICTOREEN, DWAV
      DOUBLE PRECISION     AVAL, BVAL, DIFF, C2TH
      INTEGER    I
C
C------------------------------------------------------------------------------
C
      DO 20 I = 1 , M
         DWAV = ARRAY(1,I)
         IF (FUNCT_TYPE.EQ.1) THEN
            FVECC(I) = (ARRAY(2,I)-VICTOREEN(K,A,DWAV))*ARRAY(3,I)
         ELSE IF (FUNCT_TYPE.EQ.2) THEN
            AVAL = ARRAY(2,I)
            BVAL = ARRAY(3,I)
            C2TH = ARRAY(4,I)
            DIFF = AVAL - BVAL*VICTOREEN(K,A,DWAV)**(1.0/C2TH)
            FVECC(I) = ABS(ARRAY(5,I)*DIFF)**(POWER/2.0)
            IF (FVECC(I).GT.1E17) FVECC(I) = 1E17
         END IF
 20   CONTINUE
      RETURN
      END
**==data.f  
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       DATA          **
C*************************
C
C
C Purpose: Set data parameters
C
      INTEGER FUNCTION DATA(WAVMIN,WAVMAX,IAMIN,IBMIN,ACON,BCON,ACELL,
     +                      BCELL,CCELL,CELLANG,NSPOT,FILM,ARRAY_SIZE,
     +                      WWAV,FFACTOR,WWEIGHT,AAVAL,BBVAL,IIH,IIK,
     +                      IIL,CC2TH,DDHKL,DDMIN_THRESHOLD,FITTED)
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER NSPOT, FILM, FITTED
C
C NSPOT : number of accepted singlets
C FILM : front film of current pair of films being used
C
      REAL WAVMIN, WAVMAX, ACELL, BCELL, CCELL, CELLANG(3)
      REAL IAMIN, IBMIN, ACON, BCON
C
      INTEGER NOVMAX
      PARAMETER (NOVMAX=2000)
C
      COMMON /FITTING_PARAMS/PARAMS
      DOUBLE PRECISION PARAMS(2,3,5)
C
      COMMON /BP/BOXPROF
      INTEGER BOXPROF
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
      INTEGER GET_USERS_PARAMETERS, GET_GE_DATA
      REAL STHLSQ
      DOUBLE PRECISION EVAL_VIC
      REAL WEIGHTFUN
      INTEGER ARRAY_SIZE, IPACK, NFILM, ITSPOT, NSINGLETS
      INTEGER I, IIH(ARRAY_SIZE), IH, IIK(ARRAY_SIZE)
      INTEGER RANGE, IK, IIL(ARRAY_SIZE), IL, RECORD_NUM, END_REC
      REAL WMIN, WMAX, WAV, WWAV(ARRAY_SIZE)
      REAL WWEIGHT(ARRAY_SIZE), WEIGHT, DSTAR
      REAL DDHKL(ARRAY_SIZE), DHKL, THETA, CC2TH(ARRAY_SIZE), C2TH
      REAL FFACTOR(ARRAY_SIZE), FACTOR, THFACTOR, SIN_THETA
      REAL AAVAL(ARRAY_SIZE), AVAL, BBVAL(ARRAY_SIZE), BVAL
      REAL ASIG2, BSIG2, SIG2, ASIG, BSIG, AVAL2, BVAL2
      REAL X, Y, INTENSITIES(6), SIGMAS(6)
      REAL INTENB(6), SIGB(6), INTENP(6), SIGP(6)
      INTEGER MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP
      INTEGER H_S(6), K_S(6), L_S(6)
      REAL LAMBDAS(6)
      REAL DDMIN_THRESHOLD(ARRAY_SIZE), DMIN_THRESHOLD
      DOUBLE PRECISION P, DWAV
      REAL BAND
      CHARACTER*1 CHAR
      LOGICAL LBADA, LBADB, LENT, LBAND, SPATREJ
      INTEGER VERNO
      CHARACTER*72 SPACER72
      DATA LBADA, LBADB, LENT/.TRUE., .TRUE., .TRUE./
      DATA BAND/0.50/
C
C -----------------------------------------------------------------------------
C
      DATA = 0
C
      IF (GET_USERS_PARAMETERS(BOXPROF,SPATREJ,WAVMIN,WAVMAX,IAMIN,
     +    IBMIN,ACON,BCON,ACELL,BCELL,CCELL,CELLANG,LBADA,LBADB,LENT,
     +    LBAND,BAND,FILM).NE.0) THEN
         DATA = 1
         RETURN
 
      END IF
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====== Need to have read in the film factors if using "band" option.
C
      IF (LBAND .AND. (FITTED.EQ.0)) THEN
         CALL READ_VC(PARAMS)
         FITTED = 1
      END IF
C
C====== We now read in the .ge file from unit in_data and unit 3
C====== and put only singlets into arrays.
C
      READ (IN_DATA,REC=9) SPACER72 , VERNO
      IF (VERNO.NE.2) THEN
         WRITE (STR,6030) VERNO
         CALL TWRITE(STR)
         STOP ' AFSCALE aborting'
 
      END IF
      READ (IN_DATA,REC=10) IPACK , ITSPOT
      END_REC = 10 + ITSPOT
      NSINGLETS = 0
      NSPOT = 0
      WRITE (STR,6020) ITSPOT
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
      DO 40 RECORD_NUM = 11 , END_REC
 
         MULTIPLICITY = 6
         IF (GET_GE_DATA(IN_DATA,3,RECORD_NUM,X,Y,INTENB,SIGB,INTENP,
     +       SIGP,MULTIPLICITY,NODAL,NODIDX,SPATIAL_OVERLAP,H_S,K_S,L_S,
     +       LAMBDAS,DMIN_THRESHOLD).NE.0) THEN
            DATA = 1
            RETURN
 
         END IF
 
         IF (BOXPROF.EQ.1) THEN
            DO 20 I = 1 , 6
               INTENSITIES(I) = INTENB(I)
               SIGMAS(I) = SIGB(I)
 20         CONTINUE
         ELSE
            DO 30 I = 1 , 6
               INTENSITIES(I) = INTENP(I)
               SIGMAS(I) = SIGP(I)
 30         CONTINUE
         END IF
         IF (.NOT.(SPATREJ .AND. SPATIAL_OVERLAP.EQ.1)) THEN
C
            IF (MULTIPLICITY.EQ.1) THEN
               NSINGLETS = NSINGLETS + 1
               AVAL = INTENSITIES(FILM)
               ASIG = SIGMAS(FILM)
               BVAL = INTENSITIES(FILM+1)
               BSIG = SIGMAS(FILM+1)
C
               IH = H_S(1)
               IK = K_S(1)
               IL = L_S(1)
C
               WAV = LAMBDAS(1)
C
C====== Are we rejecting unmeasurable reflexions ?
C====== If so then start rejecting!
C====== We reject intensities lower than a minimum value
C====== those values with negative or large sigmas or values
C====== outside the sigma limits.
C
               IF (LBADA .AND. (AVAL.LE.IAMIN .OR. ASIG.LE.0.0))
     +             GO TO 40
               IF (LBADB .AND. (BVAL.LE.IBMIN .OR. BSIG.LE.0.0))
     +             GO TO 40
C
               IF (LBADA .AND. 
     +             (ABS(AVAL).GT.9000.0 .OR. ABS(ASIG).GT.9000.0))
     +             GO TO 40
C
               IF (LBADB .AND. 
     +             (ABS(BVAL).GT.9000.0 .OR. ABS(BSIG).GT.9000.0))
     +             GO TO 40
C
               IF (ABS(AVAL).LT.9000.0 .AND. AVAL.LE.ACON*ASIG) GO TO 40
               IF (ABS(BVAL).LT.9000.0 .AND. BVAL.LE.BCON*BSIG) GO TO 40
C
C====== Fiddle low intensity data
C
               IF (.NOT.LBADA .AND. AVAL.LE.IAMIN) AVAL = 9999.0
               IF (.NOT.LBADB .AND. BVAL.LE.IBMIN) BVAL = 9999.0
C
C====== Check wavelength is within range
C
               IF (WAV.LT.WAVMIN .OR. WAV.GT.WAVMAX) GO TO 40
               IF ((WAV.GE.0.486 .AND. WAV.LE.0.490) .OR. 
     +             (WAV.GE.0.900 .AND. WAV.LE.0.920)) GO TO 40
C
               AVAL2 = AVAL**2
               ASIG2 = ASIG**2
               BVAL2 = BVAL**2
               BSIG2 = BSIG**2
C
               SIG2 = WEIGHTFUN(AVAL2,BVAL2,ASIG2,BSIG2)
C
               IF (LENT) THEN
                  WEIGHT = 1.0
               ELSE
                  WEIGHT = SIG2
               END IF
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(IH,IK,IL))
               DHKL = 1.0/DSTAR
               SIN_THETA = WAV*DSTAR/2.0
               THETA = ASIN(SIN_THETA)
               C2TH = COS(2.0*THETA)
               FACTOR = (AVAL/BVAL)**C2TH
C
               IF (LBAND) THEN
                  DWAV = WAV
                  P = EVAL_VIC(DWAV,RANGE,FILM)
                  THFACTOR = P
                  IF (ABS(FACTOR-THFACTOR)/THFACTOR.GT.BAND) GO TO 40
               END IF
C
               NSPOT = NSPOT + 1
               IF (NSPOT.GT.ARRAY_SIZE)
     +             STOP ' Too many singlets accepted'
               WWAV(NSPOT) = WAV
               FFACTOR(NSPOT) = FACTOR
               WWEIGHT(NSPOT) = WEIGHT
               AAVAL(NSPOT) = AVAL
               BBVAL(NSPOT) = BVAL
               IIH(NSPOT) = IH
               IIK(NSPOT) = IK
               IIL(NSPOT) = IL
               CC2TH(NSPOT) = C2TH
               DDHKL(NSPOT) = DHKL
               DDMIN_THRESHOLD(NSPOT) = DMIN_THRESHOLD
C
            END IF
         END IF
C
 40   CONTINUE
      CALL TWRITE(' ')
      WRITE (STR,6010) NSPOT , NSINGLETS
      CALL TWRITE(STR)
C
      RETURN
C
C====== Format Statements
C
 6010 FORMAT ('Total number of accepted singlets =',I6,
     +        ' out of a possible ',I6,' singlets.')
 6020 FORMAT ('===Reading',I6,' spots...')
 6030 FORMAT ('!!! The GE file version no. is ',I2,' but this is a',
     +        ' Version 2 program')
      END
**==eval_vic.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       EVAL          **
C*************************
C
C
      DOUBLE PRECISION FUNCTION  EVAL_VIC (LAMBDA, RANGE, FILM)
C
C Purpose: To apply victoreen at given lambda
C
      INCLUDE 'implicit_none.finc'
      DOUBLE PRECISION     LAMBDA
      INTEGER    RANGE, FILM
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION     PARAMS(2,3,5)
      COMMON /PARAMK/K
      DOUBLE PRECISION K
C
      COMMON /RANGES/ WAVMIN,WAVMID1,WAVMID2,WAVMAX
      DOUBLE PRECISION WAVMIN,WAVMID1,WAVMID2,WAVMAX
C
C====== LOCALS:
C
      DOUBLE PRECISION  VICTOREEN, A
      INTEGER    FIT_RANGE
C
C------------------------------------------------------------------------------
C
      IF (LAMBDA.LT.WAVMIN) THEN
         RANGE = 0
         FIT_RANGE = 1
      ELSE IF (LAMBDA.GT.WAVMAX) THEN
         RANGE = 4
         FIT_RANGE = 3
      ELSE IF (LAMBDA.LT.WAVMID1) THEN
         RANGE = 1
         FIT_RANGE = 1
      ELSE IF (LAMBDA.GE.WAVMID1 .AND. LAMBDA.LT.WAVMID2) THEN
         RANGE = 2
         FIT_RANGE = 2
      ELSE IF (LAMBDA.GE.WAVMID2) THEN
         RANGE = 3
         FIT_RANGE = 3
      END IF
C
C====== K = params(1,fit_range,film)
C
      A = PARAMS(2,FIT_RANGE,FILM)
C
      EVAL_VIC = VICTOREEN(K,A,LAMBDA)
C
      RETURN
      END
**==fit.f  
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE FIT(NSPOT,FILM,ARRAY_SIZE,WAV,RAT,SIG2,AVAL,BVAL,IH,
     +               IK,IL,C2TH,DHKL,DMIN_THRESHOLD)
C
C Purpose: Fit curve to ratio of i (filma) to i(filmb) for laue films
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER NSPOT, FILM, ARRAY_SIZE
C
      COMMON /FITTING_PARAMS/PARAMS
      DOUBLE PRECISION PARAMS(2,3,5)
C
      COMMON /RANGES/WAVMIN, WAVMID1, WAVMID2, WAVMAX
      DOUBLE PRECISION WAVMIN, WAVMID1, WAVMID2, WAVMAX
C
      INTEGER NUMSPOTMAX
      PARAMETER (NUMSPOTMAX=8000)
C
      COMMON /CALFUN_PARAMS/FUNCT_TYPE, POWER, ARRAY
      INTEGER FUNCT_TYPE
      REAL POWER
      DOUBLE PRECISION ARRAY(5,NUMSPOTMAX)
C
C FUNCT_TYPE : determines how to interpret array
C            : if funct_type = 1 then
C ARRAY(1,I) : wavelength
C ARRAY(2,I) : film to film ratio
C ARRAY(3,I) : 1/sigma2
C            : if funct_type = 2 then
C ARRAY(1,I) : wavelength
C ARRAY(2,I) : aval - intensity value on X film
C ARRAY(3,I) : bval - intensity value on Y film
C ARRAY(4,I) : c2th
C ARRAY(5,I) : weighting -- could be 1/sigma2 or 1/aval
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
      INTEGER NUMR1, NUMR2, NUMR3, I
      REAL WAV(ARRAY_SIZE), RAT(ARRAY_SIZE)
      REAL SIG2(ARRAY_SIZE)
      REAL AVAL(ARRAY_SIZE), BVAL(ARRAY_SIZE)
      INTEGER IH(ARRAY_SIZE), IK(ARRAY_SIZE)
      INTEGER IL(ARRAY_SIZE)
      REAL C2TH(ARRAY_SIZE), DHKL(ARRAY_SIZE)
      REAL DMIN_THRESHOLD(ARRAY_SIZE)
      LOGICAL REPLY, QUIT, AVAL_WEIGHT
C
C ------------------------------------------------------------------------------
C
      POWER = 2.0
C
      REPLY = .FALSE.
      CALL TASKQ('Do you want experimental fitting options',REPLY,QUIT)
      IF (REPLY) THEN
         FUNCT_TYPE = 2
         AVAL_WEIGHT = .FALSE.
         CALL TASKQ('Do you want 1/aval weighting',AVAL_WEIGHT,QUIT)
         CALL TASKF('Enter power for each term',POWER,1,QUIT)
      ELSE
         FUNCT_TYPE = 1
      END IF
C
      WRITE (STR,6010) FILM
      CALL TWRITE(STR)
C
C====== Select values from array
C
      NUMR1 = 0
C
      DO 20 I = 1 , NSPOT
C
C====== If in range 1 use this spot
C
         IF (WAV(I).LT.WAVMID1) THEN
            NUMR1 = NUMR1 + 1
            IF (NUMR1.GT.NUMSPOTMAX) STOP 'Too many spots selected'
            IF (FUNCT_TYPE.EQ.1) THEN
               ARRAY(1,NUMR1) = WAV(I)
               ARRAY(2,NUMR1) = RAT(I)
               ARRAY(3,NUMR1) = 1.0/SIG2(I)
            ELSE
               ARRAY(1,NUMR1) = WAV(I)
               ARRAY(2,NUMR1) = AVAL(I)
               ARRAY(3,NUMR1) = BVAL(I)
               ARRAY(4,NUMR1) = C2TH(I)
               ARRAY(5,NUMR1) = 1.0/SIG2(I)
               IF (AVAL_WEIGHT) ARRAY(5,NUMR1) = 1.0/AVAL(I)
            END IF
 
         END IF
C
 20   CONTINUE
C
      IF (NUMR1.GT.0) CALL VICFIT(1,NUMR1,FILM)
C
C====== Select values from array again
C
      NUMR2 = 0
C
      DO 30 I = 1 , NSPOT
C
C====== This time want range 2 spots
C
         IF (WAV(I).GE.WAVMID1 .AND. WAV(I).LT.WAVMID2) THEN
            NUMR2 = NUMR2 + 1
            IF (NUMR2.GT.NUMSPOTMAX) STOP 'Too many spots selected'
            IF (FUNCT_TYPE.EQ.1) THEN
               ARRAY(1,NUMR2) = WAV(I)
               ARRAY(2,NUMR2) = RAT(I)
               ARRAY(3,NUMR2) = 1.0/SIG2(I)
            ELSE
               ARRAY(1,NUMR2) = WAV(I)
               ARRAY(2,NUMR2) = AVAL(I)
               ARRAY(3,NUMR2) = BVAL(I)
               ARRAY(4,NUMR2) = C2TH(I)
               ARRAY(5,NUMR2) = 1.0/SIG2(I)
               IF (AVAL_WEIGHT) ARRAY(5,NUMR2) = 1.0/AVAL(I)
            END IF
 
         END IF
 
 30   CONTINUE
C
      IF (NUMR2.GT.0) CALL VICFIT(2,NUMR2,FILM)
      NUMR3 = 0
C
      DO 40 I = 1 , NSPOT
C
C====== This time want range 3 spots
C
         IF (WAV(I).GE.WAVMID2) THEN
            NUMR3 = NUMR3 + 1
            IF (NUMR3.GT.NUMSPOTMAX) STOP 'Too many spots selected'
            IF (FUNCT_TYPE.EQ.1) THEN
               ARRAY(1,NUMR3) = WAV(I)
               ARRAY(2,NUMR3) = RAT(I)
               ARRAY(3,NUMR3) = 1.0/SIG2(I)
            ELSE
               ARRAY(1,NUMR3) = WAV(I)
               ARRAY(2,NUMR3) = AVAL(I)
               ARRAY(3,NUMR3) = BVAL(I)
               ARRAY(4,NUMR3) = C2TH(I)
               ARRAY(5,NUMR3) = 1.0/SIG2(I)
               IF (AVAL_WEIGHT) ARRAY(5,NUMR3) = 1.0/AVAL(I)
            END IF
 
         END IF
C
 40   CONTINUE
C
      IF (NUMR3.GT.0) CALL VICFIT(3,NUMR3,FILM)
C
      CALL QSBELL
      CALL TREAD(STR,QUIT)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Film number',I3)
      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
      CHARACTER*80 STR
      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
      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_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
      CHARACTER*80 STR
      INTEGER    I, RANGE
      DOUBLE PRECISION     DWAV3, K, A, ARG, FF
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
**==get_users_parameters.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************************
C**       GET_USERS_PARAMETERS        **
C***************************************
C
C
      INTEGER   FUNCTION  GET_USERS_PARAMETERS(BOXPROF, SPATREJ,
     +    WAVMIN, WAVMAX, IAMIN, IBMIN, ACON, BCON, ACELL, BCELL,
     +    CCELL, CELLANG, LBADA, LBADB, LENT, LBAND, BAND, FILM)
C
C Purpose: Set data parameters
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    BOXPROF
      REAL       WAVMIN, WAVMAX, ACELL, BCELL, CCELL, CELLANG(3)
      REAL       IAMIN, IBMIN, ACON, BCON
      LOGICAL    LBADA, LBADB, LENT, LBAND, SPATREJ
      REAL       BAND
      INTEGER    FILM
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION  PARAMS(2,3,5)
      COMMON /PARAMK/ K
      DOUBLE PRECISION K
C
      INTEGER    NFILM
      PARAMETER (NFILM=6)
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
      LOGICAL    QUIT
      REAL       R4K, CP(6)
      INTEGER    IIFILM
C
C-----------------------------------------------------------------------------
C
      GET_USERS_PARAMETERS = 0
C
      BOXPROF = 2
 20   CALL TASKI('Use Box (1) or Profile (2) intensities',BOXPROF,QUIT)
      IF (BOXPROF.LE.0 .OR. BOXPROF.GT.2) GO TO 20
      CALL TASKQ('Reject spatial overlaps',SPATREJ,QUIT)
      IIFILM = FILM
      IF (IIFILM.EQ.0) IIFILM = 1
      CALL TASKI('Enter first film of pair to be scaled (A=1,B=2,etc.)',
     +           IIFILM,QUIT)
C
      IF (IIFILM.GE.1 .AND. IIFILM.LE.NFILM-1) FILM = IIFILM
      WRITE (STR,'(I1)') FILM
      CALL TWRITE(STR)
C
      CALL READ_GE_PARAMS(CP,WAVMIN,WAVMAX)
      ACELL = CP(1)
      BCELL = CP(2)
      CCELL = CP(3)
      CELLANG(1) = CP(4)
      CELLANG(2) = CP(5)
      CELLANG(3) = CP(6)
      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
      IF (ACELL.EQ.0.0 .OR. BCELL.EQ.0.0 .OR. CCELL.EQ.0.0) THEN
         GET_USERS_PARAMETERS = 1
         RETURN
      END IF
C
      CALL TASKF('Enter Ix intensity minimum',IAMIN,1,QUIT)
      CALL TASKF('Enter Iy intensity minimum',IBMIN,1,QUIT)
C
      R4K = K
      CALL TASKF('Enter constant K for Victoreen',R4K,3,QUIT)
      K = R4K
C
      CALL TASKQ('Reject non-measurable X spots',LBADA,QUIT)
      CALL TASKQ('Reject non-measurable Y spots',LBADB,QUIT)
C
      CALL TASKF('Enter Ix sigma limit',ACON,1,QUIT)
      CALL TASKF('Enter Iy sigma limit',BCON,1,QUIT)
C
      CALL TASKQ('Use unit weights',LENT,QUIT)
      CALL TASKQ('Exclude spots outside band',LBAND,QUIT)
C
      IF (LBAND) THEN
         CALL TASKF('Enter exclusion band',BAND,2,QUIT)
      END IF
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
**==heap_i.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       HEAP_I        **
C*************************
C
C
C Purpose: This subroutine is used to sort a 2d array of integer data on a
C          given key column
C
      SUBROUTINE HEAP_I(IA,IDIM1,IDIM2,NKEY,NREF)
C
C Parameters:
C
C         IA (I/O)  Integer array containing the data to be sorted with
C                   dimensions (IDIM1,IDIM2)
C      IDIM1 (I)    No. of rows in the array
C      IDIM2 (I)    No. of columns in the array (max=100)
C       NKEY (I)    Column no. for the sort key
C       NREF (I)    The no. of rows to be sorted
C
C====== Specification statements
C
      IMPLICIT INTEGER(A-Z)
      INTEGER IA(IDIM1,IDIM2),ITEMP,CUT(100)
      INTEGER SAVE(2,16)
      DATA CUTOFF/20/
C
C====== Sort data
C
      STACK = 0
      LOW = 1
      HIGH = NREF
 20   CHUNK = HIGH - LOW
      IF (CHUNK.GT.CUTOFF) THEN
         MIDDLE = (LOW+HIGH)/2
         IF (IA(LOW,NKEY).GT.IA(MIDDLE,NKEY)) THEN
            DO 30 I = 1 , IDIM2
               ITEMP = IA(LOW,I)
               IA(LOW,I) = IA(MIDDLE,I)
               IA(MIDDLE,I) = ITEMP
 30         CONTINUE
         END IF
         IF (IA(MIDDLE,NKEY).GT.IA(HIGH,NKEY)) THEN
            DO 40 I = 1 , IDIM2
               ITEMP = IA(MIDDLE,I)
               IA(MIDDLE,I) = IA(HIGH,I)
               IA(HIGH,I) = ITEMP
 40         CONTINUE
            IF (IA(LOW,NKEY).GT.IA(MIDDLE,NKEY)) THEN
               DO 50 I = 1 , IDIM2
                  ITEMP = IA(LOW,I)
                  IA(LOW,I) = IA(MIDDLE,I)
                  IA(MIDDLE,I) = ITEMP
 50            CONTINUE
            END IF
         END IF
         DO 60 I = 1 , IDIM2
            CUT(I) = IA(MIDDLE,I)
 60      CONTINUE
         UP = LOW
         DOWN = HIGH
 70      DOWN = DOWN - 1
         IF (CUT(NKEY).LT.IA(DOWN,NKEY)) GO TO 70
 80      UP = UP + 1
         IF (CUT(NKEY).GT.IA(UP,NKEY)) GO TO 80
         IF (UP.LT.DOWN) THEN
            DO 90 I = 1 , IDIM2
               ITEMP = IA(UP,I)
               IA(UP,I) = IA(DOWN,I)
               IA(DOWN,I) = ITEMP
 90         CONTINUE
            GO TO 70
         ELSE IF (UP.EQ.DOWN) THEN
            DOWN = DOWN - 1
            UP = UP + 1
         END IF
      ELSE
         IF (CHUNK.GE.2) THEN
            MM = HIGH - 1
            DO 120 M = LOW , MM
               MIN = M
               MMM = M + 1
               DO 100 N = MMM , HIGH
                  IF (IA(MIN,NKEY).GT.IA(N,NKEY)) THEN
                     MIN = N
                  END IF
 100           CONTINUE
               IF (MIN.NE.M) THEN
                  DO 110 I = 1 , IDIM2
                     ITEMP = IA(M,I)
                     IA(M,I) = IA(MIN,I)
                     IA(MIN,I) = ITEMP
 110              CONTINUE
               END IF
 120        CONTINUE
         END IF
         IF (STACK.EQ.0) RETURN
         LOW = SAVE(1,STACK)
         HIGH = SAVE(2,STACK)
         STACK = STACK - 1
         GO TO 20
      END IF
      STACK = STACK + 1
      IF (STACK.GT.16) THEN
C
C====== Error condition
C
         CALL TWRITE('** Stack limit of 16 exceeded in heap_i **')
         STOP
      ELSE
         BOTTOM = DOWN - LOW
         TOP = HIGH - UP
         IF (BOTTOM.GT.TOP) THEN
            SAVE(1,STACK) = LOW
            SAVE(2,STACK) = DOWN
            LOW = UP
         ELSE
            SAVE(1,STACK) = UP
            SAVE(2,STACK) = HIGH
            HIGH = DOWN
         END IF
         GO TO 20
      END IF
      END
**==lorentz.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       LORENTZ       **
C*************************
C
      REAL   FUNCTION  LORENTZ (THETA)
C
C Purpose: Apply Lorentz correction
C
      INCLUDE 'implicit_none.finc'
C
      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 file: ')
         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)
C
      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)
      IFAIL = 0
      CALL CCPDPN(-IN_DATA2,FNAME2,'READONLY','DU',2,IFAIL)
C
      WRITE (STR,6010) FNAME1(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 ('Input .ge1/.ge2 file root name is ',A)
      END
**==output.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       OUTPUT        **
C*************************
C
C
      SUBROUTINE  OUTPUT (ACELL, BCELL, CCELL, CELLANG)
C
C Purpose: To output reflexions suitable as input for LAUESCAL/LAUENORM
C
      INCLUDE 'implicit_none.finc'
C
      REAL       ACELL, BCELL, CCELL, CELLANG(3)
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION     PARAMS(2,3,5)
C
      COMMON /BP/ BOXPROF
      INTEGER   BOXPROF
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
      INTEGER    SCALE_SPOT
      INTEGER    IN_DATA2
      PARAMETER (IN_DATA2=3)
      CHARACTER*60  FNAM, OUT_FILE
      LOGICAL    QUIT
      REAL       SCALE_FACTOR
      INTEGER    I, J, END_REC, RECORD_NUM, NS
      INTEGER    IPACK, NFILM, ITSPOT
      INTEGER    GET_GE_DATA
      REAL       X, Y, INTENSITIES(6), SIGMAS(6),
     +           SCALED_INTENS, SCALED_SIGMA
      REAL       INTENB(6), SIGB(6), INTENP(6), SIGP(6)
      INTEGER    MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP
      INTEGER    H_S(6), K_S(6), L_S(6)
      REAL       LAMBDAS(6)
      REAL       DMIN_THRESHOLD
      INTEGER    H, K, L
      INTEGER    I_INTENS, I_SIGMA
      REAL       LAMBDA, THETA, LOR, P
      INTEGER    MAX_BAD_SPOTS, BAD_SPOT_COUNT, IFAIL
      PARAMETER (MAX_BAD_SPOTS=10)
      LOGICAL    APPLY_SIGMA, APPLY_LAMBDA, APPLY_MININT, SPATREJ
      REAL	 SIGMA_LIM
      INTEGER    MININT
      REAL       LAMBDA_LLIM, LAMBDA_HLIM
      REAL       DEV, TOT_DEV, ONE_DEV, TOT_DENOM
      INTEGER    NUM_DEV
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
C-----------------------------------------------------------------------------
C
      BAD_SPOT_COUNT = 0
      TOT_DEV = 0.0
      TOT_DENOM = 0.0
C
C====== Open output file for LAUESCAL/LAUENORM
C
 20   CALL TPROMP('Enter name for output file (default ext=.out): ')
      CALL TREAD(FNAM,QUIT)
      CALL FILEXT(FNAM,'out',OUT_FILE,I)
      CALL CCPDPN(-4,OUT_FILE,'NEW','F',0,IFAIL)
      IF (IFAIL.EQ.-1) GO TO 20
      WRITE (STR,6010) OUT_FILE(1:I)
      CALL TWRITE(STR)
C
C====== Get overall scaling factor (changed to 100. when L reciprocal)
C
      SCALE_FACTOR = 100.
      CALL TASKF('Input scaling factor',SCALE_FACTOR,2,QUIT)
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
C====== Get rejection criteria
C
      CALL TASKQ('Do you want to reject spatials',SPATREJ,QUIT)
      CALL TASKQ('Do you want to apply a sigma rejection criterion',
     +           APPLY_SIGMA,QUIT)
      IF (APPLY_SIGMA) THEN
         CALL TASKF('Enter sigma limit (reject I<x*sigma)',SIGMA_LIM,
     +              2,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
C====== Now read in singlets from .ge file
C
      READ (IN_DATA,REC=10) IPACK , ITSPOT
      END_REC = 10 + ITSPOT
      NS = 0
C
      DO 60 RECORD_NUM = 11 , END_REC
C
         MULTIPLICITY = 6
         IF (GET_GE_DATA(IN_DATA,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) THEN
            CALL TWRITE(' ')
            WRITE (STR,6050)
            CALL TWRITE(STR)
            CALL TWRITE(' ')
            CLOSE (4)
            RETURN
         END IF
         IF (.NOT.(SPATREJ .AND. SPATIAL_OVERLAP.EQ.1)) THEN
            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
C====== If this is a singlet then output
C
            IF (MULTIPLICITY.EQ.1) THEN
C
               H = H_S(1)
               K = K_S(1)
               L = L_S(1)
C
               LAMBDA = LAMBDAS(1)
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,6060) MAX_BAD_SPOTS
                     CALL TWRITE(STR)
                     CALL TWRITE(' ')
                     WRITE (STR,6070)
                     CALL TWRITE(STR)
                     CALL TWRITE(' ')
                     GO TO 70
                  END IF
C
C====== Check the criteria
C
               ELSE IF (.NOT.(APPLY_SIGMA .AND. I_INTENS.LT.
     +      (nint(SIGMA_LIM*float(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 - ignore zero sigmas
c
                        if(i_sigma.ne.0)then
C
C====== Collect stats for merging R-factor
C====== NB SCALE_SPOT has scaled up the intensities.
C
                        ONE_DEV = 0.0
                        NUM_DEV = 0
                        DO 50 I = 1 , 6
                           IF (SIGMAS(I).LT.1E6) THEN
                              DEV = ABS(INTENSITIES(I)-SCALED_INTENS)
                              IF (DEV.LT.SCALED_SIGMA*2.) THEN
                                 ONE_DEV = ONE_DEV + DEV
                                 NUM_DEV = NUM_DEV + 1
                              END IF
                           END IF
 50                     CONTINUE
                        IF (NUM_DEV.GT.0) THEN
                           TOT_DEV = TOT_DEV + ONE_DEV/NUM_DEV
                           TOT_DENOM = TOT_DENOM + SCALED_INTENS
                        END IF
C
C====== Write out to output file
C
                        NS = NS + 1
                        WRITE (4,6020) H , K , L , LAMBDA , THETA , 
     +                                 I_INTENS , I_SIGMA , X , Y
                        end if
                     END IF
                  END IF
               END IF
C
            END IF
         END IF
C
 60   CONTINUE
C
C====== Calculate merging R
C
      IF (TOT_DENOM.LE.0.0) STOP ' Bad denominator in R-merge calc.'
 70   CALL TWRITE(' ')
      CALL TWRITE(' ')
      WRITE (STR,6030) NS
      CALL TWRITE(STR)
      WRITE (STR,6040) TOT_DEV/TOT_DENOM*100.
      CALL TWRITE(STR)
C
      CLOSE (4)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Writing output to ',A)
 6020 FORMAT (3I5,2F10.5,2I10,2F8.3)
 6030 FORMAT ('Total number of spots output',I6)
 6040 FORMAT ('Merging R-factor  = ',F5.1)
 6050 FORMAT ('!!! Error in reading .ge1/.ge2 files. Aborting output.')
 6060 FORMAT ('!!! More than',I4,' bad spots detected.')
 6070 FORMAT ('    *****  ABORTING  *****')
      END
**==plot.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       PLOT          **
C*************************
C
C
      SUBROUTINE  PLOT ( FITTED, NSPOT, FILM, WAVMIN, WAVMAX,
     +     ARRAY_SIZE, WAV, RAT )
C
C Purpose: To plot film factor vs lambda points, and also a fitted curve
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FITTED, NSPOT, FILM
      DOUBLE PRECISION     WAVMIN, WAVMAX
C
C FITTED : flag indicating that data has been fitted (1=fit made)
C NSPOT : number of accepted singlets
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
      LOGICAL    QUIT
      CHARACTER*1 ANSWER
      DOUBLE PRECISION     LAMBDA, P
      REAL       D_LAMBDA, PP
      DOUBLE PRECISION     EVAL_VIC
      INTEGER    RANGE, I, IPSTYP
      INTEGER    ARRAY_SIZE
      REAL       WAV(ARRAY_SIZE), RAT(ARRAY_SIZE)
      LOGICAL    SAVPLT
      CHARACTER  PSFILE*80, TEXT*1
C
C------------------------------------------------------------------------------
C
      CALL QGPSET(2)
      CALL QGENTR
      CALL QGCLR
      CALL QGMAPC(100.0,700.0,100.0,700.0,0.0,3.0,0.0,7.0)
      CALL QGXAX1(0.0,0.0,3.0,0.0,3.0,15,1,3.0,12.0)
      CALL QGYAX1(0.0,0.0,7.0,0.0,7.0,14,1,3.0,12.0)
C
      IF (FITTED.NE.0) THEN
C
         DO 20 LAMBDA = WAVMIN , WAVMAX , 0.05
            P = EVAL_VIC(LAMBDA,RANGE,FILM)
            PP = P
            D_LAMBDA = LAMBDA
            IF (RANGE.EQ.1 .OR. RANGE.EQ.2 .OR. RANGE.EQ.3) THEN
               CALL QGSYMB(D_LAMBDA,PP,1,2.0)
            END IF
 20      CONTINUE
C
         CALL QGBELL
         CALL TREAD(ANSWER,QUIT)
C
      END IF
C
      IF (NSPOT.EQ.0) RETURN
C
      DO 30 I = 1 , NSPOT
         CALL QGSYMB(WAV(I),RAT(I),3,1.0)
 30   CONTINUE
C
      CALL QGBELL
      CALL TREAD(ANSWER,QUIT)
      CALL QGALPH
      CALL TASKQ('Do you want to save plot in a Postscript file',SAVPLT,
     +           QUIT)
      IF (SAVPLT) THEN
         CALL TPROMP('Name of Postscript file: ')
         CALL TREAD(PSFILE,QUIT)
         IF (PSFILE.NE.' ') THEN
            OPEN (UNIT=11,FILE=PSFILE)
            CALL TPROMP('Portrait (P) or landscape (L) format [P]: ')
            CALL TREAD(TEXT,QUIT)
            IPSTYP = 2
            IF (TEXT.EQ.'l' .OR. TEXT.EQ.'L') IPSTYP = 1
            CALL QGWRPS(11,IPSTYP)
         END IF
      END IF
C
      RETURN
      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: Read lambda limits and cell parameters from .GE1 file header. Expand
C          cell parameters up to complete set (depends on xtal system) and
C          convert 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
C-------------------------------------------------------------------------------
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)
C
      COMMON /IO_STREAMS/IN_TERM,OUT_TERM,IN_DATA,OUT_DATA
      INTEGER IN_TERM,OUT_TERM,IN_DATA,OUT_DATA
C
      LOGICAL CCPEXS
      INTEGER LENSTR
C
C====== LOCALS:
C
      CHARACTER*80 STR
      DOUBLE PRECISION AL150(2,3,5),AL200(2,3,5)
      CHARACTER *60 FNAM,FILENAME
      LOGICAL EXISTS, QUIT
      INTEGER COMM
      INTEGER I,J,K,TYPE,IFAIL
      DOUBLE PRECISION FILM(3),METAL(3),THICK
      REAL THICKR
      DATA FILM/1.975,.682,.2856/
C
C FILM: Holds absorbancy coefs. for film at 3 wavelength ranges
C
      DATA METAL/.0034415,.1198245,0/
C
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,
     1     2.4912,1,1.1982,1,.8018,1,2.4912,1,1.1982,1,.8018,1,7.1454,1,
     1     5.276,1,4.1812/
      DATA AL200/1,1.975,1,.682,1,.2856,1,2.6633,1,1.3703,1,.9739,1,
     1     2.6633,1,1.3703,1,.9739,1,2.6633,1,1.3703,1,.9739,1,7.1454,1,
     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(' ')
      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(' ')
      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(' ')
            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
               THICKR = THICK
               CALL TASKF('Enter foil thickness in microns',THICKR,2,
     +                    QUIT)
               THICK = THICKR
 
            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
C
         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
**==rsymm.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       RSYMM         **
C*************************
C
C
      SUBROUTINE  RSYMM (NUMSPOTS, ARRAY_SIZE, WAV, FACTOR, WEIGHT,
     +            AVAL, BVAL, IH, IK, IL, C2TH, DHKL, DMIN_THRESHOLD)
C
C Arrays WAV and FACTOR etc replace unit 9 and direct reference is made to
C them, where the previous wav etc are writen to, the variables have been
C replaced by IIH, IIK ETC
C
C Purpose: Calculate R-factor of intensities of spots hkl,-hkl for laue films
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    NUMSPOTS, ARRAY_SIZE
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
      INTEGER    INDIFA(11),NINDIFA(11),LDIFA(12),NLDIFA(12)
      INTEGER    INDIFB(11),NINDIFB(11),LDIFB(12),NLDIFB(12)
      INTEGER    INBOTA(11),LBOTA(12)
      INTEGER    INBOTB(11),LBOTB(12)
      INTEGER    IDDIFA(9),NIDDIFA(9),IDBOTA(9)
      INTEGER    IDDIFB(9),NIDDIFB(9),IDBOTB(9)
      REAL       R(12)
      INTEGER    TABLE(8000,4)
      INTEGER    ARR(8000,8)
      INTEGER    I, J, I9
      INTEGER    COMPARE_INDEX, ISUMA, NSUMA, IBOTA
      INTEGER    ISUMB, NSUMB, IBOTB
      INTEGER    IAVAL, IBVAL, IH(ARRAY_SIZE), IK(ARRAY_SIZE)
      INTEGER    ISPOT, INDBAR, NUMR1, IIH, IIK, IIL, IL(ARRAY_SIZE)
      INTEGER    IAVAL1, IAV, IDIF, IBVAL1
      REAL       AVAL(ARRAY_SIZE), BVAL(ARRAY_SIZE), FACTOR(ARRAY_SIZE)
      REAL       WAVDLIN, WAV1, WAV(ARRAY_SIZE), WEIGHT(ARRAY_SIZE)
      REAL       C2TH(ARRAY_SIZE), DHKL(ARRAY_SIZE), DDHKL, WWAV
      REAL       RAT, LAMBDA_MIS_LIMIT, RR, DMIN_THRESHOLD(ARRAY_SIZE)
      INTEGER    NUM_MISMATCHED,EMPTY_SLOT(20)
C
C EMPTY_SLOT array:
C (1) holds the number of `holes' i.e. entries in TABLE that have been deleted.
C (2) holds position of last slot in TABLE allocated
C (3--19) point to the holes - given value 9999 when hole filled.
C (20) Highest value (1) has taken so far.
C
      CHARACTER*1  REPLY
      LOGICAL    QUIT
      INTEGER   O
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
C====== Reset table
C
      DO 20 I = 1 , 8000
         TABLE(I,4) = 0
 20   CONTINUE
C
C====== Clear arrays
C
      DO 30 I = 1 , 11
         INDIFA(I) = 0
         NINDIFA(I) = 0
         INBOTA(I) = 0
 30   CONTINUE
C
      DO 40 I = 1 , 12
         LDIFA(I) = 0
         NLDIFA(I) = 0
         LBOTA(I) = 0
 40   CONTINUE
C
      ISUMA = 0
      NSUMA = 0
      IBOTA = 0
C
      DO 50 I = 1 , 11
         INDIFB(I) = 0
         NINDIFB(I) = 0
         INBOTB(I) = 0
 50   CONTINUE
C
      DO 60 I = 1 , 12
         LDIFB(I) = 0
         NLDIFB(I) = 0
         LBOTB(I) = 0
 60   CONTINUE
C
      DO 70 I = 1 , 9
         IDDIFA(I) = 0
         NIDDIFA(I) = 0
         IDBOTA(I) = 0
         IDDIFB(I) = 0
         NIDDIFB(I) = 0
         IDBOTB(I) = 0
 70   CONTINUE
C
      ISUMB = 0
      NSUMB = 0
      IBOTB = 0
C
C====== Read values into core
C
      NUMR1 = 1
C
      IF (COMPARE_INDEX.EQ.1) THEN
         DO 80 I = 1 , NUMSPOTS
            IF (IH(I).NE.0) THEN
               IAVAL = AVAL(I)
               IBVAL = BVAL(I)
               ARR(NUMR1,1) = WAV(I)*100000.0 + .5
               ARR(NUMR1,3) = IAVAL
               ARR(NUMR1,4) = IBVAL
               ARR(NUMR1,5) = IABS(IH(I))
               ARR(NUMR1,6) = IK(I)
               ARR(NUMR1,7) = IL(I)
               ARR(NUMR1,8) = DHKL(I)*100. + .5
               NUMR1 = NUMR1 + 1
               IF (NUMR1.GT.8000) STOP ' Storage overflowed in RSYMM'
            END IF
 80      CONTINUE
      ELSE IF (COMPARE_INDEX.EQ.2) THEN
         DO 90 I = 1 , NUMSPOTS
            IF (IK(I).NE.0) THEN
               IAVAL = AVAL(I)
               IBVAL = BVAL(I)
               ARR(NUMR1,1) = WAV(I)*100000.0 + .5
               ARR(NUMR1,3) = IAVAL
               ARR(NUMR1,4) = IBVAL
               ARR(NUMR1,5) = IH(I)
               ARR(NUMR1,6) = IABS(IK(I))
               ARR(NUMR1,7) = IL(I)
               ARR(NUMR1,8) = DHKL(I)*100. + .5
               NUMR1 = NUMR1 + 1
               IF (NUMR1.GT.8000) STOP ' Storage overflowed in RSYMM'
            END IF
 90      CONTINUE
      ELSE
         DO 100 I = 1 , NUMSPOTS
            IF (IL(I).NE.0) THEN
               IAVAL = AVAL(I)
               IBVAL = BVAL(I)
               ARR(NUMR1,1) = WAV(I)*100000.0 + .5
               ARR(NUMR1,3) = IAVAL
               ARR(NUMR1,4) = IBVAL
               ARR(NUMR1,5) = IH(I)
               ARR(NUMR1,6) = IK(I)
               ARR(NUMR1,7) = IABS(IL(I))
               ARR(NUMR1,8) = DHKL(I)*100. + .5
               NUMR1 = NUMR1 + 1
               IF (NUMR1.GT.8000) STOP ' Storage overflowed in RSYMM'
            END IF
 100     CONTINUE
      END IF
C
      NUMR1 = NUMR1 - 1
C
C====== Sort array
C
      CALL HEAP_I(ARR,8000,8,1,NUMR1)
C
C====== Array is now sorted
C
      NUM_MISMATCHED = 0
C
C====== Point at first elt of TABLE
C
      EMPTY_SLOT(2) = 1
C
C======  START OF ISPOT LOOP
C======  ===================
C
      DO 150 ISPOT = 1 , NUMR1
         IIH = ARR(ISPOT,5)
         IIK = ARR(ISPOT,6)
         IIL = ARR(ISPOT,7)
C
         DO 110 J = 1 , EMPTY_SLOT(2) - 1
            IF (TABLE(J,4).GT.0) THEN
               IF (TABLE(J,1).EQ.IIH) THEN
                  IF (TABLE(J,2).EQ.IIK) THEN
                     IF (TABLE(J,3).EQ.IIL) GO TO 130
                  END IF
               END IF
            END IF
 110     CONTINUE
C
C====== Indices not found - put in table
C====== first find empty slot
C
         IF (EMPTY_SLOT(1).EQ.0) THEN
C
C====== No internal holes - have to place at end.
C
            J = EMPTY_SLOT(2)
            EMPTY_SLOT(2) = EMPTY_SLOT(2) + 1
         ELSE
C
C====== Search through to find lowest empty slot
C
            J = EMPTY_SLOT(3)
            O = 3
            DO 120 I = 4 , 19
               IF (EMPTY_SLOT(I).LT.J .AND. (EMPTY_SLOT(I).GT.0)) THEN
                  J = EMPTY_SLOT(I)
                  O = I
               END IF
 120        CONTINUE
C
C====== Record that this slot in TABLE is now filled
C
            EMPTY_SLOT(O) = 9999
C
C====== Decrease count of internal holes
C
            EMPTY_SLOT(1) = EMPTY_SLOT(1) - 1
         END IF
         IF (EMPTY_SLOT(2).LE.8000) THEN
C
C====== Now put in values
C
            TABLE(J,1) = IIH
            TABLE(J,2) = IIK
            TABLE(J,3) = IIL
            TABLE(J,4) = ISPOT
            GO TO 150
         ELSE
C
            CALL TWRITE(
     +      'Table full at 8000 reflexions - partially complete R-symms'
     +      )
            GO TO 160
         END IF
C
C====== Found indices
C
 130     INDBAR = TABLE(J,4)
C
C====== Going to delete this TABLE elt creating an internal hole - inc the count
C====== of holes.
C
         EMPTY_SLOT(1) = EMPTY_SLOT(1) + 1
C
C====== Want to record, in EMPTY_SLOT, position of the deleted entry in TABLE -
C====== but first have to find an empty slot in EMPTY_SLOT!
C
         I = 2
 140     I = I + 1
C
C====== If we've freed too many holes at once (ie. if I>19) just have to forget
C====== that this slot is empty as we haven't room to store a pointer to it.
C
         IF (I.LE.19) THEN
            IF ((EMPTY_SLOT(I).GT.0) .AND. (EMPTY_SLOT(I).LT.9999))
     +          GO TO 140
            EMPTY_SLOT(I) = J
            IF (EMPTY_SLOT(1).GT.EMPTY_SLOT(20)) EMPTY_SLOT(20)
     +          = EMPTY_SLOT(1)
         END IF
C
C====== Delete this entry
C
         TABLE(J,4) = 0
         WWAV = ARR(ISPOT,1)/100000.0
         WAV1 = ARR(INDBAR,1)/100000.0
C
         IF (ABS(WWAV-WAV1).GT.LAMBDA_MIS_LIMIT) THEN
            NUM_MISMATCHED = NUM_MISMATCHED + 1
            GO TO 150
         END IF
C
         DDHKL = ARR(ISPOT,8)/100.
C
C====== Do the X film
C
         IAVAL = ARR(ISPOT,3)
         IAVAL1 = ARR(INDBAR,3)
C
         IF (IAVAL.NE.9999 .AND. IAVAL.NE.-9999) THEN
            IF (IAVAL1.NE.9999 .AND. IAVAL1.NE.-9999) THEN
C
               IAV = (IAVAL+IAVAL1)/2. + .5
               IDIF = IABS(IAVAL-IAV) + IABS(IAVAL1-IAV)
               I = IAV/200 + 1
C
C====== Do intensity ranges
C
               IF (I.GT.11) I = 11
C
               INDIFA(I) = INDIFA(I) + IDIF
               NINDIFA(I) = NINDIFA(I) + 1
               INBOTA(I) = INBOTA(I) + IAV
C
C====== Do wavelength ranges
C
               J = (WWAV-0.1)/.2 + 1
C
               IF (J.LT.1) J = 1
               IF (J.GT.12) J = 12
C
               LDIFA(J) = LDIFA(J) + IDIF
               LBOTA(J) = LBOTA(J) + IAV
               NLDIFA(J) = NLDIFA(J) + 1
C
C====== Collect Ddhkl diffs
C====== which range ?
C
               IF (DDHKL.GT.20.) THEN
                  J = 1
               ELSE IF (DDHKL.LE.20. .AND. DDHKL.GT.10.) THEN
                  J = 2
               ELSE IF (DDHKL.LE.10. .AND. DDHKL.GT.7.) THEN
                  J = 3
               ELSE IF (DDHKL.LE.7. .AND. DDHKL.GT.5.) THEN
                  J = 4
               ELSE IF (DDHKL.LE.5. .AND. DDHKL.GT.4.) THEN
                  J = 5
               ELSE IF (DDHKL.LE.4. .AND. DDHKL.GT.3.33) THEN
                  J = 6
               ELSE IF (DDHKL.LE.3.33 .AND. DDHKL.GT.2.9) THEN
                  J = 7
               ELSE IF (DDHKL.LE.2.9 .AND. DDHKL.GT.2.6) THEN
                  J = 8
               ELSE
                  J = 9
               END IF
C
               IDDIFA(J) = IDDIFA(J) + IDIF
               NIDDIFA(J) = NIDDIFA(J) + 1
               IDBOTA(J) = IDBOTA(J) + IAV
C
               ISUMA = ISUMA + IDIF
               NSUMA = NSUMA + 1
               IBOTA = IBOTA + IAV
            END IF
         END IF
C
C====== Do the Y film
C
         IBVAL = ARR(ISPOT,4)
         IBVAL1 = ARR(INDBAR,4)
C
         IF (IBVAL.NE.9999 .AND. IBVAL.NE.-9999) THEN
            IF (IBVAL1.NE.9999 .AND. IBVAL1.NE.-9999) THEN
C
               IAV = (IBVAL+IBVAL1)/2. + .5
               IDIF = IABS(IBVAL-IAV) + IABS(IBVAL1-IAV)
               I = IAV/200 + 1
C
               IF (I.GT.11) I = 11
C
               INDIFB(I) = INDIFB(I) + IDIF
               NINDIFB(I) = NINDIFB(I) + 1
               INBOTB(I) = INBOTB(I) + IAV
C
               J = (WWAV-0.1)/.2 + 1
C
               IF (J.LT.1) J = 1
               IF (J.GT.12) J = 12
C
               LDIFB(J) = LDIFB(J) + IDIF
               LBOTB(J) = LBOTB(J) + IAV
               NLDIFB(J) = NLDIFB(J) + 1
C
C====== Collect Ddhkl diffs
C====== which range ?
C
               IF (DDHKL.GT.20.) THEN
                  J = 1
               ELSE IF (DDHKL.LE.20. .AND. DDHKL.GT.10.) THEN
                  J = 2
               ELSE IF (DDHKL.LE.10. .AND. DDHKL.GT.7.) THEN
                  J = 3
               ELSE IF (DDHKL.LE.7. .AND. DDHKL.GT.5.) THEN
                  J = 4
               ELSE IF (DDHKL.LE.5. .AND. DDHKL.GT.4.) THEN
                  J = 5
               ELSE IF (DDHKL.LE.4. .AND. DDHKL.GT.3.33) THEN
                  J = 6
               ELSE IF (DDHKL.LE.3.33 .AND. DDHKL.GT.2.9) THEN
                  J = 7
               ELSE IF (DDHKL.LE.2.9 .AND. DDHKL.GT.2.6) THEN
                  J = 8
               ELSE
                  J = 9
               END IF
C
               IDDIFB(J) = IDDIFB(J) + IDIF
               NIDDIFB(J) = NIDDIFB(J) + 1
               IDBOTB(J) = IDBOTB(J) + IAV
C
               ISUMB = ISUMB + IDIF
               NSUMB = NSUMB + 1
               IBOTB = IBOTB + IAV
            END IF
         END IF
C
C====== END OF ISPOT LOOP
C====== =================
C
 150  CONTINUE
C
 160  IF (NUM_MISMATCHED.NE.1) THEN
         WRITE (STR,6010) NUM_MISMATCHED
         CALL TWRITE(STR)
      ELSE
         WRITE (STR,6020)
         CALL TWRITE(STR)
      END IF
      DO 170 I = 1 , 11
         IF (NINDIFA(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(INDIFA(I))/INBOTA(I)/2.
         END IF
 170  CONTINUE
C
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      WRITE (STR,6030)
      CALL TWRITE(STR)
      WRITE (STR,6040)
      CALL TWRITE(STR)
      WRITE (STR,6050)
      CALL TWRITE(STR)
C
      WRITE (STR,6060) (NINDIFA(I),I=1,11)
      CALL TWRITE(STR)
      WRITE (STR,6070) (R(I),I=1,11)
      CALL TWRITE(STR)
C
      CALL TWRITE(' ')
      WRITE (STR,6080)
      CALL TWRITE(STR)
      WRITE (STR,6090)
      CALL TWRITE(STR)
C
      DO 180 I = 1 , 12
         IF (NLDIFA(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(LDIFA(I))/LBOTA(I)/2.
         END IF
 180  CONTINUE
C
      WRITE (STR,6100) (NLDIFA(I),I=1,12)
      CALL TWRITE(STR)
      WRITE (STR,6110) (R(I),I=1,12)
      CALL TWRITE(STR)
C
      DO 190 I = 1 , 9
         IF (NIDDIFA(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(IDDIFA(I))/IDBOTA(I)/2.
         END IF
 190  CONTINUE
C
      CALL TWRITE(' ')
      WRITE (STR,6120)
      CALL TWRITE(STR)
      WRITE (STR,6130)
      CALL TWRITE(STR)
C
      WRITE (STR,6140) (NIDDIFA(I),I=1,9)
      CALL TWRITE(STR)
      WRITE (STR,6150) (R(I),I=1,9)
      CALL TWRITE(STR)
C
      IF (IBOTA.EQ.0) THEN
         RR = 0.0
      ELSE
         RR = FLOAT(ISUMA)/IBOTA*100.0/2.0
      END IF
C
      CALL TWRITE(' ')
      WRITE (STR,6160) RR
      CALL TWRITE(STR)
C
      DO 200 I = 1 , 11
         IF (NINDIFB(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(INDIFB(I))/INBOTB(I)/2.
         END IF
 200  CONTINUE
C
      CALL QSBELL
      CALL TREAD(STR,QUIT)
C
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      WRITE (STR,6170)
      CALL TWRITE(STR)
      WRITE (STR,6180)
      CALL TWRITE(STR)
      WRITE (STR,6190)
      CALL TWRITE(STR)
C
      WRITE (STR,6060) (NINDIFB(I),I=1,11)
      CALL TWRITE(STR)
      WRITE (STR,6070) (R(I),I=1,11)
      CALL TWRITE(STR)
C
      CALL TWRITE(' ')
      WRITE (STR,6080)
      CALL TWRITE(STR)
      WRITE (STR,6090)
      CALL TWRITE(STR)
C
      DO 210 I = 1 , 12
         IF (NLDIFB(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.0*FLOAT(LDIFB(I))/LBOTB(I)/2.0
         END IF
 210  CONTINUE
C
      WRITE (STR,6100) (NLDIFB(I),I=1,12)
      CALL TWRITE(STR)
      WRITE (STR,6110) (R(I),I=1,12)
      CALL TWRITE(STR)
C
      DO 220 I = 1 , 9
         IF (NIDDIFB(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.0*FLOAT(IDDIFB(I))/IDBOTB(I)/2.0
         END IF
 220  CONTINUE
C
      CALL TWRITE(' ')
      WRITE (STR,6120)
      CALL TWRITE(STR)
      WRITE (STR,6130)
      CALL TWRITE(STR)
C
      WRITE (STR,6140) (NIDDIFB(I),I=1,9)
      CALL TWRITE(STR)
      WRITE (STR,6150) (R(I),I=1,9)
      CALL TWRITE(STR)
C
      IF (IBOTB.EQ.0) THEN
         RR = 0.0
      ELSE
         RR = FLOAT(ISUMB)/IBOTB*100.0/2.0
      END IF
C
      CALL TWRITE(' ')
      WRITE (STR,6160) RR
      CALL TWRITE(STR)
C
      CALL QSBELL
      CALL TREAD(STR,QUIT)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('There were ',I5,' wavelength mismatches')
 6020 FORMAT ('There was only one wavelength mismatch')
 6030 FORMAT ('X film')
 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 ('R-symms on wavelength')
 6090 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')
 6100 FORMAT (12I6)
 6110 FORMAT (12F6.1)
 6120 FORMAT ('R-symms on dhkl')
 6130 FORMAT ('   >20    10     7     5     4  3.33   2.9   2.6  <2.6')
 6140 FORMAT (9I6)
 6150 FORMAT (9F6.1)
 6160 FORMAT ('Overall R-symm =',F6.1)
 6170 FORMAT ('Y film')
 6180 FORMAT ('R-symms on intensity')
 6190 FORMAT ('   200   400   600   800  1000  1200  1400  1600',
     +        ' 1800  2000 >2000')
      END
**==scale.f  
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       SCALE         **
C*************************
C
C
      SUBROUTINE SCALE(NSPOT,FILM,ARRAY_SIZE,WAV,EXPFACTOR,WEIGHT,
     +                 AVAL,BVAL,IH,IK,IL,C2TH,DHKL,DMIN_THRESHOLD)
C
C Purpose: To calculate R-factors for the film scaling factors
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER NSPOT, FILM, ARRAY_SIZE
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
      LOGICAL QUIT
      CHARACTER*80 STR
      INTEGER INDIF(11), NINDIF(11), LDIF(12), NLDIF(12)
      INTEGER LCOUNTLO(12), LCOUNTHI(12)
      INTEGER INBOT(11), LBOT(12), IDDIF(9), NIDDIF(9), IDBOT(9)
      INTEGER I, J, NS, ISUM, NSUM, IBOT, RANGE
      INTEGER IDIF
      REAL WAV(ARRAY_SIZE), EXPFACTOR(ARRAY_SIZE)
      REAL WEIGHT(ARRAY_SIZE), AVAL(ARRAY_SIZE)
      REAL BVAL(ARRAY_SIZE), SCALED_BVAL
      INTEGER IH(ARRAY_SIZE), IK(ARRAY_SIZE), IL(ARRAY_SIZE)
      REAL C2TH(ARRAY_SIZE), DHKL(ARRAY_SIZE)
      REAL DMIN_THRESHOLD(ARRAY_SIZE)
      REAL R(12), RR
      DOUBLE PRECISION DWAV, P, EVAL_VIC
      REAL BAND, FACTOR
      PARAMETER (BAND=0.20)
      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
C====== Clear arrays
C
      DO 20 I = 1 , 11
         INDIF(I) = 0
         NINDIF(I) = 0
         LDIF(I) = 0
         NLDIF(I) = 0
         LCOUNTLO(I) = 0
         LCOUNTHI(I) = 0
         INBOT(I) = 0
         LBOT(I) = 0
 20   CONTINUE
C
      DO 30 I = 1 , 9
         IDDIF(I) = 0
         NIDDIF(I) = 0
         IDBOT(I) = 0
 30   CONTINUE
C
      ISUM = 0
      NSUM = 0
      IBOT = 0
C
      DO 60 NS = 1 , NSPOT
C
         IF (ABS(AVAL(NS)).NE.9999.0 .AND. ABS(BVAL(NS)).NE.9999.0) THEN
C
            DWAV = WAV(NS)
            P = EVAL_VIC(DWAV,RANGE,FILM)
C
            IF (RANGE.EQ.1 .OR. RANGE.EQ.2 .OR. RANGE.EQ.3) THEN
C
               IF (P.GT.0.0) THEN
                  FACTOR = P**(1.0/C2TH(NS))
               ELSE
                  STOP 
     +               '!!! Attempted to calculate log of negative number'
 
               END IF
C
C====== Collect intensity diffs
C
               SCALED_BVAL = BVAL(NS)*FACTOR + 0.5
               IDIF = ABS(AVAL(NS)-SCALED_BVAL)
C
C====== Which intensity range ?
C
               I = BVAL(NS)/200.0 + 1.0
               IF (I.GT.11) I = 11
C
C====== Add up the diffs
C
               INDIF(I) = INDIF(I) + IDIF
C
C====== Keep a count
C
               NINDIF(I) = NINDIF(I) + 1
C
C====== Add up the denominator
C
               INBOT(I) = INBOT(I) + AVAL(NS)
C
C====== Collect wavelength diffs
C====== which wavlength range ?
C
               J = (WAV(NS)-0.1)/.2 + 1
C
               IF (J.LT.1) J = 1
               IF (J.GT.12) J = 12
C
               LDIF(J) = LDIF(J) + IDIF
               LBOT(J) = LBOT(J) + AVAL(NS)
               NLDIF(J) = NLDIF(J) + 1
C
               IF (ABS((EXPFACTOR(NS)-P)/P).LE.BAND) THEN
                  IF ((EXPFACTOR(NS)-P).LT.0) LCOUNTLO(J) = LCOUNTLO(J)
     +                + 1
                  IF ((EXPFACTOR(NS)-P).GE.0) LCOUNTHI(J) = LCOUNTHI(J)
     +                + 1
               END IF
C
C====== Collect dhkl diffs
C====== which range ?
C
               DO 40 J = 1 , NUM_DHKL_RANGES
                  IF (DHKL(NS).GT.DHKL_RANGES(J)) GO TO 50
 40            CONTINUE
               J = NUM_DHKL_RANGES
C
 50            IDDIF(J) = IDDIF(J) + IDIF
               NIDDIF(J) = NIDDIF(J) + 1
               IDBOT(J) = IDBOT(J) + AVAL(NS)
C
C====== Collect overall diffs
C
               ISUM = ISUM + IDIF
               NSUM = NSUM + 1
               IBOT = IBOT + AVAL(NS)
C
            END IF
         END IF
C
 60   CONTINUE
C
      DO 70 I = 1 , 11
         IF (NINDIF(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(INDIF(I))/INBOT(I)
         END IF
 
 70   CONTINUE
C
      CALL TWRITE(' ')
      CALL TWRITE(' ')
      WRITE (STR,6010)
      CALL TWRITE(STR)
      WRITE (STR,6020)
      CALL TWRITE(STR)
C
      WRITE (STR,6040) (NINDIF(I),I=1,11)
      CALL TWRITE(STR)
      WRITE (STR,6050) (R(I),I=1,11)
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
      CALL TWRITE(' ')
      WRITE (STR,6060)
      CALL TWRITE(STR)
      WRITE (STR,6070)
      CALL TWRITE(STR)
C
      DO 80 I = 1 , 12
         IF (NLDIF(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(LDIF(I))/LBOT(I)
         END IF
 
 80   CONTINUE
C
      WRITE (STR,6080) (NLDIF(I),I=1,12)
      CALL TWRITE(STR)
      WRITE (STR,6090) (R(I),I=1,12)
      CALL TWRITE(STR)
C
      CALL TWRITE('Low spots')
C
      WRITE (STR,6030) (LCOUNTLO(I),I=1,12)
      CALL TWRITE(STR)
C
      CALL TWRITE('High spots')
C
      WRITE (STR,6030) (LCOUNTHI(I),I=1,12)
      CALL TWRITE(STR)
C
      DO 90 I = 1 , 9
         IF (NIDDIF(I).EQ.0) THEN
            R(I) = 0.0
         ELSE
            R(I) = 100.*FLOAT(IDDIF(I))/IDBOT(I)
         END IF
 
 90   CONTINUE
C
      CALL TWRITE(' ')
      WRITE (STR,6100)
      CALL TWRITE(STR)
      WRITE (STR,6110)
      CALL TWRITE(STR)
C
      WRITE (STR,6120) (NIDDIF(I),I=1,9)
      CALL TWRITE(STR)
      WRITE (STR,6130) (R(I),I=1,9)
      CALL TWRITE(STR)
C
      RR = FLOAT(ISUM)/IBOT*100.0
C
      CALL TWRITE(' ')
      WRITE (STR,6140) RR
      CALL TWRITE(STR)
C
      CALL QSBELL
      CALL TREAD(STR,QUIT)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('R-factors on intensity')
 6020 FORMAT ('   200   400   600   800  1000  1200  1400  1600',
     +        ' 1800  2000 >2000')
 6030 FORMAT (12I6)
 6040 FORMAT (11I6)
 6050 FORMAT (11F6.1)
 6060 FORMAT ('R-factors on wavelength')
 6070 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')
 6080 FORMAT (12I6)
 6090 FORMAT (12F6.1)
 6100 FORMAT ('R-factors on dhkl')
 6110 FORMAT ('   >20    10     7     5     4  3.33   2.9   2.6  <2.6')
 6120 FORMAT (9I6)
 6130 FORMAT (9F6.1)
 6140 FORMAT (' Overall R-factor =',F6.1)
      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
      REAL       X, Y, INTENSITIES(6), SIGMAS(6)
      INTEGER    H, K, L
      REAL       LAMBDA
      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)
      REAL       MEAN_INTENS, MEAN_SIG, SIG1, SIG2
      REAL       SUM_OF_WEIGHTS, WEIGHTED_SUM, WEIGHTS(6)
      REAL       SUM_OF_WEIGHTS2, WEIGHTED_RESID2
      REAL       FACTOR
      INTEGER    I, J
      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,6020) SIN_THETA
         CALL TWRITE(STR)
         WRITE (STR,6030) H , K , L
         CALL TWRITE(STR)
         CALL TWRITE(' ')
         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)
         CALL TWRITE('     ** Ignoring this spot **')
         CALL TWRITE(' ')
         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 ('Warning - lambda/2d (=sin theta) =',F5.1,
     +        ' -setting theta to 90')
 6030 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
c - catch normal bad (-9999) sigmas for top film here
c
            IF (SIGMAS(1).LE.0.0) SIGMAS(1) = 1E6
C
      RETURN
      END
**==vicfit.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       VICFIT        **
C*************************
C
C
      SUBROUTINE  VICFIT (RANGE, NSPOT, FILM)
C
C Purpose: To fit Victoreen function to data
C
      INCLUDE 'implicit_none.finc'
C
      EXTERNAL CALFUN
C
      INTEGER    RANGE, NSPOT, FILM
      INTEGER    NSMAX
      PARAMETER (NSMAX=8000)
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION     PARAMS(2,3,5)
C
      INTEGER    IWORKDIM
      PARAMETER (IWORKDIM = 7 + 4*NSMAX)
C
      INTEGER    NRANGE, NFILM
      PARAMETER (NRANGE = 3,NFILM = 6)
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
      INTEGER    LIW
      PARAMETER (LIW=20)
      DOUBLE PRECISION     W(IWORKDIM), FSUMSQ
      INTEGER I
      INTEGER    IW(LIW), LW, IFLAG, REFTYP
      INTEGER    MAXFUN
      DOUBLE PRECISION     DSTEP, DAMP, ACC, FVECC(NSMAX)
      DOUBLE PRECISION     FMIN, SHFLIM
      REAL       TEMP
      LOGICAL    QUIT
      LOGICAL    LENT
      DOUBLE PRECISION     A
C
      DATA       LW,LENT /IWORKDIM,.FALSE./
C
C------------------------------------------------------------------------------
C
      A = PARAMS(2,RANGE,FILM)
C
C
      REFTYP = 1
      CALL TWRITE(' ')
      CALL TASKI('Refinement method (1=LSR, 2=Powell)',REFTYP,QUIT)
      IF (REFTYP.EQ.2) THEN
         ACC = 0.000001
         TEMP = ACC
         CALL TASKF('Acc',TEMP,-3,QUIT)
         ACC = TEMP
         CALL SSQMIN(NSPOT,1,A,ACC,50,CALFUN,FMIN,FVECC,IFLAG,W)
         IF (IFLAG.EQ.-1) THEN
            CALL TWRITE('*One or more failures in line minimisation*')
         ELSE IF (IFLAG.EQ.-2) THEN
            CALL TWRITE('*Refinement not complete after 50 iterations*')
         ELSE IF (IFLAG.GT.0) THEN
            WRITE (STR,6030) IFLAG
            CALL TWRITE(STR)
         END IF
      ELSE
         SHFLIM = 0.001
         TEMP = SHFLIM
         CALL TASKF('Shift limit for convergence',TEMP,-3,QUIT)
         SHFLIM = TEMP
         DSTEP = 1.0E-10
         TEMP = DSTEP
         CALL TASKF('Derivative calculation step',TEMP,-3,QUIT)
         DSTEP = TEMP
         DAMP = 0.5
         TEMP = DAMP
         CALL TASKF('Damping factor',TEMP,-3,QUIT)
         DAMP = TEMP
         CALL LSQMIN(NSPOT,1,A,SHFLIM,50,DSTEP,DAMP,CALFUN,FMIN,
     +               FVECC,IFLAG,W)
         IF (IFLAG.EQ.-1) THEN
            CALL TWRITE('*Refinement not complete after 50 iterations*')
         ELSE IF (IFLAG.EQ.-2) THEN
            CALL TWRITE('**Normal matrix is singular**')
         ELSE IF (IFLAG.GT.0) THEN
            WRITE (STR,6030) IFLAG
            CALL TWRITE(STR)
         END IF
      END IF
C
      PARAMS(2,RANGE,FILM) = A
C
      FSUMSQ = 0.0
      DO 20 I = 1 , NSPOT
         FSUMSQ = FSUMSQ + FVECC(I)*FVECC(I)
 20   CONTINUE
C
      CALL TWRITE(' ')
      WRITE (STR,6010) FSUMSQ
      CALL TWRITE(STR)
      WRITE (STR,6020) (PARAMS(I,RANGE,FILM),I=1,2)
      CALL TWRITE(STR)
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Final sum-of-squares  ',D10.4)
 6020 FORMAT ('Final parameters: K= ',D10.4,' a= ',D10.4)
 6030 FORMAT ('Refinement completed after',I4,' iterations')
      END
**==victoreen.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**       VICTOREEN        **
C****************************
C
C
      DOUBLE PRECISION FUNCTION  VICTOREEN (K, A, DWAV)
C
C Purpose: Applies Victoreen function given the above parameters
C          Note that sec(2*theta) correction NOT applied
C
      INCLUDE 'implicit_none.finc'
C
      DOUBLE PRECISION     K, A, DWAV
C
C====== LOCALS:
C
      CHARACTER*80 STR
      DOUBLE PRECISION     DWAV3, ARG
C
C------------------------------------------------------------------------------
C
      DWAV3 = DWAV*DWAV*DWAV
      ARG = A*DWAV3
      IF (ARG.GT.80.0) THEN
         WRITE (STR,6010) ARG
         CALL TWRITE(STR)
         ARG = 80.0
      END IF
C
      VICTOREEN = K*EXP(ARG)
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Arg =',D8.2)
      END
**==weightfun.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**       WEIGHTFUN        **
C****************************
C
C
      REAL   FUNCTION  WEIGHTFUN (AVAL2,BVAL2,ASIG2,BSIG2)
C
C Purpose: Weighting function
C
      INCLUDE 'implicit_none.finc'
C
      REAL    AVAL2, BVAL2, ASIG2, BSIG2
C
C------------------------------------------------------------------------------
C
      WEIGHTFUN = (ASIG2/AVAL2+BSIG2/BVAL2)*AVAL2/BVAL2
C
      RETURN
      END
**==write_film_factors.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************************
C**       WRITE_FILM_FACTORS        **
C*************************************
C
C
      SUBROUTINE WRITE_FILM_FACTORS
C
C Purpose: To write out Victoreen coeffs to a file on unit 10
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FILM
C
      INTEGER    NPARS,NRANGE,NFILM
      PARAMETER (NPARS=2,NRANGE=3,NFILM=6)
C
      COMMON /FITTING_PARAMS/  PARAMS
      DOUBLE PRECISION     PARAMS(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
      INTEGER LENSTR
C
C====== LOCALS:
C
      CHARACTER*80 STR
      LOGICAL QUIT
      CHARACTER*60 FNAM, FILENAME
      INTEGER    I, J, K, IFAIL
C
C------------------------------------------------------------------------------
C
 20   CALL TPROMP('File name (default ext=.vc): ')
      CALL TREAD(FNAM,QUIT)
      CALL FILEXT(FNAM,'vc ',FILENAME,I)
      WRITE (STR,6010) FILENAME(1:LENSTR(FILENAME))
      CALL TWRITE(STR)
C
      IFAIL = 1
      CALL CCPDPN(-10,FILENAME,'NEW','F',0,IFAIL)
      IF (IFAIL.EQ.-1) GO TO 20
C
      DO 50 K = 1 , NFILM - 1
         DO 40 J = 1 , NRANGE
            DO 30 I = 1 , NPARS
               WRITE (10,6020) PARAMS(I,J,K)
 30         CONTINUE
 40      CONTINUE
 50   CONTINUE
C
      CLOSE (10)
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Writing coeffs to ',A)
 6020 FORMAT (1X,G14.6)
      END