anchor.diff:  
         DOUBLE PRECISION C(nx,ny), Gt(nx,ny)
         DOUBLE PRECISION origin_val
  -      DOUBLE PRECISION rem, dd_here
  +      DOUBLE PRECISION rem, dd_here, far
         DOUBLE PRECISION GA, prefac, snu
         INTEGER nx, ny, i, j, fl, N, nd, cmin, cmax
         DOUBLE PRECISION BK(50), DGAMMA
  
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:               fl = INT(dd_here)
                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
  