module plot_twiss_data_mod
    use quick_plot
    use bmad

 interface
  subroutine qp_read_twiss_data (unit, err_flag, z, bx, by, eta) ! read in the data.
  use bmad
  implicit none
    real(rp), allocatable :: bx(:), by(:), z(:), eta(:), ax(:), ay(:), etap(:)
    integer n, ios, unit
    logical err_flag
 end subroutine
end interface

interface
  subroutine get_quad_locs(lat, zquad, xquad)  !zquad is the position, and xquad is positive or negative 1
  use bmad
  implicit none
  type(lat_struct) lat
  real(rp), allocatable :: zquad(:), xquad(:)
 end subroutine
end interface

contains
  subroutine plot_twiss_data(lat)
    integer id
    character(1) ans
    type(lat_struct)lat

    ! Generate PS and X-windows plots.
    call qp_open_page ("PS-L")  ! Tell \quickplot to generate a PS file.
    call plot_it(lat)              ! Generate the plot
    call qp_close_page        ! quick_plot.ps is the file name
    call qp_open_page ("X", id, 600.0_rp, 561.0_rp, "POINTS")
    call plot_it(lat)
    write (*, "(a28)", advance = "NO") " Hit any key to end program: "
    read(5, "(a1)") ans

  end subroutine
  !----------------------------------------------------------------------

  subroutine plot_it(lat)                             ! This generates the plot
    real(rp), allocatable :: bx(:), by(:), z(:), eta(:)
    real(rp), allocatable :: zquad(:), xquad(:)
    type(lat_struct) lat
    real(rp) x_axis_min, x_axis_max, y_axis_min, y_axis_max
    real(rp) t,h
    real(rp) xq
    integer x_places, x_divisions, y_places, y_divisions
    integer unit, ios
    integer i
    character(80) title
    logical err_flag
    namelist / parameters / title

    ! Read in the data
!    open (1, file = "plot.dat", status = "old")
!    read (1, nml = parameters)                  ! read in the parameters.
    unit=22
    call qp_read_twiss_data (unit, err_flag, z, bx, by, eta) ! read in the data.

    call get_quad_locs(lat, zquad, xquad)  !zquad is the position, and xquad is positive or negative 1
!    close (1)
    print *,' number of points = size(z))', size(z)
!    do i=1,size(z)
!     print '(4es12.4)', z(i),bx(i),by(i),eta(i)
!    end do
    title = 'Twiss Parameters'
!    title = '\ga\gb\gc\gd\ge\gf\gg\gh\gi\gj\gk\gl\gm\gn\go\gp\gq\gr\gs\gt\gu\gv\gw\gx\gy\gz'


    ! Setup the margins and page border and draw the title
    call qp_set_page_border (0.01_rp, 0.02_rp, 0.1_rp, 0.1_rp, "%PAGE")
    call qp_set_margin (0.07_rp, 0.05_rp, 0.05_rp, 0.05_rp, "%PAGE")
    call qp_draw_text (title, 0.5_rp, 0.9_rp, "%PAGE", "CT") 

    ! draw the bottom graph
    call qp_set_box (1, 1, 1, 3)
    call qp_calc_and_set_axis ("X", minval(z), maxval(z), 4, 8, "ZERO_AT_END")
    call qp_calc_and_set_axis ("Y", minval(bx), maxval(bx), 4, 8, "GENERAL")
    call qp_draw_axes ("s[m]", "\gb\dx\u[m]")
    call qp_set_symbol_attrib('dot', color = 'black', height = 10.0_rp)
    call qp_draw_data (z, bx, symbol_every = 0, draw_line=.true.)
    call qp_set_symbol_attrib('square', color = 'red', height = 30.0_rp)
!    t=0.5 !(maxval(z)-minval(z))/50.
     
    h=(maxval(bx)-minval(bx))/8.
   print *,' t=',t,' h= ', h
    do i=1,size(xquad)
      t = abs(xquad(i))/2.
      xq = sign(1.0_rp, xquad(i))
      call qp_draw_rectangle(zquad(i)-t,  zquad(i)+t,(xq-1)*h,(xq+1)*h, units='DATA/GRAPH',color='blue')
    end do
!    call qp_draw_data(zquad, xquad*maxval(bx)/10, draw_line=.false., symbol_every=1)
    call qp_save_state (.true.)

!    call qp_set_line_attrib ("PLOT", color = 'blue', style = 'dashed')
!    call qp_draw_data (z, bx, symbol_every = 5)
    call qp_restore_state

    ! draw the middle graph. 'star5_filled' is a five pointed star.
    call qp_save_state (.true.)
    call qp_set_box (1, 2, 1, 3)
    call qp_set_graph_attrib (draw_grid = .true.)
    call qp_set_symbol_attrib('dot', color = 'black', height = 10.0_rp)
 !   call qp_set_symbol_attrib ('star5_filled', height = 10.0_rp)
    call qp_calc_and_set_axis ("Y", minval(by), maxval(by), 4, 8, "GENERAL")
!    call qp_set_axis ("Y", -0.1_rp, 0.1_rp, 4, 2)
!    call qp_set_axis ('Y2', 1.0_rp, 100.0_rp, label = 'Y2 axis', &
!                                draw_numbers = .true., ax_type = 'LOG')
!    call qp_draw_axes ("\m1 \m2 \m3 \m4 \m5 \m6 \m7", "\fsLY\fn", title = "That Darn Graph")
    call qp_draw_axes ("s[m]", "\gb\dy\u[m]")
    call qp_draw_data (z, by, draw_line = .true., symbol_every = 0)
    call qp_restore_state

    ! draw the top graph. 'star5_filled' is a five pointed star.
    call qp_save_state (.true.)
    call qp_set_box (1, 3, 1, 3)
    call qp_set_graph_attrib (draw_grid = .true.)
    call qp_set_symbol_attrib ('dot', height = 10.0_rp)
    call qp_calc_and_set_axis ("Y", minval(eta), maxval(eta), 4, 8, "GENERAL") !    call qp_set_axis ("Y", -0.1_rp, 0.1_rp, 4, 2)
    call qp_draw_axes ("s[m]", "\gy[m]")
    call qp_draw_data (z, eta, draw_line = .true., symbol_every = 0)
    call qp_restore_state
  end subroutine

end module