**==aa0001.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** Program DECONV ** C*************************** C C Purpose: To deconvolute intensities for harmonics C C Input: .ge files, .vc file (film pack only), and .mtz file C ( .ge files, i.e. .ge1 and .ge2 are produced by INTLAUE C .vc (Victoreen Coefficient) files produced by AFSCALE) C .mtz file produced by SFALL) C INCLUDE 'implicit_none.finc' C INCLUDE 'overlaps.finc' C COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10) INTEGER NRANGE,ORDER REAL WAV1,WAV2,R COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C C Var to select box or profile intensities COMMON /BP/ BOXPROF INTEGER BOXPROF C C====== FUNCTIONS C INTEGER LUNSTI, LUNSTO C C====== LOCALS: C DOUBLE PRECISION PARAMS(2,3,5) INTEGER END_REC, I, END_POS INTEGER IFAIL REAL ACELL, BCELL, CCELL, CELLANG(3) REAL ITOP, IBOT REAL N_SIGS REAL WAVMIN, WAVMAX INTEGER NUM_OVERLAPS REAL OVERLAP_INTENS(6,MAX_OVERLAPS) REAL OVERLAP_SIGMAS(6,MAX_OVERLAPS) REAL OVERLAP_LAMBDAS(MAX_OVERLAPS) REAL OVERLAP_DHKLS(MAX_OVERLAPS) REAL OVERLAP_XY(2,MAX_OVERLAPS) INTEGER OVERLAP_HKLS(MAX_OVERLAPS) LOGICAL QUIT C C====== Set sensible values C DATA CELLANG /90., 90., 90./ DATA ITOP, IBOT /9000.0,100.0/ DATA N_SIGS /5/ DATA WAVMIN, WAVMAX /0.3,2.5/ C C------------------------------------------------------------------------------ C CALL TRWSET(LUNSTI(0),LUNSTO(0),8,8,0) IN_TERM = LUNSTI(0) OUT_TERM = LUNSTO(0) IN_DATA = 2 OUT_DATA = 1 IFAIL = 0 CALL CCPDPN(-8,'DECLOG','NEW','F',0,IFAIL) CALL TWRITE(' ') CALL TWRITE('Program: DECONV') CALL TWRITE('Last update: 01/07/93') CALL TWRITE('Authors: Q.Hao') CALL TWRITE(' ') C C====== Read in the Victoreen parameters C CALL READ_VC(PARAMS) C C====== Open .ge files and find last record C CALL OPEN_GE(IN_DATA,3,END_POS,END_REC) C C====== Get rejection criteria and unit cell parameters C CALL GET_REJECT_CRIT(ACELL,BCELL,CCELL,CELLANG,ITOP,IBOT,N_SIGS, + WAVMIN,WAVMAX) C C====== Use box or profile intensities? C BOXPROF = 2 30 CALL TASKI('Use box (1) or profile fitted (2) intensities', + BOXPROF,QUIT) IF (BOXPROF.LE.0 .OR. BOXPROF.GT.2) GO TO 30 C C====== Now deconvolute all the spots C CALL DECONV_SPOTS(MAX_OVERLAPS,END_REC,ACELL,BCELL,CCELL,CELLANG, + ITOP,IBOT,N_SIGS,WAVMIN,WAVMAX,PARAMS, + NUM_OVERLAPS,OVERLAP_INTENS,OVERLAP_SIGMAS, + OVERLAP_LAMBDAS,OVERLAP_HKLS,OVERLAP_DHKLS, + OVERLAP_XY) C IF (NUM_OVERLAPS.GT.0) THEN CALL OUTPUT_DECONV(NUM_OVERLAPS,OVERLAP_INTENS, + OVERLAP_LAMBDAS,OVERLAP_HKLS,OVERLAP_DHKLS, + OVERLAP_SIGMAS,OVERLAP_XY,PARAMS,ACELL, + BCELL,CCELL,CELLANG) ELSE STOP ' No overlaps found' END IF C STOP END **==apply_reject_crit.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************************ C** APPLY_REJECT_CRIT ** C************************************ C C INTEGER FUNCTION APPLY_REJECT_CRIT ( MULTIPLICITY, + INTENSITIES, SIGMAS, LAMBDAS, SPATIAL_OVERLAP, + ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX, + NUM_GOOD_FILMS ) C C Purpose: Returns 1 if spot is rejected, 0 if OK C INCLUDE 'implicit_none.finc' C INTEGER MULTIPLICITY REAL INTENSITIES(6), SIGMAS(6), LAMBDAS(MULTIPLICITY) INTEGER SPATIAL_OVERLAP REAL ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX INTEGER NUM_GOOD_FILMS C C ITOP : maximum intensity allowed C IBOT : minimum intensity allowed C 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 INTEGER I REAL WAV C C------------------------------------------------------------------------------ C APPLY_REJECT_CRIT = 0 C IF (SPATIAL_OVERLAP.EQ.0) THEN C C====== Check intensities C NUM_GOOD_FILMS = 6 DO 20 I = 1 , 6 IF (INTENSITIES(I).GT.ITOP .OR. INTENSITIES(I).LT.IBOT .OR. + INTENSITIES(I).LT.N_SIGS*ABS(SIGMAS(I))) THEN INTENSITIES(I) = -1.0*ABS(INTENSITIES(I)) NUM_GOOD_FILMS = NUM_GOOD_FILMS - 1 END IF 20 CONTINUE C C====== Check wavelengths of harmonics C DO 30 I = 1 , MULTIPLICITY C WAV = LAMBDAS(I) C IF (WAV.LT.WAVMIN .OR. WAV.GT.WAVMAX) GO TO 40 C IF ((WAV.GE.0.486 .AND. WAV.LE.0.490) .OR. + (WAV.GE.0.900 .AND. WAV.LE.0.920)) GO TO 40 C 30 CONTINUE C RETURN END IF C C====== Failed test C 40 APPLY_REJECT_CRIT = 1 RETURN END **==get_ge_data.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C****************************** C** GET_GE_DATA ** C****************************** C C INTEGER FUNCTION GET_GE_DATA (IN_DATA1, IN_DATA2, RECORD_NUM, + X, Y, INTEN_BOX_S, SIG_BOX_S, INTEN_PROF_S, SIG_PROF_S, + MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP, + H_S, K_S, L_S, LAMBDAS, DMIN_THRESHOLD ) C C Purpose: To read in data from a .ge format file C INCLUDE 'implicit_none.finc' C INTEGER IN_DATA1, IN_DATA2, RECORD_NUM REAL X, Y, INTEN_BOX_S(6), SIG_BOX_S(6), + INTEN_PROF_S(6), SIG_PROF_S(6) INTEGER MULTIPLICITY, NODAL, NODIDX, SPATIAL_OVERLAP INTEGER H_S(MULTIPLICITY), K_S(MULTIPLICITY), L_S(MULTIPLICITY) REAL LAMBDAS(MULTIPLICITY) REAL DMIN_THRESHOLD C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C INTEGER NBITST C C====== LOCALS: C INTEGER I, REC_NUM, LINK_REC, MAX_MULTI REAL BX, BY, BDMIN INTEGER BINTENBS(6), BSIGBS(6), BINTENPS(6), BSIGPS(6) REAL LAMBDA, LAMBDA1, LAMBDA2 INTEGER M_FLAGS, H, K, L, NIDX, H2, K2, L2 INTEGER IREAL REAL RINT CHARACTER*80 STR C INTEGER BUFFER(19), BUFFER2(2) INTEGER MHKL(4), HKL2(4) C EQUIVALENCE (IREAL, RINT) C EQUIVALENCE (BUFFER(1),BX) EQUIVALENCE (BUFFER(2),BY) EQUIVALENCE (BUFFER(16),LAMBDA1) EQUIVALENCE (BUFFER(18),LAMBDA2) EQUIVALENCE (BUFFER(19),BDMIN) C EQUIVALENCE (MHKL(1),M_FLAGS) EQUIVALENCE (MHKL(2),H) EQUIVALENCE (MHKL(3),K) EQUIVALENCE (MHKL(4),L) C EQUIVALENCE (HKL2(1),NIDX) EQUIVALENCE (HKL2(2),H2) EQUIVALENCE (HKL2(3),K2) EQUIVALENCE (HKL2(4),L2) C EQUIVALENCE (BUFFER2(2),LAMBDA) C C------------------------------------------------------------------------------ C GET_GE_DATA = 0 C NODAL = 0 SPATIAL_OVERLAP = 0 C MAX_MULTI = MULTIPLICITY C READ (IN_DATA1,REC=RECORD_NUM) BUFFER C X = BX Y = BY C CALL CCPI2I(BINTENBS,BUFFER(3),6,.TRUE.,.FALSE.) CALL CCPI2I(BSIGBS,BUFFER(6),6,.TRUE.,.FALSE.) CALL CCPI2I(BINTENPS,BUFFER(9),6,.TRUE.,.FALSE.) CALL CCPI2I(BSIGPS,BUFFER(12),6,.TRUE.,.FALSE.) DO 20 I = 1 , 6 INTEN_BOX_S(I) = BINTENBS(I) SIG_BOX_S(I) = BSIGBS(I) INTEN_PROF_S(I) = BINTENPS(I) SIG_PROF_S(I) = BSIGPS(I) 20 CONTINUE C CALL CCPBYI(MHKL,BUFFER(15),-4) CALL CCPBYI(HKL2,BUFFER(17),-4) H_S(1) = H K_S(1) = K L_S(1) = L H_S(2) = H2 K_S(2) = K2 L_S(2) = L2 C LAMBDAS(1) = LAMBDA1 LAMBDAS(2) = LAMBDA2 C DMIN_THRESHOLD = BDMIN C MULTIPLICITY = NBITST(M_FLAGS,0,5) IF (MULTIPLICITY.GT.MAX_MULTI) THEN WRITE (STR,6010) MULTIPLICITY - 1 , MAX_MULTI - 1 CALL TWRITE(STR) MULTIPLICITY = MAX_MULTI END IF C IF (NBITST(M_FLAGS,5,1).NE.0) THEN NODAL = 1 NODIDX = NIDX ELSE NODIDX = 0 END IF IF (NBITST(M_FLAGS,6,1).NE.0) SPATIAL_OVERLAP = 1 C IF (MULTIPLICITY.GT.2) THEN RINT = LAMBDAS(2) LINK_REC = IREAL DO 30 I = 2 , MULTIPLICITY C READ (IN_DATA2,REC=LINK_REC) BUFFER2 C LINK_REC = LINK_REC + 1 C CALL CCPBYI(MHKL,BUFFER2(1),-4) H_S(I) = H K_S(I) = K L_S(I) = L LAMBDAS(I) = LAMBDA C 30 CONTINUE END IF C RETURN C C====== Format statements C 6010 FORMAT ('Spot read that has',I3,' harmonics. Truncating ', + 'to just',I3,' harmonics.') END **==get_overlaps.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C******************************* C** GET_OVERLAPS ** C******************************* C C INTEGER FUNCTION GET_OVERLAPS ( IN_DATA1, IN_DATA2, RECORD_NUM, + END_REC, MAX_MULTI, + ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX, + X, Y, INTENSITIES, SIGMAS, NUM_GOOD_FILMS, + MULTIPLICITY, NODAL, SPATIAL_OVERLAP, + H_S, K_S, L_S, LAMBDAS, C2TH ) C C Purpose: Returns overlap data after applying rejection criteria C INCLUDE 'implicit_none.finc' C INTEGER IN_DATA1, IN_DATA2, RECORD_NUM, END_REC, MAX_MULTI C REAL ITOP, IBOT, N_SIGS REAL WAVMIN, WAVMAX C REAL X, Y, INTENSITIES(6), SIGMAS(6) INTEGER NUM_GOOD_FILMS INTEGER MULTIPLICITY, NODAL, SPATIAL_OVERLAP INTEGER H_S(MAX_MULTI), K_S(MAX_MULTI), L_S(MAX_MULTI) REAL LAMBDAS(MAX_MULTI), C2TH(MAX_MULTI) C C====== record_num : number of record to read from .ge1 file C====== max_multi : maximum multiplicity of overlap allowed C====== ITOP : maximum intensity allowed C====== IBOT : minimum intensity allowed C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C COMMON /BP/ BOXPROF INTEGER BOXPROF C C====== LOCALS: C REAL STHLSQ INTEGER GET_GE_DATA, APPLY_REJECT_CRIT INTEGER I, NODIDX REAL WAV, INTENB(6), SIGB(6), INTENP(6), SIGP(6) REAL H2, K2, L2, DSTAR, DHKL, THETA, SIN_THETA REAL DMIN_THRESHOLD C C------------------------------------------------------------------------------ C GET_OVERLAPS = 0 C C====== We want to know when mult>6 so need to use 7 here. C 20 MULTIPLICITY = 7 RECORD_NUM = RECORD_NUM + 1 IF (RECORD_NUM.GT.END_REC) THEN CALL TWRITE('Reached the end of file') GET_OVERLAPS = 1 RETURN END IF IF (GET_GE_DATA(IN_DATA1,IN_DATA2,RECORD_NUM,X,Y,INTENB,SIGB, + INTENP,SIGP,MULTIPLICITY,NODAL,NODIDX,SPATIAL_OVERLAP,H_S,K_S, + L_S,LAMBDAS,DMIN_THRESHOLD).NE.0) GO TO 20 IF (BOXPROF.EQ.1) THEN DO 30 I = 1 , 6 INTENSITIES(I) = INTENB(I) SIGMAS(I) = SIGB(I) 30 CONTINUE ELSE DO 40 I = 1 , 6 INTENSITIES(I) = INTENP(I) SIGMAS(I) = SIGP(I) 40 CONTINUE END IF C IF (MULTIPLICITY.EQ.1) GO TO 20 C C====== Now apply rejection criteria C IF (APPLY_REJECT_CRIT(MULTIPLICITY,INTENSITIES,SIGMAS,LAMBDAS, + SPATIAL_OVERLAP,ITOP,IBOT,N_SIGS,WAVMIN,WAVMAX,NUM_GOOD_FILMS) + .EQ.1) GO TO 20 C C====== If not enough good films then say so and reject C ccc IF (NUM_GOOD_FILMS.LT.MULTIPLICITY) GO TO 20 IF (NUM_GOOD_FILMS.LT.1) GO TO 20 IF (MULTIPLICITY.GT.6) GO TO 20 C C====== Calculate cos(2*theta) for the harmonics C DO 50 I = 1 , MULTIPLICITY C WAV = LAMBDAS(I) C C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2 C====== NB This function uses implicit input from SETRSL C DSTAR = 2*SQRT(STHLSQ(H_S(I),K_S(I),L_S(I))) DHKL = 1.0/DSTAR SIN_THETA = WAV*DSTAR/2.0 THETA = ASIN(SIN_THETA) C2TH(I) = COS(2.0*THETA) C 50 CONTINUE C RETURN END **==get_reject_crit.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C********************************** C** GET_REJECT_CRIT ** C********************************** C C SUBROUTINE GET_REJECT_CRIT (ACELL, BCELL, CCELL, CELLANG, + ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX) C C Purpose: To ge rejection criteria from user, as well as unit cell parameters C INCLUDE 'implicit_none.finc' C REAL ACELL, BCELL, CCELL, CELLANG(3) REAL ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX 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 REAL MILLER_INDEX, CP(6),WAV1,WAV2,R LOGICAL QUIT INTEGER I,J,NRANGE,ORDER CHARACTER*80 STR COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10) C C------------------------------------------------------------------------------ C C====== Get unit cell parameters and wavelength range from .GE1 file. C CALL READ_GE_PARAMS(CP,WAVMIN,WAVMAX) NRANGE = 2 CALL TASKI('Enter Number of Wavelength ranges',NRANGE,QUIT) DO 10 I = 1,NRANGE READ(IN_TERM,*) WAV1(I),WAV2(I),ORDER(I),(R(I,J),J=1,ORDER(I)+1) 10 CONTINUE WAVMIN=WAV1(1) WAVMAX=WAV2(NRANGE) ACELL = CP(1) BCELL = CP(2) CCELL = CP(3) DO 20 I = 1 , 3 CELLANG(I) = CP(I+3) 20 CONTINUE CALL TWRITE(' ') WRITE (STR,6010) CALL TWRITE(STR) WRITE (STR,6020) ACELL , BCELL , CCELL , CELLANG CALL TWRITE(STR) WRITE (STR,6030) WAVMIN , WAVMAX CALL TWRITE(STR) C C====== Get intensity cut-offs C CALL TASKF('Enter Itop for intensity range',ITOP,1,QUIT) CALL TASKF('Enter Ibot for intensity range',IBOT,1,QUIT) C C====== Get number of SD's before cut-off C CALL TASKF('Enter number of SD`s before cut-off',N_SIGS,1,QUIT) C C====== This concludes the interactive input of parameters C RETURN C C====== Format statements C 6010 FORMAT ('Cell Parameters: a b c alpha beta gamma') 6020 FORMAT ('=============== ',6F6.1) 6030 FORMAT ('Wavelength range ',2F5.2) END **==get_trans_factor.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*********************************** C** GET_TRANS_FACTOR ** C*********************************** C C SUBROUTINE GET_TRANS_FACTOR (PARAMS, DWAV, DC2TH, FILM, TF) C C Purpose: To calculate TRANSMISSION film factor for a given film and lambda C INCLUDE 'implicit_none.finc' C DOUBLE PRECISION PARAMS(2,3,5), DWAV, DC2TH INTEGER FILM DOUBLE PRECISION TF C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C C====== LOCALS: C INTEGER I, RANGE DOUBLE PRECISION DWAV3, K, A, ARG, FF CHARACTER*80 STR C C FF : film factor (absorption) C C------------------------------------------------------------------------------ C TF = 1.0 IF (FILM.EQ.1) RETURN C IF (DWAV.LT.0.49) THEN RANGE = 1 ELSE IF (DWAV.LT.0.92) THEN RANGE = 2 ELSE IF (DWAV.GE.0.92) THEN RANGE = 3 END IF C DWAV3 = DWAV*DWAV*DWAV C DO 20 I = 1 , FILM - 1 C K = PARAMS(1,RANGE,I) A = PARAMS(2,RANGE,I) C ARG = A*DWAV3/DC2TH IF (ARG.GT.80) THEN WRITE (STR,6010) CALL TWRITE(STR) WRITE (STR,6020) ARG CALL TWRITE(STR) CALL TWRITE(' ') ARG = 80.0 END IF FF = K*EXP(ARG) C TF = TF/FF C 20 CONTINUE C RETURN C C====== Format statements C 6010 FORMAT ('Possible overflow in calculating film factor, ', + 'K*exp(arg)') 6020 FORMAT (' with arg =',G11.4,' Truncating arg to 80.00') END **==lorentz.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************* C** LORENTZ ** C************************* C C REAL FUNCTION LORENTZ (THETA) C C Purpose: Apply Lorentz correction C INCLUDE 'implicit_none.finc' REAL THETA C C====== LOCALS: C REAL STH2 C C------------------------------------------------------------------------------ C STH2 = SIN(THETA)**2 C C====== lorentz = sth2 : changed by Pella 18.10.85 C LORENTZ = 1./STH2 RETURN END **==open_ge.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************** C** OPEN_GE ** C************************** C C SUBROUTINE OPEN_GE (IN_DATA1, IN_DATA2, END_POS, + END_REC) C C Purpose: To open .ge files using the given units and filename C Note that if the filename 'interactive' is supplied C then the user is prompted for a filename C Also note that filenames can have extensions, but these are removed C INCLUDE 'implicit_none.finc' C INTEGER IN_DATA1, IN_DATA2, END_REC INTEGER END_POS C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C LOGICAL CCPEXS C C====== LOCALS: C CHARACTER*80 STR INTEGER I,IFAIL LOGICAL EXISTS1, EXISTS2, QUIT C C GE1 : filename with .ge1 appended C GE2 : filename with .ge2 appended C C-------------------------------------------------------------------------- C C IFAIL = 0 CALL CCPDPN(-IN_DATA1,'GE1','READONLY','DU',19,IFAIL) C IFAIL = 0 CALL CCPDPN(-IN_DATA2,'GE2','READONLY','DU',2,IFAIL) C READ (UNIT=IN_DATA1,REC=10) I , END_REC IF (END_REC.LT.1) THEN READ (UNIT=IN_DATA1,REC=9) I , I , END_REC END IF END_REC = END_REC + 10 C RETURN C END **==output_deconv.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C********************************* C** OUTPUT_DECONV ** C********************************* C C SUBROUTINE OUTPUT_DECONV (NUM_OVERLAPS, OVERLAP_INTENS, + OVERLAP_LAMBDAS, OVERLAP_HKLS, OVERLAP_DHKLS, + OVERLAP_SIGMAS, + OVERLAP_XY, + PARAMS, ACELL, BCELL, CCELL, CELLANG) C C Purpose: To output deconvoluted reflections suitable as input for LAUESCAL C INCLUDE 'implicit_none.finc' C INTEGER NUM_OVERLAPS REAL OVERLAP_INTENS(6,NUM_OVERLAPS) REAL OVERLAP_SIGMAS(6,NUM_OVERLAPS) REAL OVERLAP_LAMBDAS(NUM_OVERLAPS) REAL OVERLAP_DHKLS(NUM_OVERLAPS) REAL OVERLAP_XY(2,NUM_OVERLAPS) INTEGER OVERLAP_HKLS(NUM_OVERLAPS) C DOUBLE PRECISION PARAMS(2,3,5) REAL ACELL, BCELL, CCELL, CELLANG(3) 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 C INTEGER NWORDS, NVALS REAL VALUES(10) CHARACTER*20 WORDS(10) INTEGER CLASS, MULTI INTEGER SCALE_SPOT LOGICAL QUIT REAL SCALE_FACTOR INTEGER I, J, NS, IFAIL REAL X, Y, INTENSITIES(6), SIGMAS(6) INTEGER H, K, L INTEGER I_INTENS, I_SIGMA REAL LAMBDA, THETA, LOR, P, SCALED_INTENS, SCALED_SIGMA C INTEGER MAX_BAD_SPOTS, BAD_SPOT_COUNT PARAMETER (MAX_BAD_SPOTS=10) C C C bad_spot_count: counts the number of times an integer bad_spot might occur C and if this exceeds max_bad_spots then the routine aborts C INTEGER HKL(4) C C------------------------------------------------------------------------------ C BAD_SPOT_COUNT = 0 C C====== Open output file for LAUESCALE/LAUENORM C IFAIL = 1 CALL CCPDPN(-4,'OUT_FILE','NEW','F',0,IFAIL) IF (IFAIL.EQ.-1) GO TO 220 CALL TWRITE(' ') CALL TWRITE(' ') C C====== Get class of spot to output C 30 CALL TPROMP('Enter class of spot to output (all,2,3 etc.): ') CALL TREAD(STR,QUIT) CALL strtkn(STR,' ,',10,WORDS,NWORDS,VALUES,NVALS) C IF (NVALS.GT.1) GO TO 30 IF (NWORDS.GT.1) GO TO 30 IF (NWORDS.EQ.NVALS) GO TO 30 C CLASS = -1 IF (NWORDS.EQ.1) THEN IF (WORDS(1).EQ.'all' .OR. WORDS(1).EQ.'ALL') CLASS = 0 IF (WORDS(1).EQ.'doublets' .OR. WORDS(1).EQ.'DOUBLETS') + CLASS = 2 IF (WORDS(1).EQ.'triplets' .OR. WORDS(1).EQ.'TRIPLETS') + CLASS = 3 ELSE CLASS = NINT(VALUES(1)) IF (CLASS.LT.2 .OR. CLASS.GT.6) GO TO 30 END IF IF (CLASS.EQ.-1) GO TO 30 WRITE (STR,'(''Class ='',I2)') CLASS CALL TWRITE(STR) C C====== Get overall scaling factor C SCALE_FACTOR = 100. CALL TASKF('Input scaling factor',SCALE_FACTOR,2,QUIT) C NS = 0 C C====== Use CCP4 routine SETRSL to initialise storage for STHLSQ C====== (STHLSQ is called inside scale_spot) C CALL SETRSL(ACELL,BCELL,CCELL,CELLANG(1),CELLANG(2),CELLANG(3)) C DO 50 I = 1 , NUM_OVERLAPS C C====== Is this spot in the right output class ? C CALL CCPBYI(HKL,OVERLAP_HKLS(I),-4) MULTI = HKL(4) MULTI = IABS(MULTI) IF (CLASS.EQ.0 .OR. CLASS.EQ.MULTI) THEN C H = HKL(1) K = HKL(2) L = HKL(3) C LAMBDA = OVERLAP_LAMBDAS(I) X = OVERLAP_XY(1,I) Y = OVERLAP_XY(2,I) C DO 40 J = 1 , 6 INTENSITIES(J) = OVERLAP_INTENS(J,I) SIGMAS(J) = OVERLAP_SIGMAS(J,I) 40 CONTINUE C IFAIL = SCALE_SPOT(PARAMS,SCALE_FACTOR,X,Y,INTENSITIES, + SIGMAS,H,K,L,LAMBDA,THETA,SCALED_INTENS, + SCALED_SIGMA,I_INTENS,I_SIGMA,LOR,P) IF (IFAIL.NE.0) THEN C C====== Bad spot - reject and maybe abort C IF (IFAIL.NE.2) BAD_SPOT_COUNT = BAD_SPOT_COUNT + 1 IF (BAD_SPOT_COUNT.GT.MAX_BAD_SPOTS) THEN CALL TWRITE(' ') WRITE (STR,6040) MAX_BAD_SPOTS CALL TWRITE(STR) WRITE (STR,6050) CALL TWRITE(STR) CALL TWRITE(' ') CALL TWRITE(' ') GO TO 60 END IF ELSE C C C====== Write out to output file C NS = NS + 1 WRITE (4,6010) H , K , L , LAMBDA , THETA , + I_INTENS , I_SIGMA , X , Y END IF END IF C 50 CONTINUE C 60 CALL TWRITE(' ') WRITE (STR,6030) NS CALL TWRITE(STR) C CLOSE (4) C RETURN C 220 STOP C====== Format statements C 6010 FORMAT (3I5,2F10.5,2I10,2F8.3) 6020 FORMAT ('Writing output to ',A) 6030 FORMAT ('Total number of output reflections',I6) 6040 FORMAT ('More than',I4,' bad spots detected.') 6050 FORMAT ('***** ABORTING *****') END **==polar.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************* C** POLAR ** C************************* C C REAL FUNCTION POLAR (TAU,AX,AY,THETA) C C Purpose: Apply polarization correction C INCLUDE 'implicit_none.finc' C REAL TAU, AX, AY, THETA C C====== LOCALS: C REAL PSI, C2TH2, ALPHA, P0 C C------------------------------------------------------------------------------ C IF (AY.EQ.0.) THEN C POLAR = 1.0 C ELSE C PSI = ATAN(AX/AY) C2TH2 = COS(2.0*THETA)**2 P0 = (1.0+C2TH2)/2.0 ALPHA = COS(2.0*PSI)*SIN(2*THETA)**2 ALPHA = ALPHA/(1.0+C2TH2) POLAR = P0*(1.0-TAU*ALPHA) C END IF C RETURN END **==read_ge_params.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C********************************* C** READ_GE_PARAMS ** C********************************* C C SUBROUTINE READ_GE_PARAMS (CP, WAVMIN, WAVMAX) C INCLUDE 'implicit_none.finc' C C Purpose: C Read lambda limits and cell parameters from .GE1 file header. Expand C cell parameters up to complete set (depends on xtal system) and convert C to direct cell. C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C REAL WAVE, WAVMIN, WAVMAX REAL CP(6), BUFF (9), W_FS, BUFF2(3) REAL PARAMS(6) INTEGER XALNO, I INTEGER ISYLAT(2), ISL CHARACTER*40 BUF C READ (IN_DATA,REC=3) W_FS , BUFF , CP , ISL CALL CCPI2I(ISYLAT,ISL,2,.TRUE.,.FALSE.) XALNO = ISYLAT(1) READ (IN_DATA,REC=4) BUF , BUFF2 , WAVE , WAVMIN , WAVMAX C IF (XALNO.EQ.1) THEN C C====== Triclinic C CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.2) THEN C C====== Monoclinic (unique axis b) C CP(5) = CP(4) CP(4) = 90. CP(6) = 90. CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.3) THEN C C====== Orthorhombic C CP(4) = 90. CP(5) = 90. CP(6) = 90. CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.4) THEN C C====== Tetragonal C CP(3) = CP(2) CP(2) = CP(1) CP(4) = 90. CP(5) = 90. CP(6) = 90. CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.5) THEN C C====== Hexagonal axes C CP(3) = CP(2) CP(2) = CP(1) CP(4) = 90. CP(5) = 90. CP(6) = 60. CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.6) THEN C C====== Rhombohedral axes C CP(4) = CP(2) CP(5) = CP(2) CP(6) = CP(2) CP(2) = CP(1) CP(3) = CP(1) CALL RECCEL(PARAMS,CP,WAVE) ELSE IF (XALNO.EQ.7) THEN C C====== Cubic C CP(2) = CP(1) CP(3) = CP(1) CP(4) = 90. CP(5) = 90. CP(6) = 90. CALL RECCEL(PARAMS,CP,WAVE) ELSE STOP ' ERROR IN HEADER DATA ' END IF DO 20 I = 1 , 6 CP(I) = PARAMS(I) 20 CONTINUE RETURN END **==read_vc.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************** C** READ_VC ** C************************** C C SUBROUTINE READ_VC(PARAMS) C C Purpose: To read in Victoreen parameters from unit 10 C INCLUDE 'implicit_none.finc' C DOUBLE PRECISION PARAMS(2,3,5),AL150(2,3,5),AL200(2,3,5) C COMMON /IO_STREAMS/IN_TERM,OUT_TERM,IN_DATA,OUT_DATA INTEGER IN_TERM,OUT_TERM,IN_DATA,OUT_DATA C C====== FUNCTION TYPES: C LOGICAL CCPEXS INTEGER LENSTR C C====== LOCALS: C CHARACTER*80 STR LOGICAL EXISTS, QUIT INTEGER COMM INTEGER I,J,K,TYPE,IFAIL C DATA AL150/1,1.975,1,.682,1,.2856,1,2.4912,1,1.1982,1,.8018,1, + 2.4912,1,1.1982,1,.8018,1,2.4912,1,1.1982,1,.8018,1,7.1454,1, + 5.276,1,4.1812/ DATA AL200/1,1.975,1,.682,1,.2856,1,2.6633,1,1.3703,1,.9739,1, + 2.6633,1,1.3703,1,.9739,1,2.6633,1,1.3703,1,.9739,1,7.1454,1, + 5.276,1,4.1812/ C C------------------------------------------------------------------------------ C TYPE = 1 CALL TWRITE(' ') CALL TWRITE(' ') CALL TWRITE('Options available on loading film factors:-') CALL TWRITE(' ') CALL TWRITE('0 - single image plate') CALL TWRITE('1 - AL150') CALL TWRITE('2 - AL200') CALL TWRITE('3 - enter any other .VC file') CALL TWRITE(' ') COMM = 1 20 CALL TASKI('Enter option',COMM,QUIT) IF (COMM.LT.0 .OR. COMM.GT.3) GO TO 20 C C====== Now determine which option is to be used C IF (COMM.EQ.0) THEN DO 25 K = 1 , 5 DO 24 J = 1 , 3 PARAMS(1,J,K) = 1. PARAMS(2,J,K) = 0. 24 CONTINUE 25 CONTINUE RETURN END IF IF (COMM.EQ.1) THEN DO 50 K = 1 , 5 DO 40 J = 1 , 3 DO 30 I = 1 , 2 PARAMS(I,J,K) = AL150(I,J,K) 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END IF IF (COMM.EQ.2) THEN DO 80 K = 1 , 5 DO 70 J = 1 , 3 DO 60 I = 1 , 2 PARAMS(I,J,K) = AL200(I,J,K) 60 CONTINUE 70 CONTINUE 80 CONTINUE RETURN END IF IF (COMM.EQ.3) THEN IFAIL = 1 CALL CCPDPN(-10,'VC','READONLY','F',0,IFAIL) IF (IFAIL.EQ.-1) GO TO 220 C DO 150 K = 1 , 5 DO 140 J = 1 , 3 DO 130 I = 1 , 2 READ (10,*) PARAMS(I,J,K) 130 CONTINUE 140 CONTINUE 150 CONTINUE C CLOSE (10) END IF C RETURN C 220 STOP ' .vc file not found' END **==scale_spot.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C***************************** C** SCALE_SPOT ** C***************************** C C INTEGER FUNCTION SCALE_SPOT (PARAMS, SCALE_FACTOR, + X, Y, INTENSITIES, SIGMAS, H, K, L, + LAMBDA, THETA, SCALED_INTENS, + SCALED_SIGMA, I_INTENS, + I_SIGMA, LOR,P) C C Purpose: To scale one spot combining the intensities from 6 films C INCLUDE 'implicit_none.finc' C DOUBLE PRECISION PARAMS(2,3,5) REAL SCALE_FACTOR C REAL X, Y, INTENSITIES(6), SIGMAS(6) INTEGER H, K, L REAL LAMBDA C REAL THETA, SCALED_INTENS, SCALED_SIGMA, LOR, P INTEGER I_INTENS, I_SIGMA C COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C C====== LOCALS: C CHARACTER*80 STR REAL STHLSQ REAL LORENTZ, POLAR REAL KP, T PARAMETER (KP=0.105, T=1.0) C REAL MEAN_INTENS, MEAN_SIG, SIG1, SIG2 C REAL SUM_OF_WEIGHTS, WEIGHTED_SUM, WEIGHTS(6) REAL SUM_OF_WEIGHTS2, WEIGHTED_RESID2 C REAL FACTOR C INTEGER I, J C REAL DSTAR, DHKL, SIN_THETA C C------------------------------------------------------------------------------ C SCALE_SPOT = 0 C C====== Check for overload C DO 20 J = 1 , 6 IF (ABS(INTENSITIES(J)).EQ.9999.0 .OR. ABS(SIGMAS(J)) + .EQ.9999.0) THEN INTENSITIES(J) = 0 SIGMAS(J) = 1E6 END IF 20 CONTINUE C C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2 C====== NB This function uses implicit input from SETRSL C DSTAR = 2*SQRT(STHLSQ(H,K,L)) DHKL = 1.0/DSTAR SIN_THETA = LAMBDA*DSTAR/2.0 C C====== Check that sin theta makes sense C IF (ABS(SIN_THETA).GT.1.0) THEN CALL TWRITE(' ') WRITE (STR,6030) SIN_THETA CALL TWRITE(STR) WRITE (STR,6040) H , K , L CALL TWRITE(STR) CALL TWRITE(' ') SIN_THETA = SIGN(1.0,SIN_THETA) END IF THETA = ASIN(SIN_THETA) C CALL SCALE_UP(PARAMS,INTENSITIES,SIGMAS,LAMBDA,THETA) C SUM_OF_WEIGHTS = 0.0 WEIGHTED_SUM = 0.0 C DO 30 J = 1 , 6 IF (SIGMAS(J).GE.1E6) THEN WEIGHTS(J) = 0.0 ELSE WEIGHTS(J) = 1.0/(SIGMAS(J)*SIGMAS(J)) END IF SUM_OF_WEIGHTS = SUM_OF_WEIGHTS + WEIGHTS(J) WEIGHTED_SUM = WEIGHTED_SUM + WEIGHTS(J)*INTENSITIES(J) 30 CONTINUE C IF (SUM_OF_WEIGHTS.LT.1E-6) THEN SCALE_SPOT = 2 GO TO 50 END IF C MEAN_INTENS = WEIGHTED_SUM/SUM_OF_WEIGHTS SCALED_INTENS = MEAN_INTENS SIG1 = SQRT(1.0/SUM_OF_WEIGHTS) C WEIGHTED_RESID2 = 0.0 SUM_OF_WEIGHTS2 = 0.0 C DO 40 J = 1 , 6 WEIGHTED_RESID2 = WEIGHTED_RESID2 + + ((INTENSITIES(J)-MEAN_INTENS)*WEIGHTS(J))**2 SUM_OF_WEIGHTS2 = SUM_OF_WEIGHTS2 + WEIGHTS(J)**2 40 CONTINUE C SIG2 = SQRT(WEIGHTED_RESID2/SUM_OF_WEIGHTS2) IF (SIG2.GT.SIG1) THEN MEAN_SIG = SIG2 ELSE MEAN_SIG = SIG1 END IF SCALED_SIGMA = MEAN_SIG C LOR = LORENTZ(THETA) P = POLAR(1.0,X,Y,THETA) FACTOR = SCALE_FACTOR*EXP(KP*T/COS(2.0*THETA))/(LOR*P) C MEAN_INTENS = MEAN_INTENS*FACTOR MEAN_SIG = MEAN_SIG*FACTOR C IF (MEAN_INTENS.GE.1E9 .OR. MEAN_SIG.GE.1E9) THEN CALL TWRITE(' ') WRITE (STR,6010) MEAN_INTENS CALL TWRITE(STR) WRITE (STR,6020) CALL TWRITE(STR) SCALE_SPOT = 1 ELSE I_INTENS = MEAN_INTENS I_SIGMA = MEAN_SIG END IF C 50 RETURN C C====== Format statements C 6010 FORMAT ('Intensity =',G11.4,', danger of integer overflow. ', + 'Change scaling factor.') 6020 FORMAT ('** Ignoring this spot **'/) 6030 FORMAT ('Warning - lambda/2d (=sin theta) =',F5.1, + ' -setting theta to 90') 6040 FORMAT ('(h,k,l) = (',I4,',',I4,',',I4,')') END **==scale_up.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C*************************** C** SCALE_UP ** C*************************** C C SUBROUTINE SCALE_UP (PARAMS, INTENSITIES, SIGMAS, LAMBDA, THETA) C C Purpose: Scale up film intensities and sigmas to A film C INCLUDE 'implicit_none.finc' C DOUBLE PRECISION PARAMS(2,3,5) REAL INTENSITIES(6), SIGMAS(6) REAL LAMBDA, THETA C C====== LOCALS C DOUBLE PRECISION DWAV, DC2TH, TF REAL FF INTEGER FILM C C------------------------------------------------------------------------------ C DWAV = LAMBDA DC2TH = COS(2.0*THETA) C C====== Calculate total scales C DO 20 FILM = 2 , 6 C C====== Deconv has reset bad sigmas to 1e6 C IF (SIGMAS(FILM).LT.1E6) THEN CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF) IF (TF.LT.1E-5) THEN FF = 0.0 ELSE FF = 1.0/TF END IF INTENSITIES(FILM) = INTENSITIES(FILM)*FF SIGMAS(FILM) = SIGMAS(FILM)*FF C C====== Normal bad (-9999) sigmas caught here C IF (SIGMAS(FILM).LE.0.0) SIGMAS(FILM) = 1E6 END IF 20 CONTINUE C RETURN C END **==sthlsq.f C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C************************* C** STHLSQ ** C************************* C C FUNCTION STHLSQ(IH, IK, IL ) C ============================ C INTEGER IH,IK,IL COMMON /RECPLT/ COEFHH, COEFHK, COEFHL, COEFKK, COEFKL, COEFLL C CALCULATE (SIN(THETA)/LAMBDA)**2 FROM H,K,L; COEF'S SET BY CALL TO C SETRSL : GOOD FOR ANY KIND OF AXES C C------------------------------------------------------------------------------- C STHLSQ = IH*IH*COEFHH + IH*IK*COEFHK + IH*IL*COEFHL + + IK*IK*COEFKK + IK*IL*COEFKL + IL*IL*COEFLL RETURN END **==deconv_spots.f C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C******************************** C** DECONV_SPOTS ** C******************************** C C SUBROUTINE DECONV_SPOTS (MAX_OVERLAPS, END_REC, + ACELL, BCELL, CCELL, CELLANG, + ITOP, IBOT, N_SIGS, WAVMIN, WAVMAX, + PARAMS, + NUM_OVERLAPS, OVERLAP_INTENS, + OVERLAP_SIGMAS, OVERLAP_LAMBDAS, + OVERLAP_HKLS, OVERLAP_DHKLS, + OVERLAP_XY) C C Purpose: To deconvolute intensities for harmonics C INCLUDE 'implicit_none.finc' C INTEGER MAX_OVERLAPS, END_REC C REAL ACELL, BCELL, CCELL, CELLANG(3) REAL ITOP, IBOT REAL N_SIGS REAL WAVMIN, WAVMAX C DOUBLE PRECISION PARAMS(2,3,5) C INTEGER NUM_OVERLAPS REAL OVERLAP_INTENS(6,MAX_OVERLAPS) REAL OVERLAP_SIGMAS(6,MAX_OVERLAPS) REAL OVERLAP_LAMBDAS(MAX_OVERLAPS) REAL OVERLAP_DHKLS(MAX_OVERLAPS) REAL OVERLAP_XY(2,MAX_OVERLAPS) INTEGER OVERLAP_HKLS(MAX_OVERLAPS) C C COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10) INTEGER NRANGE,ORDER REAL WAV1,WAV2,R COMMON /IO_STREAMS/ IN_TERM, OUT_TERM, IN_DATA, OUT_DATA INTEGER IN_TERM, OUT_TERM, IN_DATA, OUT_DATA C C====== LOCALS: C REAL STHLSQ INTEGER GET_OVERLAPS LOGICAL QUIT,MTZ C INTEGER NUM_SPOT_LOG,NUM_LOGGED C INTEGER OLD_OUT_TERM, RECORD_NUM, SPOT_NUM C INTEGER MAX_MULTI PARAMETER (MAX_MULTI=6) C REAL IIN(6), SIGIN(6), DELTAS(6), W_DELTAS(6) REAL TF_MAT(6,MAX_MULTI) REAL IOUT(MAX_MULTI), SIGOUT(MAX_MULTI), IFILM C DOUBLE PRECISION TF, DWAV, DC2TH INTEGER FILM, ROW INTEGER I, J, M, NUM_NEG LOGICAL IS_NEG C INTEGER H, K, L, MULTI, HKL(4) REAL DSTAR, DHKL,RESOL C REAL X, Y, INTENSITIES(6), SIGMAS(6),rat(6),ratall INTEGER NUM_GOOD_FILMS INTEGER MULTIPLICITY, NODAL, SPATIAL_OVERLAP INTEGER H_S(MAX_MULTI+1), K_S(MAX_MULTI+1), L_S(MAX_MULTI+1) REAL LAMBDAS(MAX_MULTI+1), C2TH(MAX_MULTI+1) REAL sigtemp,lamcur,rattt INTEGER IFAIL C INTEGER IPOINT(4) CHARACTER*30 LABEL(4) CHARACTER*1 CTYP(4) LOGICAL EOF REAL ADATA(4) C EQUIVALENCE (HKL(1),H) EQUIVALENCE (HKL(2),K) EQUIVALENCE (HKL(3),L) EQUIVALENCE (HKL(4),MULTI) C DATA LABEL/'H','K','L','FC'/ DATA CTYP/'H','H','H','F'/ C C------------------------------------------------------------------------------ C NUM_LOGGED = 0 NUM_SPOT_LOG = 0 CALL TASKI('How many spots do you want to output on log file', + NUM_SPOT_LOG,QUIT) DO 75 I =1,4 IPOINT(I)=-1 75 CONTINUE CALL LROPEN (3,'HKLIN',1,IFAIL) IF (IFAIL.EQ.1) STOP CALL LRASSN (3,LABEL,4,IPOINT,CTYP) C C====== Use CCP4 routine SETRSL to initialise storage for STHLSQ C CALL SETRSL(ACELL,BCELL,CCELL,CELLANG(1),CELLANG(2),CELLANG(3)) C C====== Loop through entire file, deconvoluting on the way C NUM_OVERLAPS = 0 RECORD_NUM = 10 NUM_NEG = 0 C 20 IF (GET_OVERLAPS(IN_DATA,3,RECORD_NUM,END_REC,MAX_MULTI,ITOP + ,IBOT,N_SIGS,WAVMIN,WAVMAX,X,Y,INTENSITIES,SIGMAS, + NUM_GOOD_FILMS,MULTIPLICITY,NODAL,SPATIAL_OVERLAP,H_S,K_S,L_S, + LAMBDAS,C2TH).EQ.0) THEN C C====== Write out spot info C SPOT_NUM = RECORD_NUM - 10 IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,6010) SPOT_NUM , X , Y , MULTIPLICITY , NODAL WRITE (8,6020) (INTENSITIES(I),I=1,6) , (SIGMAS(I),I=1,6) WRITE (8,6030) END IF DO 30 I = 1 , MULTIPLICITY IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,6040) H_S(I) , K_S(I) , L_S(I) , LAMBDAS(I) END IF CALL LRREFF(3,RESOL,ADATA,EOF) IF (EOF) GO TO 35 rat(i)=ADATA(4)/lamcur(lambdas(i)) if(i.eq.1)rattt=rat(i) rat(i)=rat(i)/rattt 30 CONTINUE C C====== Put the valid intensities into input vector C 35 ROW = 0 DO 50 FILM = 1 , 6 IF (INTENSITIES(FILM).GE.0.0) THEN ROW = ROW + 1 ratall=0. DO 40 M = 1 , MULTIPLICITY DWAV = LAMBDAS(M) DC2TH = C2TH(M) CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF) TF_MAT(ROW,M) = TF ratall=ratall+rat(m)*tf 40 CONTINUE IIN(ROW) = INTENSITIES(FILM)/ratall SIGIN(ROW) = 1./SIGMAS(FILM)**2 if(row.ne.1)goto 50 do 45 k=1,multiplicity iout(k)=iin(row)*rat(k) sigout(k)=1./sqrt(sigin(row)) 45 continue END IF 50 CONTINUE C C====== Display transmission factor matrix C CALL MTX_PACK(TF_MAT,TF_MAT,6,MAX_MULTI,ROW,MULTIPLICITY) IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,6090) CALL MTX_WRITE(8,TF_MAT,ROW,MULTIPLICITY) END IF C C====== Resulting deconvoluted spot in iout array C====== set reject flag if any negative components C IS_NEG = .FALSE. DO 60 I = 1 , MULTIPLICITY IF (IOUT(I).LT.0) IS_NEG = .TRUE. 60 CONTINUE C C====== Display results C IF (NUM_LOGGED.LT.NUM_SPOT_LOG) WRITE (8,6050) C DO 70 FILM = 1 , 6 IF (INTENSITIES(FILM).GE.0.0) THEN DELTAS(FILM) = -1.0*INTENSITIES(FILM) ELSE DELTAS(FILM) = -9999.0 END IF 70 CONTINUE C DO 90 I = 1 , MULTIPLICITY C C====== Record intensities, lambdas, hkls and dhkls for R-symm calculation C NUM_OVERLAPS = NUM_OVERLAPS + 1 IF (NUM_OVERLAPS.GT.MAX_OVERLAPS) THEN IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,*) ' Reached overlap limit of' , MAX_OVERLAPS END IF GO TO 130 END IF C OVERLAP_LAMBDAS(NUM_OVERLAPS) = LAMBDAS(I) C H = H_S(I) K = K_S(I) L = L_S(I) MULTI = MULTIPLICITY CALL CCPIBY(OVERLAP_HKLS(NUM_OVERLAPS),HKL,-4) C C====== Use CCP4 function STHLSQ to calculate sin(theta)**2/lambda**2 C====== NB This function uses implicit input from SETRSL C DSTAR = 2*SQRT(STHLSQ(H_S(I),K_S(I),L_S(I))) C DHKL = 1.0/DSTAR IF (DHKL.LE.0.0) THEN WRITE (8,*) ' Dhkl = ' , DHKL , ' OV=' , NUM_OVERLAPS END IF OVERLAP_DHKLS(NUM_OVERLAPS) = DHKL OVERLAP_XY(1,NUM_OVERLAPS) = X OVERLAP_XY(2,NUM_OVERLAPS) = Y C DWAV = LAMBDAS(I) DC2TH = C2TH(I) IFILM = IOUT(I) DO 80 FILM = 1 , 6 C C====== Reject spots which got bad sigmas by setting them to 0 C====== -following programs will trap this. C IF (SIGOUT(I).EQ.-9999.0) THEN OVERLAP_INTENS(FILM,NUM_OVERLAPS) = 0.0 OVERLAP_SIGMAS(FILM,NUM_OVERLAPS) = 0.0 ELSE CALL GET_TRANS_FACTOR(PARAMS,DWAV,DC2TH,FILM,TF) OVERLAP_INTENS(FILM,NUM_OVERLAPS) = TF*IFILM OVERLAP_SIGMAS(FILM,NUM_OVERLAPS) = TF*SIGOUT(I) IF (DELTAS(FILM).NE.-9999.0) THEN DELTAS(FILM) = DELTAS(FILM) + + OVERLAP_INTENS(FILM,NUM_OVERLAPS) END IF END IF 80 CONTINUE IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,6060) LAMBDAS(I) , + (OVERLAP_INTENS(J,NUM_OVERLAPS),J=1,6) WRITE (8,6070) (OVERLAP_SIGMAS(J,NUM_OVERLAPS),J=1,6) END IF 90 CONTINUE IF (NUM_LOGGED.LT.NUM_SPOT_LOG) WRITE (8,6100) + (DELTAS(I),I=1,6) DO 100 I = 1 , 6 IF (DELTAS(I).NE.-9999.0) THEN W_DELTAS(I) = DELTAS(I)/(DELTAS(I)+2*INTENSITIES(I)) ELSE W_DELTAS(I) = -9999.0 END IF 100 CONTINUE IF (NUM_LOGGED.LT.NUM_SPOT_LOG) THEN WRITE (8,6110) (W_DELTAS(I),I=1,6) WRITE (8,6080) END IF C C====== That was the last log file output for this spot so inc num_logged. C NUM_LOGGED = NUM_LOGGED + 1 C C====== Loop for more spots C C====== Reset counter if need to reject this harmonic C====== because of negative components C IF (IS_NEG) THEN NUM_OVERLAPS = NUM_OVERLAPS - MULTIPLICITY NUM_NEG = NUM_NEG + 1 ELSE C C====== Only save top film measurement C DO 120 J = 1 , MULTIPLICITY DO 110 I = 2 , 6 OVERLAP_SIGMAS(I,NUM_OVERLAPS-MULTIPLICITY+J) = -9999. 110 CONTINUE 120 CONTINUE C END IF GO TO 20 END IF C C====== Now reached end of .ge file or too many overlaps C 130 IF (NUM_LOGGED.LT.NUM_SPOT_LOG) WRITE (8,6120) NUM_OVERLAPS C IF (NUM_LOGGED.LT.NUM_SPOT_LOG) WRITE (8,6130) NUM_NEG RETURN C C====== Format statements C 6010 FORMAT (1X,130('-'),/,1X,'Spot ',I5,' X = ',G11.4,' Y = ', + G11.4,' Multiplicity = ',I2,' Nodal flag = ',I1) 6020 FORMAT (' Observed intensities on A..F films, and sigmas'/,1X,6X, + 1X,6G14.6,/,8X,6G14.6) 6030 FORMAT (1X,' h k l lambda') 6040 FORMAT (1X,3I5,5X,G11.4) 6050 FORMAT (1X,'Wavelength, Deconvoluted intensities on A..F films,', + ' and sigma') 6060 FORMAT (1X,F6.3,',',6G14.6) 6070 FORMAT (1X,'sigma= ',6G14.6) 6080 FORMAT (1X,130('.')) 6090 FORMAT (9X,'Transmission Factor Matrix') 6100 FORMAT (1X,'Deltas=',6G14.6) 6110 FORMAT (1X,'W. Del=',6G14.6) 6120 FORMAT (/' Number of reflections read and deconvoluted =',I6/) 6130 FORMAT (/' Number of rejected reflections (negative component',1X, + ')=',I6/) END **==lamcur.f function lamcur(lam) COMMON /LAUENORM/NRANGE,WAV1(8),WAV2(8),ORDER(8),R(8,10) integer kk,NRANGE,ORDER real WAV1,WAV2,lamcur,lam,lamda,r lamda=lam DO 50 I=1,8 if(lamda.LE.WAV2(I)) THEN KK=I GOTO 60 END IF 50 CONTINUE 60 lamcur=R(KK,ORDER(KK)+1) do 100 i=ORDER(KK),1,-1 lamcur=lamcur*lamda+r(kk,i) 100 continue end