program tilt_analysis

 use precision_def

 implicit none
 
  character*100 list, file_names(100)
  character*140 line
  character*7 bpm

  real(rp) sum_tilt(0:100)
  real(rp) avg_tilt(0:100)
  real(rp) sig_tilt(0:100)
  real(rp) bp_tilt, junk
  real(rp) tilt(0:100, 99)
  real(rp) sig_all_tilt
  real(rp) sum_all_tilt
  real(rp) all_tilt_avg

  integer i,j
  integer n(0:100)
  integer tot(0:100)
  integer ind
  integer ix
  integer tot_all_tilt

  sum_tilt = 0
  sig_tilt = 0
  sig_all_tilt = 0
  tot_all_tilt = 0
  sum_all_tilt = 0

  print '(a,$)',' List of files to read ? '
  accept *,list

  open(unit=1, file = trim(list), readonly)
  open(unit=73,file = 'histogram_tilts.dat')

  i=1
  do
   read(1,'(a)', end = 99)file_names(i)
   i=i+1
  end do

99 continue
close(unit=1)


  tot=0
  sig_tilt = 0.

  do j = 1,i-1
   open(unit=1, file=file_names(j), readonly)
   print '(a,a)',' file = ',file_names(j)
   do
     read(1, '(a)', end=199)line
     if(index(line,'|')/=0 .or. line == '')cycle
     call string_trim(line, line, ix)
     bpm = line(1:ix)
     call string_trim(line(ix+1:),line,ix)
     read(line(1:ix),*)bp_tilt
     call string_trim(line(ix+1:), line,ix)
     read(line(ix+1:),*)ind
     if(ind > 100)cycle
     if(abs(bp_tilt)<0.1) then
        sum_tilt(ind) = sum_tilt(ind) + bp_tilt 
        sum_all_tilt  = sum_all_tilt + bp_tilt
         tilt(ind,j) = bp_tilt
         tot(ind) = tot(ind)+1
         tot_all_tilt = tot_all_tilt +1
       else
         write(*,'(a,a7,a,es12.4)')' Tilt for ',bpm,' > 0.1rad = ', bp_tilt
     endif
     write(73,'(es12.4)')bp_tilt
   end do
 199   continue  
 end do
 close(unit=73)
 close(unit=1)

 open(unit=74,file='all_tilt.dat')
 open(unit=75,file='bpm_tilts.dat')

 all_tilt_avg = sum_all_tilt/tot_all_tilt
   do ind = 1,100
    do j=1,i-1
     if(tot(ind) <= 0)cycle
     if(abs(tilt(ind,j))>=0.1)cycle
     avg_tilt(ind) = sum_tilt(ind)/tot(ind)
     sig_tilt(ind) = sig_tilt(ind) + (tilt(ind,j)-avg_tilt(ind))**2
     sig_all_tilt = sig_all_tilt + (tilt(ind,j) - all_tilt_avg)**2
     end do
     sig_tilt(ind) = sqrt(sig_tilt(ind)/tot(ind))
     write(75,'(a17,1x,i3,1x,a1,es12.4)')' change bpm_tilt ',ind,'@',avg_tilt(ind)
     write(76,'(a17,1x,i3,1x,a1,es12.4)')' change bpm_tilt ',ind,'@',-avg_tilt(ind)
     write(74,'(i,2es12.4)')ind,avg_tilt(ind),sig_tilt(ind)
    end do
     sig_all_tilt = sqrt(sig_all_tilt/tot_all_tilt)
     write(74,'(a,a,i,a,e12.4)')'#',' Total number of tilts = ',tot_all_tilt,'   RMS = ',sig_all_tilt
     print '(a,a,i,a,e12.4)','#',' Total number of tilts = ',tot_all_tilt,'   RMS = ',sig_all_tilt
 close(unit=75)
 close(unit=74)   


  end