subroutine taylor_map(lat, taylor_order)
  use bmad

 implicit none
 type (lat_struct) lat
 type (taylor_struct) t_map(6)
 integer taylor_order
 logical err_flag
 real(rp) T(4,4,4,4)
 
 call transfer_map_calc(lat, t_map, err_flag)

 T(1,1,1,1) = taylor_coef(t_map(1),taylor_expn([1,1,1]))
 T(1,1,1,2) = taylor_coef(t_map(1),taylor_expn([1,1,2]))
 T(1,1,2,2) = taylor_coef(t_map(1),taylor_expn([1,2,2]))
 T(1,2,2,2) = taylor_coef(t_map(1),taylor_expn([2,2,2]))
 T(2,1,1,1) = taylor_coef(t_map(2),taylor_expn([1,1,1]))
 T(2,1,1,2) = taylor_coef(t_map(2),taylor_expn([1,1,2]))
 T(2,1,2,2) = taylor_coef(t_map(2),taylor_expn([1,2,2]))
 T(2,2,2,2) = taylor_coef(t_map(2),taylor_expn([2,2,2]))
 T(3,3,3,3) = taylor_coef(t_map(3),taylor_expn([3,3,3]))
 T(3,3,3,4) = taylor_coef(t_map(3),taylor_expn([3,3,4]))
 T(3,3,4,4) = taylor_coef(t_map(3),taylor_expn([3,4,4]))
 T(3,4,4,4) = taylor_coef(t_map(3),taylor_expn([4,4,4]))
 T(4,3,3,3) = taylor_coef(t_map(4),taylor_expn([3,3,3]))
 T(4,4,4,4) = taylor_coef(t_map(4),taylor_expn([3,3,4]))
 T(4,3,4,4) = taylor_coef(t_map(4),taylor_expn([3,4,4]))
  T(4,4,4,4) = taylor_coef(t_map(4),taylor_expn([4,4,4]))
 print '(a,es12.4)', 'T(1,1,1,1) = ' , T(1,1,1,1)
 print '(a,es12.4)', 'T(1,1,1,2) = ' , T(1,1,1,2)
 print '(a,es12.4)', 'T(1,1,2,2) = ' , T(1,1,2,2)
 print '(a,es12.4)', 'T(1,2,2,2) = ' , T(1,2,2,2)
 print '(a,es12.4)', 'T(2,1,1,1) = ' , T(2,1,1,1)
 print '(a,es12.4)', 'T(2,1,1,2) = ' , T(2,1,1,2)
 print '(a,es12.4)', 'T(2,1,2,2) = ' , T(2,1,2,2)
 print '(a,es12.4)', 'T(2,2,2,2) = ' , T(2,2,2,2)
 print '(a,es12.4)', 'T(3,3,3,3) = ' , T(3,3,3,3)
 print '(a,es12.4)', 'T(3,3,3,4) = ' , T(3,3,3,4)
 print '(a,es12.4)', 'T(3,3,4,4) = ' , T(3,3,4,4)
 print '(a,es12.4)', 'T(3,4,4,4) = ' , T(3,4,4,4)
 print '(a,es12.4)', 'T(4,3,3,3) = ' , T(4,3,3,3)
 print '(a,es12.4)', 'T(4,3,3,4) = ' , T(4,3,3,4)
 print '(a,es12.4)', 'T(4,3,4,4) = ' , T(4,3,4,4)
 print '(a,es12.4)', 'T(4,4,4,4) = ' , T(4,4,4,4)
 print '(a,es12.4)',' Delta nux/x^2 = (T_1111 + T_2211)/2pi = ', (T(1,1,1,1)+T(2,1,1,2))/twopi
  print '(a,es12.4)',' Delta nuy/y^2 = (T_3333 + T_4433)/2pi = ', (T(3,3,3,3)+T(4,3,3,4))/twopi
 return 
 end