      module stuff_mod
        integer, parameter :: idim=2,nmax=32000 ! idim is dimension (2 or 3), nmax is max number of atoms
        integer :: natom
	real, allocatable :: disttble(:,:),wtble(:,:)
      end module stuff_mod

      program xds_stat
!      USE stuff_mod
! compile with: f90 -O  xdsstat.f90 pck-routines.o $CCP4_LIB/src/libhtml.o -L$CCP4_LIB -lccp4f -lccp4c -o xdsstat
! maybe add -lccif if the linker complains
!
! where f90 = ifort or g95 or gfortran -fno-second-underscore ! 32bit compiler!
!
! I use
! ifort -m32 -u -warn -C -traceback -static-intel -sox xdsstat.f90 pck-routines.o \
!	/usr/local/ccp4-5.0.2.ifort/lib/src/libhtml.o \
!     /usr/local/ccp4-5.0.2.ifort/lib/libccp4f.a \
!     /usr/local/ccp4-5.0.2.ifort/lib/libccp4c.a -o xdsstat
! 
! nota bene: I had to manually compile $CCP4_LIB/src/libhtml.f because it was not compiled already
!
! 31.10.2011:
! gfortran -m32 -fno-range-check xdsstat.f90 pck-routines.f90 -o xdsstat  -L$CCP4_LIB -lccp4f -lccp4c
! seems to compile well, but has a runtime error for reading the symmetry lib
!
! 6.11.2013 this works with a 32bit CCP4-6.3.0:
! gfortran -m32 -fno-range-check lbfgs_4+8.f90 xdsstat.f90 ../pck-routines.f90 -o xdsstat -Llib32 \
! -lccp4f -lccp4c -static-libgfortran -static-libgcc
!
! 6.11.2013 this works for 64bit 6.3.0:
! gfortran -fno-range-check lbfgs_4+8.f90 xdsstat.f90 ../pck-routines.f90 -o xdsstat -Llib64 \
! -lccp4f -lccp4c -lmmdb -static-libgfortran -static-libgcc
!
! 26.3.2015 on CentOS7 after scl enable devtoolset-3 bash:
! gfortran -fno-range-check lbfgs_4+8.f90 xdsstat.f90 ../pck-routines.f90 -o xdsstat \
! -Llib64 -lccp4f -lccp4c -lmmdb -lgfortran -static-libgfortran -static-libgcc \
! /opt/rh/devtoolset-3/root/usr/lib/gcc/x86_64-redhat-linux/4.9.1/libquadmath.a 
! the last line is not needed on SL6
! for 32bit use -m32 and -Llib32
! 
! Kay Diederichs 3/2001
! to give image-wise statistics from (primarily) XDS_ASCII.HKL
!
! the program assumes that unique indices occur in groups in the input file
!
! 12.12.2007 increase maxfrm to 10000; remove long output
! 28.4.2008 fix output loggraph
! 16.2.2010 increase maxfrm to 50000
! 14.5.2010 increase maxstore from 200 to 2000
! 4.9.2010 change nint to int in "iframe=" statement
! 31.10.2011 make maxfrm arrays allocatable
! 16.8.2013 implement CC
! 9.3.2019 bugfix line 819 had been "if (iyy<1.OR.iyy>nx4) cycle" (should be ny4)
! 27.8.21 allow unlimited number of header lines (instead of 50) for XSCALE
      implicit none
      integer :: maxfrm=1000,oldmax,maxitem,maxcol,maxstore,maxsym
! maximum of maxstore equivalent reflections
! when changing maxfrm also change line 280 ! similar for maxstore
      parameter (maxitem=11,maxcol=20)
      parameter (maxstore=2000,maxsym=32)
      integer i,j,iframe,nitemfound,nitem,colindex(maxitem)
      integer, allocatable :: nref(:),misfits(:) ! nref(0:maxfrm),misfits(0:maxfrm)
      integer isymstore(maxstore),hklstore(3,maxstore),frmstore(maxstore)
      integer, allocatable :: nrfac(:),nuniq(:)  ! nrfac(0:maxfrm),nuniq(0:maxfrm)
      real refstore(maxstore),xystore(2,maxstore),rlpstore(maxstore)
      real sigstore(maxstore)
      integer peakstore(maxstore),corrstore(maxstore)
      real refmisstore(maxstore)
      integer xymisstore(2,maxstore)
      integer ispgrp,ihkl(3),jhkl(3),isym,nsym,mlaue,msymp,oldhkl(3)
      integer in1,in2,out1,lrec,ifail,nstore,ntotuniq
      real RlSymmMatrx(4,4,192)
      character line*80,spgnam*80,pgname*80
      real reflex(maxcol),sumref,total(maxitem)
      real , allocatable :: stat(:,:),rfac(:),rfacsum(:)  ! stat(maxitem,0:maxfrm),rfac(0:maxfrm),rfacsum(0:maxfrm)
      real zaehler, nenner, d,rmeas,rmeastot,rmeastotsum
      real , allocatable :: deltanofrzae(:),deltanofrnen(:) ! deltanofrzae(0:maxfrm),deltanofrnen(0:maxfrm)
      real, allocatable :: deltafrzae(:),deltafrnen(:) ! deltafrzae(0:maxfrm),deltafrnen(0:maxfrm)
      real, allocatable :: scalzae(:,:),scalnen(:,:)
      real, allocatable :: rlpzae(:,:),rlpnen(:,:)
      real, allocatable :: peakzae(:,:),peaknen(:,:)
      real, allocatable :: corrzae(:,:),corrnen(:,:)
      real, allocatable :: rfzae(:,:),rfnen(:,:)
      real, allocatable :: friedz(:,:),friedn(:,:)
      integer, parameter :: maxbin=12  ! if maxbin is changed two formats also have to be changed 
      real :: peak(50:100,maxbin)=0.,peaksum(50:100,maxbin)=0.
      integer :: npeak(50:100,maxbin)=0, ier
      integer*2, allocatable :: img2(:,:), img3(:,:)
      integer*2, allocatable :: ndet(:,:)
      real, allocatable :: miszae(:,:), misnen(:,:)
      
      integer, allocatable :: deltanofried(:),deltafried(:) ! deltanofried(0:maxfrm),deltafried(0:maxfrm)
      real, allocatable :: deltanominator(:),deltanumerator(:) ! deltanominator(0:maxfrm),deltanumerator(0:maxfrm)
      
      integer, allocatable :: delsymcnt(:,:)  ! delsymcnt(0:maxfrm,0:maxsym)
      real, allocatable :: delsymnom(:,:),delsymden(:,:) ! delsymnom(0:maxfrm,0:maxsym),delsymden(0:maxfrm,0:maxsym)
      real dsd(maxsym)
      integer dsc(maxsym)
      
      integer hklcnt(-512:512,3)
      real hkln(-512:512,3),hkld(-512:512,3)
      
      integer, allocatable :: deltacounter(:),tmp_i(:),tmp_i2(:,:)  ! deltacounter(0:maxfrm)
      real, allocatable :: tmp_r(:),tmp_r2(:,:),tmp_r3(:,:,:)
