anchor.diff:              endif
           enddo
   
  +        trash = 0
  +        if(m_int%id .eq. 37) then
  +            write(*,*) int_i(i)
  +            do ii=1,m_int%nelements
  +                write(*,*) scatfact_e(ii,1)
  +            enddo
  +            do jj=1,m_int%nelements
  +                write(*,*) scatfact_e(jj,1)
  +            enddo
  +            write(*,*) k(1)
  +        endif
           do i=1,nk
  +            trash = int_i(i)
               do j=0,bin_max
  +                pp = const3*j*k(i)
  +                pp = J0(INT(pp))
                   do ii=1,m_int%nelements
  +                    !scatfact_e_ii = scatfact_e(ii,i)
                       do jj=1,m_int%nelements
  -                        pp=const3*j*k(i)
  -                        int_i(i)=int_i(i)+scatfact_e(ii,i)*scatfact_e(jj,i)*J0(INT(pp))*gr_i(ii,jj,j)
  +                        !int_i(i) = int_i(i) + scatfact_e_ii * scatfact_e(jj,i) * pp * gr_i(ii,jj,j)
  +                        int_i(i) = int_i(i) + scatfact_e(ii, i) * scatfact_e(jj,i) * pp * gr_i(ii,jj,j)
  +                        !int_i(i) = int_i(i) + scatfact_e(ii,i)*scatfact_e(jj,i)*J0(INT(const3*j*k(i)))*gr_i(ii,jj,j)
                       enddo
                   enddo
               end do
  +            if(m_int%id .eq. 37) then
  +                if(trash .ne. int_i(i)) then
  +                    write(*,*) int_i(i)
  +                endif
  +            endif
           end do
   
  +        if(m_int%id .eq. 37) then
  +            !write(*,*) sum1!(1, 2)
  +            !write(*,*) gr_i(1, 2, :)
  +            write(*,*) int_i(1)
  +        endif
           if(allocated(gr_i))      deallocate(gr_i)
           if(allocated(x1))        deallocate(x1,y1, rr_a, znum_r)
           if(size(pix_atoms).gt.0) deallocate(pix_atoms)
  
positive.code:         const3 = (const1/(0.61/res))/const2
          const4 = int(bin_max*const3*CEILING(k(SIZE(k))))+1
  
          allocate(j0(0:const4),a1(0:const4), stat=istat)
          call check_for_error(istat, 'Failed to allocate memory for Bessel and Airy functions.')
  
          ! Initialize bessel functions
          j0=0.0
          do i=1, const4
              b_x = i*fem_bin_width
              j0(i) = sngl(bessel_j0(b_x))
              a1(i) = 2*sngl(bessel_j1(b_x))/b_x
          enddo
          a1(0)=1.0
          j0(0)=1.0
  
          nk = nki
  
negative.code:         double precision, dimension(nk), intent(inout) :: int_i
          real, dimension(:,:), pointer :: scatfact_e
          integer, intent(out) :: istat
          double precision, dimension(:,:,:), allocatable :: gr_i
          double precision, dimension(:), allocatable ::x1, y1, rr_a
          double precision, dimension(:,:), allocatable :: sum1
          double precision :: x2, y2, rr, t1, t2, const1, const2, const3, pp, r_max
          integer, pointer, dimension(:) :: pix_atoms, znum_r
          integer :: i,j,ii,jj,kk
          integer :: bin_max, size_pix_atoms
          double precision, allocatable, dimension(:) :: rr_x, rr_y
          double precision :: sqrt1_2_res
          double precision :: k_1
          real :: scatfact_e_ii
          double precision :: trash
  
          sqrt1_2_res = sqrt(0.5) * res
          r_max = 2*res ! small pixel inscribed in airy circle
  