program plot_da

 use precision_def
 use quick_plot
 use da_data_mod

implicit none

interface
subroutine read_dafile(iunit, dadata)
 use precision_def
 use da_data_mod
type (da_data_struct) dadata

integer iunit
end subroutine read_dafile
end interface


type init_struct
  real(rp) emittance,beta_x,beta_y,alpha_x,alpha_y
  integer type
  character*100 file
end type

type (da_data_struct) da_data(3), da_refdata(3), dadata
type (qp_point_struct) origin
type (qp_line_struct) line(3)
type (qp_symbol_struct) symbol(3)
type (init_struct) dr, ler, def, undulator, pretz
character*100 text(3) , text_ellipse(3)
character*100 text_frf  

character*18 def_prefix/'da'/, prefix, ref_prefix/'da_100turns'/
character*1 suffix(3)/'1','2','3'/
integer ios, iunit/1/
integer i,j, id
integer n,ix
integer mult/1/
integer iu, iostat

  character*100 file_name, plot_file
  character*3 device, ans 
  character*7 ctune(3)/'Q\dx\u','Q\dy\u','Q\dz\u'/
  character*100 string
  character*100 directory

  real(rp) xlen/800./, ylen/540./
  real(rp) xmax/0./,ymax/0./
  real(rp)  ellipse_x(100), ellipse_y(100)
  real(rp) rgamma_x, rgamma_y
  real(rp) dy

  namelist/input/ def

  dr%emittance=1.e-6; dr%beta_x=19.2241; dr%beta_y=73.4635; dr%alpha_x = -0.1355; dr%alpha_y=-0.0349; dr%type=1
  ler%emittance = 1.3e-9; ler%beta_x = 41.2; ler%beta_y=15.7; ler%alpha_x=0.0357; ler%alpha_y=0.0940; ler%type=2
  undulator%emittance = 5.4e-8;    undulator%beta_x =6.76;    undulator%beta_y = 15.02;    undulator%alpha_x=-0.094; undulator%alpha_y=0.3; undulator%type=2
  pretz%emittance = 1.545e-7;    pretz%beta_x =3.0276;    pretz%beta_y = 34.1638;    pretz%alpha_x=0.1524; pretz%alpha_y=-0.4391; pretz%type=3
!  dr%file = '/Users/dlr10/development/bmad_dist/plot_da/dr.dat'
!  ler%file = '/Users/dlr10/development/bmad_dist/plot_da/ler.dat'
!  undulator%file = '/Users/dlr10/development/bmad_dist/plot_da/undulator.dat'
!  pretz%file = '/Users/dlr10/development/bmad_dist/plot_da/pretz.dat'
  dr%file = '/home/dlr/development9_linux/plot_da/dr.dat'
  ler%file = '/home/dlr/development9_linux/plot_da/ler.dat'
  undulator%file = '/home/dlr/development9_linux/plot_da/undulator.dat'
  pretz%file = '/home/dlr/development9_linux//plot_da/pretz.dat'

  call read_string(' File prefix ? (Def = "da") ',def_prefix, prefix)

  print '(a,$)',' ilc damping ring, ler, undulator, pretz params (i,l, und, pretz) ? '
  read(*,'(a)')string
  if(string(1:1) == 'i')def = dr
  if(string(1:1) == 'l')def = ler
  if(string(1:1) == 'u')def = undulator
  if(string(1:1) == 'p')def = pretz
  iu = lunget()
  print *,' file ',def%file
  open(unit=iu,file = def%file, status='old', IOSTAT = ios)
  if(ios /=0) then
      print *,' Use default values'
     else
      print *,' Use most recent values '
      read (iu, nml = input)
  endif
  close(iu)
  
    
  print '(a,e12.4,a,$)', ' emittance = ', def%emittance, ' '
  read(*,'(a)') string
  call string_trim(string, string, ix)
  if(ix /= 0)read(string,*)def%emittance


  print '(a,f10.2,a,$)', ' beta_x = ', def%beta_x, ' '
  read(*,'(a)') string
  call string_trim(string, string, ix)
  if(ix /= 0)read(string,*)def%beta_x

  print '(a,f10.2,a,$)', ' beta_y = ', def%beta_y, ' '
  read(*,'(a)') string
  call string_trim(string, string, ix)
  if(ix /= 0)read(string,*)def%beta_y

  print '(a,f10.2,a,$)', ' alpha_x = ', def%alpha_x, ' '
  read(*,'(a)') string
  call string_trim(string, string, ix)
  if(ix /= 0)read(string,*)def%alpha_x

  print '(a,f10.2,a,$)', ' alpha_y = ', def%alpha_y, ' '
  read(*,'(a)') string
  call string_trim(string, string, ix)
  if(ix /= 0)read(string,*)def%alpha_y

  rgamma_x=sqrt((1+def%alpha_x**2)/def%beta_x)
  rgamma_y=sqrt((1+def%alpha_y**2)/def%beta_y)
  print *,' rgamma_x = ', rgamma_x,' rgamma_y = ', rgamma_y

  open(unit=iu,file = def%file)
   write(iu,'(a)')'&input'
   write(iu,'(a,es12.4)')'def%emittance = ',def%emittance
   write(iu,'(a,es12.4)')'def%beta_x = ',def%beta_x 
   write(iu,'(a,es12.4)')'def%beta_y = ',def%beta_y 
   write(iu,'(a,es12.4)')'def%alpha_x = ',def%alpha_x 
   write(iu,'(a,es12.4)')'def%alpha_y = ',def%alpha_y 
   write(iu,'(a)')'/'
 close(iu)

 ios=0 
