! KD 21.3.2018
! check A,B,C,F,I,R,H centering of XDS_ASCII.HKL
! see http://reference.iucr.org/dictionary/Integral_reflection_conditions
!
! ifort -c symops.f90 resolution.f90
! ifort -check all,noarg_temp_created -traceback -static-intel -sox checkcentering.f90 resolution.o symops.o -o checkcentering


MODULE reflection_mod
      INTEGER, ALLOCATABLE :: hkl(:,:)
      REAL, ALLOCATABLE :: intensity(:)
      INTEGER nref
END MODULE reflection_mod

PROGRAM checkcentering
      USE reflection_mod, ONLY: hkl,nref
      USE resolution_mod, ONLY: a,b,c,alf,bet,gam,getres,setbin,getlim,getbin
      IMPLICIT NONE
      INTEGER, PARAMETER :: nbin=10
      INTEGER :: hp(21)=(/143,144,145,149,150,151,152,153,154,168,169,170,171,172,173,&
                         177,178,179,180,181,182/)
      INTEGER :: ispgr,nc(2,nbin,8),i,j,itype,ibin,ncsum(2,8)
      REAL :: xc(2,nbin,8),s2,resol,resmin,resmax,reslim(nbin),xcsum(2,8)
      REAL, PARAMETER :: cut=0.5,deg2rad=3.141592654/180.
      CHARACTER :: centr*2(8)=(/'A ','B ','C ','F ','I ','Rr','Ro','H '/)

    1 FORMAT(a)

      WRITE(*,1)'Checking integral reflection conditions for XDS_ASCII.HKL - see Chapter 2.1.3.13 of International Tables for'
      WRITE(*,1)'Crystallography, Volume A, and http://reference.iucr.org/dictionary/Reflection_conditions (open access).'
      WRITE(*,1)'This program currently only implements checking of A B C F I R H centering of data processed in the P lattice.'
      WRITE(*,*)
      
      CALL read_reflection_file(ispgr,a,b,c,alf,bet,gam)

! find min and max resolution
      resmin=HUGE(resmin)
      resmax=0.
      DO i=1,nref
        CALL getres(hkl(1,i),hkl(2,i),hkl(3,i),s2,resol)
        IF (resol<resmin) resmin=resol
        IF (resol>resmax) resmax=resol      
      END DO
      CALL setbin(nbin,resmax,resmin,1)
      CALL getlim(reslim)
!      WRITE(*,'(a,99f7.3)')'resolution limits:',reslim

! create tables
      xc=0.
      nc=0
      DO i=1,nref
        CALL getres(hkl(1,i),hkl(2,i),hkl(3,i),s2,resol)
        CALL getbin(s2,ibin)
        CALL centered(i,hkl(1,i),hkl(2,i),hkl(3,i),ibin,nbin,xc,nc)
      END DO

      DO i=1,8
        xcsum(:,i)=SUM(xc(:,:,i),DIM=2)
        ncsum(:,i)=SUM(nc(:,:,i),DIM=2) 
        CALL get_ratio(xcsum(:,i),ncsum(:,i))
      END DO
      
      WRITE(*,*)
      WRITE(*,1)'Average intensity <I>, relative to mean, of reflections fulfilling ("present")'
      WRITE(*,1)'or violating ("absent") a given condition, and their numbers #:'
      WRITE(*,1)'centering type       A        B        C        F        I    R_revers  R_obvers    H'
      WRITE(*,1)'                                             all of'
      WRITE(*,1)'condition:         k+l=2n   h+l=2n   h+k=2n   A,B,C  h+k+l=2n k-h+l=3n  h-k+l=3n  h-k=3n'
      WRITE(*,'(a,8f9.2)')'<I> present:  ',xcsum(1,:)
      WRITE(*,'(a,8f9.2)')'<I> absent :  ',xcsum(2,:)
      WRITE(*,'(a,8i9)')  '# present:    ',ncsum(1,:)
      WRITE(*,'(a,8i9)')  '# absent :    ',ncsum(2,:)
      WRITE(*,*)
      WRITE(*,1)'Result of analysis:'
      IF (ANY(ncsum(2,:)==0)) THEN
        WRITE(*,1) 'One or more of the "absent" possibilities has zero reflections. This means that the'
        STOP 'data were processed in a centered setting. However, the program expects primitive (P) data.'
      END IF

      itype=MINLOC(xcsum(2,:),DIM=1)    ! find centering with lowest "absent" intensity
      IF (xcsum(2,4)<=cut) itype=4     ! F centering includes A, B, C centering
! early stop if no centering:
      IF (xcsum(2,itype)>cut) THEN
        WRITE(*,'(a,f5.2)') 'None of the centering possibilities has average "absent" intensity <',cut
        STOP 'Thus, no centering is suspected.'
      END IF
