program resonantextraction

   use bmad
    use beam_mod
     implicit none

      type (lat_struct), target :: lat
      type (branch_struct), pointer :: to_branch
      type (ele_struct), pointer :: rex, esrefpt, magrefpt, test, mq1, mq2, mq3, mq4, to_ele
      type (ele_struct), pointer :: einj, pinj, eext1, eext2, pext1, pext2, bump1, bump2, bump3
      type (ele_pointer_struct), allocatable :: eles(:), elesa(:), elesb(:), elesc(:), elesd(:), elese(:)
      type (ele_pointer_struct), allocatable :: elesf(:), elesg(:), elesh(:), elesi(:), elesj(:), elesk(:)
      type (ele_pointer_struct), allocatable :: elesl(:), elesm(:), elesn(:)
      type (beam_init_struct) init
      type (beam_struct) beam
      type (coord_struct), allocatable :: outgoing(:,:),trajectory(:)

      integer, allocatable :: transferturn(:)
      integer ix, nargs
      integer i, j, k, m, n, p, q, r, m1, n1, p1, q1, r1
      integer lun
      integer nturns, turnstep, startpoint
      integer n_loc
      integer ib_to, ie_to
      integer bumpinitturn, bumpramplength, rexinitturn, rexramplength

      real rexinit, rexmax, kqhinit, kqh1, kqh2, kqh3, kqh4
      real initialwait, ramp1, ramp2, ramp3, ramp4, ramp5, finalwait
      real esseptinner, esseptouter, magseptinner, magseptouter, esseptkick, esouterwall, magouterwall
      real magseptx, magseptpx, xadj, pxadj
      real kick1, kick2, kick3, initkick1, initkick2, initkick3
      real turnreal

      character*140 lat_file
      character*120 line, lineout
      character*10 iteration

      logical err

      bmad_com%auto_bookkeeper = .false.

  !ready an input file from the command line
      nargs = cesr_iargc()
      if(nargs == 1)then
        call cesr_getarg(1,lat_file)
        print *, 'Using ', trim(lat_file)
       else
        lat_file = 'bmad.'
        print '(a,$)',' Lattice file name ? (default= synch.lat) '
        read(5,'(a)') line
        call string_trim(line, lineout, ix)
        lat_file = lineout
        if(ix == 0) lat_file = '/home/jdp279/Documents/Track6610/Outgoing/darkphoton/files/synch.lat'
        print *, ' lat_file = ', lat_file
      endif

  !run parameters
      turnstep = 40       ! turns between beam profile outputs
      nturns = 4080       ! should be a multiple of turnstep
      startpoint = 2440   ! should be a multiple of turnstep and less than initialwait. Use to start run in the middle.
                          ! For now, depending on what parts are skipped, need to manually change the seven !** values below.

      initialwait = 2525
      ramp1 = 250
      ramp2 = 100
      ramp3 = 750
      ramp4 = 150
      ramp5 = 300
      finalwait = 5      ! summing these should equal nturns

      kqhinit = 0.40096   ! initial main quadrupole strength
      kqh1 = 0.4470       ! after ramp 1
      kqh2 = 0.4495       ! after ramp 2
      kqh3 = 0.4575       ! after ramp 3
      kqh4 = 0.4575       ! after ramp 4

      rexinit = 5.0       !** 0
      rexinitturn = 500
      rexramplength = 1000
      rexmax = 5.0

  !beam parameters
      init%n_particle=10000 ! number of particles to track
      init%a_emit=1.68E-7   !** 6e-6
      init%b_emit=2.8E-8    !** 1e-6
      init%sig_e=3.3E-4     !** 2e-3
      init%n_bunch=1

  !septum parameters
      esseptinner=0.0310
      esseptouter=0.0311
      esouterwall=0.040
      esseptkick=0.00017

      magseptinner=0.0247
      magseptouter=0.0257
      magouterwall=0.040
      magseptx=0.0271
      magseptpx=0.00303
      xadj=-0.20
      pxadj=-0.037

  !bump parameters

      initkick1=0.00121      !** 0
      initkick2=0.00057      !** 0
      initkick3=-0.00178      !** 0
      bumpinitturn=500
      bumpramplength=1000
      kick1= 0.00121
      kick2= 0.00057
      kick3= -0.00178

  !set up lattice parameters

      call bmad_parser (lat_file, lat)  ! parse lattice file

    !shortcut references for lattice elements.
      !extraction elements
      call lat_ele_locator('REX',lat,elesa,n_loc,err)
      rex => elesa(1)%ele

      call lat_ele_locator('ESSEPT',lat,elesb,n_loc,err)
      esrefpt => elesb(1)%ele

      call lat_ele_locator('MAGSEPT',lat,elesc,n_loc,err)
      magrefpt => elesc(1)%ele

      call lat_ele_locator('Q01V',lat,elesd,n_loc,err)
      mq1 => elesd(1)%ele
      mq2 => elesd(2)%ele

      call lat_ele_locator('Q02H',lat,elese,n_loc,err)
      mq3 => elese(1)%ele
      mq4 => elese(2)%ele

      !preexisting septa
      call lat_ele_locator('ELINJSEP',lat,elesf,n_loc,err)
      einj => elesf(1)%ele

      call lat_ele_locator('POSINJSEP',lat,elesg,n_loc,err)
      pinj => elesg(1)%ele

      call lat_ele_locator('ELEXT1SEP',lat,elesh,n_loc,err)
      eext1 => elesh(1)%ele

      call lat_ele_locator('ELEXT2SEP',lat,elesi,n_loc,err)
      eext2 => elesi(1)%ele

      call lat_ele_locator('POSEXT1SEP',lat,elesj,n_loc,err)
      pext1 => elesj(1)%ele

      call lat_ele_locator('POSEXT2SEP',lat,elesk,n_loc,err)
      pext2 => elesk(1)%ele

      !bump kickers
      call lat_ele_locator('BUMP1',lat,elesl,n_loc,err)
      bump1 => elesl(1)%ele

      call lat_ele_locator('BUMP2',lat,elesm,n_loc,err)
      bump2 => elesm(1)%ele

      call lat_ele_locator('BUMP3',lat,elesn,n_loc,err)
      bump3 => elesn(1)%ele


      !set initial element strengths
      mq1%value(k1$)=-kqhinit
      mq2%value(k1$)=-kqhinit
      mq3%value(k1$)=kqhinit
      mq4%value(k1$)=kqhinit
      call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
      call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
      call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
      call set_flags_for_changed_attribute(mq4,mq4%value(k1$))

      rex%a_pole(2)=rexinit
      call set_flags_for_changed_attribute(rex,rex%a_pole(2))

      bump1%value(kick$)=initkick1
      bump2%value(kick$)=initkick2
      bump3%value(kick$)=initkick3
      call set_flags_for_changed_attribute(bump1,bump1%value(kick$))
      call set_flags_for_changed_attribute(bump2,bump2%value(kick$))
      call set_flags_for_changed_attribute(bump3,bump3%value(kick$))

      call lattice_bookkeeper(lat)
      call lat_make_mat6(lat,-1)        ! compute transfer matrices for all of the elements in the lattice

      call twiss_at_start(lat)          ! compute twiss parameters at the starting point
      call twiss_propagate_all(lat)     ! propagate twiss parameters to all elements in the lattice
      call writeinfo(lat)               ! write initial twiss parameters to file

      lun = lunget()
      open(unit=lun, file='placements.dat')
      write(lun,'(a15,i4)')'rex = ', rex%ix_ele
      write(lun,'(a15,i4)')'esrefpt = ', esrefpt%ix_ele
      write(lun,'(a15,i4)')'magrefpt = ', magrefpt%ix_ele
      write(lun,'(a15,i4)')'einj = ', einj%ix_ele
      write(lun,'(a15,i4)')'pinj = ', pinj%ix_ele
      write(lun,'(a15,i4)')'eext1 = ', eext1%ix_ele
      write(lun,'(a15,i4)')'eext2 = ', eext2%ix_ele
      write(lun,'(a15,i4)')'pext1 = ', pext1%ix_ele
      write(lun,'(a15,i4)')'pext2 = ', pext2%ix_ele
      write(lun,'(a15,i4)')'bump1 = ', bump1%ix_ele
      write(lun,'(a15,i4)')'bump2 = ', bump2%ix_ele
      write(lun,'(a15,i4)')'bump3 = ', bump3%ix_ele
      close(unit=lun)


  !locate forking structures

      ib_to=nint(magrefpt%value(ix_to_branch$))
      ie_to=nint(magrefpt%value(ix_to_element$))

      to_branch => lat%branch(ib_to)
      to_ele => to_branch%ele(ie_to)

      allocate(outgoing(init%n_particle,0:to_branch%n_ele_track))
      allocate(trajectory(0:to_branch%n_ele_track))
      allocate(transferturn(init%n_particle))  !this index is used somewhat generally to track the status of particles.
      do j=1,init%n_particle
        transferturn(j)=-1
      end do

  !track beam through lattice. (The proper functioning of this section depends on the order of elements in the lattice file!)

      call init_beam_distribution(lat%ele(0),lat%param,init,beam)
      do i=(startpoint/turnstep), (nturns/turnstep)-1
        do k=0, turnstep-1
          print *, 'turn', (i*turnstep)+k, 'quad = ', mq4%value(k1$), 'bump = ', bump1%value(kick$), 'rex = ', rex%a_pole(2)

          call track_beam(lat,beam,lat%ele(0),rex,err)

          !checks for 'aperture violations' over the last full turn, ignoring extracted particles.
          m=0
          do j=1, init%n_particle
            if(transferturn(j)==-1)then
              if(beam%bunch(1)%particle(j)%state == alive$)then
               else
                transferturn(j)=(i*turnstep)+k+(nturns*4)
                m=m+1
              end if
            end if
          end do
          if(m>1)then
            print *, m, 'particles hit beam pipe!'
           else if(m==1)then
            print *, m, 'particle hit beam pipe!'
          end if

          !writes beam profile at sextupole
          if(k==0)then
            write(iteration,'(I4.4)')(i*turnstep)+k
            lun = lunget()
            open(unit=lun, file='rexturn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          call track_beam(lat,beam,rex,bump1,err)

          !writes beam profile at bump 1
          if(k==0)then
            write(iteration,'(I4.4)')(i*turnstep)+k
            lun = lunget()
            open(unit=lun, file='bump1turn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          call track_beam(lat,beam,bump1,esrefpt,err)

          !electrostatic septum checks and kick
          m=0
          n=0
          p=0
          do j=1, init%n_particle
            if(transferturn(j)==-1)then
              if(beam%bunch(1)%particle(j)%vec(1) > esseptinner)then
                if(beam%bunch(1)%particle(j)%vec(1) > esseptouter)then
                  if(beam%bunch(1)%particle(j)%vec(1) > esouterwall)then
                    transferturn(j) = (i*turnstep)+k+(nturns*2)
                    m=m+1
                   else                
                    beam%bunch(1)%particle(j)%vec(2) = beam%bunch(1)%particle(j)%vec(2)+esseptkick
                    n=n+1
                  end if
                 else
                  transferturn(j) = (i*turnstep)+k+nturns
                  p=p+1
                end if
              end if
            end if
          end do
          if(m>1)then
            print *, m, 'particles hit outer wall (essept)!'
           else if(m==1)then
            print *, m, 'particle hit outer wall (essept)!'
          end if
          if(n>1)then
            print *, n, 'particles kicked'
           else if(n==1)then
            print *, n, 'particle kicked'
          end if
          if(p>1)then
            print *, p, 'particles hit electrostatic septum!'
           else if(p==1)then
            print *, p, 'particle hit electrostatic septum!'
          end if

          !writes beam profile at electrostatic septum
          if(k==0)then
            write(iteration,'(I4.4)')(i*turnstep)+k
            lun = lunget()
            open(unit=lun, file='esturn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          call track_beam(lat,beam,esrefpt,magrefpt,err)

          !writes beam profile at magnetic septum
          if(k==0)then
            lun = lunget()
            open(unit=lun, file='magturn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          !magnetic septum checks and extraction
          m=0
          n=0
          p=0
          do j=1, init%n_particle
            if(transferturn(j)==-1)then
              if(beam%bunch(1)%particle(j)%vec(1) > magseptinner)then
                if(beam%bunch(1)%particle(j)%vec(1) > magseptouter)then
                  if(beam%bunch(1)%particle(j)%vec(1) > magouterwall)then
                    transferturn(j) = (i*turnstep)+k+(nturns*2)
                    m=m+1
                   else                
                    outgoing(j,ie_to)=beam%bunch(1)%particle(j)
                    outgoing(j,ie_to)%vec(1)=outgoing(j,ie_to)%vec(1)-(magseptx+xadj*(mq3%value(k1$)-(kqh4+kqh1)/2))
                    outgoing(j,ie_to)%vec(2)=outgoing(j,ie_to)%vec(2)-(magseptpx+pxadj*(mq3%value(k1$)-(kqh4+kqh1)/2))
                    transferturn(j) = (i*turnstep)+k
                    n=n+1
                  end if
                 else
                  transferturn(j)=(i*turnstep)+k+(nturns*3)
                  p=p+1
                end if
              end if
            end if
          end do
          if(m>1)then
            print *, m, 'particles hit outer wall (magsept)!'
           else if(m==1)then
            print *, m, 'particle hit outer wall (magsept)!'
          end if
          if(n>1)then
            print *, n, 'particles transferred'
           else if(n==1)then
            print *, n, 'particle transferred'
          end if
          if(p>1)then
            print *, p, 'particles hit magnetic septum!'
           else if(p==1)then
            print *, p, 'particle hit magnetic septum!'
          end if

          call track_beam(lat,beam,magrefpt,bump2,err)

          !writes beam profile at bump 2
          if(k==0)then
            write(iteration,'(I4.4)')(i*turnstep)+k
            lun = lunget()
            open(unit=lun, file='bump2turn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          call track_beam(lat,beam,bump2,bump3,err)

          !writes beam profile at bump 3
          if(k==0)then
            write(iteration,'(I4.4)')(i*turnstep)+k
            lun = lunget()
            open(unit=lun, file='bump3turn'//trim(iteration)//'.dat')
            write(lun,'(a10,6a12)')'particle','x','px','y','py','z','delta'
            do j=1,init%n_particle
              if(transferturn(j)==-1)then
                write(lun,'(i10,6es12.4)')j, beam%bunch(1)%particle(j)%vec(1:6)
              end if
            end do
            close(unit=lun)
          end if

          call track_beam(lat,beam,bump3,lat%ele(lat%n_ele_track),err)

          !set element strengths for next turn, write twiss parameters during extraction, and write chromaticities.

          if (((i*turnstep)+k) > (bumpinitturn-1) .AND. ((i*turnstep)+k) < (bumpinitturn+bumpramplength)) then
            bump1%value(kick$)=bump1%value(kick$)+ kick1/bumpramplength
            bump2%value(kick$)=bump2%value(kick$)+ kick2/bumpramplength
            bump3%value(kick$)=bump3%value(kick$)+ kick3/bumpramplength
            call set_flags_for_changed_attribute(bump1,bump1%value(kick$))
            call set_flags_for_changed_attribute(bump2,bump2%value(kick$))
            call set_flags_for_changed_attribute(bump3,bump3%value(kick$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > (rexinitturn-1) .AND. ((i*turnstep)+k) < (rexinitturn+rexramplength)) then
            rex%a_pole(2)=rex%a_pole(2)+ rexmax/rexramplength
            call set_flags_for_changed_attribute(rex,rex%a_pole(2))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > initialwait-1 .AND. ((i*turnstep)+k) < (initialwait+ramp1)) then
            mq1%value(k1$)=mq1%value(k1$) -((kqh1-kqhinit)/ramp1)
            mq2%value(k1$)=mq2%value(k1$) -((kqh1-kqhinit)/ramp1)
            mq3%value(k1$)=mq3%value(k1$) +((kqh1-kqhinit)/ramp1)
            mq4%value(k1$)=mq4%value(k1$) +((kqh1-kqhinit)/ramp1)
            call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
            call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
            call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
            call set_flags_for_changed_attribute(mq4,mq4%value(k1$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > (initialwait+ramp1-1) .AND. ((i*turnstep)+k) < (initialwait+ramp1+ramp2)) then
            mq1%value(k1$)=mq1%value(k1$) -((kqh2-kqh1)/ramp2)
            mq2%value(k1$)=mq2%value(k1$) -((kqh2-kqh1)/ramp2)
            mq3%value(k1$)=mq3%value(k1$) +((kqh2-kqh1)/ramp2)
            mq4%value(k1$)=mq4%value(k1$) +((kqh2-kqh1)/ramp2)
            call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
            call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
            call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
            call set_flags_for_changed_attribute(mq4,mq4%value(k1$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > (initialwait+ramp1+ramp2-1) .AND. ((i*turnstep)+k) < (initialwait+ramp1+ramp2+ramp3)) then
            mq1%value(k1$)=mq1%value(k1$) -((kqh3-kqh2)/ramp3)
            mq2%value(k1$)=mq2%value(k1$) -((kqh3-kqh2)/ramp3)
            mq3%value(k1$)=mq3%value(k1$) +((kqh3-kqh2)/ramp3)
            mq4%value(k1$)=mq4%value(k1$) +((kqh3-kqh2)/ramp3)
            call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
            call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
            call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
            call set_flags_for_changed_attribute(mq4,mq4%value(k1$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > (initialwait+ramp1+ramp2+ramp3-1) .AND.&
          ((i*turnstep)+k) < (initialwait+ramp1+ramp2+ramp3+ramp4)) then
            mq1%value(k1$)=mq1%value(k1$) -((kqh4-kqh3)/ramp4)
            mq2%value(k1$)=mq2%value(k1$) -((kqh4-kqh3)/ramp4)
            mq3%value(k1$)=mq3%value(k1$) +((kqh4-kqh3)/ramp4)
            mq4%value(k1$)=mq4%value(k1$) +((kqh4-kqh3)/ramp4)
            call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
            call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
            call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
            call set_flags_for_changed_attribute(mq4,mq4%value(k1$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) > (initialwait+ramp1+ramp2+ramp3+ramp4-1) .AND.&
         ((i*turnstep)+k) < (initialwait+ramp1+ramp2+ramp3+ramp4+ramp5)) then
            mq1%value(k1$)=mq1%value(k1$) -((kqhinit-kqh4)/ramp5)
            mq2%value(k1$)=mq2%value(k1$) -((kqhinit-kqh4)/ramp5)
            mq3%value(k1$)=mq3%value(k1$) +((kqhinit-kqh4)/ramp5)
            mq4%value(k1$)=mq4%value(k1$) +((kqhinit-kqh4)/ramp5)
            call set_flags_for_changed_attribute(mq1,mq1%value(k1$))
            call set_flags_for_changed_attribute(mq2,mq2%value(k1$))
            call set_flags_for_changed_attribute(mq3,mq3%value(k1$))
            call set_flags_for_changed_attribute(mq4,mq4%value(k1$))
            call lat_make_mat6(lat,-1)
            call lattice_bookkeeper(lat)
          end if

          if (((i*turnstep)+k) == 3610) then
            call twiss_at_start(lat)
            call twiss_propagate_all(lat)
            call writeinfo2(lat)
          end if

          !adiabatic damping from energy change
          turnreal = REAL((i*turnstep)+k)
          do j=1, init%n_particle
            beam%bunch(1)%particle(j)%vec(1) = &
            beam%bunch(1)%particle(j)%vec(1)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
            beam%bunch(1)%particle(j)%vec(2) = &
            beam%bunch(1)%particle(j)%vec(2)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
            beam%bunch(1)%particle(j)%vec(3) = &
            beam%bunch(1)%particle(j)%vec(3)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
            beam%bunch(1)%particle(j)%vec(4) = &
            beam%bunch(1)%particle(j)%vec(4)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
            beam%bunch(1)%particle(j)%vec(5) = &
            beam%bunch(1)%particle(j)%vec(5)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
            beam%bunch(1)%particle(j)%vec(6) = &
            beam%bunch(1)%particle(j)%vec(6)*(1+((sin(turnreal/1052.0))/(2104.0*(cos(turnreal/1052.0)-1.0)-103.5)))
          end do

        end do
      end do

  !track outgoing particles to target

      lun = lunget()
      open(unit=lun, file='outgoing.dat')
      write(lun,'(a10,7a12)')'particle','transferturn','x','px','y','py','z','delta'
      do j=1,init%n_particle
        if(-1<transferturn(j) .AND. transferturn(j)<nturns)then
          trajectory=outgoing(j,0:to_branch%n_ele_track)
          call track_all(lat,trajectory,ib_to)
          write(lun,'(2i10,6es12.4)')j, transferturn(j), trajectory(to_branch%n_ele_track)%vec(1:6)
        end if
      end do
      close(unit=lun)

  !write run statistics to file

      lun = lunget()
      open(unit=lun, file='runstats.dat')
      write(lun,'(a10,5a15)')'turn','transferred','hitesseptum','hitmagseptum','hitouterwall','hitbeampipe'
      m1=0
      n1=0
      p1=0
      q1=0
      r1=0
      do i=0, nturns-1
        m=0
        n=0
        p=0
        q=0
        r=0
        do j=1, init%n_particle
          if(transferturn(j)==i)then
            m=m+1
           else if(transferturn(j)==i+nturns)then
            n=n+1
           else if(transferturn(j)==i+(nturns*2))then
            p=p+1
           else if(transferturn(j)==i+(nturns*3))then
            q=q+1
           else if(transferturn(j)==i+(nturns*4))then
            r=r+1
          end if
        end do
        write(lun,'(i10,5i15)')i,m,n,q,p,r
        m1=m1+m
        n1=n1+n
        p1=p1+p
        q1=q1+q
        r1=r1+r
      end do
      write(lun,'(a10,5i15)')'total',m1,n1,q1,p1,r1
      close(unit=lun)

 end