C C C Unix version: Wed Oct 2 08:37:52 BST 1991 C ----------------------------------------------- C C This version - DATE: Thu Oct 24 10:26:37 BST 1991 C C Changes: C C 1) Invert image left to right so it is in the standard orientation for C MOSFLM. A.G.W.L. 1/10/91 C 2) Comment out call to CCPFYP 1/10/91 (For VMS) C C ================ PROGRAM IMSTILLS C ================ C C C This version specific for EMBL Hamburg image plates. C C This program is part the Imperial College oscillation film C measurement package. The program is used together with IDXREF C to determine the orientation of a crystal from "still" photographs. C This program looks for spots on a film and outputs a list of C coordinates that can then read by IDXREF. Peter Brick C C Original version October 1985 C Key-worded version with width checking and no constants file. C July 1987 C Radial Background following EMBL/BASEL version added November 1988, C A.G.W. Leslie C C Modified to add "pseudo fiducials" for consistency with film images. C This, of course, changes the format of the SPOTS file, so this must C be used in conjunction with the new version of REFIX. C A.G.W. Leslie 11/1/91 C C C Changes to finding radial background, now print values also C Check record number when finding background so that it does not C try to read strips outside the scanned image. C Set pixels with values gt 32767 to 32767 rather than 1600 C Convert keywords to upper case C Add keyword DLAB to set 1200*1200 image size with 1 header record, C and subroutine RHEAD to read (and ignore) header for Mar scanners. C Default beam centre is that for F1 native data. C Test input IXLEN, IYLEN is LE IXLENGTH,IYLENGTH C AGWL 11/9/91 C C C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN,CUTWYMAX,CUTWYMIN,RMAX, + RMIN,SPSIGFACT,XSPLIT,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,NBINR,NBINT,NFILM, + NPIXMAX,NPIXMIN LOGICAL LIST,PRINT C .. C .. Arrays in Common .. REAL PHI INTEGER BGOD INTEGER*2 IMAGE LOGICAL RDSTRIP CHARACTER IMFILE*60 C .. C .. Local Scalars .. INTEGER I,ISAVE LOGICAL EOF C .. C .. External Subroutines .. EXTERNAL FINDBG,OPENODS,PICKSPOTS,READDAT C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /FNAMES/IMFILE(3) COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN,IMAGE(IXLENGTH*IYLENGTH), + RDSTRIP(IYLENGTH),IPOINT,IHD COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX,NBINR,NBINT, + IBINLIM,PHI(3),NPIXMIN,NPIXMAX,CUTPIXMIN,CUTPIXMAX, + CUTWXMIN,CUTWXMAX,CUTWYMIN,CUTWYMAX,LIST C .. save C C *********** CALL CCPFYP C *********** C C EOF = .FALSE. C C WRITE (6,FMT=6000) 6000 FORMAT (20X,'PROGRAM STILLS') NFILM = 0 C C 10 ISAVE = NFILM + 1 C C ************* CALL READDAT(EOF) C ************* C DO 20 I = ISAVE,NFILM IFILM = I C C---- Write terminator for previous still C IF (IFILM.GT.1) THEN WRITE (2,FMT='(1X,5F8.2)') - 99.0,-99.0,-99.0,-99.0,-99.0 END IF C C ****************** CALL OPENODS(IMFILE(I)) C ****************** C WRITE (6,FMT=6002) I,PHI(I) 6002 FORMAT (' ',15X,'PROCESSING STILL ',I1,' AT A PHI OF',F8.3) C C call findfids(ier) C if(ier.ne.0) stop C C ****** CALL FINDBG C ****** C C if (radx.or.radxy) call findbgx C C ********** CALL PICKSPOTS C ********** C 20 CONTINUE C IF (NFILM.LT.3 .AND. .NOT.EOF) GO TO 10 C C---- Write terminator for last film C WRITE (2,FMT='(1X,5F8.2)') - 999.0,-999.0,-999.0,-999.0,-999.0 C C STOP END SUBROUTINE READDAT(EOF) C ======================== C C C C---- Read in data cards and set defaults C C C C»»» C»»» C C C C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalar Arguments .. LOGICAL EOF C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN, + CUTWYMAX,CUTWYMIN,PIXEL,RMAX,RMIN,SPSIGFACT, + THRDIR,THRFID,XC,XSPLIT,YC,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,N1OD, + NBINR,NBINT,NFILM,NPIXMAX,NPIXMIN LOGICAL LIST,PRINT C .. C .. Arrays in Common .. REAL PHI INTEGER BGOD INTEGER*2 IMAGE LOGICAL RDSTRIP CHARACTER IMFILE*60 C .. C .. Local Scalars .. REAL STEP,TEMP,XC1,XC2,XC3,XC4,XC5,YC1,YC2,YC3,YC4, + YC5 INTEGER I,IXLEN1,IXLEN2,IXLEN3,IXLEN4,IYLEN1,IYLEN2, + IYLEN3,IYLEN4,LUNIN,LUNOUT,NFILE,NTOK LOGICAL SITE CHARACTER K1*4,K2*4,KEY*4,LINE*400 C .. C .. Local Arrays .. REAL VALUE(200) INTEGER IBEG(200),IDEC(200),IEND(200),ITYP(200) C .. C .. External Subroutines .. EXTERNAL CCPUPC,KEYERR,KEYNUM,PARSER,PRINTDAT C .. C .. Intrinsic Functions .. INTRINSIC NINT C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /COORDS/XC,YC,PIXEL COMMON /FNAMES/IMFILE(3) COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN, + IMAGE(IXLENGTH*IYLENGTH),RDSTRIP(IYLENGTH), + IPOINT,IHD COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX, + NBINR,NBINT,IBINLIM,PHI(3),NPIXMIN,NPIXMAX, + CUTPIXMIN,CUTPIXMAX,CUTWXMIN,CUTWXMAX, + CUTWYMIN,CUTWYMAX,LIST COMMON /TSTILL/N1OD,THRFID,THRDIR C .. save C .. Data statements .. C C Local variables C C C---- add XC4,YC4 for the new black scanner 4/1/91. C These values are after the camera alignment 15/12/90 C XC5,YC5 are for F1 data at DLAB C DATA XC1/100.12/,YC1/99.71/,PIXEL/0.187/, + IXLEN1/1068/,IYLEN1/1068/,XC2/108.627/, + YC2/108.192/,IXLEN2/1160/,IYLEN2/1160/, + XC3/111.08/,YC3/112.11/,IXLEN3/1186/, + IYLEN3/1186/,IXLEN4/1200/,IYLEN4/1200/, + XC4/87.65/,YC4/88.94/,XC5/92.69/,YC5/91.33/ DATA RMAX/90.0/,RMIN/10.0/,XSPLIT/0.3/,YSPLIT/0.3/ DATA IBINLIM/50/,NBINR/4/,NBINT/8/ DATA SPSIGFACT/10/,NPIXMIN/5/,NPIXMAX/250/ DATA CUTPIXMIN/0.5/,CUTPIXMAX/10.0/ DATA CUTWXMIN/0.5/,CUTWXMAX/5.0/ DATA CUTWYMIN/0.5/,CUTWYMAX/5.0/ DATA LUNIN/5/,LUNOUT/6/ DATA SITE/.FALSE./ C .. C C IHD = 0 YOFF = 0.0 C C ******************************************************** 10 CALL PARSER(LUNIN,LUNOUT,LINE,IBEG,IEND,ITYP,VALUE,IDEC,NTOK) C ******************************************************** C C---- End of file? C IF (NTOK.EQ.-1) GO TO 90 C C---- First 4 chars of 1st token C KEY = LINE(IBEG(1) :IEND(1)) C C---- Convert to upper case C C *********** CALL CCPUPC(KEY) C *********** C C---- PHI C === C IF (KEY.EQ.'PHI') THEN DO 20 I = 2,NTOK CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) NFILM = NFILM + 1 IF (NFILM.GT.3) THEN WRITE (6,FMT='(1X,''*** Only 3 Stills allowed ***'')') STOP END IF PHI(NFILM) = VALUE(I) 20 CONTINUE C C---- FILE C ==== C ELSE IF (KEY.EQ.'FILE') THEN DO 30 I = 2,NTOK NFILE = NFILE + 1 IF (NFILE.GT.3) THEN WRITE (6,FMT='(1X,''*** Only 3 Stills allowed ***'')') STOP END IF IMFILE(NFILE) = LINE(IBEG(I) :IEND(I)) 30 CONTINUE C C---- CENTer C ====== C ELSE IF (KEY.EQ.'CENT') THEN CALL KEYNUM(2,2,LINE,IBEG,IEND,ITYP,NTOK) XC = VALUE(2) YC = VALUE(3) C C---- PIXEl size of one pixel in mm C =============================== C ELSE IF (KEY.EQ.'PIXE') THEN CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) PIXEL = VALUE(2) C C---- XLENgth pixels in x direction C ================================== C ELSE IF (KEY.EQ.'XLEN') THEN CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) IXLEN = NINT(VALUE(2)) IF (IXLEN.GT.IXLENGTH) THEN WRITE (6,FMT=6000) IXLEN,IXLENGTH 6000 FORMAT (//1X, + 'Requested IXLEN of',I6,' is greater than maximum allowed (',I6, + ')',/1X,'Recompile program after changing parameter IXLENGTH') STOP END IF C C----- YLENgth pixels in y direction C ================================== C ELSE IF (KEY.EQ.'YLEN') THEN CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) IYLEN = NINT(VALUE(2)) IF (IYLEN.GT.IYLENGTH) THEN WRITE (6,FMT=6002) IYLEN,IYLENGTH 6002 FORMAT (//1X, + 'Requested IYLEN of',I6,' is greater than maximum allowed (',I6, + ')',/1X,'Recompile program after changing parameter IYLENGTH') STOP END IF C C---- SITE... can be DLAB, EMBL, LMB C ================================ C C for EMBL, subkeywords are SCR1,SCR2,SCR3,MAR for different C scanners. DLAB is assumed to be Mar. C C ELSE IF (KEY.EQ.'SITE') THEN SITE = .TRUE. C C---- Get site (second keyword) C K1 = LINE(IBEG(2) :IEND(2)) C C---- Convert to upper case C C *********** CALL CCPUPC(K1) C *********** C IF (K1.EQ.'EMBL') THEN C C---- Get scanner (third keyword) C IF (NTOK.LT.3) THEN WRITE (6,FMT=6004) 6004 FORMAT (//1X, + 'ERROR, MUST GIVE SCANNER TYPE AS SCR1,SCR2,SCR3 or MAR') STOP END IF K2 = LINE(IBEG(3) :IEND(3)) C C---- Convert to upper case C C *********** CALL CCPUPC(K2) C *********** C C IMAGE (VERSION 1.0) C IF (K2.EQ.'SCR1') THEN XC = XC1 YC = YC1 IXLEN = IXLEN1 IYLEN = IYLEN1 PIXEL = 0.187 C C SCRAP (VERSION 2.0) C ELSE IF (K2.EQ.'SCR2') THEN XC = XC2 YC = YC2 IXLEN = IXLEN2 IYLEN = IYLEN2 PIXEL = 0.187 C C SCRAP (VERSION 3.0) C ELSE IF (K2.EQ.'SCR3') THEN XC = XC3 YC = YC3 IXLEN = IXLEN3 IYLEN = IYLEN3 PIXEL = 0.187 C C---- Mar Research scanner C ELSE IF (K2.EQ.'MAR') THEN IHD = 1 IXLEN = IXLEN4 IYLEN = IYLEN4 PIXEL = 0.15 ELSE WRITE (6,FMT=6006) 6006 FORMAT (1X,'SCANNER NOT RECOGNISED',/1X, + 'Must be SCR1,SCR2',',SCR3 or MAR') STOP END IF C C LMB SCANNER C ELSE IF (K1.EQ.'LMB') THEN XC = XC4 YC = YC4 PIXEL = 0.150 IF (NTOK.EQ.4) THEN C C---- Allow for user specified direct beam coordinates C XC = VALUE(3) YC = VALUE(4) END IF IXLEN = IXLEN3 IYLEN = IYLEN3 C C---- Mar scanner at DLAB C ELSE IF (K1.EQ.'DLAB') THEN XC = XC5 YC = YC5 PIXEL = 0.15 IF (NTOK.EQ.4) THEN C C---- Allow for user specified direct beam coordinates C XC = VALUE(3) YC = VALUE(4) END IF IHD = 1 IXLEN = IXLEN4 IYLEN = IYLEN4 ELSE WRITE (6,FMT=6008) K1 6008 FORMAT (1X,'SITE',A,' NOT RECOGNISED',/1X, + 'Must be one of: EMBL,DLAB,LMB') STOP END IF C C ORIEntation C C ELSE IF (KEY.EQ.'ORIE') THEN C IF(NTOK.LT.2) GOTO 71 C K1 = LINE(IBEG(2):IEND(2)) C C ORIEntation ROTAted C C IF(K1.EQ.'ROTA') THEN C DO II = 1,3 C TEMP = XYFID(1,II) C XYFID(1,II) = XYFID(2,II) C XYFID(2,II) = -TEMP C ENDDO C C ORIEntation STANdard- the default C C ELSE IF(K1.EQ.'STAN') THEN C C ORIEntation ?????? C C ELSE C CALL KEYERR(I,0,LINE,IBEG,IEND,ITYP) C STOP C ENDIF C 40 CONTINUE C C---- BACKground C =========== C ELSE IF (KEY.EQ.'BACK') THEN I = 1 42 I = I + 1 IF(I.GT.NTOK) GOTO 44 ! Skip if no more tokens on line K1 = LINE(IBEG(I):IEND(I)) C C BACKground SIZE C c IF(K1.EQ.'SIZE') THEN c I = I + 1 c CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) c RBG = VALUE(I) C C BACKground XOFF C c ELSE IF(K1.EQ.'XOFF') THEN c I = I + 1 c CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) c XOFF = VALUE(I) cC C BACKground YOFF C IF(K1.EQ.'YOFF') THEN I = I + 1 CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) YOFF = VALUE(I) C C BACKground RATIo C c ELSE IF(K1.EQ.'RATI') THEN c I = I + 1 c CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) c RATIO = VALUE(I) C C BACKground STEPS C c ELSE IF(K1.EQ.'STEP') THEN c I = I + 1 c CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) c NSTEP = VALUE(I) C C BACKground PRINT C C ELSE IF(K1.EQ.'PRIN') THEN C PRINT=.TRUE. C C BACKground RADIal X, Y or XY C c ELSE IF(K1.EQ.'RADI') THEN c I = I + 1 c STR=LINE(IBEG(I):IEND(I)) c RADX=(STR.EQ.'X ') cc RADY=(STR.EQ.'Y ') c RADXY=(STR.EQ.'XY') C C BACKground FIXEd C C ELSE IF(K1.EQ.'FIXE') THEN C FIXED=.TRUE. C I = I + 1 C CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) C XBG = VALUE(I) C I = I + 1 C YBG = VALUE(I) C C BACKground ?????? C ELSE CALL KEYERR(I,0,LINE,IBEG,IEND,ITYP) STOP ENDIF C GOTO 42 44 CONTINUE C C C----- Now for the spot search and selection parameters C C---- RADIus C ====== C ELSE IF (KEY.EQ.'RADI') THEN IF (NTOK.EQ.3) THEN CALL KEYNUM(2,2,LINE,IBEG,IEND,ITYP,NTOK) RMAX = VALUE(2) RMIN = VALUE(3) IF (RMAX.LT.RMIN) THEN TEMP = RMAX RMAX = RMIN RMIN = TEMP END IF ELSE CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) RMAX = VALUE(2) END IF C C---- SPLIt C ===== C ELSE IF (KEY.EQ.'SPLI') THEN IF (NTOK.EQ.3) THEN CALL KEYNUM(2,2,LINE,IBEG,IEND,ITYP,NTOK) XSPLIT = VALUE(2) YSPLIT = VALUE(3) ELSE CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) XSPLIT = VALUE(2) YSPLIT = XSPLIT END IF C C---- BINS C ==== C ELSE IF (KEY.EQ.'BINS') THEN I = 1 50 I = I + 1 C C---- Skip if no more tokens on line C IF (I.GT.NTOK) GO TO 60 K1 = LINE(IBEG(I) :IEND(I)) C C---- Convert to upper case C C *********** CALL CCPUPC(K1) C *********** C C---- BINS OSC or STILL C ================= C IF (K1.EQ.'OSC ') THEN IBINLIM = 100 NBINR = 4 ELSE IF (K1.EQ.'STIL') THEN IBINLIM = 3 NBINR = 2 C C---- BINS SPOTs C ========== C ELSE IF (K1.EQ.'SPOT') THEN I = I + 1 CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) IBINLIM = NINT(VALUE(I)) C C---- BINS RADIal C =========== C ELSE IF (K1.EQ.'RADI') THEN I = I + 1 CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) NBINR = NINT(VALUE(I)) C C---- BINS THETa C ========== C ELSE IF (K1.EQ.'THET') THEN I = I + 1 CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) NBINT = NINT(VALUE(I)) C C---- BINS ?????? C =========== C ELSE CALL KEYERR(I,0,LINE,IBEG,IEND,ITYP) STOP END IF C GO TO 50 60 CONTINUE C C---- THREshold C ========== C ELSE IF (KEY.EQ.'THRE') THEN CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) SPSIGFACT = VALUE(2) C C---- LIST reflections C ================ C ELSE IF (KEY.EQ.'LIST') THEN LIST = .TRUE. C C----- SELEct C ======= C ELSE IF (KEY.EQ.'SELE') THEN I = 1 70 I = I + 1 C C---- Skip if no more tokens on line C IF (I.GT.NTOK) GO TO 80 K1 = LINE(IBEG(I) :IEND(I)) C C---- Convert to upper case C C *********** CALL CCPUPC(K1) C *********** C C---- SELEct PIXEl C ============= C IF (K1.EQ.'PIXE') THEN I = I + 1 CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) CUTPIXMIN = VALUE(I) CUTPIXMAX = VALUE(I+1) IF (CUTPIXMIN.GT.CUTPIXMAX) THEN TEMP = CUTPIXMAX CUTPIXMAX = CUTPIXMIN CUTPIXMIN = TEMP END IF I = I + 1 C C---- SELEct X C ======== C ELSE IF (K1.EQ.'X') THEN I = I + 1 CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) CUTWXMIN = VALUE(I) CUTWXMAX = VALUE(I+1) IF (CUTWXMIN.GT.CUTWXMAX) THEN TEMP = CUTWXMAX CUTWXMAX = CUTWXMIN CUTWXMIN = TEMP END IF I = I + 1 C C---- SELEct Y C ======== C ELSE IF (K1.EQ.'Y') THEN I = I + 1 CALL KEYNUM(2,I,LINE,IBEG,IEND,ITYP,NTOK) CUTWYMIN = VALUE(I) CUTWYMAX = VALUE(I+1) IF (CUTWYMIN.GT.CUTWYMAX) THEN TEMP = CUTWYMAX CUTWYMAX = CUTWYMIN CUTWYMIN = TEMP END IF I = I + 1 C C---- SELEct MINPixel C ================ C ELSE IF (K1.EQ.'MINP') THEN I = I + 1 CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) NPIXMIN = NINT(VALUE(I)) C C----- SELEct MAXPixel C =============== C ELSE IF (K1.EQ.'MAXP') THEN I = I + 1 CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) NPIXMAX = NINT(VALUE(I)) C C---- SELEct ?????? C ============== C ELSE CALL KEYERR(I,0,LINE,IBEG,IEND,ITYP) STOP END IF C GO TO 70 80 CONTINUE C C---- NOW - Print out all parameters C =============================== C ELSE IF (KEY.EQ.'NOW ') THEN CALL PRINTDAT C C---- GO C === C ELSE IF (KEY.EQ.'GO ') THEN GO TO 100 C ELSE WRITE (6,FMT=6010) 6010 FORMAT (' ERROR - LINE IGNORED ') END IF GO TO 10 90 CONTINUE EOF = .TRUE. C 100 CONTINUE C C ITHRFID = NINT( THRFID * N1OD ) C ITHRDIR = NINT( THRDIR * N1OD ) c IF (RBG.EQ.-99.0) THEN c IF (FIXED) THEN c RBG=40.0 c ELSE c RBG=MAX(0.5,0.5*RMAX/NSTEP) c END IF c END IF c STEP=rmax/nstep C STEP = 1 C IF (.NOT.SITE) THEN WRITE (6,FMT=6012) 6012 FORMAT (//1X,'***** ERROR *****',/1X, + 'A SITE card must be given (DLAB,EMBL or LMB)') STOP END IF C C ************** IF (.NOT.EOF) CALL PRINTDAT C ************** C RETURN END REAL FUNCTION DOCALC(VALUE,OPER,VALUE0) C ================================== C C C C Do simple arithmetic on two arguments C C .. Scalar Arguments .. REAL VALUE,VALUE0 INTEGER OPER C .. C .. Local Scalars .. INTEGER DIVIDE,EENOTN,MINUS,MULT,PLUS C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Data statements .. DATA PLUS/1/,MINUS/2/,MULT/3/,DIVIDE/4/,EENOTN/5/ C .. C C IF (OPER.EQ.PLUS) DOCALC = VALUE0 + VALUE IF (OPER.EQ.MINUS) DOCALC = VALUE0 - VALUE IF (OPER.EQ.MULT) DOCALC = VALUE0*VALUE IF (OPER.NE.DIVIDE) GO TO 10 IF (VALUE.NE.0.0) DOCALC = VALUE0/VALUE IF (VALUE.EQ.0.0) DOCALC = 0.0 10 CONTINUE IF (OPER.NE.EENOTN) GO TO 20 IF (ABS(VALUE).LE.76.0) DOCALC = VALUE0* (10.0**VALUE) IF (ABS(VALUE).GT.76.0) DOCALC = 0.0 20 CONTINUE RETURN END SUBROUTINE FINDBG C ================= C C C C C---- This subroutine calculates a radial background from rmin to rmax C and returns BGOD containing the background values. C C C C---- convert mm dimensions to pixels C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) INTEGER NPIX PARAMETER (NPIX=51) C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN,CUTWYMAX,CUTWYMIN, + PIXEL,RMAX,RMIN,SPSIGFACT,XC,XSPLIT,YC,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,LIST,NBINR,NBINT, + NFILM,NPIXMAX,NPIXMIN LOGICAL PRINT C .. C .. Arrays in Common .. REAL PHI INTEGER BGOD INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Local Scalars .. INTEGER I,I1,I2,IBA,II,IRMAX,IRMIN,IXCEN,IYCEN,J,N,NBACK,NLI,NOFF, + NPT,NRJ,NRJLAST,SUM,IYOFF LOGICAL FIRSTTIME C .. C .. Local Arrays .. REAL RAD(800),SIGAV(800) INTEGER AV(800),IBGAV(800),IS(NPIX),ISIGBGAV(800),NREJ(800), + NREJAV(800),PIX(NPIX) C .. C .. External Subroutines .. EXTERNAL GETOD,RDBLK,SORTUP4 C .. C .. Intrinsic Functions .. INTRINSIC MIN,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /COORDS/XC,YC,PIXEL COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN,IMAGE(IXLENGTH*IYLENGTH), + RDSTRIP(IYLENGTH),IPOINT,IHD COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX,NBINR,NBINT, + IBINLIM,PHI(3),NPIXMIN,NPIXMAX,CUTPIXMIN,CUTPIXMAX, + CUTWXMIN,CUTWXMAX,CUTWYMIN,CUTWYMAX,LIST C .. save IXCEN = XC/PIXEL IYCEN = YC/PIXEL IRMIN = RMIN/PIXEL IRMAX = RMAX/PIXEL C C---- Radial background calculation for image plates C C Background calculated from NPIX pixels in strips in the area C from ixcen+rmin to ixcen+rmax C WRITE (6,FMT=6000) NPIX,YOFF 6000 FORMAT (1X,'Calculating radial background..... ',/,1X,'The ', + 'background is radial along the X (horizontal) direction', + /,1X,'The average of',I4,' pixels, excluding outliers,', + ' is used',/,1X,'The scan is offset by',F4.1,'mm ', + 'from the image centre in Y to avoid any',/,1X, + 'backstop shadow') C NOFF = NPIX/2 IYOFF = YOFF/PIXEL C C---- Loop over strips C DO 50 I = IXCEN + IRMIN,IXCEN + IRMAX II = I - IXCEN BGOD(II) = 0.0 C C---- test limits of I C IF ((I.LT.1) .OR. (I.GT.IXLEN)) GO TO 50 C C ******** CALL RDBLK(I) C ******** C FIRSTTIME = .TRUE. N = 0 SUM = 0 NRJ = 0 10 CONTINUE C C---- Loop over 51 pixels C DO 20 J = IYCEN - NOFF + IYOFF,IYCEN + NOFF + IYOFF C C ************* CALL GETOD(J,IBA) C ************* C IF (FIRSTTIME) THEN N = N + 1 PIX(N) = IBA ELSE Cal if (iba.le.(av(ii)+sigav(ii)*2))then IF (IBA.LE. (AV(II)+SIGAV(II)*1.5)) THEN C C---- Check if too big, ie probably spot C SUM = SUM + IBA N = N + 1 ELSE NRJ = NRJ + 1 END IF END IF 20 CONTINUE C C IF (FIRSTTIME) THEN C C---- Sort pixel values into ascending order C C ******************** CALL SORTUP4(NPIX,PIX,IS) C ******************** C C---- Set the background to the mean of the NPIX/3 smallest values, C excluding zeros C NBACK = 0 SUM = 0.0 C C DO 30 J = 1,NPIX IBA = PIX(IS(J)) IF (IBA.EQ.0) GO TO 30 NBACK = NBACK + 1 SUM = SUM + IBA IF (NBACK.GE.NPIX/3) GO TO 40 30 CONTINUE C C WRITE (6,FMT=6002) NPIX/3,I 6002 FORMAT (1X,'*** ERROR ****',/1X,'Less than',I3,' pixels with', + ' non-zero values for stripe',I6) STOP 40 AV(II) = SUM/NBACK ELSE IF (N.EQ.0) THEN AV(II) = 0 ELSE C C---- Average background C AV(II) = SUM/N END IF END IF SIGAV(II) = SQRT(REAL(AV(II))) C C type *,'ii,av,sig',ii,av(ii),sigav(ii) C C IF (FIRSTTIME) THEN C WRITE(6,2222) PIX 6004 FORMAT (1X,11I6) FIRSTTIME = .FALSE. SUM = 0 N = 0 NRJ = 0 NRJLAST = 0 GO TO 10 ELSE C C---- At end, calculate threshold for plot C IF (NRJ.EQ.NRJLAST) THEN C C type *,'ii, final av,sig',ii,av(ii),sigav(ii) C BGOD(II) = NINT((SIGAV(II)*SPSIGFACT)+AV(II)) NREJ(II) = NRJ ELSE SUM = 0.0 N = 0 NRJLAST = NRJ NRJ = 0 GO TO 10 END IF END IF 50 CONTINUE WRITE (6,FMT=6006) 6006 FORMAT (1X,'Radial Background calculation complete') C IF (PRINT) THEN C C---- Find average over 4 pixels (0.6mm) at a time C NPT = 0 C C---- Loop over strips C DO 60 I = IXCEN + IRMIN,IXCEN + IRMAX,4 NPT = NPT + 1 II = I - IXCEN IBGAV(NPT) = (AV(II+1)+AV(II)+AV(II+2)+AV(II+3))*0.25 ISIGBGAV(NPT) = (SIGAV(II+1)+SIGAV(II)+SIGAV(II+2)+SIGAV(II+3))* + 0.25 RAD(NPT) = II*PIXEL NREJAV(NPT) = (NREJ(II+1)+NREJ(II)+NREJ(II+2)+NREJ(II+3))*0.25 60 CONTINUE C C type *,'npt',npt C NLI = NPT/12 + 1 I1 = -11 C C DO 70 I = 1,NLI I1 = I1 + 12 I2 = I1 + 11 I2 = MIN(I2,NPT) WRITE (6,FMT=6008) (RAD(J),J=I1,I2) WRITE (6,FMT=6010) (IBGAV(J),J=I1,I2) WRITE (6,FMT=6012) (ISIGBGAV(J),J=I1,I2) WRITE (6,FMT=6014) (NREJAV(J),J=I1,I2) 70 CONTINUE C C 6008 FORMAT (1X,'RADIUS',12F6.1) 6010 FORMAT (1X,'BACKG ',12I6) 6012 FORMAT (1X,'SIGMA ',12I6) 6014 FORMAT (1X,'NREJ ',12I6) C C RETURN END SUBROUTINE FINDSPOTS(XMM,YMM,RMM,ITHRESH,WXMM,WYMM,MINPIX,NSPOT, + XSPOT,YSPOT,RSPOT,FSPOT,NPIX,IWXSPOT, + IWYSPOT,FID) C =========================================================== C C C C---- This version includes spot widths in x and y C C Scans film for spots within circle of radius RMM centered on XMM,YMM C Minimum dimensions of spot are given by WXMM and WYMM C Pixels in spot have to be above ITHRESH C C Routine returns NSPOT,XSPOT,YSPOT,RSPOT,FSPOT,NPIX,IWXSPOT,IWYSPOT C C This version has been made considerably more complicated in C an attempt to avoid problems with the masking tape being C treated as fiducials. C C The algorithm started off similar to that used by Bob Sweet. C C IMPLICIT NONE C C C---- Maximum active spots C C C---- Maximum number spots C C C Dummy arguments C C C C Local variables C C C .. Parameters .. INTEGER MXCENT PARAMETER (MXCENT=250) INTEGER MXSPOT PARAMETER (MXSPOT=5000) INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalar Arguments .. REAL RMM,WXMM,WYMM,XMM,YMM INTEGER ITHRESH,MINPIX,NSPOT LOGICAL FID C .. C .. Array Arguments .. REAL FSPOT(MXSPOT),RSPOT(MXSPOT),XSPOT(MXSPOT), + YSPOT(MXSPOT) INTEGER IWXSPOT(MXSPOT),IWYSPOT(MXSPOT),NPIX(MXSPOT) C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN, + CUTWYMAX,CUTWYMIN,PIXEL,RMAX,RMIN,SPSIGFACT, + XC,XSPLIT,YC,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,NBINR, + NBINT,NFILM,NPIXMAX,NPIXMIN LOGICAL LIST,PRINT C .. C .. Arrays in Common .. REAL PHI INTEGER BGOD INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Local Scalars .. REAL FACT,FINT,FITEMP,FXL,STEP,SX,SY,TSUMI, + TSUMYI,XSQ,YSQ INTEGER I,IB,IBA,ICENT,IDEL,IE,ILEN,IR,ITEMP,ITMP, + IX,IXCT,IXH,IXL,IXSPLIT,IY,IYH,IYL,IYMID, + IYPT,IYSPLIT,J,JDX,JDY,NCENT,NP,NPS,R1 LOGICAL ENDSTRIPE,LINE C .. C .. Local Arrays .. REAL SUMI(MXCENT),SUMXI(MXCENT),SUMYI(MXCENT) INTEGER IFLAG(MXCENT),IPT(MXCENT),IXCENT(MXCENT), + IXMAX(MXCENT),IXMIN(MXCENT),IYBEG(MXCENT), + IYCENT(MXCENT),IYEND(MXCENT),IYMAX(MXCENT), + IYMIN(MXCENT),NPIXCENT(MXCENT) C .. C .. External Subroutines .. EXTERNAL GETOD,RDBLK C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,NINT,REAL,SQRT C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /COORDS/XC,YC,PIXEL COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN, + IMAGE(IXLENGTH*IYLENGTH),RDSTRIP(IYLENGTH), + IPOINT,IHD COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX, + NBINR,NBINT,IBINLIM,PHI(3),NPIXMIN,NPIXMAX, + CUTPIXMIN,CUTPIXMAX,CUTWXMIN,CUTWXMAX, + CUTWYMIN,CUTWYMAX,LIST C .. save C C STEP = 1 C C---- number of pixels in a background step C NPS = NINT(STEP/PIXEL) C C---- conversion of pixels to background steps C FACT = PIXEL/STEP C C type *,'ithresh,nps,fact',ithresh,nps,fact C C---- Zero array, used as flag C DO 10 I = 1,MXCENT NPIXCENT(I) = 0 10 CONTINUE C C---- Calculate limits in pixel units C IX = NINT(XMM/PIXEL) IY = NINT(YMM/PIXEL) C C---- Limits on X C IR = NINT(RMM/PIXEL) IXL = MAX(1,IX-IR) IXH = MIN(IXLEN,IX+IR) C C---- Limits on spot dimensions for split spots C IXSPLIT = NINT(WXMM/PIXEL) IYSPLIT = NINT(WYMM/PIXEL) C NSPOT = 0 NCENT = 0 C C---- Outer loop - Loop over stripes/rows C WRITE (6,FMT=6000) 6000 FORMAT (2X,'Starting spot search...') C C DO 100 IXCT = IXL,IXH JDX = IXCT - IX XSQ = JDX*JDX C C *********** CALL RDBLK(IXCT) C *********** C FXL = REAL(IXCT) C C---- Determine limits of this stripe C IDEL = NINT(SQRT(REAL(IR*IR-JDX*JDX))) IYL = MAX(1,IY-IDEL) IYH = MIN(IYLEN,IY+IDEL) C C if (.not.fid) type *,'stripe, ylim',ixct,iyl,iyh C LINE = .FALSE. C C---- Loop over pixels along stripe and C find lines of pixels above threshold C At end of line, add line into active spot or start a new spot C NP = -1 C C DO 50 IYPT = IYL,IYH NP = NP + 1 C C **************** CALL GETOD(IYPT,IBA) C **************** C JDY = IYPT - IY YSQ = JDY*JDY R1 = NINT(SQRT(XSQ+YSQ)) C C IF (R1.EQ.0) THEN ITHRESH = 0 ELSE ITHRESH = BGOD(R1) END IF C C ITEMP = IBA C C IF (.NOT.LINE) THEN C C---- No line now or before C IF (ITEMP.LE.ITHRESH) THEN GO TO 50 ELSE C C---- (.not.line.and.itemp.gt.ithresh) ! New line C LINE = .TRUE. IB = IYPT FITEMP = REAL(ITEMP-ITHRESH) TSUMYI = REAL(IYPT)*FITEMP TSUMI = FITEMP END IF C C----! (itemp.le.ithresh) C ELSE C C---- ! (line) C ENDSTRIPE = (IYPT.EQ.IYH) C C---- ! End of line C IF (ITEMP.LE.ITHRESH .OR. ENDSTRIPE) THEN LINE = .FALSE. IE = IYPT - 1 ILEN = IE - IB + 1 IYMID = NINT(TSUMYI/TSUMI) C C---- Inner loop - for this line loop over active spots/centroids C Complicated test to allow for split spots. C C Line is considered to be part of spot if either C a) current line overlaps with line from previous stripe C b) or centroid of line (IYMID) lies C within IYSPLIT of current spot centroid C DO 20 ICENT = 1,NCENT J = IPT(ICENT) ITMP = ABS(IYMID-IYCENT(J)) C C---- Skip if not overlapping and not within IYSPLIT C IF (ABS((IYEND(J)+IYBEG(J))- (IE+IB)).GT. + ABS((IYEND(J)-IYBEG(J))+ (IE-IB)) .AND. + ITMP.GT.IYSPLIT) GO TO 20 C C SUMXI(J) = SUMXI(J) + TSUMI*FXL SUMYI(J) = SUMYI(J) + TSUMYI SUMI(J) = SUMI(J) + TSUMI IXCENT(J) = NINT(SUMXI(J)/SUMI(J)) IYCENT(J) = NINT(SUMYI(J)/SUMI(J)) NPIXCENT(J) = NPIXCENT(J) + ILEN IYBEG(J) = IB IYEND(J) = IE IXMAX(J) = IXCT IYMIN(J) = MIN(IB,IYMIN(J)) IYMAX(J) = MAX(IB,IYMAX(J)) C C---- Flag spot if line overlaps search area limits C IF (IB.EQ.IYL .OR. ENDSTRIPE) IFLAG(J) = 1 C GO TO 50 C C---- End of inner loop C 20 CONTINUE C C---- No active spot for current line - start of new spot C NCENT = NCENT + 1 IF (NCENT.GT.MXCENT) THEN WRITE (6,FMT= + '('' **Too many active spots, Maximum is:'',I6)') + MXCENT STOP END IF C C---- Find empty slot in which to store info about this spot C DO 30 I = 1,MXCENT IF (NPIXCENT(I).EQ.0) THEN IPT(NCENT) = I GO TO 40 END IF 30 CONTINUE C C 40 CONTINUE C C J = IPT(NCENT) NPIXCENT(J) = ILEN SUMXI(J) = TSUMI*FXL SUMYI(J) = TSUMYI SUMI(J) = TSUMI IYBEG(J) = IB IYEND(J) = IE IXCENT(J) = NINT(SUMXI(J)/SUMI(J)) IYCENT(J) = NINT(SUMYI(J)/SUMI(J)) IYMIN(J) = IB IYMAX(J) = IE IXMIN(J) = IXCT IXMAX(J) = IXCT C C---- Flag spot if strip overlaps search area C IF (IB.EQ.IYL .OR. IXCT.EQ.IXL) THEN IFLAG(J) = 1 ELSE IFLAG(J) = 0 END IF C C---- (line.and.itemp.gt.ithresh) ! ongoing line C ELSE FITEMP = REAL(ITEMP-ITHRESH) TSUMYI = REAL(IYPT)*FITEMP + TSUMYI TSUMI = TSUMI + FITEMP END IF C C----! (itemp.le.ithresh) C END IF C C----! (.not.line) C 50 CONTINUE C C---- ! End of loop over pixels C C C C Loop over active spots and store those that have finished C C---- Skip if no active spots C IF (NCENT.EQ.0) GO TO 100 ICENT = 0 60 ICENT = ICENT + 1 J = IPT(ICENT) C C---- Spot is considered to be finished if C a) No pixels for spot found on last stripe C b) and current stripe is at least IXSPLIT from centroid C C Skip if active or within IXSPLIT C IF (IXCT-IXCENT(J).LE.IXSPLIT .OR. IXCT.EQ.IXMAX(J)) GO TO 90 C C---- Reject spot if it contains less than MINPIX pixels C IF (NPIXCENT(J).LT.MINPIX) GO TO 70 C C---- Reject spot if it overlaps limits of search area C IF (IFLAG(J).EQ.1) GO TO 70 FINT = SUMI(J) C C---- Transform to mm C SX = SUMXI(J)*PIXEL/FINT SY = SUMYI(J)*PIXEL/FINT C C---- Store spot if centroid lies within RMM of search box centre C R1 = SQRT((XMM-SX)**2+ (YMM-SY)**2) IF (R1.GT.RMM) GO TO 70 NSPOT = NSPOT + 1 C C type *,'nspot',nspot C IF (NSPOT.GT.MXSPOT) THEN WRITE (6,FMT='('' **Too many spots, Maximum is:'',I6)') MXSPOT STOP END IF C C XSPOT(NSPOT) = SX YSPOT(NSPOT) = SY FSPOT(NSPOT) = FINT RSPOT(NSPOT) = R1 NPIX(NSPOT) = NPIXCENT(J) IWXSPOT(NSPOT) = (IXMAX(J)-IXMIN(J)+1) IWYSPOT(NSPOT) = (IYMAX(J)-IYMIN(J)+1) C C---- Spot stored or rejected C Set flag and shuffle IPT down to fill gap in array C 70 NPIXCENT(J) = 0 NCENT = NCENT - 1 C C---- Skip if last spot C IF (ICENT.GT.NCENT) GO TO 100 C C DO 80 I = ICENT,NCENT IPT(I) = IPT(I+1) 80 CONTINUE C C 90 IF (ICENT.LT.NCENT) GO TO 60 100 CONTINUE C C---- End of loop on stripes/rows/Y C C RETURN END SUBROUTINE GETOD(I,IBA) C ======================== C C C C C---- Get the I'th value for the current (IPOINT'th) C strip from IMAGE C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalar Arguments .. INTEGER I,IBA C .. C .. Scalars in Common .. INTEGER IFILM,IHD,IPOINT,IXLEN,IYLEN,NFILM C .. C .. Arrays in Common .. INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Common blocks .. COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN, + IMAGE(IXLENGTH*IYLENGTH),RDSTRIP(IYLENGTH), + IPOINT,IHD C .. save IBA = IMAGE(IPOINT+I-1) RETURN END SUBROUTINE KEYERR(I,MODE,LINE,IBEG,IEND,ITYP) C ============================================= C C Print warning when token not of correct type. C C Variables used by PARSER C C .. Scalar Arguments .. INTEGER I,MODE CHARACTER LINE*400 C .. C .. Array Arguments .. INTEGER IBEG(200),IEND(200),ITYP(200) C .. C .. Local Arrays .. CHARACTER TYPE(3)*12 C .. C .. Data statements .. DATA TYPE/'alphanumeric','numeric ', + 'quoted '/ C .. C IF (MODE.EQ.0) THEN WRITE (6,FMT=6000) LINE(IBEG(I) :IEND(I)) 6000 FORMAT (' ** ERROR : Key word ',A,' not recognized and has, th', + 'erefore been ignored') C ELSE C WRITE (6,FMT=6002) LINE(IBEG(I) :IEND(I)),TYPE(ITYP(I)),TYPE(I) 6002 FORMAT (' ** ERROR: Token ',A,' is ',A,' while a ',A,' token w', + 'as expected') C END IF RETURN END SUBROUTINE KEYNUM(N,NSTART,LINE,IBEG,IEND,ITYP,NTOK) C ==================================================== C C Check that correct number of numbers are present C C Variables used by PARSER C C C .. Scalar Arguments .. INTEGER N,NSTART,NTOK CHARACTER LINE*400 C .. C .. Array Arguments .. INTEGER IBEG(200),IEND(200),ITYP(200) C .. C .. Local Scalars .. INTEGER I C .. C .. External Subroutines .. EXTERNAL KEYERR C .. DO 10 I = NSTART,NSTART + N - 1 IF (I.GT.NTOK) THEN WRITE (6,FMT=6000) I - NSTART,N 6000 FORMAT (' *** TOO FEW NUMBERS - ',I3,' FOUND WHEN ',I3,' EXP', + 'ECTED') STOP END IF IF (ITYP(I).NE.2) THEN CALL KEYERR(I,2,LINE,IBEG,IEND,ITYP) STOP END IF 10 CONTINUE C RETURN END C C C LOGICAL FUNCTION LITEND(idum) C ========================= C C integer idum C C---- Check endedness, Returns TRUE if little endian (VAX, FX2800, C ULTRIX, CONVEX) C C FALSE if big endian (IBM,IRIS,ESV) C C C .. Local Scalars .. INTEGER I,JDO C .. C .. Local Arrays .. BYTE B(4) C .. C .. Equivalences .. EQUIVALENCE (I,B(1)) C .. C C---- Initialise B C DO 10 JDO = 1,4 B(JDO) = 0 10 CONTINUE C C I = 1 C C IF (B(1).NE.0) THEN LITEND = .TRUE. ELSE LITEND = .FALSE. END IF C C END SUBROUTINE OPENODS(IMFILE) C ========================== C C---- Opens the file of optical densities and positions file C at ods of the first film to be processed C c modified at HH for files from Image plate scanner - c These files are sequential, unformatted, with 16 bit unsigned c integer values, and we read the whole image in here. c images are always in FILM: C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalar Arguments .. CHARACTER IMFILE* (*) C .. C .. Scalars in Common .. INTEGER IFILM,IHD,IPOINT,IXLEN,IYLEN,NFILM C .. C .. Arrays in Common .. INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Local Scalars .. INTEGER I,IBAD,ICURR,II,IJ,INOD,J,NCH,NPIXP1,IYLENM1, + IERR,IBYTES CHARACTER ERRSTR*200,HANDLE*5 cc cc- lines for Bigendian byte swapping BYTE BWORK(2),BTEMP INTEGER*2 IWORK,bloody_extra_2_bytes,AWORK(IYLENGTH) cc- endof mod cc LOGICAL BIGEND,EFILE C .. C .. Local Arrays .. INTEGER IHEAD(IYLENGTH/2) C .. C .. External Functions .. LOGICAL LITEND,VAXVMS INTEGER LENSTR EXTERNAL LITEND,LENSTR,VAXVMS C .. C .. External Subroutines .. EXTERNAL UGERR,CCPERR,QOPEN,QMODE,QREAD,QCLOSE,UBYTES, + RHEAD C .. C .. Common blocks .. COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN, + IMAGE(IXLENGTH*IYLENGTH),RDSTRIP(IYLENGTH), + IPOINT,IHD C .. SAVE C .. C .. Equivalences .. EQUIVALENCE (IWORK,BWORK) C .. C C ICURR = 0 cc cc- lines for Bigendian byte swapping cc BIGEND = .TRUE. C C ********** IF (LITEND()) BIGEND = .FALSE. C ********** C IF(BIGEND) + WRITE(6,*) ' Byte swapping of image will be done' C C WRITE(6,6000) IYLEN,IXLEN 6000 FORMAT(' IYLEN = ',I8, ' IXLEN = ',I8) C C cc- end of Bigendian mod cc C INOD = 4 WRITE (6,FMT=6010) IMFILE(1:LENSTR(IMFILE)) 6010 FORMAT (1X,'IMAGE FILE: ',A) C C ********************** CALL UBYTES (IBYTES,HANDLE) C ********************** C IRECLEN = IYLEN*2 IF (IYLEN.EQ.1186) IRECLEN = 1187*2 IF (HANDLE.EQ.'WORDS') IRECLEN = IRECLEN/2 cc IF (VAXVMS()) THEN OPEN(UNIT=INOD,FILE=IMFILE,STATUS='OLD',FORM='UNFORMATTED', * READONLY,BLOCKSIZE=59904,ERR=9877) ELSE C C ****************************** CALL QOPEN (INOD,IMFILE,'READONLY') CALL QMODE (INOD,1,NCMITM) C ****************************** C END IF cc C C---- Now read the image into the IMAGE array, first initialising C Rdstrip array to false. C DO 10 I = 1,IXLEN RDSTRIP(I) = .FALSE. 10 CONTINUE C C---- Read header if present C IF (IHD.NE.0) THEN IF (VAXVMS()) THEN C C ***************************** CALL RHEAD(INOD,IHD,IHEAD,IYLEN/2) C ***************************** C ELSE C C **************************** CALL QREAD(INOD,IHEAD,IYLEN,IERR) C **************************** C IF (IERR.NE.0) THEN CALL UGERR (-IERR,ERRSTR) CALL CCPERR (1, ERRSTR) ENDIF ENDIF ENDIF C C IBAD = 0 C C---- Mar scanner puts out the image inverted left to right compared with C what MOSFLM expects, so correct for that by reading the image C inverted here C NPIXP1 = IXLEN*IYLEN + 1 IYLENM1 = IYLEN - 1 C C IF (VAXVMS()) THEN C C DO 300 J = 1,IXLEN I = NPIXP1 - J*IYLEN II = I + IYLENM1 C C IF (IYLEN.EQ.1186) THEN READ (INOD) (IMAGE(IJ),IJ=I,II),bloody_extra_2_bytes ELSE READ (INOD) (IMAGE(IJ),IJ=I,II) END IF C C---- check for bad values C DO 200 IJ = I,II IF (IMAGE(IJ).LT.0) THEN C C---- ie 16th bit set C IMAGE(IJ) = 32767 IBAD = IBAD + 1 END IF 200 CONTINUE C C RDSTRIP(J) = .TRUE. 300 CONTINUE C C ELSE C C DO 30 J = 1,IXLEN I = NPIXP1 - J*IYLEN C C ******************************* CALL QREAD(INOD,IMAGE(I),IYLEN,IERR) C ******************************* C IF (IERR.NE.0) THEN C C ************************* CALL UGERR (-IERR,ERRSTR) CALL CCPERR (1, ERRSTR) C ************************* C ENDIF C C IF (IYLEN.EQ.1186) C C *************************************** + CALL QREAD(INOD,bloody_extra_2_bytes,1,IERR) C *************************************** C C---- check for bad values C DO 20 IJ = I,IYLEN+I-1 C C---- Test for BigEndianness Before testing for negative Image(value) C IF (BIGEND) THEN C C---- Ok now swap bytes, as read as an Integer*2 C must equivalence C IWORK = IMAGE(IJ) BTEMP = BWORK(1) BWORK(1) = BWORK(2) BWORK(2) = BTEMP IMAGE(IJ) = IWORK ENDIF C C IF (IMAGE(IJ).LT.0) THEN C C---- ie 16th bit set C IMAGE(IJ) = 32767 IBAD = IBAD + 1 ENDIF 20 CONTINUE RDSTRIP(J) = .TRUE. 30 CONTINUE C C ************* CALL QCLOSE (INOD) C ************* C END IF C C IF (IBAD.GT.0) WRITE (6,FMT=6020) IBAD RETURN C C 9877 CONTINUE WRITE (6,FMT=6002) IMFILE(1:LENSTR(IMFILE)) 6002 FORMAT (//1X,'**** ERROR ****',/1X,'IMAGE FILE ',A,' DOES NOT ', + 'EXIST') STOP C C---- Format Statements C 6020 FORMAT (' This file contains ',i4, ' pixels with value > 32768') C END SUBROUTINE PARSE(LINE,IBEG,IEND,ITYP,FVALUE,IDEC,NTOK) C ====================================================== C C C Finds the tokens in the character string LINE. C Based on Mike Levitt's routine of the same name. C Peter Brick 1-June-1983 C C Modifications C******No***Date*********Purpose**************** C C C Tokens are strings separated by any of the delimiters = , ( ) or C blank. The delimiters are ignored when they occur within a quoted C string. C C If an unquoted exclamation mark (!) occurs in a line then it C and the rest of the line will be skipped. C On entry LINE should contain the line to be parsed C C NTOK number of tokens found in LINE C IBEG(I) start position of the i'th token in LINE C IEND(I) end position of the i'th token in LINE C ITYP(I) = 1 for alphanumeric tokens C = 2 for numeric tokens C = 3 for quoted tokens C FVALUE(I) contains value of numeric tokens C IDEC(I) token length if alphanumeric C if token numeric it contains 100 if decimal point present C + number of places after the decimal point C C Numeric tokens can be simple two argument signed arithmetic C expressions or can contain the E notation (but cannot both contain C two arguments and use the E notation). C Valid examples are: C 10, 10., 1.4-0.6, -1.32*-1.67, +1.45/-1, 1.7E-3 C C .. Scalar Arguments .. INTEGER NTOK CHARACTER LINE* (*) C .. C .. Array Arguments .. REAL FVALUE(40) INTEGER IBEG(40),IDEC(40),IEND(40),ITYP(40) C .. C .. Local Scalars .. REAL F10,SIGN,SIGN0,VALUE,VALUE0 INTEGER I,IDOT,J,LENLIM,NPLACE,OPER LOGICAL BNUM,NUMBER,OPRATR,QUOTE,TOKEN,TQUOTE CHARACTER OLDQUT,LETQT*2,DELIM*7,DIGS*16 C .. C .. Local Arrays .. INTEGER ISGN(2) C .. C .. External Functions .. REAL DOCALC EXTERNAL DOCALC C .. C .. Intrinsic Functions .. INTRINSIC INDEX,LEN C .. C .. Data statements .. DATA ISGN/1,-1/ DATA LETQT/'''"'/ DATA DELIM/' = , ()'/ DATA DIGS/'0123456789+-*/E.'/ C .. C C LENLIM = LEN(LINE) NTOK = 1 TOKEN = .FALSE. VALUE = 0.0 OPRATR = .TRUE. IDOT = 0 SIGN = 1.0 OPER = 0 OLDQUT = ' ' QUOTE = .FALSE. TQUOTE = .FALSE. NUMBER = .FALSE. BNUM = .FALSE. NPLACE = 0 C C*********Start of main character counter loop C DO 30 I = 1,LENLIM C Skip rest of LINE if ! IF (.NOT.QUOTE .AND. LINE(I:I).EQ.'!') GO TO 40 C C Check for quotation mark C J = INDEX(LETQT,LINE(I:I)) IF (J.GT.0) THEN IF (OLDQUT.EQ.' ') THEN OLDQUT = LETQT(J:J) QUOTE = .TRUE. ELSE IF (OLDQUT.EQ.LETQT(J:J)) THEN OLDQUT = ' ' QUOTE = .FALSE. END IF GO TO 30 END IF C C Check for a delimiter C J = INDEX(DELIM,LINE(I:I)) C IF (.NOT. (.NOT.QUOTE) .OR. (J.EQ.0) .AND. + (I.LE.LENLIM)) GO TO 10 C C ie if( quote.OR.( (J.eq.0).and.(I.le.lenlim) ) ) C IF (.NOT. (TOKEN)) GO TO 30 C C We have a token and a delimiter C IF (NUMBER) THEN C ITYP(NTOK) = 2 IEND(NTOK) = I - 1 FVALUE(NTOK) = VALUE*SIGN IF (OPER.GT.0) FVALUE(NTOK) = DOCALC(FVALUE(NTOK),OPER, + SIGN0*VALUE0) IDEC(NTOK) = 100*IDOT + NPLACE C ELSE C IF (TQUOTE) THEN ITYP(NTOK) = 3 IEND(NTOK) = I - 2 ELSE ITYP(NTOK) = 1 IEND(NTOK) = I - 1 END IF IDEC(NTOK) = IEND(NTOK) - IBEG(NTOK) + 1 END IF C NTOK = NTOK + 1 TOKEN = .FALSE. VALUE = 0.0 OPRATR = .TRUE. IDOT = 0 SIGN = 1.0 OPER = 0 TQUOTE = .FALSE. NUMBER = .FALSE. BNUM = .FALSE. GO TO 30 C 10 CONTINUE C C ie if ( quote .OR. (token .AND. (.NOT.number)) ) C IF (.NOT. (.NOT.QUOTE) .OR. (.NOT. (.NOT.TOKEN)) .AND. + ((.NOT. (NUMBER).AND. (.NOT.BNUM)))) GO TO 20 C C Check for digit or +,-,E,/ etc C J = INDEX(DIGS,LINE(I:I)) IF (J.EQ.0) THEN NUMBER = .FALSE. ELSE IF (J.LE.10) THEN NUMBER = .TRUE. BNUM = .FALSE. IF (IDOT.EQ.0) THEN VALUE = (J-1) + VALUE*10 ELSE IF (IDOT.EQ.1) THEN VALUE = (J-1)*F10 + VALUE F10 = F10*0.1 NPLACE = NPLACE + 1 END IF OPRATR = .FALSE. C ELSE IF (OPRATR .AND. ((J.EQ.11).OR. (J.EQ.12))) THEN C C Opratr is true and we have a + or - C BNUM = .TRUE. OPRATR = .FALSE. SIGN = ISGN(J-10) C ELSE IF (J.GE.11 .AND. J.LE.15) THEN C C We have an operator C NUMBER = .TRUE. BNUM = .FALSE. IF (OPRATR .OR. (OPER.NE.0)) NUMBER = .FALSE. C C - store 1st number C VALUE0 = VALUE SIGN0 = SIGN OPER = J - 10 VALUE = 0.0 SIGN = 1.0 IDOT = 0 OPRATR = .TRUE. C ELSE IF (J.EQ.16) THEN C C We have a decimal point C NUMBER = .TRUE. BNUM = .FALSE. IDOT = IDOT + 1 NPLACE = 0 F10 = 0.1 IF (IDOT.EQ.2) NUMBER = .FALSE. OPRATR = .FALSE. C END IF C 20 CONTINUE IF (.NOT.TOKEN) THEN TOKEN = .TRUE. IBEG(NTOK) = I IF (QUOTE) TQUOTE = .TRUE. END IF C 30 CONTINUE C 40 CONTINUE C C C End of line acts as a terminator C C IF (TOKEN .AND. NUMBER) THEN C ITYP(NTOK) = 2 FVALUE(NTOK) = VALUE*SIGN IF (OPER.GT.0) FVALUE(NTOK) = DOCALC(FVALUE(NTOK),OPER, + SIGN0*VALUE0) IDEC(NTOK) = 100*IDOT + NPLACE C ELSE IF (TOKEN .AND. (TQUOTE.AND. (.NOT.QUOTE))) THEN ITYP(NTOK) = 3 IEND(NTOK) = LENLIM - 1 IDEC(NTOK) = IEND(NTOK) - IBEG(NTOK) + 1 ELSE IF (TOKEN .AND. (.NOT.TQUOTE)) THEN ITYP(NTOK) = 1 IEND(NTOK) = LENLIM IDEC(NTOK) = IEND(NTOK) - IBEG(NTOK) + 1 ELSE NTOK = NTOK - 1 END IF RETURN END SUBROUTINE PARSER(INFILE,IPRINT,LINEX,IBEGX,IENDX,ITYPX,VALUEX, + IDECX,NTOKX) C =========================================================== C C This routine reads in lines from unit INFILE, reflects them to C unit number IPRINT and parses the lines by calling PARSE. C The lines can be continued by using either a minus sign "-" C or an ampersand "&" as a token and then continuing the record on the C following line. C See PARSE for a description of the subroutine arguments. C C C The following 3 lines needed for this routine C C C Following 3 lines needed for the PARSE subroutine. C C C .. Scalar Arguments .. INTEGER INFILE,IPRINT,NTOKX CHARACTER LINEX*400 C .. C .. Array Arguments .. REAL VALUEX(200) INTEGER IBEGX(200),IDECX(200),IENDX(200),ITYPX(200) C .. C .. Local Scalars .. INTEGER ICX,ILAST,ILINX,INCR,ITOK,MAXLIN,NTOK LOGICAL CONTIN CHARACTER LINE*80 C .. C .. Local Arrays .. REAL VALUE(40) INTEGER IBEG(40),IDEC(40),IEND(40),ITYP(40) C .. C .. External Subroutines .. EXTERNAL PARSE C .. C .. Data statements .. DATA MAXLIN/5/ C .. C NTOKX = 0 ILINX = 1 INCR = 0 ICX = 0 C 10 READ (INFILE,FMT=6000,END=30) LINE 6000 FORMAT (A) C WRITE (IPRINT,FMT=6002) LINE 6002 FORMAT (' ===> ',A) C CALL PARSE(LINE,IBEG,IEND,ITYP,VALUE,IDEC,NTOK) C IF (NTOK.EQ.0) GO TO 10 C C--Check if this LINE is to be continued C via "-" or "&" as last token in line C IF (LINE(IBEG(NTOK) :IEND(NTOK)).EQ.'-' .OR. + LINE(IBEG(NTOK) :IEND(NTOK)).EQ.'&') THEN CONTIN = .TRUE. NTOK = NTOK - 1 ELSE CONTIN = .FALSE. END IF C C--Copy line arrays to record arrays C DO 20 ITOK = 1,NTOK NTOKX = NTOKX + 1 IBEGX(NTOKX) = IBEG(ITOK) + INCR IENDX(NTOKX) = IEND(ITOK) + INCR ITYPX(NTOKX) = ITYP(ITOK) VALUEX(NTOKX) = VALUE(ITOK) IDECX(NTOKX) = IDEC(ITOK) 20 CONTINUE C C--Concatenate lines C IF (CONTIN) THEN ILAST = IEND(NTOK+1) - 1 ELSE ILAST = IEND(NTOK) END IF C IF (ICX.EQ.0) THEN LINEX = LINE(1:ILAST) ELSE LINEX = LINEX(:ICX)//LINE(1:ILAST) END IF ICX = ICX + ILAST C IF (.NOT.CONTIN) RETURN C ILINX = ILINX + 1 INCR = IEND(NTOK+1) + INCR - 1 IF (ILINX.LE.MAXLIN) GO TO 10 C C-- Too many lines C WRITE (IPRINT,FMT=6004) 6004 FORMAT (' --> *** TOO MANY CONTINUED LINES - MAXIMUM OF 5 ***') STOP C 30 CONTINUE NTOKX = -1 RETURN END SUBROUTINE PICKSPOTS C ===================== C C C C---- This is the subroutine that controls the spot search. C The coordinates of the spots on the still are written C out to a file which can then be read in by IDXREF. C C C---- Maximum number spots C C C C Local variables C C .. Parameters .. INTEGER MXSPOT PARAMETER (MXSPOT=5000) INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN,CUTWYMAX,CUTWYMIN, + PIXEL,RMAX,RMIN,SPSIGFACT,XC,XSPLIT,YC,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,NBINR,NBINT,NFILM, + NPIXMAX,NPIXMIN LOGICAL LIST,PRINT C .. C .. Arrays in Common .. REAL FSPOT,PHI,RSPOT,XSPOT,YSPOT INTEGER BGOD,IORDER,IWXSPOT,IWYSPOT,NPIX INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Local Scalars .. REAL FRACT,PI,RMM,RSPOTMAX,THETA,WBINR2,WBINT,WXMAX,WXMIN,WXMM, + WYMAX,WYMIN,WYMM,XF,XF1,XMM,XS,YF,YF1,YMM,YS INTEGER I,IFRIG,IR,IT,ITHRESH,J,MEDWXSPOT,MEDWYSPOT,MINPIX, + NREJBIGX,NREJBIGY,NREJSMALLX,NREJSMALLY,NSPOT,NWRITE LOGICAL FID C .. C .. Local Arrays .. INTEGER ICOUNT(16,5) C .. C .. External Subroutines .. EXTERNAL FINDSPOTS,SORTUP,SORTUPINT C .. C .. Intrinsic Functions .. INTRINSIC ATAN,ATAN2,INT,MAX,MIN C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /COORDS/XC,YC,PIXEL COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN,IMAGE(IXLENGTH*IYLENGTH), + RDSTRIP(IYLENGTH),IPOINT,IHD COMMON /SCRATCH/XSPOT(MXSPOT),YSPOT(MXSPOT),RSPOT(MXSPOT), + FSPOT(MXSPOT),NPIX(MXSPOT),IWXSPOT(MXSPOT),IWYSPOT(MXSPOT), + IORDER(MXSPOT) COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX,NBINR,NBINT, + IBINLIM,PHI(3),NPIXMIN,NPIXMAX,CUTPIXMIN,CUTPIXMAX, + CUTWXMIN,CUTWXMAX,CUTWYMIN,CUTWYMAX,LIST C .. save C .. Data statements .. DATA FRACT/0.5/ C .. C C FID = .FALSE. PI = ATAN(1.0)*4.0 C C---- Zero counter for spots in each bin C DO 20 I = 1,NBINR DO 10 J = 1,NBINT ICOUNT(J,I) = 0 10 CONTINUE 20 CONTINUE C C---- Determine orientation of film on drum of scanner. C C If X coord of fid 1 is positive then film has been placed on drum C in the same orientation as film was on camera and we need to frig the C axes to get coords into the standard frame used by the IDXREF etc. C C IF(XYFID(1,1).GT.0.0) THEN C IFRIG = 1 C ELSE C IFRIG = 0 C ENDIF C IFRIG = 1 C RMM = RMAX WXMM = XSPLIT WYMM = YSPLIT C if (fixed) ithresh = bgody(-1) + spsigfact * bgsigy(-1) XMM = XC YMM = YC MINPIX = NPIXMIN C C *********************************************************** CALL FINDSPOTS(XMM,YMM,RMM,ITHRESH,WXMM,WYMM,MINPIX,NSPOT,XSPOT, + YSPOT,RSPOT,FSPOT,NPIX,IWXSPOT,IWYSPOT,FID) C *********************************************************** C IF (NSPOT.EQ.0) THEN WRITE (6,FMT='('' ** No spots found on still!!'')') STOP END IF C C---- Open output SPOT file and write the camera constants and C fiducial coordinates. All coordinates in mms C IF (IFILM.EQ.1) THEN LDUMK = 80 IFAILK = 0 CALL CCPDPN(2,'SPOTS','NEW','F',LDUMK,IFAILK) cc-vms open (2,file='spots',status='new',carriagecontrol='list') END IF C C---- Write out "pseudo fiducials" for image plates C First the camera constants (zero), then give fiducials 1 and 3 C the direct beam coordinates as specified by XC,YC C XF = 0.0 YF = 0.0 WRITE (2,FMT='(1X,2F10.3)') XF,YF C C----- Note that for the direct beam coordinates we C need to transform the input XC,YC to the MOSFLM frame C XF1 = YC YF1 = 2*XC - XC WRITE (2,FMT='(1X,2F10.3)') XF1,YF1 WRITE (2,FMT='(1X,2F10.3)') XF,YF1 WRITE (2,FMT='(1X,2F10.3)') XF1,YF1 C C Camera constants C C IF(XBEAM.NE.0.0.AND.YBEAM.NE.0.0) THEN C CCX = XBEAM - XC C CCY = YBEAM - YC C IF(IFRIG.EQ.1) THEN C XF = CCY C YF = -CCX C ELSE C XF = CCX C YF = CCY C ENDIF C ELSE C XF = 0.0 C YF = 0.0 C ENDIF C C WRITE(2,FMT='(1X,2F10.3)') XF,YF C C Fiducial coords C C DO 80 J = 1,3 C IF(IFRIG.EQ.1) THEN C XF = FPOSN(2,J) C YF = -FPOSN(1,J) + 2.0*XC C ELSE C XF = FPOSN(1,J) C YF = FPOSN(2,J) C ENDIF C WRITE(2,FMT='(1X,2F10.3)') XF,YF C 80 CONTINUE C C Test for acceptable spots C C We try to reject spots from bits of dirt, hair, streaks etc. C To do this we find the median widths and then reject those C spots that are either too small or too big relative to the C median spot size. C C---- Sort spots on width in x C C ******************************** CALL SORTUPINT(NSPOT,IWXSPOT,IORDER) C ******************************** C C---- Determine median x width and use it to set cutoffs C MEDWXSPOT = IWXSPOT(IORDER(INT(NSPOT*0.6))) C C---- Sort spots on width in Y C C ******************************** CALL SORTUPINT(NSPOT,IWYSPOT,IORDER) C ******************************** C C---- Determine median Y width and use it to set cutoffs C MEDWYSPOT = IWYSPOT(IORDER(INT(NSPOT*0.6))) C WRITE (6,FMT=6000) MEDWXSPOT,MEDWYSPOT 6000 FORMAT (//1X,'Median spot size is',I4,' x',I4,' pixels') C WXMIN = MEDWXSPOT*CUTWXMIN WXMAX = MEDWXSPOT*CUTWXMAX C WYMIN = MEDWYSPOT*CUTWYMIN WYMAX = MEDWYSPOT*CUTWYMAX C NREJBIGX = 0 NREJSMALLX = 0 C NREJBIGY = 0 NREJSMALLY = 0 C C---- Sort all spots on intensity C C ************************** CALL SORTUP(NSPOT,FSPOT,IORDER) C ************************** C C---- Find actual maximum radius for binning purposes C RSPOTMAX = 0.0 C C DO 30 I = 1,NSPOT IF (IWXSPOT(I).LT.WXMIN) GO TO 30 IF (IWXSPOT(I).GT.WXMAX) GO TO 30 IF (IWYSPOT(I).LT.WYMIN) GO TO 30 IF (IWYSPOT(I).GT.WYMAX) GO TO 30 RSPOTMAX = MAX(RSPOTMAX,RSPOT(I)) 30 CONTINUE C C---- Calculate widths of bins C WBINR2 = (RSPOTMAX**2)/NBINR WBINT = 2*PI/NBINT NWRITE = 0 IF (LIST) WRITE (6,FMT=6002) 6002 FORMAT (/6X,'X',7X,'Y',8X,'R',4X,'Intensity',2X,'Pixels WX ', + ' WY') C C---- Loop over spots starting with the strongest C Write out the strongest IBINLIM spots in each bin C DO 40 I = NSPOT,1,-1 J = IORDER(I) IF (RSPOT(J).LT.RMIN) GO TO 40 IF (IWXSPOT(J).LT.WXMIN) THEN NREJSMALLX = NREJSMALLX + 1 IF (LIST) WRITE (6,FMT=6004) (XSPOT(J)-XC), (YSPOT(J)-YC), + RSPOT(J),FSPOT(J),NPIX(J),IWXSPOT(J),IWYSPOT(J) 6004 FORMAT (1X,2F8.2,2X,F7.2,2X,F7.0,2X,I6,2X,2I6,' REJECTED') GO TO 40 ELSE IF (IWYSPOT(J).LT.WYMIN) THEN NREJSMALLY = NREJSMALLY + 1 IF (LIST) WRITE (6,FMT=6004) (XSPOT(J)-XC), (YSPOT(J)-YC), + RSPOT(J),FSPOT(J),NPIX(J),IWXSPOT(J),IWYSPOT(J) GO TO 40 ELSE IF (IWXSPOT(J).GT.WXMAX) THEN NREJBIGX = NREJBIGX + 1 IF (LIST) WRITE (6,FMT=6004) (XSPOT(J)-XC), (YSPOT(J)-YC), + RSPOT(J),FSPOT(J),NPIX(J),IWXSPOT(J),IWYSPOT(J) GO TO 40 ELSE IF (IWYSPOT(J).GT.WYMAX) THEN NREJBIGY = NREJBIGY + 1 IF (LIST) WRITE (6,FMT=6004) (XSPOT(J)-XC), (YSPOT(J)-YC), + RSPOT(J),FSPOT(J),NPIX(J),IWXSPOT(J),IWYSPOT(J) GO TO 40 END IF C C---- Find radial bin C IR = MIN(NBINR,INT(RSPOT(J)**2/WBINR2+1.0)) C C---- Find theta bin C XS = XSPOT(J) - XC YS = YSPOT(J) - YC THETA = ATAN2(YS,XS) IT = MIN(NBINT,INT((THETA+PI)/WBINT+1.0)) C C---- Increment counter C ICOUNT(IT,IR) = ICOUNT(IT,IR) + 1 C IF (ICOUNT(IT,IR).LE.IBINLIM) THEN NWRITE = NWRITE + 1 IF (IFRIG.EQ.1) THEN XF = YSPOT(J) YF = -XSPOT(J) + 2.0*XC ELSE XF = XSPOT(J) YF = YSPOT(J) END IF C WRITE (2,FMT='(1X,4F8.3,F10.1)') XF,YF,FRACT,PHI(IFILM), + FSPOT(J) C IF (LIST) WRITE (6,FMT=6006) (XSPOT(J)-XC), (YSPOT(J)-YC), + RSPOT(J),FSPOT(J),NPIX(J),IWXSPOT(J),IWYSPOT(J) 6006 FORMAT (1X,2F8.2,2X,F7.2,2X,F7.0,2X,I6,2X,2I6) END IF 40 CONTINUE C WRITE (6,FMT=6008) NSPOT,NREJSMALLX,NREJSMALLY,NREJBIGX,NREJBIGY, + NWRITE 6008 FORMAT (/2X, + 'Total spots found:',I4,/2X, + 'Number rejected as too small on X:',I4,/2X, + 'Number rejected as too small on Y:',I4,/2X, + 'Number rejected as too big on X:',I4,/2X, + 'Number rejected as too big on Y:',I4,/2X, + 'Number written to file:',I4) C WRITE (6,FMT=6010) 6010 FORMAT (/13X,'Number of spots found in each bin') C C DO 50 J = 1,NBINR WRITE (6,FMT='(12X,12I4)') (ICOUNT(I,J),I=1,NBINT) 50 CONTINUE C C RETURN END SUBROUTINE PRINTDAT C ==================== C C C C---- Print out values used by program C C C C C---- Coords in mm of film center wrt bottom left corner C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalars in Common .. REAL CUTPIXMAX,CUTPIXMIN,CUTWXMAX,CUTWXMIN,CUTWYMAX,CUTWYMIN, + PIXEL,RMAX,RMIN,SPSIGFACT,XC,XSPLIT,YC,YSPLIT,YOFF INTEGER IBINLIM,IFILM,IHD,IPOINT,IXLEN,IYLEN,NBINR,NBINT,NFILM, + NPIXMAX,NPIXMIN LOGICAL LIST,PRINT C .. C .. Arrays in Common .. REAL PHI INTEGER BGOD INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Common blocks .. COMMON /BACKGR/BGOD(800),PRINT,YOFF COMMON /COORDS/XC,YC,PIXEL COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN,IMAGE(IXLENGTH*IYLENGTH), + RDSTRIP(IYLENGTH),IPOINT,IHD COMMON /SPOTS/XSPLIT,YSPLIT,SPSIGFACT,RMIN,RMAX,NBINR,NBINT, + IBINLIM,PHI(3),NPIXMIN,NPIXMAX,CUTPIXMIN,CUTPIXMAX, + CUTWXMIN,CUTWXMAX,CUTWYMIN,CUTWYMAX,LIST C .. save WRITE (6,FMT=6000) XC,YC 6000 FORMAT (/1X,'Coords in mm of image center ',2F7.2/,1X,'These ', + 'should be wrt the lower left hand corner of the image', + /,1X,'Note that the display programs display the image ', + 'looking from the source',/,1X,'towards the detector', + /,1X,'MOSFLM works with a view towards the source, so', + ' the X coordinate of the',/,1X,'beam centre (found from', + ' the display program or WAX) is changed internally to', + /,1X,'(width of image - input coordinate) to allow for', + ' this') XC = (IXLEN*PIXEL - XC) C C---- Scanner-increment/pixel size in mms C WRITE (6,FMT=6002) PIXEL 6002 FORMAT (1X,'Scanner pixel size in mm:',F7.4) C C---- Number of pixels in the X (along drum) C and Y (around drum) direction C WRITE (6,FMT=6004) IXLEN,IYLEN 6004 FORMAT (/1X,'Number of pixels in the X and Y directions:',2I5) C C---- Now for numbers that control the spot search procedure C Spot size in X and Y directions, C threshold (in sigmas) for spot C WRITE (6,FMT=6006) XSPLIT,YSPLIT,SPSIGFACT 6006 FORMAT (/1X,'Split spot distance in X and Y:',2F5.2,/1X, + 'Standard deviations above background for spot threshold OD:', + F5.2) C C---- Radii limits on image for spots, C number of bins across or down image and C maximum number of spots in a bin C WRITE (6,FMT=6008) RMIN,RMAX,NBINR,NBINT,IBINLIM,NPIXMIN,NPIXMAX, + CUTPIXMIN,CUTPIXMAX 6008 FORMAT (1X, + 'Minimum and maximum radii for spot search:',2F7.2,/1X, + 'Number of annuli and angular subdivisions for spot selection:', + 2I5,/1X, + 'Maximum number of spots in a bin:',I4,/1X, + 'Minimum and maximum pixels in a spot:',2I6,/1X, + 'Factors specifying minimum and maximum spot sizes',2X, + 'as a function',/1X,'of the median size:',2F7.2) C RETURN END SUBROUTINE RDBLK(IBLK) C ====================== C C C---- Check if record IBLK from image file has already been read into C IMAGE, and if not, read it. Length of record comes from /SCN/. C In Image plate version length is NWORD, not NBYTE C C C Note that if all image is not stored in memory (INCORE=FALSE) C C C C .. Parameters .. INTEGER IXLENGTH PARAMETER (IXLENGTH=1200) INTEGER IYLENGTH PARAMETER (IYLENGTH=1200) C .. C .. Scalar Arguments .. INTEGER IBLK C .. C .. Scalars in Common .. INTEGER IFILM,IHD,IPOINT,IXLEN,IYLEN,NFILM C .. C .. Arrays in Common .. INTEGER*2 IMAGE LOGICAL RDSTRIP C .. C .. Common blocks .. COMMON /ODREC/NFILM,IFILM,IXLEN,IYLEN, + IMAGE(IXLENGTH*IYLENGTH),RDSTRIP(IYLENGTH), + IPOINT,IHD C .. save IPOINT = (IBLK-1)*IYLEN + 1 RETURN END C C C C ======================================= SUBROUTINE RHEAD(INOD,IHD,IHEAD,IYLEN2) C ======================================= C C C C----- This will require a change for BigEndian machines C C .. Scalar Arguments .. INTEGER IHD,INOD,IYLEN2 C .. C .. Array Arguments .. INTEGER IHEAD(IYLEN2) C .. C .. Local Scalars .. INTEGER I C .. DO 10 I = 1,IHD READ (INOD) IHEAD 10 CONTINUE RETURN END C C SUBROUTINE GETOD(I) C C COMMON /STRIPE/ BVAL(4),BBA(2560) C BYTE BVAL,BBA C C BVAL(1)=BBA(I) C RETURN C END C C SUBROUTINE SORTUP(N,A,IN) C ========================= C C---- Ref: Comm. ACM VOL.12 #3 MARCH 1969, R.C.SINGLETON C C Routine returns order of A in IN C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL A(N) INTEGER IN(N) C .. C .. Local Scalars .. INTEGER I,IJ,J,K,L,M,T,TT C .. C .. Local Arrays .. INTEGER IL(16),IU(16) C .. DO 10 I = 1,N IN(I) = I 10 CONTINUE C M = 1 I = 1 J = N C 20 IF (I.GE.J) GO TO 100 30 K = I IJ = (I+J)/2 T = IN(IJ) IF (A(IN(I)).LE.A(T)) GO TO 40 IN(IJ) = IN(I) IN(I) = T T = IN(IJ) 40 L = J IF (A(IN(J)).GE.A(T)) GO TO 70 IF (A(IN(J)).LT.A(IN(I))) GO TO 50 IN(IJ) = IN(J) IN(J) = T T = IN(IJ) GO TO 70 C 50 IN(IJ) = IN(I) IN(I) = IN(J) IN(J) = T T = IN(IJ) GO TO 70 C 60 IN(L) = IN(K) IN(K) = TT 70 L = L - 1 IF (A(IN(L)).GT.A(T)) GO TO 70 TT = IN(L) 80 K = K + 1 IF (A(IN(K)).LT.A(T)) GO TO 80 IF (K.LE.L) GO TO 60 IF ((L-I).LE. (J-K)) GO TO 90 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 110 90 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 110 C 100 M = M - 1 IF (M.EQ.0) GO TO 140 I = IL(M) J = IU(M) 110 IF ((J-I).GE.11) GO TO 30 IF (I.EQ.1) GO TO 20 I = I - 1 120 I = I + 1 IF (I.EQ.J) GO TO 100 T = IN(I+1) IF (A(IN(I)).LE.A(T)) GO TO 120 K = I 130 IN(K+1) = IN(K) K = K - 1 IF (A(T).LT.A(IN(K))) GO TO 130 IN(K+1) = T GO TO 120 C 140 RETURN END SUBROUTINE SORTUP4(N,A,IN) C ========================== C C---- Ref: Comm. ACM VOL.12 #3 MARCH 1969, R.C.SINGLETON C C---- Routine returns order of A in IN C C C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. INTEGER A(N),IN(N) C .. C .. Local Scalars .. INTEGER I,IJ,J,K,L,M INTEGER*2 T,TT C .. C .. Local Arrays .. INTEGER*2 IL(16),IU(16) C .. C C DO 10 I = 1,N IN(I) = I 10 CONTINUE C C M = 1 I = 1 J = N C 20 IF (I.GE.J) GO TO 80 30 K = I IJ = (I+J)/2 T = IN(IJ) C C IF (A(IN(I)).GT.A(T)) THEN IN(IJ) = IN(I) IN(I) = T T = IN(IJ) END IF C C L = J IF (A(IN(J)).GE.A(T)) GO TO 50 C C IF (A(IN(J)).GE.A(IN(I))) THEN IN(IJ) = IN(J) IN(J) = T T = IN(IJ) GO TO 50 END IF C IN(IJ) = IN(I) IN(I) = IN(J) IN(J) = T T = IN(IJ) GO TO 50 C 40 IN(L) = IN(K) IN(K) = TT 50 L = L - 1 IF (A(IN(L)).GT.A(T)) GO TO 50 TT = IN(L) 60 K = K + 1 IF (A(IN(K)).LT.A(T)) GO TO 60 IF (K.LE.L) GO TO 40 IF ((L-I).LE. (J-K)) GO TO 70 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 90 70 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 90 C 80 M = M - 1 IF (M.EQ.0) RETURN I = IL(M) J = IU(M) 90 IF ((J-I).GE.11) GO TO 30 IF (I.EQ.1) GO TO 20 I = I - 1 100 I = I + 1 IF (I.EQ.J) GO TO 80 T = IN(I+1) IF (A(IN(I)).LE.A(T)) GO TO 100 K = I 110 IN(K+1) = IN(K) K = K - 1 IF (A(T).LT.A(IN(K))) GO TO 110 IN(K+1) = T GO TO 100 C C C END SUBROUTINE SORTUPINT(N,A,IN) C ============================= C C C Ref: Comm. ACM VOL.12 #3 MARCH 1969, R.C.SINGLETON C C Routine returns order of A in IN C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. INTEGER A(N),IN(N) C .. C .. Local Scalars .. INTEGER I,IJ,J,K,L,M,T,TT C .. C .. Local Arrays .. INTEGER IL(16),IU(16) C .. DO 10 I = 1,N IN(I) = I 10 CONTINUE C M = 1 I = 1 J = N C 20 IF (I.GE.J) GO TO 100 30 K = I IJ = (I+J)/2 T = IN(IJ) IF (A(IN(I)).LE.A(T)) GO TO 40 IN(IJ) = IN(I) IN(I) = T T = IN(IJ) 40 L = J IF (A(IN(J)).GE.A(T)) GO TO 70 IF (A(IN(J)).LT.A(IN(I))) GO TO 50 IN(IJ) = IN(J) IN(J) = T T = IN(IJ) GO TO 70 C 50 IN(IJ) = IN(I) IN(I) = IN(J) IN(J) = T T = IN(IJ) GO TO 70 C 60 IN(L) = IN(K) IN(K) = TT 70 L = L - 1 IF (A(IN(L)).GT.A(T)) GO TO 70 TT = IN(L) 80 K = K + 1 IF (A(IN(K)).LT.A(T)) GO TO 80 IF (K.LE.L) GO TO 60 IF ((L-I).LE. (J-K)) GO TO 90 IL(M) = I IU(M) = L I = K M = M + 1 GO TO 110 90 IL(M) = K IU(M) = J J = L M = M + 1 GO TO 110 C 100 M = M - 1 IF (M.EQ.0) GO TO 140 I = IL(M) J = IU(M) 110 IF ((J-I).GE.11) GO TO 30 IF (I.EQ.1) GO TO 20 I = I - 1 120 I = I + 1 IF (I.EQ.J) GO TO 100 T = IN(I+1) IF (A(IN(I)).LE.A(T)) GO TO 120 K = I 130 IN(K+1) = IN(K) K = K - 1 IF (A(T).LT.A(IN(K))) GO TO 130 IN(K+1) = T GO TO 120 C 140 RETURN END