! in case of possible centering:
      WRITE(*,'(a,f7.3,a,f5.2,a,a,a)') 'Relative intensity=',xcsum(2,itype),' <',cut,' (an arbitrary cutoff) for ',TRIM(centr(itype)),' centering'
      WRITE(*,'(a,a,a)')'The data appear to be (pseudo-)centered. Analysing ',TRIM(centr(itype)),' centering by resolution:'
      WRITE(*,1)'reslim    rel.intensity     reflections'
      WRITE(*,1)'        present  absent  present  absent'
      DO i=1,nbin
        CALL get_ratio(xc(:,i,itype),nc(:,i,itype))
        WRITE(*,'(f7.3,2f8.3,1x,2i8)')reslim(i),xc(:,i,itype),nc(:,i,itype)
      END DO
      WRITE(*,*)
      WRITE(*,1)'If IDXREF indexes SPOT.XDS in a P lattice and you suspect it is actually centered based'
      WRITE(*,1)'on the analysis above, you may inspect the "absent" reflections in SPOT.XDS with e.g.:'
      WRITE(*,*)
      IF (itype==1) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($6+$7)%2\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==2) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$7)%2\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==3) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$6)%2\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==4) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (\!(($5+$6)%2==0 && ($5+$7)%2==0 && ($6+$7)%2==0)) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==5) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$6+$7)%2\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==6) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($6-$5+$7)%3\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==7) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5-$6+$7)%3\!=0) print $0}'' SPOT.XDS | less'
      ELSE IF (itype==8) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5-$6)%3\!=0) print $0}'' SPOT.XDS | less'
      END IF
      WRITE(*,*)
      WRITE(*,1)'To make IDXREF index the data in a centered lattice, filter out those reflections from SPOT.XDS'
      WRITE(*,1)'that violate the condition, and re-run IDXREF with SPACE_GROUP_NUMBER=0. Filtering means e.g:'
      WRITE(*,*)
      WRITE(*,1)'cp -p SPOT.XDS SPOT.XDS.orig'
      IF (itype==1) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($6+$7)%2==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==2) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$7)%2==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==3) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$6)%2==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==4) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$6)%2==0 && ($5+$7)%2==0 && ($6+$7)%2==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==5) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5+$6+$7)%2==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==6) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($6-$5+$7)%3==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==7) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5-$6+$7)%3==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      ELSE IF (itype==8) THEN
        WRITE(*,1)'awk ''\!/ 0 0 0/{if (($5-$6)%3==0) print $0}'' SPOT.XDS.orig >SPOT.XDS'
      END IF
! some cases where the cell changes are explicitly told to the user (complete the list as needed!):
      IF ((ispgr==3.OR.ispgr==4).AND.itype==5) THEN
        WRITE(*,*)
        bet=bet*deg2rad
        s2=sqrt(a*a+c*c+2*a*c*cos(bet))
        bet=acos(c*cos(bet)/s2) / deg2rad
        WRITE(*,1)'Or, to force IDXREF indexing in C2 instead of I-centered P2(1) (as found here), run with'
        WRITE(*,'(a,3f7.1,3f6.1)')'SPACE_GROUP_NUMBER=5 and UNIT_CELL_CONSTANTS=',s2,b,a,alf,bet,gam
      ELSE IF (ANY(hp==ispgr).AND.itype==8) THEN
        WRITE(*,*)
        WRITE(*,1)'Or, to force IDXREF indexing in (shorter-axis) hP instead of H-centered hP (as found here),'
        WRITE(*,'(a,i0,a,3f7.1,3f6.1)')'run with SPACE_GROUP_NUMBER= ',ispgr,' and UNIT_CELL_CONSTANTS=',&
                                                a/SQRT(3.),b/SQRT(3.),c,alf,bet,gam
      ELSE IF ((ispgr==16.OR.ispgr==17.OR.ispgr==18.OR.ispgr==19).AND.itype==1) THEN
        WRITE(*,*)
        WRITE(*,1)'Or, to force IDXREF indexing in oC instead of A-centered oP (as found here),'
        WRITE(*,'(a,3f7.1,3f6.1)')'run with SPACE_GROUP_NUMBER= 21 or 22, and UNIT_CELL_CONSTANTS=',&
                                                b,c,a,alf,bet,gam
      ELSE IF ((ispgr==16.OR.ispgr==17.OR.ispgr==18.OR.ispgr==19).AND.itype==2) THEN
        WRITE(*,*)
        WRITE(*,1)'Or, to force IDXREF indexing in oC instead of B-centered oP (as found here),'
        WRITE(*,'(a,3f7.1,3f6.1)')'run with SPACE_GROUP_NUMBER= 21 or 22, and UNIT_CELL_CONSTANTS=',&
                                                c,a,b,alf,bet,gam
      END IF
      
