! compute electron and positron trajectories, including synchrotron radiation
! bunches collide at L0 and L3

 program speedoflight

   use bmad
   use reverse_mod

      type (lat_struct) lat, lat_rev
      type (ele_struct) ele 
      type (coord_struct), allocatable:: co(:),co_pos(:), co_ele(:), co_pos_a(:)
      type (rad_int_all_ele_struct) rad_int
      type (normal_modes_struct) mode
      integer ix, i, j
      integer nargs,iargc
      integer iu, det_num     
      integer ix_cache
      integer lun
!      integer, parameter :: rf_w1$=1, rf_w2$=2, rf_e1$ =3, rf_e2$ = 4
      character*140 lat_file
      character*120 line
      logical err/.false./
      real(rp) mc2, gamma,Npart/1.e10/,dqh,dqv,sigx,sigy
      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
      lat_file = 'bmad.'
      print '(a37,$)',' Lattice file name ? (default= bmad.) '
      read(5,'(a)') line
       call string_trim(line, line, ix)
       lat_file = line
       if(ix == 0) lat_file = 'bmad.'
       print *, ' lat_file = ', lat_file
      endif

       call bmad_parser (lat_file, lat)
       call lat_reverse(lat, lat_rev)
       call reallocate_coord(co,lat%n_ele_max)
       call reallocate_coord(co_pos,lat%n_ele_max)
       call reallocate_coord(co_pos_a,lat%n_ele_max)
       call reallocate_coord(co_ele, lat_rev%n_ele_max)


       bmad_com%radiation_damping_on = .false.
       bmad_com%radiation_fluctuations_on = .false.

       call closed_orbit_calc(lat,co,i_dim=6,direction = 1, ix_branch=0, err_flag=err)
       call twiss_at_start(lat)
       call twiss_propagate_all(lat)
       ix_cache=0
       call radiation_integrals (lat, co, mode, ix_cache, 0, rad_int)

       mc2 = mass_of(positron$)
       gamma = lat%ele(0)%value(e_tot$)/mc2
       print '(a,es12.4,a,es12.4)',' mc2 = ', mc2, ' gamma = ',gamma
       print '(a,es12.4,a,es12.4,a,es12.4)',' Nparticles = ', Npart,' emit_x = ',mode%a%emittance, ' delta e/e = ',mode%sige_e
       do i=1,lat%n_ele_track
        if(trim(lat%ele(i)%name)=='IP_L0' .or. trim(lat%ele(i)%name)=='IP_L3')then
         sigx = sqrt(mode%a%emittance * lat%ele(i)%a%beta + (lat%ele(i)%x%eta*mode%sige_e)**2)
         sigy = sqrt(mode%a%emittance/10 * lat%ele(i)%b%beta)
         dqh = Npart*classical_radius_factor/twopi/gamma * lat%ele(i)%a%beta/sigx*(sigx+sigy) 
         dqv = Npart*classical_radius_factor/twopi/gamma * lat%ele(i)%b%beta/sigy*(sigx+sigy) 
         print '(a,a,es12.4,a,es12.4)',lat%ele(i)%name,' dqv = ',dqv,' dqh =',dqh
        endif
       end do
       bmad_com%radiation_damping_on = .true.
       bmad_com%radiation_fluctuations_on = .true.

       call closed_orbit_calc(lat,co_pos,i_dim=6,direction = 1, ix_branch=0, err_flag=err)
       call closed_orbit_calc(lat_rev,co_ele,i_dim=6,direction = 1, ix_branch=0, err_flag=err)

      lun=lunget()
      open(unit=lun,file='electron_positron_closed_orbits.dat')
       write(lun,'(36x,a29,1x,a29,1x,a29)')'no synch radiation','positrons','electrons'
       write(lun,'(a16,1x,a10,1x,10a12)')'ele','bpm','s','co x','co y','co delta','pos x','pos y','pos delta','ele x','ele y','ele delta'
       do i=1,lat%n_ele_track
        j = lat%n_ele_track+1-i
        print '(i10,1x,i10,1x,a,a)',i,j,lat%ele(i)%name, lat_rev%ele(j)%name
        if(index(lat%ele(i)%name,'DET')/=0 .and. (index(lat%ele(i)%name(7:7),'W')/= 0 .or. index(lat%ele(i)%name(7:7),'E')/= 0))then
         read(lat%ele(i)%name(5:6),'(i2)')det_num
         if(index(lat%ele(i)%name(7:7),'E')/=0)det_num = 99-det_num
         write(lun,'(a16,1x,i10,1x,13es12.4)')trim(lat%ele(i)%name),det_num, lat%ele(i)%s,co(i)%vec(1),co(i)%vec(3),co(i)%vec(6), co_pos(i)%vec(1), co_pos(i)%vec(3),co_pos(i)%vec(6),&
                           co_ele(j)%vec(1), co_ele(j)%vec(3),co_ele(j)%vec(6)
        endif
       end do
      close(unit=lun)

      do i=1,lat%n_ele_track
       do j=1,lat_rev%n_ele_track
         if(lat%ele(i)%name == lat_rev%ele(j)%name .and. index(lat%ele(i)%name,'DET')/= 0)then
          write(12,'(a16,1x,13es12.4)')lat%ele(i)%name, lat%ele(i)%s,co(i)%vec(1),co(i)%vec(3),co(i)%vec(6), co_pos(i)%vec(1), co_pos(i)%vec(3),co_pos(i)%vec(6),&
                           co_ele(j)%vec(1), co_ele(j)%vec(3),co_ele(j)%vec(6)
         endif
       end do
     end do  


end