**==aa0001.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**  Program GENLAUE    **
C*************************
C
C
C Program: GENLAUE
C
C Purpose: Allows user to enter and refine the experimental parameters
C          that govern a Laue pattern, and also produces the generate file.
C
C The Program has several steps...
C
C      1. Read in initial parameters from .gen file
C      2. Predict pattern based on given parameters
C      3. Find energy and spatial overlaps in prediction, i.e. find nodals
C      4. Find nodals on film and match with predicted nodal spots
C      5. Using the information from step 4 calculate refined parameters
C      6. Generate refined prediction
C      7. Find energy and spatial overlaps in refined prediction
C      8. Do a vector plot showing deviation of predicted from observed spots
C      9. If user is not satisfied then goto 4
C      10. Write out .ge1/2 files
C      11. Write out refined parameters to a new .gen file
C      12. Exit
C
C------------------------------------------------------------------------------
C
      INCLUDE 'implicit_none.finc'
C
C======  INCLUDE this file 'cos it's different for small & large versions.
C
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'afilmspots.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'gen_file_header.finc'
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'refined_film_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'observed_spots.finc'
      INCLUDE 'backgr.finc'
      INCLUDE 'ramask.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      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
      INTEGER    LUNSTI, LUNSTO
C
C====== LOCALS:
C
      LOGICAL EOF
      CHARACTER*80 STR
      CHARACTER*20 TRM
      INTEGER      ITTYP
      CHARACTER*40  FILM_NAME, GE_NAME, FNAME
      INTEGER END_POS
      REAL       WAVSTD, WAVMIN, WAVMAX
      REAL       DMIN, FILM_RAD
      LOGICAL    QUIT, REFINE, HAPPY, MAP, NO_FILM, RERUN, ANOTHER
      LOGICAL    REGEN, DONE_A
      REAL       DELTA, EPSILON, XSTART, XEND, BOX_SIZE
      REAL       X_CEN_FS(6), Y_CEN_FS(6), W_FS(6)
      REAL       X_REJECT, Y_REJECT
      INTEGER    I, FILM_NUM, LAMBDA_BINS(16)
      INTEGER    IQUIT, IFAIL
      INTEGER    NUM_OLD_OBS, NUM_REF_SPOTS, MININT
      REAL       THICK
      INTEGER    IMARK
C
      DATA   FILM_NUM /1/
      DATA   THICK /.22/
      DATA   REGEN, DONE_A /.TRUE., .FALSE./
C
C------------------------------------------------------------------------------
C
C
C======  Initialise graphics & terminal streams
C
      CALL TRWSET(LUNSTI(0),LUNSTO(0),8,8,0)
      ITTYP = 0
      CALL UGTENV('TERM',TRM)
      IF (TRM(1:5).EQ.'xterm') ITTYP = 1
      CALL QGINIT(0,0,ITTYP,LUNSTI(0),LUNSTO(0))
      IN_TERM = LUNSTI(0)
      OUT_TERM = LUNSTO(0)
      IN_DATA = 2
      OUT_DATA = 1
C
C======  Log file, please
C
      MAX_REFLEX = MAX_REF
      IPLATE = .FALSE.
      IFAIL = 0
      CALL CCPDPN(-8,'GENLOG','NEW','F',0,IFAIL)
C
      CALL TWRITE(' ')
      CALL TWRITE('      GENLAUE: Laue Image Analysis Software')
      CALL TWRITE('      -------------------------------------')
      CALL TWRITE('     S. Zurek, I.J. Clifton and P.A. Machin')
      CALL TWRITE('           SERC Daresbury Laboratory')
      CALL TWRITE(' ')
      CALL TWRITE('       Latest version October 3''rd, 1992')
      CALL TWRITE('   including elliptical masking and IP software')
      CALL TWRITE('   J.W. Campbell, T.J. Greenhough & A.K. Shrive')
      CALL TWRITE(' ')
C
C====== Get starting parameters
C
      MAX_OBSERVED = MAX_OBS
      CALL GET_GEN_PARAMTERS(WAVSTD,WAVMIN,WAVMAX,DMIN,FILM_RAD,
     +                       X_REJECT,Y_REJECT)
C
C======= Initialise fiducial parameters to zero
C
      DO 20 I = 1 , 6
         X_CEN_FS(I) = 0.0
         Y_CEN_FS(I) = 0.0
         W_FS(I) = 0.0
 20   CONTINUE
C
C====== Initialize film parameters to be refined
C
      DO 30 I = 1 , 6
         X_C_S(I) = X_C
         Y_C_S(I) = Y_C
         W_C_S(I) = W_C
         TWISTS(I) = TWIST
         TILTS(I) = TILT
         BULGES(I) = BULGE
         Y_SCALES(I) = Y_SCALE
         CTOF_S(I) = CTOF
 30   CONTINUE
C
C====== Predict pattern
C
      DELTA = 0.3
      EPSILON = 0.05
      CALL CALC_LAUE(WAVSTD,WAVMIN,WAVMAX,DMIN,FILM_RAD,DELTA,EPSILON,
     +               LAMBDA_BINS,FILM_NUM)
C
C====== Does user want to skip refinement ?
C
      REFINE = .TRUE.
      CALL TASKQ('Do you want refinement option',REFINE,EOF)
      IF (.NOT.REFINE) GO TO 70
C
C====== Start of refinement loop
C====== ++++++++++++++++++++++++
C
 40   NO_FILM = .FALSE.
      CALL TASKQ('Do you want to refine with a saved file of matches',
     +           NO_FILM,EOF)
C
C====== Locate nodals on film and match with prediction
C
      IF (.NOT.NO_FILM) THEN
         CALL GET_FILM_INFO(FILM_NAME,FILM_NUM)
         IF (FILM_NUM.EQ.1) DONE_A = .TRUE.
      ELSE
C
C====== If no film, treat as standard 2400sq 50u
C
         RASTER_SIZE = 0.05
         X_RASTERS = 2400
         Y_RASTERS = 2400
         X_CEN_F = 1200
         Y_CEN_F = 1200
         FILM_NUM = 1
      END IF
C
 50   CALL MATCH_NODALS(DELTA,FILM_NAME,BOX_SIZE,NUM_OLD_OBS,
     +                  NUM_REF_SPOTS,X_REJECT,Y_REJECT,NO_FILM)
C
C====== Now refine parameters
C
      CALL REFINE_GEN_PARAMETERS(WAVSTD,WAVMIN,WAVMAX,FILM_RAD,
     +                           NUM_REF_SPOTS,FILM_NUM)
C
C====== Plot vector shifts from observed to predicted position
C
 60   IF (NO_FILM) THEN
         MAP = .FALSE.
      ELSE
         MAP = .TRUE.
         CALL TASKQ('Do you want a plot of the difference map',MAP,EOF)
      END IF
      IF (MAP) THEN
         CALL PLOT_DIFF(FILM_NAME,DELTA,BOX_SIZE,NUM_OLD_OBS)
      END IF
C
C====== Give option to re-run
C
      RERUN = .TRUE.
      CALL TASKQ('Do you want to re-run refinement immediately',RERUN,
     +           EOF)
      IF (RERUN) THEN
      CALL REFINE_GEN_PARAMETERS(WAVSTD,WAVMIN,WAVMAX,FILM_RAD,
     +                           NUM_REF_SPOTS,FILM_NUM)
         GO TO 60
      END IF
C
C====== Generate refined prediction
C
      IF (REGEN) THEN
         CALL CALC_LAUE(WAVSTD,WAVMIN,WAVMAX,DMIN,FILM_RAD,DELTA,
     +                  EPSILON,LAMBDA_BINS,FILM_NUM)
      END IF
C
C====== Put fiducial parameters into this film's slot
C====== and convert to mm. units.
C
      X_CEN_FS(FILM_NUM) = X_CEN_F/MM_RAST_X
      Y_CEN_FS(FILM_NUM) = Y_CEN_F/MM_RAST_Y
      W_FS(FILM_NUM) = W_F
C
      IF (.NOT.DONE_A) THEN
         CALL TWRITE(
     +   ' WARNING - you have not yet defined the A-film spot positions'
     +   )
         HAPPY = .FALSE.
      ELSE
         HAPPY = .TRUE.
      END IF
      CALL TASKQ('Are you happy with the refinement for this film',
     +           HAPPY,EOF)
      IF (.NOT.HAPPY) THEN
         CALL TASKQ('Do you want to refine with a saved file of matches'
     +              ,NO_FILM,EOF)
         GO TO 50
      ELSE IF (.NOT.IPLATE) THEN
         CALL TASKQ('Do you want to refine with another film',ANOTHER,
     +              EOF)
         IF (ANOTHER) THEN
            CALL TASKF('Enter c-to-f increment to next film (mm)',THICK,
     +                 3,QUIT)
            CTOF = CTOF + THICK
            STREAK = .FALSE.
            GO TO 40
         END IF
      END IF
C
C====== End of refinement loop
C====== ----------------------
C
C
C====== Get filename from user and open the relevant files
C
 70   CALL TPROMP('Enter root name for output .gen/.ge1/.ge2 files: ')
      CALL TREAD(FNAME,EOF)
      IF (FNAME(1:1).EQ.' ') STOP
      CALL FILEXT(FNAME,'gen',GE_NAME,END_POS)
C
C====== Recalculate pattern for A film if lower films refined
C====== (Note: Storing A film coords in xf_a, yf_a was not adequate)
C
      IF (FILM_NUM.NE.1) THEN
         CTOF = CTOF_S(1)
         FILM_NUM = 1
         CALL CALC_LAUE(WAVSTD,WAVMIN,WAVMAX,DMIN,FILM_RAD,DELTA,
     +        EPSILON,LAMBDA_BINS,FILM_NUM)
      END IF
C
C====== Now write out refined parameters
C
      CALL WRITE_GEN_PARAMETERS(GE_NAME(:END_POS),WAVSTD,WAVMIN,WAVMAX,
     +                          DMIN,FILM_RAD,X_REJECT,Y_REJECT)
C
C====== And finally write out to a .ge1/2 file
C
C
C====== Set a sensible default for "min intensity to measure on subs. films"
C
      MININT = 100
      IF (.NOT.IPLATE) CALL TASKI(
     +         'Enter Min I to by-pass measurements on subsequent films'
     +         ,MININT,EOF)
      IF (GE_NAME.NE.' ') THEN
         CALL WRITE_GE_FILE(GE_NAME(:END_POS),NUM_REFLEX,H_S,K_S,L_S,
     +                      M_FLAGS,LAMBDAS,XF_A,YF_A,DMIN_THRESHOLDS2,
     +                      TITLE,CRYSTAL_ID,CELL_PARAMS,
     +                      NUM_CELL_PARAMS,XTAL_SYSTEM,LATTICE,X_AXIS,
     +                      Z_AXIS,PHIS,WAVSTD,WAVMIN,WAVMAX,DMIN,
     +                      FILM_RAD,SPINDLE,DELTA,EPSILON,MININT,
     +                      RASTER_SIZE,X_CEN_FS,Y_CEN_FS,W_FS,X_C_S,
     +                      Y_C_S,W_C_S,TWISTS,TILTS,BULGES,Y_SCALES,
     +                      CTOF_S,LAMBDA_BINS,X_REJECT,Y_REJECT)
      END IF
C
      CALL QGCLOS
      CALL TWRITE('*Genlaue completed*')
      STOP
      END
**==a_calc.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       A_CALC        **
C*************************
C
C
      SUBROUTINE  A_CALC (ORIENTATION, PHIS, A_MAT)
