C
C PROGRAM 'FILMDSEQ'
C ==================
C
C PROGRAM TO FIND SPOTS ON A FILM IMAGE
C
C
C FILMDSEQ HISTORY:
C
C     REWORKED BY J.W.CAMPBELL & P.D.CARR, 1989
C
C     INCLUDES DISPLAY OF A FILM IMAGE 
C     7/6/84 DEVELOPED FOR MMH SMALL SPOTS
C
C     CALCULATE A RADIAL BACKGROUND CURVE
C     AND USE A THRESHOLD ABOVE THIS FOR DISPLAY
C
C     MODIFIED FOR SEQUENTIAL READ 28/6/84
C
C     SPOTS OPTIONALLY SAVED IN BUFFER AND OUTPUT TO FILE
C     WITH MAX INTENSITY , NO. OF VALUES, APPROX INTEG VALUES
C
c     image plate and qscreen graphics 11/12/91 pdc
C
C Parameter for maximum no. of spots which can be stored
      PARAMETER (MAXSPT=20000)
C
      common /fiducials/ x_cen_f, y_cen_f, w_f, x_c, y_c, w_c,
     1                   fidx, fidy, fid_pos_x, fid_pos_y
      real               x_cen_f, y_cen_f, w_f, x_c, y_c, w_c
C     Units:             rasters  rasters  rasters       rasters
      real               fidx(3), fidy(3), fid_pos_x(3), fid_pos_y(3)
      integer max_rasters
      parameter (max_rasters=9600)
C Size in pixels of scanned area of film.
      integer x_rasters, y_rasters
C Multiply by mm_rast to convert mm to pixels.
C   Units: mm           pix/mm   mm
      real raster_size, mm_rast, fid_box_size
      common /sizes/ raster_size, x_rasters, y_rasters, mm_rast,
     1 fid_box_size
      common /ip/ iplate,nchar_raster
      integer iplate,nchar_raster
      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
      COMMON /CHANEL/ ITERM ,IPRINT,IHELP ,IFILEA,IFILEB,ISPARE,
     1                IOUT  ,IRUNFL,ILRNFL,ISTATF,IUNIT8,INTEST
      INTEGER         ITERM ,IPRINT,IHELP ,IFILEA,IFILEB,ISPARE,
     1                IOUT  ,IRUNFL,ILRNFL,ISTATF,IUNIT8,INTEST
C
      COMMON /PGSTAT/ NOVICE,ABORT ,ALARM ,OVERID,CONTIN,PASWRD,
     1                ADUMY1,ADUMY2,ADUMY3,ADUMY4,ADUMY5,ADUMY6
      LOGICAL         NOVICE,ABORT ,ALARM ,OVERID,CONTIN,PASWRD,
     1                ADUMY1,ADUMY2,ADUMY3,ADUMY4,ADUMY5,ADUMY6
C
      CHARACTER*80 FNAME,FILOUT,STRING
      character*150 longstring
      COMMON/CDEN/IBB(max_rasters)
      COMMON/XCONST/NX,NY,ICX,ICY,NREC,NYREC,ISTART,IYST
      COMMON/CBACK/NBG,IBGR(30),IBGY(30),IBGV(30)
      COMMON/CSPOT/NSPOT,SPOT(2,MAXSPT),ISPOT(5,MAXSPT),IFLAG
      integer*4 itmp(max_rasters/2)
      LOGICAL BAD
      LOGICAL VEC
      INTEGER film_unit
      CHARACTER*1 IYN
      CHARACTER*2 ANS
      REAL FIDXC1(3),FIDYC1(3),FIDXC2(3),FIDYC2(3)
      DATA FIDXC1/-5000,-5000,5000/
      DATA FIDYC1/-4000,4000,4000/
      DATA FIDXC2/-4000,-4000,4000/
      DATA FIDYC2/5000,-5000,-5000/
C
C=====Initialisations
C
      in_term=5
      out_term=6
      ITERM=5
      IPRINT=6
      INTEST=99
      iplate=0
      nchar_raster=1
      film_unit=4
      NSPOT=0
      NBG=20
      DO 10 I=1,3
      FIDX(I)=FIDXC1(I)
      FIDY(I)=FIDYC1(I)
10    CONTINUE
C
C======INITIALISE SCANNER VARIABLES
C
      NX=2400
      NY=2400
      IFLAG=1
C
C=====DIALOGUE WITH USER RE PARAMETERS
      write(out_term,2000)
2000  format (' >>> Film (F) or Image Plate (I): ',$)
      read(in_term,1002)ans
       if(ans.eq.'I'.or.ans.eq.'i') then 
        iplate=1
        nchar_raster=2
       endif
100   WRITE(6,1001)
1001  FORMAT(' >>> Enter filename of image : ',$)
      READ(5,1002)FNAME
      NRECST=1
