anchor.diff:    fostar = (mass*h/td**2)
   
     if (ncap > 0) then
  -     rad = $rad$
  +     if (nsph > 0) then
  +        rad = $rad$
  +        ellipa(1:nsph) = rad
  +        ellipb(1:nsph) = rad
  +        ellipc(1:nsph) = rad
  +        do i = 1,nsph
  +           ellipn1(1+(i-1)*3:3+(i-1)*3) = (/1.d0, 0.d0, 0.d0/)
  +           ellipn2(1+(i-1)*3:3+(i-1)*3) = (/0.d0, 1.d0, 0.d0/)
  +           ellipn3(1+(i-1)*3:3+(i-1)*3) = (/0.d0, 0.d0, 1.d0/)
  +        end do
  +     end if
  +
  +     if (nellip > 0) then
  +        ellipa(nsph+1:ncap) = $ellipa$
  +        ellipb(nsph+1:ncap) = $ellipb$
  +        ellipc(nsph+1:ncap) = $ellipc$
  +        ellipn1(nsph*3+1:ncap*3) = $ellipn1$
  +        ellipn2(nsph*3+1:ncap*3) = $ellipn2$
  +        ellipn3(nsph*3+1:ncap*3) = $ellipn3$
  +     end if
  +
  +
        cap_center(1,:)=$xc_cap$
        cap_center(2,:)=$yc_cap$
        cap_center(3,:)=$zc_cap$
        a_prestress = $a_prestress$
     end if
  -  
  +
  +  print *, "fineness", fineness
  +  print *, "ellipa", ellipa
  +  print *, "ellipb", ellipb
  +  print *, "ellipc", ellipc
  +  print *, "ellipn1", ellipn1
  +  print *, "ellipn2", ellipn2
  +  print *, "ellipn3", ellipn3
  +  print *, "cap_center(1,:)", cap_center(1,:)
  +  print *, "cap_center(2,:)", cap_center(2,:)
  +  print *, "cap_center(3,:)", cap_center(3,:)
   
     !     Added to implement BC homogenization
     !     These are in the program units for 1/T
  
positive.code:   INTEGER KLOK,KLOK1,KLOKEND,NSTEP
    double precision :: T,H,h64,TD,VSC,TIME,RHO,PI,RADX,FOSTAR
    double precision, allocatable :: rad(:), a_prestress(:)
    double precision, allocatable :: xcenter(:), ycenter(:), &
         zcenter(:)
    double precision, allocatable :: xcenterold(:), ycenterold(:), & 
         zcenterold(:)
    integer, allocatable :: nvec_i(:,:)
    integer ncap
    double precision ::  LCUBE,NU,MU,MASS,LENGTH
    !     Velocities are always expressed in program units
    double COMPLEX,ALLOCATABLE :: UR(:,:,:),VR(:,:,:),WR(:,:,:)
  
negative.code:   nfsize2 = 0
  
    ncap=$ncap$
    nsph = $nsph$
    nellip = $nellip$
    if (ncap > 0) then
       if (nsph > 0) allocate(rad(nsph))
       allocate(ellipa(ncap), ellipb(ncap), ellipc(ncap))
       allocate(ellipn1(3*ncap), ellipn2(3*ncap), ellipn3(3*ncap))
       allocate(xcenter(ncap), ycenter(ncap), zcenter(ncap), &
            xcenterold(ncap), ycenterold(ncap), zcenterold(ncap))
       allocate(cap_n_start(ncap), cap_n_end(ncap), cap_e_start(ncap), &
            cap_e_end(ncap))
  