C
C Purpose: Calculates the master matrix, a_mat, from the crystal system,
C          crstal alignment and cell parameters.
C          This is then premultiplied by the phi matrices.
C Optimisation: Don`t re-calculate if orientation & phis are the same from
C               call to call.
C
      INCLUDE 'implicit_none.finc'
C
      REAL       ORIENTATION(3,3)
      REAL       PHIS(3),  A_MAT(3,3)
C
C====== LOCALS:
C
      REAL    OLD_ORIENT(3,3), OLD_PHIS(3), OLD_A_MAT(3,3)
      LOGICAL CALLED_ONCE
      REAL    PC(3,3)
      INTEGER   I, J
      SAVE    OLD_ORIENT, OLD_PHIS, CALLED_ONCE, OLD_A_MAT
      DATA    CALLED_ONCE/.FALSE./
C
C------------------------------------------------------------------------------
C
      IF (CALLED_ONCE) THEN
C
C====== See if phis have changed - return if not
C
         DO 20 I = 1 , 3
            IF (PHIS(I).NE.OLD_PHIS(I)) GO TO 70
 20      CONTINUE
C
C====== Better check orientation as well
C
         DO 40 I = 1 , 3
            DO 30 J = 1 , 3
               IF (ORIENTATION(I,J).NE.OLD_ORIENT(I,J)) GO TO 70
 30         CONTINUE
 40      CONTINUE
C
C====== Return old a_mat
C
         DO 60 I = 1 , 3
            DO 50 J = 1 , 3
               A_MAT(I,J) = OLD_A_MAT(I,J)
 50         CONTINUE
 60      CONTINUE
         RETURN
      END IF
C
 70   CALLED_ONCE = .TRUE.
      DO 90 I = 1 , 3
         OLD_PHIS(I) = PHIS(I)
         DO 80 J = 1 , 3
            OLD_ORIENT(I,J) = ORIENTATION(I,J)
 80      CONTINUE
 90   CONTINUE
C
C====== Define phi rotation matrix
C
      CALL MTX_ROTATE(PHIS,PC)
C
C====== Apply phi rotations
C
      CALL MTX_MULT(PC,ORIENTATION,A_MAT,3,3,3)
C
C====== Save old a_mat
C
      DO 110 I = 1 , 3
         DO 100 J = 1 , 3
            OLD_A_MAT(I,J) = A_MAT(I,J)
 100     CONTINUE
 110  CONTINUE
      RETURN
      END
**==bkset.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       BKSET         **
C*************************
C
C
      SUBROUTINE BKSET(IX)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose : For a given record set yvalues where background
C           changes (in ibgy) and store the background to
C           use in ibgv. note there are usually an odd
C           number of regions, except when ix=x_cen_f
C
      INTEGER IX
C
      INCLUDE 'backgr.finc'
C
      COMMON /CMBACK/IBGY,IBGV
      INTEGER IBGY(30),IBGV(30)
C
      INCLUDE 'fiducials.finc'
C
      COMMON /PLOTS/SCAN_START_X,SCAN_START_Y,IDX,IDY
      INTEGER IDX,IDY,SCAN_START_X,SCAN_START_Y
C
C Author: IJC
C
C====== LOCALS:
C
      INTEGER RAD,RAD_MAX,Y_BIN,I,NUM_SEEN
      REAL Y_DIFF
C
C-------------------------------------------------------------------------------
C
      RAD_MAX = Y_RANGE/2
      DO 20 I = 1 , 30
         IBGY(I) = 10000
         IBGV(I) = BACKGR(1)
 20   CONTINUE
      RAD = RAD_MAX - N_PER_BIN
      I = 0
C
C====== See if we cut the next background region
C
 30   IF (IABS(IX-XCEN).LT.RAD) THEN
         I = I + 1
         Y_DIFF = RAD**2 - (IX-XCEN)**2
         Y_BIN = YCEN - SQRT(Y_DIFF)
         IBGY(I) = Y_BIN
         IBGV(I) = BACKGR(I)
C
C====== Seen all the regions?
C
         IF (I.LT.N_BACKGR) THEN
            RAD = RAD - N_PER_BIN
            GO TO 30
         END IF
      END IF
 
      NUM_SEEN = I
      IF (NUM_SEEN.LT.1) RETURN
C
      IBGV(NUM_SEEN+1) = BACKGR(I+1)
C
C====== Apply symmetry
C
      DO 40 I = 1 , NUM_SEEN
         IBGV(NUM_SEEN+I+1) = IBGV(NUM_SEEN-I+1)
         IBGY(NUM_SEEN+I) = 2*YCEN - IBGY(NUM_SEEN-I+1)
 40   CONTINUE
      RETURN
      END
**==c_of_g.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       C_OF_G        **
C*************************
C
C
      SUBROUTINE  C_OF_G (FILM_UNIT, XC, YC, X_SIZE, Y_SIZE, THRESH,
     +                    CGX, CGY, NUM_OVER)
C
C Purpose: Finds c of g of intensities lying in box centred on (xc, yc)
C          of dimensions 2*x_size+1 and 2*y_size+1,
C          allowing for a constant intensity threshold
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FILM_UNIT
C
C====== Integers because in SCANNER space
C
      INTEGER    XC, YC, X_SIZE, Y_SIZE
      INTEGER    THRESH, NUM_OVER
      REAL       CGX, CGY
C
C====== NUM_OVER: number of spots in the box over the threshold.
C
      INCLUDE 'sizes.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      INCLUDE 'ramask.finc'
C
C====== LOCALS:
C
      LOGICAL    ERR
      CHARACTER*80 STR
      INTEGER    IX, IY, NSTRIP
      INTEGER    ISTRIP(MAX_RASTERS)
      INTEGER    INTENS, I
      REAL       SUMX, SUMY, SUM_INTENS, DX, DY
      INTEGER IMARK
C
C------------------------------------------------------------------------------
C
c
      NUM_OVER = 0
      SUMX = 0.0
      SUMY = 0.0
      SUM_INTENS = 0.0
C
      NSTRIP = XC - X_SIZE
      IF (XC-X_SIZE.LT.1 .OR. XC+X_SIZE.GT.X_RASTERS .OR. 
     +    YC-Y_SIZE.LT.1 .OR. YC+Y_SIZE.GT.Y_RASTERS) THEN
         WRITE (STR,6010) XC , YC
         CALL TWRITE(STR)
         NUM_OVER = 0
         CGX = XC
         CGY = YC
         RETURN
      END IF
      IF (IPLATE) THEN
         CALL IMPPOS(FILM_UNIT,NSTRIP,OUT_TERM,ERR)
      ELSE
         CALL FLMPOS(FILM_UNIT,NSTRIP,OUT_TERM,ERR)
      END IF
C
      DO 30 IX = -X_SIZE , X_SIZE
         IF (IPLATE) THEN
            CALL IMPLRD(FILM_UNIT,ISTRIP,Y_RASTERS,OUT_TERM,ERR)
         ELSE
            CALL FILMRD(FILM_UNIT,ISTRIP,Y_RASTERS,OUT_TERM,ERR)
         END IF
         NSTRIP = NSTRIP + 1
C
         DO 20 IY = -Y_SIZE , Y_SIZE
C
C====== Call radmask for streaked spots
C
            IMARK = 1
            IF (STREAK) THEN
               CALL RADMASK(1,XC,YC,NSTRIP,IY+YC,9111,IMARK)
            END IF
            IF (IMARK.EQ.1) THEN
               INTENS = ISTRIP(YC+IY)
               INTENS = MAX(0,INTENS-THRESH)
               IF (INTENS.GT.0) NUM_OVER = NUM_OVER + 1
               SUMX = SUMX + IX*INTENS
               SUMY = SUMY + IY*INTENS
               SUM_INTENS = SUM_INTENS + INTENS
            END IF
 20      CONTINUE
 
 30   CONTINUE
C
      IF (SUM_INTENS.EQ.0.0) THEN
         DX = 0
         DY = 0
      ELSE
         DX = SUMX/SUM_INTENS
         DY = SUMY/SUM_INTENS
      END IF
C
      CGX = XC + DX
      CGY = YC + DY
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('**Attempt to read beyond scanned film** xc =',I5,' yc =',
     +        I5)
      END
**==calc_backgr.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       CALC_BACKGR        **
C******************************
C
C
      SUBROUTINE  CALC_BACKGR (FILM_UNIT)
C
C Purpose: Estimates a radial background curve
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FILM_UNIT
C
      INCLUDE 'backgr.finc'
      INCLUDE 'sizes.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      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    ERR
      CHARACTER*80 STR
      INTEGER    I, FREQ(1000), IVAL(1000), NFREQ, IBIN
      INTEGER    X, Y
      INTEGER    BIN, MAX_IN_BIN
      INTEGER    N_INTENS, INTENS, SUM_INTENS(N_BINS), SUM
      INTEGER NUM_PLOT_X, NUM_PLOT_Y, IDY
      PARAMETER (NUM_PLOT_X=700,NUM_PLOT_Y=700)
      INTEGER    ISTRIP(MAX_RASTERS)
C
C------------------------------------------------------------------------------
C
C====== Calculate the approx centre used for background estimation
C
      XCEN = X_RASTERS/2
      YCEN = Y_RASTERS/2
      NFREQ = 256
      IF (IPLATE) NFREQ = 1000
C
      IDY = Y_RASTERS/NUM_PLOT_Y
      Y_RANGE = IDY*NUM_PLOT_Y
      Y_START = YCEN - Y_RANGE/2
      N_PER_BIN = Y_RANGE/N_BINS
C
      CALL TWRITE('=== Estimating background...')
C
      MAX_IN_BIN = 0.7*N_PER_BIN + 1
      DO 20 BIN = 1 , N_BINS
         SUM_INTENS(BIN) = 0
 20   CONTINUE
C
      IF (IPLATE) THEN
         CALL IMPPOS(FILM_UNIT,XCEN-20,OUT_TERM,ERR)
      ELSE
         CALL FLMPOS(FILM_UNIT,XCEN-20,OUT_TERM,ERR)
      END IF
 
      DO 80 X = XCEN - 20 , XCEN + 20
         IF (IPLATE) THEN
            CALL IMPLRD(FILM_UNIT,ISTRIP,Y_RASTERS,OUT_TERM,ERR)
         ELSE
            CALL FILMRD(FILM_UNIT,ISTRIP,Y_RASTERS,OUT_TERM,ERR)
         END IF
         Y = Y_START
C
         DO 70 BIN = 1 , N_BINS
            DO 30 I = 1 , NFREQ
               FREQ(I) = 0
               IVAL(I) = 0
 30         CONTINUE
            DO 40 I = 1 , N_PER_BIN
               INTENS = ISTRIP(Y)
               IF (IPLATE) THEN
                  IBIN = INTENS/1000 + 1
                  IF (IBIN.GT.NFREQ) IBIN = NFREQ
               ELSE
                  IBIN = INTENS + 1
               END IF
               FREQ(IBIN) = FREQ(IBIN) + 1
               IVAL(IBIN) = IVAL(IBIN) + INTENS
               Y = Y + 1
 40         CONTINUE
            N_INTENS = 0
            SUM = 0
            DO 50 I = 1 , NFREQ
               N_INTENS = N_INTENS + FREQ(I)
               SUM = SUM + IVAL(I)
               IF (N_INTENS.GT.MAX_IN_BIN) GO TO 60
 50         CONTINUE
 60         SUM_INTENS(BIN) = SUM_INTENS(BIN) + SUM/N_INTENS
 70      CONTINUE
 80   CONTINUE
C
      DO 90 BIN = 1 , N_BINS
C
C====== 41 is range -20..20
C
         SUM_INTENS(BIN) = SUM_INTENS(BIN)/41
 90   CONTINUE
C
C====== Adjust edges
C
      DO 100 I = 1 , N_BINS/2
         BACKGR(I) = MAX(SUM_INTENS(I),SUM_INTENS(N_BINS+1-I))
 100  CONTINUE
 
      BACKGR(1) = BACKGR(2)
      I = N_BINS/2
      BACKGR(I) = BACKGR(I-1) + IABS(BACKGR(I-1)-BACKGR(I-2))
 
      WRITE (STR,6010) (BACKGR(I),I=1,N_BINS/2)
      CALL TWRITE(STR)
      RETURN
C
C====== Format statements
C
 6010 FORMAT (20I4)
      END
**==calc_laue.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       CALC_LAUE     **
C*************************
C
C
      SUBROUTINE  CALC_LAUE (WAVSTD, WAVMIN, WAVMAX,
     +                       DMIN, FILM_RAD, DELTA, EPSILON,
     +                       LAMBDA_BINS, FILM_NUM)
C
C Purpose: To generate a Laue pattern and calculate overlaps and nodals
C
      INCLUDE 'implicit_none.finc'
C
      REAL       WAVSTD, WAVMIN, WAVMAX
      REAL       DMIN, FILM_RAD, DELTA, EPSILON
      INTEGER    LAMBDA_BINS(16), FILM_NUM
C
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'afilmspots.finc'
C
C====== LOCALS:
C
      REAL       XSTART, XEND
      INTEGER    I
      LOGICAL    QUIT
C
C------------------------------------------------------------------------------
C
C====== Allow user to change dmin if need be
C
      CALL TASKF('Dmin for this prediction',DMIN,2,QUIT)
C
C====== Generate new prediction
C
      CALL GENERATE_REFLEXIONS(WAVSTD,WAVMIN,WAVMAX,DMIN,FILM_RAD,
     +                         XSTART,XEND)
C
C====== Save A film (x,y) to go into .ge1/.ge2 file.
C
      IF (FILM_NUM.EQ.1) THEN
         DO 20 I = 1 , NUM_REFLEX
            XF_A(I) = XF_S(I)
            YF_A(I) = YF_S(I)
 20      CONTINUE
      END IF
C
      CALL LAMBDA_HIST(NUM_REFLEX,LAMBDAS,WAVSTD,WAVMIN,WAVMAX,
     +                 LAMBDA_BINS)
C
C====== Find new overlap spots
C
      CALL TASKF('Enter delta - the spatial overlap criterion (mm)',
     +           DELTA,2,QUIT)
      CALL TASKF(
     +        'Enter epsilon - the spatial overlap closeness limit (mm)'
     +        ,EPSILON,2,QUIT)
 
      CALL FIND_OVERLAPS(DELTA,EPSILON,XSTART,XEND)
 
      RETURN
      END
**==calc_xym.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       CALC_XYM      **
C*************************
C
C
      SUBROUTINE  CALC_XYM (H, K, L, XM, YM)
C
C Purpose: Calculates (xm,ym) given h,k,l values and current parameter values
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    H, K, L
      REAL       XM, YM
C
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
C
C====== LOCALS:
C
      REAL       A(3,3), SPIN(3,3)
      REAL       X, Y, Z
      REAL       XF, YF
      REAL       X1, Y1, Z1
      REAL       X_CEN, Y_CEN, OMEGA, Y_SCRAST
      REAL       COS_OMEGA, SIN_OMEGA
      REAL       RF, K_FACTOR, DEG_TO_RAD
C
C------------------------------------------------------------------------------
C
      DEG_TO_RAD = 3.1415926536/180.0
 
      X_CEN = X_CEN_F + X_C
      Y_CEN = Y_CEN_F + Y_C
      OMEGA = W_F + W_C
 
      COS_OMEGA = COS(OMEGA*DEG_TO_RAD)
      SIN_OMEGA = SIN(OMEGA*DEG_TO_RAD)
 
      Y_SCRAST = MM_RAST_Y
 
      CALL MAKE_ORIENT(CELL_PARAMS,XTAL_SYSTEM,X_AXIS,Z_AXIS,
     +                 ORIENTATION)
 
      CALL A_CALC(ORIENTATION,PHIS,A)
C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
      X1 = A(1,1)*H + A(1,2)*K + A(1,3)*L
      Y1 = A(2,1)*H + A(2,2)*K + A(2,3)*L
      Z1 = A(3,1)*H + A(3,2)*K + A(3,3)*L
C
C====== Apply spindle rotation
C
      X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
      Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
      Z = Z1
 
      RF = (X*X+Y*Y+Z*Z)/(2*X)
 
      XF = CTOF*Y/(X-RF)
      YF = CTOF*Z/(X-RF)
C
C====== K_factor = 1.0 + tilt*xf + twist*yf + bulge*sqrt(xf*xf + yf*yf)
C
      K_FACTOR = 1.0
 
      XM = X_CEN + K_FACTOR*(XF*COS_OMEGA-YF*SIN_OMEGA)*MM_RAST_X
      YM = Y_CEN + K_FACTOR*(XF*SIN_OMEGA+YF*COS_OMEGA)*Y_SCRAST
 
      RETURN
      END
**==calfun.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       CALFUN        **
C*************************
C
C
      SUBROUTINE  CALFUN (M, N, FVECC, XC)
C
C Purpose: Harwell routine uses this in refining parameters
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    M, N
      DOUBLE PRECISION     XC(N), FVECC(M)
C
      COMMON /CALFUN_PARAMS/ SCALING, XC_START, FIRST
      DOUBLE PRECISION      SCALING(20), XC_START(20)
      LOGICAL    FIRST
 
      COMMON /RFLG/ REFLAG(20), RFCL
      INTEGER       REFLAG
      LOGICAL       RFCL
      SAVE /RFLG/
C
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
C
C====== LOCALS:
C
      INTEGER    I, J
      LOGICAL    REFINE_CELL
      INTEGER    H, K, L
      DOUBLE PRECISION     A(3,3), SPIN(3,3)
      DOUBLE PRECISION     X, Y, Z
      DOUBLE PRECISION     XF, YF
      DOUBLE PRECISION     XM, YM
      DOUBLE PRECISION     X1, Y1, Z1
      DOUBLE PRECISION     X_CEN, Y_CEN, OMEGA, Y_SCRAST
      DOUBLE PRECISION     COS_OMEGA, SIN_OMEGA
      DOUBLE PRECISION     RF, K_FACTOR, DEG_TO_RAD
      DOUBLE PRECISION     PC(3,3), C(3,3), CX(3,3), CY(3,3), CZ(3,3)
      DOUBLE PRECISION     CPX, CPY, CPZ, SPX, SPY, SPZ
C
C====== Save in common to make sure its preserved for successive calls
C
      DOUBLE PRECISION     LOC_ORIENT(3,3)
      COMMON/SSAVE/ LOC_ORIENT
C
      DOUBLE PRECISION     RL(3), CC(3), SS(3)
      INTEGER    SI, SJ, SK, ILOOP, JLOOP
      DOUBLE PRECISION     COS_PSI, COS_RHO, A_AXIS
      DOUBLE PRECISION     COS_ALPHA, SIN_ALPHA, FUN
      INTEGER    IANG
C
C====== Local double precision copies of the params being refined.
C
      DOUBLE PRECISION     LOC_CELL_PARAMS(6)
      DOUBLE PRECISION     LOC_PHIS(3)
      DOUBLE PRECISION     LOC_CTOF
      DOUBLE PRECISION     LOC_X_C
      DOUBLE PRECISION     LOC_Y_C
      DOUBLE PRECISION     LOC_W_C
C
C------------------------------------------------------------------------------
C
      DEG_TO_RAD = 3.1415926536/180.0
 
      REFINE_CELL = RFCL
C
C====== Unpack parameters and correct for scaling
C
      DO 20 I = 1 , 3
         IF (REFLAG(I).GT.0) THEN
            LOC_PHIS(I) = XC(REFLAG(I))*SCALING(REFLAG(I)) + XC_START(I)
         ELSE
            LOC_PHIS(I) = XC_START(I)
         END IF
 20   CONTINUE
      IF (REFLAG(4).GT.0) THEN
         LOC_CTOF = XC(REFLAG(4))*SCALING(REFLAG(4)) + XC_START(4)
      ELSE
         LOC_CTOF = XC_START(4)
      END IF
      IF (REFLAG(5).GT.0) THEN
         LOC_X_C = XC(REFLAG(5))*SCALING(REFLAG(5)) + XC_START(5)
      ELSE
         LOC_X_C = XC_START(5)
      END IF
      IF (REFLAG(6).GT.0) THEN
         LOC_Y_C = XC(REFLAG(6))*SCALING(REFLAG(6)) + XC_START(6)
      ELSE
         LOC_Y_C = XC_START(6)
      END IF
      IF (REFLAG(7).GT.0) THEN
         LOC_W_C = XC(REFLAG(7))*SCALING(REFLAG(7)) + XC_START(7)
      ELSE
         LOC_W_C = XC_START(7)
      END IF
 
      IF (REFINE_CELL) THEN
         DO 30 I = 8 , 7 + NUM_CELL_PARAMS
            IF (REFLAG(I).GT.0) THEN
               LOC_CELL_PARAMS(I-7) = XC(REFLAG(I))*SCALING(REFLAG(I))
     +                                + XC_START(I)
            ELSE
               LOC_CELL_PARAMS(I-7) = XC_START(I)
            END IF
 30      CONTINUE
      ELSE
         DO 40 I = 1 , 6
            LOC_CELL_PARAMS(I) = CELL_PARAMS(I)
 40      CONTINUE
      END IF
C
C====== Now prepare to evaluate residual
C
      X_CEN = X_CEN_F + LOC_X_C
      Y_CEN = Y_CEN_F + LOC_Y_C
      OMEGA = W_F + LOC_W_C
 
      COS_OMEGA = COS(OMEGA*DEG_TO_RAD)
      SIN_OMEGA = SIN(OMEGA*DEG_TO_RAD)
 
      Y_SCRAST = MM_RAST_Y
 
      IF (REFINE_CELL .OR. FIRST) THEN
C
C====== In-line, R*8 version of MAKE_ORIENT--------------------
C====== i, j, k have special meanings in this piece of code.
C====== i refers to the crystal rotation axis, Z
C
         FIRST = .FALSE.
         I = IABS(Z_AXIS)
         SI = 1
         IF (Z_AXIS.LT.0) SI = -1
C
C====== j refers to X_axis, the x-ray beam axis, X
C
         J = IABS(X_AXIS)
         SJ = 1
         IF (X_AXIS.LT.0) SJ = -1
C
C====== k is therefore the other axis, Y
C
         K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
         IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
            SK = SI*SJ
         ELSE
            SK = -SI*SJ
         END IF
C
C====== Zero orientation matrix
C
         DO 60 ILOOP = 1 , 3
            DO 50 JLOOP = 1 , 3
               LOC_ORIENT(ILOOP,JLOOP) = 0.
 50         CONTINUE
 60      CONTINUE
C
C====== Convert cell angles to radians
C
         IANG = 4
C
C====== Check if rhombohedral
C
         IF (XTAL_SYSTEM.EQ.6) IANG = 2
C
         DEG_TO_RAD = 3.1415926536/180.0
         LOC_CELL_PARAMS(IANG) = LOC_CELL_PARAMS(IANG)*DEG_TO_RAD
         LOC_CELL_PARAMS(IANG+1) = LOC_CELL_PARAMS(IANG+1)*DEG_TO_RAD
         LOC_CELL_PARAMS(IANG+2) = LOC_CELL_PARAMS(IANG+2)*DEG_TO_RAD
C
C====== Now go and do specific stuff for each crystal system
C
         IF (XTAL_SYSTEM.EQ.2) THEN
C
C====== Monoclinic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            CC(1) = 0.0
            CC(2) = COS(LOC_CELL_PARAMS(4))
            CC(3) = 0.0
 
            SS(1) = 1.0
            SS(2) = SIN(LOC_CELL_PARAMS(4))
            SS(3) = 1.0
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*CC(I)
            LOC_ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         ELSE IF (XTAL_SYSTEM.EQ.3) THEN
C
C====== Orthorhombic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            LOC_ORIENT(1,J) = SJ*RL(J)
            LOC_ORIENT(2,K) = SK*RL(K)
            LOC_ORIENT(3,I) = SI*RL(I)
         ELSE IF (XTAL_SYSTEM.EQ.4) THEN
C
C====== Tetragonal
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(1)
            RL(3) = LOC_CELL_PARAMS(2)
 
            LOC_ORIENT(1,J) = SJ*RL(J)
            LOC_ORIENT(2,K) = SK*RL(K)
            LOC_ORIENT(3,I) = SI*RL(I)
         ELSE IF (XTAL_SYSTEM.EQ.5) THEN
C
C====== Hexagonal
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(1)
            RL(3) = LOC_CELL_PARAMS(2)
 
            CC(1) = 0.0
            CC(2) = 0.0
            CC(3) = 0.5
 
            SS(1) = 1.0
            SS(2) = 1.0
            SS(3) = SQRT(3.0)/2.0
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*CC(I)
            LOC_ORIENT(2,K) = SK*RL(K)*SS(I)*SS(J)
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         ELSE IF (XTAL_SYSTEM.EQ.6) THEN
C
C====== Rhombohedral
C
            A_AXIS = LOC_CELL_PARAMS(1)
            COS_ALPHA = COS(LOC_CELL_PARAMS(2))
            SIN_ALPHA = SIN(LOC_CELL_PARAMS(2))
            FUN = SQRT(1.-3.*COS_ALPHA**2+2.*COS_ALPHA**3)
 
            LOC_ORIENT(1,J) = SJ*A_AXIS*SIN_ALPHA
            LOC_ORIENT(1,K) = SJ*A_AXIS*COS_ALPHA*(1.-COS_ALPHA)
     +                        /SIN_ALPHA
            LOC_ORIENT(2,K) = SK*A_AXIS*FUN/SIN_ALPHA
            LOC_ORIENT(3,I) = SI*A_AXIS
            LOC_ORIENT(3,J) = SI*A_AXIS*COS_ALPHA
            LOC_ORIENT(3,K) = SI*A_AXIS*COS_ALPHA
         ELSE IF (XTAL_SYSTEM.EQ.7) THEN
C
C====== Cubic
C
            A_AXIS = LOC_CELL_PARAMS(1)
 
            LOC_ORIENT(1,J) = SJ*A_AXIS
            LOC_ORIENT(2,K) = SK*A_AXIS
            LOC_ORIENT(3,I) = SI*A_AXIS
         ELSE
C
C====== Triclinic
C
            RL(1) = LOC_CELL_PARAMS(1)
            RL(2) = LOC_CELL_PARAMS(2)
            RL(3) = LOC_CELL_PARAMS(3)
 
            CC(1) = COS(LOC_CELL_PARAMS(4))
            CC(2) = COS(LOC_CELL_PARAMS(5))
            CC(3) = COS(LOC_CELL_PARAMS(6))
 
            SS(1) = SIN(LOC_CELL_PARAMS(4))
            SS(2) = SIN(LOC_CELL_PARAMS(5))
            SS(3) = SIN(LOC_CELL_PARAMS(6))
C
C====== See Buerger's Xray Crystallography for the next formulae, p348
C
            COS_PSI = (CC(I)-CC(K)*CC(J))/SS(K)
            COS_RHO = SQRT(1.0-CC(I)**2-CC(J)**2-CC(K)**2+2.0*CC(I)
     +                *CC(J)*CC(K))/SS(K)
 
            LOC_ORIENT(1,J) = SJ*RL(J)*SS(K)
            LOC_ORIENT(1,K) = SJ*RL(K)*COS_PSI
            LOC_ORIENT(2,K) = SK*RL(K)*COS_RHO
            LOC_ORIENT(3,I) = SI*RL(I)
            LOC_ORIENT(3,J) = SI*RL(J)*CC(K)
            LOC_ORIENT(3,K) = SI*RL(K)*CC(J)
         END IF
C
C====== End of code from MAKE_ORIENT--------------------
C
      END IF
C
C====== In-line version of A_CALC - made DOUBLE PRECISION throughout.---------
C
C
C====== Define phi rotation matrix
C
      CPX = COS(LOC_PHIS(1)*DEG_TO_RAD)
      SPX = SIN(LOC_PHIS(1)*DEG_TO_RAD)
      CPY = COS(LOC_PHIS(2)*DEG_TO_RAD)
      SPY = SIN(LOC_PHIS(2)*DEG_TO_RAD)
      CPZ = COS(LOC_PHIS(3)*DEG_TO_RAD)
      SPZ = SIN(LOC_PHIS(3)*DEG_TO_RAD)
 
      DO 80 I = 1 , 3
         DO 70 J = 1 , 3
            CX(I,J) = 0.
 70      CONTINUE
 80   CONTINUE
      CX(1,1) = 1.0
      CX(2,2) = CPX
      CX(2,3) = -SPX
      CX(3,2) = SPX
      CX(3,3) = CPX
 
      DO 100 I = 1 , 3
         DO 90 J = 1 , 3
            CY(I,J) = 0.
 90      CONTINUE
 100  CONTINUE
      CY(1,1) = CPY
      CY(1,3) = SPY
      CY(2,2) = 1.0
      CY(3,1) = -SPY
      CY(3,3) = CPY
 
      DO 120 I = 1 , 3
         DO 110 J = 1 , 3
            CZ(I,J) = 0.
 110     CONTINUE
 120  CONTINUE
      CZ(1,1) = CPZ
      CZ(1,2) = -SPZ
      CZ(2,1) = SPZ
      CZ(2,2) = CPZ
      CZ(3,3) = 1.0
 
      DO 150 I = 1 , 3
         DO 140 J = 1 , 3
            C(I,J) = 0.0
            DO 130 K = 1 , 3
               C(I,J) = C(I,J) + CY(I,K)*CX(K,J)
 130        CONTINUE
 140     CONTINUE
 150  CONTINUE
 
      DO 180 I = 1 , 3
         DO 170 J = 1 , 3
            PC(I,J) = 0.0
            DO 160 K = 1 , 3
               PC(I,J) = PC(I,J) + CZ(I,K)*C(K,J)
 160        CONTINUE
 170     CONTINUE
 180  CONTINUE
C
C====== Apply phi rotations
C
      DO 210 I = 1 , 3
         DO 200 J = 1 , 3
            A(I,J) = 0.0
            DO 190 K = 1 , 3
               A(I,J) = A(I,J) + PC(I,K)*LOC_ORIENT(K,J)
 190        CONTINUE
 200     CONTINUE
 210  CONTINUE
C
C====== End of text included from a_calc --------------------
C
C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
 
      DO 220 I = 1 , M
         J = REF_NODALS(I)
         H = NODAL_H_S(J)
         K = NODAL_K_S(J)
         L = NODAL_L_S(J)
 
         X1 = A(1,1)*H + A(1,2)*K + A(1,3)*L
         Y1 = A(2,1)*H + A(2,2)*K + A(2,3)*L
         Z1 = A(3,1)*H + A(3,2)*K + A(3,3)*L
C
C====== Apply spindle rotation
C
         X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
         Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
         Z = Z1
 
         RF = (X*X+Y*Y+Z*Z)/(2*X)
 
         XF = LOC_CTOF*Y/(X-RF)
         YF = LOC_CTOF*Z/(X-RF)
C
C====== K_factor = 1.0 + tilt*xf + twist*yf + bulge*sqrt(xf*xf + yf*yf)
C
         K_FACTOR = 1.0
 
         XM = X_CEN + K_FACTOR*(XF*COS_OMEGA-YF*SIN_OMEGA)*MM_RAST_X
         YM = Y_CEN + K_FACTOR*(XF*SIN_OMEGA+YF*COS_OMEGA)*Y_SCRAST
 
         FVECC(I) = SQRT(((XM_OBS(J)-XM)/MM_RAST_X)**2+((YM_OBS(J)-YM)/
     +              MM_RAST_Y)**2)
 
 220  CONTINUE
 
      RETURN
      END
**==displaue.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       DISPLAUE      **
C*************************
C
C
      SUBROUTINE DISPLAUE(FILM_NAME,FILM_UNIT,NUM_CURS,NO_FILM)
C
      INCLUDE 'implicit_none.finc'
C
C
C Purpose : program to display a film image
C           7/6/84 developed for mmh small spots
C           20/2/85  adapted to display laue film and calc nodal spots
C
C           calculate a radial background curve
C           and use a threshold above this for display
C
C           generate file read for nodal spots
C           which are then output to the screen
C           the cursor can then be used to match nodal spots
C           with spots on the film
C           for output for subsequent GENLAUE refinement
C
C           modified to include fiducial constants 29/9/85
C
C
C Author : PAM
C
      CHARACTER * (*)  FILM_NAME
      INTEGER FILM_UNIT,NUM_CURS
      LOGICAL NO_FILM
C
      INCLUDE 'backgr.finc'
      INCLUDE 'sizes.finc'
 
      COMMON /PLOTS/SCAN_START_X,SCAN_START_Y,IDX,IDY
      INTEGER IDX,IDY,SCAN_START_X,SCAN_START_Y
      COMMON /CMBACK/IBGY,IBGV
      INTEGER IBGY(30),IBGV(30)
      COMMON /NODAL/NTOT,NPLOTTED,NODE
      INTEGER NTOT,NPLOTTED,NODE(5,1000)
C
      INCLUDE 'fiducials.finc'
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      COMMON /OMEGA/COSOM,SINOM
      REAL COSOM,SINOM
      INTEGER LUNSTO
C
C====== LOCALS:
C
      INTEGER IBB(MAX_RASTERS)
      CHARACTER *1 ANS
      CHARACTER*40 FILNAM
      CHARACTER *40 IMG_FILE
      LOGICAL QUIT,IMSAV,IMREP
      LOGICAL ERR
      LOGICAL DRAW_VEC
      INTEGER I,J,JJ,K,KK,NUM_PLOT_X,NUM_PLOT_Y,NREC,NYREC
      INTEGER NEXT_BACK_SEG,COFG
      INTEGER JTHRESH,ITHRESH,ITYPE,IFAIL
      REAL OMEGA,OMRAD,X,Y
      INTEGER IRECL,IER,IBLK,CUR_BACK_SEG,VEC_BEGIN
      INTEGER IADDR,IADR,IVAL,VEC_END
C
C-----------------------------------------------------------------------
C
      NUM_PLOT_X = 700
      NUM_PLOT_Y = 700
      IF (NO_FILM) THEN
         CALL READ_SPOTS(NUM_CURS)
         RETURN
 
      END IF
C
C====== Initialise scanner variables
C
      DO 20 I = 1 , MAX_RASTERS
         IBB(I) = 0
 20   CONTINUE
C
C====== Initialise plot variables
C
      IDX = X_RASTERS/NUM_PLOT_X
      IDY = Y_RASTERS/NUM_PLOT_Y
      IF (IDX.NE.IDY) THEN
         NUM_PLOT_X = MIN(X_RASTERS,Y_RASTERS) / MAX(IDX,IDY)
         NUM_PLOT_Y = NUM_PLOT_X
         IDX = MAX(IDX,IDY)
         IDY = IDX
      END IF
      NREC = NUM_PLOT_X*IDX
      NYREC = NUM_PLOT_Y*IDY
      JTHRESH = 50
      IF (IPLATE) JTHRESH = 500
      ITYPE = 1
 
      OMEGA = W_F + W_C
      OMRAD = OMEGA*3.141592653/180.
      COSOM = COS(OMRAD)
      SINOM = SIN(OMRAD)
 
      SCAN_START_X = X_CEN_F - NREC/2 + 1
      SCAN_START_Y = Y_CEN_F - NYREC/2 + 1
C
C====== Enter spot plot mode
C
 30   CALL TASKI('Type of spots to plot (1=nodal, 2=all)',ITYPE,QUIT)
      IF (ITYPE.LT.1 .OR. ITYPE.GT.2) GO TO 30
C
 40   CALL TWRITE('Do you want to select spots by')
      CALL TWRITE('1= Cursor coordinates')
      CALL TWRITE('or')
      CALL TWRITE('2= Centre of gravity of box centred at cursor')
      COFG = 2
      CALL TASKI(' ',COFG,QUIT)
      IF (COFG.LT.1 .OR. COFG.GT.2) GO TO 40
C
      IMREP = .TRUE.
      CALL TASKQ('Do you want to replot an image',IMREP,QUIT)
      IF (IMREP) THEN
         CALL TPROMP('Image file name (default ext=.img): ')
         CALL TREAD(FILNAM,QUIT)
         CALL FILEXT(FILNAM,'img',IMG_FILE,I)
         CALL QGENTR
         CALL QGCLR
         CALL QGPIXC
         CALL TWRITE('1. nodal spot 2. assoc spot on film')
         CALL TWRITE('Type e to finish, a to abort')
         CALL REPLOT(60,IMG_FILE,LUNSTO(0))
         GO TO 100
      END IF
      CALL TASKI('Enter above background threshold',JTHRESH,QUIT)
C
C====== Output message re cursor input
C
      CALL TLOG(.FALSE.)
      CALL TWRITE('Note when plot of film and nodal spots is done,')
      CALL TWRITE('use the cursor to input pairs of')
      CALL TWRITE('1. nodal spot    and   2. associated spot on film')
      CALL TWRITE('Type e to finish, a to abort')
      CALL TWRITE('! You need > 20 spots for refinement !')
      CALL TLOG(.TRUE.)
      CALL QSBELL
      IMSAV = .FALSE.
      CALL TASKQ('Do you want to save film image',IMSAV,QUIT)
      IF (IMSAV) THEN
         CALL TPROMP('Image file name (default ext=.img): ')
         CALL TREAD(FILNAM,QUIT)
         CALL FILEXT(FILNAM,'img',IMG_FILE,I)
         CALL STASAV(61,IMG_FILE)
      END IF
      CALL TLOG(.FALSE.)
      CALL TPROMP('Type return to continue and WAIT for plot...')
      CALL TREAD(ANS,QUIT)
      CALL TLOG(.TRUE.)
C
C===== Reopen film file film_unit
C
      IF (IPLATE) THEN
         IRECL = 2*Y_RASTERS
         CALL IMPOPN(-FILM_UNIT,FILM_NAME,IRECL,SCAN_START_X,OUT_TERM,
     +               ERR)
      ELSE
         IRECL = Y_RASTERS
         CALL FLMOPN(-FILM_UNIT,FILM_NAME,IRECL,SCAN_START_X,OUT_TERM,
     +               ERR)
      END IF
C
C====== Initialise plotter package
C
      CALL QGENTR
      CALL QGCLR
      CALL QGPIXC
      CALL TWRITE('1. nodal spot 2. assoc spot on film')
      CALL TWRITE('Type e to finish, a to abort')
C
C====== Do the display
C
      IBLK = SCAN_START_X
      DO 90 J = 1 , NUM_PLOT_X
C
C====== Set background divisions in y for this scan
C
         CALL BKSET(IBLK)
 
         DO 80 JJ = 1 , IDX
            ITHRESH = IBGV(1) + JTHRESH
            NEXT_BACK_SEG = 2
            CUR_BACK_SEG = IBGY(1)
            IF (IPLATE) THEN
               CALL IMPLRD(FILM_UNIT,IBB,Y_RASTERS,6,ERR)
            ELSE
               CALL FILMRD(FILM_UNIT,IBB,Y_RASTERS,6,ERR)
            END IF
            IBLK = IBLK + 1
C
C====== Loop over record and plot
C
            DRAW_VEC = .FALSE.
            IADDR = SCAN_START_Y
C
C====== Find start of background
C
 50         IF (IADDR.GE.CUR_BACK_SEG) THEN
               CUR_BACK_SEG = IBGY(NEXT_BACK_SEG)
               ITHRESH = IBGV(NEXT_BACK_SEG) + JTHRESH
               NEXT_BACK_SEG = NEXT_BACK_SEG + 1
               GO TO 50
 
            END IF
 
            DO 70 K = 1 , NUM_PLOT_Y
               IVAL = 0
C
C====== Set ival to max. of the (idy) scan pts at this plot pt.
C
               DO 60 KK = 1 , IDY
                  IF (IBB(IADDR).GT.IVAL) IVAL = IBB(IADDR)
                  IADDR = IADDR + 1
 60            CONTINUE
 
               IF (IADDR.GE.CUR_BACK_SEG) THEN
                  ITHRESH = IBGV(NEXT_BACK_SEG) + JTHRESH
                  CUR_BACK_SEG = IBGY(NEXT_BACK_SEG)
                  NEXT_BACK_SEG = NEXT_BACK_SEG + 1
               END IF
 
               IF (IVAL.GT.ITHRESH .AND. .NOT.DRAW_VEC) THEN
                  VEC_BEGIN = K
                  DRAW_VEC = .TRUE.
                  CALL QGMOVE(FLOAT(J),FLOAT(VEC_BEGIN))
                  CALL PUTSAV
               END IF
 
               IF (IVAL.LE.ITHRESH .AND. DRAW_VEC) THEN
                  VEC_END = (K-1)
                  DRAW_VEC = .FALSE.
                  CALL QGDRAW(FLOAT(J),FLOAT(VEC_END))
                  CALL PUTSAV
               END IF
 
 70         CONTINUE
 
            IF (DRAW_VEC) THEN
               CALL QGDRAW(FLOAT(J),FLOAT(NUM_PLOT_Y))
               CALL PUTSAV
            END IF
 80      CONTINUE
 
 90   CONTINUE
C
C====== End of film plot
C
      IF (IMSAV) CALL FINSAV
      IF (IPLATE) THEN
         CALL IMPCLS(FILM_UNIT)
      ELSE
         CALL FLMCLS(FILM_UNIT)
      END IF
C
C====== Find nodal spots and plot on same scale
C
 100  CALL PLOTGE(ITYPE)
C
C====== Cursor in nodal spots
C
      IF (IPLATE) THEN
         IRECL = 2*Y_RASTERS
         CALL IMPOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      ELSE
         IRECL = Y_RASTERS
         CALL FLMOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      END IF
      CALL NODEIN(NUM_CURS,COFG)
      IF (IPLATE) THEN
         CALL IMPCLS(FILM_UNIT)
      ELSE
         CALL FLMCLS(FILM_UNIT)
      END IF
      RETURN
      END
**==dump_spots.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**       DUMP_SPOTS       **
C****************************
C
C
      SUBROUTINE DUMP_SPOTS (NUM_CURS)
C
 
C
C Purpose : Write out the spots cursored in by the user
C
      INTEGER NUM_CURS
C
C Author : IJC
C
      INCLUDE 'observed_spots.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'gen_file_header.finc'
      INCLUDE 'sizes.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      INTEGER I, DUMP_FIL
      REAL X, Y
      PARAMETER (DUMP_FIL=15)
      CHARACTER*40 FILNAM
      CHARACTER*40 DUMP_FILE
C
C-----------------------------------------------------------------------
C
      CALL TPROMP('Dump file name (default ext=.spots): ')
      CALL TREAD(FILNAM,QUIT)
      CALL FILEXT(FILNAM,'spots',DUMP_FILE,I)
      OPEN (FILE=DUMP_FILE,STATUS='NEW',UNIT=DUMP_FIL)
      WRITE (DUMP_FIL,6010) CRYSTAL_ID , TITLE
      WRITE (DUMP_FIL,6020) NUM_CURS
      DO 20 I = 1 , NUM_CURS
C
C====== Convert from scanner units to mm.s
C
         X = (XM_OBS(I)-X_CEN_F)/MM_RAST_X
         Y = (YM_OBS(I)-Y_CEN_F)/MM_RAST_Y
         WRITE (DUMP_FIL,6030) NODAL_H_S(I) , NODAL_K_S(I) , 
     +                         NODAL_L_S(I) , X , Y
 20   CONTINUE
      CLOSE (DUMP_FIL)
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,A6/1X,A20)
 6020 FORMAT (' Number=',I4)
 6030 FORMAT (' h k l  ',3I4,T40,'x y  ',2F6.1)
      END
**==factor_index.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**       FACTOR_INDEX        **
C*******************************
C
C
      SUBROUTINE FACTOR_INDEX (H, K, L)
C
C Purpose: To factor out the HCF of 3 indices.
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER   H, K, L
C
C====== LOCALS:
C
      INTEGER   HCF, HCF1, HCF2
C
C------------------------------------------------------------------------------
C
      HCF1 = HCF(H,K)
      HCF2 = HCF(HCF1,L)
      H = H/HCF2
      K = K/HCF2
      L = L/HCF2
      RETURN
      END
**==find_fid.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       FIND_FID      **
C*************************
C
C
      SUBROUTINE  FIND_FID (FILM_UNIT, ROUGH_X, ROUGH_Y, FIDX, FIDY)
C
C Purpose: This routine is supposed to locate a fiducial accurately given
C          its rough position
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FILM_UNIT
      REAL       ROUGH_X, ROUGH_Y
      REAL       FIDX, FIDY
C
C Note: if the routine fails to locate the fiducial within 5mm then
C       fidx and fidy are both -1
C
      INCLUDE 'sizes.finc'
      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       MAX_X_SIZE, MAX_SHIFT, MAX_SHIFT2
      INTEGER    XC, YC, X_SIZE, Y_SIZE
      INTEGER    THRESH, NUM_OVER
      REAL       CGX, CGY
      LOGICAL    TRY_AGAIN, QUIT, LOW_THR_FIRST, FIRST_TRY
C
C------------------------------------------------------------------------------
C
C====== Max size is twice starting size
C====== NB fid box is actually fid_box_size/2
C
      MAX_X_SIZE = FID_BOX_SIZE*MM_RAST_X
      MAX_SHIFT = FID_BOX_SIZE/4.*MM_RAST_X
      MAX_SHIFT2 = MAX_SHIFT*MAX_SHIFT
C
C====== Assume background is 125 initially
C
      THRESH = 125
      TRY_AGAIN = .FALSE.
      LOW_THR_FIRST = .FALSE.
      FIRST_TRY = .TRUE.
 
      XC = ROUGH_X
      YC = ROUGH_Y
 20   X_SIZE = FID_BOX_SIZE/2.*MM_RAST_X
      Y_SIZE = FID_BOX_SIZE/2.*MM_RAST_Y
 30   CALL C_OF_G(FILM_UNIT,XC,YC,X_SIZE,Y_SIZE,THRESH,CGX,CGY,NUM_OVER)
C
      IF (NUM_OVER.GT.0) THEN
C
C====== If we have moved far from the rough place, re-position the box
C
         IF (((CGX-XC)**2+(CGY-YC)**2).GT.MAX_SHIFT2) THEN
 
            XC = CGX
            YC = CGY
            GO TO 20
 
         END IF
C
C====== OK found a spot that hasn't shifted far from box centre
C
         FIDX = CGX
         FIDY = CGY
         RETURN
      END IF
C
C====== If flat film, try lowering thresh first.
C====== Only ask question first time round
C
      IF (FIRST_TRY .AND. .NOT.LOW_THR_FIRST) THEN
         LOW_THR_FIRST = .TRUE.
         CALL TASKQ(
     +           'Couldn`t find fiducial - try lowering threshold first'
     +           ,LOW_THR_FIRST,QUIT)
      END IF