120   NX=2400
      NY=2400
      WRITE (6,1003)
1003  FORMAT(' >>> No. of x rasters, no. of y rasters : ',$)        
      READ(5,*) NX,NY
      x_cen_f = float(nx/2)
      y_cen_f = float(ny/2)
140   CONTINUE
      OPEN(UNIT=film_unit,FILE=FNAME,ACCESS='DIRECT',
     +RECL=(NY*nchar_raster)/4)
      CLOSE(film_unit)
      ICX=NX/2
      ICY=NY/2
      DO 150 I=1,NY
      IBB(I)=0
150   CONTINUE
1002  FORMAT(A)
       if (iplate.eq.0) then
180     WRITE (6,1004)
1004    FORMAT (' >>> Camera type (C1 or C2) or I (to input fids):',$)
        READ(5,'(A)')ANS
          IF (ANS.EQ.'I'.OR.ANS.EQ.'i') then
           WRITE(6,3010)
3010       format(' INPUT FIDX(1),FIDY(1),FIDX(2),FIDY(2),',
     +     'FIDX(3),FIDY(3) in 10 micron units:')
           READ(5,*) FIDX(1),FIDY(1),FIDX(2),FIDY(2),
     +          FIDX(3),FIDY(3)
          ELSEIF (ANS.EQ.'C1'.OR.ANS.EQ.'c1') then
           ICAM = 1
            DO 200 I=1,3
             FIDX(I)=FIDXC1(I)
             FIDY(I)=FIDYC1(I)
200         CONTINUE
          ELSE 
           ICAM = 2
            DO 220 I=1,3
             FIDX(I)=FIDXC2(I)
             FIDY(I)=FIDYC2(I)
220         CONTINUE
          ENDIF
         
C
250      fid_box_size=8.0
         WRITE(6,1010)
1010      FORMAT(' >>> Enter fiducials search box size in mm :',$)
          READ(5,*)fid_box_size
       endif
C
C=====Prepare data in suitable form for find_fiducials routines
C
260   continue
      WRITE (6,1011)
1011  FORMAT(' >>> Raster size(microns) : ',$)
      READ(5,*)raster_size
      IR10=NINT(10.0*raster_size)
      raster_size=IR10/10000.0
      mm_rast=1.0/raster_size
      x_rasters=NX
      y_rasters=NY
c
       if(iplate.eq.0) then
        DO 280 I=1,3
         fid_pos_x(I) = (FIDX(I)/100.0 *mm_rast) + ICX
         fid_pos_y(I) = (FIDY(I)/100.0 *mm_rast) + ICY
280     CONTINUE
       OPEN(UNIT=film_unit,FILE=FNAME,ACCESS='DIRECT',STATUS=
     * 'OLD',FORM='UNFORMATTED',RECL=(NX*nchar_raster)/4)
       CALL FINDFD (film_unit)
       CLOSE (film_unit)
      endif
C
300   WRITE(6,1005)
1005  FORMAT(' >>> Enter crystal to image dist. in mm: ',$)
      READ(5,*)CTFILM
      WRITE(6,1006)
1006  FORMAT(' >>> Enter above background threshold for spots: ',$)
      READ(5,*)JTHRESH
      WRITE(6,1007)
1007  FORMAT(' >>> Do you want to save spots in buffer [y]: ',$)
      READ(5,1002)IYN
      IF (IYN.EQ.'N'.OR.IYN.EQ.'n') IFLAG=0
C
C=====CHECK FOR BAD REGION OF FILM
C
      WRITE(6,1008)
1008  FORMAT(' >>> Is there a bad region? [n]: ',$)
      READ(5,1002)IYN
      IXMIN=10000
       IF(IYN.EQ.'Y'.OR.IYN.EQ.'y') THEN
        WRITE(6,1009)
1009    FORMAT(' >>> Enter xmin,xmax,ymin,ymax: ',$)
        READ(5,*)IXMIN,IXMAX,IYMIN,IYMAX
       ENDIF
350   CONTINUE
      AMARG=8.0
      WRITE(6,1012)
1012  FORMAT(' >>> Enter width of border to exclude (multiple of',
     *' 2.0 mm) in mm : ',$)
      READ(5,*)AMARG 
      AMARG=2.0*NINT(AMARG/2.0)
      MARGIN=NINT(AMARG*mm_rast)
c
c=====INITIALISE PLOT VARIABLES AND SCAN RANGE VARIABLES
c
       if(nx.ne.ny)then
        nmax=max0(nx,ny)
        idx=nmax/600
        idy=idx
       else 
        idx=nx/600
        idy=ny/600
       endif
      NPX=600-margin/idx
      NPY=600-margin/idy
      NMX=margin/idx
      NMY=NMX
      NREC=600*IDX-2*MARGIN
      NYREC=600*IDY-2*MARGIN
      ISTART=ICX-NREC/2+1
      IYST=ICY-NYREC/2+1
