!+
! Subroutine radiation_integrals_custom (lat, ir, orb, rad_int1, err_flag)
!
! Dummy routine for custom elements. Will generate an error if called.
! A valid radiation_integrals_custom is needed only if the 
! radiation_integrals routine is being used.
!
! General rule: Your code may NOT modify any argument that is not listed as
! an output agument below.
!
! Modules needed:
!   use rad_int_common
!
! Input:
!   lat         -- lat_struct: Lattice with the custom element.
!   ir          -- integer: lat%ele(ir) is the custom element.
!   orb(:)      -- coord_struct: Orbit around which integrals are to be evaluated.
!
! Output:
!   rad_int1    -- rad_int1_struct: Structure for storing the results.
!   err_flag    -- logical: Set true if there is an error. False otherwise.
!-

subroutine radiation_integrals_custom (lat, ir, orb, rad_int1, err_flag)

use bmad_interface, dummy => radiation_integrals_custom

implicit none

type (lat_struct) lat
type (coord_struct) orb(0:)
type (rad_int1_struct) rad_int1
integer ir
logical err_flag
character(32) :: r_name = 'radiation_integrals_custom'

!These are my additions
type (ele_struct), pointer :: ele
real(rp) L0, K_und, B, vbar, gamma_rel
real(rp) lambda, rho, curly_H_a, curly_H_b, r0
real(rp) p0

!get info on element in question
ele => lat%ele(ir)

p0 = orb(ir-1)%p0c/c_light*e_charge !momentum of electron in SI units
L0 = ele%value(n_period$)*ele%value(l_period$) !undulator length
lambda = ele%value(l_period$) !undulator wavelength
B = ele%value(b_max$) !peak field (T)
gamma_rel = sqrt((orb(ir-1)%p0c*(1.+orb(ir-1)%vec(6)))**2+m_electron**2)/m_electron !relativistic gamma

K_und = e_charge*B*lambda/(twopi*m_electron*e_charge/c_light) !undualtor paramater
vbar = c_light*(1.-1./2./gamma_rel/gamma_rel*(1.+K_und*K_und)) !average longitudinal velocity (m/s)

rho = p0/(e_charge*B) !radius of curavature

r0 = e_charge*B*c_light**2/(vbar*gamma_rel*m_electron*e_charge)*(lambda/twopi)**2

curly_H_a = ele%a%gamma*ele%a%eta**2 + 2.*ele%a%alpha*ele%a%eta*ele%a%etap + ele%a%beta*ele%a%etap**2 !dispersion invariants in a,b modes
curly_H_b = ele%b%gamma*ele%b%eta**2 + 2.*ele%b%alpha*ele%b%eta*ele%b%etap + ele%b%beta*ele%b%etap**2

!these are actually important
rad_int1%i0 = L0*gamma_rel/rho
rad_int1%i1 = -L0*r0/rho
rad_int1%i2 = L0/rho**2
rad_int1%i3 = L0/rho**3
rad_int1%i4a = L0*(-1./rho**2*r0/rho/2. - 25./32.*(lambda/twopi)**2*(e_charge*B/p0)**4)
rad_int1%i4b = L0*(-1./rho**2*r0/rho/2. - 25./32.*(lambda/twopi)**2*(e_charge*B/p0)**4)
rad_int1%i4z = rad_int1%i4a + rad_int1%i4b
rad_int1%i5a = L0*curly_H_a/rho**3
rad_int1%i5b = L0*curly_H_b/rho**3
rad_int1%i6b = L0*ele%b%beta/rho**3

!unsure about whether we really need the below, but can include them anyway...
rad_int1%lin_i2_E4 = L0*gamma_rel**4/rho**2
rad_int1%lin_i3_E7 = L0*gamma_rel**7/rho**3
rad_int1%lin_i5a_E6 = L0*gamma_rel**6/rho**3*curly_H_a
rad_int1%lin_i5b_E6 = L0*gamma_rel**6/rho**3*curly_H_b

!unsure about how to handle the below
rad_int1%lin_norm_emit_a = 0 !running sum from start of branch
rad_int1%lin_norm_emit_b = 0
rad_int1%n_steps = 0 !for use in qromb


err_flag = .false.


end subroutine