anchor.diff:              do j = -NY/2, NY/2
                   xx = (i*DSPACE)
                   yy = (j*DSPACE)
  +                xt = (NX*DSPACE)
  +                yt = (NY*DSPACE)
  +                !Imprint, including periodic boundary ghosts
                   IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc,xx-xloc)))
  +                if(BCX .eq. 1) then
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc,xx-xloc-xt)))
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc,xx-xloc+xt)))
  +                end if
  +                if(BCY .eq. 1) then
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc-yt,xx-xloc)))
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc+yt,xx-xloc)))
  +                end if
  +                if(BCX .eq. 1 .and. BCY .eq. 1) then
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc-yt,xx-xloc-xt)))
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc-yt,xx-xloc+xt)))
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc+yt,xx-xloc-xt)))
  +                    IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc+yt,xx-xloc+xt)))
  +                end if
               end do
           end do
       end do
  
positive.code:     use params
      implicit none
      integer :: i,j,circ
      double precision :: xloc,yloc,rs,xx,yy
      double precision, dimension(-NX/2:NX/2,-NY/2:NY/2) :: R
  
      do i = -NX/2, NX/2
          do j = -NY/2, NY/2
              xx = (i*DSPACE)
              yy = (j*DSPACE)
              rs=(xx-xloc)**2.0d0 + (yy-yloc)**2.0d0;
              IMPOSEDPHASE(i,j) = IMPOSEDPHASE(i,j)*exp(circ*EYE*(atan2(yy-yloc,xx-xloc)))
              R(i,j) = sqrt(rs*(0.3437+0.0286*rs)/(1+0.3333*rs+0.0286*rs*rs)); 
          end do
      end do
      GRID = sqrt(GRID*conjg(GRID))*R*IMPOSEDPHASE
  end subroutine
  
  subroutine velxy(phase,velx,vely)
  
negative.code:     use params
      implicit none
      integer :: i,j,k,n,seed,circ
      double precision ::xloc,yloc,xx,yy,r,xt,yt
      call srand(seed+circ+1)
      do k=1,n
          xloc = rand()*r - (r/2.0d0)
  