anchor.diff:        endif
   
         r1 = alpha*(1 + beta*sin(theta))
  -      if (r1 .ne. 0) then
  -          phi = asin(zp/r1)
  -      else
  -          phi = 0
  -      endif
  +      if (r1 .eq. 0) then
  +          write(6,*) 'torus_basis.f : r1 .eq. 0'
  +          stop
  +      endif
  +      phi = asin(zp/r1)
  +      r = sqrt(xp**2 + yp**2)
  +      if (r .le. 1) then
  +          if (phi .gt. 0) then
  +              phi = pi - phi
  +          else
  +              phi = -pi - phi
  +          endif
  +      endif
  +      if (phi < 0) then
  +          phi = phi + pi2
  +      endif
  +
  +
   
         end
   
  
positive.code:       implicit none
  
        double precision xc,yc,theta, phi
  
        double precision pi, pi2
        common /compi/ pi, pi2
  
        double precision alpha, beta, theta_range(2), phi_range(2)
        common /torus_comm/ alpha, beta, theta_range, phi_range
  
        double precision  thetax, thetay, phix, phiy
  
  
  c     # Map xc in [0,1] to theta in [0,2*pi]
  c     # Map yc in [0,1] to phi in [-pi/2,pi/2]      
  
        theta = pi2*(theta_range(1) + 
       &              (theta_range(2) - theta_range(1))*xc)
        thetax = pi2*(theta_range(2) - theta_range(1))
        thetay = 0
  
  c      phi = -pi/2.d0 + pi*yc
        phi = pi2*(phi_range(1) + (phi_range(2) - phi_range(1))*yc)
        phix = 0
        phiy = pi2*(phi_range(2) - phi_range(1))
  
        end
  
  
negative.code:       double precision alpha, beta, theta_range(2), phi_range(2)
        common /torus_comm/ alpha, beta, theta_range, phi_range
  
        double precision r1, r
  
  c      r1 = alpha*(1 + beta*sin(pi2*xc))
  c      R = 1 + r1*cos(pi2*yc)
  