!      real, allocatable :: big_dist(:,:,:)
      integer nrefsum,misfitsum,nrfactot,nuniqtot,idiff,k,mi,ma,l,colcount
      character colitem(maxitem)*11
      logical first,hkleq
      integer nx,ny,zmin,zmax,ix,iy,iz,nx4,ny4,ixx,iyy,nmist
      real abszae(25,25,10),absnen(25,25,10)
      integer iargc
      real res,reslo,reshi,sthlsq,a,b,c,al,be,ga
! interesting columns in XDS_ASCII.HKL:
      data colitem/'H','K','L','ZD','IOBS','SIGMA(IOBS)','PEAK','CORR',  &
                   'XD','YD','RLP'/
      data first/.true./
      
! INTEGER arrays, rank 1
      allocate(nref(0:maxfrm))
      allocate(misfits(0:maxfrm))
      allocate(nrfac(0:maxfrm))
      allocate(nuniq(0:maxfrm))
      allocate(deltanofried(0:maxfrm))
      allocate(deltafried(0:maxfrm))
      allocate(deltacounter(0:maxfrm))
! INTEGER arrays, rank 2
      allocate(delsymcnt(0:maxfrm,0:maxsym)) 
! REAL arrays, rank 1
      allocate(rfac(0:maxfrm))
      allocate(rfacsum(0:maxfrm))
      allocate(deltanofrzae(0:maxfrm))
      allocate(deltanofrnen(0:maxfrm))
      allocate(deltafrzae(0:maxfrm))
      allocate(deltafrnen(0:maxfrm))
      allocate(deltanominator(0:maxfrm))
      allocate(deltanumerator(0:maxfrm))
! REAL arrays, rank 2
      allocate(stat(maxitem,0:maxfrm))
      allocate(delsymnom(0:maxfrm,0:maxsym))
      allocate(delsymden(0:maxfrm,0:maxsym))
! REAL array, rank 3
!      allocate(big_dist(6,maxfrm,maxfrm))
      
      write(0,*) 'XDSSTAT version 27.8.2021'
      write(0,*) 'The program can be run with two commandline ',         &
           'arguments: the lower and upper limit of the resolution range'
      write(0,*)'e.g. xdsstat 100 6   will use 100-6A data'
      write(0,*)' '
      write(0,*)'When using results from XDSSTAT please cite'
      write(0,*)'Diederichs, K.: Some aspects of quantitative analysis'
      write(0,*)'and correction of radiation damage.'
      write(0,*)'Acta Cryst (2006) D62:96-101'
      write(0,*)'Thank you!' 
      
      if (iargc()==2) then
! read hi and lo resolution limits from command line
        call getarg(1,line)
        print*,'arg1= ',line(1:len_trim(line))
        read(line,*) reslo
        call getarg(2,line)
        print*,'arg2= ',line(1:len_trim(line))
        read(line,*) reshi
        if (reslo<reshi) then ! reverse if given in wrong order
          a=reslo
          reslo=reshi
          reshi=a
        endif
      else
        reslo=99999.
        reshi=0.
      endif
      hklcnt=0
      hkln=0.
      hkld=0.
      abszae=0.
      absnen=0.
      
      delsymcnt = 0
      delsymnom = 0.
      delsymden = 0.
      deltanominator = 0.
      deltanumerator = 0.
      deltacounter   = 0
      deltafried     = 0
      deltanofried   = 0
      deltanofrzae   = 0
      deltanofrnen   = 0
      deltafried     = 0
      deltafrzae     = 0
      deltafrnen     = 0
!      big_dist       = 0
      stat = 0.
      total = 0.
      nref = 0
      nrfac = 0
      nuniq = 0
      nuniqtot = 0
      rmeastot = 0.
      rmeastotsum = 0.
      nrfactot = 0
      misfits = 0
      colindex = 0
      rfac = 0.
      rfacsum = 0.
      lrec = 0
      ifail = 0
      nrefsum = 0
      misfitsum = 0
      sumref = 0.
      nstore = 0
      nmist = 0
      out1 = 6
      in1  = 8
      in2  = 7
      nx = 0
      ny = 0
      zmin = 0
      zmax = 0
      

! initialize CCP4 and HTML output
      call ccpfyp()
      call ccp4h_init()
      call ccp4h_pre_beg()
!
      write(0,*) ' '
      write(0,*) 'advice: run the program e.g. as ...'
      write(0,*) 'xdsstat > XDSSTAT.LP'
      write(0,*) ' '
      write(0,*)'inputfile (from XDS/XSCALE)? (just hit ENTER for XDS_ASCII.HKL)'
      read(5,'(a)') line
      if (line=='') then
        line='XDS_ASCII.HKL'
      else
! get rid of leading blanks
!      line = adjustl(line)
  10    if (line(1:1)==' ') then
          line(1:)=line(2:)
          goto 10
        endif
      endif
      write(out1,'(a,a,/)') 'reading from ',line
! I couldn't get ccpdpn to work for me. KD 3/01
!      call ccpdpn(in2,line(1:lentrim(line)),'READONLY','F',lrec,ifail)
      open(in2,file=line,status='OLD',iostat=ifail)
      if (ifail > 0) stop 'file not found'
!      print*,line(1:lentrim(line)),' opened on unit ',in2
!   
! read header      
      nitemfound=0
      i=0
      do ! i=1,50   no limit 27.8.21   
        read (in2,'(a80)') line
        IF (line(1:1)=='!') i=i+1
        write(out1,*) line
