**==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 C C where: C 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 name of crystal system - only first three letters C are used, e.g. tri for triclinic C C a single letter for the type of unit cell. Default P C C 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 similar form to 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