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'

! loop over lines of INTEGRATE.HKL
do
if (line(1:21)=='!UNIT_CELL_CONSTANTS=') then
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
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
```