Aniso cutoff

A program to cut your data anisotropically is given below. To use it, save it as aniso_cutoff.f90 and compile&link with e.g.

gfortran -C aniso_cutoff.f90 -o aniso_cutoff

Then run it with ./aniso_cutoff

! program to impose elliptic resolution cutoff
! Kay Diederichs 2010-04-16
! input file: INTEGRATE.HKL
! output file: INTEGRATE.HKL.modified       - to be used for the CORRECT step
! user input: resolution cutoffs (in A) along a*, b* and c* (3 numbers)
!
module cell_mod
     real as,bs,cs,cosgs,cosbs,cosas
contains
subroutine getcell(a,b,c,alf,bet,gam)
     implicit none
     real, intent(in) :: a,b,c,alf,bet,gam
     real cosa,cosb,cosg,sina,sinb,sing,vol

     cosa=cos(alf/57.2958)
     cosb=cos(bet/57.2958)
     cosg=cos(gam/57.2958)
     sina=sin(alf/57.2958)
     sinb=sin(bet/57.2958)
     sing=sin(gam/57.2958)
     cosas=(cosb*cosg-cosa)/(sinb*sing)
     cosbs=(cosa*cosg-cosb)/(sina*sing)
     cosgs=(cosa*cosb-cosg)/(sina*sinb)
     vol=a*b*c*sqrt(1.-cosa**2-cosb**2-cosg**2+2*cosa*cosb*cosg)
     as=b*c*sina/vol
     bs=a*c*sinb/vol
     cs=a*b*sing/vol
end subroutine getcell
!
subroutine getres(ih,ik,il,resol)
     implicit none
     integer, intent(in) :: ih,ik,il
     real, intent(out) :: resol
     real s2,stol

     s2=((ih*as)**2+(ik*bs)**2+(il*cs)**2+2.*ih*ik*as*bs*cosgs         &
                   +2.*ih*il*as*cs*cosbs+2.*ik*il*bs*cs*cosas)/4.
     stol = sqrt(s2)
     resol=1./(2.*stol)
end subroutine getres
end module cell_mod
!
program main
     USE cell_mod
     implicit none
     real acut,bcut,ccut,a,b,c,alf,bet,gam,resol
     integer :: ih,ik,il,n1=0,n2=0
     character line*160

! open files
     open(1,file='INTEGRATE.HKL',status='OLD')
     open(2,file='INTEGRATE.HKL.modified',status='UNKNOWN')

     print*,'what is the resolution cutoff along a*, b*, c* (in Angstrom)?'
     print*,'for trigonal, hexagonal and tetragonal spacegroups, the first 2 numbers should be the same'
     print*,'for cubic spacegroups, all 3 numbers should be the same'
     read*,acut,bcut,ccut

! loop over lines of INTEGRATE.HKL
     do
       read(1,'(a)',end=99) line
       if (line(1:21)=='!UNIT_CELL_CONSTANTS=') then
         read(line(22:),*) a,b,c,alf,bet,gam
         print '(a,6f11.3)','!UNIT_CELL_CONSTANTS=',a,b,c,alf,bet,gam
         a=a/acut
         b=b/bcut
         c=c/ccut
         call getcell(a,b,c,alf,bet,gam)
       endif
       if (line(1:1)=='!') then
         write(2,'(a)') line(:len_trim(line))
       else
         read(line,*) ih,ik,il
         call getres(ih,ik,il,resol)
         if (resol>=1.) then
           write(2,'(a)') line(:len_trim(line))
           n2=n2+1
         endif
         n1=n1+1
       endif
     end do

! finish
99   print*,n1,'reflections found in INTEGRATE.HKL'
     print*,n2,'reflections written to INTEGRATE.HKL.modified'
     print*,'you should now enter these two commands:'
     print*,'mv INTEGRATE.HKL INTEGRATE.HKL.original'
     print*,'mv INTEGRATE.HKL.modified INTEGRATE.HKL'
     print*,'and then run the CORRECT job of XDS'
     close(1)
     close(2)
end program main