linse/sics
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
f3853c20f0868e2396766978712d5bab75206705
sics/difrac/profil.f
cvs ff5e8cf0b2 - Improved centering in DIFRAC 
- Fixed a bug in UserWait
- Improved scan message in scancom
- Added zero point correction in lin2ang
- fixed an issue with uuencoded messages
2000-04-06 12:18:53 +00:00

476 lines
18 KiB
Fortran
Raw Blame History

C-----------------------------------------------------------------------
C Subroutine to find peak limits in the profile
C-----------------------------------------------------------------------
SUBROUTINE PROFIL
INCLUDE 'COMDIF'
DIMENSION SMOOTH(500),IHTAGS(4),ID(500)
C-----------------------------------------------------------------------
C
C Explanation of symbols used for Profile Analysis
C
C PLIM Probability lower limit. Arbitrarily 0.01
C A12 Ratio of Int(alpha1)/Int(alpha2) (1.8)
C NWIND No. of profile points in the test window. (6)
C IDEL No. of pts in the profile + 1
C COUNT Sum of all profile points
C FRAC Ratio of 0.5*Background time/Peak time
C SIGLIM Inet significance limit
C CON No. of profile pts per degree scan (STEPDG). This
C gives a power of two steps for all speeds.
C SPEED Scan speed in degs per min.
C D12 Alpha1 to Alpha2 seperation in degrees
C ITYPE Scan type indicator. 0 or 1 2Theta; 2 or 3 Omega
C AS Scan before Alpha1 in degrees
C CS Scan after Alpha2 in degrees
C ACOUNT(I) Array of profile intensity values
C IWARN Warning flag from the measuring routine. 0 = OK.
C IPRFLG Profile analysis indicator. 0 = Do; 1 = Dont.
C BGRD1 Low angle background, taken for FRAC*Peak-time
C BGRD2 High angle background, as BGRD1.
C RSW PDP8E Read switch register routine.
C RSW(N,J) Reads 1-bit switch N into J
C STEPOF Fraction of As (and Cs) to step off from Alpha1
C (and Alpha2) before starting the profile analysis
C
C-----------------------------------------------------------------------
DATA PLIM/0.01/,A12/1.8/
IF ((IPRFLG .NE. 0 .AND. KI(1:1) .NE. 'G') .OR.
$ IDEL .LT. 10) RETURN
C-----------------------------------------------------------------------
C Results are sent to the printer for either :--
C 1. Individual measurements, i.e. not part of GO; or
C 2. Part of GO and Switch 4 is set to 1.
C-----------------------------------------------------------------------
IF (KI(1:1) .EQ. 'G') THEN
CALL RSW (4,ILPT)
ELSE
ILPT = 0
ENDIF
A1 = A12/(A12 + 1.0)
A2 = 1.0 - A1
NWIND = 6
ILOW = 1
NP = IDEL - 1
IHIGH = NP
SUM = COUNT
FRAC1 = FRAC
SIGLIM = 2.0
CON = STEPDG
RD12 = CON*D12
IF (ITYPE .EQ. 2 .OR. ITYPE .EQ. 3) RD12 = RD12*0.5
NXPTS = 1000.0/((AS + CS)*CON + RD12)
ID12 = RD12 + 0.5
ITOL = CON*(AS + CS)/8.0
MAXI = 0
C-----------------------------------------------------------------------
C Do not try to process peak top measurements
C-----------------------------------------------------------------------
IF (ITYPE .GE. 4) RETURN
CMAX = 0
DO 100 I = 1,NP
ACT = ACOUNT(I+1)
IF (CMAX .LT. ACT) CMAX = ACT
100 ACOUNT(I) = ACT
C-----------------------------------------------------------------------
C Smooth the profile (5-point average)
C-----------------------------------------------------------------------
SMOOTH(1) = (ACOUNT(1) + ACOUNT(2) + ACOUNT(3))/3.0
SMOOTH(2) = (3.0*SMOOTH(1) + ACOUNT(4))/4.0
DO 110 I = 3,NP-2
SMOOTH(I) = (ACOUNT(I-2) + ACOUNT(I-1) + ACOUNT(I) +
$ ACOUNT(I+1) + ACOUNT(I+2))/5.0
110 CONTINUE
SMOOTH(NP-1) = (SMOOTH(NP-2)*5 - ACOUNT(NP-4))/4.0
SMOOTH(NP) = (4.0*SMOOTH(NP-1) - ACOUNT(NP-3))/3.0
C-----------------------------------------------------------------------
C Test if peak is OK from MESINT or profile not needed
C-----------------------------------------------------------------------
IF (IWARN .NE. 0) GO TO 240
C-----------------------------------------------------------------------
C Work out Inet and sigma(Inet)
C-----------------------------------------------------------------------
BTOT = (BGRD1 + BGRD2)/(FRAC + FRAC)
BKN = BTOT/NP
TOP = COUNT - BTOT
BOT = SQRT(COUNT + BTOT/(FRAC + FRAC))
C-----------------------------------------------------------------------
C If GO mode and no profile analysis print results for non-standards
C-----------------------------------------------------------------------
IF (IPRFLG .NE. 0 .AND. KI(1:1) .EQ. 'G') THEN
IF (ISTAN .EQ. 0) THEN
ITOP = TOP + 0.5
IBOT = BOT + 0.5
IF (TOP .LE. SIGLIM*BOT) THEN
IF (ILPT .EQ. 0) WRITE (LPT,10000) IH,IK,IL,ITOP,IBOT,NREF
WRITE (COUT,10000) IH,IK,IL,ITOP,IBOT,NREF
ELSE
IF (NATT .NE. 0) THEN
IF (ILPT .EQ. 0)
$ WRITE (LPT,15100) IH,IK,IL,NATT,ITOP,IBOT,NREF
WRITE (COUT,15100) IH,IK,IL,NATT,ITOP,IBOT,NREF
ELSE
IF (ILPT .EQ. 0)
$ WRITE (LPT,15200) IH,IK,IL,ITOP,IBOT,NREF
WRITE (COUT,15200) IH,IK,IL,ITOP,IBOT,NREF
ENDIF
ENDIF
CALL GWRITE (ITP,' ')
GO TO 240
ENDIF
ENDIF
C-----------------------------------------------------------------------
C Test if peak is considered significant and print if not
C-----------------------------------------------------------------------
IF (TOP .LE. SIGLIM*BOT) IWARN = 1
IF (IWARN .NE. 0) THEN
ITOP = TOP + 0.5
IBOT = BOT + 0.5
IF (ILPT .EQ. 0) WRITE (LPT,10000) IH,IK,IL,ITOP,IBOT,NREF
WRITE (COUT,10000) IH,IK,IL,ITOP,IBOT,NREF
CALL GWRITE (ITP,' ')
ENDIF
IF (IWARN .NE. 0) GO TO 240
C-----------------------------------------------------------------------
C Profile is OK and significant. Print smoothed profile if Demo
C-----------------------------------------------------------------------
IF (KI .EQ. 'DE') THEN
WRITE (COUT,11000)
CALL GWRITE (ITP,' ')
WRITE (COUT,12000) (SMOOTH(J),J = 1,NP)
CALL GWRITE (ITP,' ')
ENDIF
C-----------------------------------------------------------------------
C Test that there are no funny bumps in the profile, by ensuring that
C the max of the peak is near the correct position.
C MAXA is the calculated position of the alpha peak
C MAXI is the intensity weighted maximum
C-----------------------------------------------------------------------
MAXA = AS*CON + RD12*A2 + 0.5
SUMI = 0
SUMNI = 0
DO 120 N = 1,NP
D = SMOOTH(N) - BKN
SUMI = SUMI + D
SUMNI = SUMNI + N*D
120 CONTINUE
MAXI = 0.5 + SUMNI/SUMI
C-----------------------------------------------------------------------
C Allow for a variable acceptance window
C-----------------------------------------------------------------------
CALL RSW(8,I)
ITOL = 5*I + ITOL
CALL RSW(7,I)
ITOL = 10*I + ITOL
CALL RSW(6,I)
ITOL = 20*I + ITOL
IF (ABS(MAXI-MAXA) .GT. ITOL) THEN
IF (TOP .GT. 2.0*SIGLIM*BOT) THEN
IWARN = 2
WRITE (COUT,14000) IH,IK,IL,MAXI,MAXA,BGRD1,COUNT,BGRD2
CALL GWRITE (ITP,' ')
IF (ILPT .EQ. 0)
$ WRITE (LPT,14000) IH,IK,IL,MAXI,MAXA,BGRD1,COUNT,BGRD2
ELSE
WRITE (COUT,14100) IH,IK,IL,MAXI,MAXA,BGRD1,COUNT,BGRD2
CALL GWRITE (ITP,' ')
IF (ILPT .EQ. 0)
$ WRITE (LPT,14100) IH,IK,IL,MAXI,MAXA,BGRD1,COUNT,BGRD2
ENDIF
GO TO 240
ENDIF
C-----------------------------------------------------------------------
C The profile is suitable for analysis to find the limits
C J1 is the beginning of the low angle search
C J2 is the beginning of the high angle search
C-----------------------------------------------------------------------
C J1 = MAXI - STEPOF*CON*AS - A2*ID12
C J2 = MAXI + STEPOF*CON*CS + A1*ID12
J1 = MAXI - ((STEPOF*AS)/STEP) - A2*ID12
J2 = MAXI + ((STEPOF*CS)/STEP) + A1*ID12
IF (J1 .LE. NWIND .OR. J2 .GE. NP-NWIND) THEN
ILOW = 1
IHIGH = NP
GO TO 210
ENDIF
C-----------------------------------------------------------------------
C Find the low angle limit by moving down from J1
C Set the window width to 0.67*0.67*CNT/5
C Find how many of the next NWIND values are in the window and if more
C than half are in the window, switch on the detector PROB.
C-----------------------------------------------------------------------
J = J1
LIM = J - 1
IFLAG = 0
PROB = 1
DO 160 I = NWIND,LIM
CNT = SMOOTH(J)
W = 0.08978*CNT
SUM = 0
DO 150 KK = J-NWIND,J-1
DIFF = CNT - SMOOTH(KK)
DC = DIFF*DIFF
IF (DC .LT. W) SUM = SUM + 1
150 CONTINUE
IF (SUM .GE. NWIND/2) IFLAG = 1
IF (IFLAG .NE. 0) THEN
PROB = PROB*(NWIND - SUM)/NWIND
IF (PROB .LE. PLIM) GO TO 170
ENDIF
J = J - 1
160 CONTINUE
170 ILOW = J-NWIND
IF (ILOW .LE. 0) ILOW = 1
C-----------------------------------------------------------------------
C Do the same for the high angle side
C-----------------------------------------------------------------------
J = J2
LIM = J + 1
IFLAG = 0
PROB = 1
DO 190 I = LIM,IDEL-NWIND
CNT = SMOOTH(J)
W = 0.08978*CNT
SUM = 0
DO 180 KK = J+1,J+NWIND
DIFF = CNT - SMOOTH(KK)
DC = DIFF*DIFF
IF (DC .LT. W) SUM = SUM + 1
180 CONTINUE
IF (SUM .GE. NWIND/2) IFLAG = 1
IF (IFLAG .NE. 0) THEN
PROB = PROB*(NWIND - SUM)/NWIND
IF (PROB .LE. PLIM) GO TO 200
ENDIF
J = J + 1
190 CONTINUE
200 IHIGH = J + NWIND
IF (IHIGH .GT. NP) IHIGH = NP
C-----------------------------------------------------------------------
C Now work out the net count & esd for profile between
C ILOW & IHIGH, using BGRD1 & BGRD2 plus pts between 1 to ILOW
C and IHIGH to NP for the background
C Revised EJG Aug 94 to allow for sloping backgrounds better
C-----------------------------------------------------------------------
210 NPK = IHIGH - ILOW + 1
B1 = BGRD1
IF (ILOW .GT. 1) THEN
DO 220 I = 1,ILOW-1
B1 = B1 + ACOUNT(I)
220 CONTINUE
ENDIF
FRAC1 = (FRAC*NP + ILOW - 1)/NPK
PEAK = 0.0
DO 225 I = ILOW,IHIGH
PEAK = PEAK + ACOUNT(I)
225 CONTINUE
B2 = BGRD2
IF (IHIGH .LT. NP) THEN
DO 230 I = IHIGH+1,NP
B2 = B2 + ACOUNT(I)
230 CONTINUE
ENDIF
FRAC2 = (FRAC*NP + NP - IHIGH)/NPK
BTOT = 0.5*(B1/FRAC1 + B2/FRAC2)
TOP1 = PEAK - BTOT
BOT1 = SQRT(PEAK + 0.25*(B1/(FRAC1*FRAC1) + B2/(FRAC2*FRAC2)))
FRAC1 = 0.5*(FRAC1 + FRAC2)
BGRD1 = BTOT*FRAC1
SUM = PEAK
BGRD2 = BGRD1
C-----------------------------------------------------------------------
C Print Inet and sigma(Inet) for non-standards in GO mode
C-----------------------------------------------------------------------
IF (KI(1:1) .EQ. 'G' .AND. ISTAN .EQ. 0) THEN
ITOP = TOP1 + 0.5
IBOT = BOT1 + 0.5
IF (TOP .LE. SIGLIM*BOT) THEN
IF (ILPT .EQ. 0) WRITE (LPT,10000) IH,IK,IL,ITOP,IBOT,NREF
WRITE (COUT,10000) IH,IK,IL,ITOP,IBOT,NREF
ELSE
IF (NATT .NE. 0) THEN
IF (ILPT .EQ. 0)
$ WRITE (LPT,15100) IH,IK,IL,NATT,ITOP,IBOT,NREF
WRITE (COUT,15100) IH,IK,IL,NATT,ITOP,IBOT,NREF
ELSE
IF (ILPT .EQ. 0)
$ WRITE (LPT,15200) IH,IK,IL,ITOP,IBOT,NREF
WRITE (COUT,15200) IH,IK,IL,ITOP,IBOT,NREF
ENDIF
ENDIF
CALL GWRITE (ITP,' ')
ENDIF
240 CALL RSW(9,JSW)
C------- always write profile at TRICS!
C IF (JSW .NE. 0 .and. istan .ne. 0) CALL PRFWRT (NP)
CALL PRFWRT (NP)
C-----------------------------------------------------------------------
C Prepare the profile for display on the c.r.t. if wanted
C Code below here is not needed for profile analysis
C The display is 10-bits * 10-bits
C If this reflection is to be plotted, the scaling is done in the
C display routine itself as the profile is developed.
C If the last reflection is to be plotted, the scaling is done here
C and an origin offset is added. Scaling is to a max of 1000 in each
C direction and the packing is
C 4096*scaled-counts + scaled-width + 4096*1024
C The marks are shifted by 100 points.
C
C SR 0 = 0 for normal display; = 1 for profile display
C-----------------------------------------------------------------------
CALL RSW (1,I)
IF (I .NE. 0) THEN
C-----------------------------------------------------------------------
C SR 1 = 0 not this time; = 1 for last reflection
C-----------------------------------------------------------------------
N = NXPTS
C-----------------------------------------------------------------------
C SR 2 = 0 for raw counts; = 1 for smoothed counts
C-----------------------------------------------------------------------
CALL RSW (2,J)
C-----------------------------------------------------------------------
C Insert marks at ILOW,IHIGH and ALPHA1 obs and calc positions
C-----------------------------------------------------------------------
IHTAGS(1) = AS * CON
IHTAGS(2) = AS * CON
IF (IWARN .NE. 1) IHTAGS(1) = MAXI - A2*ID12
IHTAGS(3) = ILOW
IHTAGS(4) = IHIGH
IF (J .NE. 0) THEN
CALL PTPREP (NP,SMOOTH,IHTAGS)
ELSE
CALL PTPREP (NP,ACOUNT,IHTAGS)
ENDIF
ENDIF
CALL RSW (3,J)
IF (J .EQ. 1) THEN
C-----------------------------------------------------------------------
C Dump the difference profile for Ladge
C-----------------------------------------------------------------------
C ic = 0
C do 1000 i = 1,np
C j = acount(i) + 0.5
C id(i) = j - ic
C ic = j
C 1000 continue
C WRITE (LPT,17100) (id(I),I=1,NP)
C17100 format (10(3x,z4))
WRITE (LPT,17000) (acount(I),I=1,NP)
WRITE (LPT,17000) (SMOOTH(I),I=1,NP)
ENDIF
RETURN
10000 FORMAT (3I4,2X,I7,'(',I4,') ',I5,' **')
11000 FORMAT (/,' The Profile counts are:')
12000 FORMAT (1X,10F7.0)
14000 FORMAT (3I4,' Max Profile',I4,', Alpha',I5,3F7.0)
14100 FORMAT (3I4,' Max Profile',I4,', Alpha',I4,3F7.0,' Weak Peak')
15000 FORMAT (3I4,F5.0,F7.0,F5.0,3I4,5F7.0/1X,F5.0,F8.4,2F6.0)
15100 FORMAT (3I4,I2,I7,'(',I4,') ',I5)
15200 FORMAT (3I4,2X,I7,'(',I4,') ',I5)
17000 FORMAT (1X,10F7.0)
END
C-----------------------------------------------------------------------
C Write a profile on unit 7 (32 4-byte variables per record) :--
C Each reflection is written as several records.
C Record 1:
C Bytes Symbol Contents
C 1 to 12 IH IK IL h, k, l 4 bytes each
C 13 to 16 NP2 number of pts in profile + 1000*std #
C 17 to 20 ILOW the point number on the low angle side
C 1 if no analysis
C 21 to 24 IHIGH the point number on the high angle side
C NP if no analysis
C 25 to 28 FRAC1 b/P time ratio (0.1 if no analysis)
C 29 to 32 IB1 Low angle background
C 31 to 36 ICOUNT Sum of all NP profile points
C 37 to 40 IB2 High angle background
C 41 to 28 44 profile points - 32000 (2 bytes each)
C
C Record 2 on:
C 1 to 128 64 profile points
C-----------------------------------------------------------------------
SUBROUTINE PRFWRT (NP)
INCLUDE 'COMDIF'
INTEGER*2 IPTS(500)
EQUIVALENCE (ACOUNT(501),IPTS(1))
NP2 = NP2 + 1000*NN
IB1 = BGRD1
ICOUNT = COUNT
IB2 = BGRD2
NREC = (NP + 20 + 63)/64 - 1
DO 100 I = 1,NP
IPTS(I) = ACOUNT(I) - 32000
100 CONTINUE
IPR = IOUNIT(7)
IDREC = 32*IBYLEN
STATUS = 'DO'
CALL IBMFIL (PRNAME, IPR,IDREC,STATUS,IERR)
NPR = NPR + 1
WRITE (IPR,REC=NPR) IH,IK,IL,NP2,ILOW,IHIGH,FRAC1,IB1,ICOUNT,IB2,
$ (IPTS(J),J=1,44)
IF (NREC .NE. 0) THEN
J1 = 45
DO 110 I = 1,NREC
J2 = J1 + 63
NPR = NPR + 1
WRITE (IPR,REC=NPR) (IPTS(J),J=J1,J2)
J1 = J2 + 1
110 CONTINUE
ENDIF
CALL IBMFIL (PRNAME,-IPR,IDREC,STATUS,IERR)
RETURN
END
C-----------------------------------------------------------------------
C Routine to write the binary stored profiles to an ASCII file
C The format of the ASCII file for each reflection is :--
C Line 1
C h,k,l, Npts, Ilow, Ihigh, Frac, Ib1, Icount, Ib2
C ( 3I4, 3I5, F8.5, I6, I7, I6)
C NREC lines of IPTS (10I6)
C-----------------------------------------------------------------------
SUBROUTINE PROFAS
INCLUDE 'COMDIF'
DIMENSION JPTS(500)
INTEGER*2 IPTS(500)
CHARACTER ASPROF*40
EQUIVALENCE (ACOUNT(501),IPTS(1)),(ACOUNT(1001),JPTS(1))
IPR = IOUNIT(7)
IDREC = 32*IBYLEN
CALL IBMFIL (PRNAME, IPR,IDREC,'DO',IERR)
IAS = IOUNIT(8)
WRITE (COUT,10000)
ASPROF = 'DONT DO IT'//' '
CALL ALFNUM (ASPROF)
IF (ASPROF .EQ. ' ') ASPROF = 'PROFL7.ASC'
CALL IBMFIL (ASPROF, IAS,IDREC,'SU',IERR)
NPR = 0
100 NPR = NPR + 1
READ (IPR,REC=NPR,IOSTAT=I)
$ IH,IK,IL,NP2,ILOW,IHIGH,FRAC,IB1,ICOUNT,IB2,(IPTS(J),J=1,52)
IF (I .EQ. 0) THEN
NP = NP2 - 1000*(NP2/1000)
NREC = (NP + 20 + 63)/64 - 1
IF (NREC .GT. 0) THEN
J1 = 45
DO 110 I = 1,NREC
J2 = J1 + 63
NPR = NPR + 1
READ (IPR,REC=NPR) (IPTS(J),J=J1,J2)
J1 = J2 + 1
110 CONTINUE
ENDIF
DO 120 I = 1,NP
JPTS(I) = IPTS(I) + 32000
120 CONTINUE
WRITE (IAS,11000) IH,IK,IL,NP2,ILOW,IHIGH,FRAC,IB1,ICOUNT,IB2
WRITE (IAS,12000) (JPTS(I),I=1,NP)
GO TO 100
ENDIF
CALL IBMFIL (PRNAME,-IPR,IDREC,'DO',IERR)
CALL IBMFIL (ASPROF,-IAS,IDREC,'SU',IERR)
KI = ' '
RETURN
10000 FORMAT (' Type the name of the ASCII file (PROFL7.ASC) ',$)
11000 FORMAT (3I4,3I5,F8.5,I6,I7,I6)
12000 FORMAT (10I6)
END
Reference in New Issue View Git Blame Copy Permalink