C   
C======CALCULATE RADIAL BACKGROUND CURVE
C
      NRECST=ICX-20
      OPEN(UNIT=film_unit,FILE=FNAME,RECL=(NY*nchar_raster)/4,
     +access='DIRECT')
      CALL BACK(MARGIN,itmp)
C
C=====INITIALISE TERMINAL
C
      IXTER=0
      WRITE(6,3000)
3000  FORMAT(' IS THIS A UNIX XTERM WINDOW [N]: ')
      READ(5,1002)IYN
      IF(IYN.EQ.'Y'.OR.IYN.EQ.'y') IXTER=1
      CALL QGINIT(0,0,IXTER,5,6)
C
C=====CALLS TO INITIALISE GRAPHICS
C
      CALL QGENTR
C
C=====DO THE DISPLAY
C
      K1=1
      K2=NPY
      NOMIT=0
      JVAL=0
      BAD=.FALSE.
      IBLK=ISTART
      DO 500 J=1,NPX
C
C=====SET BACKGROUND DIVISIONS IN Y FOR THIS SCAN
C
      CALL BKSET(IBLK)
C     CHECK FOR BAD REGION
      IF(BAD)GO TO 380
      IF(IBLK.LE.IXMIN.OR.IBLK.GT.IXMAX)GO TO 400
C     JUST REACHED BAD REGION OF FILM
      BAD=.TRUE.
      K1=(IYMIN-IYST)/IDY
      IF(K1.LT.1)K1=1
      K2=(IYMAX-IYST)/IDY
      IF(K2.GT.NPY)K2=NPY
      NOMIT=IYMIN-IYST
      IF(NOMIT.LT.0)NOMIT=0
      GO TO 400
C
C=====END OF BAD PATCH?
C
380   IF(IBLK.LE.IXMAX)GO TO 400
      BAD=.FALSE.
      K1=1
      K2=NPY
      NOMIT=0
400   CONTINUE
      DO 480 JJ=1,IDX
      ITHRESH=IBGV(1)+JTHRESH
      NB=2
      IBY=IBGY(1)
      CALL READBK(film_unit,IBLK,itmp)
      IBLK=IBLK+1
C
C=====LOOP OVER RECORD AND PLOT
C
      VEC=.FALSE.
      IADDR=IYST+NOMIT
C
C=====FIND START OF BACKGROUND
C
      IF(NOMIT.EQ.0)GO TO 420
410   IF(IADDR.LT.IBY)GO TO 420
      IBY=IBGY(NB)
      ITHRESH=IBGV(NB)+JTHRESH
      NB=NB+1
      GO TO 410
420   CONTINUE
      DO 460 K=K1,K2
      IVAL=0
      DO 430 KK=1,IDY
      IF(IBB(IADDR).GT.IVAL)IVAL=IBB(IADDR)
      IADDR=IADDR+1
C   C   IBB(IADDR)=IBB(IADDR)-ITHRESH
430   CONTINUE
      IF(IADDR.LT.IBY)GO TO 440
      ITHRESH=IBGV(NB)+JTHRESH
      IBY=IBGY(NB)
      NB=NB+1
C
440   IF(IVAL.LE.ITHRESH)GO TO 450
      IF(VEC)JVAL=MAX0(JVAL,IVAL)
      IF(VEC)GO TO 460
      IB=K
      L1=IADDR-2
      VEC=.TRUE.
      rx=float(j+nmx)
      ry=float(ib+nmy)
      call qgmove(rx,ry)
      JVAL=IVAL
      GO TO 460
450   IF(.NOT.VEC)GO TO 460
      IE=(K-1)
      VEC=.FALSE.
      rx=float(j+nmx)
      ry=float(ie+nmy)
      call qgdraw(rx,ry)
      CURRY=(IE+IB)/2.
      CURRY=CURRY*IDY+IYST
      ICURRX=IBLK-1
      KVAL=JVAL-ITHRESH
      JVAL=0
      L2=IADDR
      IF(IFLAG.EQ.1)CALL OUTSPOT(ICURRX,CURRY,KVAL,L1,L2,ITHRESH)
460   CONTINUE
C
470   IF(.NOT.VEC)GO TO 480
      rx=float(j+nmx)
      ry=float(nmy+npy)
      call qgdraw(rx,ry)
480   CONTINUE
C
500   CONTINUE
C
C=====OUTPUT SPOTS IF REQUIRED
C
      IF(IFLAG.EQ.0) go to 600
      chsize=15.0
      write(string,1120)
1120  format(' do you want to cursor in beam centre (n): ')
      call qgmove(5.0,650.0)
      call qgtext(5.0,750.0,chsize,string)
      read(5,1002)iyn