print*,i,trim(line)
        if (line(1:37)=='!NUMBER_OF_ITEMS_IN_EACH_DATA_RECORD=') then
          read(line(38:),'(i4)') nitem
          write(out1,*)'number of items according to header is ',nitem
        elseif (line(1:21)=='!UNIT_CELL_CONSTANTS=') then
          read(line(22:),*) a,b,c,al,be,ga   ! '(3f9.3,3f8.3)'
          write(*,'(a,6f9.3)')'unit cell parameters are',a,b,c,al,be,ga
          call setrsl(a,b,c,al,be,ga)  ! prepare to calculate resolution
        elseif (line(1:6)=='!ITEM_') then
          do j=1,maxitem
            if (line(7:6+lentrim(colitem(j)))==colitem(j)) then
              read(line(8+lentrim(colitem(j)):),'(i2)') colindex(j)
              if (colindex(j) > maxcol) stop 'maxcol'
            endif
          end do
          nitemfound=nitemfound+1
        elseif (line=='!END_OF_HEADER') then
          goto 50
        elseif (line(1:4)=='!NX=') then
          read(line( 5:),'(i6)') nx
          read(line(16:),'(i6)') ny
          write(out1,*)'NX, NY =',nx,ny
          nx4 = nx/4
          ny4 = ny/4
          allocate(img2(nx4,ny4),img3(nx4,ny4))
          allocate(scalzae(nx4,ny4),scalnen(nx4,ny4))
          scalzae=0.
          scalnen=0.
          allocate(rlpzae(nx4,ny4),rlpnen(nx4,ny4))
          rlpzae=0.
          rlpnen=0.
          allocate(peakzae(nx4,ny4),peaknen(nx4,ny4))
          peakzae=0.
          peaknen=0.
          allocate(corrzae(nx4,ny4),corrnen(nx4,ny4))
          corrzae=0.
          corrnen=0.
          allocate(rfzae(nx4,ny4),rfnen(nx4,ny4))
          rfzae=0.
          rfnen=0.
          allocate(friedz(nx4,ny4),friedn(nx4,ny4))
          friedz=0.
          friedn=0.
          allocate(ndet(nx4,ny4),miszae(nx4,ny4),misnen(nx4,ny4))
          ndet = 0
          miszae = 0.
          misnen = 0.
        elseif (line(1:12)=='!DATA_RANGE=') then
          read(line(13:),'(2i8)') zmin,zmax
          write(out1,*) 'data_range=',zmin,zmax
        elseif (line(1:21)=='!SPACE_GROUP_NUMBER=') then
          read(line(21:),'(i5)') ispgrp
          print*,'the space group number is', ispgrp
          call ccpdpn(in1, 'SYMOP','READONLY','F',lrec,ifail)
          call msyget(in1,ispgrp,nsym,RlSymmMatrx)
          close(in1)
          call asuset(spgnam,ispgrp,pgname,nsym,RlSymmMatrx,msymp,         &
     &                mlaue,.false.)
! previous line: changing false to true is supposed to make asusset 
! print out the symops but this does not seem to work
        endif
      end do
      stop 'no !END_OF_HEADER line found?!'
!
! come here after reading header
!
  50  write(out1,*)'found ',i,' header lines with ',nitemfound,' items'
      if (nitem /= nitemfound) stop 'nitem /= nitemfound'
      if (colindex(4)==0) stop 'no ZD column given => no listing'
      write(out1,*)'the interesting columns and their items are:'
      do i=1,maxitem
      write(out1,'(i11,2x,a11)') colindex(i),colitem(i)
      end do
      if (nitem > maxcol) stop 'nitem > maxcol'
!     
! reflection-scanning loop
!      
      ntotuniq = 1 ! 1 because update_stuff is called in the end, too
 100  read(in2,*,err=999,end=998) (reflex(i),i=1,nitem)
      ihkl(1) = nint(reflex(colindex(1)))
      ihkl(2) = nint(reflex(colindex(2)))
      ihkl(3) = nint(reflex(colindex(3)))
      res = 1./sqrt(4.*sthlsq(ihkl(1),ihkl(2),ihkl(3)))
      if (res > reslo .OR. res < reshi) goto 100
! look into ZD column to find which frame this refl belongs to:
      iframe=int(reflex(colindex(4))) ! changed from nint() 4.9.2010
! NB int(-0.2)=int(0.2)=0
      iframe=iframe+1  ! 4.9.2010
      if (iframe > maxfrm) then
              oldmax=maxfrm
              maxfrm=max(maxfrm*2,iframe+999)
              print*,'so many frames ... increasing maxfrm to',maxfrm 
! INTEGER arrays, rank 1
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=nref;        tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,nref)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=misfits;     tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,misfits)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=nrfac;       tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,nrfac)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=nuniq;       tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,nuniq)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=deltanofried;tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,deltanofried)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=deltafried;  tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,deltafried)
      allocate(tmp_i(0:maxfrm));tmp_i(0:oldmax)=deltacounter;tmp_i(oldmax+1:)=0;call move_alloc(tmp_i,deltacounter)
! INTEGER arrays, rank 2
      allocate(tmp_i2(0:maxfrm,0:maxsym));tmp_i2(0:oldmax,:)=delsymcnt;tmp_i2(oldmax+1:,:)=0;call move_alloc(tmp_i2,delsymcnt)
! REAL arrays, rank 1
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=rfac          ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,rfac)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=rfacsum       ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,rfacsum)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltanofrzae  ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltanofrzae)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltanofrnen  ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltanofrnen)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltafrzae    ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltafrzae)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltafrnen    ;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltafrnen)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltanominator;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltanominator)
      allocate(tmp_r(0:maxfrm));tmp_r(0:oldmax)=deltanumerator;tmp_r(oldmax+1:)=0.;call move_alloc(tmp_r,deltanumerator)

! REAL arrays, rank 2
      allocate(tmp_r2(maxitem,0:maxfrm));tmp_r2(:,0:oldmax)=stat 
      tmp_r2(:,oldmax+1:)=0.;call move_alloc(tmp_r2,stat)
      allocate(tmp_r2(0:maxfrm,0:maxsym));tmp_r2(0:oldmax,:)=delsymnom
      tmp_r2(oldmax+1:,:)=0.;call move_alloc(tmp_r2,delsymnom)
      allocate(tmp_r2(0:maxfrm,0:maxsym));tmp_r2(0:oldmax,:)=delsymden
      tmp_r2(oldmax+1:,:)=0.;call move_alloc(tmp_r2,delsymden)      

