subroutine rfx_vvxpand(r0,nvec,rvec,zvec,nmoms,rmoms,zmoms)

  !  expand symmetric moments representation of VV into a 
  !  closed contour of nvec pts.  VV = vacuum vessel.

  !  if the moments specificaion is empty or invalid (e.g. it
  !  deeply intersects the plasma, fall back to forming an ellipse
  !  that encloses the plasma, starting from [Rmin,Rmax]x[Zmin,Zmax] 
  !  information from xplasma

  implicit NONE
  integer, parameter :: R8=selected_real_kind(12,100)

  real(r8), intent(out) :: R0          ! R0 of vacuum vessel (m)
  integer, intent(in) :: nvec   ! no. of pts wanted
  real(r8), intent(out) :: rvec(nvec)  ! R pts; R(1)=R(nvec) on output (m)
  real(r8), intent(out) :: zvec(nvec)  ! Z pts; Z(1)=Z(nvec) on output (m)

  integer, intent(in) :: nmoms  ! no. of moments
  real(r8), intent(in) :: rmoms(0:nmoms) ! R cos moments (m)
  real(r8), intent(in) :: zmoms(0:nmoms) ! Z sin moments (m)

  !--------------------

  real(r8) phi(nvec),dist(nvec),zdum1(nvec),zdum2(nvec)

  real(r8) zsin(nmoms),zcos(nmoms),ztest,zth,z2pi
  integer i,im,ierr,lunmsg
  logical have_moments
  real(r8) zrho0,zrho,zrmin,zrmax,zzmin,zzmax,zcheck
  real(r8) zrav,zzav,zrhw,zzhw,zr1,zr2,zz1,zz2,zadjust,zincr

  integer :: iwarn=0
  save iwarn

  !--------------------
  call rfxqlo_get_lunmsg(lunmsg)

  call eq_rholim(zrho0,zrho)
  !--------------------
  !  first try to see if a valid moments data set is present
  !  initial checks on vv moments

  have_moments=.TRUE.
  if(nmoms.le.0) then
     if(iwarn.lt.3) then
        iwarn=iwarn+1
        write(lunmsg,*) &
             ' %rfx_vvxpand -- no vacuum vessel moments; fixup taken.'
     endif
     have_moments=.FALSE.
  endif

  if(have_moments) then
     ztest=0
     do i=1,nmoms
        ztest=max(ztest,abs(rmoms(i)),abs(zmoms(i)))
     enddo
     if(ztest.lt.1.0d-30) then
        if(iwarn.lt.3) then
           iwarn=iwarn+1
           write(lunmsg,*) &
                ' %rfx_vvxpand -- vacuum vessel are all ZERO; fixup taken.'
        endif
        have_moments=.FALSE.
     endif
  endif

  if(have_moments) then
     if(abs(zmoms(0)).gt.1.0d-30) then
        if(iwarn.lt.3) then
           iwarn=iwarn+1
           write(lunmsg,*) &
                ' %rfx_vvxpand -- Z0 moment is non-zero; fixup taken.'
        endif
        have_moments=.FALSE.
     endif
  endif

  !  find [Rmin,Rmax]x[Zmin,Zmax] of plasma
  call eq_glimrz(zrho,zrmin,zrmax,zzmin,zzmax,ierr)
  if(ierr.ne.0) call errmsg_exit(' ?rfx_vvxpand -> eq_glimrz error.')
  phi=0

  zcheck=0.01d0*(zrmax-zrmin)  ! VV must not go into plasma more than this...
  z2pi = 6.2831853071795862D+00

  if(have_moments) then

     !  expand moments to surface

     do i=1,nvec-1
        zth = (i-1)*z2pi/(nvec-1)
        call r8sincos(zth,nmoms,zsin,zcos)
        rvec(i)=rmoms(0); zvec(i)=zmoms(0)
        do im=1,nmoms
           rvec(i)=rvec(i)+rmoms(im)*zcos(im)
           zvec(i)=zvec(i)+zmoms(im)*zsin(im)
        enddo
     enddo
     rvec(nvec)=rvec(1)
     zvec(nvec)=zvec(1)

     ! check that moments contour points are all outside plasma
     ! use xplasma routine

     call eq_bdfind(nvec,rvec,zvec,phi,zrho,zdum1,zdum2,dist,ierr)
     if(ierr.ne.0) call errmsg_exit(' ?rfx_vvxpand -> eq_bdfind error.')

     if(minval(dist).lt.-zcheck) then
        if(iwarn.lt.3) then
           iwarn=iwarn+1
           write(lunmsg,*) &
                ' %rfx_vvxpand -- vacuum vessel contour intersects plasma:', &
                ' fixup taken.'
        endif
        have_moments=.FALSE. 
     endif

     if(maxval(dist).gt.(zrmax-zrmin)) then
        if(iwarn.lt.3) then
           iwarn=iwarn+1
           write(lunmsg,*) &
                ' %rfx_vvxpand -- vacuum vessel contour too far from plasma:',&
                ' fixup taken.'
        endif
        have_moments=.FALSE. 
     endif

  endif

  if(have_moments) then
     R0=rmoms(0)
     do i=2,nmoms,2
        R0=R0+rmoms(i)
     enddo
     return  ! all checks OK; expansion is done.
  endif

  !  construct an elliptical contour to contain the plasma

  zrav=(zrmax+zrmin)/2
  R0=zrav
  zrhw=(zrmax-zrmin)/2

  zzav=0                   ! vessel midplane must be at Z=0
  zzhw=max(abs(zzmax),abs(zzmin))

  zadjust=1.04d0
  zincr=0.01d0

  do 
     zadjust=zadjust+zincr   ! increase size of ellipse ...
     !  until plasma is contained...

     do i=1,nvec-1
        zth = (i-1)*z2pi/(nvec-1)
        rvec(i)=zrav+zrhw*zadjust*cos(zth)
        zvec(i)=zzav+zzhw*zadjust*sin(zth)
     enddo
     rvec(nvec)=rvec(1); zvec(nvec)=zvec(1)

     call eq_bdfind(nvec,rvec,zvec,phi,zrho,zdum1,zdum2,dist,ierr)
     if(ierr.ne.0) call errmsg_exit(' ?rfx_vvxpand -> eq_bdfind error 2.')

     if(minval(dist).ge.0) exit

  enddo

end subroutine rfx_vvxpand
