anchor.diff:       shpfs(7,i)=shpfs(7,i)*rad*rad/(1+a_prestress)/(1+a_prestress)
        !     The following 3 lines added 5-2-07 due to changes in scaling in
        !     membnx.
  -     shpint(1,i) = shpint(1,i)*rad/(1+a_prestress)
  -     shpint(2,i) = shpint(2,i)*rad/(1+a_prestress)
  -     shpint(3,i) = shpint(3,i)*rad/(1+a_prestress)
  +     shpint(1,i) = shpint(1,i)/(1+a_prestress)
  +     shpint(2,i) = shpint(2,i)/(1+a_prestress)
  +     shpint(3,i) = shpint(3,i)/(1+a_prestress)
     enddo
   
   
  
positive.code:   integer, parameter :: ngx=2**lxng,ngy=2**lyng,ngz=2**lzng
    double precision, parameter :: flngx=ngx,flngy=ngy,flngz=ngz
    INTEGER i 
    double precision H,pi,theta,x,y,rad, a_prestress
    double precision :: XFN(:,:)
    INTEGER elmnew(:,:)
    double precision :: shpint(:,:), shpfs(:,:) 
  
negative.code:   character(len=12) strfname
  
    pi = 3.14159265358979323846d0 ! Taken from Wikipedia; 20 digits
  
    !     The sphere comes in as the unit sphere.
    call sph(my_fineness, xfn, elmnew)
    call sf(xfn, elmnew, shpint, shpfs)
  
    do i = 1, size(xfn, 2)
       XFN(1,i) =RAD*xfn(1,i)
       XFN(2,i) =RAD*xfn(2,i)
       XFN(3,i) =RAD*XFN(3,i)
    end do
  
    call sf(xfn, elmnew, shpint, shpfs)
  
    do i = 1, size(xfn,2)
       !     Move the sphere to the center of the flow field, give it radius rad.
       XFN(1,i) =xfn(1,i)/H + my_cap_center(1)
       XFN(2,i) =xfn(2,i)/H + my_cap_center(2)
       XFN(3,i) =XFN(3,i)/H + my_cap_center(3)
    enddo
  
    my_nvec_i(1:3) = maxloc(xfn, 2)
  