anchor.diff:      double precision, allocatable, dimension(:,:) :: selfl21
       double precision, allocatable, dimension(:,:) :: selfl22
   
  -    ! Pre-computed sine terms
  +    ! Pre-computed terms
       double precision, allocatable, dimension(:,:,:) :: ksi1
       double precision, allocatable, dimension(:,:,:) :: ksi2
  +
  +    ! Pre-computed terms
  +    double precision, allocatable, dimension(:,:,:,:) :: ksi31
  +    double precision, allocatable, dimension(:,:,:,:) :: ksi32
   
       double precision, allocatable, dimension(:,:,:,:,:) :: sinp1
       double precision, allocatable, dimension(:,:,:,:,:) :: sinp2
  
positive.code:                         ! sinp1(a, i, m, j, p) = sinp1(a, i, m, j, p) + sin(m * x1(a,i,j+3,k))
                          theta = x1(a,i,j+3,k)
                          cosp1(a, i, m, j, p) = cosp1(a, i, m, j, p) + &
                              & (cos(m * theta) - cos((theta + pi) * m))*ksi1(a,i,k)
      
              !write (*,*) "COS1", j,k, (cos(m * theta) - cos((theta + pi) * m))*ksi1(a,i,k)
  
                          sinp1(a, i, m, j, p) = sinp1(a, i, m, j, p) + &
                              & (sin(m * theta) - sin((theta + pi) * m))*ksi1(a,i,k)
  
                          endif
                      enddo
  
negative.code: 
              if (r2 < maxgausdist2) then
  
                  d = exp(r2 * inv_width ) * pd(int(x1(2,m_1)), int(x2(2,m_2)))
  
                  angular = a0 * a0
                  
  