c
      if(iyn.eq.'y'.or.iyn.eq.'Y') then
       CALL QGCURS(x,y,char)
       x_cen_f=x*idx
       y_cen_f=y*idy
       write(string,1126) x_cen_f,y_cen_f
       call qgtext(5.0,730.0,chsize,string)
      else
       continue
      endif
c
1126   format(' new centre at ',2f10.2)
c
      write(string,1101)
1101  format(' do you want to output spots (y): ')
      call qgtext(5.0,710.0,chsize,string)
      read(5,1002) iyn
      if(iyn.eq.'n'.or.iyn.eq.'N') goto 600
C
C     CHECK FOR OUTPUT FILE
C
520   WRITE(string,1102)
1102  FORMAT(' Enter name of output file: ')
      call qgtext(5.0,690.0,chsize,string)
      READ(5,1002)FILOUT
      call qgtext(250.0,690.0,chsize,FILOUT)
c
530   CONTINUE
      OPEN(UNIT=8,FILE=FILOUT,STATUS='NEW')
      CALL CLEAN
      DO 550 I=1,NSPOT
      IF(ISPOT(4,I).GT.0)GO TO 540
c      WRITE(6,1103)I,(ISPOT(JJJ,I),JJJ=1,6)
      ISPOT(4,I)=1
1103  FORMAT(7I8)
540   CONTINUE
C
C=====ADD BACK VIEWING THRESHOLD INTO INTEGRATED INTENSITY
C
      ISPOT(5,I)=ISPOT(5,I)+ISPOT(4,I)*JTHRESH
550   CONTINUE
C
C     OUTPUT TITLE
      WRITE(longstring,1105) NSPOT,jthresh,fname
      call qgtext(5.0,670.0,chsize,longstring)
1105  FORMAT(1x,i6,' spots (threshold=', i4,') read from ',a )
c
      WRITE(8,1106)CTFILM
1106  FORMAT(1X,F9.3)
      DO 560 I=1,NSPOT
      XX=(SPOT(1,I)-x_cen_f)/mm_rast
      YY=(SPOT(2,I)-y_cen_f)/mm_rast
      WRITE(8,1108)XX,YY,ISPOT(5,I),ISPOT(4,I)
1108  FORMAT(1X,2F9.3,2I6)
560   CONTINUE
      CLOSE(UNIT=8)
      call qgtmod
 600  CALL QGCLOS
      CLOSE(UNIT=2)
      STOP
      END
C
C
C**********************************************************
      SUBROUTINE BACK(MARGIN,itmp)
C
C     USE 40 RECORDS AT CENTRAL X TO GET
C     A RADIAL (DIAMETER REALLY) BACKGROUND CURVE
C     THE NBG (20) SYMMETRIC VALUES ARE SAVED IN IBGR
C     THE CENTRE AND EXTREME VALUES ARE BASED ON NEIGHBOURS
C
      integer max_rasters
      integer film_unit
      parameter (max_rasters=9600)
      COMMON/CBACK/NBG,IBGR(30),IBGY(30),IBGV(30)
      COMMON/CDEN/IBB(max_rasters)
      COMMON/XCONST/NX,NY,ICX,ICY,NREC,NYREC,ISTART,IYST
      common /ip/iplate,nchar_raster
      integer iplate,nchar_raster
      integer*4 itmp(ny*nchar_raster/4)
      DIMENSION ISUM(30)
      CHARACTER*1 IYN
      film_unit=4
      DO 10 I=1,30
      ISUM(I)=0
   10 CONTINUE
      NPERB=(NYREC-2*MARGIN)/NBG
      IBLK=ICX-20
C     SKIP OVER RECORDS TILL NEAR MIDDLE
C
      iadr=iblk
      DO 40 I=1,40
      CALL READBK(film_unit,IADR,itmp)
C
      L=IYST
      DO 30 J=1,NBG
      JSUM=0
      DO 20 II=1,NPERB
      JSUM=JSUM+IBB(L)
      L=L+1
   20 CONTINUE
      ISUM(J)=ISUM(J)+JSUM/NPERB
   30 CONTINUE
C
      iadr=iadr+1
   40 CONTINUE
      DO 42 II=1,NBG
      ISUM(II)=ISUM(II)/40
   42 CONTINUE
      WRITE(6,1010)(ISUM(J),J=1,NBG)
1010  FORMAT(/,' BACKGROUND VALUES ACROSS FILM:',/,' ',10I5,/,' ',10I5)
C
      NN=NBG
      N=NBG/2
      DO 50 I=1,N
      JSUM=MAX0(ISUM(I),ISUM(NN))
      IBGR(I)=JSUM
      IBGR(NN)=JSUM
      NN=NN-1
   50 CONTINUE
C
C     ADJUST EDGES
      IBGR(1)=IBGR(2)
      IBGR(NBG)=IBGR(NBG-1)
