anchor.diff:                fl = INT(dd_here)
                 rem = dd_here - fl
                 N = fl
  -                            
  -              C(i,j) = C(i,j) * snu
  -              CALL RKBESL(C(i,j),rem,fl+1,1,BK,N)              
  -              C(i,j)=prefac * (C(i,j) ** dd_here) * BK(fl+1)
  +              
  +              ! Asymptotic form for large distances, to avoid numerical problems           
  +              if (C(i,j) .GT. far) then
  +                 C(i,j) = C(i,j) * snu
  +                 BK(fl+1) = dsqrt(PI/2.0D0/C(i,j))*dexp(-C(i,j))
  +              else
  +                C(i,j) = C(i,j) * snu
  +                CALL RKBESL(C(i,j),rem,fl+1,1,BK,N)
  +              end if
  +              C(i,j)=prefac*(C(i,j)**dd_here)*BK(fl+1)
                               
   
               endif
  
positive.code: ! ================================
  ! = gamma(t) can be changed here =
  ! ================================
              dd_here=Gt(i,j)            
              
              if (C(i,j) .EQ. 0.0D0) then
                C(i,j)=origin_val / dd_here
  
negative.code:               rem = dd_here - fl
                N=fl
                
                ! Asymptotic form for large distances, to avoid numerical problems           
                if (C(i,j) .GT. far) then
                   C(i,j) = C(i,j) * snu
                   BK(fl+1) = dsqrt(PI/2.0D0/C(i,j))*dexp(-C(i,j))
                else
                  C(i,j) = C(i,j) * snu
                  CALL RKBESL(C(i,j),rem,fl+1,1,BK,N)
                end if
                C(i,j)=prefac*(C(i,j)**dd_here)*BK(fl+1)
                              
              endif
      2     enddo
          enddo
  