! REAL arrays, rank 3
!      allocate(tmp_r3(6,maxfrm,maxfrm));tmp_r3=0.;tmp_r3(:,:oldmax,:oldmax)=big_dist;call move_alloc(tmp_r3,big_dist)

      endif
      if (iframe < 0) then
        iframe=0
        print*,'reflection with surprisingly low z-value:',reflex(colindex(4))
      endif
! ... and update the statistics for this frame (6 is sigma(Iobs))
      if (reflex(colindex(6)) < 0.) then
        misfits(iframe) = misfits(iframe) + 1
        misfitsum  = misfitsum  + 1
      else
        nref(iframe) = nref(iframe) + 1
        nrefsum = nrefsum + 1
        do j=5,8
          if (colindex(j) > 0) then
            stat(j,iframe) = stat(j,iframe) + reflex(colindex(j))
            total(j)  = total(j)  + reflex(colindex(j))
          endif
        end do
      endif
! get unique hkl
      call asuput(ihkl,jhkl,isym)
      if (first) then
        oldhkl = jhkl
        first = .false.
      endif
      if (.NOT.hkleq(jhkl,oldhkl)) then
! this is a different unique reflection than what we looked at so far
        call update_stuff
        ntotuniq = ntotuniq + 1
      endif
      if (colindex(9)>0) then
        ix = max(1,min(nx/4,nint(reflex(colindex( 9))/4)))
        iy = max(1,min(ny/4,nint(reflex(colindex(10))/4)))
      else
        ix = 0
      endif
      if (reflex(colindex(6)) >= 0.) then
        nstore = nstore + 1
        if (nstore.gt.maxstore) stop 'maxstore=2000 too low. recompile!'
        refstore(nstore) = reflex(colindex(5)) ! intensity
        sigstore(nstore) = reflex(colindex(6)) ! sigma                  
        frmstore(nstore) = iframe                   ! frame number
        isymstore(nstore)= isym
        hklstore(:,nstore) = ihkl
        if (colindex(9)>0) then
          xystore(1,nstore)= reflex(colindex(9))    ! XD
          xystore(2,nstore)= reflex(colindex(10))   ! YD
        endif
        if (colindex(11)>0) rlpstore(nstore) = reflex(colindex(11))
        if (colindex(7)>0)  peakstore(nstore) = nint(reflex(colindex(7)))
        if (colindex(8)>0)  corrstore(nstore) = nint(reflex(colindex(8)))
        sumref = sumref + reflex(colindex(5))
      else  ! also take misfits into account
        nmist = nmist + 1
        if (nmist.gt.maxstore) stop 'maxstore too low. fix & recompile!'
        refmisstore(nmist) = reflex(colindex(5)) ! intensity
        if (colindex(9)>0) then
          xymisstore(1,nmist)= ix   ! XD
          xymisstore(2,nmist)= iy   ! YD
        endif
      endif
      if (colindex(9)>0.AND.NX>0) then
        ndet(ix,iy) = ndet(ix,iy) + 1
      endif
      goto 100
! end of reflection-scanning loop

998   stop 'invalid file or !END_OF_DATA line is missing'
999   backspace(in2)
      read(in2,'(a)') line
      if (line/='!END_OF_DATA') then
        write(*,*) line
        stop 'read error'
      endif
! treat last reflection
      if (nstore>0) call update_stuff
! make output such that when viewed in netscape we get a Java plot window !
      call ccp4h_graph_beg(800,600)
      if (colindex(7)==0.or.colindex(8)==0) then
! no peak and corr plot
        write(out1,5)
      else
        write(out1,2)
      endif
    2 format(/,                                                          &
     & ' $TABLE: Statistics of data set(s):',/,                          &
     & ' $GRAPHS',/,                                                     &
     & ' :reflections, misfits, contributing_to_Rmeas, unique',/,        &
     & ' :A:1,2,3,10,11:',/,                                             &
     & ' :mean Intensity and sigma, and Intensity/sigma',/,              &
     & ' :A:1,4,5,6:',/,                                                 &
     & ' :mean fraction of theoretical intensity, and correlation',      &
     & ' with standard profiles',/,                                      &
     & ' :A:1,7,8:',/,                                                   &
     & ' :individual Rmeas of all frames',/,                             &
     & ' :A:1,9:',/,                                                     &
     & ' $$',/                                                           &
     & ' Frame #refs #misfits  Iobs     sigma   Iobs/sigma ',            &
     &      '  Peak      Corr     Rmeas  #rmeas #unique  dummy $$',      &
     &/,' $$',/)
    5 format(/,                                                          &
     & ' $TABLE: Statistics of data set(s):',/,                          &
     & ' $GRAPHS',/,                                                     &
     & ' :individual Rmeas of all frames',/,                             &
     & ' :A:1,9:',/,                                                     &
     & ' :reflections, misfits, contributing_to_Rmeas, unique',/,        &
     & ' :A:1,2,3,10,11:',/,                                             &
     & ' :mean Intensity and sigma, and Intensity/sigma',/,              &
     & ' :A:1,4,5,6:',/,                                                 &
     & ' $$',/                                                           &
     & ' Frame #refs #misfits  Iobs     sigma   Iobs/sigma ',            &
     &      '  Peak      Corr     Rmeas  #rmeas #unique  dummy $$',      &
     &/,' $$',/)
      do i=0,maxfrm
        if (nref(i)>0.or.misfits(i)>0) then
        if (nref(i)>0) then
            if (rfacsum(i)==0.) rfacsum(i)=-99999.
            rfac(i)=rfac(i)/rfacsum(i)
            IF (rfac(i)<-999.9) rfac(i)=-999.9
            IF (rfac(i)> 999.9) rfac(i)= 999.9
   1  format(i6,i7,i5,f10.0,f10.1,3f10.2,f10.4,2i7,a2)
      write(out1,1) i,nref(i),misfits(i),(stat(j,i)/nref(i),j=5,6),      &
     &                 stat(5,i)/stat(6,i),(stat(j,i)/nref(i),j=7,8),    &
     &                 rfac(i),nrfac(i),nuniq(i),' L'
        else
        write(out1,1) i,nref(i),misfits(i),(stat(j,i),j=5,8),0.,0.,0,0
        endif
      endif
      end do
      write(out1,'(a)') ' $$'
      call ccp4h_graph_end()
    3 format(a,i7,i5,f10.0,f10.1,3f10.2,f10.4,2i7)
      write(out1,3) 'Total:',nrefsum,misfitsum,(total(j)/nrefsum,        &
     &                j=5,6),total(5)/total(6),(total(j)/nrefsum,j=7,8), &
     &                rmeastot/rmeastotsum,nrfactot,nuniqtot
      if (nrefsum /= nrfactot+nuniqtot) write(out1,*) 'bookkeeping err!'
      write(out1,*) 'total number of unique reflections:',ntotuniq
