program field_emitter

use field_emitter_mod
use beam_mod
use wall3d_mod

implicit none

type (lat_struct), target :: lat
type (branch_struct), pointer :: branch
type (ele_struct), pointer ::  eles(:)
type (ele_struct) :: ele
type (coord_struct) :: orb0
type (coord_struct), allocatable :: closed_orb(:), particles(:),&
all_particles(:)
real(rp) ::  beta_fn, A_fn, work_f,q, t, t_a, t_b, perp(3), angle, s_init
type(field_emitter_struct), allocatable :: FE_list(:)
integer ::  i_dim, n_t, n, j,  p, n_contour_points, contour_file_num, o_stat
integer :: tot_particles, upto_particles, n_particles, namelist_file, n_char, open_status, unit_spec, unit_num
integer:: ix_FE, ix_ele, wall_contour_num
character(40) :: x1, spec_file, FE_num_fmt, filename, contour_filename
character(100) :: lat_name, lat_path, base_name, in_file
character(30), parameter :: r_name = 'field_emitter'
logical :: separate_files, write_info
namelist / field_emitter_params / &
    lat_name, spec_file, separate_files, n_contour_points, write_info

!------------------------------------------

print *, ' '
print '(a85)',' ___    ___       __       ___         ___ ___  ___  __  '
print '(a85)','|__  | |__  |    |  \     |__   |\/| |  |   |  |__  |__) '
print '(a85)','|    | |___ |___ |__/ ___ |___  |  | |  |   |  |___ |  \ '
print *, ' '
print *, ' '


 !Defaults for namelist
lat_name = 'lat.bmad'
spec_file = 'field_emitters.dat'
n_contour_points = 0
write_info = .false.

!Read namelist
in_file = 'field_emitter.in'
if (command_argument_count() > 0) call get_command_argument(1, in_file)

namelist_file = lunget()
print *, 'Opening: ', trim(in_file)
print *, ' '
open (namelist_file, file = in_file, status = "old",  iostat=open_status)
if (open_status /= 0) then
  print *, 'File missing: ', in_file
  print *, 'Stopping'
  stop
endif
read (namelist_file, nml = field_emitter_params)
close (namelist_file)

if (separate_files) then
  print*, 'Separate_files = true'
  print *, 'All field emitters will be written into different files.'
else
  print *, 'Separate_files = false'
  print *, 'All field emitters will be written into the same file.'
endif

print *, ' '

!Trim filename
n_char= SplitFileName(lat_name, lat_path, base_name) 

print *, 'Parsing ', trim(adjustl(lat_name)), '..'
!Parse Lattice
call bmad_parser (lat_name, lat)
!branch => lat%branch(0)

! Check for absolute time tracking. If not, abort!
if (.not. bmad_com%absolute_time_tracking ) then
  call out_io (s_error$, r_name, 'absolute time tracking must be set to True.')
  stop
endif

!Get number of field emitters and allocate FE_list accordingly
print *, 'Opening: ', trim(adjustl(spec_file))
print *, ' '
print *, 'Finding total number of field emitters..'
unit_spec = lunget()

o_stat = 0
open(unit_spec, file = spec_file, iostat = o_stat)

if ( o_stat /= 0 ) then
  print *, 'Data file ', spec_file, ' not found!'
  print *, 'Stopping.'
  stop
endif

read(unit_spec, '(I8)') n
print *, 'Field emitters found: ', n
allocate(FE_list(n))
!Leave spec_file open so read_FE picks up at line 2

!Populate the coord_struct array FE%particles with proper x, y, s_rel, t,
!and charge values and loop over all field emitters to FE_list.
print *, 'Calculating magnitude of charge for each field emitter..'
do ix_FE=1, ubound(FE_list, 1)
  call read_FE(unit_spec, lat, FE_list(ix_FE))
enddo
close(unit_spec)
!separate_file = .true./.false. is an option in the namelist that, if set to true, writes the!data of each particle to a separate file in the current directory of the form
!field_emitter[field emitter number].out. If separate_files is set to false, the data is written into field_emitters.out


if (separate_files) then
  do ix_FE=1, ubound(FE_list, 1)
    if (.not. allocated(FE_list(ix_FE)%particles)) cycle
    unit_num = lunget()
    FE_num_fmt = '(I8)' 
    write(x1, FE_num_fmt) ix_FE 
    filename = 'field_emitter'//trim(adjustl(x1))//'.dat'
    open(unit_num, file = filename)
    write(unit_num, '(I8)') size(FE_list(ix_FE)%particles)
    call particle_writer(FE_list(ix_FE)%particles, unit_num)
    close(unit_num)
    print *, 'Particle data written: ', filename
    if (.not. write_info) cycle
    filename = 'field_emitter'//trim(adjustl(x1))//'.info'
    unit_num = lunget()
    open(unit_num, file = filename)
    call write_fe_info(FE_list(ix_FE), unit_num)
    close(unit_num)
    print *, 'Field emitter info written: ', filename
  enddo
else
  upto_particles = 0
  do ix_FE=1, ubound(FE_list, 1)
    n_particles = particle_counter(FE_list(ix_FE)%particles)
    tot_particles = n_particles + upto_particles
    upto_particles = tot_particles
  enddo
  
  unit_num = lunget()
  filename = 'field_emitter.dat'
  open(unit_num, file = filename)
  write(unit_num, '(I8)') tot_particles
  do ix_FE=1, ubound(FE_list, 1)
    call particle_writer(FE_list(ix_FE)%particles, unit_num)
  enddo
  print *, 'Particle data written: ', filename
  print *, 'Total particles=', tot_particles
close(unit_num)
endif

print *, ' '
print *, ' '



if (n_contour_points == 0) stop

print *, 'n_contour_points is nonzero.' 
print *, 'Running wall_contour..'

angle = 0
s_init=0
contour_file_num = lunget()
print *, 'Contour points being found at angle ', angle
contour_filename = 'wall_points.contour'
open(contour_file_num, file = contour_filename) 
write(contour_file_num, '(4A25,$)') 'x', 'y', 's', 'E_max','ele_name'
do p=1, ubound(lat%ele, 1)
  ele = lat%ele(p)
  if (ele%key /= lcavity$ .and. ele%key /= rfcavity$) then
    cycle
  else
  print *,'Finding contour points for ', trim(adjustl(ele%name))
!  s_init = ele%s - ele%value(l$)
  call wall_contour(lat%ele(p), angle, contour_file_num, n_contour_points, s_init)
  endif
enddo
print *, 'Contour points found.'
print *, 'Points written in: ', contour_filename
close(contour_file_num)


end program