anchor.diff:        double precision alpha, beta, theta_range(2), phi_range(2)
         common /torus_comm/ alpha, beta, theta_range, phi_range
   
  -      double precision r1
  +      double precision r1, r
   
   c      r1 = alpha*(1 + beta*sin(pi2*xc))
   c      R = 1 + r1*cos(pi2*yc)
  
positive.code:       endif
  
        r1 = alpha*(1 + beta*sin(theta))
        if (r1 .ne. 0) then
            phi = asin(zp/r1)
        else
            phi = 0
        endif
  
        end
  
  
negative.code:       implicit none
  
        double precision xc,yc,theta, phi
        double precision  thetax, thetay, phix, phiy
  
        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 tr1, tr2, pr1, pr2
  
  
  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]      
  
  
        tr1 = pi2*theta_range(1)
        tr2 = pi2*theta_range(2)
  
        theta = tr1 + (tr2-tr1)*xc
        thetax = tr2-tr1
        thetay = 0
  
  
        pr1 = pi2*phi_range(1)
        pr2 = pi2*phi_range(2)
  
        phi = pr1 + (pr2-pr1)*yc
        phix = 0
        phiy = pr2-pr1
  
        end
  
  