C
C====== Next thing to try - enlarge box slightly
C====== If box gets to twice the starting size then abort
C
 40   IF (.NOT.(FIRST_TRY .OR. TRY_AGAIN)) THEN
         TRY_AGAIN = .TRUE.
         CALL TWRITE('*Couldn''t find fiducial*')
         CALL TASKQ('Try expanding box and then lowering threshold',
     +              TRY_AGAIN,QUIT)
         IF (.NOT.TRY_AGAIN) THEN
            WRITE (STR,6010) ROUGH_X , ROUGH_Y
            CALL TWRITE(STR)
            FIDX = ROUGH_X
            FIDY = ROUGH_Y
            RETURN
         END IF
      END IF
C
C====== Expand box second try if not too big
C
      IF (.NOT.LOW_THR_FIRST .AND. X_SIZE.LE.MAX_X_SIZE) THEN
C
C======  OK. - expAND the box by 20% - and make sure it does get bigger
C======        else infinite loop.
C
         X_SIZE = X_SIZE*1.2 + 1
         Y_SIZE = Y_SIZE*1.2 + 1
         WRITE (STR,6020) X_SIZE*2.0
         CALL TWRITE(STR)
         GO TO 30
      END IF
C
C====== Try lowering threshold (again)
C
      IF (THRESH.GT.75) THEN
         THRESH = THRESH - 25
         WRITE (STR,6030) THRESH
         CALL TWRITE(STR)
         GO TO 30
      ELSE
         WRITE (STR,6040) THRESH
         CALL TWRITE(STR)
         THRESH = THRESH - 25
         CALL TASKI('Enter threshold or 0 to give up',THRESH,QUIT)
         IF (THRESH.GT.0) GO TO 30
         IF (FIRST_TRY) THEN
            FIRST_TRY = .FALSE.
C
C======= Reset threshold for use with expanded box
C======= It's going to be lowered by 25 immediately
C
            THRESH = 150
            GO TO 40
         ELSE
C
C====== Failed to locate fiducial - give up!
C
            FIDX = -1.0
            FIDY = -1.0
            RETURN
         END IF
      END IF
C
C====== Format statements
C
 6010 FORMAT ('Using ESTIMATE fiducial position of',2F8.1)
 6020 FORMAT ('Expanding box to ',I4,' rasters across')
 6030 FORMAT ('Setting threshold to ',I3)
 6040 FORMAT ('STILL no density found above ',I4)
C
      END
**==find_fiducials.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*********************************
C**       FIND_FIDUCIALS        **
C*********************************
C
C
      SUBROUTINE  FIND_FIDUCIALS (FILM_UNIT)
C
C Purpose: Locates fiducials on film and determines fiducial parameters
C          i.e. x_cen_f, y_cen_f and w_f
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    FILM_UNIT
C
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
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
      REAL       X_CENTRE, Y_CENTRE
      REAL       RAD_TO_DEG
      PARAMETER (RAD_TO_DEG = 360./(2.*3.1415926))
C
C------------------------------------------------------------------------------
C
      WRITE (STR,6030)
      CALL TWRITE(STR)
C
      DO 20 I = 1 , 3
C
         CALL FIND_FID(FILM_UNIT,FID_POS_X(I),FID_POS_Y(I),FIDX(I),
     +                 FIDY(I))
         IF (FIDX(I).EQ.-1.0) THEN
            WRITE (STR,6100) I
            CALL TWRITE(STR)
            STOP
         END IF
C
 20   CONTINUE
C
      WRITE (STR,6010) FIDX(1) , FIDY(1)
      CALL TWRITE(STR)
      WRITE (STR,6020) FIDX(2) , FIDY(2)
      CALL TWRITE(STR)
      WRITE (STR,6020) FIDX(3) , FIDY(3)
      CALL TWRITE(STR)
C
      X_CEN_F = (FIDX(1)+FIDX(3))/2.0
      Y_CEN_F = (FIDY(1)+FIDY(3))/2.0
C
      W_F = (FIDY(3)-FIDY(2))/(FIDX(3)-FIDX(2))
      W_F = ATAN(W_F)*RAD_TO_DEG
C
      WRITE (STR,6080) RASTER_SIZE*1000. , X_CEN_F , Y_CEN_F
      CALL TWRITE(STR)
      WRITE (STR,6090) W_F
      CALL TWRITE(STR)
C
C====== See if we can find centre by using same procedure for
C====== finding fiducials
C====== Only do this if x_c and y_c are zero
C
      IF (X_C.EQ.0.0 .AND. Y_C.EQ.0.0) THEN
 
         CALL FIND_FID(FILM_UNIT,X_CEN_F,Y_CEN_F,X_CENTRE,Y_CENTRE)
         IF (X_CENTRE.GT.0) THEN
            WRITE (STR,6040) X_CENTRE , Y_CENTRE
            CALL TWRITE(STR)
            X_C = X_CENTRE - X_CEN_F
            Y_C = Y_CENTRE - Y_CEN_F
            WRITE (STR,6050)
            CALL TWRITE(STR)
            WRITE (STR,6060) RASTER_SIZE*1000.
            CALL TWRITE(STR)
            WRITE (STR,6070) X_C , Y_C
            CALL TWRITE(STR)
         ELSE
            WRITE (STR,6110)
            CALL TWRITE(STR)
         END IF
      END IF
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Fiducials found at:  (',F7.1,',',F7.1,')')
 6020 FORMAT ('                     (',F7.1,',',F7.1,')')
 6030 FORMAT ('Locating fiducials...')
 6040 FORMAT ('Beam centre at:      (',F7.1,',',F7.1,')')
 6050 FORMAT ('Fiducial constants as determined from beam centre:')
 6060 FORMAT ('(units are in ',F5.1,' micron scanner steps)')
 6070 FORMAT ('x_c =',G14.6,',    y_c =',G14.6)
 6080 FORMAT ('Fiducial centre (',F5.1,'u units) ',2G14.6)
 6090 FORMAT ('Omega as determined from fiducials (degrees)',G14.6)
 6100 FORMAT ('!!! Failed to locate fiducial ',I1)
 6110 FORMAT ('!!! Failed to locate beam centre')
      END
**==find_nodals.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       FIND_NODALS       **
C*****************************
C
C
      SUBROUTINE  FIND_NODALS
C
C Purpose: To identify and mark spots as nodals
C
      INCLUDE 'implicit_none.finc'
C
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
 
      INTEGER NBITST
 
C
C====== LOCALS:
C
      CHARACTER*80 STR
      INTEGER    HCF
      INTEGER    MAX_HKL_LIMIT, DIM, MIN_NUM_REFLEX, OFFSET
      PARAMETER (MIN_NUM_REFLEX=50)
      PARAMETER (MAX_HKL_LIMIT=12)
      PARAMETER (DIM=2*MAX_HKL_LIMIT+1)
      PARAMETER (OFFSET=MAX_HKL_LIMIT+1)
      INTEGER    NODALS(DIM, DIM, DIM)
      INTEGER    I, J, NUM_NODALS, NUM_NODAL_REFLEX, SEL_NODAL_HKL
      INTEGER    H, K, L, HCF1, HCF2
      LOGICAL    QUIT
      INTEGER    NODAL_HIST(20)
      INTEGER    NODAL_INDEX, NOD_IDX, NNOD_SPA
C
C------------------------------------------------------------------------------
C
      DO 20 I = 1 , MAX_HKL_LIMIT
         NODAL_HIST(I) = 0
 20   CONTINUE
C
C====== Clear array of pointers
C
      DO 50 K = 1 , DIM
         DO 40 J = 1 , DIM
            DO 30 I = 1 , DIM
               NODALS(I,J,K) = 0
 30         CONTINUE
 40      CONTINUE
 50   CONTINUE
C
C====== Find the nodals
C
      NUM_NODAL_REFLEX = 0
      NNOD_SPA = 0
      DO 60 I = 1 , NUM_REFLEX
 
         H = H_S(I)
         K = K_S(I)
         L = L_S(I)
 
         NOD_IDX = NODAL_INDEX(H,K,L)
         IF (NOD_IDX.LE.MAX_HKL_LIMIT) THEN
C
C====== Check if it was a spatial
C
            IF (NBITST(M_FLAGS(I),6,1).EQ.0) THEN
C
C====== Found a nodal!
C
               CALL FACTOR_INDEX(H,K,L)
               NUM_NODAL_REFLEX = NUM_NODAL_REFLEX + 1
               CALL STBITS(M_FLAGS(I),5,1,1)
               IF (NODALS(H+OFFSET,K+OFFSET,L+OFFSET).EQ.0) THEN
C
C====== Inc the histogram (but count distinct nodals)
C======                              --------
C
                  NODAL_HIST(NOD_IDX) = NODAL_HIST(NOD_IDX) + 1
               END IF
               NODALS(H+OFFSET,K+OFFSET,L+OFFSET) = I
            ELSE
               NNOD_SPA = NNOD_SPA + 1
            END IF
         END IF
 60   CONTINUE
      WRITE (STR,6020) NNOD_SPA
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      WRITE (STR,6040)
      CALL TWRITE(STR)
      WRITE (STR,6050)
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
C====== Calc cumulative totals
C
      DO 70 I = 1 , MAX_HKL_LIMIT - 1
         NODAL_HIST(I+1) = NODAL_HIST(I+1) + NODAL_HIST(I)
 70   CONTINUE
 
      DO 80 I = 1 , MAX_HKL_LIMIT
         WRITE (STR,6030) NODAL_HIST(I) , I
         CALL TWRITE(STR)
 80   CONTINUE
      SEL_NODAL_HKL = MAX_HKL_LIMIT
 90   CALL TASKI('Nodal spot selection index',SEL_NODAL_HKL,QUIT)
      IF (SEL_NODAL_HKL.GT.20) GO TO 90
C
C====== Now assign found nodal spots to the ref_nodals vector
C
      NUM_NODALS = 0
      DO 120 K = OFFSET - SEL_NODAL_HKL , OFFSET + SEL_NODAL_HKL
         DO 110 J = OFFSET - SEL_NODAL_HKL , OFFSET + SEL_NODAL_HKL
            DO 100 I = OFFSET - SEL_NODAL_HKL , OFFSET + SEL_NODAL_HKL
               IF (NODALS(I,J,K).NE.0) THEN
                  NUM_NODALS = NUM_NODALS + 1
                  IF (NUM_NODALS.LE.MAX_OBS) THEN
                     REF_NODALS(NUM_NODALS) = NODALS(I,J,K)
                  END IF
               END IF
 100        CONTINUE
 110     CONTINUE
 120  CONTINUE
 
      NUM_OBS = MIN(NUM_NODALS,MAX_OBS)
 
      WRITE (STR,6010) NUM_NODALS
      CALL TWRITE(STR)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('...',I6,' nodal spots')
 6020 FORMAT (I5,' spatials were candidate nodals')
 6030 FORMAT (I5,' spots selected at nodal index of ',I3)
 6040 FORMAT ('Histogram of nodal spots')
 6050 FORMAT ('========================')
      END
**==find_overlaps.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**       FIND_OVERLAPS       **
C*******************************
C
C
      SUBROUTINE  FIND_OVERLAPS (DELTA, EPSILON, XSTART, XEND)
C
C Purpose: Sort reflexions to find and mark energy and spatial overlaps
C
      INCLUDE 'implicit_none.finc'
C
      REAL       DELTA, EPSILON
      REAL       XSTART, XEND
C
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
 
      INTEGER NBITST
C
C====== LOCALS:
C
      CHARACTER*80 STR
      INTEGER    I, J, IS, JS
      LOGICAL    QUIT
      REAL       OV_DELTA, OV_DELTA2
      REAL       XMIN, XMAX, X_STEP, XF
      REAL       X1, Y1, X2, Y2
      REAL       DELTA2, EPSILON2, D2, X_DIFF, Y_DIFF
      INTEGER    CURRENT_REFLEX, REFLEX1, REFLEX2
      INTEGER    MAX_SPOTS_IN_STRIP
      PARAMETER (MAX_SPOTS_IN_STRIP=10000)
      INTEGER    N, NUM_STRIPS, STRIP, NUM_SPAT_OV_REFLEX,
     + NUM_CLOS_OV_REFLEX
      INTEGER    STRIP_Y(MAX_SPOTS_IN_STRIP)
      INTEGER    STRIP_REFLEX(MAX_SPOTS_IN_STRIP)
      INTEGER    STRIP_REFLEX1(MAX_SPOTS_IN_STRIP)
      INTEGER    STRIP_INDEX(MAX_SPOTS_IN_STRIP)
C
C------------------------------------------------------------------------------
C
      WRITE (STR,6010)
      CALL TWRITE(STR)
C
C====== Adjust delta to a minimum of 10 micron
C
      DELTA = MAX(DELTA,0.01)
      DELTA2 = DELTA**2
      EPSILON2 = EPSILON**2
C
      OV_DELTA = 0.01
C------>  CALL  TASKI('Lambda overlap limit',OV_DELTA,QUIT)
      OV_DELTA2 = OV_DELTA*OV_DELTA
C
      NUM_SPAT_OV_REFLEX = 0
      NUM_CLOS_OV_REFLEX = 0
C
C====== First, set all spots to be singlets (and clear flags)
C
      DO 20 I = 1 , NUM_REFLEX
         M_FLAGS(I) = 1
 20   CONTINUE
C
C====== To allow for fact that some strips will contain more spots
C====== than others, fill strips half full on average.
C====== Besides this, it seems to be necessary to take a few extra strips.
C====== (Currently 2 extra)
C====== Changed again for Janos  28 Feb 87 (/3 rather than /2)
C
      NUM_STRIPS = NUM_REFLEX/(MAX_SPOTS_IN_STRIP/3.0) + 1 + 1
C
C====== Loop through the strips
C
      X_STEP = (XEND-XSTART)/NUM_STRIPS
      XMIN = XSTART - DELTA - X_STEP
      XMAX = XMIN + X_STEP + 2.0*DELTA
C
      DO 70 STRIP = 1 , NUM_STRIPS
         N = 0
         XMIN = XMIN + X_STEP
         XMAX = XMAX + X_STEP
         CURRENT_REFLEX = 0
C
C====== Read through list of sorted reflexions
C
 30      CURRENT_REFLEX = CURRENT_REFLEX + 1
         IF (CURRENT_REFLEX.GT.NUM_REFLEX) THEN
C
C======  O.K. now we have finished reading...
C======  ...so sort the spots within the strip
C
            CALL INDXIV(STRIP_Y,N,STRIP_INDEX)
            DO 40 I = 1 , N
               STRIP_REFLEX(I) = STRIP_REFLEX1(STRIP_INDEX(I))
 40         CONTINUE
C
C====== Loop over each spot within the strip
C
            IF (N.GT.1) THEN
               DO 60 I = 1 , N - 1
                  DO 50 J = I + 1 , N
                     REFLEX1 = STRIP_REFLEX(I)
                     REFLEX2 = STRIP_REFLEX(J)
                     X1 = XF_S(REFLEX1)
                     Y1 = YF_S(REFLEX1)
                     X2 = XF_S(REFLEX2)
                     Y2 = YF_S(REFLEX2)
C
C====== If y difference exceeds delta then finish spot i.
C
                     Y_DIFF = ABS(Y1-Y2)
                     IF (Y_DIFF.GT.DELTA) GO TO 60
C
C====== If y difference is possible then check x difference.
C
                     X_DIFF = ABS(X1-X2)
                     IF (X_DIFF.LE.DELTA) THEN
C
C====== If both are possible then check precise difference.
C
                        D2 = X_DIFF**2 + Y_DIFF**2
                        IF (D2.LE.DELTA2) THEN
C
C====== An overlap of some sort ... set both non-measurable and count.
C
                           IF (D2.LE.OV_DELTA2) THEN
C
C====== Wavelength overlap of unknown multiplicity, indicate by zero
C====== in multiplicity field,
C====== preserve spatial overlap information
C
                              CALL STBITS(M_FLAGS(REFLEX1),0,5,0)
                              CALL STBITS(M_FLAGS(REFLEX2),0,5,0)
C
                           ELSE IF (D2.LE.EPSILON2) THEN
C
C====== Some spot with spatial overlap, but too close to separate.
C====== Preserve wavelength overlap information
C====== But do not count twice!
C
                              IF (NBITST(M_FLAGS(REFLEX1),7,1).EQ.0)
     +                            THEN
                                 CALL STBITS(M_FLAGS(REFLEX1),6,2,3)
                                 NUM_CLOS_OV_REFLEX = NUM_CLOS_OV_REFLEX
     +                              + 1
                              END IF
C
                              IF (NBITST(M_FLAGS(REFLEX2),7,1).EQ.0)
     +                            THEN
                                 CALL STBITS(M_FLAGS(REFLEX2),6,2,3)
                                 NUM_CLOS_OV_REFLEX = NUM_CLOS_OV_REFLEX
     +                              + 1
                              END IF
C
                           ELSE
C
C====== Some spot with spatial overlap, suitable for separation.
C
                              IF (NBITST(M_FLAGS(REFLEX1),6,1).EQ.0)
     +                            THEN
                                 CALL STBITS(M_FLAGS(REFLEX1),6,1,1)
                                 NUM_SPAT_OV_REFLEX = NUM_SPAT_OV_REFLEX
     +                              + 1
                              END IF
C
                              IF (NBITST(M_FLAGS(REFLEX2),6,1).EQ.0)
     +                            THEN
                                 CALL STBITS(M_FLAGS(REFLEX2),6,1,1)
                                 NUM_SPAT_OV_REFLEX = NUM_SPAT_OV_REFLEX
     +                              + 1
                              END IF
C
                           END IF
                        END IF
                     END IF
C
 50               CONTINUE
C
C====== End of loop through the strip of spots
C
 60            CONTINUE
            END IF
         ELSE
            XF = XF_S(CURRENT_REFLEX)
            IF (XF.GE.XMIN .AND. XF.LE.XMAX) THEN
C
C====== Add spot to strip, whilst checking for too many spots in strip
C
               N = N + 1
               IF (N.GT.MAX_SPOTS_IN_STRIP) THEN
                  WRITE (STR,6040) MAX_SPOTS_IN_STRIP , STRIP , 
     +                             NUM_STRIPS
                  CALL TWRITE(STR)
                  WRITE (STR,6050)
                  CALL TWRITE(STR)
                  STOP
               END IF
C
               STRIP_Y(N) = YF_S(CURRENT_REFLEX)*1E2
C
C====== Loop for more reflexions
C
               STRIP_REFLEX1(N) = CURRENT_REFLEX
            END IF
            GO TO 30
         END IF
C
C====== End of loop over the strips
C
 70   CONTINUE
C
      WRITE (STR,6020) NUM_SPAT_OV_REFLEX , DELTA
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      WRITE (STR,6030) NUM_CLOS_OV_REFLEX , EPSILON
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
C====== Now see which are nodals
C
      CALL FIND_NODALS
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('===Finding nodals and spatial overlaps...')
 6020 FORMAT ('...',I6,' reflexions spatially overlapped for delta =',
     +        F5.2,' mm.')
 6030 FORMAT (' &',I6,' reflexions too close to measure for epsilon =',
     +        F5.2,' mm.')
 6040 FORMAT ('**OVERFLOW** More than',I6,' spots in strip',I4/' (of ',
     +        I3,' strips)')
 6050 FORMAT ('!!! Error in program module: FIND_OVERLAPS')
      END
**==generate_reflexions.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************************
C**       GENERATE_REFLEXIONS       **
C*************************************
C
C
      SUBROUTINE  GENERATE_REFLEXIONS (WAVSTD, WAVMIN, WAVMAX,
     +            DMIN, FILM_RAD, XSTART, XEND)
C
C Purpose: Generate reflexions, producing hkl, lambda and xf,yf
C
      INCLUDE 'implicit_none.finc'
C
      REAL       WAVSTD, WAVMIN, WAVMAX
      REAL       DMIN, FILM_RAD
      REAL       XSTART, XEND
C
      INCLUDE 'crystal_params.finc'
      INCLUDE 'cell_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'flmconsts.finc'
      INCLUDE 'scnlim.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
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
      INTEGER    H, K, L
      INTEGER    NEXT_H, NEXT_K, NEXT_L
      REAL       DSTAR2, DSTAR_MAX2, DSTAR2_K_MAX, DSTAR2_L_MAX
      LOGICAL    FOUND_ANY_L, FOUND_ANY_K, PLOT, ABSENT
      REAL       X_CEN, Y_CEN
      REAL       A(3,3), SPIN(3,3)
      REAL       X, Y, Z, LAMBDA
      REAL       XF, YF
      REAL       RAD2, FILM_RAD2
      REAL       X1, Y1, Z1
      REAL       RMIN, RMAX, RMIN2, RMAX2
      REAL       RF, FACTOR
      INTEGER    TRIES
      REAL       DEG_TO_RAD
      REAL       XTEM,YTEM,XP_SC,YP_SC,P_SC
C
C------------------------------------------------------------------------------
C
      DEG_TO_RAD = 3.1415926536/180.0
C
      PLOT = .FALSE.
      CALL TASKQ('Do you want calculated pattern displayed on screen',
     +           PLOT,QUIT)
      WRITE (STR,6030)
      CALL TWRITE(STR)
C
C====== Initialise graphics if required
C
      IF (PLOT) THEN
         X_CEN = X_RASTERS/2.0
         Y_CEN = Y_RASTERS/2.0
         XP_SC = X_CEN/MM_RAST_X
         YP_SC = Y_CEN/MM_RAST_Y
         P_SC = AMAX1(XP_SC,YP_SC) + 3.0
         CALL PREPARE_PLOT(-P_SC,P_SC,-P_SC,P_SC)
