Subroutine write_header(lat,nbranch, twiss, kicker_params,tlength, pz, initial_offsets)
 use bmad
 use muon_mod
 use parameters_bmad

 implicit none
  type (lat_struct) lat
  type (g2twiss_struct) twiss
  type (initial_offsets_struct) initial_offsets
  type (kicker_params_struct) kicker_params
  real(rp) epsx, epsy, tlength, pz
  real(rp) angle
  integer j, nbranch,i
  integer n
  integer lun

  ! Print the values read from the namelist
  lun = lunget()
  print *,' write to file:', trim(directory)//'/'//'header.dat'
  open(unit=lun, file=trim(directory)//'/'//'header.dat')
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a)')        "@Beam::Longitudinal"
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,'(a,6a12)')   "    ", "tlength(ns)", "dP/P(%)", "What", "else", "???"
  WRITE(lun,'(a)')        "----------------------------------------------------------------------------"
  WRITE(lun,'(a,6f12.6)') "    ", tlength*1.e9, pz*1.e2
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a)')        "@Beam::Transverse(InflectorMidpoint)"
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,'(a,6a12)')   "    ", "beta(m)", "alpha ", "eta(m)", "etap  ", "phi(rad)", "gamma(1/m)"
  WRITE(lun,'(a)')        "----------------------------------------------------------------------------"
  WRITE(lun,'(a,6f12.6)') " X |", twiss%betax, twiss%alphax, twiss%etax, twiss%etapx, twiss%phix, twiss%gammax
  WRITE(lun,'(a,6f12.6)') " Y |", twiss%betay, twiss%alphay, twiss%etay, twiss%etapy, twiss%phiy, twiss%gammay
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a)')        "@Offsets(InflectorExit)"
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,'(a,6a12)')   "    ","x(m)", "y(m)", "t(ns)", "x'(rad)", "y'(rad)", "dP/P"
  WRITE(lun,'(a)')        "----------------------------------------------------------------------------"
  WRITE(lun,'(a,6es12.4)')"    ",initial_offsets%x_mean, initial_offsets%y_mean, initial_offsets%tmean, initial_offsets%pxmean, initial_offsets%pymean, initial_offsets%pzmean
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a)')        "@Kickers"
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,'(a,6a12)')   "   |", "By(Gauss)", "twidth(ns)", "What", "else", "???"
  WRITE(lun,'(a)')        "----------------------------------------------------------------------------"
  WRITE(lun,'(a,6f12.1)') " K1|", kicker_params%kicker_field(1)*1.e4, kicker_params%kick_width(1)*1.e9
  WRITE(lun,'(a,6f12.1)') " K2|", kicker_params%kicker_field(2)*1.e4, kicker_params%kick_width(2)*1.e9
  WRITE(lun,'(a,6f12.1)') " K3|", kicker_params%kicker_field(3)*1.e4, kicker_params%kick_width(3)*1.e9
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a)')        "@Quads"
  WRITE(lun,'(a)')        "============================================================================"
  WRITE(lun,'(a,6a12)')   "   |", "Field index", "?", "What", "else", "???"
  WRITE(lun,'(a)')        "----------------------------------------------------------------------------"
  WRITE(lun,'(a,6f12.1)') " Q1|", quad_params%long_quad_field_index(1)
  WRITE(lun,'(a,6f12.1)') " Q2|", quad_params%long_quad_field_index(2)
  WRITE(lun,'(a,6f12.1)') " Q3|", quad_params%long_quad_field_index(3)
  WRITE(lun,'(a,6f12.1)') " Q4|", quad_params%long_quad_field_index(4)
  WRITE(lun,'(a)')        "============================================================================"
  do j = 0,nbranch
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,*)
  WRITE(lun,'(a,i10)')        "@Lattice branch = ", j
  WRITE(lun,'(a)')        "================================================================================"
  WRITE(lun,'(7a16,2a22)')     "ELEMENT NAME", "r(m)", "s(m)", "angle(deg)", "total(deg)", "field_index","B_field","Tracking method","Field_calc_method"
  WRITE(lun,'(a)')        "--------------------------------------------------------------------------------"
  angle = 0
  DO i=0, lat%branch(j)%n_ele_track
     angle = angle + lat%branch(j)%ele(i)%value(angle$)
    WRITE(lun,'(a16,6f16.3,2a22)') lat%branch(j)%ele(i)%name, lat%branch(j)%ele(i)%value(rho$), lat%branch(j)%ele(i)%s, &
         lat%branch(j)%ele(i)%value(angle$)*180./pi, angle*180./pi, value_of_attribute(lat%branch(j)%ele(i), 'FIELD_INDEX'), &
         lat%branch(j)%ele(i)%value(b_field$), trim(tracking_method_name(lat%branch(j)%ele(i)%tracking_method)), trim(field_calc_name(lat%branch(j)%ele(i)%field_calc))
  END DO
  end do !branch
  n = lat%branch(nbranch)%n_ele_track
  WRITE(lun,'(a)')        "--------------------------------------------------------------------------------"
  WRITE(lun,'(a,2(3x,a,f8.5))') 'Fractional Tune:     ', 'Qx =', lat%a%tune/twopi , 'Qy =', lat%b%tune/twopi
  WRITE(lun,'(a,5(3x,a,f8.5))') 'RING PARAMS:     ', 'Qx =', lat%branch(nbranch)%ele(n)%a%phi/twopi , 'Qy =', lat%branch(nbranch)%ele(n)%b%phi/twopi, &
    'betax =', lat%branch(nbranch)%ele(n)%a%beta, 'betay =', lat%branch(nbranch)%ele(n)%b%beta, 'field_index =', (lat%branch(nbranch)%ele(n)%b%phi/twopi) **2
  WRITE(*,'(a,2(3x,a,f8.5))') 'ELE(0)%beta:     ',    'betax =', lat%branch(nbranch)%ele(0)%a%beta, 'betay =', lat%branch(nbranch)%ele(0)%b%beta
  WRITE(*,'(a)')        "--------------------------------------------------------------------------------"
  WRITE(*,'(a,2(3x,a,f8.5))') 'Fractional Tune:     ', 'Qx =', lat%a%tune/twopi , 'Qy =', lat%b%tune/twopi
  WRITE(*,'(a,5(3x,a,f8.5))') 'RING PARAMS:     ', 'Qx =', lat%branch(nbranch)%ele(n)%a%phi/twopi , 'Qy =', lat%branch(nbranch)%ele(n)%b%phi/twopi, &
    'betax =', lat%branch(nbranch)%ele(n)%a%beta, 'betay =', lat%branch(nbranch)%ele(n)%b%beta, 'field_index =', (lat%branch(nbranch)%ele(n)%b%phi/twopi) **2
  WRITE(*,'(a,2(3x,a,f8.5))') 'ELE(0)%beta:     ',    'betax =', lat%branch(nbranch)%ele(0)%a%beta, 'betay =', lat%branch(nbranch)%ele(0)%b%beta
! , &
!    " Q'x =",chrom_x," Q'y=",chrom_y
  WRITE(lun,'(a)')        "================================================================================"
  WRITE(lun,*)
  WRITE(lun,*)
  close(unit=lun)
 return
end subroutine write_header