C     ADJUST CENTRE
      IBGR(N)=IBGR(N-1)+IABS(IBGR(N-1)-IBGR(N-2))
      IBGR(N+1)=IBGR(N)
      WRITE(6,1011)(IBGR(J),J=1,NBG)
1011  FORMAT(/,' SYMMETRICAL, EDGE ADJUSTED BACKGROUND VALUES:',
     * /,' ',10I5,/,' ',10I5)
      IBGR(NBG+1)=IBGR(NBG)
      WRITE (6,1012)
      READ (5,1013) IYN
1012  FORMAT (' >>> Do you wish to override these values [n]: ',$)
1013  FORMAT (A)
      IF (IYN.NE.'Y'.AND.IYN.NE.'y') RETURN
      WRITE (6,1014) NBG
1014  FORMAT (' >>> Enter',I3,' background values:')
      READ (5,*) (IBGR(I),I=1,NBG)
      IBGR(NBG+1)=IBGR(NBG)
      RETURN
      END
C
C
C**************************************************
C     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.E=ICX
C
      SUBROUTINE BKSET(IX)
      COMMON/CBACK/NBG,IBGR(30),IBGY(30),IBGV(30)
      COMMON/XCONST/NX,NY,ICX,ICY,NREC,NYREC,ISTART,IYST
      DATA NPERB/105/
C
      IBGY(1)=10000
      IBGV(1)=IBGR(1)
      JX=IABS(ICX-IX)
      NST=105*9
C
      DO 50 I=1,30
      IF(NST.LT.JX)GO TO 60
      DIF=NST*NST-JX*JX
      IBGY(I)=ICY-SQRT(DIF)
      IBGV(I)=IBGR(I)
      NST=NST-NPERB
   50 CONTINUE
C
   60 N=I-1
      IF(N.LE.0)GO TO 80
      NN=2*N
      IF(IBGY(N).EQ.ICY)NN=NN-1
      IBGY(NN+1)=10000
      DO 70 I=1,N
      IBGY(NN)=2*ICY-IBGY(I)
      IBGV(NN+1)=IBGV(I)
      NN=NN-1
   70 CONTINUE
      IBGV(NN+1)=IBGR(N+1)
   80 CONTINUE
      N=2*N+1
      RETURN
      END
C********************************************************
      SUBROUTINE OUTSPOT(ICURRX,CURRY,JVAL,L1,L2,ITHRESH)
C
C     COLLECT UP ARRAY OF SPOTS
C     SPOT(2,N) HAS CURRECT XC,YC FOR SPOT
C     ISPOT(2,N) HAS NO OF ROWS OF X IN SPOT,LAST X
C     ( 0 MEANS SPOT IS FINISHED)
C     ISPOT(3,N) MAX INTENSITY 
C     ISPOT(4,N) NO. OF POINTS IN SPOT
C     ISPOT(5,N) INTEG INTENSITY, CORRECTED FOR BACKGROUND
C
C Parameter for maximum no. of spots which can be stored
      PARAMETER (MAXSPT=20000)
      integer max_rasters
      parameter (max_rasters=9600)
C
C
      COMMON/CSPOT/NSPOT,SPOT(2,MAXSPT),ISPOT(6,MAXSPT),IFLAG
      COMMON/CDEN/IBB(max_rasters)
      DATA ISF/1/
      DATA DELY/5./
C
      L=ISF-1
C
   10 L=L+1
      IF(L.GT.NSPOT)GO TO 40
      IF(ISPOT(2,L).EQ.0)GO TO 16
      DIF=ABS(SPOT(2,L)-CURRY)
      IF(DIF.LE.DELY)GO TO 20
C
C     SPOT DOESNT MATCH ,
C     SWITCH OFF IF Y<CURRY
      IF(CURRY.LT.SPOT(2,L))GO TO 30
      IF(ISPOT(2,L).EQ.ICURRX)GO TO 10
   15 ISPOT(2,L)=0
   16 IF(L.EQ.ISF)ISF=ISF+1
      GO TO 10
C
C     CONTINUE AN EXISTING SPOT
   20 N=ISPOT(1,L)
      SPOT(1,L)=(N*SPOT(1,L)+ICURRX)/(N+1)
      SPOT(2,L)=(N*SPOT(2,L)+CURRY)/(N+1)
      ISPOT(1,L)=N+1
      IF(ISPOT(2,L).NE.ICURRX-1)GO TO 99
   22 ISPOT(2,L)=ICURRX
      ISPOT(3,L)=MAX0(ISPOT(3,L),JVAL)
      GO TO 90
C
C     SPOT GREATER IN Y
   30 IF(ispot(2,L).EQ.ICURRX-1)GO TO 10
      GO TO 15
