module muon_interface
interface
  subroutine create_phase_space(lat, nmuons, create_new_distribution, new_file, muon_file, &
        epsdistr, epsx, epsy, tdistr, tlength, tsigma, pzdistr, pz, pzsigma, twiss2, muons, twiss_ref,mat, mat_inv, nmuon_first)
    USE nr
!    use bmad
    use muon_mod

   IMPLICIT NONE
  type (lat_struct) lat
  integer nmuons
  integer nmuon_first
  character*16 epsdistr, tdistr, pzdistr, inf_aperture, twiss_ref
  character*120 new_file, muon_file
  real(rp) tlength, tsigma, pz, pzsigma, epsx, epsy
  real(rp), optional :: mat_inv(6,6),mat(6,6)
  logical inf_end_us, inf_end_ds, create_new_distribution
  logical energyloss
  type (muon_struct), allocatable :: muons(:)
  type (g2twiss_struct) twiss2
  end subroutine
end interface

interface
  subroutine compute_moments(muons, turn,moment, nmuons, nmu, unit, ele_name, ix)
  use muon_mod
  implicit none

  real(rp) turn
  integer n,i,nmuons, nmu, ix
  integer unit
  type (muon_struct), allocatable :: muons(:)
  type (g2moment_struct) moment
  character*(*) ele_name
  end subroutine
end interface

interface
subroutine compute_spins(muons,turn,nmuons,unit,ele_name,ix,ele_s)

  USE bmad
  USE muon_mod
!  USE spin_mod
  IMPLICIT NONE

  integer  unit
  integer  n,i,j,nmuons,ix
  real(rp) turn,ele_s
  
  type (muon_struct), dimension(nmuons) :: muons   !defined this way to avoid some execptions
  character*(*) ele_name
 end subroutine
end interface

INTERFACE
  SUBROUTINE compute_beam_params(muons,verbosity, where)
    USE muon_mod
    IMPLICIT NONE
    ! Subroutine arguments
    TYPE (muon_struct), ALLOCATABLE, INTENT(IN) :: muons(:)
    INTEGER, OPTIONAL :: verbosity ! controls how much information is printed
    character*(*), optional:: where
  END SUBROUTINE compute_beam_params
END INTERFACE

interface
 subroutine collect_times(muons, nturns, NN, bin_max)
  use muon_mod
  use parameters_bmad
  implicit none

  real(rp), allocatable :: NN(:)
  integer nturns, bin_max
  type (muon_struct), allocatable :: muons(:)
  end subroutine
end interface

interface
 subroutine write_phase_space(nmuons, muons, string, tot, turn, BetaCrossE_start_time)

  use muon_mod
  implicit none
  integer nmuons
  integer i
  integer tot
  integer turn
  character*(*) string
  type (muon_struct), allocatable :: muons(:)
  real(rp), optional :: BetaCrossE_start_time
 end subroutine
end interface

interface
 subroutine inflector_end_scatter(nmuons, muons)

   use muon_mod
   use parameters_bmad
   implicit none
   integer unit, nmuons
   integer i, tot
   type (muon_struct), allocatable :: muons(:)

 end subroutine
end interface

interface
  function fnalw(center,width)
    USE parameters_bmad ! fnalw_integral table
    IMPLICIT NONE

    REAL(rp)             :: fnalw                       ! quantity to determine
    REAL(rp), INTENT(IN) :: center, width               ! center and width of distribution [a.u.]
    INTEGER,  PARAMETER  :: nbins=100                   ! number of bins in fnalw_integral table
    REAL(rp) rand, di, frac_low, frac_high, i_interp    ! helper variables
    INTEGER  i
  end function fnalw
end interface

interface
subroutine read_phase_space(unit, nmuons, muons, string, tot, nmuon_first)
 use muon_mod
 implicit none
 integer unit, nmuons
 integer nmuon_first
 integer i, tot
 character*(*) string
 type (muon_struct), allocatable :: muons(:)
 end subroutine
end interface

interface
 subroutine twiss_propagate_cm(lat, nbranch)
  use bmad
  implicit none
  type (lat_struct) lat
  integer nbranch
 end subroutine
end interface

interface
 subroutine read_Vmuons(unit, nmuons, muons,toff, string, tot)
  use muon_mod
  implicit none
  integer unit, nmuons
  integer i, tot
  character*(*) string
  type (muon_struct), allocatable :: muons(:)
  real(rp) toff
 end subroutine
