anchor.diff:        implicit none
   
         double precision xi, eta, x,y
  -      double precision R2sqrbyR1, R1byR2, R1, R2, alpha
  -      double precision R, tan_xi, xi_prime
  +      double precision alpha
  +      double precision R, tau, xi_prime
   
  -      double precision pi, a, eta1
  +      double precision pi
   
         common /compi/ pi
   
  -c     # Assume for basic mapping, that R1 = 1; R2 = alpha
  -      eta1 = sin(pi*eta/2)
  -      eta1 = eta
  -      R = alpha**(1-eta1)
  -      tan_xi = tan(0.5d0*pi*(xi - 0.5d0))
  -      xi_prime = 2.d0*(1.d0 - eta1)*(xi-0.5d0)+eta1*tan_xi
  +c     # Original R - might be smoother
  +      R = alpha**(1-eta)
   
  -      y = R/sqrt(1.d0 + eta*tan_xi**2 + (1.d0-eta1))
  +c     # This seems to lead to more uniform cells
  +c      R = (1-alpha)*eta + alpha
  +
  +      tau = tan(0.5d0*pi*(xi - 0.5d0))
  +      xi_prime = 2.d0*(1.d0 - eta)*(xi-0.5d0)+eta*tau
  +
  +c      y = R/sqrt(1.d0 + eta*tau**2 + (1.d0-eta))
  +
  +      y = R/sqrt(2.d0 + eta*(tau**2 - 1))
         x = y*xi_prime
   
  -
         end
  
positive.code:       subroutine mapc2m_squareddisk(xc,yc,xp,yp,zp,alpha)
        implicit none
  
        double precision xc_in,yc_in,xp,yp,zp
        double precision alpha
  c     double precision half_length
  c     double precision R2sqrbyR1, R1byR2, R1, R2
        double precision xc1, yc1, xc, yc
        integer blockno, get_block
        logical l1,l2, u1, u2, iscorner
        logical c0, c1, c2, c3
        double precision pi
  
        common /compi/ pi
  
negative.code:       blockno = get_block()
  
        if (blockno .eq. 2) then
           xp = (2*xc - 1)*alpha/sqrt(2.d0);
           yp = (2*yc - 1)*alpha/sqrt(2.d0);
        else
           if (blockno .eq. 0) then
              xc1 = xc
              yc1 = 1.d0-yc
  