C
C     ADD IN A NEW SPOT
   40 NSPOT=NSPOT+1
      IF(NSPOT.LE.MAXSPT)GO TO 45
      CALL CLEAN
      DO 41 II=1,NSPOT
      IF(ISPOT(2,II).NE.0)GO TO 42
   41 CONTINUE
      II=NSPOT+1
   42 ISF=II

      IF(NSPOT.LE.MAXSPT)GO TO 45
      IFLAG=0
   44 RETURN
   45 SPOT(1,NSPOT)=ICURRX
      SPOT(2,NSPOT)=CURRY
       ISPOT(1,NSPOT)=1
      ISPOT(2,NSPOT)=ICURRX
      ISPOT(3,NSPOT)=JVAL
      ISPOT(4,NSPOT)=0
      ISPOT(5,NSPOT)=0
      ISPOT(6,NSPOT)=ITHRESH
      L=NSPOT
C
   90 CONTINUE
C     ADD SPOT TO GET INTEGRATED VALUE
      ISUM=0
      JSUM=0
      DO 96 I=L1,L2
      JUNK=IBB(I)-ITHRESH
      IF(JUNK.LE.0)GO TO 96
      ISUM=ISUM+1
      JSUM=JSUM+JUNK
   96 CONTINUE
      ISPOT(4,L)=ISPOT(4,L)+ISUM
      ISPOT(5,L)=ISPOT(5,L)+JSUM
      RETURN
C   99 WRITE(6,100)ICURRX,CURRY,L,SPOT(1,L),SPOT(2,L),
   99 CONTINUE
C     1 ISPOT(1,L),ISPOT(2,L)
  100 FORMAT(1X,I5,1X,F8.2,I5,1X,2F8.2,2I5)
      GO TO 22
      END
C
C*********************************************************
      SUBROUTINE CLEAN
C
C Parameter for maximum no. of spots which can be stored
      PARAMETER (MAXSPT=20000)