C
C====== Plot fiducials
C
         IF (NUM_FID.GT.0) THEN
           IF(.NOT.IPLATE) THEN
            CALL QGSYMB((FID_POS_X(1)-X_CEN)/MM_RAST_X,
     +                  (FID_POS_Y(1)-Y_CEN)/MM_RAST_Y,1,3.0)
            CALL QGSYMB((FID_POS_X(2)-X_CEN)/MM_RAST_X,
     +                  (FID_POS_Y(2)-Y_CEN)/MM_RAST_Y,1,3.0)
            CALL QGSYMB((FID_POS_X(3)-X_CEN)/MM_RAST_X,
     +                  (FID_POS_Y(3)-Y_CEN)/MM_RAST_Y,1,3.0)
           END IF
            CALL QGSYMB(0.,0.,1,3.0)
         END IF
      END IF
C
C====== Set up spindle rotation matrix
C
      SPIN(1,1) = COS(SPINDLE*DEG_TO_RAD)
      SPIN(1,2) = -SIN(SPINDLE*DEG_TO_RAD)
      SPIN(2,1) = -SPIN(1,2)
      SPIN(2,2) = SPIN(1,1)
C
      RMIN = WAVSTD/WAVMAX
      RMAX = WAVSTD/WAVMIN
      RMIN2 = RMIN*RMIN
      RMAX2 = RMAX*RMAX
C
      FILM_RAD2 = FILM_RAD*FILM_RAD
C
      DSTAR_MAX2 = (WAVSTD/DMIN)**2
C
      CALL A_CALC(ORIENTATION,PHIS,A)
C
      NUM_REFLEX = 0
      TRIES = 0
      XSTART = 1E6
      XEND = -1E6
C
C====== Start h loop
C
      H = 0
      NEXT_H = 1
C
C====== Start k loop
C
 20   K = 0
      NEXT_K = 1
      FOUND_ANY_K = .FALSE.
      DSTAR2_K_MAX = 1000.0
C
C====== Start l loop - skip 000
C
C
C====== Check h & k for systematic absences
C
 30   IF (ISYS.GT.0) THEN
         CALL SYSAHK(H,K,ABSENT)
         IF (ABSENT) THEN
            K = K + NEXT_K
            GO TO 30
         END IF
      END IF
C
      L = 0
      NEXT_L = 1
      IF (H.EQ.0 .AND. K.EQ.0) L = 1
      FOUND_ANY_L = .FALSE.
      DSTAR2_L_MAX = 1000.0
C
C====== Next l
C
 40   TRIES = TRIES + 1
C
C====== Check h, k & l for systematic absences
C
      IF (ISYS.GT.0) THEN
         CALL SYSABS(H,K,L,ABSENT)
         IF (ABSENT) THEN
            L = L + NEXT_L
            GO TO 40
         END IF
      END IF
C
C====== Check if within resolution
C====== Calculate reciprocal lattice co-ords from hkl values
C
      X1 = A(1,1)*H + A(1,2)*K + A(1,3)*L
      Y1 = A(2,1)*H + A(2,2)*K + A(2,3)*L
      Z1 = A(3,1)*H + A(3,2)*K + A(3,3)*L
C
      DSTAR2 = X1*X1 + Y1*Y1 + Z1*Z1
C
      IF (DSTAR2.LE.DSTAR_MAX2) THEN
C
C====== Apply spindle rotation
C
         X = SPIN(1,1)*X1 + SPIN(1,2)*Y1
         Y = SPIN(2,1)*X1 + SPIN(2,2)*Y1
         Z = Z1
C
C====== Convert x,y,z to xf,yf
C
C====== Is r.l. point between min and max spheres ?
C
         IF ((X+RMIN)**2+Y*Y+Z*Z.GE.RMIN2) THEN
            IF ((X+RMAX)**2+Y*Y+Z*Z.LE.RMAX2) THEN
C
C====== x must be negative, but greater than rf
C
               RF = DSTAR2/(2*X)
               IF (X.LT.0.00001 .AND. X.GT.RF) THEN
C
C====== Project onto film
C
                  XF = CTOF*Y/(X-RF)
                  YF = CTOF*Z/(X-RF)
C
C====== Check that point is inside the circle requested by the user
C
	TEMCX = X_CEN_F
	TEMCY = Y_CEN_F
	IF(.NOT.CENT_IN) THEN
		TEMCX = X_RASTERS/2.0
		TEMCY = Y_RASTERS/2.0
	ENDIF
	XTEM = X_C + TEMCX + XF*MM_RAST_X
	YTEM = Y_C + TEMCY + YF*MM_RAST_Y
C
                  IF (ABS(XF).LE.FILM_RAD) THEN
                     IF (ABS(YF).LE.FILM_RAD) THEN
                        RAD2 = XF*XF + YF*YF
                        IF (RAD2.LE.FILM_RAD2) THEN
	IF(XTEM.GE.XLOW.AND.XTEM.LE.XHIGH.AND.
     +     YTEM.GE.YLOW.AND.YTEM.LE.YHIGH) THEN
C
C====== O.K. so present on film
C
                           NUM_REFLEX = NUM_REFLEX + 1
                           IF (NUM_REFLEX.GT.MAX_REFLEX) THEN
                              NUM_REFLEX = MAX_REFLEX
                              WRITE (STR,6040)
                              CALL TWRITE(STR)
                              GO TO 50
                           END IF
C
                           H_S(NUM_REFLEX) = H
                           K_S(NUM_REFLEX) = K
                           L_S(NUM_REFLEX) = L
                           XF_S(NUM_REFLEX) = XF
                           YF_S(NUM_REFLEX) = YF
                           LAMBDAS(NUM_REFLEX) = -WAVSTD/RF
                           DMIN_THRESHOLDS2(NUM_REFLEX) = DSTAR2
C
                           IF (PLOT) CALL QGSYMB(XF,YF,0,0.15)
C
                           IF (XF.LT.XSTART) XSTART = XF
                           IF (XF.GT.XEND) XEND = XF
	END IF
                        END IF
                     END IF
                  END IF
               END IF
            END IF
         END IF
C
         FOUND_ANY_L = .TRUE.
         FOUND_ANY_K = .TRUE.
C
C====== Increment l and continue
C
         L = L + NEXT_L
         GO TO 40
C
      END IF
C
C====== hkl outside resolution ... test if any hits for this l so far
C
      IF (FOUND_ANY_L) THEN
C
C====== yes ... reverse l direction
C
         IF (NEXT_L.EQ.1) THEN
            NEXT_L = -1
            L = -1
            GO TO 40
C
C====== or increment k if already reversed l
C
         ELSE
            K = K + NEXT_K
            GO TO 30
         END IF
C
      END IF
C
C====== no  ... check if getting closer ... save and continue if so
C
      IF (DSTAR2.LT.DSTAR2_L_MAX) THEN
         DSTAR2_L_MAX = DSTAR2
         L = L + NEXT_L
         GO TO 40
      END IF
C
C====== getting further away ... reverse l and continue
C
      IF (NEXT_L.EQ.1) THEN
         NEXT_L = -1
         L = -1
         GO TO 40
      END IF
C
C====== no hits for this l so far ... test if any for this k so far
C
      IF (FOUND_ANY_K) THEN
C
C====== yes ... reverse k direction
C
         IF (NEXT_K.EQ.1) THEN
            NEXT_K = -1
            K = -1
            GO TO 30
C
C====== or bump h if already reversed k
C
         ELSE
            H = H + NEXT_H
            GO TO 20
         END IF
C
      END IF
C
C====== no  ... check if getting closer ... save and continue if so
C
      IF (DSTAR2.LT.DSTAR2_K_MAX) THEN
         DSTAR2_K_MAX = DSTAR2
         K = K + NEXT_K
         GO TO 30
      END IF
C
C====== getting further away ... reverse k and continue
C
      IF (NEXT_K.EQ.1) THEN
         NEXT_K = -1
         K = -1
         GO TO 30
      END IF
C
C====== check if we can go in other h direction
C
      IF (NEXT_H.EQ.1) THEN
         NEXT_H = -1
         H = -1
         GO TO 20
      END IF
C
C====== no hits at all for this h value so stop searching further
C
 50   IF (PLOT) THEN
         CALL QGBELL
         READ (IN_TERM,'(A)')
         CALL QGALPH
      END IF
      WRITE (STR,6010) NUM_REFLEX , H , K , L
      CALL TWRITE(STR)
      WRITE (STR,6020) TRIES
      CALL TWRITE(STR)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('...',I9,' reflections generated at ',3I4)
 6020 FORMAT ('out of',I9,' hkl values tried.')
 6030 FORMAT ('===Generating reflexions...')
 6040 FORMAT ('!!! Maximum number of reflexions reached - terminating')
      END
**==get_film_info.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C********************************
C**       GET_FILM_INFO        **
C********************************
C
C
      SUBROUTINE  GET_FILM_INFO (FILM_NAME, FILM_NUM)
C
C Purpose: To get various parameters in preparation for matching
C
      INCLUDE 'implicit_none.finc'
C
      CHARACTER*(*)  FILM_NAME
      INTEGER    FILM_NUM
C
      INCLUDE 'fiducials.finc'
      INCLUDE 'backgr.finc'
      INCLUDE 'sizes.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      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 ERR
      CHARACTER*6   FILMS
      CHARACTER*1   FILM_LETTER
      CHARACTER*60 PATH
      CHARACTER*40 NAME
      CHARACTER*12 TYPE
      CHARACTER*6 VERS
      CHARACTER*60 LOGNAM
      CHARACTER*80 FILTMP
      INTEGER    FILM_UNIT
      INTEGER    I
      INTEGER L
      INTEGER IRECL, IFAIL
      INTEGER LENSTR
      LOGICAL    QUIT, EXISTS
      LOGICAL CCPEXS
      DATA   FILMS /'ABCDEF'/
C
C------------------------------------------------------------------------------
C
      FILM_UNIT = 4
C
C====== Get film number
C
      IF (IPLATE) THEN
         FILM_NUM = 1
      ELSE
C
C====== Get film number
C
 20      FILM_NUM = MIN(MAX(FILM_NUM,1),6)
         FILM_LETTER = FILMS(FILM_NUM:FILM_NUM)
         WRITE (STR,6010) FILM_LETTER
         CALL TPROMP(STR)
         CALL TREAD(FILM_LETTER,QUIT)
         IF (FILM_LETTER.EQ.' ') THEN
            FILM_LETTER = FILMS(FILM_NUM:FILM_NUM)
         END IF
         CALL CCPUPC(FILM_LETTER)
         IF (FILM_LETTER.LT.'A' .OR. FILM_LETTER.GT.'F') GO TO 20
         FILM_NUM = ICHAR(FILM_LETTER) - ICHAR('A') + 1
      END IF
C
C====== Get film filename and open file
C
 30   CALL TPROMP('Film file name (default ext=.dat): ')
      CALL TREAD(FILTMP,QUIT)
      CALL FILEXT(FILTMP,'dat',FILM_NAME,I)
C
C====== Assume default device is film: or from env variable FILM
C
      CALL CCPPSF(FILM_NAME,PATH,NAME,TYPE,VERS)
      IF (PATH.EQ.' ') THEN
         CALL UGTENV('FILM',LOGNAM)
         IF (LOGNAM.NE.' ') THEN
            FILTMP = FILM_NAME
            FILM_NAME = LOGNAM
            L = LENSTR(FILM_NAME)
            FILM_NAME(L+1:) = FILTMP
         END IF
      END IF
      EXISTS = CCPEXS(FILM_NAME)
C
      IF (.NOT.EXISTS) THEN
         CALL TWRITE('!!! File not found')
         GO TO 30
      END IF
C
      IF (IPLATE) THEN
         IRECL = 2*Y_RASTERS
         CALL IMPOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      ELSE
         IRECL = Y_RASTERS
         CALL FLMOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      END IF
C
C====== Find fiducials and set fiducial parameters
C
      IF (.NOT.CENT_IN.AND..NOT.IPLATE) CALL FIND_FIDUCIALS(FILM_UNIT)
C
C====== Establish radial background intensity
C
      CALL CALC_BACKGR(FILM_UNIT)
C
      IF (IPLATE) THEN
         CALL IMPCLS(FILM_UNIT)
      ELSE
         CALL FLMCLS(FILM_UNIT)
      END IF
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Enter film to be refined by letter, A to F [',A1,']: ')
      END
**==get_gen_paramters.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C************************************
C**       GET_GEN_PARAMTERS        **
C************************************
C
C
      SUBROUTINE  GET_GEN_PARAMTERS (WAVSTD, WAVMIN, WAVMAX,
     +                   DMIN, FILM_RAD, X_REJECT, Y_REJECT)
C
      INCLUDE 'implicit_none.finc'
C
      REAL    WAVSTD, WAVMIN, WAVMAX
      REAL    DMIN, FILM_RAD
      REAL    X_REJECT, Y_REJECT
C
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'flmconsts.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'gen_file_header.finc'
      INCLUDE 'scnlim.finc'
C
      CHARACTER*1 TYP
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      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 IFAIL
      LOGICAL QUIT
      CHARACTER*120 STR
      CHARACTER*40  FNAME, GEN_FILE
      INTEGER    I, J, NVALS
      LOGICAL    ERROR, EXISTS
      LOGICAL    CCPEXS
      REAL       VALUES(10)
C
C------------------------------------------------------------------------------
C
 20   WRITE (STR,6120)
      CALL TPROMP(STR)
      CALL TREAD(FNAME,QUIT)
C
      IF (FNAME(1:1).NE.' ') THEN
         CALL FILEXT(FNAME,'gen',GEN_FILE,I)
         EXISTS = CCPEXS(GEN_FILE)
         IF (.NOT.EXISTS) THEN
            CALL TWRITE('!!! Cannot find file')
            GO TO 20
         END IF
         CALL TWRITE(' ')
         WRITE (STR,6130)
         CALL TPROMP(STR)
         CALL TREAD(TYP,QUIT)
         IF (TYP.EQ.'i' .OR. TYP.EQ.'I') IPLATE = .TRUE.
         IF (IPLATE) THEN
            CALL TWRITE(' ')
            CALL TWRITE('Processing image plate data')
         END IF
         IFAIL = 0
         CALL CCPDPN(-9,GEN_FILE,'READONLY','F',0,IFAIL)
C
C====== Read in `nearly-constants' from .GEN file or LAUEFLM.DAT
C
         READ (9,6140) TITLE
         WRITE (STR,6150) TITLE
         CALL TWRITE(STR)
         READ (9,6180) CRYSTAL_ID
         WRITE (STR,6190) CRYSTAL_ID
         CALL TWRITE(STR)
         CALL GET_ORIENT(9,ERROR,ORIENTATION,UNIT_CELL_INFO)
         IF (ERROR) GO TO 30
         WRITE (STR,6150) UNIT_CELL_INFO
         CALL TWRITE(STR)
         WRITE (STR,'(3F12.7)') (ORIENTATION(1,I),I=1,3)
         CALL TWRITE(STR)
         WRITE (STR,'(3F12.7)') (ORIENTATION(2,I),I=1,3)
         CALL TWRITE(STR)
         WRITE (STR,'(3F12.7)') (ORIENTATION(3,I),I=1,3)
         CALL TWRITE(STR)
         READ (9,END=40,FMT=*) PHIS
         WRITE (STR,'(3F12.4)') PHIS
         CALL TWRITE(STR)
         READ (9,END=50,FMT=*) WAVSTD , WAVMIN , WAVMAX
         WRITE (STR,'(3F12.4)') WAVSTD , WAVMIN , WAVMAX
         CALL TWRITE(STR)
         READ (9,END=60,FMT=*) CTOF , DMIN , FILM_RAD
         WRITE (STR,'(3F12.4)') CTOF , DMIN , FILM_RAD
         CALL TWRITE(STR)
         READ (9,END=70,FMT=*) SPINDLE
         WRITE (STR,'(F12.4)') SPINDLE
         CALL TWRITE(STR)
         READ (9,END=80,FMT=*) X_C , Y_C , W_C
         WRITE (STR,'(3F12.7)') X_C , Y_C , W_C
         CALL TWRITE(STR)
         READ (9,END=90,FMT=*) TWIST , TILT , BULGE
         WRITE (STR,'(3F12.7)') TWIST , TILT , BULGE
         CALL TWRITE(STR)
         READ (9,END=100,FMT=*) Y_SCALE
         WRITE (STR,'(F12.7)') Y_SCALE
         CALL TWRITE(STR)
         READ (9,END=110,FMT=*) X_REJECT , Y_REJECT
         WRITE (STR,'(2F12.7)') X_REJECT , Y_REJECT
         CALL TWRITE(STR)
         REWIND 9
         CALL READCONSTS(9)
         CLOSE (9)
C
C====== Film constants
C
         WRITE (STR,'(I3,6F9.1)') NUM_FID , 
     +                            (FID_POS_X10U(I),FID_POS_Y10U(I),I=1,
     +                            3)
         CALL TWRITE(STR)
         WRITE (STR,'(F6.1,F6.2,I5,3F7.3)') FILM_WIDTH , FID_BOX_SIZE , 
     +          N1OD , BASEOD , G1OD , CURVE
         CALL TWRITE(STR)
         WRITE (STR,'(2I6,4F10.3)') X_RASTERS , Y_RASTERS , 
     +                             FILMPLOT_SCALE , X_CEN_F , Y_CEN_F , 
     +                             RASTER_SIZE*1000.0
         CALL TWRITE(STR)
         WRITE(STR,'(4F6.0)') XLOW,XHIGH,YLOW,YHIGH
         CALL TWRITE(STR)
         GO TO 120
C
      END IF
C
      WRITE (STR,6160)
      CALL TPROMP(STR)
      CALL TREAD(TITLE,QUIT)
      WRITE (STR,6170)
      CALL TPROMP(STR)
      CALL TREAD(CRYSTAL_ID,QUIT)
C
 30   WRITE (STR,6010)
      CALL TWRITE(STR)
      WRITE (STR,6020)
      CALL TWRITE(STR)
      WRITE (STR,6030)
      CALL TPROMP(STR)
      CALL GET_ORIENT(IN_TERM,ERROR,ORIENTATION,UNIT_CELL_INFO)
      IF (ERROR) GO TO 30
C
 40   WRITE (STR,6040)
      PHIS(1) = 0
      PHIS(2) = 0
      PHIS(3) = 0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) PHIS(1) = VALUES(1)
      IF (NVALS.GE.2) PHIS(2) = VALUES(2)
      IF (NVALS.GE.3) PHIS(3) = VALUES(3)
      READ (IN_TERM,*) PHIS
C
 50   WRITE (STR,6050)
      WAVSTD = 0.0
      WAVMIN = 0.0
      WAVMAX = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) WAVSTD = VALUES(1)
      IF (NVALS.GE.2) WAVMIN = VALUES(2)
      IF (NVALS.GE.3) WAVMAX = VALUES(3)
C
 60   WRITE (STR,6060)
      CTOF = 0.0
      DMIN = 0.0
      FILM_RAD = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) CTOF = VALUES(1)
      IF (NVALS.GE.2) DMIN = VALUES(2)
      IF (NVALS.GE.3) FILM_RAD = VALUES(3)
C
 70   WRITE (STR,6070)
      SPINDLE = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) SPINDLE = VALUES(1)
C
 80   WRITE (STR,6080)
      X_C = 0.0
      Y_C = 0.0
      W_C = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) X_C = VALUES(1)
      IF (NVALS.GE.2) Y_C = VALUES(2)
      IF (NVALS.GE.3) W_C = VALUES(3)
C
 90   WRITE (STR,6090)
      TWIST = 0.0
      TILT = 0.0
      BULGE = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) TWIST = VALUES(1)
      IF (NVALS.GE.2) TILT = VALUES(2)
      IF (NVALS.GE.3) BULGE = VALUES(3)
C
 100  WRITE (STR,6100)
      Y_SCALE = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) Y_SCALE = VALUES(1)
C
 110  WRITE (STR,6110)
      X_REJECT = 0.0
      Y_REJECT = 0.0
      CALL TPROMP(STR)
      CALL TREAD(STR,QUIT)
      CALL STRVAL(STR,' ,',1,VALUES,NVALS)
      IF (NVALS.GE.1) X_REJECT = VALUES(1)
      IF (NVALS.GE.2) Y_REJECT = VALUES(2)
      CALL READCONSTS(0)
C
C====== Convert x_c, y_c to scanner space
C
 120  X_C = X_C*MM_RAST_X
      Y_C = Y_C*MM_RAST_Y
      RETURN
C
 6010 FORMAT ('Enter recip. unit cell parameters (A) & (degrees)',
     +        ' and orientation info')
 6020 FORMAT ('Eg. ".01 62. RHOMBO P +ast -bst"')
 6030 FORMAT ('(If lattice type omitted P assumed): ')
 6040 FORMAT ('Enter phi angles (degrees): ')
 6050 FORMAT ('Enter wavelength standard, min and max (A):')
 6060 FORMAT ('Enter ctof (mm), dmin (A), film radius (mm):')
 6070 FORMAT ('Enter spindle angle (degrees):')
 6080 FORMAT ('Enter x_c (mm), y_c (mm), w_c (degrees):')
 6090 FORMAT ('Enter twist, tilt, bulge (in units of equation):')
 6100 FORMAT ('Enter Y scaling factor:')
 6110 FORMAT ('Enter central rejection box dimensions, x_size, y_size',
     +        ' (mm):')
 6120 FORMAT ('Enter root name for .gen file: ')
 6130 FORMAT ('Processing film (f) or image plate data (i) [f]: ')
 6140 FORMAT (A76)
 6150 FORMAT (1X,A76)
 6160 FORMAT ('Enter title: ')
 6170 FORMAT ('Enter crystal identifier: ')
 6180 FORMAT (A6)
 6190 FORMAT (1X,A6)
      END
**==get_orient.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       GET_ORIENT        **
C*****************************
C
C
      SUBROUTINE  GET_ORIENT (UNIT, ERROR, ORIENTATION, UNIT_CELL_INFO)
C
C Purpose: To derive orientation matrix from cell parameters - which are
C          to be read from the specified unit
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    UNIT
      LOGICAL    ERROR
      REAL       ORIENTATION(3,3)
      CHARACTER*80  UNIT_CELL_INFO
C
      INCLUDE 'cell_params.finc'
C
C====== LOCALS:
C
      INTEGER       NVALS, NWORDS
      CHARACTER*80  STR
      CHARACTER*32  WORDS(10)
      CHARACTER*40  BLANKS40
      INTEGER      EXPECTED_NVALS(7)
      CHARACTER*3  SYSTEM_NAME
      CHARACTER*5  SYSTEMS(7)
      INTEGER      I, AXES(2), ISIGN
      CHARACTER*1  CHAR
      CHARACTER*2  AXIS_NAMES(7)
      CHARACTER*1 LATTICE_TYPES(7)
      REAL   VALUES(10)
      DATA   BLANKS40  /'                                        '/
C
      DATA   SYSTEMS  /'TRICL',
     1                 'MONOC',
     2                 'ORTHO',
     3                 'TETRA',
     4                 'HEXAG',
     5                 'RHOMB',
     6                 'CUBIC'/
C
      DATA   EXPECTED_NVALS  /6, 4, 3, 2, 2, 2, 1/
      DATA   AXIS_NAMES /'-c', '-b', '-a', ' ?', '+a', '+b', '+c'/
      DATA   LATTICE_TYPES /'P', 'A', 'B', 'C', 'I', 'F', 'R'/
C
C------------------------------------------------------------------------------
C
C====== Expect input to be of form:
C
C     <cell_parameters> <crystal_system> <lattice> <spindle_axis> <beam_axis>
C
C where:
C     <cell_parameters> are up to six numbers, and are in the order
C                       astar, bstar, cstar, alpha, beta, gamma
C                       Only supply as many as needed e.g. for rhombohedral
C                       just give   astar, alpha.
C
C     <crystal_system>  name of crystal system - only first three letters
C                       are used, e.g. tri for triclinic
C
C     <lattice>         a single letter for the type of unit cell. Default P
C
C     <beam_axis>       a word is needed here - only first two characters used
C                       i.e.  can start with  +a , +b , +c , -a , -b , -c
C                       specifies alignment of reciprocal lattice wrt.
C                       beam axis (i.e. X-axis)
C
C     <spindle_axis>    similar form to <beam_axis> but this time specifies
C                       alignment wrt. spindle axis (i.e. Z-axis)
C
C An e.g. of a complete input line is:
C
C     0.123 0.456 0.0789  ortho a +astar -cstar
C
C but could be shortened to
C
C     0.123 0.456 0.0789  ort a a -c
C
C Upper or lower case can be used and intermixed freely
C
C------------------------------------------------------------------------------
C
      ERROR = .FALSE.
C
      READ (UNIT,'(A)') STR
      CALL STRTKN(STR,' ,',10,WORDS,NWORDS,VALUES,NVALS)
      IF (NVALS.GT.0) THEN
         DO 20 I = 1 , NVALS
            CELL_PARAMS(I) = VALUES(I)
 20      CONTINUE
      END IF
C
      IF (NWORDS.LE.4) THEN
         IF (NWORDS.EQ.3) THEN
            WORDS(4) = WORDS(3)
            WORDS(3) = WORDS(2)
            WORDS(2) = 'P'
         END IF
C
         CALL CCPUPC(WORDS(1))
         CALL CCPUPC(WORDS(2))
         CALL CCPUPC(WORDS(3))
         CALL CCPUPC(WORDS(4))
C
         SYSTEM_NAME = WORDS(1)(1:3)
         DO 30 I = 1 , 7
            IF (SYSTEM_NAME.EQ.SYSTEMS(I)(1:3)) GO TO 40
 30      CONTINUE
      END IF
      GO TO 60
C
 40   XTAL_SYSTEM = I
      IF (NVALS.EQ.EXPECTED_NVALS(I)) THEN
C
         DO 50 I = 1 , 2
            IF (WORDS(I+2)(1:1).EQ.'-') THEN
               ISIGN = -1
               CHAR = WORDS(I+2)(2:2)
            ELSE IF (WORDS(I+2)(1:1).EQ.'+') THEN
               ISIGN = 1
               CHAR = WORDS(I+2)(2:2)
            ELSE
               CALL TWRITE('Unspecified direction for axis. Assuming +')
               ISIGN = 1
               CHAR = WORDS(I+2)(1:1)
            END IF
            AXES(I) = 0
            IF (CHAR.EQ.'A') AXES(I) = 1
            IF (CHAR.EQ.'B') AXES(I) = 2
            IF (CHAR.EQ.'C') AXES(I) = 3
            AXES(I) = AXES(I)*ISIGN
            IF (AXES(I).EQ.0) GO TO 60
 50      CONTINUE
         Z_AXIS = AXES(1)
         X_AXIS = AXES(2)
C
C======Set up arrays for systematic absences.
C
C====== r only applies to rhombohedral indexed on hexagonal lattice
C
         IF (WORDS(2).EQ.'R' .AND. XTAL_SYSTEM.NE.5) WORDS(2) = 'P'
         IF (WORDS(2).EQ.'P') THEN
            LATTICE = 1
            ISYS = 0
            KSYS(1) = 0
            KSYS(2) = 0
            KSYS(3) = 0
         ELSE IF (WORDS(2).EQ.'A') THEN
            LATTICE = 2
            ISYS = 2
            KSYS(2) = 1
            KSYS(3) = 1
         ELSE IF (WORDS(2).EQ.'B') THEN
            LATTICE = 3
            ISYS = 2
            KSYS(1) = 1
            KSYS(3) = 1
         ELSE IF (WORDS(2).EQ.'C') THEN
            LATTICE = 4
            ISYS = 2
            KSYS(1) = 1
            KSYS(2) = 1
         ELSE IF (WORDS(2).EQ.'I') THEN
            LATTICE = 5
            ISYS = 2
            KSYS(1) = 1
            KSYS(2) = 1
            KSYS(3) = 1
         ELSE IF (WORDS(2).EQ.'F') THEN
            LATTICE = 6
            ISYS = 4
         ELSE IF (WORDS(2).EQ.'R') THEN
            LATTICE = 7
            ISYS = 3
            KSYS(1) = -1
            KSYS(2) = 1
            KSYS(3) = 1
         ELSE
            GO TO 60
         END IF
C
         CALL MAKE_ORIENT(CELL_PARAMS,XTAL_SYSTEM,X_AXIS,Z_AXIS,
     +                    ORIENTATION)
C
         UNIT_CELL_INFO = BLANKS40//BLANKS40
         IF (XTAL_SYSTEM.EQ.1) THEN
            WRITE (UNIT_CELL_INFO(1:60),6010) (CELL_PARAMS(I),I=1,6)
         ELSE IF (XTAL_SYSTEM.EQ.2) THEN
            WRITE (UNIT_CELL_INFO(1:60),6020) (CELL_PARAMS(I),I=1,4)
         ELSE IF (XTAL_SYSTEM.EQ.3) THEN
            WRITE (UNIT_CELL_INFO(1:60),6030) (CELL_PARAMS(I),I=1,3)
         ELSE IF (XTAL_SYSTEM.EQ.4) THEN
            WRITE (UNIT_CELL_INFO(1:60),6040) (CELL_PARAMS(I),I=1,2)
         ELSE IF (XTAL_SYSTEM.EQ.5) THEN
            WRITE (UNIT_CELL_INFO(1:60),6050) (CELL_PARAMS(I),I=1,2)
         ELSE IF (XTAL_SYSTEM.EQ.6) THEN
            WRITE (UNIT_CELL_INFO(1:60),6060) (CELL_PARAMS(I),I=1,2)
         ELSE IF (XTAL_SYSTEM.EQ.7) THEN
            WRITE (UNIT_CELL_INFO(1:60),6070) CELL_PARAMS(1)
         ELSE
            CALL TWRITE('!!! Illegal Crystal System Code')
            STOP
         END IF
         UNIT_CELL_INFO = UNIT_CELL_INFO(1:60)
     +                    //' '//SYSTEMS(XTAL_SYSTEM)
     +                    //' '//LATTICE_TYPES(LATTICE)
     +                    //' '//AXIS_NAMES(AXES(1)+4)
     +                    //'st '//AXIS_NAMES(AXES(2)+4)//'st'
C
         NUM_CELL_PARAMS = NVALS
C
         RETURN
      END IF
C
 60   ERROR = .TRUE.
      RETURN
C
C====== Format statements
C
 6010 FORMAT (3F10.7,3F10.5)
 6020 FORMAT (3F10.7,F10.5)
 6030 FORMAT (3F10.7)
 6040 FORMAT (2F10.7)
 6050 FORMAT (2F10.7)
 6060 FORMAT (F10.7,F10.5)
 6070 FORMAT (F10.7)
      END
**==get_raster_size.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************************
C**       GET_RASTER_SIZE        **
C**********************************
C
C
      SUBROUTINE  GET_RASTER_SIZE(RASTER_SIZE)
C
C Purpose: Match value of raser_size on input against the range of values
C          supported, and then give the user a change to over-ride.
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER NUM_POSS_VALS
      PARAMETER (NUM_POSS_VALS=4)
      REAL RASTER_SIZE
      REAL POSSIBLE_VALS(NUM_POSS_VALS)
      REAL RANGE_LO(NUM_POSS_VALS), RANGE_HI(NUM_POSS_VALS)
      INTEGER I
      LOGICAL QUIT
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      DATA POSSIBLE_VALS/12.5E-3,25.E-3,50.E-3,100.E-3/
      DATA RANGE_LO/10.E-3,20.E-3,40.E-3,90.E-3/,
     + RANGE_HI/15.E-3,30.E-3,60.E-3,110.E-3/
C
C------------------------------------------------------------------------------
C
      IF (.NOT.IPLATE) THEN
         DO 20 I = 1 , NUM_POSS_VALS
            IF (RASTER_SIZE.GE.RANGE_LO(I) .AND. 
     +          RASTER_SIZE.LE.RANGE_HI(I)) THEN
               RASTER_SIZE = POSSIBLE_VALS(I)
               GO TO 30
            END IF
 20      CONTINUE
         CALL TWRITE('Raster size calculated from film dimensions')
         CALL TWRITE('doesn''t correspond to a standard size')
         CALL QSBELL
      END IF
C
C======  Convert to um units during interaction
C
 30   RASTER_SIZE = RASTER_SIZE*1000.
      CALL TASKF('Enter raster size (microns)',RASTER_SIZE,1,QUIT)
      RASTER_SIZE = RASTER_SIZE/1000.
      RETURN
      END
**==hcf.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       HCF           **
C*************************
C
C
      INTEGER    FUNCTION  HCF (II,JJ)
C
C Purpose: To calculate h.c.f. of parameters
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    II, JJ
C
C====== LOCALS:
C
      INTEGER    I, J, ISAV, K
C
C------------------------------------------------------------------------------
C
      I = II
      J = JJ
C
      IF (I.EQ.J) THEN
         HCF = I
         RETURN
      END IF
C
      IF (I.LT.J) THEN
         ISAV = I
         I = J
         J = ISAV
      END IF
C
      IF (J.EQ.0) THEN
         HCF = I
         RETURN
      END IF
C
 20   K = MOD(I,J)
C
      IF (K.EQ.0) THEN
         HCF = J
         RETURN
      END IF
C
      I = J
C
      J = K
      GO TO 20
C
      END
**==lambda_hist.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       LAMBDA_HIST       **
C*****************************
C
C
      SUBROUTINE  LAMBDA_HIST (NUM_REFLEX, LAMBDAS,
     +                         WAVSTD, WAVMIN, WAVMAX,
     +                         N_BIN)
C
C Purpose: Calculate and display wavelength histogram
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER    NUM_REFLEX
      REAL       LAMBDAS(NUM_REFLEX)
      REAL       WAVSTD, WAVMIN, WAVMAX
      INTEGER    N_BIN(16)
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, J
      REAL       BIN_CUTOFF(16), LAMBDA
C
C------------------------------------------------------------------------------
C
C====== Set up bins
C
      DO 20 I = 1 , 16
         BIN_CUTOFF(I) = WAVMIN + I*(WAVMAX-WAVMIN)/15.0
         N_BIN(I) = 0
 20   CONTINUE
      BIN_CUTOFF(16) = 100.0
C
C====== Work out histogram
C
      DO 50 J = 1 , NUM_REFLEX
C
         LAMBDA = LAMBDAS(J)
         DO 30 I = 1 , 16
            IF (LAMBDA.LT.BIN_CUTOFF(I)) GO TO 40
 30      CONTINUE
         I = 16
C
 40      N_BIN(I) = N_BIN(I) + 1
C
 50   CONTINUE
C
C====== Display results
C
      CALL TWRITE(' ')
      CALL TWRITE('Wavelength  :   Number')
      CALL TWRITE(' ')
C
      DO 60 I = 1 , 16
         WRITE (STR,6010) BIN_CUTOFF(I) , N_BIN(I)
         CALL TWRITE(STR)
 60   CONTINUE
C
      CALL TWRITE(' ')
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,F9.3,I12)
      END
**==make_orient.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       MAKE_ORIENT        **
C******************************
C
C
      SUBROUTINE  MAKE_ORIENT (CELL_PARAMS, XTAL_SYSTEM, X_AXIS,
     +               Z_AXIS, ORIENTATION)
C
C Purpose: Makes orientation matrix given recip. cell parameters,
C          the crystal system number (1..7), and alignment of
C          recip. cell lattice wrt. beam and spindle axes.
C Optimisation: Don`t re-calculate if cell params are the same from
C               call to call.
C !!!! Warning !!!! Assumes that cell params will change if one
C                   of xtal-system,...etc do.
C
      INCLUDE 'implicit_none.finc'
