subroutine write_loop (vparam_id, vparam_label, vparam_col_ind, nmuons,kicker_params, moment, tot_inf_scatter, tot, twiss, initial_offsets, &
  xmax, ymax, inflector_angle,lat)

  use bmad
  use muon_mod
  use muon_interface
  use parameters_bmad

  implicit none
  type (lat_struct) lat
  type (g2twiss_struct) twiss
  type (g2moment_struct) moment
  type (initial_offsets_struct) initial_offsets
  type (kicker_params_struct) kicker_params

  integer, save :: lun32
  integer vparam_id, vparam_col_ind, nmuons
  integer tot_inf_scatter, tot
  character*20 vparam_label 
  logical first/.true./
  real(rp) xmax, ymax
  real(rp) inflector_angle

   if(first)then
    lun32 = lunget()
    open(unit = lun32, file = trim(directory)//'/'//'VparamDependence.dat')
    write(lun32,'(a,i4)')'vparam_id =', vparam_id
    write(lun32,'(a,a1,a,a1)')'set xlabel ','"',vparam_label,'"'
    write(lun32,'(a,i5)')'column1 =',vparam_col_ind
    write(lun32,'(a,i10)')'nmuons =', nmuons
    write(lun32,'(a)')'exit'
    write(lun32,'(/,a,/)')'#'
    write(lun32,'(16a16,9a12,a12,a12, a12, a12,a12,a12,a12,a12,a12, a12, 2a12)')  &
   'K1[Guass]-1', 'K2[Gauss]-2', 'K3[Gauss]-3', &
   'depth<x>-4',' depth<xp>-5','depth<y>-6','depth<yp>-7','depth<t>-8','depth<pz>-9', &
   'depth_SIG(1,1)-10','depth_SIG(3,3)-11','depth_SIG(5,5)-12','depth_SIG(6,6)-13', &
   'n muons-14',' after_inf-15','remaining-16','beta x-17','eta x-18','kicker_field-19','beta y-20',' kickWid-21','dtRise-22','dtFall-23','x init','px init', &
   'inf_field','pz init','etap','alphax','Kick_tStart','Quad [n]',' xmax ',' ymax ', 'inf angle','dp/p','x tune','y tune'
   first=.false.
  endif !open lun32

  write(lun32,'(3es16.3,6es16.4,4es16.4,3i16,2es12.4,7es12.4, 2es12.4, 5es12.4, 6es12.4, es12.4, 2es12.4)') &
  kicker_params%kicker_field(1) * 1.e4, & ! Gauss
  kicker_params%kicker_field(2) * 1.e4, & ! Gauss
  kicker_params%kicker_field(3) * 1.e4, & ! Gauss
  (moment%max_ave-moment%min_ave)                          *1.e3, & ! mm or mrad
  (sqrt(moment%max_sigma(1,1))-sqrt(moment%min_sigma(1,1)))*1.e3, & ! mm or mrad
  (sqrt(moment%max_sigma(3,3))-sqrt(moment%min_sigma(3,3)))*1.e3, & ! mm or mrad
  (sqrt(moment%max_sigma(5,5))-sqrt(moment%min_sigma(5,5)))*1.e3, & ! mm or mrad
  (sqrt(moment%min_sigma(6,6))), & ! mm or mrad
  nmuons,   tot_inf_scatter, tot, &
  twiss%betax, twiss%etax, maxval(kicker_params%kicker_field(1:3))*10000, twiss%betay, kicker_params%kick_width(1)*1.e9,&
  kicker_params%dtRise(1)*1.e9, kicker_params%dtFall(1)*1.e9, &
  initial_offsets%x_mean*1.e3, initial_offsets%pxmean*1.e3, inflector_field, initial_offsets%pzmean*1.e2, twiss%etapx, &
  twiss%alphax, kicker_tStart(1)*1.e9, quad_params%long_quad_field_index(1), xmax, ymax, inflector_angle, moment%ave(6), &
  lat%a%tune/twopi, lat%b%tune/twopi

  return
 end