program test

 use bmad
  implicit none

 integer i, number/100/, nquads/10/, n
 integer one/1/,zero/0/,minus_one/-1/
 integer ix
 real(rp) radius/100/,dt, t(0:110)
 real(rp) x(0:110),y(0:110)
 real(rp) theta, rquad, rbend
 real(rp) arc_length, width, a, b

!create ring

 dt = twopi/number 
 theta = -dt
 n=0
 do while (theta <= twopi)
  theta = theta+dt
  a = radius*cos(theta)
  b = radius*sin(theta)
  x(n) = radius*cos(theta)
  y(n) = radius*sin(theta)
  t(n) = theta
 write(11,'(2f12.4)')x(n), y(n)
  n=n+1
end do

 rquad = 2*twopi/nquads/5*radius
 rbend = twopi/nquads/radius

 width = 4
 do i=1,nquads
  ix = (i-1)*number/nquads
  rquad=2*twopi/nquads/5*radius
  arc_length = twopi/8
  call focus_lens(rquad*1.5,arc_length/1.5, width,x(ix),y(ix),t(ix),one)
  call focus_lens(rquad*3,arc_length/3, width,x(ix+2),y(ix+2),t(ix+2),one)
  ix = ix + number/nquads/2+1
  arc_length = twopi/(2*nquads+3)
  call focus_lens(radius,arc_length,width, x(ix),y(ix),t(ix),zero)
  one = -one
end do

end

subroutine focus_lens(radius,arc_length,width, x0,y0,tilt_angle, fdfb)

  use bmad
  use plots_interface

  implicit none
  real(rp) size, x0,y0, tilt_angle, arc, radius, darc, theta
  real(rp) arc_length, width
  real(rp) x(0:110), y(0:110), t(0:110)
  real(rp) tilt_offset
 
  integer i,n, fdfb
  integer p

!  arc = twopi/8
  arc=arc_length
  darc= 2*arc/20
  theta = -arc 
  n=-1
print '(/,a)',' radius, arc_length, width, x0, y0'
print '(5es12.4)', radius, arc_length, width, x0, y0

 if(fdfb /= 0)then
 if(fdfb>0)then !convex lens
    do while(theta <= arc)
   n=n+1
   x(n) = radius * cos(theta) - radius*cos(arc)
   y(n) = radius* sin(theta)
   theta = theta + darc
  end do
   x(2*n+1) = x(0)
   y(2*n+1) = y(0)
   p =2*n+1
 endif

 if(fdfb<0)then  !concave lens
  do while(theta <= arc)
   n=n+1
   x(n) = radius * cos(theta) - radius*cos(arc)
   x(n) = x(n) - 1.25*radius*(1-cos(arc))
   y(n) = radius* sin(theta)
   theta = theta + darc
  end do
 endif

! left half
print *,' n= ', n
  do i=1,n+1
   x(i+n) = -x(n-i+1)
   y(i+n) = y(n-i+1)
  end do 
   x(2*n+2) = x(0)
   y(2*n+2) = y(0)
  p = 2*n +2
tilt_offset = twopi/4
endif

if(fdfb==0)then !bend
  darc=arc_length/20.
  theta = -arc_length/2
  n=-1
  do while (theta <= arc_length/2)
   n=n+1
   x(n) = (radius+width)*cos(theta) - radius
   y(n) = (radius+width)*sin(theta)
   t(n) = theta
   theta = theta + darc
  end do
  do i=1,n+1
   x(i+n) = (radius-width)*cos(t(n-i+1))-radius
   y(i+n) = (radius-width)*sin(t(n-i+1))
  end do
   x(2*n+2) = x(0)
   y(2*n+2) = y(0)
   p = 2*n+2
tilt_offset = 0
endif  


do i=0,2*n
  print '(i10,2es12.4)',i,x(i),y(i)
end do
print *
  call rotate(tilt_angle+tilt_offset, x,y,p)
do i=0,p
  print '(i10,2es12.4)',i,x(i),y(i)
end do

   write(11,'(/,a,/)')'#'

  write(11,'(2f12.4)')(x(i)+x0,y(i)+y0,i=0,p)
return
end

subroutine rotate(theta,x,y,m)
use precision_def

real(rp) theta, x(0:), y(0:), a(0:m), b(0:m)
integer m,i

do i=0,m
  a(i) = x(i)*cos(theta)-y(i)*sin(theta)
  b(i) = y(i)*cos(theta)+x(i)*sin(theta)
  print '(4es12.4)',x(i),y(i),a(i),b(i)
end do


x(0:m) = a(0:m)
y(0:m) = b(0:m)

return
end