end interface

interface
 subroutine read_Dmuons(lunit, nmuons, muons,toff, string, tot, nmuon_first)
  use muon_mod
  implicit none
  integer lunit, nmuons, nmuon_first
  integer i, tot
  character*(*) string
  type (muon_struct), allocatable :: muons(:)
  real(rp) toff
 end subroutine
end interface

interface
 subroutine inflector_scatter(nmuons, muons, inf_end_us, inf_end_ds,energy_loss)
  use parameters_bmad
  use muon_mod
  use materials_mod
  IMPLICIT NONE
  type (muon_struct), allocatable :: muons(:)
  integer i,nmuons
  logical inf_end_us, inf_end_ds, energy_loss
 end subroutine
end interface

interface
 subroutine check_inflector_aperture(nmuons,muons, inf_aperture, lost_at_inflector)
  use parameters_bmad
  use muon_mod
  use materials_mod
  IMPLICIT NONE
  type (muon_struct), allocatable :: muons(:)
  integer i,nmuons, lost_at_inflector
  logical withinInflectorAperture
  character*16 inf_aperture
 end subroutine
end interface

interface
 subroutine create_inv_taylor_ele(ele)
  use bmad
  implicit none
  type (ele_struct) ele
  real(rp) m(6,6)
 end subroutine
end interface

interface
 subroutine compute_emittance_beta(nmuons,muons, twiss1, epsx,epsy,averages)
  use bmad
  use muon_mod
  implicit none
  type (muon_struct), allocatable :: muons(:) 
  type (averages_struct) averages
  type (g2twiss_struct) twiss1
  real(rp) epsx, epsy
  integer nmuons, i
 end subroutine
end interface

interface
  SUBROUTINE fields(ele,x,y,s,t,field)
  use bmad
  use parameters_bmad
  IMPLICIT NONE
  type (ele_struct) ele
  type (em_field_struct) field
  real(rp) x,y,s,t
  end subroutine fields
end interface

interface
 subroutine get_field(y,map,B,b1,b2,out_of_range)
   use magfield
   implicit none
   logical out_of_range
   type(magfield_struct), allocatable :: map(:,:,:)
   real(rp) x(3), B(3) ,b1(3,3),b2(3,3)
   real(rp) y(3)
 end subroutine
end interface

interface
 subroutine optimize_incident(start_orb,lat, inflector_end_target, opt_orb, inflector_angle)
  use bmad
  use parameters_bmad
  implicit none
  type (lat_struct) lat
  type (coord_struct) start_orb, opt_orb
  type (initial_offsets_struct) inflector_end_target
  real(rp) inflector_angle
 end subroutine
end interface

interface
 subroutine optimize_twiss(lat, twiss, twiss_opt)
  use bmad
  use muon_mod
  implicit none
  type (lat_struct) lat
  type (g2twiss_struct) twiss, twiss_opt
 end subroutine
end interface

interface
  subroutine set_kicker_params(lat, kicker_params)
  use bmad
  use muon_mod
  use parameters_bmad

 implicit none
  type (lat_struct), target::lat
  type (kicker_params_struct) kicker_params
 end subroutine
end interface

interface
  subroutine MatchToRing(lat,initial_offsets,kicker_params,twiss_match)

  use bmad
  use muon_mod
  use parameters_bmad

  implicit none

  type(lat_struct), target:: lat 
  type(kicker_params_struct) kicker_params
  type (g2twiss_struct) twiss_match
  type (coord_struct), allocatable :: co(:), to_orbit(:)
  type (initial_offsets_struct) initial_offsets
 end subroutine
end interface

interface
 subroutine focus_ele(ele,x,y,t,field)
 use bmad
 implicit none
 type(ele_struct) ele
 type(em_field_struct) field
 real(rp) x,y,t
 end subroutine
end interface

interface
  subroutine step_floor_around_ring(lat,j, nbranch, co_start, vec_init, s_start)
    use bmad
    implicit none

    type (lat_struct) lat
    type (coord_struct) co_start
!    type (coord_struct), allocatable :: co(:)
    real(rp) vec_init(3)
    real(rp) s_start      
    integer j, nbranch
  end subroutine
end interface

