! creates lst of quads, sextupoles, dipoles, wigglers,
! with gaussian distribution of offsets in format
! that can be called by bmad.

 program quad_offsets

   use bmad
   use nr

      type (lat_struct) lat
      type (ele_struct) ele 

      real(rp) offset
      real(rp) offset_mag
      real(rp) dist
      real(rp) rms_quad_offset, rms_quad_tilt, rms_sext_offset, rms_wiggler_offset
      real(rp) rms_dipole_roll
      real(rp) tilt_mag, tilt, roll_mag, roll 
      real(rp) idum/0/      
      real(rp) sum/0./,sumqoff/0./, sumqtilt/0./, sumsext/0./, sumroll/0./, sumwig/0./

      integer ix, i, j, m, ii, n_turns, index, ixx, k
      integer nargs,iargc
      integer iu,iu1
      integer nbend/0/, nquad/0/, nsext/0/, nwig/0/
      integer time(8)
      integer key

      character*140 lat_file
      character*120 line, last_line, changes(1:3)/3*' '/
      character*16 name, last_name

      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)

       print '(a,$)',' RMS quad vertical offset ? (microns) '
       read(5,*)rms_quad_offset
       rms_quad_offset = rms_quad_offset * 1.e-6

       print '(a,$)',' RMS quad tilt ? (micro-radians) '
       read(5,*)rms_quad_tilt
       rms_quad_tilt = rms_quad_tilt * 1.e-6

       print '(a,$)',' RMS sextupole vertical offset ? (microns) '
       read(5,*)rms_sext_offset
       rms_sext_offset = rms_sext_offset * 1.e-6

       print '(a,$)',' RMS dipole roll ? (micro-radians) '
       read(5,*)rms_dipole_roll
       rms_dipole_roll = rms_dipole_roll * 1.e-6

       print '(a,$)',' RMS wiggler vertical offset ? (microns) '
       read(5,*)rms_wiggler_offset
       rms_wiggler_offset = rms_wiggler_offset * 1.e-6

       if (idum < 1) then
          call date_and_time(values=time)
          idum = time(8)
       end if

      do i=1,int(idum)
      call gasdev_s(idum)
      sum = sum+idum**2
      end do
      print *,' call gasdev ', i-1, ' times', ' rms = ', sqrt(sum/(i-1))
      sum=0.

!       call ran_seed(time(8))
       print *,' idum = ', idum
       

       last_name = " "
       do i=1,lat%n_ele_max
        key = lat%ele(i)%key
        if(key /= sextupole$ .and. key /= wiggler$ .and. key /= quadrupole$ .and. key /= sbend$)cycle
        if(index(lat%ele(i)%name,'SK')/= 0)cycle
        ix=index(lat%ele(i)%name,'#')
          if(ix /= 0 .and. lat%ele(i)%name(1:ix-1) == last_name(1:ix-1))cycle
          if(ix ==0)ix = len(trim(lat%ele(i)%name))+1
          last_name = lat%ele(i)%name

        if(lat%ele(i)%key == quadrupole$)then

         call gasdev_s(idum)
         nquad=nquad+1

         offset_mag = min(abs(idum),3.)*rms_quad_offset
         offset = sign(offset_mag,idum)
!         print *,' idum= ', idum

         sumqoff = sumqoff+offset**2

         write(11,'(a,f12.8)')lat%ele(i)%name(1:ix-1)//'[y_offset]:=',offset
         if(rms_quad_tilt /= 0.)then
           call gasdev_s(idum)
           tilt_mag = min(abs(idum),3.)*rms_quad_tilt
           tilt = sign(tilt_mag,idum)
           sumqtilt=sumqtilt+tilt**2
           write(11,'(a,f12.8)')lat%ele(i)%name(1:ix-1)//'[tilt]:=',tilt
         endif
        endif

        if(lat%ele(i)%key == sextupole$ .and. rms_sext_offset /= 0.)then
         call gasdev_s(idum)
         nsext=nsext+1
         offset_mag = min(abs(idum),3.)*rms_sext_offset
         offset = sign(offset_mag,idum)
         sumsext=sumsext+offset**2
!         print *,' idum= ', idum
          write(11,'(a,f12.8)')lat%ele(i)%name(1:ix-1)//'[y_offset]:=',offset
        endif

        if(lat%ele(i)%key == sbend$ .and. rms_dipole_roll /= 0.)then
         nbend=nbend+1
         call gasdev_s(idum)
         roll_mag = min(abs(idum),3.)*rms_dipole_roll
         roll = sign(roll_mag,idum)
         sumroll=sumroll+roll**2
!         print *,' idum= ', idum
          write(11,'(a,f12.8)')lat%ele(i)%name(1:ix-1)//'[roll]:=',roll
        endif

        if(lat%ele(i)%key == wiggler$ .and. rms_wiggler_offset /= 0.)then
         nwig=nwig+1
         call gasdev_s(idum)
         offset_mag = min(abs(idum),3.)*rms_wiggler_offset
         offset = sign(offset_mag,idum)
         sumwig=sumwig+offset**2
!         print *,' idum= ', idum
          write(11,'(a,f12.8)')lat%ele(i)%name(1:ix-1)//'[y_offset]:=',offset
        endif

       end do
      if(nbend>0)print '(i4,a,a,e12.4)',nbend,' bends', ' rms roll = ', sqrt(sumroll/nbend)
      if(nquad>0)print '(i4,a,a,e12.4)',nquad,' quads', ' rms offset = ', sqrt(sumqoff/nquad)
      if(nquad>0)print '(i4,a,a,e12.4)',nquad,' quads', ' rms tilt = ', sqrt(sumqtilt/nquad)
      if(nsext>0)print '(i4,a,a,e12.4)',nsext,' sextupoles', ' rms offset = ', sqrt(sumsext/nsext)
      if(nwig>0)print '(i4,a,a,e12.4)',nwig,' wigglers', ' rms offset = ', sqrt(sumwig/nwig)
end