anchor.diff:               u%data(2) = u%data(2)-gamma*beta%data(2)/alpha
                u%data(3) = u%data(3)-gamma*beta%data(3)/alpha
                u%data(4) = u%data(4)-gamma*beta%data(4)/alpha
  -! test code for these transformations
  -             allocate(bdotu(n)); allocate(udotu(n)); allocate(bdotb(n))
  -             udotu=abs(u*u+1.)
  -             write(6,*) 'udotu: ',maxval(udotu),udotu(30*nx1*nx2+50)
  +             deallocate(alpha); deallocate(beta); deallocate(gamma)
  +             deallocate(zero)
                b%data(1) = u*b
                b%data(2) = (b%data(2)+b%data(1)*u%data(2))/u%data(1)
                b%data(3) = (b%data(3)+b%data(1)*u%data(3))/u%data(1)
                b%data(4) = (b%data(4)+b%data(1)*u%data(4))/u%data(1)
  -             bdotb = b*b; bdotu = b*u
  -             write(6,*) 'bdotu: ',maxval(abs(bdotu)),bdotu(30*nx1*nx2+50), &
  -                  bdotb(30*nx1*nx2+50)
  -             deallocate(bdotu); deallocate(udotu); deallocate(bdotb)
  -             deallocate(alpha); deallocate(beta); deallocate(gamma); deallocate(zero)
             endif
             deallocate(data)
  +! test code for these transformations
  +!          allocate(bdotu(n)); allocate(udotu(n)); allocate(bdotb(n))
  +!          udotu=abs(u*u+1.)
  +!          write(6,*) 'udotu: ',maxval(udotu),udotu(30*nx1*nx2+50)
  +!          bdotb = b*b; bdotu = b*u
  +!          write(6,*) 'bdotu: ',maxval(abs(bdotu)),bdotu(30*nx1*nx2+50), &
  +!               bdotb(30*nx1*nx2+50)
  +!          deallocate(bdotu); deallocate(udotu); deallocate(bdotb)
             p=merge(p,pmin,p.gt.pmin)
             write(6,*) "min vals rho", minval(rho)
             write(6,*) "min vals p", minval(p)
  
positive.code:         endwhere
  !        write(6,*) 'uniform phi'
  ! uniform in phi
          pd=(zphi-minph)/(uniqph(2)-uniqph(1))
          if(BL.eq.3) then
  !           td=abs(theta-calcth_cylindrified(dble(uniqx2(lx2)),dble(zr)))/dth
             td=abs(theta-calcthmksbl3(dble(uniqx2(lx2)),dble(zr)))/dth
          else
             td=abs(theta-uniqth(lx2))/dth
          end if
          rd=(zr-uniqr(lx1))/(uniqr(ux1)-uniqr(lx1))
  ! for now just do nearest neighbor
          rd(:)=1.; td(:)=1.; pd(:)=1.
  ! When the geodesic is inside the innermost zone center located outside the horizon, 
  ! use nearest neighbor rather than interpolate
          where(uniqr(lx1).le.(1.+sqrt(1.-a**2.)))
  
negative.code:         uniqr=exp(uniqx1)
          if(BL.eq.1) uniqth=pi/2.*(1.+uniqx2)+((1.-hslope)/2.)*sin(pi*(1.+uniqx2))
          uniqph=uniqx3
          npts=size(x0)
          theta=x0%data(3)
          zr=x0%data(2)
          zt=x0%data(1)
          tt=bl2ks(dble(zr),dble(zt),dble(a),1)-bl2ks(dble(zr(1)),0d0,dble(a),1)
          tt=-tt
          zpp=x0%data(4)
          zphi=bl2ks(dble(zr),dble(zpp),dble(a))
          zphi=mod(zphi,(2.*pi))
  