do i=1,3
  file_name = trim(prefix)//suffix(i)//'.dat'
  open(unit=iunit, file = file_name, status='old')
  if(ios /= 0)cycle
  dadata = da_data(i)
!  print *,'1 filename  = ',file_name
  call read_dafile(iunit, dadata)
  da_data(i) = dadata
  close(unit=1)
  n=da_data(i)%npoints
  da_data(i)%x(1:n) = da_data(i)%x(1:n)*1000 *rgamma_x
  da_data(i)%y(1:n) = da_data(i)%y(1:n)*1000 *rgamma_y

end do
  
  call read_string(' Reference file prefix ? (Def = "da") ',ref_prefix, prefix)

da_refdata(1)%npoints=0
do i=1,3
  file_name = trim(prefix)//suffix(i)//'.dat'
  open(unit=iunit, file = file_name, status='old', iostat=ios)
  if(ios /= 0)exit
  call read_dafile(iunit, da_refdata(i))
  close(unit=1)
  n=da_refdata(i)%npoints
  da_refdata(i)%x(1:n) = da_refdata(i)%x(1:n)*1000*rgamma_x
  da_refdata(i)%y(1:n) = da_refdata(i)%y(1:n)*1000*rgamma_y

end do
  if(def%type == 2)mult=10
  if(def%type == 3)mult= 5
  do i=1,100
  ellipse_x(i) = mult*sqrt(def%emittance) * cos(pi*i*0.01) * 1000
  ellipse_y(i) = mult*sqrt(def%emittance) * sin(pi*i*0.01) * 1000
  end do




  device = 'X'