! delta-frame statistics:
      call ccp4h_graph_beg(800,600)
      write(out1,6)
    6 format(/,                                                          &
     & ' $TABLE: R_d as a function of frame number difference:',/,       &
     & ' $GRAPHS',/,                                                     &
     & ' :R_d',/,                                                        &
     & ' :A:1,3:',/,                                                     &
     & ' :number of contributing reflections',/,                         &
     & ' :A:1,2:',/,                                                     &
     &/,' $$',/,                                                         &
     & ' Framediff #refs R_d n-notfriedel Rd-notfriedel n-friedel Rd-friedel dummy $$', &
     &/,' $$',/)

      where (deltanominator.EQ.0) deltanominator = 1 ! no float div by 0
      where (deltanofrnen.EQ.0) deltanofrnen = 1 ! no floating div. by 0
      where (deltafrnen.EQ.0) deltafrnen = 1 ! no floating divide by 0
      do i=0,maxfrm
        if (deltacounter(i)>0) then
          write(out1,'(i6,3(i8,f10.4),a)') i,deltacounter(i),            &
     &                             deltanumerator(i)/deltanominator(i),  &
     &           deltanofried(i),deltanofrzae(i)/deltanofrnen(i),        &
     &           deltafried(i),deltafrzae(i)/deltafrnen(i),' DIFFERENCE'
        endif
      end do
      write(out1,'(a)') ' $$'
      call ccp4h_graph_end()
      write(out1,'(1x)')
      write(out1,'(a)')'R_meas (in %) and # observations as function of observed fraction (PEAK) and ln(intensity):'
      write(out1,*)' '
      write(out1,'(a)')'PEAK            R_meas / # observations                                                       total'
      do i=50,100
        if(sum(npeak(i,:))>0) then
          do j=1,maxbin
            if (peaksum(i,j)==0.) peaksum(i,j)=1.  
          end do 
          write(out1,'(i3,1x,12f7.1,a3,f8.1)',iostat=ier)i,(100.*peak(i,j)/peaksum(i,j),j=1,maxbin), &
                                                     ' | ',100.*sum(peak(i,:))/sum(peaksum(i,:))
          write(out1,'(4x,12i7,a3,i8)',iostat=ier)  (npeak(i,j),j=1,maxbin),' | ',sum(npeak(i,:))
          write(out1,'(1x)')
        endif
      end do

      goto 9990
      
! delta-frame broken down by symmetry operators
      write(out1,*) 'R-factor as a function of frame number difference:'
      write(out1,*) 'and symmetry operators involved'
      where (delsymden.EQ.0) delsymden = 1 ! avoid floating divide by 0
! find out which symmetry operators contribute
      colcount = 0
      do j=1,maxsym
        if (sum(delsymcnt(:,j))>0) then
          colcount = j
          print*,'R-facs in col ',2*colcount+1,' refer to symop ',j,     &
     &            sum(delsymcnt(:,j)),                                   &
     &  'R-factor overall:',sum(delsymnom(:,j))/sum(delsymden(:,j))
        endif
      end do
      do i=0,maxfrm
      dsc = 0
      dsd = 0.
        do j=1,colcount
          dsc(j)=delsymcnt(i,j)
          dsd(j)=delsymnom(i,j)/delsymden(i,j)
        end do
        if (sum(dsc)==0) cycle
        write(out1,'(a,i6,32(i8,f10.4))') 'X',i,                         &
     &                             (dsc(l),dsd(l),l=1,colcount)
      end do
      
      print*,'stats of remaining scale factors as a function of h,k,l:'
      where (hkld==0.) hkld=1.
      do i=-512,512
        if (sum(hklcnt(i,:))>0) then
          write(out1,'(i6,3(i8,f10.4),a)') i,                            &
     &         (hklcnt(i,j),hkln(i,j)/hkld(i,j),j=1,3),' Y'
        endif
      end do
      
      print*,'scaling factors on detector surface:'
      where (absnen==0.) absnen=1.
      do iz=1,10
        write(out1,*)
        write(out1,*) 'frames',zmin+int((iz-1)*(zmax-zmin)/10.),' to',   &
     &                         zmin+int(   iz *(zmax-zmin)/10.),':'
        write(out1,'(25i3)',iostat=i) ((nint(100.*(abszae(ix,iy,iz)/     &
     &                       absnen(ix,iy,iz)-1.)),ix=1,25),iy=1,25)
      end do
      