C
      REAL       CELL_PARAMS(6)
      INTEGER    XTAL_SYSTEM, X_AXIS, Z_AXIS
      REAL       ORIENTATION(3,3)
C
C Notes: Generalised matrix builder
C
C    The reciprocal cell parameters come in cell_params(1..6) but only the
C    unique ones ... cell_params(1) is astar, cell_params(2) is alphastar
C    for rhombohedral for example.
C    orientation returns the orientation matrix
C      ... row 1 is components of reciprocal cell along X,
C      ... row 2 along Y, row 3 along Z.
C    Convention for crystal system:
C         1 - triclinic    2 - monoclinic    3 - orthorhombic
C         4 - tetragonal   5 - hexagonal     6 - rhombohedral
C         7 - cubic
C    Conventions for alignment are:
C       1..astar,   2..bstar,   3..cstar
C      -1..-astar, -2..-bstar, -3..-cstar
C    Thus if the astar axis lies along the X-axis in the direction of +x
C    then X_axis = 1, and similarly for Z_axis.
C
C
C====== LOCALS:
C
      REAL       OLD_CELL_PARAMS(6), OLD_ORIENT(3,3)
      LOGICAL    CALLED_ONCE
      SAVE       OLD_CELL_PARAMS, CALLED_ONCE, OLD_ORIENT
      REAL       RL(3), CC(3), SS(3)
      INTEGER    I, J, K, SI, SJ, SK, ILOOP
      REAL       COS_PSI, COS_RHO, A, C, S, FUN
      INTEGER    IANG
      REAL       DEG_TO_RAD
      REAL       ALPHA, BETA, GAMMA
      DATA       CALLED_ONCE/.FALSE./
C
C------------------------------------------------------------------------------
C
      IF (CALLED_ONCE) THEN
C
C====== See if params have changed - just return old orient if they haven`t
C
         DO 20 I = 1 , 6
            IF (CELL_PARAMS(I).NE.OLD_CELL_PARAMS(I)) GO TO 50
 20      CONTINUE
         DO 40 I = 1 , 3
            DO 30 J = 1 , 3
               ORIENTATION(I,J) = OLD_ORIENT(I,J)
 30         CONTINUE
 40      CONTINUE
         RETURN
      END IF
C
 50   CALLED_ONCE = .TRUE.
      DO 60 I = 1 , 6
         OLD_CELL_PARAMS(I) = CELL_PARAMS(I)
 60   CONTINUE
C
      DO 80 I = 1 , 3
         DO 70 J = 1 , 3
            ORIENTATION(I,J) = 0.
 70      CONTINUE
 80   CONTINUE
C
C====== i, j, k have special meanings in this code.
C====== i refers to the crystal rotation axis, Z
C
      I = IABS(Z_AXIS)
      SI = 1
      IF (Z_AXIS.LT.0) SI = -1
C
C====== j refers to X_axis, the x-ray beam axis, X
C
      J = IABS(X_AXIS)
      SJ = 1
      IF (X_AXIS.LT.0) SJ = -1
C
C====== k is therefore the other axis, Y
C
      K = 6 - I - J
C
C====== The sign of k relies on i,j,k and X,Y,Z being right-handed
C
      IF (K.EQ.J+1 .OR. K.EQ.J-2) THEN
         SK = SI*SJ
      ELSE
         SK = -SI*SJ
      END IF
C
C====== Zero orientation matrix
C
      CALL MTX_ZERO(ORIENTATION,3,3)
C
C====== Convert cell angles to radians
C
      IANG = 4
C
C====== Check if rhombohedral
C
      IF (XTAL_SYSTEM.EQ.6) IANG = 2
C
      DEG_TO_RAD = 3.1415926536/180.0
      ALPHA = CELL_PARAMS(IANG)*DEG_TO_RAD
      BETA = CELL_PARAMS(IANG+1)*DEG_TO_RAD
      GAMMA = CELL_PARAMS(IANG+2)*DEG_TO_RAD
C
C====== Now go and do specific stuff for each crystal system
C
      IF (XTAL_SYSTEM.EQ.2) THEN
C
C====== Monoclinic
C
         RL(1) = CELL_PARAMS(1)
         RL(2) = CELL_PARAMS(2)
         RL(3) = CELL_PARAMS(3)
C
         CC(1) = 0.0
         CC(2) = COS(ALPHA)
         CC(3) = 0.0
C
         SS(1) = 1.0
         SS(2) = SIN(ALPHA)
         SS(3) = 1.0
C
         ORIENTATION(1,J) = SJ*RL(J)*SS(K)
         ORIENTATION(1,K) = SJ*RL(K)*CC(I)
         ORIENTATION(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENTATION(3,I) = SI*RL(I)
         ORIENTATION(3,J) = SI*RL(J)*CC(K)
         ORIENTATION(3,K) = SI*RL(K)*CC(J)
      ELSE IF (XTAL_SYSTEM.EQ.3) THEN
C
C====== Orthorhombic
C
         RL(1) = CELL_PARAMS(1)
         RL(2) = CELL_PARAMS(2)
         RL(3) = CELL_PARAMS(3)
C
         ORIENTATION(1,J) = SJ*RL(J)
         ORIENTATION(2,K) = SK*RL(K)
         ORIENTATION(3,I) = SI*RL(I)
      ELSE IF (XTAL_SYSTEM.EQ.4) THEN
C
C====== Tetragonal
C
         RL(1) = CELL_PARAMS(1)
         RL(2) = CELL_PARAMS(1)
         RL(3) = CELL_PARAMS(2)
C
         ORIENTATION(1,J) = SJ*RL(J)
         ORIENTATION(2,K) = SK*RL(K)
         ORIENTATION(3,I) = SI*RL(I)
      ELSE IF (XTAL_SYSTEM.EQ.5) THEN
C
C====== Hexagonal
C
         RL(1) = CELL_PARAMS(1)
         RL(2) = CELL_PARAMS(1)
         RL(3) = CELL_PARAMS(2)
C
         CC(1) = 0.0
         CC(2) = 0.0
         CC(3) = 0.5
C
         SS(1) = 1.0
         SS(2) = 1.0
         SS(3) = SQRT(3.0)/2.0
C
         ORIENTATION(1,J) = SJ*RL(J)*SS(K)
         ORIENTATION(1,K) = SJ*RL(K)*CC(I)
         ORIENTATION(2,K) = SK*RL(K)*SS(I)*SS(J)
         ORIENTATION(3,I) = SI*RL(I)
         ORIENTATION(3,J) = SI*RL(J)*CC(K)
         ORIENTATION(3,K) = SI*RL(K)*CC(J)
      ELSE IF (XTAL_SYSTEM.EQ.6) THEN
C
C====== Rhombohedral
C
         A = CELL_PARAMS(1)
         C = COS(ALPHA)
         S = SIN(ALPHA)
         FUN = SQRT(1.-3.*C**2+2.*C**3)
C
         ORIENTATION(1,J) = SJ*A*S
         ORIENTATION(1,K) = SJ*A*C*(1.-C)/S
         ORIENTATION(2,K) = SK*A*FUN/S
         ORIENTATION(3,I) = SI*A
         ORIENTATION(3,J) = SI*A*C
         ORIENTATION(3,K) = SI*A*C
      ELSE IF (XTAL_SYSTEM.EQ.7) THEN
C
C====== Cubic
C
         A = CELL_PARAMS(1)
C
         ORIENTATION(1,J) = SJ*A
         ORIENTATION(2,K) = SK*A
         ORIENTATION(3,I) = SI*A
      ELSE
C
C====== Triclinic
C
         RL(1) = CELL_PARAMS(1)
         RL(2) = CELL_PARAMS(2)
         RL(3) = CELL_PARAMS(3)
C
         CC(1) = COS(ALPHA)
         CC(2) = COS(BETA)
         CC(3) = COS(GAMMA)
C
         SS(1) = SIN(ALPHA)
         SS(2) = SIN(BETA)
         SS(3) = SIN(GAMMA)
C
C====== See Buerger's Xray Crystallography for the next formulae, p348
C
         COS_PSI = (CC(I)-CC(K)*CC(J))/SS(K)
         COS_RHO = SQRT(1.0-CC(I)**2-CC(J)**2-CC(K)**2+2.0*CC(I)*CC(J)
     +             *CC(K))/SS(K)
C
         ORIENTATION(1,J) = SJ*RL(J)*SS(K)
         ORIENTATION(1,K) = SJ*RL(K)*COS_PSI
         ORIENTATION(2,K) = SK*RL(K)*COS_RHO
         ORIENTATION(3,I) = SI*RL(I)
         ORIENTATION(3,J) = SI*RL(J)*CC(K)
         ORIENTATION(3,K) = SI*RL(K)*CC(J)
      END IF
C
      DO 100 I = 1 , 3
         DO 90 J = 1 , 3
            OLD_ORIENT(I,J) = ORIENTATION(I,J)
 90      CONTINUE
 100  CONTINUE
      RETURN
      END
**==match_nodals.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**       MATCH_NODALS        **
C*******************************
C
C
      SUBROUTINE  MATCH_NODALS (DELTA, FILM_NAME,
     +                          BOX_SIZE,
     +                          NUM_OLD_OBS, NUM_REF_SPOTS,
     +                          X_REJECT, Y_REJECT, NO_FILM)
C
C Purpose: To pair up spots in prediction with those on film
C
      INCLUDE 'implicit_none.finc'
C
      CHARACTER*40  FILM_NAME
      REAL       DELTA
      REAL       BOX_SIZE
      INTEGER    NUM_OLD_OBS, NUM_REF_SPOTS
      REAL       X_REJECT, Y_REJECT
      LOGICAL    NO_FILM
C
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
      INCLUDE 'backgr.finc'
      INCLUDE 'ramask.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      COMMON /RMS/ OLD_RMS
      REAL   OLD_RMS
 
      INTEGER NBITST
C
C====== LOCALS:
C
      LOGICAL ERR
      CHARACTER*80 STR
      INTEGER    FILM_UNIT
      INTEGER    I, J, IREF, H, K, L, IFAIL
      INTEGER    HKL_KEY(MAX_OBS), WORK1(MAX_OBS), WORK2(MAX_OBS)
      INTEGER    NUM_OVER, NUM_CURS
      INTEGER    THRESH
      REAL       SPOT_DIAM, TOO_FAR
      LOGICAL    QUIT, DO_DISPLAUE, EXISTS
      CHARACTER*40  MATCHES_FILE
      REAL       X_CEN, Y_CEN, OMEGA, Y_SCRAST
      REAL       COS_OMEGA, SIN_OMEGA
      REAL       K_FACTOR, DEG_TO_RAD
      REAL       XF, YF
      INTEGER    BACKGROUND, IRAD
      INTEGER    XC, YC, X_SIZE, Y_SIZE
      REAL       CGX, CGY, XM, YM
      INTEGER    CUR_OBS, NUM_STARRED
      REAL       SUMSQ, DIST2, RMS, DX, DY
      REAL       PLOT_SIZE_X, PLOT_SIZE_Y
      INTEGER    IRECL
      INTEGER    IMASK
      REAL CMASKM, BMINORM
      DATA     SPOT_DIAM /0.3/
C
C------------------------------------------------------------------------------
C
      CALL TASKQ('Radially streaked spots?',STREAK,QUIT)
      IF (STREAK) THEN
         CALL TASKF('Enter streak length (mm). ',CMASKM,3,QUIT)
         CALL TASKF('Enter streak width (mm). ',BMINORM,3,QUIT)
         CALL TASKF('Enter bg border width (rasters)',BORDER,1,QUIT)
         CALL TASKF('Enter constant K ',KMASK,2,QUIT)
         BOX_SIZE = CMASKM + 2.0/MM_RAST_X + BORDER/MM_RAST_X
         SPOT_DIAM = BMINORM
         TOO_FAR = (BOX_SIZE-SPOT_DIAM)/2.0
         CMASK = MM_RAST_X*CMASKM/2.0
         BMINOR = MM_RAST_X*BMINORM/2.0
         DMASK = CTOF*MM_RAST_X
      END IF
C
      IF (NO_FILM) THEN
C
C====== If not using a film, have to use a saved file inside DISPLAUE.
C
         DO_DISPLAUE = .TRUE.
      ELSE
         DO_DISPLAUE = .TRUE.
         CALL TASKQ('Do you want to match spots on a threshold plot',
     +              DO_DISPLAUE,QUIT)
      END IF
C
      IF (DO_DISPLAUE) THEN
         FILM_UNIT = 4
         CALL DISPLAUE(FILM_NAME,FILM_UNIT,NUM_CURS,NO_FILM)
C
C====== Just do an auto match-up if too few spots
C
         IF (NUM_CURS.GE.1) THEN
            DO 20 CUR_OBS = 1 , NUM_CURS
               H = NODAL_H_S(CUR_OBS)
               K = NODAL_K_S(CUR_OBS)
               L = NODAL_L_S(CUR_OBS)
               HKL_KEY(CUR_OBS) = IABS(H) + IABS(K) + IABS(L)
               REF_NODALS(CUR_OBS) = CUR_OBS
 20         CONTINUE
            CUR_OBS = NUM_CURS
            GO TO 40
         END IF
      END IF
C
C====== Get box size multiple
C
      IF (.NOT.STREAK) THEN
         CALL TASKF('Enter spot diameter in mm.',SPOT_DIAM,2,QUIT)
         IF (SPOT_DIAM.GT.DELTA) THEN
            WRITE (STR,6080)
            CALL TWRITE(STR)
            WRITE (STR,6090) DELTA
            CALL TWRITE(STR)
         END IF
C
         BOX_SIZE = 2.0*DELTA - SPOT_DIAM
         IF (BOX_SIZE.LT.SPOT_DIAM) BOX_SIZE = SPOT_DIAM
         WRITE (STR,6100) BOX_SIZE , IFIX(BOX_SIZE*MM_RAST_X)
         CALL TWRITE(STR)
         CALL TASKF('Enter box size in mm.',BOX_SIZE,2,QUIT)
C
C====== Too_far (mm) is max. shift possible with all the spot inside box.
C
         TOO_FAR = (BOX_SIZE-SPOT_DIAM)/2.0
      END IF
C
C====== Get threshold
C
      THRESH = 20
      CALL TASKI('Enter threshold above background for c_of_g calc',
     +           THRESH,QUIT)
C
C====== Open film file
C
      FILM_UNIT = 4
      IF (IPLATE) THEN
         IRECL = 2*Y_RASTERS
         CALL IMPOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      ELSE
         IRECL = Y_RASTERS
         CALL FLMOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      END IF
C
      WRITE (STR,6020)
      CALL TWRITE(STR)
C
      CUR_OBS = 0
C
C====== Prepare for transform from recip. space to scanner space
C
      X_CEN = X_CEN_F + X_C
      Y_CEN = Y_CEN_F + Y_C
      OMEGA = W_F + W_C
C
      DEG_TO_RAD = 3.1415926536/180.0
      COS_OMEGA = COS(OMEGA*DEG_TO_RAD)
      SIN_OMEGA = SIN(OMEGA*DEG_TO_RAD)
C
      Y_SCRAST = MM_RAST_Y
C
      SUMSQ = 0.0
C
      DO 30 I = 1 , NUM_OBS
C
         X_SIZE = MM_RAST_X*BOX_SIZE/2.0
         Y_SIZE = MM_RAST_Y*BOX_SIZE/2.0
C
         IREF = REF_NODALS(I)
C
C====== Reject if spatially overlapped
C
         IF (NBITST(M_FLAGS(IREF),6,2).EQ.0) THEN
C
            XF = XF_S(IREF)
            YF = YF_S(IREF)
C
C====== K_factor = 1.0 + tilt*xf + twist*yf + bulge*sqrt(xf*xf + yf*yf)
C
            K_FACTOR = 1.0
C
            XM = X_CEN + K_FACTOR*(XF*COS_OMEGA-YF*SIN_OMEGA)*MM_RAST_X
            YM = Y_CEN + K_FACTOR*(XF*SIN_OMEGA+YF*COS_OMEGA)*Y_SCRAST
C
C====== Reject if in backstop region
C
            IF (ABS(XM-X_CEN_F).GE.X_REJECT*10.0) THEN
               IF (ABS(YM-Y_CEN_F).GE.Y_REJECT*10.0) THEN
C
C====== Round to nearest integer to be in scanner space
C
                  XC = XM + 0.5
                  YC = YM + 0.5
C
C====== Overestimate background
C
                  IRAD = SQRT(FLOAT((XC-XCEN)**2+(YC-YCEN)**2))
                  J = 10 - (IRAD/N_PER_BIN)
                  J = MAX(1,MIN(10,J))
                  BACKGROUND = BACKGR(J) + THRESH
C
                  CALL C_OF_G(FILM_UNIT,XC,YC,X_SIZE,Y_SIZE,BACKGROUND,
     +                        CGX,CGY,NUM_OVER)
                  IF (NUM_OVER.GT.0) THEN
                     DIST2 = (XM-CGX)**2 + (YM-CGY)**2
C
C====== Check each coord to see if some of spot is outside c of g box.
C
                     IF (ABS(XM-CGX).GT.(TOO_FAR*MM_RAST_X)) THEN
C
C====== Because some of spot is outside box, expect true c of g to be
C====== further towards edge of box. An approximation of distance
C====== to move is the amount of spot that`s been cut off.
C====== ie. (cgx-xm) + spot_diam/2 - box_size/2 (put into right units)
C
                        XC = CGX + (CGX-XM)
     +                       + ((SPOT_DIAM-BOX_SIZE)/2.*MM_RAST_X)
                     ELSE
                        XC = CGX
                     END IF