do 
 
 if(device(1:1) == 'X')call qp_open_page('X', id, xlen, ylen,'POINTS')
 if(device(1:2) == 'PS') call qp_open_page('PS-L',plot_file='da.ps')
 if(device(1:3) == 'GIF') call qp_open_page('GIF-L',plot_file='da.gif')

  call qp_set_page_border (0.02_rp, 0.04_rp, 0.04_rp, 0.04_rp, '%PAGE')
  call qp_set_margin (0.035_rp, 0.02_rp, 0.035_rp, 0.015_rp, '%PAGE')

   call qp_set_box (1,1,1,1)

 xmax = (maxval(da_data(1)%x(1:da_data(1)%npoints)))
 xmax = ((xmax/1)*1+1)
 ymax = (maxval(da_data(1)%y(1:da_data(1)%npoints)))
 ymax = ((ymax/1)*1+1)
 call qp_calc_and_set_axis ('X', -xmax,xmax, 4, 8, 'GENERAL')
 call qp_calc_and_set_axis ('Y', 0.0_rp,ymax, 4, 8, 'GENERAL')
 call qp_set_text_attrib('MAIN_TITLE',  height=25.0_rp)
 call qp_set_text_attrib('GRAPH_TITLE', height=20.0_rp)
 call qp_set_text_attrib('LEGEND', color=black$,height=15.0_rp)
 call qp_set_text_attrib('AXIS_NUMBERS', color=black$,height=20.0_rp)
 call qp_set_text_attrib('AXIS_LABEL', color=black$,height=20.0_rp)
 call qp_set_text_attrib('TEXT', color=black$,height=20.0_rp)
 call qp_set_line_attrib('AXIS',pattern=1, width=4,color=black$)  
 call qp_draw_axes ("horizontal amplitude [mm]", "vertical amplitude [mm]",trim(da_data(1)%lattice), draw_grid = .false.)


 do j=1,3
   n=da_data(j)%npoints
 call qp_set_symbol_attrib (type=j+3, color = j, fill_pattern=3, height = 20.0_rp)
  call qp_set_line_attrib('PLOT',pattern=1, width=3,color=j)  
  call qp_draw_data (da_data(j)%x(1:n), da_data(j)%y(1:n),draw_line = .true., symbol_every = 1)
  line(j)%pattern = 1  
  line(j)%width=3
  line(j)%color = j
  symbol(j)%type = j+3
  write(string,'(f5.1)')da_data(j)%delta_e    
!  print *,da_data(j)%delta_e, 'string = ', trim(string)
  text(j) = '\gd\fnE/E = '//trim(string)//' %'
if(da_refdata(1)%npoints /= 0) &
  call qp_draw_data (da_refdata(j)%x(1:n), da_refdata(j)%y(1:n),draw_line = .true., symbol_every = 0)
  
  end do

 origin%x= 0.70_rp
 origin%y= 0.98_rp
 origin%units = "%GRAPH"
 print *,' origin = ',origin
 
 call qp_draw_curve_legend (x_origin=origin%x,y_origin=origin%y, line= line, symbol=symbol, text= text, units='%GRAPH')

  call qp_draw_text (trim(da_data(1)%data_date), 0.01_rp, 1.02_rp,"%GRAPH",height=10.0_rp)
  call qp_draw_text (trim(directory), 0.79_rp, 1.02_rp,"%GRAPH",height=10.0_rp)
! print tunes and energy key
 do j=1,3
  write(string,'(f6.3)')da_data(1)%tune(j)  
  text(j) = ctune(j)//' ='//trim(string)

  call qp_set_line_attrib('PLOT',pattern=2,width=3,color=j+3)  
  call qp_draw_data(ellipse_x(1:100)*j,ellipse_y(1:100)*j, draw_line=.true., symbol_every=0, clip=.true.)
  line(j)%pattern = 2
  line(j)%color = j+3

 end do

 call qp_draw_text_legend (text(1:3), 0.10_rp, 0.98_rp, "%GRAPH")

 text_ellipse(1) = '\gs\fn\dinj\u'
 text_ellipse(2) = '2\gs\fn\dinj\u'
 text_ellipse(3) = '3\gs\fn\dinj\u'

if(def%type == 2)then
 text_ellipse(1) = '10\gs\fn'
 text_ellipse(2) = '20\gs\fn'
 text_ellipse(3) = '30\gs\fn'
endif
if(def%type == 3)then
 text_ellipse(1) = '5\gs\fn'
 text_ellipse(2) = '10\gs\fn'
 text_ellipse(3) = '15\gs\fn'