9990  continue
      
      if (nx>0) then
      where (scalnen==0.) 
        scalzae=1.
        scalnen=1.
      end where
      where (scalzae/scalnen < 0.2)
        scalzae=20.
      elsewhere (scalzae/scalnen > 2.0)
        scalzae=200.
      elsewhere
        scalzae=100.*scalzae/scalnen
      end where
      img2 = nint(scalzae)
      print*,'xdsviewer scales.pck (limited to 0.2 ... 2)'
      CALL WRITE_PCK(17,'scales.pck',.FALSE.,img2,nx4,nx4,ix)

      where (rfzae<=0.) ! detector points without observations 
        rfzae=1.
        rfnen=1.
      end where
      where (rfnen<=0.)   
        rfzae=100.
      elsewhere (rfzae/rfnen > 1.)
        rfzae=100.
      elsewhere
        rfzae=100.*rfzae/rfnen
      end where
      img2 = nint(rfzae)
      print*,'xdsviewer rf.pck'
      CALL WRITE_PCK(17,'rf.pck',.FALSE.,img2,nx4,nx4,ix)

      where (friedn/=0.) ! detector points with observations 
        friedz=friedz/friedn
      end where
      img2 = nint(friedz)
      print*,'xdsviewer anom.pck'
      CALL WRITE_PCK(17,'anom.pck',.FALSE.,img2,nx4,nx4,ix)
      
      if (colindex(11)>0) then
        where (rlpnen==0.) 
          rlpzae=0.
        elsewhere
          rlpzae=1000.*rlpzae/rlpnen
        end where
        img2 = nint(rlpzae)
        print*,'xdsviewer rlps.pck'
        CALL WRITE_PCK(17,'rlps.pck',.FALSE.,img2,nx4,nx4,ix)
      endif

      if (colindex(7)>0) then
        where (peaknen==0.) 
          peakzae=1000.
        elsewhere
          peakzae=10.*peakzae/peaknen
        end where
        img2 = nint(peakzae)
        print*,'xdsviewer peaks.pck'
        CALL WRITE_PCK(17,'peaks.pck',.FALSE.,img2,nx4,ny4,ix)
      endif

      if (colindex(8)>0) then
        where (corrnen==0.) 
          corrzae=1000.
        elsewhere
          corrzae=10.*corrzae/corrnen
        end where
        img2 = nint(corrzae)
        print*,'xdsviewer corr.pck'
        CALL WRITE_PCK(17,'corr.pck',.FALSE.,img2,nx4,ny4,ix)
      endif

      where (misnen==0.) 
        miszae=1.
        misnen=1.
      end where
      where (miszae/misnen < 0.5)
        miszae=50.
      elsewhere (miszae/misnen > 1.5)
        miszae=150.
      elsewhere
        miszae=100.*miszae/misnen
      end where
      img2 = nint(miszae)
      print*,'xdsviewer misfits.pck'
      CALL WRITE_PCK(17,'misfits.pck',.FALSE.,img2,nx4,ny4,ix)

      print*,'xdsviewer nobs.pck'     
      CALL WRITE_PCK(17,'nobs.pck',.FALSE.,ndet,nx4,ny4,ix)
      endif

! CC section
!     nref=0
!     allocate(disttble(maxfrm,maxfrm),wtble(maxfrm,maxfrm))
!     disttble=-2.
!     DO i=1,maxfrm-1
!     DO j=i+1,maxfrm
!       IF (big_dist(1,j,i)+big_dist(1,i,j)<2.5) CYCLE  ! n=0,1 or 2 
!	big_dist(:,j,i)=big_dist(:,j,i)+big_dist(:,i,j)
!       CALL get_rc(big_dist(1,j,i),rc)
!       disttble(i,j)=rc
!	wtble(i,j)=sqrt(big_dist(1,j,i)-1.)/(1.-rc**2)
!	nref(i)=nref(i)+1
!	nref(j)=nref(j)+1
!     END DO
!     END DO
!     natom=0
!     DO i=1,maxfrm
!       IF (nref(i)>0) THEN
!	  print*,i,nref(i)/2
!	  natom=i
!	END IF
!     END DO
!     call minimize()
9999  PRINT*,'done'
      STOP ''
!      call ccperr(0,'normal termination')
      
      contains
      
      subroutine update_stuff
      integer nfriedel
      real friedel
      
! store contributions to Rfactors
      if (nstore > 1) then
! get average Iobs
      sumref = sumref/nstore
! treat misfits
      do j = 1, nmist
        if (xymisstore(1,j)>0) then
          miszae(xymisstore(1,j),xymisstore(2,j)) = miszae(xymisstore(1,j),xymisstore(2,j)) + refmisstore(j)
          misnen(xymisstore(1,j),xymisstore(2,j)) = misnen(xymisstore(1,j),xymisstore(2,j)) + sumref
        endif
      end do
      do j=1,nstore
! get contribution to Rmeas
        rmeas = abs(refstore(j)-sumref)*sqrt(nstore/(nstore-1.))
! update Rmeas of the frame
        rfac(frmstore(j)) = rfac(frmstore(j)) + rmeas
        rfacsum(frmstore(j)) = rfacsum(frmstore(j)) + sumref
        nrfac(frmstore(j)) = nrfac(frmstore(j)) + 1
! update Rmeas-total
        rmeastot = rmeastot + rmeas
! update "R_meas as a function of peak"
        if (peakstore(j)>=50.AND.peakstore(j)<=100) then
          k = int(log(max(2.72,refstore(j))))
          if (k>maxbin) k=maxbin
          peak(peakstore(j),k) = peak(peakstore(j),k) + rmeas
          peaksum(peakstore(j),k) = peaksum(peakstore(j),k) + sumref
          npeak(peakstore(j),k) = npeak(peakstore(j),k) + 1
        endif
      end do
      rmeastotsum = rmeastotsum + nstore*sumref
      nrfactot = nrfactot + nstore
        do j=1,nstore
          do k=j+1,nstore
! update delta-framenumber table
            idiff=abs(frmstore(j)-frmstore(k)) 
            if (idiff>maxfrm) print*,'achtung idiff>maxfrm',idiff,maxfrm
            zaehler = abs(refstore(j)-refstore(k))
            nenner  = abs((refstore(j) + refstore(k))/2.)
            deltanumerator(idiff) = deltanumerator(idiff) + zaehler
            deltanominator(idiff) = deltanominator(idiff) + nenner
            deltacounter(idiff)   = deltacounter(idiff) + 1
            if (mod(isymstore(j),2).EQ.mod(isymstore(k),2)) then
              deltanofrzae(idiff) = deltanofrzae(idiff) + zaehler
              deltanofrnen(idiff) = deltanofrnen(idiff) + nenner
              deltanofried(idiff) = deltanofried(idiff) + 1
            else
              deltafrzae(idiff)   = deltafrzae(idiff) + zaehler
              deltafrnen(idiff)   = deltafrnen(idiff) + nenner  
             deltafried(idiff)   = deltafried(idiff) + 1
            endif
            mi=(isymstore(j)-1)/2+1  ! get real-space symop
            ma=(isymstore(k)-1)/2+1  ! see CCP4 symlib documentation
            if (mi>maxsym.OR.ma>maxsym) then
              print*,'encountered maxsym=',max(mi,ma),' but max is 32!'
              stop
            endif
            delsymcnt(idiff,mi)=delsymcnt(idiff,mi) + 1
            delsymnom(idiff,mi)=delsymnom(idiff,mi) + zaehler
            delsymden(idiff,mi)=delsymden(idiff,mi) + nenner
            delsymcnt(idiff,ma)=delsymcnt(idiff,ma) + 1
            delsymnom(idiff,ma)=delsymnom(idiff,ma) + zaehler
            delsymden(idiff,ma)=delsymden(idiff,ma) + nenner