C
                     IF (ABS(YM-CGY).GT.(TOO_FAR*MM_RAST_Y)) THEN
                        YC = CGY + (CGY-YM)
     +                       + ((SPOT_DIAM-BOX_SIZE)/2.*MM_RAST_Y)
                     ELSE
                        YC = CGY
                     END IF
C
C====== Re-centre if we didn`t get a true centre of gravity.
C
                     IF (ABS(XM-CGX).GT.(TOO_FAR*MM_RAST_X) .OR. 
     +                   ABS(YM-CGY).GT.(TOO_FAR*MM_RAST_Y)) THEN
C
C====== True c of g should be within a spot radius of here, but
C====== make box 1.5 times bigger just in case.
C
                        X_SIZE = SPOT_DIAM/2.0*1.5*MM_RAST_X
                        Y_SIZE = SPOT_DIAM/2.0*1.5*MM_RAST_Y
                        CALL C_OF_G(FILM_UNIT,XC,YC,X_SIZE,Y_SIZE,
     +                              BACKGROUND,CGX,CGY,NUM_OVER)
                        DIST2 = (XM-CGX)**2 + (YM-CGY)**2
                     END IF
                  END IF
C
                  IF (NUM_OVER.GT.0) THEN
C
                     CUR_OBS = CUR_OBS + 1
                     XM_OBS(CUR_OBS) = CGX
                     YM_OBS(CUR_OBS) = CGY
                     REF_NODALS(CUR_OBS) = IREF
                     H = H_S(IREF)
                     K = K_S(IREF)
                     L = L_S(IREF)
                     NODAL_H_S(CUR_OBS) = H
                     NODAL_K_S(CUR_OBS) = K
                     NODAL_L_S(CUR_OBS) = L
                     HKL_KEY(CUR_OBS) = IABS(H) + IABS(K) + IABS(L)
C
                     CALL CALC_XYM(H,K,L,XM,YM)
                     SUMSQ = SUMSQ + DIST2
C
                  END IF
               END IF
            END IF
         END IF
C
 30   CONTINUE
C
      IF (IPLATE) THEN
         CALL IMPCLS(FILM_UNIT)
      ELSE
         CALL FLMCLS(FILM_UNIT)
      END IF
C
C====== Give rms for all spots
C
      RMS = SQRT(SUMSQ/CUR_OBS)/MM_RAST_X
      CALL TWRITE(
     +        '********************************************************'
     +        )
      WRITE (STR,6050) RMS , CUR_OBS
      CALL TWRITE(STR)
      CALL TWRITE(
     +        '********************************************************'
     +        )
C
C====== Get number of spots to refine against
C
 40   WRITE (STR,6030) CUR_OBS
      CALL TWRITE(STR)
      NUM_REF_SPOTS = CUR_OBS
      CALL TASKI('Enter number of spots to be used in refinement',
     +           NUM_REF_SPOTS,QUIT)
      IF (NUM_REF_SPOTS.LT.20) THEN
         WRITE (STR,6110)
         CALL TWRITE(STR)
         NUM_REF_SPOTS = CUR_OBS
      END IF
      NUM_REF_SPOTS = MIN(NUM_REF_SPOTS,CUR_OBS)
C
C====== Index reflections on hkl
C
      CALL INDXIV(HKL_KEY,CUR_OBS,REF_NODALS)
C
C====== Set up for plot
C
      PLOT_SIZE_X = X_RASTERS
      PLOT_SIZE_Y = Y_RASTERS
      CALL PREPARE_PLOT(1.,PLOT_SIZE_X,1.,PLOT_SIZE_Y)
      DO 50 I = 1 , 3
         CALL QGSYMB(FIDX(I),FIDY(I),4,2.0)
 50   CONTINUE
C
      NUM_STARRED = 0
      SUMSQ = 0.0
C
      DO 60 I = 1 , NUM_REF_SPOTS
C
         J = REF_NODALS(I)
C
         H = NODAL_H_S(J)
         K = NODAL_K_S(J)
         L = NODAL_L_S(J)
         CALL CALC_XYM(H,K,L,XM,YM)
C
C====== Check how near we are to the spot on the film
C
         DIST2 = (XM-XM_OBS(J))**2 + (YM-YM_OBS(J))**2
         IF (DIST2.GT.(MM_RAST_X*DELTA)**2) THEN
C
C====== More than delta ? - yes, then flag by drawing a star
C
            CALL QGSYMB(XM_OBS(J),YM_OBS(J),1,3.0)
            CALL QGSYMB(XM_OBS(J),YM_OBS(J),2,3.0)
            NUM_STARRED = NUM_STARRED + 1
C
         ELSE
C
            CALL QGSYMB(XM_OBS(J),YM_OBS(J),3,1.0)
C
         END IF
C
C====== Exaggerate vector x 5
C
         DX = 5.0*(XM-XM_OBS(J))
         DY = 5.0*(YM-YM_OBS(J))
         CALL QGMOVE(XM_OBS(J),YM_OBS(J))
         CALL QGDRAW(XM_OBS(J)+DX,YM_OBS(J)+DY)
         SUMSQ = SUMSQ + DIST2
C
 60   CONTINUE
C
      CALL QGBELL
      CALL TREAD(STR,QUIT)
      CALL QGALPH
      RMS = SQRT(SUMSQ/NUM_REF_SPOTS)/MM_RAST_X
C
C====== Save RMS for comparison after ref`mnt
C
      OLD_RMS = RMS
C
      WRITE (STR,6040) NUM_STARRED
      CALL TWRITE(STR)
C
C====== Give rms for only those spots to be used in refinement
C
      WRITE (STR,6060) RMS , NUM_REF_SPOTS
      CALL TWRITE(STR)
C
      WRITE (STR,6010) CUR_OBS , NUM_OBS
      CALL TWRITE(STR)
C
      NUM_OBS = CUR_OBS
      NUM_OLD_OBS = NUM_OBS
C
      IF (NUM_OBS.LT.20) THEN
         WRITE (STR,6070)
         CALL TWRITE(STR)
      END IF
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('...',I6,' spots located on film out of',I6,
     +        ' predicted and stored')
 6020 FORMAT ('===Matching up predicted and observed spots...')
 6030 FORMAT ('...',I6,' spots located')
 6040 FORMAT ('...',I6,' spots were starred in difference map')
 6050 FORMAT ('*****   rms deviation:',F7.3,' mm. for ',I3,
     +        ' spots    *****')
 6060 FORMAT ('rms deviation: ',F7.3,' mm.  for the ',I3,
     +        ' spots to be used in refinement')
 6070 FORMAT ('** Less than 20 spots located for refinement **')
 6080 FORMAT ('** Spot diameter greater than delta **')
 6090 FORMAT ('Currently delta =',F5.2,' mm.')
 6100 FORMAT ('Suggested box size - ',F5.2,' mm.  (ie. ',I3,
     +        ' rasters across)')
 6110 FORMAT ('Can''t refine with less than 20 spots - ',
     +        'so using all spots')
      END
**==nodal_index.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C******************************
C**       NODAL_INDEX        **
C******************************
C
C
      INTEGER   FUNCTION NODAL_INDEX (H, K, L)
C
C Purpose: To calculate the "nodal spot selection index" at which a
C          nodal would be selected.
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER   H, K, L
C
C====== LOCALS:
C
      INTEGER   HCF, HCF1, HCF2
      INTEGER   H1, K1, L1
      INTEGER   HI, LO
C
C------------------------------------------------------------------------------
C
      H1 = H
      K1 = K
      L1 = L
      HCF1 = HCF(H1,K1)
      HCF2 = HCF(HCF1,L1)
      H1 = H1/HCF2
      K1 = K1/HCF2
      L1 = L1/HCF2
C
      HI = MAX0(H1,K1,L1)
      LO = MIN0(H1,K1,L1)
      NODAL_INDEX = MAX0(HI,-LO)
C
      RETURN
      END
**==nodein.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       NODEIN        **
C*************************
C
      SUBROUTINE NODEIN(NUM_CURS,COFG)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose : use the cursor to input positions in pairs
C           1 the nodal spot (for hkl)
C           2 the associated spot on the film (for x,y)
C
C     If aborted by user return -1 in num_curs
C
      INTEGER NUM_CURS,COFG
C
      INCLUDE 'backgr.finc'
      INCLUDE 'observed_spots.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'ramask.finc'
C
      COMMON /NODAL/NTOT,NPLOTTED,NODE
      INTEGER NTOT,NPLOTTED,NODE(5,1000)
      COMMON /PLOTS/SCAN_START_X,SCAN_START_Y,IDX,IDY
      INTEGER IDX,IDY,SCAN_START_X,SCAN_START_Y
C
C Author : PAM
C
C====== LOCALS:
C
      LOGICAL QUIT
      CHARACTER*80 STR
      INTEGER I,II,IBACK,IX,IY
      REAL XX,YY,XM,YM,XP,YP
      INTEGER DIF,DELX,DELY,IRAD
      CHARACTER*1 KCH
      INTEGER FILM_UNIT,NUM_OVER
      LOGICAL DO_DUMP,REPEAT
      INTEGER IXBOX,IYBOX
      DATA DIF/5/
C
C-----------------------------------------------------------------------
C
      FILM_UNIT = 4
C
C====== Get nodal spots
C
      NUM_CURS = 0
      REPEAT = .FALSE.
C
C====== Now input from the cursor
C====== first a nodal spot
C
 20   IF (REPEAT) THEN
C
C====== Repeat previous spot indices- inc count of spots
C
         NUM_CURS = NUM_CURS + 1
         IF (NUM_CURS.EQ.(NUM_CURS/5)*5) THEN
            WRITE (STR,6010) NUM_CURS
            CALL TLOG(.FALSE.)
            CALL TWRITE(STR)
            CALL TWRITE(' ')
            CALL TLOG(.TRUE.)
         END IF
C
         NODAL_H_S(NUM_CURS) = NODAL_H_S(NUM_CURS-1)
         NODAL_K_S(NUM_CURS) = NODAL_K_S(NUM_CURS-1)
         NODAL_L_S(NUM_CURS) = NODAL_L_S(NUM_CURS-1)
 
      ELSE
         CALL QGCURS(XM,YM,KCH)
         CALL QGMPXY(XM,YM,XP,YP)
         IX = NINT(XP)
         IY = NINT(YP)
C
C====== See if user wants to end or abort
C
         IF (KCH.EQ.'E' .OR. KCH.EQ.'e') GO TO 50
         IF (KCH.EQ.'A' .OR. KCH.EQ.'a') THEN
C
C====== "A" - abort
C
            CALL QGALPH
            NUM_CURS = -1
            RETURN
 
         END IF
 
         DO 30 I = 1 , NPLOTTED
            DELX = NODE(4,I) - IX
            IF (ABS(DELX).LE.DIF) THEN
               DELY = NODE(5,I) - IY
               IF (ABS(DELY).LE.DIF) THEN
C
C====== Spot OK - inc count of spots
C
                  IF (NUM_CURS.GE.MAX_OBS) GO TO 50
                  NUM_CURS = NUM_CURS + 1
                  IF (NUM_CURS.EQ.(NUM_CURS/5)*5) THEN
                     WRITE (STR,6010) NUM_CURS
                     CALL TLOG(.FALSE.)
                     CALL TWRITE(STR)
                     CALL TWRITE(' ')
                     CALL TLOG(.TRUE.)
                  END IF
C
                  NODAL_H_S(NUM_CURS) = NODE(1,I)
                  NODAL_K_S(NUM_CURS) = NODE(2,I)
                  NODAL_L_S(NUM_CURS) = NODE(3,I)
                  GO TO 40
               END IF
            END IF
 
 30      CONTINUE
C
C====== No match
C
         WRITE (STR,6020)
         CALL TLOG(.FALSE.)
         CALL TWRITE(STR)
         CALL TLOG(.TRUE.)
         GO TO 20
 
      END IF
 
C
C====== Now input observed spot position
C
 40   CALL QGBELL
      CALL QGCURS(XM,YM,KCH)
      CALL QGMPXY(XM,YM,XP,YP)
      IX = NINT(XP)
      IY = NINT(YP)
      IF (KCH.EQ.'r' .OR. KCH.EQ.'R') REPEAT = .TRUE.
      IF (KCH.EQ.'s' .OR. KCH.EQ.'S') REPEAT = .FALSE.
      IF (COFG.EQ.2) THEN
C
C====== Box input mode.
C
C====== Calculate background using radial curve and spot position
C====== change to scanner units
C
         XX = IX*IDX + SCAN_START_X
         YY = IY*IDY + SCAN_START_Y
C
C====== NB Background is centred on (xcen,ycen)
C
         IRAD = SQRT((XX-XCEN)**2+(YY-YCEN)**2)
         II = N_BACKGR - (IRAD/N_PER_BIN)
 
C
C====== Overestimate background
C
         IBACK = BACKGR(II) + 10
         IX = XX
         IY = YY
         IF (.NOT.STREAK) THEN
            CALL C_OF_G(FILM_UNIT,IX,IY,20,20,IBACK,XM_OBS(NUM_CURS),
     +                  YM_OBS(NUM_CURS),NUM_OVER)
C
         ELSE
            IXBOX = 2.0*CMASK + 2.0 + BORDER
            IYBOX = IXBOX
            CALL C_OF_G(FILM_UNIT,IX,IY,IXBOX,IYBOX,IBACK,
     +                  XM_OBS(NUM_CURS),YM_OBS(NUM_CURS),NUM_OVER)
         END IF
         IF (NUM_OVER.LE.0) THEN
            WRITE (STR,6030)
            CALL TLOG(.FALSE.)
            CALL TWRITE(STR)
            CALL TLOG(.TRUE.)
            GO TO 40
 
         END IF
 
      ELSE
C
C====== Simple input mode.
C
C====== Change to scanner units
C
         XM_OBS(NUM_CURS) = IX*IDX + SCAN_START_X
         YM_OBS(NUM_CURS) = IY*IDY + SCAN_START_Y
 
      END IF
      GO TO 20
C
C====== End of input
C
 50   IF (NUM_CURS.LT.20) THEN
         CALL QGTMOD
         WRITE (STR,6040) NUM_CURS
         CALL TLOG(.FALSE.)
         CALL TWRITE(STR)
         CALL TLOG(.TRUE.)
 
      END IF
 
      CALL QGBELL
      CALL TREAD(STR,QUIT)
      CALL QGALPH
C
C====== Dump of spot positions
C
      CALL TASKQ('Dump spots',DO_DUMP,QUIT)
      IF (DO_DUMP) CALL DUMP_SPOTS(NUM_CURS)
      RETURN
C
C====== Format statements
C
 6010 FORMAT (T66,I2,'spots')
 6020 FORMAT ('no match , try again ')
 6030 FORMAT ('Flat film, try again')
 6040 FORMAT ('WARNING - refinement may not be stable with',' only ',I3,
     +        ' spots.')
      END
**==plot_diff.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************
C**       PLOT_DIFF        **
C****************************
C
C
      SUBROUTINE  PLOT_DIFF (FILM_NAME, DELTA, BOX_SIZE,
     + NUM_OLD_OBS)
C
C Purpose: Plot vector differences between observed and predicted spot positions
C
      INCLUDE 'implicit_none.finc'
C
      CHARACTER*40  FILM_NAME
      REAL       DELTA, BOX_SIZE
      INTEGER    NUM_OLD_OBS
C
      INCLUDE 'observed_spots.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'backgr.finc'
      INCLUDE 'sizes.finc'
C
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      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, ERR
      CHARACTER*1 ANSWER
      CHARACTER*80 STR
      INTEGER    FILM_UNIT
      INTEGER    I, H, K, L, J, NUM_BOXED
      REAL       XM, YM
      INTEGER    BACKGROUND, IRAD, NUM_OVER
      INTEGER    XC, YC, X_SIZE, Y_SIZE
      REAL       CGX, CGY, RMS, SUMSQ, DIST2
      REAL       DX, DY
      REAL       PLOT_SIZE_X, PLOT_SIZE_Y
      INTEGER    IRECL, IFAIL
C
C------------------------------------------------------------------------------
C
      PLOT_SIZE_X = X_RASTERS
      PLOT_SIZE_Y = Y_RASTERS
      CALL PREPARE_PLOT(1.,PLOT_SIZE_X,1.,PLOT_SIZE_Y)
C
C====== Draw in fiducials
C
      DO 20 I = 1 , 3
         CALL QGSYMB(FIDX(I),FIDY(I),1,3.0)
 20   CONTINUE
C
      X_SIZE = MM_RAST_X*BOX_SIZE/2.0
      Y_SIZE = MM_RAST_Y*BOX_SIZE/2.0
C
C====== Open film file
C
      FILM_UNIT = 4
      IF (IPLATE) THEN
         IRECL = 2*Y_RASTERS
         CALL IMPOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      ELSE
         IRECL = Y_RASTERS
         CALL FLMOPN(-FILM_UNIT,FILM_NAME,IRECL,1,OUT_TERM,ERR)
      END IF
C
      NUM_BOXED = 0
      SUMSQ = 0.0
C
      DO 30 I = 1 , NUM_OLD_OBS
C
         H = NODAL_H_S(I)
         K = NODAL_K_S(I)
         L = NODAL_L_S(I)
         CALL CALC_XYM(H,K,L,XM,YM)
C
         DIST2 = (XM-XM_OBS(I))**2 + (YM-YM_OBS(I))**2
         SUMSQ = SUMSQ + DIST2
C
         XC = XM + 0.5
         YC = YM + 0.5
C
         IRAD = SQRT(FLOAT((XC-XCEN)**2+(YC-YCEN)**2))
         J = 10 - (IRAD/N_PER_BIN)
         J = MAX(1,MIN(10,J))
         BACKGROUND = BACKGR(J) + 10
C
         CALL C_OF_G(FILM_UNIT,XC,YC,X_SIZE,Y_SIZE,BACKGROUND,CGX,CGY,
     +               NUM_OVER)
C
C====== Check if we are still near the old spot on the film
C
         IF (((CGX-XM_OBS(I))/MM_RAST_X)**2+((CGY-YM_OBS(I))/MM_RAST_Y)
     +       **2.GT.(DELTA)**2) THEN
C
C====== NO!! - then flag by drawing a box centred on old position
C
            CALL QGSYMB(XM_OBS(I),YM_OBS(I),3,3.0)
            NUM_BOXED = NUM_BOXED + 1
C
         ELSE
C
            CALL QGSYMB(XM_OBS(I),YM_OBS(I),3,1.0)
C
         END IF
C
C====== Exaggerate vector x 5
C
         DX = 5.0*(XM-XM_OBS(I))
         DY = 5.0*(YM-YM_OBS(I))
         CALL QGMOVE(XM_OBS(I),YM_OBS(I))
         CALL QGDRAW(XM_OBS(I)+DX,YM_OBS(I)+DY)
C
 30   CONTINUE
C
      IF (IPLATE) THEN
         CALL IMPCLS(FILM_UNIT)
      ELSE
         CALL FLMCLS(FILM_UNIT)
      END IF
C
C====== Pause for appreciation
C
      CALL QGBELL
      CALL TREAD(ANSWER,QUIT)
      CALL QGALPH
C
      RMS = SQRT(SUMSQ/NUM_OLD_OBS)/MM_RAST_X
C
      WRITE (STR,6010) NUM_BOXED
      CALL TWRITE(STR)
      CALL TWRITE(
     +        '********************************************************'
     +        )
      WRITE (STR,6020) RMS , NUM_OLD_OBS
      CALL TWRITE(STR)
      CALL TWRITE(
     +        '********************************************************'
     +        )
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('...',I6,' spots were boxed in difference map')
 6020 FORMAT ('*****   rms deviation:',F7.3,' mm. for ',I3,
     +        ' spots    *****')
      END
**==plotge.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       PLOTGE        **
C*************************
C
C
      SUBROUTINE PLOTGE(ITYPE)
C
      INCLUDE 'implicit_none.finc'
C
C Purpose : read spots from generated array, plot all or nodals according
C           to itype.
C            Convert coords from mm, central origin to plotter space
C            ie 700 by 700 origin bottom left
C            NODE stores plotter units.
C
      INTEGER ITYPE
C
      INCLUDE 'fiducials.finc'
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'film_params.finc'
      COMMON /NODAL/NTOT,NPLOTTED,NODE
      INTEGER NTOT,NPLOTTED,NODE(5,1000)
      COMMON /PLOTS/SCAN_START_X,SCAN_START_Y,IDX,IDY
      INTEGER IDX,IDY,SCAN_START_X,SCAN_START_Y
      COMMON /OMEGA/COSOM,SINOM
      REAL COSOM,SINOM
C
C Author : PAM
C
C====== LOCALS:
C
      INTEGER LX,LY,I,IH,IK,IL
      REAL X_CEN,Y_CEN,XM,YM,X,Y
      INTEGER JX,JY,K_FACTOR
      INTEGER NUM_TO_PLOT,IREF
      DATA LX,LY/3,3/
C
C-----------------------------------------------------------------------
C
      NPLOTTED = 0
      NTOT = 0
      X_CEN = X_CEN_F + X_C
      Y_CEN = Y_CEN_F + Y_C
C
C====== This factor used to apply twist, tilt, bulge.
C
      K_FACTOR = 1.0
C
C====== Loop over reflexions. If plotting nodals, find spots with
C====== the "ref_nodals" array of pointers.
C
      IF (ITYPE.EQ.1) THEN
         NUM_TO_PLOT = NUM_OBS
 
      ELSE
         NUM_TO_PLOT = NUM_REFLEX
      END IF
 
      DO 20 I = 1 , NUM_TO_PLOT
         IF (ITYPE.EQ.1) THEN
            IREF = REF_NODALS(I)
 
         ELSE
            IREF = I
         END IF
 
         X = XF_S(IREF)
         Y = YF_S(IREF)
C
C====== Transform from mms to scanner space
C
         X = X*MM_RAST_X
         Y = Y*MM_RAST_Y
C
C====== Adjust wrt centre and omega, convert to mm.s
C
         XM = K_FACTOR*(X*COSOM-Y*SINOM)
         YM = K_FACTOR*(X*SINOM+Y*COSOM)
         JX = (XM+X_CEN-SCAN_START_X)/IDX + 1
         JY = (YM+Y_CEN-SCAN_START_Y)/IDY + 1
         IH = H_S(IREF)
         IK = K_S(IREF)
         IL = L_S(IREF)
         NPLOTTED = NPLOTTED + 1
         IF (NPLOTTED.GT.1000) GO TO 30
 
         NODE(1,NPLOTTED) = IH
         NODE(2,NPLOTTED) = IK
         NODE(3,NPLOTTED) = IL
         NODE(4,NPLOTTED) = JX
         NODE(5,NPLOTTED) = JY
C
C====== Note that calculated spot positions are stored in
C====== plotter space units
C
C
C====== Then output cross at point
C
         CALL QGSYMB(FLOAT(JX),FLOAT(JY),1,3.0)
 
 20   CONTINUE
 30   RETURN
      END
**==prepare_plot.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*******************************
C**       PREPARE_PLOT        **
C*******************************
C
C
      SUBROUTINE  PREPARE_PLOT (X1, X2, Y1, Y2)
C
C Purpose: To set up VT100 screen for plotting pattern
C
      INCLUDE 'implicit_none.finc'
C
      REAL    X1, X2, Y1, Y2
C
C------------------------------------------------------------------------------
C
      CALL QGENTR
      CALL QGCLR
      CALL QGMAPC(0.0,767.0,0.0,767.0,X1,X2,Y1,Y2)
      CALL QGMOVE(X1,Y1)
      CALL QGDRAW(X1,Y2)
      CALL QGDRAW(X2,Y2)
      CALL QGDRAW(X2,Y1)
      CALL QGDRAW(X1,Y1)
      RETURN
      END
**==radmask.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**        RADMASK          **
C*****************************
C
C
C  Program designed to create a mask for radial streaks
C  This version for genlaue includes commons from there for
C  xcentre and ycentre
C  Not really similar in many ways to INTLAUE, since here we
C  test a single point xpnt, ypnt related to a centre
C  ixmin,iymin. INTLAUE generates the whole mask
C
      SUBROUTINE RADMASK(IRECM,IXMIN,IYMIN,XPNT,YPNT,ISPOTR,IMARK)

      INTEGER IRECM,IXMIN,IYMIN,XPNT,YPNT,ISPOTR,IMARK
      INTEGER XPN,YPN
C
C
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'ramask.finc'
C
C-------------------------------------------------------------------------------
C
C====== Incoming ixmin,iymin,xpnt,ypnt are in rasters, positive for irecm=1
C====== (this is also the case for genlaue. Note that the call from
C====== cgmask in seekrs has irecm = 1 and ispotr set to 9111 to
C====== indicate that mask needs rearranging as in genlaue)
C====== convert to rasters about centre using fact (0.2 for 50micron)
C====== or in case of genlaue, using mm_rast
C
C
      IXM = ANINT(FLOAT(IXMIN)-X_CEN_F-X_C)
      IYM = ANINT(FLOAT(IYMIN)-Y_CEN_F-Y_C)
      XPN = ANINT(FLOAT(XPNT)-X_CEN_F-X_C)
      YPN = ANINT(FLOAT(YPNT)-Y_CEN_F-Y_C)
 
C
C====== cmask is constant eta*D or FULL spot size (mm) for k=0
C====== bminor is full width in mm
C====== dmask is cf in mm
C====== kmask is scale: 0 for constant, 1 for MMH
C====== all read in in match_nodals and converted to rasters
C====== via
C======     bminor = bminor/2.0*mm_rast
C======     cmask = cmask/2.0*mm_rast
C======     dmask = ctof*mm_rast
C
C
      PI = 4.0*ATAN(1.0)
C
      RADIUS = SQRT(FLOAT(IXM**2)+FLOAT(IYM**2))
      AMAJOR = CMASK*(1.+KMASK*RADIUS/DMASK)
      IF (IXM.EQ.0) THEN
         THETA = PI/2.0
         IF (IYM.LT.0) THEN
            THETA = 1.5*PI
         END IF
         GO TO 20
      END IF
C
      THETA = ATAN(FLOAT(IYM)/FLOAT(IXM))
      IF (IXM.LT.0) THEN
         THETA = THETA + PI
      ELSE IF (IYM.LT.0) THEN
         THETA = THETA + 2.0*PI
      END IF
 20   IMARK = 0
