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

  implicit none

  integer np
  real(rp) r0/44.45/,deep_but_r0/46.482/,shallow_but_r0/44.45/,flat/27.94/,but_diam/19.0/ 
  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/
  real(rp) t1,t2
  real(rp) xmin, xmax
  real(rp) wire_radius
  real(rp) piby4, theta_but, theta_flat_shallow, theta_flat_deep, dep, theta_flat
  real(rp) x,y
  real(rp) tmin
  real(rp) r_deep, r_small

  integer i,j, jmin, jbut

  real(rp) twopi/6.28318530718/

  theta_but = asin(but_diam/r0)/2
  theta_flat = asin(flat/r0/2)

  if(2*(np/2) == np) then
   print *,' NP = ',np,' NP must be odd '
   stop
  endif  

   print *,' size of grid = ', size(grid)
  print *,' np = ', np
  print *, ' wire_radius = ', wire_radius

  deltat = twopi/np

  do i=1,np
   do j=1,4
    grid(i)%point(j) = .false.
   end do
  end do

  do i=1,np
    theta = i*deltat
     tmin = 99.
     do j =0,8
       if(abs(theta- j*twopi/8) < tmin)then
         tmin = abs(theta-j*twopi/8)
         jmin = j
       endif
     end do
     jbut = jmin+1
     if((jbut/2)*2 == jbut)grid(i)%point(jbut/2) = .true.
     if (jmin == 2 .or. jmin == 6)then
         r_deep = deep_but_r0
         r_small = r0
        else
         r_deep = r0
         r_small = r0
     endif 

     call notch(theta,jmin*twopi/8, r_small, flat, r_deep, x,y)

    grid(i)%xs = x
    grid(i)%ys = y
    grid(i)%ts = theta
    grid(i)%V = 0.
    
  end do

  return
end

subroutine notch(theta, theta_ref, r_small, width,r_deep, x,y)

use precision_def
implicit none

real(rp) theta, theta_ref, r_small, width, r_deep, x, y
real(rp) theta_small, theta_big, dt
real(rp) xout, yout, M11, M12, M21, M22

  theta_small = 2*atan(width/2/r_deep)
  theta_big = 2*atan(width/2/r_small)

  dt = theta - theta_ref

   if(abs(dt) < theta_small)then
     x=r_deep/cos(theta)
     y = r_deep * tan(theta)
   endif
   if(abs(dt) > theta_small .and. abs(dt) < theta_big)then
     y = width/2
     x = y/tan(theta)
   endif

   M11 = cos(theta_ref)
   M12 = sin(theta_ref)
   M21 = -M12
   M22= M11

   xout = M11 * x + M12 * y
   yout = M12 * x + M22 * y

   x = xout
   y = yout

  return
 end