interface
  subroutine track_write_backleg_trajectory(lat,co)
  use bmad
  implicit none
  type(lat_struct) lat
  type(coord_struct), allocatable:: co(:)
  end subroutine
end interface

interface
  subroutine write_single_particle_tbt(i,muons,Q, nturns, ref_time)
       use bmad
       use muon_mod
       implicit none
       integer i
       type(muon_struct), allocatable :: muons(:)
       integer nturns
       real(rp) Q(3)
       real(rp) ref_time
  end subroutine
end interface

interface
  subroutine compute_moments_vs_time(nturns, turn,nmuons, muons, string, coef_save, chisq_save, ndf_save, test)
   use muon_mod
   implicit none
   integer  nmuons
   integer i, turn
   integer nturns
   integer, optional :: ndf_save   
   character*(*) string
   integer, optional :: test

   type (muon_struct), allocatable :: muons(:)
   real(rp) time_bin_width
   real(rp), optional :: coef_save(4)
   real(rp), optional :: chisq_save
   
  end subroutine
end interface

interface
  subroutine write_lost_muons(n,s, orbit, ele_name,i)
   use bmad
   use parameters_bmad
   implicit none
   type (coord_struct) orbit
   integer n,i
   real(rp) s
   character*16 ele_name
  end subroutine
end interface

interface
       subroutine first_turn(lat,spin_tracking_on, from_orbit, circumference, make_movie, start_tracking_at_inflector_exit)
         use bmad
         implicit none
         type (lat_struct) lat
         type (coord_struct), allocatable :: from_orbit(:)
         real(rp) circumference
         logical spin_tracking_on, make_movie
         logical start_tracking_at_inflector_exit
 end subroutine
end interface

interface
       subroutine write_single_particle_by_element(i, ele_name, coord, s_frac,s_tot, Q) 
        use bmad
        use parameters_bmad
        implicit none
        type(coord_struct) coord
        integer i
        real(rp) s_frac, s_tot  ! fractional turn, total turns
        real(rp) Q(3)
        character*16 ele_name
 end subroutine
end interface


interface
 subroutine set_steering(lat, name, kick)
  use bmad
  use muon_mod

 implicit none
  type (lat_struct), target::lat
  character*16 name(:)
  real(rp) kick(:)

 end subroutine
end interface

interface
subroutine set_a_quad(ele, verbose)
  use bmad
  use muon_mod
  use parameters_bmad
  implicit none
  type (ele_struct), target :: ele
  logical verbose, err
end subroutine
end interface

interface
subroutine generate_time_long_dist(tdistr,muons, nmuons, tlength, tsigma)

  use muon_mod
  use parameters_bmad

  implicit none
  type (muon_struct), allocatable :: muons(:)
  character*16 tdistr
  integer nmuons
  integer ix, pulse
  real(rp) tlength, tsigma
end subroutine
end interface

interface
subroutine set_dipole_params(lat,nbranch, DeltaB_onB)
use bmad
use parameters_bmad

  type (lat_struct), target::lat
  type(ele_struct) ele
  integer nbranch
  real(rp) DeltaB_onB
end subroutine
end interface

interface
subroutine track_all_int_efield (lat, orbit,BetaCrossE_start_time,  sumBetaCrossE, ix_branch, track_state,  err_flag, orbit0, verbose)

 use bmad_interface
 use muon_mod

 implicit none

 type (lat_struct), target :: lat
 type (coord_struct), allocatable, target :: orbit(:)
 type (coord_struct), optional, allocatable, target :: orbit0(:)
 type (sumBetaCrossE_struct), allocatable :: sumBetaCrossE(:)
 integer, optional :: ix_branch, track_state
 real(rp) BetaCrossE_start_time
 logical, optional :: err_flag
 logical, optional :: verbose
