program sextupole_geometry

use bmad
implicit none

real(rp) rin/3./,rout/6./,w/2./, r/8./
real(rp) theta, gamma,delta
real(rp) x(100),y(100), rad1, rad2
integer i,j

i=0
theta = -twopi/4

 i=i+1
 x(i) = r * cos(theta)
 y(i)= r * sin(theta)

 i=i+1
 x(i) = rin * cos(theta)
 y(i) = rin * sin(theta)

 i=i+1
 delta = atan(w/2./rin)
 rad1 = rin**2+(w/2)**2
 x(i) =  sqrt(rad1)* cos(theta+delta)
 y(i) =  sqrt(rad1)* sin(theta+delta)

 i=i+1
 gamma = atan(w/2,rout)
 rad2 = rout**2+(w/2)**2
 x(i) = sqrt(rad2)*cos(theta+gamma)
 y(i) = sqrt(rad2)*sin(theta+gamma)

 theta=-twopi/12
 do while (theta <= twopi/12)
 
 i=i+1
 gamma = atan(w/2,rout)
 x(i) = sqrt(rad2)*cos(theta-gamma)
 y(i) = sqrt(rad2)*sin(theta-gamma)
  
 i=i+1
 delta = atan(w/2./rin)
 x(i) = sqrt(rad1)*cos(theta-delta)
 y(i) = sqrt(rad1)*sin(theta-delta)

 i=i+1
 delta = atan(w/2./rin)
 x(i) = sqrt(rad1)*cos(theta+delta)
 y(i) = sqrt(rad1)*sin(theta+delta)

 i=i+1
 gamma = atan(w/2,rout)
 x(i) = sqrt(rad2)*cos(theta+gamma)
 y(i) = sqrt(rad2)*sin(theta+gamma)

 theta = theta + twopi/6
end do

 theta = twopi/4
 
 i=i+1 
 gamma = atan(w/2,rout)
 x(i) = sqrt(rad2)*cos(theta-gamma)
 y(i) = sqrt(rad2)*sin(theta-gamma)


 i=i+1
 delta = atan(w/2./rin)
 x(i) = sqrt(rad1)*cos(theta-delta)
 y(i) = sqrt(rad1)*sin(theta-delta)

 i=i+1
 x(i) = rin * cos(theta)
 y(i) = rin * sin(theta)

 i=i+1
 x(i) = r*cos(theta)
 y(i) = r* sin(theta)

 do j=1,i
 write(11,'(i10,1x,2es12.4)')j,x(j),y(j)
 end do
end program