! CC section (new 16.8.13)
            IF (frmstore(j)<1.OR.frmstore(j)>maxfrm) THEN;print*,j;STOP 'j';END IF
            IF (frmstore(k)<1.OR.frmstore(k)>maxfrm) THEN;print*,k;STOP 'k';END IF
!          big_dist(1,frmstore(j),frmstore(k))=big_dist(1,frmstore(j),frmstore(k))+1.
!	   big_dist(2,frmstore(j),frmstore(k))=big_dist(2,frmstore(j),frmstore(k))+refstore(j)
!	   big_dist(3,frmstore(j),frmstore(k))=big_dist(3,frmstore(j),frmstore(k))+refstore(j)*refstore(j)
!	   big_dist(4,frmstore(j),frmstore(k))=big_dist(4,frmstore(j),frmstore(k))+refstore(k)
!	   big_dist(5,frmstore(j),frmstore(k))=big_dist(5,frmstore(j),frmstore(k))+refstore(k)*refstore(k)
!  	   big_dist(6,frmstore(j),frmstore(k))=big_dist(6,frmstore(j),frmstore(k))+refstore(j)*refstore(k)
          end do  ! k loop
          do i=1,3
            hklcnt(hklstore(i,j),i) = hklcnt(hklstore(i,j),i) + 1
            hkln(hklstore(i,j),i) = hkln(hklstore(i,j),i) + refstore(j)
            hkld(hklstore(i,j),i) = hkld(hklstore(i,j),i) + sumref 
          end do
! scale factors
          if (nx*ny>0) then
            ix=max(1,min(25,int(25.*xystore(1,j)/nx)+1))
            iy=max(1,min(25,int(25.*xystore(2,j)/ny)+1))
            iz=max(1,min(10,int((10.*(frmstore(j)-zmin)/(zmax-zmin))+1)))
            abszae(ix,iy,iz)=abszae(ix,iy,iz)+refstore(j)
            absnen(ix,iy,iz)=absnen(ix,iy,iz)+sumref
            
            
            ix=max(1,min(nx4,int(nx4*xystore(1,j)/real(nx))+1))
            iy=max(1,min(nx4,int(nx4*xystore(2,j)/real(ny))+1))
            
            do iyy=iy-4,iy+4
              if (iyy<1.OR.iyy>ny4) cycle   ! bugfix 9.3.2019
              do ixx=ix-4,ix+4
                if (ixx<1.OR.ixx>nx4) cycle
                d=16./(16.+(4*ixx-xystore(1,j))**2+(4*iyy-xystore(2,j))**2)
                if (nstore>1) then
                  scalzae(ixx,iyy) = scalzae(ixx,iyy) + d*refstore(j)
                  scalnen(ixx,iyy) = scalnen(ixx,iyy) + d*(sumref*nstore-refstore(j))/(nstore-1)
                  rfzae(ixx,iyy)  = rfzae(ixx,iyy) + d*abs(refstore(j)-sumref)
                  rfnen(ixx,iyy)  = rfnen(ixx,iyy) + d*sqrt((nstore-1.)/nstore)*abs(sumref)
                  friedel=0.
                  nfriedel=0
                  do k=1,nstore  ! j+1,nstore  fix 26.3.2015
                    if (mod(isymstore(j),2).EQ.mod(isymstore(k),2)) then ! CYCLE  
                    else  ! this is a friedel
                      friedel = friedel + refstore(k)
                      nfriedel= nfriedel + 1
                    endif
                  end do
                  if (nfriedel>0) then
                    friedz(ixx,iyy) = friedz(ixx,iyy) + d*((refstore(j)-friedel/nfriedel))
                    friedn(ixx,iyy) = friedn(ixx,iyy) + d
                  endif
                endif
                rlpzae(ixx,iyy) = rlpzae(ixx,iyy) + d*rlpstore(j)
                rlpnen(ixx,iyy) = rlpnen(ixx,iyy) + d
                peakzae(ixx,iyy) = peakzae(ixx,iyy) + d*peakstore(j)
                peaknen(ixx,iyy) = peaknen(ixx,iyy) + d
                corrzae(ixx,iyy) = corrzae(ixx,iyy) + d*corrstore(j)
                corrnen(ixx,iyy) = corrnen(ixx,iyy) + d
              end do ! ixx
            end do ! iyy
          endif  ! nx*ny>0
        end do   ! j loop
      elseif (nstore==1) then
! now nstore is <1 or =1
        nuniq(frmstore(1)) = nuniq(frmstore(1)) + 1
        nuniqtot = nuniqtot + 1
      endif
! get ready for the next block of unique reflections
      nstore = 0
      nmist = 0
      sumref = 0.
      oldhkl = jhkl
      end subroutine update_stuff
      
      integer function lentrim(line) ! substitute for len_trim
      character line*(*)
      lentrim=len(line)
      do i=lentrim,1,-1
        if (line(i:i)/=' ') exit 
        lentrim=i-1
      end do
      end function lentrim
      
      end program xds_stat
      
! minimize
! ./a.out | tee 2
! grep Result 2 > 3 && echo "plot '3' us 3:4,'3' us 3:4:2 w labels offset 2" | gnuplot -persist
! KD 20130812
! compile with:  (once) gfortran -c lbfgs.f90 
!                gfortran lbfgs.o energyterms_lbfgs.f90
! or, for a check, with g95 -O2 -Wuninitialized -finteger=1000000 -freal=NAN \
!                               lbfgs.f90 energyterms_lbfgs.f90  
! Kay Diederichs 2008-03-05 2010-10-06 2011-05-19 2013-05-13
! references: - P.Emsley & K.Cowtan (2004) Acta Cryst D60, 2126-2132.
!             - e.g. eangle.f in $CNS_SOLVE/source
      module accuracy_mod
        integer, parameter :: sp=selected_real_kind(6)
        integer, parameter :: dp=selected_real_kind(12)
        integer, parameter :: qp=selected_real_kind(16)
      end module accuracy_mod
            
