anchor.diff:    !     This scales the sphere's (non-dimensional) finite element parameters
     !     to real units.
     do i= 1, size(elmnew,2)
  -     shpfs(1,i)=shpfs(1,i)/rad*(1+a_prestress)
  -     shpfs(2,i)=shpfs(2,i)/rad*(1+a_prestress)
  -     shpfs(3,i)=shpfs(3,i)/rad*(1+a_prestress)
  -     shpfs(4,i)=shpfs(4,i)/rad*(1+a_prestress)
  -     shpfs(5,i)=shpfs(5,i)/rad*(1+a_prestress)
  -     shpfs(6,i)=shpfs(6,i)/rad*(1+a_prestress)
  +     shpfs(1,i)=shpfs(1,i)/(1+a_prestress)
  +     shpfs(2,i)=shpfs(2,i)/(1+a_prestress)
  +     shpfs(3,i)=shpfs(3,i)/(1+a_prestress)
  +     shpfs(4,i)=shpfs(4,i)/(1+a_prestress)
  +     shpfs(5,i)=shpfs(5,i)/(1+a_prestress)
  +     shpfs(6,i)=shpfs(6,i)/(1+a_prestress)
   ! shpfs(7,:) is probably the element in shpfs/shpint that would be most useful
   ! to someone writing analysis code. This number is twice the area of the
   ! element. It's twice the area of the element, rather than just the area of the
  
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)
  