program eta_test

   use bmad
   use touschek_mod
   use ibs_mod
   use rad_int_common

      implicit none

      type (lat_struct) lat
      type (coord_struct), allocatable :: orbit(:)
      type (normal_modes_struct) mode, inmode, ibsmode
      type (rad_int_common_struct) rad_int

      integer ix, i, j, m, ii, n_turns, index, ixx
      integer nargs,iargc
      integer n
      integer l,k
      integer ix_ele
      integer ix_cache/0/
      integer i_dim/6/
      integer nwig, ndipole

      real(rp) I_Amps
      real(rp) f_rev/390100./,T1
      real(rp) synch_tune, Qz
      real(rp) initial_blow_up/3./
      real(rp) ratio
      real(rp) i2,i4,i5, i2w,i4w,i5w
      real(rp) Cq/3.83e-13/, gamma_f, mc2

      character(60) lat_file
      character(120) line, last_line
      character*4 formula
   
      real q_over
      real U_zero
      real E_zero
      real alpha_mom
      real k_harm
      real f_func
      real circum
      real app_max

      real aperture
      real count
      real ratio_input

      nargs = cesr_iargc()
      if(nargs == 1)then
         call cesr_getarg(1,lat_file)
         print *, 'Using ', trim(lat_file)
       else
      lat_file = 'bmad.'
      type '(a,$)',' Lattice file name ? (default= bmad.) '
      read(5,'(a)') line
       call string_trim(line, line, ix)
       lat_file = line
       if(ix == 0) lat_file = 'bmad.'
       type *, ' lat_file = ', lat_file
      endif

       call bmad_parser (lat_file, lat)

       print*, "Ratio? "
       read(*,*) ratio_input
       count = 1.
       do k=1,4
          aperture = (count/1000.)+.006
          
          count = count+1
          do n=3,9
             synch_tune = 5851.5+(1950.5*n)
!       synch_tune = 22710.
!       Qz=synch_tune/390100.
!       lat%z%tune = -Qz*twopi
       
       Qz = synch_tune/390100
       lat%z%tune = -Qz*twopi
       

 !      print*, "Aperture? "
  !     read(*,*) aperture

!      do while(line(1:1) /= 'G')
!10    print '(a, $)', '  element change or GO> '
!      read(5, '(a)',err=10) line
!     
!!$  ix = index(line, '!')
!!$  if (ix /= 0) line = line(:ix-1)        ! strip off comments
!!$
!!$  call str_upcase(line, line)
!!$  call string_trim(line, line, ix)
!!$
!!$  if (ix == 0) then       ! nothing typed. do the same thing
!!$      line = last_line
!!$  endif
!!$
!!$   last_line = line
!!$
!!$   call str_upcase(line,line)
!!$
!!$   if(line(1:1) /= 'G') call find_change( line, lat)
!!$   end do
!!$   


       call set_z_tune(lat)

       call reallocate_coord(orbit,lat%n_ele_max)
     call twiss_at_start(lat)
     call twiss_propagate_all(lat)
     call closed_orbit_calc (lat,orbit,i_dim)
     call radiation_integrals(lat,orbit,mode, ix_cache)
     call element_locator('RF_W1',lat,ix_ele)

     E_zero = lat%ele(lat%n_ele_track)%value(E_TOT$)
     circum = lat%ele(lat%n_ele_track)%s
     U_zero = mode%e_loss
     alpha_mom = mode%synch_int(1)/circum
     k_harm = floor((circum/c_light)*5e8)
     
     q_over = (4*lat%ele(ix_ele)%value(voltage$))/U_zero
     f_func = 2*(sqrt(q_over*q_over-1)-acos(1/q_over))
     app_max = sqrt((U_zero*f_func)/(twopi/2*alpha_mom*k_harm*E_zero))
     

   mc2 = mass_of (lat%param%particle)
   gamma_f = lat%ele(lat%n_ele_track)%value(E_TOT$) / mc2

   i2=0.
   i4=0.
   i5=0.
   i2w=0.
   i4w=0.
   i5w=0.
   nwig=0
   ndipole = 0

!   do i =1,lat%n_ele_track
!  write(14,'(1x,a12,3e12.4)')lat%ele(i)%name, rad_int%i2(i),rad_int%i4a(i),rad_int%i5a(i)
!  if(lat%ele(i)%key == wiggler$)then
!  write(16,'(1x,a12,3e12.4)')lat%ele(i)%name, rad_int%i2(i),rad_int%i4a(i),rad_int%i5a(i)
!   i2w = i2w + rad_int%i2(i)
!   i4w = i4w + rad_int%i4a(i)
!   i5w = i5w + rad_int%i5a(i)
!   if(rad_int%i2(i) /= 0.)nwig=nwig+1
!   else
!   i2 = i2 + rad_int%i2(i)
!   i4 = i4 + rad_int%i4a(i)
!   i5 = i5 + rad_int%i5a(i)
!   if(rad_int%i2(i) /= 0.)ndipole = ndipole + 1  
!  endif