C
C====== ods are stored as if film scanned backwards (IDIR = -1)
C====== An old throwback from interactive work.
C====== RMASK must therefore be stored the same way.
C====== Printing of masks, spots, profiles etc by RASPLOT, PROSCL
C====== and here are now the same. Eventually need to use IDIR since
C====== not sure which way spot refinement works.
C
      IMARK = 0
      XPRIME = (FLOAT(XPN)-IXM)*COS(THETA)
      XPRIME = XPRIME + (FLOAT(YPN)-IYM)*SIN(THETA)
      YPRIME = (FLOAT(YPN)-IYM)*COS(THETA)
      YPRIME = YPRIME - (FLOAT(XPN)-IXM)*SIN(THETA)
      ABG = AMAJOR + 1.0 + BORDER
      BBG = BMINOR + 1.0 + BORDER
      ELLIPB = (XPRIME**2/ABG**2) + (YPRIME**2/BBG**2)
      IF (ELLIPB.LE.1.0) THEN
         AIG = AMAJOR + 1.0
         BIG = BMINOR + 1.0
         ELLIPI = (XPRIME**2/AIG**2) + (YPRIME**2/BIG**2)
         IF (ELLIPI.LE.1.0) THEN
            ELLIP = (XPRIME**2/AMAJOR**2) + (YPRIME**2/BMINOR**2)
            IF (ELLIP.LE.1.0) THEN
C
C====== PEAK return imark = 1 to indicate include in sums
C
               IMARK = 1
            END IF
         END IF
      END IF
C
      RETURN
      END
**==read_spots.f  
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       READ_SPOTS        **
C*****************************
C
C
      SUBROUTINE READ_SPOTS (NUM_CURS)
C
C Purpose : read spots in from a dump file
C
C
C Author : IJC
C
      INTEGER NUM_CURS
C
      INCLUDE 'observed_spots.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'gen_file_header.finc'
      INCLUDE 'sizes.finc'
C
C====== LOCALS:
C
      LOGICAL QUIT
      CHARACTER*80 STR
      INTEGER I, DUMP_FIL
      REAL X, Y
      PARAMETER (DUMP_FIL=15)
      CHARACTER*40 FILNAM
      CHARACTER*40 DUMP_FILE
      CHARACTER XTAL_ID*6, TTL*20
      LOGICAL CCPEXS
      LOGICAL EXISTS
C
C-----------------------------------------------------------------------
C
 20   CALL TPROMP('Dump file name (default ext=.spots): ')
      CALL TREAD(FILNAM,QUIT)
      CALL FILEXT(FILNAM,'spots',DUMP_FILE,I)
      EXISTS = CCPEXS(DUMP_FILE)
      IF (.NOT.EXISTS) THEN
         CALL TWRITE('!!! Cannot find file')
         GO TO 20
      END IF
      OPEN (FILE=DUMP_FILE,STATUS='OLD',UNIT=DUMP_FIL)
      READ (DUMP_FIL,6010) XTAL_ID , TTL
      IF (XTAL_ID.NE.CRYSTAL_ID) THEN
         WRITE (STR,6020) XTAL_ID
         CALL TWRITE(STR)
         WRITE (STR,6030) CRYSTAL_ID
         CALL TWRITE(STR)
      END IF
      READ (DUMP_FIL,6040) NUM_CURS
      NUM_CURS = MIN(NUM_CURS,MAX_OBS)
      DO 30 I = 1 , NUM_CURS
         READ (DUMP_FIL,6050) NODAL_H_S(I) , NODAL_K_S(I) , NODAL_L_S(I)
     +                        , X , Y
C
C====== Convert to scanner units, centre
C
         XM_OBS(I) = (X*MM_RAST_X) + X_CEN_F
         YM_OBS(I) = (Y*MM_RAST_Y) + Y_CEN_F
 30   CONTINUE
      CLOSE (DUMP_FIL)
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,A6/1X,A20)
 6020 FORMAT ('Warning - crystal id in dump file ',A6)
 6030 FORMAT ('doesn`t match the id in the generate file ',A6,'.')
 6040 FORMAT (8X,I4)
 6050 FORMAT (8X,3I4,T40,5X,2F6.1)
      END
**==readconsts.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*****************************
C**       READCONSTS        **
C*****************************
C
C
      SUBROUTINE  READCONSTS(IUNIT)
C
C Purpose: Reads `nearly-constant' parameters from the .GEN file
C          (on unit IUNIT) or alternatively `LAUEFLM(.DAT)'.
C          Don't try .GEN file if IUNIT<=0.
C          NB These params are also used by INTLAUE and so some aren't
C          needed in GENLAUE.
C
      INCLUDE 'implicit_none.finc'
C
      INCLUDE 'sizes.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'flmconsts.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'scnlim.finc'
C
      INTEGER IUNIT
C
C====== LOCALS:
C
      INTEGER IFAIL
      INTEGER I, XCEN, YCEN
      CHARACTER*80 DUMMY, STR
      REAL VALUES(10)
      INTEGER NVALS
      REAL RASTSIZ
      COMMON/IPLATE/ IPLATE
      LOGICAL IPLATE
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
C------------------------------------------------------------------------------
C
C====== Just try LAUEFLM.DAT if unit no. <= 0
C
      IF (IUNIT.GT.0) THEN
C
C====== First try and read from .GEN file - skip 11 records.
C
         DO 20 I = 1 , 11
            READ (9,6010) DUMMY
 20      CONTINUE
         READ (IUNIT,END=30,FMT=*) NUM_FID , 
     +                             (FID_POS_X10U(I),FID_POS_Y10U(I),I=1,
     +                             3)
         READ (IUNIT,*) FILM_WIDTH , FID_BOX_SIZE , N1OD , BASEOD , 
     +                  G1OD , CURVE
         READ (IUNIT,'(A)') DUMMY
         CALL STRVAL(DUMMY,' ,',10,VALUES,NVALS)
         X_CEN_F = 0.0
         Y_CEN_F = 0.0
         RASTSIZ = 0.0
         XLOW = 1.0
         YLOW = 1.0
         IF (NVALS.GE.1) X_RASTERS = NINT(VALUES(1))
         IF (NVALS.GE.2) Y_RASTERS = NINT(VALUES(2))
         XHIGH = X_RASTERS
         YHIGH = Y_RASTERS
         IF (NVALS.GE.3) FILMPLOT_SCALE = VALUES(3)
         IF (NVALS.GE.4) X_CEN_F = VALUES(4)
         IF (NVALS.GE.5) Y_CEN_F = VALUES(5)
         IF (IPLATE) THEN
            IF(X_CEN_F.EQ.0) X_CEN_F = VALUES(1)/2.0
            IF(Y_CEN_F.EQ.0) Y_CEN_F = VALUES(2)/2.0
         ENDIF
         IF (NVALS.GE.6) RASTSIZ = VALUES(6)
         IF(NVALS.GE.7) XLOW = VALUES(7)
         IF(NVALS.GE.8) XHIGH = VALUES(8)
         IF(NVALS.GE.9) YLOW = VALUES(9)
         IF(NVALS.GE.10) YHIGH = VALUES(10)
         GO TO 40
      END IF
C
C====== Data missing from .GEN file - try to read LAUEFLM(.DAT)
C
 30   WRITE (STR,6020)
      CALL TWRITE(STR)
      IFAIL = 0
      CALL CCPDPN(IN_DATA,'LAUEFLM','READONLY','F',0,IFAIL)
      READ (IN_DATA,*) NUM_FID , (FID_POS_X10U(I),FID_POS_Y10U(I),I=1,3)
      READ (IN_DATA,*) FILM_WIDTH , FID_BOX_SIZE , N1OD , BASEOD , 
     +                 G1OD , CURVE
      READ (IN_DATA,'(A)') DUMMY
      CALL STRVAL(DUMMY,' ,',6,VALUES,NVALS)
      X_CEN_F = 0.0
      Y_CEN_F = 0.0
      RASTSIZ = 0.0
      IF (NVALS.GE.1) X_RASTERS = NINT(VALUES(1))
      IF (NVALS.GE.2) Y_RASTERS = NINT(VALUES(2))
      IF (NVALS.GE.3) FILMPLOT_SCALE = VALUES(3)
      IF (NVALS.GE.4) X_CEN_F = VALUES(4)
      IF (NVALS.GE.5) Y_CEN_F = VALUES(5)
      IF (IPLATE) THEN
         IF(X_CEN_F.EQ.0) X_CEN_F = VALUES(1)/2.0
         IF(Y_CEN_F.EQ.0) Y_CEN_F = VALUES(2)/2.0
      ENDIF
      IF (NVALS.GE.6) RASTSIZ = VALUES(6)
      IF(NVALS.GE.7) XLOW = VALUES(7)
      IF(NVALS.GE.8) XHIGH = VALUES(8)
      IF(NVALS.GE.9) YLOW = VALUES(9)
      IF(NVALS.GE.10) YHIGH = VALUES(10)
      CLOSE (IN_DATA)
C
 40   IF (NUM_FID.NE.0 .AND. NUM_FID.NE.3) THEN
         CALL TWRITE('!!! Genlaue can only handle 3 fiducials')
         STOP
      END IF
      CENT_IN = .FALSE.
      IF (X_CEN_F.GT.0.0 .AND. Y_CEN_F.GT.0.0) CENT_IN = .TRUE.
C
C====== Default for film width is 125mm (standard film)
C
      RASTER_SIZE = RASTSIZ/1000.0
      IF (FILM_WIDTH.EQ.0.0) FILM_WIDTH = 125.0
      IF (RASTER_SIZE.LE.0) THEN
         RASTER_SIZE = FILM_WIDTH/X_RASTERS
C
C====== Set raster size to the nearest value that GENLAUE supports.
C
         CALL GET_RASTER_SIZE(RASTER_SIZE)
         MM_RAST_X = 1./RASTER_SIZE
         MM_RAST_Y = MM_RAST_X*Y_SCALE
      ELSE
         MM_RAST_X = 1./RASTER_SIZE
         MM_RAST_Y = MM_RAST_X*Y_SCALE
      END IF
C
C====== Approx centre of film
C
      XCEN = X_RASTERS/2
      YCEN = Y_RASTERS/2
      IF (NUM_FID.EQ.0) THEN
         XCEN = X_CEN_F
         YCEN = Y_CEN_F
      END IF
C
C====== Now convert fid posns from 10u units to rasters.
C
      DO 50 I = 1 , 3
         FID_POS_X(I) = (FID_POS_X10U(I)/100.*MM_RAST_X) + XCEN
         FID_POS_Y(I) = (FID_POS_Y10U(I)/100.*MM_RAST_Y) + YCEN
C
C====== Initialize these in case we don't look at a film
C
         FIDX(I) = FID_POS_X(I)
         FIDY(I) = FID_POS_Y(I)
 50   CONTINUE
      REWIND (IUNIT)
C
C====== Format statements
C
 6010 FORMAT (A80)
 6020 FORMAT ('Warning - data missing in .gen file'/
     +        '  -trying to read from LAUEFLM')
      END
**==refine_gen_parameters.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C****************************************
C**       REFINE_GEN_PARAMETERS        **
C****************************************
C
C
      SUBROUTINE  REFINE_GEN_PARAMETERS (WAVSTD, WAVMIN, WAVMAX,
     +                        FILM_RAD, NUM_REF_SPOTS, FILM_NUM)
C
C Purpose: To refine parameters for generate_reflexions
C
C Notes: the fiducial constant w_c is refined only on lower (B,C,...)
C        films. The 3 phis are refined on A film only.
C        The fiducial constants x_c, y_c are refined on all films-
C        although they are probably determined well by the position of
C        the main beam.
C        ctof, twist, tilt, bulge and y_scale are not refined.
C        Cell parameters may be optionally refined.
C
      INCLUDE 'implicit_none.finc'
C
      REAL       WAVSTD, WAVMIN, WAVMAX
      REAL       FILM_RAD
      INTEGER    NUM_REF_SPOTS, FILM_NUM
C
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'refined_film_params.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'reflexions_def.finc'
      INCLUDE 'reflexions.finc'
      INCLUDE 'observed_spots.finc'
C
      COMMON /CALFUN_PARAMS/ SCALING, XC_START, FIRST
      DOUBLE PRECISION      SCALING(20), XC_START(20)
      LOGICAL    FIRST
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
C
      REAL       OLD_RMS
      COMMON /RMS/ OLD_RMS
 
      COMMON /RFLG/ REFLAG(20), RFCL
      INTEGER       REFLAG
      LOGICAL       RFCL
      SAVE /RFLG/

      EXTERNAL CALFUN
C
C====== LOCALS
C
      CHARACTER*80 STR
      LOGICAL    RFN
      INTEGER    I, J, IFAIL, IFLAG, REFTYP, N, M
      INTEGER    H, K, L
      REAL       XM, YM, XM_O, YM_O
      INTEGER    MAX_M, MAX_N, LW1
      PARAMETER (MAX_N=13, MAX_M=MAX_OBS)
      PARAMETER (LW1=5*MAX_N+2*MAX_N*MAX_N+2*MAX_M*MAX_N+2*MAX_M)
      DOUBLE PRECISION     XC(20), FSUMSQ, RMS
      REAL       XY_C_PROMPT, TEMP
      DOUBLE PRECISION     W(LW1)
      INTEGER    MAXFUN
      DOUBLE PRECISION     DSTEP, ACC, FVECC(MAX_M), DAMP, SHFLIM
      DOUBLE PRECISION     FMIN
      LOGICAL    LOCSCH
      INTEGER    IPRINT, INTYPE, MAXCAL, IBOUND, LH, ISTATE(MAX_N)
      DOUBLE PRECISION ETA, XTOL, STEPMX, FEST, DELTA(MAX_N)
      DOUBLE PRECISION BL(MAX_N), BU(MAX_N), HESL(MAX_N*(MAX_N-1)/2)
      DOUBLE PRECISION HESD(MAX_N), F, G
      LOGICAL    QUIT, REFINE_CELL
      REAL   MM_TEMP
      CHARACTER*7  PROMPT(MAX_N)
      DATA  PROMPT /'phi x', 'phi y', 'phi z', 'c to f',
     1   'x_c', 'y_c', 'w_c',
     2   '1st', '2nd', '3rd', '4th', '5th', '6th' /
C
C------------------------------------------------------------------------------
C
      M = NUM_REF_SPOTS
C
      REFINE_CELL = .FALSE.
      CALL TASKQ('Do you want to refine cell parameters',REFINE_CELL,
     +           QUIT)
C
      RFCL = REFINE_CELL
      IF (REFINE_CELL) THEN
         CALL TLOG(.FALSE.)
         CALL TWRITE('Warning - you must leave at least 1 cell length')
         CALL TWRITE('fixed for proper convergence')
         CALL TLOG(.TRUE.)
      END IF
C
C====== We refine different subsets of the params on different films.
C
C
C====== Set up scaling constants
C
      DO 20 I = 1 , 3
         SCALING(I) = 1
 20   CONTINUE
C
C====== SCALING(4) was 1e2 at one stage
C
      SCALING(4) = 1
      SCALING(5) = 1
      SCALING(6) = 1
      SCALING(7) = 1
      DO 30 I = 8 , 13
         SCALING(I) = 1E-3
 30   CONTINUE
C
C====== Zero vector of offsets
C
      DO 40 I = 1 , 20
         XC(I) = 0.0
 40   CONTINUE
C
C====== Transfer parameters which may be refined into XC_START vector
C
      DO 50 I = 1 , 3
         XC_START(I) = PHIS(I)
 50   CONTINUE
      XC_START(4) = CTOF
      XC_START(5) = X_C
      XC_START(6) = Y_C
      XC_START(7) = W_C
C
      IF (REFINE_CELL) THEN
         DO 60 I = 8 , 7 + NUM_CELL_PARAMS
            XC_START(I) = CELL_PARAMS(I-7)
 60      CONTINUE
      END IF
 
      N = 0
C
C====== Get starting values from user - defaults are the current values
C
C
C====== On A film refine phis but not omega, on others refine w but not phis
C
      DO 70 I = 1 , 7
         REFLAG(I) = 0
         IF (FILM_NUM.EQ.1) THEN
            IF (I.EQ.7) GO TO 70
         ELSE IF (I.LT.4) THEN
            GO TO 70
         END IF
         IF (I.EQ.5 .OR. I.EQ.6) THEN
C
C====== Have to fiddle x_c, y_c into mm.
C
            MM_TEMP = MM_RAST_X
            IF (I.EQ.6) MM_TEMP = MM_RAST_Y
            XY_C_PROMPT = XC_START(I)/MM_TEMP
            CALL TASKF('Enter starting value for '//PROMPT(I),
     +                 XY_C_PROMPT,3,QUIT)
            XC_START(I) = XY_C_PROMPT*MM_TEMP
         ELSE
            TEMP = XC_START(I)
            CALL TASKF('Enter starting value for '//PROMPT(I),TEMP,3,
     +                 QUIT)
            XC_START(I) = TEMP
         END IF
         RFN = .TRUE.
         CALL TASKQ('Refine this parameter',RFN,QUIT)
         IF (RFN) THEN
            N = N + 1
            REFLAG(I) = N
         END IF
 70   CONTINUE
