subroutine ptc_calc_fps_param (lat, map_order, param)

  use bmad
  use ptc_interface_mod
  use ptc_layout_mod
  use pointer_lattice, latptc => lat, dummy => pi, dummy1 =>twopi
  use c_TPSA, only: operator(.sub.)

  implicit none
  
  type(lat_struct), target :: lat
  integer map_order
  
  type(probe) xs0 
  type(probe_8) xs
  type(layout), pointer :: ring
  type(fibre), pointer:: p, ptc_fibre
  type(internal_state), target :: state

  type(c_damap) id, one_turn_map
  type(c_normal_form) normal_form

  real(rp) closed(6)
  
  !resonant map calculation
  type(c_damap) n_c
  type(c_universal_taylor) h_res
  integer i, j, fact
  complex(rp) del, alpha_xx_half, gamma
  
  real(rp) g, phi0

  real(rp) n_res/3.0/
  real(rp) param(3)
  
  param(:) = 0.0_rp
  
  !
 
  use_info = .true.
  use_quaternion=.true.
  
  n_cai=-i_                             !  J= x_1 x_2 (h+ h-)    h+=(x - i px)/sqrt(2)
 
  call ptc_ini_no_append
  
  switch_to_drift_kick=.false.
  check_excessive_cutting=.false.
  
  call lat_to_ptc_layout (lat) 
  ring => lat%branch(0)%ptc%m_t_layout
  p=>ring%start

  closed=0
  my_estate=>state
  state=only_4d0 ! care only (x,px,y,py), ignore (z,pz). In PTC (5,6) defined as (e,t)
  
  call init(state,map_order,0)
  call find_orbit_x(closed,state,1.e-7_rp,fibre1=p)

  !!  In PTC all objects that contain a TPSA variable 
  !!  or might contain one (..._8)  must be allocated and killed on exit

  call alloc(xs)
  call alloc(id,one_turn_map)
  call alloc(normal_form)  
  call alloc(n_c)
  
  xs0=closed 
  id=1       
  xs=xs0+id  ! FPP -> PTC, closed orbit added

  call propagate(xs,state,fibre1=p)
  one_turn_map=xs ! PTC -> FPP 
  
! get normal form
  call c_normal(one_turn_map, normal_form)

! check whether it is 3nux or 4nux
  if ( abs(normal_form%tune(1)*3 - nint(normal_form%tune(1)*3)) .lt. 0.08) then
     n_res = 3.0_rp
  elseif( abs(normal_form%tune(1)*4 - nint(normal_form%tune(1)*4)) .lt. 0.1) then
     n_res = 4.0_rp
  end if

 ! get resonant map for the coefficients
  normal_form%nres=(c_%no+1)/n_res
  do i=1, (c_%no+1)/n_res
     normal_form%m(1,i)=n_res*i ! keep only the 3vx or 4vx term
  enddo

  call c_normal(one_turn_map,normal_form, canonize=.true.)
  call clean(normal_form%h_nl,normal_form%h_nl,prec=1.d-10)
  
  
  fact=nint(normal_form%tune(1)*n_res)  ! determine near which resonance line (1/n or (n-1)/n)
  
  normal_form%h_l=0
  normal_form%h_l%v(1)=( i_*twopi*fact/n_res)*dz_c(1) 
  normal_form%h_l%v(2)=(-i_*twopi*fact/n_res)*dz_c(2)
  normal_form%h_l%v(3)=( i_*twopi*normal_form%tune(2))*dz_c(3) 
  normal_form%h_l%v(4)=(-i_*twopi*normal_form%tune(2))*dz_c(4) 

  normal_form%h_l=ci_phasor()*normal_form%h_l
  n_c=normal_form%atot**(-1)*one_turn_map*normal_form%atot
  id=exp(normal_form%h_l)
  n_c=id*n_c
  n_c=ci_phasor()*n_c*c_phasor()

  normal_form%h=c_logf_spin(n_c)  ! get Harmitonian

  call d_field_for_demin(normal_form%h, h_res)
  call clean(h_res, h_res, prec=1.d-5)

  del=h_res .sub. [1,1,0,0]
  alpha_xx_half= h_res .sub. [2,2,0,0]

  if (n_res .eq. 3.0_rp) then
     gamma= h_res .sub. [3,0,0,0]
  elseif (n_res .eq. 4.0_rp) then
     gamma= h_res .sub. [4,0,0,0]
  end if
  phi0=atan2(aimag(gamma),real(gamma) )
  g=abs(gamma)
  
  param(1) = real(alpha_xx_half)
  param(2) = real(gamma)
  param(3) = aimag(gamma) 

!  write(6,'(a,1es15.5,a,2es15.5)') "phi0: ", phi0, ' g: ', real(gamma), aimag(gamma)
  
  call kill(xs)
  call kill(id,one_turn_map)
  call kill(normal_form)  
  call kill(n_c)

end subroutine ptc_calc_fps_param