!   enddo
!   write(15,'(1x,a9,i3,a11,3e12.4)')'  sum of ',ndipole,' dipoles = ',i2,i4,i5
!   write(15,'(1x,a9,i3,a12,3e12.4)')'  sum of ', nwig,' wigglers = ',i2w,i4w,i5w
!   write(15,'(1x,a13,3e12.4)')' emittance = ', Cq*gamma_f**2*i5/(i2-i4),Cq*gamma_f**2*(i5+i5w)/(i2+i2w-i4-i4w), mode%a%emittance
!     mode%b%emittance = mode%a%emittance * 0.01
   write(15,'(1x,a13,3e12.4)')' emittance = ', mode%a%emittance
     mode%b%emittance = mode%a%emittance * 0.01

     print '(a24,e12.4)',' horizontal emittance = ',mode%a%emittance
     print '(a22,e12.4)',' vertical emittance = ',mode%b%emittance
     print '(a11,e12.4)',' sigE/E = ',mode%sigE_E
     print '(a16,e12.4)',' bunch length = ',mode%sig_z

     if (app_max > aperture)then
          mode%pz_aperture = aperture
     else
        mode%pz_aperture = app_max
     endif

     print '(a19,e12.4,a7,e12.4,a6,e12.4,a3)',' Energy aperture = ', mode%pz_aperture, '  Qz = ', lat%z%tune/twopi,' fz = ',lat%z%tune/twopi*390.1,' kHz'
!     write(13, '(a19,e12.4,a7,e12.4,a6,e12.4,a3)')' Energy aperture = ', mode%pz_aperture, '  Qz = ', lat%z%tune/twopi,' fz = ',lat%z%tune/twopi*390.1,' kHz'
     call element_locator('RF_W1',lat,ix_ele)
     print '(a24,e12.4,a3)',' Accelerating voltage = ', 4*lat%ele(ix_ele)%value(voltage$)/1.e6,' MV'
!    write(13, '(a24,e12.4,a3)')' Accelerating voltage = ', 4*lat%ele(ix_ele)%value(voltage$)/1.e6,' MV'
!    print '(2a15)','   Current (A) ',' Lifetime(min) '
     print*, 'file:', 1000+k*100+n
  do j =1,20
   mode%b%emittance = mode%a%emittance * 0.001*j 

!   print '(/)'
   print '(a23,1x,i,1x,e12.4)',' #  mode%b%emittance = ',j,mode%b%emittance

   write(k*100+n, '(/)')
   write(k*100+n, '(/,a23,1x,i,1x,e12.4)')' #  mode%b%emittance = ',j,mode%b%emittance
   write(k*100+n, '(/)')
!   write(13, '(/)')
   
!   write(13,'(/,a23,1x,i,1x,e12.4)')' # mode%b%emittance = ',j,mode%b%emittance
!   write(13, '(/)')

   write(1000+k*100+n, '(/)')
   write(1000+k*100+n, '(/,a23,1x,i,1x,e12.4)')' #  mode%b%emittance = ',j,mode%b%emittance
   write(1000+k*100+n, '(/)')

   do i = 1,20
     I_Amps = 0.00025*i
     lat%param%n_part = I_Amps/f_rev/1.6e-19
     ibsmode = mode
 !     call touschek_lifetime (ibsmode, T1, lat, orbit)

!     print '(1x,e12.4,4x,e12.4)',I_Amps, T1/60.
!     write(13, '(1x,e12.4,4x,e12.4)')I_Amps, T1/60.


     inmode=mode
     formula = 'cimp'
     ratio = ratio_input
     
    
     print '(2(a20,e12.4))',' mode%a%emittance = ', mode%a%emittance,' mode%b%emittance = ', mode%b%emittance
     print '(2(a15,e12.4))','  mode%sig_z = ',mode%sig_z,' mode%sigE_E = ', mode%sigE_E

     call ibsequilibrium2(lat,inmode, ibsmode, formula, ratio, initial_blow_up)
!     print '(5e12.4)', I_Amps,mode%a%emittance, ibsmode%a%emittance, &
!                              mode%b%emittance, ibsmode%b%emittance
     write(1000+k*100+n,'(5e12.4)') I_Amps,mode%a%emittance, ibsmode%a%emittance, &
                              mode%b%emittance, ibsmode%b%emittance
     
      call touschek_lifetime (ibsmode, T1, lat, orbit)

!     print '(1x,e12.4,4x,e12.4)',I_Amps, T1/60.
     write(k*100+n, '(1x,e12.4,4x,e12.4)')I_Amps, T1/60.

   end do

  end do
end do
end do
     end