C
      COMMON/CSPOT/NSPOT,SPOT(2,MAXSPT),ISPOT(5,MAXSPT),IFLAG
      IC=0
      DO 10 I=1,NSPOT
      IF(ISPOT(1,I).GT.1)GO TO 15
      IF(ISPOT(2,I).GT.0)GO TO 15
      IF(ISPOT(3,I).LE.20)GO TO 10
   15 IC=IC+1
      IF(I.EQ.IC)GO TO 10
      DO 20 J=1,2
   20 SPOT(J,IC)=SPOT(J,I)
      DO 22 J=1,5
   22 ISPOT(J,IC)=ISPOT(J,I)
   10 CONTINUE
      NSPOT=IC
      RETURN
      END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE  FINDFD (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
C ***Modified version of the routine find-fiducials from GENLAUE***
C
      implicit none
C
      integer film_unit, max_rasters
      common /fiducials/ x_cen_f, y_cen_f, w_f, x_c, y_c, w_c,
     1                   fidx, fidy, fid_pos_x, fid_pos_y
      real               x_cen_f, y_cen_f, w_f, x_c, y_c, w_c
C     Units:             rasters  rasters  rasters       rasters
      real               fidx(3), fidy(3), fid_pos_x(3), fid_pos_y(3)
      parameter (max_rasters=9600)
C Size in pixels of scanned area of film.
      integer x_rasters, y_rasters
C Multiply by mm_rast to convert mm to pixels.
C   Units: mm           pix/mm   mm
      real raster_size, mm_rast, fid_box_size
      common /sizes/ raster_size, x_rasters, y_rasters, mm_rast,
     1 fid_box_size
      common /ip/ iplate,nchar_raster
      integer iplate,nchar_raster
      common /io_streams/ in_term, out_term, in_data, out_data
      integer             in_term, out_term, in_data, out_data
C
C in_term : main input stream (usually the terminal =5)
C out_term : main output (printer) stream (usually =6)
C in_data : input stream for data (to be written to files)
C out_data : output stream for data (to be written to files)
C
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local Variables:
C
      integer  i
      real       rad_to_deg
      parameter (rad_to_deg = 360./(2.*3.1415926))
C
C------------------------------------------------------------------------------
C
      write (out_term, 2000)
      do i=1,3
         call  ffid (film_unit, fid_pos_x(i), fid_pos_y(i),
     1                   fidx(i), fidy(i) )
         if  (fidx(i).eq.-1)  then
            write (out_term, 9000) i
            stop
         end if
      end do
C
      write (out_term, 1000)  fidx(1), fidy(1),
     1                        fidx(2), fidy(2),
     2                        fidx(3), fidy(3)
      x_cen_f = (fidx(1) + fidx(3)) / 2.0
      y_cen_f = (fidy(1) + fidy(3)) / 2.0
      w_f = (fidy(3) - fidy(2)) / (fidx(3) - fidx(2))
      w_f = atan(w_f) * rad_to_deg
C
      write(out_term,4000)  raster_size*1000., x_cen_f, y_cen_f, w_f
      return
C
 1000 format(' Fiducials found at:  (',f7.1,',',f7.1,')',
     1   2(/,'                      (',f7.1,',',f7.1,')'))
 2000 format(' === Locating fiducials')
 4000 format(' Fiducial centre (',f5.1,'u units) ', 2g14.6/' Omega',
     1       ' as determined from fiducials (degrees)'/g14.6)
 9000 format(1x,' !!! Failed to locate fiducial ',i1,' - Program',
     1' terminated')
      end
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE ffid(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
C ***Modified version of subroutine find_fid from GENLAUE***
C
      implicit none
C
      integer  film_unit
      real     rough_x, rough_y
      real     fidx, fidy
      CHARACTER*1 IYN
C
C Note: if the routine fails to locate the fiducial within 5mm then
C       fidx and fidy are both -1
C
      integer max_rasters
      parameter (max_rasters=9600)
C Size in pixels of scanned area of film.
      integer x_rasters, y_rasters
C Multiply by mm_rast to convert mm to pixels.
C   Units: mm           pix/mm   mm
      real raster_size, mm_rast, fid_box_size
      common /sizes/ raster_size, x_rasters, y_rasters, mm_rast,
     1 fid_box_size
      common /io_streams/ in_term, out_term, in_data, out_data
      integer             in_term, out_term, in_data, out_data
C
C in_term : main input stream (usually the terminal =5)
C out_term : main output (printer) stream (usually =6)
C in_data : input stream for data (to be written to files)
C out_data : output stream for data (to be written to files)
C
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local Variables:
C
      real*4     max_x_size, max_shift, max_shift2
C
      integer*4  xc, yc, x_size, y_size
      integer*4  thresh, num_over
      real*4     cgx, cgy
      logical    try_again
      byte strip(max_rasters)
C
C------------------------------------------------------------------------------
C max size is twice starting size
C NB fid box is actually fid_box_size/2
      max_x_size = fid_box_size *mm_rast
      max_shift = fid_box_size/4. *mm_rast
      max_shift2 = max_shift * max_shift
      thresh = 125
      try_again = .false.
C
      xc = rough_x
      yc = rough_y
 1    x_size = fid_box_size/2. *mm_rast
      y_size = x_size
C
C========  Assume background is 25 - there is no reasoning behind this value
C
 10   call  COFG (film_unit, xc, yc, x_size, y_size, thresh,
     1              cgx, cgy, num_over,strip)
C
C========  If flat film, then enlarge box slighty
C          if box gets to 10mm x 10mm then abort
C                        = 200 x 200
C
      if  (num_over.eq.0)  then
C
         if (.not. try_again) then
           try_again = .true.
c next 4 lines replace call to ask (CCP4)
        WRITE(OUT_TERM,200)
200   FORMAT(' Couldnt find fiducial - try with bigger box/',
     +'lower threshold?')
        READ(IN_TERM,210)iyn
210     format(a)
        IF (IYN.EQ.'n'.OR.IYN.EQ.'N') then
            write (out_term,120) rough_x, rough_y
 120        format(' Using ESTIMATED fiducial pos. of ', 2f8.1)
            fidx = rough_x
            fidy = rough_y
            return
           end if
         end if
         if  (x_size.gt.max_x_size)  then
C
C========  Can`t find spot in big box - try lowering threshold (weak fiducial?)
C
            if (thresh.gt.75) then
              thresh = thresh - 25
              write (out_term, 100) thresh
  100         format (' Could not find fiducial - lowering',
     1         ' threshold to ',i3)
              goto 10
            else
              write (out_term,110) thresh
 110          format(' STILL no density found above ', i4)
              thresh = thresh - 25
              write(out_term,115)
115           format(1x,'Enter threshold or 0 to give up',$)
              read(in_term,210)thresh
              if (thresh.gt.0) goto 10
            end if
C
C========  Failed to locate fiducial - give up!
C
            fidx = -1
            fidy = -1
            return
         else
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 = x_size
            goto 10
C
         end if
C
      end if
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
C
         xc = cgx
         yc = cgy
         goto 1
C
      end if
C
C========  OK found a spot that hasn't shifted far from box centre
C
      fidx = cgx
      fidy = cgy
C
      return
      end

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE  COFG (film_unit, xc, yc, x_size, y_size, thresh,
     1                    cgx, cgy, num_over,strip)
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
C ***Modified version of subroutine c_of_g from GENLAUE***
c
C only used for film not image plates
c
      implicit none
C
      integer*4  film_unit
      integer*4  xc, yc, x_size, y_size 
      integer*4  thresh, num_over
      real*4     cgx, cgy
C
C     num_over: number of spots in the box over the threshold.
C
      integer max_rasters
      parameter (max_rasters=9600)
C Size in pixels of scanned area of film.
      integer x_rasters, y_rasters
C Multiply by mm_rast to convert mm to pixels.
C   Units: mm           pix/mm   mm
      real raster_size, mm_rast, fid_box_size
      common /sizes/ raster_size, x_rasters, y_rasters, mm_rast,
     1 fid_box_size
      common /ip/iplate,nchar_raster
      integer iplate,nchar_raster
      common /io_streams/ in_term, out_term, in_data, out_data
      integer           in_term, out_term, in_data, out_data
C
C in_term : main input stream (usually the terminal =5)
C out_term : main output (printer) stream (usually =6)
C in_data : input stream for data (to be written to files)
C out_data : output stream for data (to be written to files)
C
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local Variables:
C
      integer*4  ix, iy, nstrip
c      byte       strip(max_rasters)
c debug 
      byte  strip(y_rasters)
c end debug
      integer*4  intens, i,ios
      real*4     sumx, sumy, sum_intens, dx, dy
C
C------------------------------------------------------------------------------
CMods: IJC Feb 86
      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.
     1     yc-y_size.lt.1 .or. yc+y_size.gt.y_rasters)  then
         write (out_term, 9000) xc, yc
         num_over = 0
         cgx = 0.0
         cgy = 0.0
         return
      end if
C
      do ix=-x_size,x_size
C
         read (film_unit, rec=nstrip) (strip(i), i=1, y_rasters)
c        READ(film_unit,rec=nstrip,IOSTAT=ios,err=100) strip
c         read (film_unit,rec=nstrip,err=200) (strip(i),i=
c     +                                     1,y_rasters/nchar_raster)
200      nstrip = nstrip + 1
C
         do  iy=-y_size,y_size
C
C###            intens = jzext(strip(yc+iy))
C### DL mod. replaces jzext call
            intens = strip(yc+iy)
            if (intens.lt.0) intens = 256 + intens

            intens = max(0, intens - thresh)
C
            if (intens.gt.0) num_over = num_over+1
            sumx = sumx + ix * intens
            sumy = sumy + iy * intens
            sum_intens = sum_intens + intens
C
         end do
C
      end do
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
100   WRITE(6,101) ios
101   format(' read failed iostat = ',i6)
 9000 format(' !!! Attempted to read beyond scanned film',/,
     1       '  xc =',i4,', yc=',i4)
      end
c
c
c
      SUBROUTINE READBK(IUNIT,IBLK,itmp)
c read and unpack NY integer values from .od file and 
c return in IBB
c
c replaces msbr_r and iunpack
c
      integer max_rasters
      PARAMETER (max_rasters=9600)
      COMMON/CDEN/IBB(max_rasters)
      COMMON/XCONST/NX,NY,ICX,ICY,NREC,NYREC,ISTART,IYST
      integer nx,ny
      common /ip/iplate,nchar_raster
      integer iplate,nchar_raster
      integer*4 itmp(ny*nchar_raster/4),IBYTE(4),i2byte(2)
      INTEGER*4 JUNK,junk1,junk2,junk3
c      DATA MASKA /B'11111111000000000000000000000000'/
c      DATA MASKB /B'00000000111111110000000000000000'/
c      DATA MASKC /B'00000000000000001111111100000000'/
c      DATA MASKD /B'00000000000000000000000011111111'/
      data maska /'ff000000'x/
      data maskb /'00ff0000'x/
      data maskc /'0000ff00'x/
      data maskd /'000000ff'x/
c
       READ(IUNIT,rec=IBLK,IOSTAT=ios,err=300)ITMP
c
40    if (iplate.eq.0) then
        DO 100 K=1,ny*nchar_raster/4
         JUNK=IAND(MASKA,ITMP(K))
         IBYTE(1)=ISHFT(JUNK,-24)
         JUNK=IAND(MASKB,ITMP(K))
         IBYTE(2)=ISHFT(JUNK,-16)
         JUNK=IAND(MASKC,ITMP(K))
         IBYTE(3)=ISHFT(JUNK,-8)
         IBYTE(4)=IAND(MASKD,ITMP(K))
          DO 50 IB=1,4
           NCOUNT=(4*K)-4+IB
           IBB(NCOUNT)=IBYTE(IB)
50        CONTINUE
100     CONTINUE
C
      elseif (iplate.eq.1) then
       do 200 j=1,ny*nchar_raster/4
        junk1=iand(maska,itmp(j))
        junk2=iand(maskb,itmp(j))
        junk3=ior(junk1,junk2)
        i2byte(1)=ishft(junk3,-16)
        junk1=iand(maskc,itmp(j))
        junk2=iand(maskd,itmp(j))
        junk3=ior(junk1,junk2)
        i2byte(2)=iand(junk3,itmp(j))
         do 250 ib=1,2
          ncount=(2*j)-2+ib
          ibb(ncount)=i2byte(ib)
250      continue
200     continue
       endif
      goto 350
c
300   write(6,310)ios
310   format('read error ,ios= ',i5)
c
350   RETURN
      END