end subroutine
end interface

  interface
        subroutine write_BetaCrossE_turn(n, co, sumBetaCrossE, muons, BetaCrossE_start_time)

        use bmad
        use parameters_bmad
        
        implicit none
        type (coord_struct) co
        type(sumBetaCrossE_struct)sumBetaCrossE
        type (muon_struct) muons
        integer n
        real(rp) sum_vec(3), sum_vec_pitch(3), path, pathx
        real(rp) sum_time
        real(rp) BetaCrossE_start_time
      end subroutine write_BetaCrossE_turn
   end interface
   
   interface
          subroutine create_and_propagate_electron(lat,nbranch, co_muon_end, co_electron, n, turn, nmuons, make_movie)
               use bmad
               implicit none
               type (lat_struct) , target ::lat
               type (coord_struct) co_muon_end, co_electron
               integer n, turn, nmuons
               integer nbranch, track_state
               logical make_movie, err_flag
             end subroutine create_and_propagate_electron
      end interface
      interface 
       subroutine decay_info(turn,n,ecoord,coord, muons)
         use bmad
         use parameters_bmad
         use muon_mod
         IMPLICIT NONE
         real(rp) turn
         integer  n             !which muon
         type (coord_struct) :: coord,ecoord ! for muon and electron
         type (muon_struct), optional :: muons
      end subroutine decay_info
    end interface

    interface
       subroutine bad_resistor_voltages(time, q1Lbottom, q1Ltop)
         use precision_def
         use sim_utils
         implicit none
         real(rp) time  !seconds
         real(rp) q1Lbottom, q1ltop
       end subroutine bad_resistor_voltages
    end interface

    interface
      subroutine kick_to_pulse_time (kick_time, pulse_time, delay, j)
        use precision_def
        implicit none
        real(rp) kick_time, delay
        real(rp), allocatable :: pulse_time(:)
        integer j,i
       end subroutine
    end interface

    interface
      subroutine combine_kick_and_pulse(kick_and_pulse, pulse_trace, kick_trace, pulse_delay, kick_delay, pulse_points, kick_points)

        use precision_def
        use sim_utils
        use muon_mod
        implicit none
        character*290 new_string,  pulse_trace, kick_trace
        type (kick_and_pulse_struct), allocatable :: kick_and_pulse(:)
        real(rp) delay, pulse_delay, kick_delay
        integer kick_points, pulse_points        
      end subroutine
    end interface

    interface
      subroutine write_amp_phase_momentum(npoints,ix_slice,index, eta, beta, revolution, b_2, b_3, pulse,kick_time,time_step_kick, slice_ave, width, xinf, bxpinf, dxpinf, n_momenta, all_done)

        use precision_def
        use sim_utils
        implicit none
        real(rp) index, eta, slice_ave, width, xinf,  bxpinf, dxpinf, beta,  pulse,time_step_kick, kick_time
        real(rp) revolution, omega_c, xp0
        real(rp) b_2, b_3 ! sextupole and quadrupole terms in units of the field index
        real(rp) Qx, Qx_0
        real(rp) momentum_step
        integer n_momenta, i, ix_slice,npoints
        logical all_done
     end subroutine
   end interface

  interface
   subroutine get_momentum_range(kick_tot, kick_to_gauss, B_to_kick, eta, beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, &
        Aperture, Bmin, Bmax, sigma_xp,delta_p, dp_max, dp_min, angle_correction, &
        delta_central, Nxp,B,kick_tot_sin_to_angle_beta, bxpinf, xpinf2_avg, xp_limit)    
    use precision_def
    implicit none
    real(rp) kick_tot, Aperture, kick_to_gauss, B_to_kick, beta, kick_tot_cos, kick_tot_sin, xinf, xpinf, Bmin, Bmax, sigma_xp,B
    real(rp) delta_p, dp_max, dp_min, angle_correction,kick_tot_sin_to_angle_beta
    real(rp) Nxp, delta_central
    real(rp) bxpinf, eta, xpinf2_avg, xp_limit
   end subroutine
  end interface

  interface
   subroutine opt_inf_angle(lat, co,nbranch, track_state, err_flag, kicker_params)
    use bmad
    use muon_mod

   implicit none

   type (lat_struct), target :: lat
   type (coord_struct), allocatable, target :: co(:)
   type (kicker_params_struct) kicker_params, kicker_params_0
   integer nbranch, track_state, nturns
   logical err_flag
  end subroutine

  end interface

  interface
       subroutine injection_channel(lat,spin_tracking_on, from_orbit, circumference, make_movie)
         use bmad
         use parameters_bmad
         implicit none
         type (lat_struct), target :: lat
         type (ele_struct) ele_at_s, ele_offset
         type (coord_struct), allocatable :: from_orbit(:), to_orbit(:)
         type (coord_struct) orb_at_s, co
         type (em_field_struct) field
         type (g2twiss_struct) twiss, twiss2
         real(rp) circumference
         logical spin_tracking_on, make_movie
       end subroutine injection_channel
    end interface

    interface
       subroutine bin_spin_angle_by_time( turn, unit, nmuons,muons, dphi_dp, dphi_dp_err)
         use muon_mod
         use parameters_bmad
         implicit none
         type (muon_struct), allocatable :: muons(:)

         real(rp), save :: p_min, p_max, p_over
         real(rp) pvec(3), svec(3),sdotp,sdotp_norm, spin_phase

         real(rp), save, allocatable :: spin_phase_sum(:,:),spin_phase_prod(:,:)
         real(rp), allocatable :: avg_phase(:), p(:), rms(:)
         real(rp), optional,  allocatable :: dphi_dp(:),dphi_dp_err(:)
         real(rp) siga, sigb, a,b, chi2,q
         real(rp) turn

         integer unit, nmuons
       end subroutine bin_spin_angle_by_time
    end interface

    interface
       subroutine time_bins(co, n_ele_track, nmuons, time_binning_start, ele)
         use bmad
         use parameters_bmad
         implicit none
         type (coord_struct), allocatable :: co(:)
         type (ele_struct), allocatable :: ele(:)
         integer bin,j,n_ele_track
         integer nmuons
         real(rp) time_binning_start
       end subroutine time_bins
    end interface

    interface
      subroutine manage_polarization(nmuons,muons, gamma)
        use muon_mod
        use parameters_bmad
        implicit none
        type(muon_struct), allocatable :: muons(:)
        integer nmuons
        real(rp) gamma
      end subroutine manage_polarization
    end interface

    interface
       subroutine RightOrLeft(pvec,svec,sxp_pp)
         use sim_utils
         implicit none
         real(rp) svec(3), pvec(3)
         real(rp) sxp_p(3), sxp_pp(3)
       end subroutine RightOrLeft
     end interface

     interface
        subroutine write_phase_space_at_time(co, spin_phase,write_time, bin, write_bins, close_files)
          use bmad
          use bmad_struct
          use bmad_interface
          use parameters_bmad
          implicit none
          type (coord_struct) co
          real(rp) spin_phase
          real(rp) write_time
          integer, allocatable :: write_bins(:)
          integer bin
          logical close_files
        end subroutine write_phase_space_at_time
     end interface

     interface
        subroutine get_time_dep_moments(n,moments, running_moment)
          use precision_def
          use muon_mod