C
      IF (REFINE_CELL) THEN
         DO 80 I = 8 , 7 + NUM_CELL_PARAMS
            REFLAG(I) = 0
            TEMP = XC_START(I)
            CALL TASKF('Enter starting value for '//PROMPT(I)
     +                 //'cell parameter',TEMP,3,QUIT)
            XC_START(I) = TEMP
            RFN = .TRUE.
            IF (I.EQ.8) RFN = .FALSE.
            CALL TASKQ('Refine this parameter',RFN,QUIT)
            IF (RFN) THEN
               N = N + 1
               REFLAG(I) = N
            END IF
 80      CONTINUE
      END IF
      WRITE (STR,6080) N , M
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      CALL TWRITE('===Starting refinement...')
      FIRST = .TRUE.
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(M,N,XC,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,6010) 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(M,N,XC,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,6010) IFLAG
            CALL TWRITE(STR)
         END IF
      END IF
C
      RMS = FMIN/M
      RMS = DSQRT(RMS)
C
C====== Copy refined parameters back, scaling appropriately
C
      DO 90 I = 1 , 3
         IF (REFLAG(I).GT.0) THEN
            PHIS(I) = XC(REFLAG(I))*SCALING(REFLAG(I)) + XC_START(I)
         ELSE
            PHIS(I) = XC_START(I)
         END IF
 90   CONTINUE
      IF (REFLAG(4).GT.0) THEN
         CTOF = XC(REFLAG(4))*SCALING(REFLAG(4)) + XC_START(4)
      ELSE
         CTOF = XC_START(4)
      END IF
      IF (REFLAG(5).GT.0) THEN
         X_C = XC(REFLAG(5))*SCALING(REFLAG(5)) + XC_START(5)
      ELSE
         X_C = XC_START(5)
      END IF
      IF (REFLAG(6).GT.0) THEN
         Y_C = XC(REFLAG(6))*SCALING(REFLAG(6)) + XC_START(6)
      ELSE
         Y_C = XC_START(6)
      END IF
      IF (REFLAG(7).GT.0) THEN
         W_C = XC(REFLAG(7))*SCALING(REFLAG(7)) + XC_START(7)
      ELSE
         W_C = XC_START(7)
      END IF
C
      IF (REFINE_CELL) THEN
         DO 100 I = 8 , 7 + NUM_CELL_PARAMS
            IF (REFLAG(I).GT.0) THEN
               CELL_PARAMS(I-7) = XC(REFLAG(I))*SCALING(REFLAG(I))
     +                            + XC_START(I)
            ELSE
               CELL_PARAMS(I-7) = XC_START(I)
            END IF
 100     CONTINUE
      END IF
C
C====== Also remember film parameters for later output to .ge file
C
      CTOF_S(FILM_NUM) = CTOF
      X_C_S(FILM_NUM) = X_C
      Y_C_S(FILM_NUM) = Y_C
      W_C_S(FILM_NUM) = W_C
C
C====== Not refining these parameters
C
      TWISTS(FILM_NUM) = TWIST
      TILTS(FILM_NUM) = TILT
      BULGES(FILM_NUM) = BULGE
      Y_SCALES(FILM_NUM) = Y_SCALE
C
C====== Output params, converting x_c, y_c to mms.
C
      WRITE (STR,6020) (PHIS(I),I=1,3)
      CALL TWRITE(STR)
      WRITE (STR,6030) CTOF
      CALL TWRITE(STR)
      WRITE (STR,6040) (X_C/MM_RAST_X) , (Y_C/MM_RAST_Y) , W_C
      CALL TWRITE(STR)
      WRITE (STR,6050) (CELL_PARAMS(I),I=1,NUM_CELL_PARAMS)
      CALL TWRITE(STR)
      WRITE (STR,6060) RMS , OLD_RMS
      CALL TWRITE(STR)
      WRITE (STR,6070) M
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
C====== Save RMS for comparison with subsequent ref`mnt.
C
      OLD_RMS = RMS
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Refinement completed after',I4,' iterations')
 6020 FORMAT ('Phi''s (deg):',3G14.6)
 6030 FORMAT ('C to f (mm.):',G14.6)
 6040 FORMAT ('X_c, y_c, w_c:',2G14.6,' (mm.) ',G14.6,' (deg)')
 6050 FORMAT ('Cell parameters (A):',6G14.6)
 6060 FORMAT ('Rms deviation: ',F7.3,' mm. (from ',F7.3,')')
 6070 FORMAT ('...for a total of',I4,' spots')
 6080 FORMAT (' Refining ',I2,' parameters against',I4,' spots...')
      END
**==sysabs.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       SYSABS        **
C*************************
C
C
        SUBROUTINE SYSABS (IH,IK,IL,ABSNT)
C
      INCLUDE 'cell_params.finc'
C
        INTEGER H(3)
        LOGICAL ABSNT
C
C-------------------------------------------------------------------------------
C
      H(1) = IH
      H(2) = IK
      H(3) = IL
      IF (ISYS.EQ.4) THEN
 
         ABSNT = .TRUE.
         DO 20 I = 1 , 2
            IF (MOD(H(I)+H(I+1),2).NE.0) RETURN
 20      CONTINUE
         ABSNT = .FALSE.
         RETURN
      ELSE
C
C====== This version for a,b,c,r,i & f lattice absences
C
         ABSNT = .FALSE.
         J = 0
         DO 30 I = 1 , 3
            J = J + H(I)*KSYS(I)
 30      CONTINUE
         IF (MOD(J,ISYS).NE.0) ABSNT = .TRUE.
         RETURN
      END IF
      END
**==sysahk.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       SYSAHK        **
C*************************
C
C
        SUBROUTINE SYSAHK(IH,IK,ABSNT)
C
C   This variant of sysabs checks to see if a given h,k is absent,
C   whatever the l value
C
C
      INCLUDE 'cell_params.finc'
C
        LOGICAL ABSNT
C
C-------------------------------------------------------------------------------
C
      ABSNT = .FALSE.
      IF ((ISYS.EQ.2 .OR. ISYS.EQ.4) .AND. KSYS(3).EQ.0) THEN
         IF (MOD(IH+IK,2).NE.0) ABSNT = .TRUE.
      END IF
      RETURN
      END
**==write_ge_file.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C********************************
C**       WRITE_GE_FILE        **
C********************************
C
C
      SUBROUTINE  WRITE_GE_FILE (GE_NAME,
     +      NUM_REFLEX, H_S, K_S, L_S, M_FLAGS, LAMBDAS, XF_S, YF_S,
     +      DMIN_THRESHOLDS2,
     +      TITLE, CRYSTAL_ID,
     +      CELL_PARAMS, NUM_CELL_PARAMS,
     +      XTAL_SYSTEM, LATTICE, X_AXIS, Z_AXIS,
     +      PHIS, WAVSTD, WAVMIN, WAVMAX, DMIN, FILM_RAD, SPINDLE,
     +      DELTA, EPSILON, MININT, RASTER,
     +      X_CEN_FS, Y_CEN_FS, W_FS,
     +      X_C_S, Y_C_S, W_C_S, TWISTS, TILTS, BULGES,
     +      Y_SCALES, CTOF_S,
     +      LAMBDA_BINS, X_REJECT, Y_REJECT)
C
C Purpose: To write out a skeleton .ge1/2 file
C
      INCLUDE 'implicit_none.finc'
C
      CHARACTER*(*) GE_NAME
      INTEGER    NUM_REFLEX
      INCLUDE 'reflexions_def.finc'
      INTEGER    M_FLAGS(*)
      REAL       LAMBDAS(NUM_REFLEX)
      REAL       XF_S(NUM_REFLEX), YF_S(NUM_REFLEX)
      REAL       DMIN_THRESHOLDS2(NUM_REFLEX)
      CHARACTER*76  TITLE
      CHARACTER*6   CRYSTAL_ID
      REAL          CELL_PARAMS(6)
      INTEGER       NUM_CELL_PARAMS, XTAL_SYSTEM, LATTICE,
     +              X_AXIS, Z_AXIS
      REAL          PHIS(3), WAVSTD, WAVMIN, WAVMAX
      REAL          DMIN, FILM_RAD, SPINDLE
      REAL          DELTA, EPSILON
      INTEGER       MININT
      REAL          RASTER
      REAL          X_CEN_FS(6), Y_CEN_FS(6), W_FS(6)
      REAL          X_C_S(6), Y_C_S(6), W_C_S(6)
      REAL          TWISTS(6), TILTS(6), BULGES(6)
      REAL          Y_SCALES(6), CTOF_S(6)
      INTEGER       LAMBDA_BINS(16)
      REAL          X_REJECT, Y_REJECT
C
      COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
      INTEGER             IN_TERM, OUT_TERM, IN_DATA, OUT_DATA
 
      INTEGER NBITST
C
C====== LOCALS:
C
      CHARACTER*80 STR
      CHARACTER*40  FNAME_GE1, FNAME_GE2
      LOGICAL    QUIT, WR_SPAT, WR_CLOS, WR_BACK
      INTEGER    NODAL_INDEX
      INTEGER    CURRENT_REFLEX, I, J, PACK_ID
      INTEGER    LINK_REC
      INTEGER    BINTENS(6), BSIGS(6), PINTENS(6), PSIGS(6)
      REAL       BX, BY, BDMIN2
      REAL       LAMBDA1, LAMBDA2
      INTEGER    M_FLAG, H, K, L, NODIDX, H2, K2, L2
      INTEGER    IREAL
      REAL       RINT
      INTEGER    IH,IK,IL
      INTEGER    NUM_OVERLAPS
      INTEGER    GE1_REC, GE2_REC
      INTEGER    MULTIPLICITY, MULTI, IX, IY, HKLPAK
      INTEGER    IPACK, NFILM
      INTEGER    MULTI_BINS(21), NUM_NODALS
      INTEGER    NUM_SPOTS, NUM_SPATIALS, NUM_MULTIPLETS, NUM_BOTH
      INTEGER    BUFFER(19), BUFFER2(2), MM
      INTEGER    MHKL(4), HKL2(4)
      INTEGER    NODAL_BINS(12)
      EQUIVALENCE (IREAL, RINT)
      EQUIVALENCE (BUFFER(1),BX)
      EQUIVALENCE (BUFFER(2),BY)
      EQUIVALENCE (BUFFER(16),LAMBDA1)
      EQUIVALENCE (BUFFER(18),LAMBDA2)
      EQUIVALENCE (BUFFER(19),BDMIN2)
C
      EQUIVALENCE (MHKL(1),M_FLAG)
      EQUIVALENCE (MHKL(2),H)
      EQUIVALENCE (MHKL(3),K)
      EQUIVALENCE (MHKL(4),L)
C
      EQUIVALENCE (HKL2(1),NODIDX)
      EQUIVALENCE (HKL2(2),H2)
      EQUIVALENCE (HKL2(3),K2)
      EQUIVALENCE (HKL2(4),L2)
C
      EQUIVALENCE (BUFFER2(2),LAMBDA1)
C
      INCLUDE 'novmax_def.finc'
      INTEGER    OVERLAP_SORT(NOVMAX)
      REAL       OVERLAP_X(NOVMAX)
      REAL       OVERLAP_Y(NOVMAX)
      INTEGER    OVERLAP_INDEX(NOVMAX)
      INTEGER    OVERLAP_HKL(NOVMAX)
      REAL       OVERLAP_LAMBDA(NOVMAX)
      REAL       OVERLAP_DMIN2(NOVMAX)
      INTEGER    IFAIL, NSING, NSINGOV
      DATA  BINTENS, BSIGS, PINTENS, PSIGS /6*-9999, 6*-9999,
     1                6*-9999, 6*-9999/
C
C------------------------------------------------------------------------------
C
      FNAME_GE1 = GE_NAME//'.ge1'
      FNAME_GE2 = GE_NAME//'.ge2'
      CALL CCPII2(BUFFER(3),BINTENS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(6),BSIGS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(9),PINTENS,6,.TRUE.,.FALSE.)
      CALL CCPII2(BUFFER(12),PSIGS,6,.TRUE.,.FALSE.)
      WR_SPAT = .FALSE.
      CALL TASKQ('Do you want to write spatials to the file',WR_SPAT,
     +           QUIT)
      WR_CLOS = .FALSE.
      WR_BACK = .FALSE.
      IF (X_REJECT.GT.0 .OR. Y_REJECT.GT.0) THEN
         CALL TASKQ('Do you want to write central spots to the file',
     +              WR_BACK,QUIT)
      ELSE
         WR_BACK = .TRUE.
      END IF
C
      IFAIL = 1
      CALL CCPDPN(2,FNAME_GE1,'NEW','DU',19,IFAIL)
      IFAIL = 1
      CALL CCPDPN(3,FNAME_GE2,'NEW','DU',2,IFAIL)
      IF (IFAIL.LT.0) THEN
         CLOSE (UNIT=2)
         CLOSE (UNIT=3)
      END IF
C
C====== Zero histograms
C
      NUM_OVERLAPS = 0
      NUM_SPOTS = 0
      NUM_SPATIALS = 0
      NUM_MULTIPLETS = 0
      NUM_BOTH = 0
      NUM_NODALS = 0
      DO 20 I = 1 , 21
         MULTI_BINS(I) = 0
 20   CONTINUE
      NSING = 0
      NSINGOV = 0
      DO 30 I = 1 , 12
         NODAL_BINS(I) = 0
 30   CONTINUE
C
C====== Now write to units 2 and 3
C
      GE1_REC = 10
      GE2_REC = 0
C
      DO 40 CURRENT_REFLEX = 1 , NUM_REFLEX
C
         BX = XF_S(CURRENT_REFLEX)
         BY = YF_S(CURRENT_REFLEX)
C
C====== Check central box
C
         IF (WR_BACK .OR. (ABS(BX).GT.X_REJECT .OR. ABS(BY).GT.Y_REJECT)
     +       ) THEN
C
 
            M_FLAG = M_FLAGS(CURRENT_REFLEX)
            MULTIPLICITY = NBITST(M_FLAG,0,5)
            IF ((NBITST(M_FLAG,6,1).EQ.0 .OR. WR_SPAT) .AND. 
     +          (NBITST(M_FLAG,7,1).EQ.0 .OR. WR_CLOS)) THEN
C
C====== Get hkl,lambda for at least the singlet
C
               H = H_S(CURRENT_REFLEX)
               K = K_S(CURRENT_REFLEX)
               L = L_S(CURRENT_REFLEX)
C
               LAMBDA1 = LAMBDAS(CURRENT_REFLEX)
C
               BDMIN2 = DMIN_THRESHOLDS2(CURRENT_REFLEX)
C
               IF (MULTIPLICITY.EQ.1) THEN
C
C====== Now write out .ge1 record for the singlet
C
                  GE1_REC = GE1_REC + 1
                  IF (NBITST(M_FLAG,5,1).NE.0) THEN
                     NUM_NODALS = NUM_NODALS + 1
                     IH = H
                     IK = K
                     IL = L
                     NODIDX = NODAL_INDEX(IH,IK,IL)
                     NODAL_BINS(NODIDX) = NODAL_BINS(NODIDX) + 1
                  ELSE
                     NODIDX = 0
                  END IF
C
                  CALL CCPIBY(BUFFER(15),MHKL,-4)
                  CALL CCPIBY(BUFFER(17),HKL2,-4)
                  WRITE (2,REC=GE1_REC) BUFFER
C
                  MULTI_BINS(1) = MULTI_BINS(1) + 1
                  NSING = NSING + 1
                  IF (NBITST(M_FLAG,6,1).NE.0)
     +                NUM_SPATIALS = NUM_SPATIALS + 1
C
               ELSE
C
                  IF (MULTIPLICITY.NE.0) STOP 'write_ge_file: multi<>0'
C
C====== This is a multiplet so add to overlap list
C
                  NUM_OVERLAPS = NUM_OVERLAPS + 1
                  IF (NUM_OVERLAPS.GT.NOVMAX) THEN
                     CALL TWRITE('!!! Too many overlaps - stop')
                     STOP
                  END IF
C
                  IH = H
                  IK = K
                  IL = L
                  CALL FACTOR_INDEX(IH,IK,IL)
                  OVERLAP_SORT(NUM_OVERLAPS) = 40000*(IH+100)
     +               + 200*(IK+100) + IL + 100
                  OVERLAP_X(NUM_OVERLAPS) = BX
                  OVERLAP_Y(NUM_OVERLAPS) = BY
                  CALL CCPIBY(OVERLAP_HKL(NUM_OVERLAPS),MHKL,-4)
                  OVERLAP_LAMBDA(NUM_OVERLAPS) = LAMBDA1
                  OVERLAP_DMIN2(NUM_OVERLAPS) = BDMIN2
C
               END IF
            END IF
         END IF
C
 40   CONTINUE
C
C====== Now sort out overlap list
C
      IF (NUM_OVERLAPS.EQ.0) GO TO 90
      CALL INDXIV(OVERLAP_SORT,NUM_OVERLAPS,OVERLAP_INDEX)
C
      J = 1
C
C====== Determine multiplicity by finding out how many
C====== consecutive reflexions have the same factored-out indices.
C
 50   MULTI = 1
      HKLPAK = OVERLAP_SORT(OVERLAP_INDEX(J))
C
 60   IF (OVERLAP_SORT(OVERLAP_INDEX(J+MULTI)).EQ.HKLPAK .AND. 
     +    J+MULTI.LE.NUM_OVERLAPS) THEN
         MULTI = MULTI + 1
         IF (MULTI.LE.30) GO TO 60
         WRITE (STR,6010) MULTI
         CALL TWRITE(STR)
      END IF
C
      BX = OVERLAP_X(OVERLAP_INDEX(J))
      BY = OVERLAP_Y(OVERLAP_INDEX(J))
C
      CALL CCPBYI(MHKL,OVERLAP_HKL(OVERLAP_INDEX(J)),-4)
      M_FLAG = M_FLAG + MULTI
      LAMBDA1 = OVERLAP_LAMBDA(OVERLAP_INDEX(J))
C
C====== We need lowest Dmin2 value calculated
C====== this corresponds to transition from singlet to doublet in a doublet
C====== or from not being on film to being present for a singlet
C
      BDMIN2 = 1E6
      DO 70 I = 1 , MULTI
         IF (OVERLAP_DMIN2(OVERLAP_INDEX(J-1+I)).LT.BDMIN2)
     +       BDMIN2 = OVERLAP_DMIN2(OVERLAP_INDEX(J-1+I))
 70   CONTINUE
C
      CALL CCPBYI(HKL2,OVERLAP_HKL(OVERLAP_INDEX(J+1)),-4)
      LAMBDA2 = OVERLAP_LAMBDA(OVERLAP_INDEX(J+1))
C
C====== Update histogram
C
      I = MIN(MULTI,21)
      MULTI_BINS(I) = MULTI_BINS(I) + 1
      IF (MULTI.EQ.1) NSINGOV = NSINGOV + 1
      IF (MULTI.GT.1) NUM_MULTIPLETS = NUM_MULTIPLETS + 1
      IF (NBITST(M_FLAG,5,1).NE.0) THEN
         NUM_NODALS = NUM_NODALS + 1
         IH = H
         IK = K
         IL = L
         NODIDX = NODAL_INDEX(IH,IK,IL)
         NODAL_BINS(NODIDX) = NODAL_BINS(NODIDX) + 1
      ELSE
         NODIDX = 0
      END IF
      IF (NBITST(M_FLAG,6,1).NE.0) THEN
         NUM_SPATIALS = NUM_SPATIALS + 1
         IF (MULTI.GT.1) NUM_BOTH = NUM_BOTH + 1
      END IF
C
C====== Now write out .ge1 record for the multiplet
C
      GE1_REC = GE1_REC + 1
      CALL CCPIBY(BUFFER(15),MHKL,-4)
      CALL CCPIBY(BUFFER(17),HKL2,-4)
      IF (MULTI.LE.2) THEN
C
         WRITE (2,REC=GE1_REC) BUFFER
C
      ELSE
C
         LINK_REC = GE2_REC + 1
         IREAL = LINK_REC
         LAMBDA2 = RINT
         WRITE (2,REC=GE1_REC) BUFFER
C
C====== Write out consecutive records to .ge2 file
C
         DO 80 I = 2 , MULTI
            CALL CCPBYI(MHKL,OVERLAP_HKL(OVERLAP_INDEX(J+I-1)),-4)
            LAMBDA1 = OVERLAP_LAMBDA(OVERLAP_INDEX(J+I-1))
            GE2_REC = GE2_REC + 1
            CALL CCPIBY(BUFFER2(1),MHKL,-4)
            WRITE (3,REC=GE2_REC) BUFFER2
 80      CONTINUE
C
      END IF
C
C====== More multiplets ?
C
      J = J + MULTI
      IF (J.LE.NUM_OVERLAPS) GO TO 50
C
C====== Finally write out header info, number of spots output
C
 90   IPACK = 0
      NFILM = 0
      NUM_SPOTS = GE1_REC - 10
      WRITE (2,REC=9) IPACK , NFILM , NUM_SPOTS
C
C======= Display histogram results
C
      CALL TWRITE(' ')
      WRITE (STR,6020)
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      WRITE (STR,6030)
      CALL TWRITE(STR)
      CALL TWRITE(' ')
      DO 100 I = 1 , 10
         WRITE (STR,6040) I , MULTI_BINS(I) , I + 10 , MULTI_BINS(I+10)
         CALL TWRITE(STR)
 100  CONTINUE
      WRITE (STR,6050) MULTI_BINS(21)
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
      WRITE (STR,6060) NUM_SPOTS
      CALL TWRITE(STR)
      WRITE (STR,6070) NUM_MULTIPLETS
      CALL TWRITE(STR)
      WRITE (STR,6080) NUM_SPATIALS
      CALL TWRITE(STR)
      WRITE (STR,6090) NUM_BOTH
      CALL TWRITE(STR)
      WRITE (STR,6100) NUM_NODALS
      CALL TWRITE(STR)
      CALL TWRITE(' ')
C
C====== Calc cumulative totals
C
      DO 110 I = 1 , 11
         NODAL_BINS(I+1) = NODAL_BINS(I+1) + NODAL_BINS(I)
 110  CONTINUE
C
C====== Now write out header!
C
C
      CALL TASKI('Pack ID',PACK_ID,QUIT)
C
      CALL WRITE_GE_HEADER(2,TITLE,CRYSTAL_ID,CELL_PARAMS,
     +                     NUM_CELL_PARAMS,XTAL_SYSTEM,LATTICE,X_AXIS,
     +                     Z_AXIS,PHIS,WAVSTD,WAVMIN,WAVMAX,DMIN,
     +                     FILM_RAD,SPINDLE,DELTA,EPSILON,MININT,RASTER,
     +                     X_CEN_FS,Y_CEN_FS,W_FS,X_C_S,Y_C_S,W_C_S,
     +                     TWISTS,TILTS,BULGES,Y_SCALES,CTOF_S,
     +                     MULTI_BINS,LAMBDA_BINS,NODAL_BINS,PACK_ID,
     +                     NUM_SPOTS,NUM_MULTIPLETS,NUM_SPATIALS,
     +                     NUM_BOTH,NUM_NODALS)
C
      CLOSE (2)
      CLOSE (3)
      RETURN
C
C====== Format statements
C
 6010 FORMAT ('Large multiplicity in writegefile ',I4)
 6020 FORMAT ('              Histogram of output spots')
 6030 FORMAT ('Multiplicity  -  Frequency',13X,
     +        'Multiplicity  -  Frequency')
 6040 FORMAT (I8,I17,13X,I8,I17)
 6050 FORMAT ('     >20',I17)
 6060 FORMAT ('Summary:',I6,' spots output, of which')
 6070 FORMAT (8X,I6,' were multiplets,')
 6080 FORMAT (8X,I6,' were spatial overlaps,')
 6090 FORMAT (8X,I6,' were both, and')
 6100 FORMAT (8X,I6,' were nodals.')
      END
**==write_ge_header.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**********************************
C**       WRITE_GE_HEADER        **
C**********************************
C
C
      SUBROUTINE  WRITE_GE_HEADER (GE1_UNIT,
     +            TITLE, CRYSTAL_ID,
     +            CELL_PARAMS, NUM_CELL_PARAMS,
     +            XTAL_SYSTEM, LATTICE, X_AXIS, Z_AXIS,
     +            PHIS, WAVSTD, WAVMIN, WAVMAX, DMIN, FILM_RAD, SPINDLE,
     +            DELTA, EPSILON, MININT, RASTER,
     +            X_CEN_FS, Y_CEN_FS, W_FS,
     +            X_C_S, Y_C_S, W_C_S, TWISTS, TILTS, BULGES,
     +            Y_SCALES, CTOF_S,
     +            MULTI_BINS, LAMBDA_BINS, NODAL_BINS,
     +            PACK_ID, NSPOTS, NMULTIPLETS,
     +            NSPATOVS, NBOTH, NNODALS)
C
C Purpose: To write out a COMPREHENSIVE header for a .ge1/2 file
C
      INCLUDE 'implicit_none.finc'
C
      INTEGER       GE1_UNIT
      CHARACTER*76  TITLE
      CHARACTER*6   CRYSTAL_ID
      REAL          CELL_PARAMS(6)
      INTEGER       NUM_CELL_PARAMS, XTAL_SYSTEM, LATTICE,
     +              X_AXIS, Z_AXIS
      REAL          PHIS(3), WAVSTD, WAVMIN, WAVMAX
      REAL          DMIN, FILM_RAD, SPINDLE
      REAL          DELTA, EPSILON
      INTEGER       MININT
      REAL          RASTER
      REAL          X_CEN_FS(6), Y_CEN_FS(6), W_FS(6)
      REAL          X_C_S(6), Y_C_S(6), W_C_S(6)
      REAL          TWISTS(6), TILTS(6), BULGES(6)
      REAL          Y_SCALES(6), CTOF_S(6)
      INTEGER       MULTI_BINS(21), LAMBDA_BINS(16), NODAL_BINS(12)
      INTEGER       PACK_ID, NSPOTS,NMULTIPLETS,NSPATOVS,NBOTH,NNODALS
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    EOF
      INTEGER    I
      INTEGER    BUFFER(19)
      REAL       X_C_MM(6), Y_C_MM(6)
      INTEGER    BOX(5)
      REAL       CTOF_RATS(6)
      LOGICAL    USE22, USEREF
      INTEGER    ISYLAT(2), ISL
      REAL       THICK
      INTEGER    QUIT
      DATA       BUFFER/19*0/
C
C------------------------------------------------------------------------------
C
      WRITE (GE1_UNIT,REC=1) TITLE
      WRITE (GE1_UNIT,REC=2) CRYSTAL_ID(1:6) , X_CEN_FS , Y_CEN_FS , 
     +                       (W_FS(I),I=1,5)
C
      ISYLAT(1) = XTAL_SYSTEM
      ISYLAT(2) = LATTICE
      CALL CCPII2(ISL,ISYLAT,2,.TRUE.,.FALSE.)
      WRITE (GE1_UNIT,REC=3) W_FS(6) , (BUFFER(I),I=1,9) , CELL_PARAMS , 
     +                       ISL , X_AXIS , Z_AXIS
C
C====== Initialise integ. box params
C
      BOX(1) = DELTA/RASTER + 1
C
C====== Force box size to be odd.
C
      IF (BOX(1)*2/2.EQ.BOX(1)) BOX(1) = BOX(1) + 1
      BOX(2) = BOX(1)
C
C====== Cut off a corner 1/4 way along.
C
      BOX(3) = BOX(1)/4
      BOX(4) = 1
      BOX(5) = BOX(4)
      WRITE (GE1_UNIT,REC=4) RASTER , BOX , (BUFFER(I),I=1,4) , PHIS , 
     +                       WAVSTD , WAVMIN , WAVMAX , DMIN , 
     +                       FILM_RAD , SPINDLE
C
C====== Convert centres to mm units from scanner units
C
      DO 20 I = 1 , 6
         X_C_MM(I) = X_C_S(I)*RASTER
         Y_C_MM(I) = Y_C_S(I)*RASTER
 20   CONTINUE
      WRITE (GE1_UNIT,REC=5) X_C_MM , Y_C_MM , W_C_S , DELTA
C
C====== Calculate cryst-film ratios either from user`s given film thickness
C====== or from refined c-to-f's.
C
      CTOF_RATS(1) = 1.
      USEREF = .FALSE.
      USE22 = .TRUE.
      CALL TASKQ('Use individual refined crystal-film dists',USEREF,EOF)
      IF (.NOT.USEREF) CALL TASKQ(
     +                     'Use a film thickness of .22mm for all films'
     +                     ,USE22,EOF)
      IF (USEREF) THEN
         DO 30 I = 2 , 6
            CTOF_RATS(I) = CTOF_S(I)/CTOF_S(1)
 30      CONTINUE
      ELSE IF (USE22) THEN
         DO 40 I = 2 , 6
            CTOF_RATS(I) = CTOF_RATS(1) + (I-1)*.22/CTOF_S(1)
 40      CONTINUE
      ELSE
         THICK = .22
         DO 50 I = 2 , 6
            CALL TASKF('Enter each interfilm thickness(mm)',THICK,2,EOF)
            CTOF_RATS(I) = CTOF_RATS(I-1) + THICK/CTOF_S(1)
 50      CONTINUE
      END IF
      WRITE (GE1_UNIT,REC=6) TWISTS , TILTS , BULGES , CTOF_S(1)
      WRITE (GE1_UNIT,REC=7) Y_SCALES , CTOF_RATS , MININT , EPSILON
C
      WRITE (GE1_UNIT,REC=8) (MULTI_BINS(I),I=1,19)
      WRITE (GE1_UNIT,REC=9) (MULTI_BINS(I),I=20,21) , 
     +                       (LAMBDA_BINS(I),I=1,16) , 2
C
      WRITE (GE1_UNIT,REC=10) PACK_ID , NSPOTS , NMULTIPLETS , 
     +                        NSPATOVS , NBOTH , NNODALS , NODAL_BINS
C
      RETURN
      END
**==write_gen_parameters.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************************
C**       WRITE_GEN_PARAMETERS       **
C**************************************
C
C
      SUBROUTINE WRITE_GEN_PARAMETERS (GE_NAME, WAVSTD, WAVMIN,
     +                  WAVMAX, DMIN, FILM_RAD, X_REJECT, Y_REJECT)
C
C Purpose: To write out refined gen parameters to a .gen file
C
      INCLUDE 'implicit_none.finc'
C
      CHARACTER*(*) GE_NAME
      REAL    WAVSTD, WAVMIN, WAVMAX
      REAL    DMIN, FILM_RAD
      REAL    X_REJECT, Y_REJECT
C
      INCLUDE 'cell_params.finc'
      INCLUDE 'crystal_params.finc'
      INCLUDE 'film_params.finc'
      INCLUDE 'flmconsts.finc'
      INCLUDE 'fiducials.finc'
      INCLUDE 'refined_film_params.finc'
      INCLUDE 'sizes.finc'
      INCLUDE 'gen_file_header.finc'
      INCLUDE 'scnlim.finc'
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 EOF
      CHARACTER*80 STR
      CHARACTER*40 FNAME, FTEMP
      INTEGER      I, J, IFAIL
      REAL         FIDX10U(3), FIDY10U(3), XCEN, YCEN, XXC, YYC
C
C------------------------------------------------------------------------------
C
      FNAME = GE_NAME//'.gen'
 20   IFAIL = 1
      CALL CCPDPN(9,FNAME,'NEW','F',0,IFAIL)
      IF (IFAIL.LT.0) THEN
         WRITE (STR,6140)
         CALL TPROMP(STR)
         CALL TREAD(FTEMP,EOF)
         CALL FILEXT(FTEMP,'gen',FNAME,I)
         GO TO 20
      END IF
      WRITE (9,6010) TITLE
      WRITE (9,6020) CRYSTAL_ID
C
C====== UNIT_CELL_INFO already contains info - don't overwrite it
C
      IF (XTAL_SYSTEM.EQ.1) THEN
         WRITE (UNIT_CELL_INFO(1:60),6030) (CELL_PARAMS(I),I=1,6)
      ELSE IF (XTAL_SYSTEM.EQ.2) THEN
         WRITE (UNIT_CELL_INFO(1:60),6040) (CELL_PARAMS(I),I=1,4)
      ELSE IF (XTAL_SYSTEM.EQ.3) THEN
         WRITE (UNIT_CELL_INFO(1:60),6050) (CELL_PARAMS(I),I=1,3)
      ELSE IF (XTAL_SYSTEM.EQ.4) THEN
         WRITE (UNIT_CELL_INFO(1:60),6060) (CELL_PARAMS(I),I=1,2)
      ELSE IF (XTAL_SYSTEM.EQ.5) THEN
         WRITE (UNIT_CELL_INFO(1:60),6070) (CELL_PARAMS(I),I=1,2)
      ELSE IF (XTAL_SYSTEM.EQ.6) THEN
         WRITE (UNIT_CELL_INFO(1:60),6080) (CELL_PARAMS(I),I=1,2)
      ELSE IF (XTAL_SYSTEM.EQ.7) THEN
         WRITE (UNIT_CELL_INFO(1:60),6090) CELL_PARAMS(1)
      ELSE
         CALL TWRITE('!!! Illegal crystal system code')
         STOP
      END IF
      WRITE (9,6100) UNIT_CELL_INFO
      WRITE (9,*) PHIS
      WRITE (9,*) WAVSTD , WAVMIN , WAVMAX
      WRITE (9,*) CTOF_S(1) , DMIN , FILM_RAD
      WRITE (9,*) SPINDLE
C
C====== Convert scanner units to mm.s
C
      WRITE (9,*) X_C_S(1)/MM_RAST_X , Y_C_S(1)/MM_RAST_Y , W_C_S(1)
      WRITE (9,*) TWISTS(1) , TILTS(1) , BULGES(1)
      WRITE (9,*) Y_SCALES(1)
      WRITE (9,*) X_REJECT , Y_REJECT
C
C====== Convert fid posns to original units, uncentred.
C
      XCEN = X_RASTERS/2
      YCEN = Y_RASTERS/2
      DO 30 I = 1 , 3
         FIDX10U(I) = (FIDX(I)-XCEN)/MM_RAST_X*100.
         FIDY10U(I) = (FIDY(I)-YCEN)/MM_RAST_Y*100.
 30   CONTINUE
      WRITE (9,6110) NUM_FID , (FIDX10U(I),FIDY10U(I),I=1,3)
      WRITE (9,6120) FILM_WIDTH , FID_BOX_SIZE , N1OD , BASEOD , G1OD , 
     +               CURVE
      IF (CENT_IN) THEN
         XXC = X_CEN_F
         YYC = Y_CEN_F
      ELSE
         XXC = 0.0
         YYC = 0.0
      END IF
      WRITE (9,6130) X_RASTERS , Y_RASTERS , FILMPLOT_SCALE , XXC , 
     +               YYC , RASTER_SIZE*1000.0,
     +               XLOW,XHIGH,YLOW,YHIGH
      CLOSE (9)
C
      RETURN
C
C====== Format statements
C
 6010 FORMAT (1X,A76)
 6020 FORMAT (A6)
 6030 FORMAT (3F10.7,3F10.5)
 6040 FORMAT (3F10.7,F10.5)
 6050 FORMAT (3F10.7)
 6060 FORMAT (2F10.7)
 6070 FORMAT (2F10.7)
 6080 FORMAT (F10.7,F10.5)
 6090 FORMAT (F10.7)
 6100 FORMAT (1X,A80)
 6110 FORMAT (I2,3(2F8.1))
 6120 FORMAT (F6.1,F4.1,I4,3E10.3)
 6130 FORMAT (2I5,F4.1,3F9.3,4F7.1)
 6140 FORMAT ('Root name for .ge1/.ge2/.gen files: ')
      END