END PROGRAM checkcentering
!
SUBROUTINE centered(i,h,k,l,ibin,nbin,xc,nc)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: i,h,k,l,ibin,nbin
      REAL, INTENT(INOUT) :: xc(2,nbin,8)
      INTEGER, INTENT(INOUT) :: nc(2,nbin,8)
      
      CALL assign(i,MOD(k+l,2)==0,  xc(:,ibin,1),nc(:,ibin,1)) ! A 
      CALL assign(i,MOD(h+l,2)==0,  xc(:,ibin,2),nc(:,ibin,2)) ! B 
      CALL assign(i,MOD(h+k,2)==0,  xc(:,ibin,3),nc(:,ibin,3)) ! C 
      CALL assign(i,MOD(k+l,2)+MOD(h+l,2)+MOD(h+k,2)==0,xc(:,ibin,4),nc(:,ibin,4)) ! F 
      CALL assign(i,MOD(h+k+l,2)==0,xc(:,ibin,5),nc(:,ibin,5)) ! I 
      CALL assign(i,MOD(k-h+l,3)==0,xc(:,ibin,6),nc(:,ibin,6)) ! R reverse 
      CALL assign(i,MOD(h-k+l,3)==0,xc(:,ibin,7),nc(:,ibin,7)) ! R obverse 
      CALL assign(i,MOD(h-k,3)==0,  xc(:,ibin,8),nc(:,ibin,8)) ! H 
END SUBROUTINE centered

SUBROUTINE assign(i,truth,means,counts)
      USE reflection_mod, ONLY: intensity
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: i
      LOGICAL, INTENT(IN) :: truth
      REAL, INTENT(INOUT) :: means(2)
      INTEGER, INTENT(INOUT) :: counts(2)
      
      IF (truth) THEN
        means(1)=means(1)+intensity(i)
        counts(1)=counts(1)+1
      ELSE
        means(2)=means(2)+intensity(i)
        counts(2)=counts(2)+1
      END IF
END SUBROUTINE assign

SUBROUTINE get_ratio(means,counts)
      IMPLICIT NONE
      REAL, INTENT(INOUT) :: means(2)
      INTEGER, INTENT(INOUT) :: counts(2)
      
      WHERE (counts>0) means=means/counts
      IF (SUM(means)>0.) means=2*means/SUM(means)
END SUBROUTINE get_ratio

SUBROUTINE read_reflection_file(ispgr,a,b,c,alf,bet,gam)
      USE symops_mod, ONLY: getops,asuput,op,ncenop
      USE reflection_mod, ONLY: hkl,intensity,nref
      IMPLICIT NONE
      INTEGER, INTENT(OUT) :: ispgr
      REAL, INTENT(OUT) :: a,b,c,alf,bet,gam
      INTEGER, ALLOCATABLE :: t_hkl(:,:)
      REAL, ALLOCATABLE :: t_intensity(:)
      INTEGER :: ih,ik,il,oldhkl,jh,jk,jl,ihkl,iop,nobs,isize,jhold,jkold,jlold
      REAL :: fsq,sigma,xint,siga
      CHARACTER (LEN=80) :: line
      LOGICAL :: last

 1000 FORMAT(A)
      last=.FALSE.
      nobs=0
      nref=0
      xint=0.
      siga=0.
      ALLOCATE(hkl(3,1000),intensity(1000))

      OPEN(1,FILE='XDS_ASCII.HKL',ACTION='READ')    
! header of XDS_ASCII.HKL
      DO
        READ(1,1000) line
        IF (line(1:20)=='!SPACE_GROUP_NUMBER=') then
          READ(line(21:),*) ispgr
          CALL getops(ispgr)
          WRITE(*,1000) TRIM(line)
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
          WRITE(*,1000) TRIM(line)
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO

 ! read reflections 
      DO
        READ(1,1000) line
        IF (line(1:12) == '!END_OF_DATA') last=.TRUE. 
        IF (.NOT.last) THEN
          READ(line,*)ih,ik,il,fsq,sigma
!          IF (sigma<0.) CYCLE
          CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,SIZE(op,3)/ncenop)
          nobs=nobs+1
          IF (nobs==1) THEN
            oldhkl=ihkl  ! make sure old hkl is sensibly initialized
            jhold=jh
            jkold=jk
            jlold=jl
          END IF
        END IF
        IF (ihkl /= oldhkl.OR.last) THEN  ! process unique reflection
          IF (siga>0.) THEN
            nref=nref+1
            IF (nref>SIZE(intensity)) THEN  ! make more space
              isize=SIZE(intensity)
              ALLOCATE(t_hkl(3,2*isize))
              ALLOCATE(t_intensity(2*isize))
              t_hkl(:,:isize)=hkl
              t_intensity(:isize)=intensity
              CALL MOVE_ALLOC(t_hkl,hkl)
              CALL MOVE_ALLOC(t_intensity,intensity)
            END IF
            hkl(1,nref)=jhold
            hkl(2,nref)=jkold
            hkl(3,nref)=jlold
            intensity(nref)=xint/siga
          END IF
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          xint=0.
          siga=0.
          oldhkl=ihkl
          jhold=jh
          jkold=jk
          jlold=jl
        END IF  ! process unique reflection
! store this observation
        xint=xint + fsq/sigma**2
        siga=siga + 1./sigma**2
      END DO   ! the loop reading the XDS_ASCII.HKL file 
      WRITE(*,'(2(a,i8))')'nobs =',nobs,' nunique=',nref
END SUBROUTINE read_reflection_file
