program test

 use bmad
 implicit none
 real(rp) x,y,width/0.05/,offset/0.05/, r0/0.001/, V0/-25000./
 real(rp) E(2)
 integer nwires/50/
integer i
 
 x=0.
 do i=-50,50
 y = i*0.001

  call esquad_field(V0,x,y,width,offset,nwires,E)
 print '(2es12.4)',E
 write(12, '(2es12.4,2x,2es12.4)') x,y,E
end do
end

! program to compute efield of electrostatic quadrupole



subroutine esquad_field(V0,x,y,width,offset,n,E)

 use bmad
 use nr

 implicit none

 type plate_struct
  real(rp) x,y
 end type

 type (plate_struct), allocatable, save :: left_plate(:), top_plate(:)
 real(rp) width, offset
 real(rp), allocatable, save :: Mpp(:,:), Mtr(:,:), M00(:,:), T(:,:),T_inv(:,:), Mbr(:,:)
 real(rp), allocatable, save :: Q(:), V(:)
 real(rp), allocatable :: Vcalc(:)
 real(rp) r0/0.0001/
 real(rp) V0
 real(rp) delta
 real(rp) x,y,E(2)
 real(rp) dx, dy, r2
 integer n
 integer i,j
 logical recompute/.true./

if(recompute)then
 recompute = .false.
 allocate(left_plate(-n:n), top_plate(-n:n))
 allocate(Mpp(-n:n,-n:n), Mtr(-n:n,-n:n),M00(-n:n,-n:n), Mbr(-n:n,-n:n))
 allocate(T(-n:n, -n:n), T_inv(-n:n, -n:n))
 allocate(Q(-n:n), V(-n:n))
 allocate(Vcalc(-n:n))
 delta = width/(2*n+1)

do i=-n,n
 left_plate(i)%x = -offset
 left_plate(i)%y = i*delta
 top_plate(i)%x = i*delta
 top_plate(i)%y = offset
 V(i) = V0
end do

!do i=-n,n
!  print '(11es12.4)',left_plate(i)%y
!end do

do i=-n,n
  do j=-n,n
    Mpp(i,j) = 0.5*log((left_plate(i)%x + left_plate(j)%x)**2 + (left_plate(i)%y - left_plate(j)%y)**2/r0**2)
    Mtr(i,j) = 0.5*log((-left_plate(i)%x  -top_plate(j)%x)**2 + (left_plate(i)%y - top_plate(j)%y)**2/r0**2)
    Mbr(i,j) = 0.5*log((-left_plate(i)%x  -top_plate(j)%x)**2 + (left_plate(i)%y + top_plate(j)%y)**2/r0**2)
    if(i == j) M00(i,j) = log(r0/r0)
    if(i/=j)M00(i,j) = 0.5*log((left_plate(i)%y - left_plate(j)%y)**2/r0**2)
 end do
end do

!print *
!do i=-n,n
! print '(11es12.4)',Mpp(i,-n:n)
!end do

T = Mpp-Mtr-Mbr + M00
call mat_inverse(T, T_inv)

Q = matmul(T_inv,V)
Vcalc = matmul(T,Q)

print *
do i=-n,n
 print '(11es12.4)',V(i)
end do

!print *
!do i=-n,n
! print '(11es12.4)',T(i,-n:n)
!end do


!print *
!print '(11es12.4)', Q
!write(11,'(es12.6)')Q

endif ! Q calculated

! compute field

E=0.
do i = -n,n
 do j=1,4
  if(j==1)then
   dx = x-left_plate(i)%x
   dy = y-left_plate(i)%y
  elseif(j==2)then
   dx = x-top_plate(i)%x
   dy = y-top_plate(i)%y
 elseif(j==3)then
   dx = x+left_plate(i)%x
   dy = y-left_plate(i)%y
  elseif(j==4)then
   dx = x-top_plate(i)%x
   dy = y+top_plate(i)%y
 endif
  r2 = dx**2+dy**2
  E(1) = E(1) + (-1)**(j+1)*Q(i)*dx/r2
  E(2) = E(2) + (-1)**(j+1)*Q(i)*dy/r2
 end do
end do

return
end