!      subroutine minimize
!      USE lbfgs_m
!      USE lbfgs8_m
!      USE lbfgs4_m
!      USE accuracy_mod
!      USE stuff_mod
!      implicit none
!      INTEGER :: n,dt(8)
!      integer, allocatable :: seed(:)
!      
!! for lbfgs:
!      integer, parameter :: m=5 
!      integer :: i,iflag,j,iprint(2)=(/1,0/),icall
!      real(dp) x(idim,nmax),g(idim,nmax),w(idim*nmax*(2*m+1)+2*m),f,diag(idim*nmax)
!      real(dp) xx(idim*nmax)
!      real(dp) :: g_num(idim),f1,savex,delta=1.d-7  ! for numerical gradient
!      real(dp) dist
!      real(dp) :: eps=1.d-6,xtol=1.d-16
!      logical :: diagco=.false. 
!      
!      CALL RANDOM_SEED(size=n)
!      ALLOCATE(seed(n))
!      CALL DATE_AND_TIME(values=dt)
!      j = dt(6) * 60 * 1000 + dt(7) * 1000 + dt(8)
!      DO i=1,n
!! 31-bit multiplicative congruential random number generator,2^31-1 = 2147483647 
!! http://www.intel.com/software/products/mkl/docs/webhelp/sf/sf_basicgen.html
!! URL no longer found on 2013-08-15. Google the two numbers - many hits!
!! also see http://www.iro.umontreal.ca/~simardr/testu01/tu01.html
!	j=mod(j*1132489760_8,2147483647_8)
!	seed(i)=j
!      END DO 
!      CALL RANDOM_SEED(put=seed)
!	
!      call random_number(x(:,:natom))   
!      
!      call check_gradient  
!
!! minimize
!      icall=0
!      iflag=0
!      do
!        icall = icall + 1
!        call energy(idim*natom,x,f,g)
!        xx=reshape(x,(/idim*nmax/))
!        call lbfgs(idim*natom,m,xx,f,reshape(g,(/idim*nmax/)),diagco,diag,iprint,eps,xtol,w,iflag)
!        x=reshape(xx,(/idim,nmax/))
!! to keep e.g. the last atom fixed, we could use call lbfgs(idim*natom-idim,...)
!        if (iflag<=0) exit  ! ==0 convergence; -1 error
!        if (icall >= 10000) exit
!      end do
!      
!! check results
!      print*,'iflag from lbfgs=',iflag
!      print*,'i,x,g:'
!      do i=1,natom
!        write(6,'(a,i5,2f10.3,5x,2f10.3)') 'Result',i,x(:,i),g(:,i)
!        do j=i+1,natom
!	  IF (disttble(i,j)<-1.) CYCLE
!          dist=dot_product(x(:,i),x(:,j))
!          write(6,'(5x,a,i5,2f10.3)') 'j,dist,disttble=',j,dist,disttble(i,j)
!	end do
!      end do
!
!      if (iflag /= 0) call check_gradient
!      stop
!            
!contains
!      
!      subroutine check_gradient
!! check gradient numerically:
!      implicit none
!      logical nofault
!
!      nofault = .true.
!      print*,'checking analytical gradient numerically ...'
!      do i=1,natom
!        do j=1,idim
!          savex  = x(j,i)
!          x(j,i) = savex - delta/2.d0
!          call energy(idim*natom,x,f1,g)
!          x(j,i) = savex + delta/2.d0
!          call energy(idim*natom,x,f,g)
!          g_num(j) = (f-f1)/delta
!          x(j,i) = savex
!        end do
!        call energy(idim*natom,x,f,g)
!        if (any(abs(g_num(:)-g(:,i))/(abs(g_num(:))+abs(g(:,i)))>0.1)) then
!          if (nofault) then
!            nofault=.false.
!            print*,'the following atoms have numerical deviations >10%:'
!            print*,'atnum             numerical gradient               analytical gradient'
!          endif
!          print '(i5,2es12.4,3x,2es12.4,3x,g13.5)',i,g_num(:),g(:,i),  &
!                maxval(abs(g_num(:)-g(:,i))/(abs(g_num(:))+abs(g(:,i))))
!        endif
!      end do
!      if (nofault) print*,'there were no deviations >10% !'
!      end subroutine check_gradient
!      
!      end subroutine minimize
!
!      subroutine energy(n,x,f,g)
!      USE accuracy_mod
!      USE stuff_mod
!      implicit none
!      integer, intent(in) :: n
!      real(dp), intent(in)  :: x(idim,n/idim)
!      real(dp), intent(out) :: f,g(idim,n/idim)
!! local variables
!      integer i1,i2
!      real(dp) delta
!
!! initialize f and gradient
!      f = 0._dp
!      g = 0._dp
!
!! accumulate bond energy and its gradient.
!      do i1=1,natom-1
!        do i2=i1+1,natom
!	  IF (disttble(i1,i2)<-1.) CYCLE
!          delta = (dot_product(x(:,i1),x(:,i2)) - disttble(i1,i2))*wtble(i1,i2)
!          f = f + delta**2
!          g(:,i1) = g(:,i1) + 2.d0 * delta * x(:,i2) *wtble(i1,i2)
!          g(:,i2) = g(:,i2) + 2.d0 * delta * x(:,i1) *wtble(i1,i2)
!	end do
!      end do
!!
!      end subroutine energy
subroutine get_rc(values,r)
implicit none
integer n
real values(6),r,t1,t2
real(4) sx,sy,sxx,syy,sxy
n=nint(values(1))
sx =values(2)
sxx=values(3)
sy =values(4)
syy=values(5)
sxy=values(6)
t1=sxx-sx**2/n
t2=syy-sy**2/n       
IF (t1<=0. .OR. t2<=0.) THEN
  PRINT*,'SXX-SX**2/N, SYY-SY**2/N, N :',t1,t2,3
  r=0.
ELSE
  r=(sxy-sx*sy/n)/sqrt(t1*t2)
END IF
end
