subroutine generalize_cyl_exp(azimuthal_exp_z, azimuthal_exp_r, r,phi,z, field)
 use bmad
 
 implicit none
 type(em_field_struct) field
 real(rp) r,phi,z
 real(rp) a(0:200), b(0:200), c(0:200), d(0:200), x0(0:200)
 real(rp), save :: AA(0:200), BB(0:200), CC(0:200), DD(0:200)
 real(rp) first_zero_at/8.1/
 real(rp) Br,Bz,Bphi,x,jn,jp, jb(0:200)
 real(rp), save::k(0:200)
 real(rp) r0/7.112/
 character *120 azimuthal_exp_z, azimuthal_exp_r
 integer, save :: nmodes
 integer m
 logical first/.true./

 if(first)then
  call read_fourier(azimuthal_exp_z,nmodes,a,b)
  call read_fourier(azimuthal_exp_r,nmodes,c,d)
  print '(a,a)',' read ',azimuthal_exp_z
  print '(a,a)',' read ',azimuthal_exp_r
  nmodes = 10
  call compute_bessel_zeros(x0, nmodes)
print '(a)',' Initialize generalize_cyl_exp'
write(96, '(a)')' Initialize generalize_cyl_exp'
  k(0:nmodes) = x0(0:nmodes)/first_zero_at
  do m=0,nmodes
   jb(m) = bessel_jn(m,k(m)*r0)
   AA(m) = a(m)/jb(m)/k(m)
   BB(m) = b(m)/jb(m)/k(m)
   jp = m/r0/k(m)*jb(m)-bessel_jn(m+1,k(m)*r0)
   CC(m) = c(m)/jp/k(m)
   DD(m) = d(m)/jp/k(m)
 print '(i10,1x,4es12.4)',m, a(m), b(m), c(m), d(m)
 write(6, '(i10,1x,4es12.4)')m, a(m), b(m), c(m), d(m)
  end do


 first=.false.
endif


 Br=0
 Bphi=0
 Bz=0
! print *
 do m=0,nmodes
  x = k(m)*r
  jn = bessel_jn(m,x)
  jp = m/x*jn-bessel_jn(m+1,x)
  Bz = Bz + jn * k(m) &
       * (cosh(k(m)*z)* (AA(m)*sin(m*phi)+BB(m)*cos(m*phi)) - sinh(k(m)*z)*(CC(m)*sin(m*phi)+DD(m)*cos(m*phi)) )
  Br = Br + jp * k(m) &
       * (sinh(k(m)*z)* (AA(m)*sin(m*phi)+BB(m)*cos(m*phi)) + cosh(k(m)*z)*(CC(m)*sin(m*phi)+DD(m)*cos(m*phi)))
  Bphi = Bphi + jn /r * m &
       * (sinh(k(m)*z)* (AA(m)*cos(m*phi)-BB(m)*sin(m*phi)) +cosh(k(m)*z)*(CC(m)*cos(m*phi)-DD(m)*sin(m*phi)))
!  print *,m,Bz
 end do
 field%B(1) = Br
 field%B(2) = Bz
 field%B(3) = Bphi

return

end subroutine generalize_cyl_exp

subroutine compute_bessel_zeros(x0,nmodes)
 use bmad
 use nr
 implicit none
 real(rp) a(0:200), b(0:200), c(0:200), d(0:200), x0(0:200), x(0:200)
 real(rp) jx(0:200), jxp(0:200), delta(0:200), dx(0:200)
 integer nmodes, i
 integer m/0/

delta(:) = 0.00001
jx(0:nmodes) = 100.
! compute first zeros of J_n
  x0(0) = 2.405
  x0(1) = 3.832
  x0(2) = 5.136
  x0(3) = 6.3802
  x0(4) = 7.5883
  x0(5) = 8.7715

  forall (i=6:nmodes)x0(i) = (i-0.5)*twopi/4  !first guess
! print *
!  do i=1,nmodes
!     print '(i10,2es12.4)',i,x0(i), bessel_jn(i,x0(i)) 
!     jx(i) = bessel_jn(i,x0(i)) 
!  end do
!  write(11,'(i10,1x,22es12.4)')m,x0(0:10), jx(0:10)

 do while (any(abs(jx(0:nmodes)) > 0.001)) 
  m=m+1
  
  do  i=0,nmodes
   jx(i) = bessel_jn(i,x0(i))
   jxp(i) = (bessel_jn(i,x0(i)+delta(i)) - jx(i))/delta(i)
   dx(i) = -jx(i)/ jxp(i)
   x0(i) = x0(i) + dx(i)/2
  end do
!  print *
!   do i=1,nmodes
!      print '(i10,2es12.4)',i,x0(i), bessel_jn(i,x0(i)) 
!      jx(i) = bessel_jn(i,x0(i))
!   end do
!   write(11,'(i10,1x,22es12.4)')m,x0(0:10), jx(0:10)
end do
end subroutine compute_bessel_zeros



 subroutine read_fourier(azimuthal_exp,nmodes,a,b)
  use precision_def
  implicit none
  real(rp) a(0:200), b(0:200)
  real(rp) x, ab(0:400) 
  integer lun, lunget
  integer mode, ix, nmodes
  character*120 azimuthal_exp
  character*100 string, word, equal, value

  a(:)=0
  b(:)=0
  lun = lunget()
  open(unit=lun, file = azimuthal_exp, status = 'old')
  do while(.true.)
    read(lun, '(a)', end=99) string

    call string_trim(string, word, ix)  
    if(word(1:1) /= 'p')cycle

    call string_trim(word(ix+1:),equal, ix)
    if(equal(1:1) /= '=')cycle
!    print *, 'string = ', trim(string), '   word = ', word(1:ix), ' equal = ', equal(1:1)
    call string_trim(equal(ix+1:), value, ix)
    read(value(1:ix),*)x
   
    if(word(3:4) == '  ')then
        read(word(2:2),*)mode
      else if(word(4:4) == ' ')then
        read(word(2:3),*)mode
      else
       read(word(2:4),*)mode
    endif
!    print *,' word(1:4) =', word(1:4),' mode = ', mode, x
    if(2*(mode/2) == mode .and. mode /= 0)then
      ix=mode/2
      a(ix) = x
     else
      ix = (mode+1)/2
      b(ix)=x
    endif
    ab(mode) = x
!    print '(i10,es12.4)', mode, ab(mode)

 end do
99 continue
close(unit = lun)

!do ix=0,(mode+1)/2
    nmodes = ix
!   print '(i10,2es12.4)', ix, a(ix), b(ix)
!end do

return
end