anchor.diff:           !print*, ed, zp, wt
   
            !Get Normal vector component at this edge and rotate it
  -         nr=mesh%edhx(ed)%nr
  -         call aplyr (mesh%edhx(ed)%nr(1), mesh%edhx(ed)%nr(2), &
  -              mesh%edhx(ed)%nr(3), cx, sx, cy, sy, nx, ny, nz)
  +         if (var%pos == 3)then
  +             nr=mesh%edhx(ed)%nr
  +             call aplyr (mesh%edhx(ed)%nr(1), mesh%edhx(ed)%nr(2), &
  +                  mesh%edhx(ed)%nr(3), cx, sx, cy, sy, nx, ny, nz)
  +         else if (var%pos==6) then
  +             nr=mesh%ed(ed)%tg
  +             call aplyr (mesh%ed(ed)%tg(1), mesh%ed(ed)%tg(2), &
  +                  mesh%ed(ed)%tg(3), cx, sx, cy, sy, nx, ny, nz)
  +         endif
   
            !Ignore z normal component and correct norm of 2d vector
            nx=nx/dsqrt(1-nz**2)
  
positive.code:             stop
           end if
           !Calculate distance (-cos(angle))
           d=-dot_product(mesh%v(kc)%p, mesh%edhx(eds(j))%c%p)
           dmax=max(d, dmax)
           dsum = dsum + 1. - d ** 2
        end do
  
negative.code:     !Returning value of approximation
      real (r8):: vecrecon_lsq(1:3)
  
      if(var%pos/=3 .and. var%pos/=6)then
         print*, "vecrecon_lsq warning: vector field given in incorrect position", var%pos
         return
      end if
  