!          use muon_interface, dummy => get_time_dep_moments
          implicit none
          integer n,i,j,idum
          TYPE (g2moment_struct), allocatable :: moments(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width
          type (running_moment_struct), allocatable :: running_moment(:)

        end subroutine get_time_dep_moments
     end interface

     interface
        subroutine fit_average_x(n,moments, running_moment, coef, chisq, ndf)
          use precision_def
          use muon_mod
!          use muon_interface, dummy => get_time_dep_moments
          use parameters_bmad
          implicit none
  
          TYPE (g2moment_struct), allocatable :: moments(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width
          type (running_moment_struct), allocatable :: running_moment(:)

          integer n,i,j,idum
          real(rp) coef(4)
          real(rp) chisq
          integer ndf
          
        end subroutine fit_average_x
     end interface

     interface
        subroutine time_dep_quad_field (start_orb, ele, quad_number, err_flag)

          use bmad !, except_dummy => track1_preprocess
          use quad_scrape_parameters
          use parameters_bmad
          use muon_mod
          implicit none

          type (coord_struct) :: start_orb, orb
          type (ele_struct), target :: ele
          type (ele_struct), pointer :: ele_init
          type (ele_struct), pointer :: lord

          logical err_flag, err
          integer quad_number          
        end subroutine time_dep_quad_field
     end interface

     interface
         subroutine get_xavg( n,moments, running_moment, fit_number)
           use precision_def
           use muon_mod
           implicit none
           
           TYPE (g2moment_struct), allocatable :: moments(:) !moment structure for each time bin, 1000 turns => an array with 1000*149e-9/bin_width
           type (running_moment_struct), allocatable :: running_moment(:)

           integer n,i,j,idum
           integer lun
           integer fit_number
         end subroutine get_xavg
      end interface

   end module muon_interface