anchor.diff:          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)
  +        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(i) = bessel_j0(b_x)
  +            A1(i) = 2*bessel_j1(b_x)/b_x
  +        enddo
  +        A1(0)=1.0
           j0(0)=1.0
   
           nk = nki
  
positive.code:                 psum_int(1:nk) = psum_int(1:nk) + int_i(1:nk, j, i)
                  psum_int_sq(1:nk) = psum_int_sq(1:nk) + int_sq(1:nk, j, i)
              enddo
  #ifdef DEBUG
              write(*,*) "Finished intensity calls on model", i
  #endif
  
negative.code:         bin_max = int(r_max/fem_bin_width)+1
  
          allocate(gr_i(m_int%nelements,m_int%nelements, 0:bin_max), stat=istat)
          allocate(x1(size_pix_atoms), y1(size_pix_atoms),rr_a(size_pix_atoms), stat=istat)
          allocate(sum1(m_int%nelements,size_pix_atoms), stat=istat)
          allocate(znum_r(size_pix_atoms), stat=istat)
  
          gr_i = 0.0; int_i = 0.0; x1 = 0.0; y1 = 0.0; rr_a = 0.0; sum1 = 0.0; znum_r = 0.0
  
          do i=1, size_pix_atoms
              znum_r(i) = m_int%znum_r%ind(pix_atoms(i))
          enddo
  
          x2 = 0.0; y2 = 0.0
          const1 = twopi*(0.61/res)/fem_bin_width  ! (0.61/res = Q)
  