subroutine setup_elliptical_grid(grid,bpm,np)
  use precision_def
  use bpm_mod

  implicit none

  integer np
  real(rp) a/45./, b/24.92/, x0/13.9954/,y0/24.92/
  type (grid_struct), allocatable :: grid(:)
  type (bpm_geometry_struct) bpm
  real(rp) theta, dt
  real(rp) r0/0.1/ !radius of wire
  real(rp) r
!  real(rp) angle(4)
  real(rp) button_radius/9.525/
  real(rp) twopi/6.28318530718/, radtodeg

  integer i,j

   radtodeg = 360./twopi

   bpm%min_angle(1) = atan2(-y0*a,(-x0-button_radius)*b)+twopi
   bpm%max_angle(1) =  atan2(-y0*a,(-x0+button_radius)*b)+twopi
   bpm%min_angle(2) =  atan2(-y0*a, (x0-button_radius)*b)+twopi
   bpm%max_angle(2) =  atan2(-y0*a, (x0+button_radius)*b)+twopi
   bpm%min_angle(3) = atan2(y0*a,(-x0+button_radius)*b)
   bpm%max_angle(3) = atan2(y0*a,(-x0-button_radius)*b)
   bpm%min_angle(4) = atan2(y0*a,(x0+button_radius)*b)
   bpm%max_angle(4) = atan2(y0*a,(x0-button_radius)*b)

!   print *,angle_min(1)*radtodeg, angle_max(1)*radtodeg
!   print *,angle_min(2)*radtodeg, angle_max(2)*radtodeg
!   print *,angle_min(3)*radtodeg, angle_max(3)*radtodeg
!   print *,angle_min(4)*radtodeg, angle_max(4)*radtodeg

  dt = twopi/(np-1)
  
!  set up surface as ellipse
   print *,' size of grid = ', size(grid)

  do i=2,np
   theta = dt*(i-1)
   grid(i)%xs = a*cos(theta)
   grid(i)%ys = b*sin(theta)
   grid(i)%ts = theta
   grid(i)%V = 0.
   do j=1,4
     grid(i)%point(j) = .false.
     if(theta > bpm%min_angle(j) .and. theta < bpm%max_angle(j))then
      grid(i)%point(j) = .true.
!      print '(i,3e12.4)',j,xs(i),ys(i),theta*radtodeg
     endif
   end do
  end do

  return
end

  subroutine setup_extrusion_grid(grid,bpm,np, wire_radius)
  use precision_def
  use bpm_mod

  implicit none

  integer np
  real(rp) r0/75.0062/,center/50.0126/,flat/45.0088/ !, x0/12.98/,y0/24.92/
  type (grid_struct), allocatable :: grid(:)
  type (bpm_geometry_struct) bpm

  real(rp) theta, L_top, L_side, theta_0, deltat, deltay
  real(rp) button_radius/9.525/,x0/13.9954/, y0
  real(rp) t1,t2
  real(rp) xmin, xmax
  real(rp) wire_radius

  integer i,j
  integer np_top, np_side
  integer last
  integer np_1

  real(rp) twopi/6.28318530718/, radtodeg

  theta = asin(flat/r0)
  L_top = 2*r0*theta
  L_side = 2*(r0*cos(theta)-center)

  if(2*(np/2) == np) then
   print *,' NP = ',np,' NP must be odd '
   stop
  endif  
  np_1  = np-1
  np_top = L_top/(L_top+L_side) * np_1
  np_side = np_1-np_top
  if(2*(np_top/2) /= np_top)then
    np_top = np_top + 1
    np_side = np_side - 1
  endif
  deltat = L_top/(np_top/2)/r0
  wire_radius = L_top/(np_top/2)/2

   print *,' size of grid = ', size(grid)
  print *,' np = ', np
  print *,' L_top=', L_top
  print *,' L_side=', L_side
  print *, ' np_top = ', np_top
  print *, ' np_side = ', np_side
  print *, ' wire_radius = ', wire_radius

  do i=1,np
   do j=1,4
    grid(i)%point(j) = .false.
   end do
  end do
  t1 = asin((x0+button_radius)/sqrt(r0**2+(x0+button_radius)**2))
  t2 = asin((x0-button_radius)/sqrt(r0**2+(x0-button_radius)**2))
  xmax=r0*sin(t1)
  xmin=r0*sin(t2)
  bpm%min_angle(4)= twopi/4-t1
  bpm%max_angle(4)= twopi/4-t2
  bpm%min_angle(3)= twopi/4+t2
  bpm%max_angle(3)= twopi/4+t1
  bpm%min_angle(2)= 3*twopi/4+t2
  bpm%max_angle(2)= 3*twopi/4+t1
  bpm%min_angle(1)= 3*twopi/4-t1
  bpm%max_angle(1) = 3*twopi/4-t2
   print '(a8,a20,2e12.4)','   top  ',' bpm%..angle(3) = ', bpm%min_angle(3), bpm%max_angle(3)
   print '(a8,a20,2e12.4)','   top  ',' bpm%.._angle(4) = ', bpm%min_angle(4), bpm%max_angle(4)
   print '(a8,a20,2e12.4)',' bottom ',' bpm%.._angle(1) = ', bpm%min_angle(1), bpm%max_angle(1)
   print '(a8,a20,2e12.4)',' bottom ',' bpm%..._angle(2) = ', bpm%min_angle(2), bpm%max_angle(2)

! first do top
  theta_0 = -theta
  last = 1
  do i=last+1,np_top/2+last
    theta = theta_0+(i-last)*deltat
    grid(i)%xs = -r0*sin(theta)
    grid(i)%ys = r0*cos(theta)-center
    grid(i)%ts = theta +twopi/4
    grid(i)%V = 0.
    if(xmin < grid(i)%xs .and. grid(i)%xs < xmax)then
       grid(i)%point(4) = .true.
    endif
    if(-xmax < grid(i)%xs .and. grid(i)%xs < -xmin)then
       grid(i)%point(3) = .true.
    endif
    
  end do

! bottom
  last = np_top/2+last
  do i=last+1,np_top/2+last
    theta = theta_0+(i-last)*deltat + twopi/2
    grid(i)%xs = -r0*sin(theta)
    grid(i)%ys = r0*cos(theta)+center
    grid(i)%ts = theta + twopi/4
    grid(i)%V = 0.
    if(-xmax < grid(i)%xs .and. grid(i)%xs < -xmin)then
       grid(i)%point(1) = .true.
    endif 
    if(xmin < grid(i)%xs .and. grid(i)%xs < xmax)then
       grid(i)%point(2) = .true.
    endif

  end do

! left
   deltay = L_side/(np_side/2)
   y0 = L_side/2
   last = np_top/2+last
   do i=last+1,np_side/2+last
     grid(i)%xs = -flat
     grid(i)%ys = y0-(i-last)*deltay
     grid(i)%V = 0.
     grid(i)%ts = twopi/2-atan2(grid(i)%ys,flat)
      if(grid(i)%ts < 0.) grid(i)%ts = grid(i)%ts + twopi
   end do

! right
   last = np_side/2+last
!   print *, last, np_side/2
   do i=last+1,np_side/2+last
     grid(i)%xs = flat
!    print *,i,grid(i)%xs
     grid(i)%ys = -y0+(i-last)*deltay
     grid(i)%V = 0.
     grid(i)%ts = atan2(grid(i)%ys,flat)
      if(grid(i)%ts < 0.) grid(i)%ts = grid(i)%ts + twopi
   end do

  if(np /= np_side/2+last) then
    print *,' SETUP_GRID: np /= last point in grid '
    stop
  endif
  return
end