endif

 origin%y=0.86_rp
 call qp_draw_curve_legend(x_origin=origin%x, y_origin=origin%y,line = line, text=text_ellipse, units='%GRAPH')

 text(1) = 'horiz amp = x\gg\dh\u\fn\u1/2'
 text(2) = 'vert amp = y\gg\dv\u\fn\u1/2'
 call qp_draw_text_legend (text(1:1),0.1_rp,0.8_rp,"%GRAPH")
 call qp_draw_text_legend (text(2:2),0.1_rp,0.75_rp,"%GRAPH")
! print *,da_data(1)%delta_frf
 write(string,'(f10.1)')da_data(1)%delta_frf
! print *,string
 text_frf = '\gd\fnf\dRF\u = '//trim(string)//' Hz'
  if(da_data(1)%delta_frf /= 0)call qp_draw_text (trim(text_frf), 0.1_rp, 0.65_rp,"%GRAPH",height=15.0_rp)


 if(device /= 'X')then
   call qp_close_page
   exit
 endif
 print '(a,$)' ,' Postscript (p) or Gif (g) '
 call read_string('Postscript (p) or Gif (g) ? ', 'no',ans)
 call upcase_string(ans)
 if(index(ans,'N') /= 0) exit
 if(index(ans,'P')/= 0)device = 'PS' 
 if(index(ans,'G')/= 0)device = 'GIF' 
 print '(a,$)',' Directory ? '
 read(5,'(a)')directory
end do

end


!------------------------------------------------
subroutine read_dafile(iunit, dadata)
 use precision_def
 use da_data_mod
 implicit none

type (da_data_struct) dadata

character*100 string, ctune(3),crf
character*1 plane
character*16 lost_at
integer n,ios,ix1,ix2
integer iunit
integer i
integer ix/0/
logical data_flag

data_flag=.false.
n=0
  do
   read(iunit, '(a)', iostat=ios)string
    if(ios ==-1)exit
    if(.not.data_flag)then
    ix1=index(string,"`")
    ix2=index(string(ix1+1:),"'")
    if(ix1 /=0)then
!      print '(a,a)','subroutine read_datafile: string = ', string
      if(index(string,'Lattice')/= 0)dadata%lattice = string(ix1+1:ix1+ix2-1)
      if(index(string,'data_date')/= 0)dadata%data_date = string(ix1+1:ix1+ix2-1)
      if(index(string,'Q_x')/= 0)ctune(1) = string(ix1+1:ix1+ix2-1)
      if(index(string,'Q_y')/= 0)ctune(2) = string(ix1+1:ix1+ix2-1)
      if(index(string,'Q_z')/= 0)ctune(3) = string(ix1+1:ix1+ix2-1)
      if(index(string,'Delta_fRF')/=0)crf=string(ix1+1:ix1+ix2-1)
    endif

    if(index(string,'dE_E')/=0)then
     ix1 = index(string, '=')
     call from_string(string(ix1+1:), dadata%delta_e)
     dadata%delta_e = dadata%delta_e * 100
    endif
    if(index(string,'!')/=0)data_flag = .true.
   else
    n=n+1
       read(string(1:30),*)dadata%x(n), dadata%y(n), dadata%turn(n)
   endif
  end do
  dadata%npoints=n

  do i=1,3
   print '(a)',ctune(i)
   call from_string(ctune(i), dadata%tune(i))
  end do
   call from_string(crf, dadata%delta_frf)
 return
 end subroutine

 subroutine from_string(string, x)

 use precision_def
  implicit none

  character*(*) string
  character*40 word
  real(rp) x
 
  integer ix1, ix2
  integer exp

  ix1=index(string,'=')
  ix2 = index(string,'E')

  word = string(ix1+1:)
  if(ix2 /= 0) then
   read(string(ix2+1:),*)exp
   read(word,*)x
!   print *,' exp = ',exp, ' x =', x
!   x= x*10.**exp
  else
   read(word, *)x
 endif

  return
end subroutine