!........................................................................
!
! Subroutine : MINIMIZE_MOMENTS(LAT, MOMNT, SEXT_VAR, CON)
!
! Description: Compute sextupoles that minimize moments, and energy and 
!              amplitude dependence
!
! Arguments  : 
!     INPUT:
!        LAT     -- lat_struct : Lat
!        MOMNT   -- Moment_struct: momnt         
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
!             
! save sextupole values
! enforce sextupole symmetry?
! find matching sextupoles
! setup useable_ele array
! setup matrix
! find which sextupoles to use
!------------------------------------------------------------
! sext_res
! sext moments
! k2^4 term
! sext symmetry
! sext antisymmetry
! sign_symmetry
! chrom
! tonality
! coupling_a_real
! coupling_a_image
! coupling_b_real
! coupling_b_image
! x plane delta_beta
! y plane delta_beta
! x plane dbeta_dpretz
! y plane dbeta_dpretz
! x plane dbeta_dcos
! y plane dbeta_dcos
! x plane dbeta_dsin
! y plane dbeta_dsin
!--------------------------------------------------------------
! find optimal solution
! find better sextupole distribution
! update new results
!----------------------------------------------------------------------------
! write to file
! calculate new chi2
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!
!........................................................................
!
!
! $Log$
! Revision 1.13  2007/01/30 16:15:13  dcs
! merged with branch_bmad_1.
!
! Revision 1.9  2006/11/16 18:55:45  mjf7
! Bmad change in ele_struct name length caused bugs.
!
! Revision 1.8  2006/02/08 19:09:46  mjf7
! Replaced obsolete bmad subroutine calls with new versions, and found variables which needed a save attribute.
!
! Revision 1.7  2004/11/08 19:16:29  dlr
! replace n_ele_maxx with n_ele_max
!
! Revision 1.6  2004/05/04 15:24:23  dcs
! ele%position -> ele%floor.
!
! Revision 1.5  2004/01/30 19:15:42  dlr
! correct calculation of electron positron differences
! correct sextupole optimization to deal properly with both beams
! increase gjdet array size for making scmating knobs
! add new constraints eletron emittance
!
! Revision 1.4  2003/09/30 12:51:39  dcs
! Update to correspond with new floor structure.
!
! Revision 1.3  2003/07/08 19:27:54  mjf7
!
!
! Modified all subroutines to correctly use allocatable lat elements. n_ele_maxx is no longer global, but a member variable of the lat struct. -mjf
!
! Revision 1.2  2003/04/30 17:14:52  cesrulib
! dlr's changes since last import
!
! Revision 1.1.1.1  2002/12/13 19:23:29  cesrulib
! import bmadz
!
!
!........................................................................
!


!subroutine minimize_moments (lat_pos, lat_ele, momnt_positron, momnt_electron, sext_var, con)
subroutine minimize_moments (lat, momnt, sext_var, con)

  use moments_com_mod
  use nr

  use bmad


  use nonlin_mod
  use constraints_mod

  implicit none

  type (lat_struct)  lat !, lat_pos, lat_ele
  type (lat_struct), save :: lat_old
!  type (moment_struct) momnt, momnt_positron, momnt_electron
  type (moment_struct) momnt
  type (indep_var_struct)  sext_var
  type (constraint_struct)  con
  type (sextupole_res_struct), save ::  res, res_e

  integer n_sext, n_var, i, i2, jsp, n_res, j, k, key
  integer n1, n2, n3, n4, n5, n6, n7, n8
  integer nv_use, ix2, n_row
  integer nv, nn, ix_w, ix_e
  integer use_to_i(275)
  integer ix, east_sext(49), west_sext(49), west_to_east(275)
  integer ix_sp(8000)
  integer ix_ele_sp(lat%n_ele_max), ix_species

  real(rp) weight(8000), sig(8000), covar(275, 275), alpha(275, 275)
  real(rp) a_lambda, wgt, l1, l2
  real(rp) w(275), v(275,275), b(8000), c(8000), x(275), bb(8000)
  real(rp) R0(8000,275), fixed_sext(275)
  real(rp) old_chi2, chi2, phi_res, res_factor, jz_factor, tune_vec(3)

  character*40 name
  character*60 fmt1, fmt2

  logical eastwest_sym, useable_ele(lat%n_ele_max)
  logical, save :: been_here = .false.
  external eval_sex_moments


!dlr 4/12/03  momnt = momnt_positron
!dlr 4/12/03  lat = lat_pos

! save sextupole values
   lat_old = lat

! enforce sextupole symmetry?

  eastwest_sym = con%sextupole_symmetry

! find matching sextupoles

  west_sext = 0
  east_sext = 0
  do i = 1, momnt%n_sext
    name = momnt%name(i)
    if (name(7:7) == 'W' .and. name(1:3) == 'SEX') then
      read (name(5:6), *) ix
      west_sext(ix) = i
    elseif (name(7:7) == 'E' .and. name(1:3) == 'SEX') then
      read (name(5:6), *) ix
      east_sext(ix) = i
    endif
  enddo

  west_to_east = 0
  where (west_sext /= 0 .and. east_sext /= 0) &
                                  west_to_east(west_sext(:)) = east_sext(:)

  if (eastwest_sym) then
    n_var = count(west_to_east /= 0)
  else
    n_var = momnt%n_sext
  endif

! setup useable_ele array

  do k = 1, lat%n_ele_track
    key= lat%ele(k)%key
    if (key == quadrupole$ .or. key == rfcavity$ .or.  &
                                                     key ==  marker$) then
      useable_ele(k) = .true.
    else
      useable_ele(k) = .false.
    endif
  enddo

! setup matrix

  b(:) = 0.0
  row(:,:) = 0.0
  n_row = 0
  k2_4_row = 0

! find which sextupoles to use

  n_sext = momnt%n_sext

  do i = 1, n_sext
    momnt%use(i) = .false.
!dlr 4/12/03    momnt_positron%use(i) = .false.
!dlr 4/12/03    momnt_electron%use(i) = .false.
    do j = 1, sext_var%n_var
      if (momnt%name(i) == sext_var%v(j)%ele_name)then
          momnt%use(i) = .true.
!dlr 4/12/03          momnt_positron%use(i) = .true.
!dlr 4/12/03          momnt_electron%use(i) = .true.
      endif
    end do
  end do


! loop first positron moments and then electron moments

!dlr 4/12/03 this loop
! electron_positron_loop : do ix_species = 1, 2
!  if(ix_species == 1)then
!    lat = lat_pos
!    momnt = momnt_positron
!  else
!   lat = lat_ele
!   momnt = momnt_electron
!  endif

    

!------------------------------------------------------------
! sext_res

  do jsp = 1, con%n_constraint

    if (con%c(jsp)%weight == 0) cycle

    if (con%c(jsp)%variable == sext_resonances$) then

      call lat_geometry (lat)

      if (lat%z%tune == 0) then
        print *, 'WARNING FROM MINIMIZE_MOMENTS: SYNCH TUNE IS ZERO!'
      endif

      do k = 1, lat%n_ele_track
        lat%ele(k)%z%phi = lat%ele(k)%floor%theta * lat%z%tune
      enddo

      ix = momnt%ix(1)
      call sextupole_resonance_calc (lat%ele(ix), res)

      wgt = con%c(jsp)%weight
      do k = 1, res%n_res
        tune_vec = (/ lat%a%tune, lat%b%tune, lat%z%tune /)
        phi_res = dot_product(res%x(k)%iq, tune_vec) / 2
        res_factor = max (abs(sin(phi_res)), sin(con%del_q_res_min/2))
        jz_factor = (sqrt(con%jz_ja_betaz_ratio))**res%x(k)%ij(3)
        res%x_weight(k) = jz_factor * wgt / res_factor

        phi_res = dot_product (res%y(k)%iq, tune_vec) / 2
        res_factor = max (abs(sin(phi_res)), sin(con%del_q_res_min/2))
        jz_factor = (sqrt(con%jz_ja_betaz_ratio))**res%y(k)%ij(3)
        res%y_weight(k) = jz_factor * wgt / res_factor
      enddo

      n_res = res%n_res

      n1 = n_row + 1
      n2 = n_row + n_res
      weight(n1:n2)   = res%x_weight(:)
      ix_sp(n1:n2) = jsp
      n_row = n2

      n3 = n_row + 1
      n4 = n_row + n_res
      weight(n3:n4)   = res%x_weight(:)
      ix_sp(n3:n4) = jsp
      n_row = n4

      n5 = n_row + 1
      n6 = n_row + n_res
      weight(n5:n6)   = res%y_weight(:)
      ix_sp(n5:n6) = jsp
      n_row = n6

      n7 = n_row + 1
      n8 = n_row + n_res
      weight(n7:n8)   = res%y_weight(:)
      ix_sp(n7:n8) = jsp
      n_row = n8

      do i = 1, n_sext

        ix = momnt%ix(i)
        call sextupole_resonance_calc (lat%ele(ix), res)

        if (eastwest_sym) then
          if (west_to_east(i) == 0) exit
          i2 = momnt%ix(west_to_east(i))
          call sextupole_resonance_calc (lat%ele(i2), res_e)
          row(n1:n2, i) = res%a%cos + res_e%a%cos
          row(n3:n4, i) = res%a%sin + res_e%a%sin
          row(n5:n6, i) = res%b%cos + res_e%b%cos
          row(n7:n8, i) = res%b%sin + res_e%b%sin
        else
          row(n1:n2, i) = res%a%cos
          row(n3:n4, i) = res%a%sin
          row(n5:n6, i) = res%b%cos
          row(n7:n8, i) = res%b%sin
        endif

      enddo

    endif

! sext moments

    if (con%c(jsp)%variable == sext_moments$) then
      n1 = n_row + 1
      n2 = n_row + 10
      weight(n1:n2) = con%c(jsp)%weight
      ix_sp(n1:n2) = jsp
      call row_stuffit (row(n_row+1, :), momnt%m1)
      call row_stuffit (row(n_row+2, :), momnt%m2)
      call row_stuffit (row(n_row+3, :), momnt%m3)
      call row_stuffit (row(n_row+4, :), momnt%m4)
      call row_stuffit (row(n_row+5, :), momnt%m5)
      call row_stuffit (row(n_row+6, :), momnt%m6)
      call row_stuffit (row(n_row+7, :), momnt%m7)
      call row_stuffit (row(n_row+8, :), momnt%m8)
      call row_stuffit (row(n_row+9, :), momnt%m9)
      call row_stuffit (row(n_row+10, :), momnt%m10)
      n_row = n_row + 10
    endif

! k2^4 term

    if (con%c(jsp)%variable == sex_k2_4$) then
      call next_row (jsp)
      k2_4_row = n_row
    endif

! sext symmetry

    if (con%c(jsp)%variable == sext_symmetry$) then
      if (.not. eastwest_sym) then  ! only do for no symmetry
        do i = 1, n_sext
          if (west_to_east(i) /= 0) then
            call next_row (jsp)
            row(n_row, i) = 1.0
            row(n_row, west_to_east(i)) = -1.0
          endif
        enddo
      endif
    endif

! sext antisymmetry

    if (con%c(jsp)%variable == sext_antisymmetry$) then
      if (.not. eastwest_sym) then  ! only do for no symmetry
        do i = 1, n_sext
          if (west_to_east(i) /= 0) then
            call next_row (jsp)
            row(n_row, i) = 1.0
            row(n_row, west_to_east(i)) = 1.0
          endif
        enddo
      endif
    endif

! sign_symmetry

    if (con%c(jsp)%variable == sign_symmetry$) then
      if (.not. eastwest_sym) then  ! only do for no symmetry
        call next_row (jsp)
        do i = 1, n_sext
          ix = momnt%ix(i)
          if (con%c(jsp)%where1_ix <= ix .and. &
                                  ix <= con%c(jsp)%where2_ix) then

            if(abs(momnt%chrom_x(i)) > abs(momnt%chrom_y(i)))then
              row(n_row,i) = 1./momnt%s_h
              row(n_row,n_sext+1-i) = 1./momnt%s_h
             else
              row(n_row,i) = -1./momnt%s_v
              row(n_row,n_sext+1-i) = -1./momnt%s_v
            endif
          endif
        end do
      endif
    endif

! chrom

    if (con%c(jsp)%variable == chrom$) then
      call next_row (jsp)
      b(n_row) = con%c(jsp)%target_value
      if (con%c(jsp)%plane == x_plane$) then
        if (con%compensate_natural) b(n_row) = b(n_row) - momnt%chrom_x_0
        call row_stuffit (row(n_row, :), momnt%chrom_x)
      else
        if (con%compensate_natural) b(n_row) = b(n_row) - momnt%chrom_y_0
        call row_stuffit (row(n_row, :), momnt%chrom_y)
      endif
    endif

! tonality

    if (con%c(jsp)%variable == tonality$) then
      call next_row (jsp)
      b(n_row) = con%c(jsp)%target_value / 390.1  ! convert from kHz
      if (con%c(jsp)%plane == x_plane$) then
        call row_stuffit (row(n_row, :), momnt%tonality_x)
      else
        call row_stuffit (row(n_row, :), momnt%tonality_y)
      endif
    endif

  enddo  ! jsp

! coupling_a_real

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == coupling_a_real$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%a_real(j,:))
    endif
  enddo

! coupling_a_image

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == coupling_a_image$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%a_image(j,:))
    endif
  enddo

! coupling_b_real

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == coupling_b_real$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%b_real(j,:))
    endif
  enddo

! coupling_b_image

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == coupling_b_image$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%b_image(j,:))
    endif
  enddo

! x plane delta_beta

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == delta_beta$ &
                        .and. con%c(jsp)%plane == x_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), momnt%dbeta_x(j,:)*con%delta_e)
      if (con%compensate_natural) b(n_row) = -momnt%dbeta_x_0(j) * con%delta_e
    endif
  enddo

! y plane delta_beta

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == delta_beta$ &
                        .and. con%c(jsp)%plane == y_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), momnt%dbeta_y(j,:)*con%delta_e)
      if (con%compensate_natural) b(n_row) = -momnt%dbeta_y_0(j) * con%delta_e
    endif
  enddo

! x plane dbeta_dpretz

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dpretz$ &
                        .and. con%c(jsp)%plane == x_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), momnt%dbeta_dpretz_x(j,:))
       b(n_row) = con%c(ix_ele_sp(j))%target_value
      call next_row (ix_ele_sp(j))   ! this is for electrons
      call row_stuffit (row(n_row, :), -1*momnt%dbeta_dpretz_x(j,:))
       b(n_row) = -con%c(ix_ele_sp(j))%target_value
    endif
  enddo

! y plane dbeta_dpretz

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dpretz$ &
                        .and. con%c(jsp)%plane == y_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), momnt%dbeta_dpretz_y(j,:))
       b(n_row) = con%c(ix_ele_sp(j))%target_value
      call next_row (ix_ele_sp(j))          ! this is for electrons
      call row_stuffit (row(n_row, :), -1 * momnt%dbeta_dpretz_y(j,:))
       b(n_row) = -con%c(ix_ele_sp(j))%target_value
    endif
  enddo

! x plane dbeta_dcos

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dcos$ &
                        .and. con%c(jsp)%plane == x_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%dbeta_dcos_x(j,:))
    endif
  enddo

! y plane dbeta_dcos

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dcos$ &
                        .and. con%c(jsp)%plane == y_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%dbeta_dcos_y(j,:))
    endif
  enddo

! x plane dbeta_dsin

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dsin$ &
                        .and. con%c(jsp)%plane == x_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%dbeta_dsin_x(j,:))
    endif
  enddo

! y plane dbeta_dsin

  ix_ele_sp = 0
  do jsp = 1, con%n_constraint
   if (con%c(jsp)%weight == 0) cycle
    if (con%c(jsp)%variable == dbeta_dsin$ &
                        .and. con%c(jsp)%plane == y_plane$) then
      do j = con%c(jsp)%where1_ix, con%c(jsp)%where2_ix
        ix_ele_sp(j) = jsp
      end do
    endif
  enddo
  do j = 1, lat%n_ele_track
    if (useable_ele(j) .and. ix_ele_sp(j) /= 0) then
      call next_row (ix_ele_sp(j))
      call row_stuffit (row(n_row, :), -momnt%dbeta_dsin_y(j,:))
    endif
  enddo

  if (n_row > 8000) then
    print *,' ERROR IN MINIMIZE_MOMENTS: NUMBER OF ROW GREATER THAN 8000'
    call err_exit
  endif

!enddo electron_positron_loop

!--------------------------------------------------------------
! find optimal solution

  c(1:n_row) = matmul(row(1:n_row, 1:n_var), momnt%k2l(1:n_var))
  c = (c - b) * weight
  chi2 = dot_product(c(1:n_row), c(1:n_row))
  if (k2_4_row /= 0) chi2 = chi2 + weight(k2_4_row)**2 * sum(momnt%k2l**4)

! find better sextupole distribution

  if (con%n_loops >= 0) then

    if (k2_4_row /= 0) then

    sig(1:n_row) = 1 / weight(1:n_row)
    a_lambda = -1
    old_chi2 = chi2

    nn = n_row
    do i = 1, 10
      if (i > 3 .and. chi2 > 0.95*old_chi2) a_lambda = 0
      call mrqmin (x(1:nn), b(1:nn), sig(1:nn), momnt%k2l, momnt%use, &
                             covar, alpha, chi2, eval_sex_moments, a_lambda)
      if (a_lambda == 0) exit
      old_chi2 = chi2
    enddo

    else  ! k2^4 not used

      nv_use = 0

      fixed_sext = 0
      where (.not. momnt%use) fixed_sext = momnt%k2l

    if(con%compensate_natural)then
      c(1:n_row) = matmul(row(1:n_row, 1:n_var), fixed_sext(1:n_var))
    else
      c(1:n_row) = 0.
    endif

      bb(:) = (b(:)-c(:)) * weight(:)

      do i=1, n_var
        if (momnt%use(i)) then
          nv_use = nv_use+1
          use_to_i(nv_use) = i
          r0(:,nv_use) = row(:,i) * weight(:)
        endif
      end do

      nn = n_row; nv = nv_use
      call svdcmp (R0(1:nn,1:nv), w(1:nv), v(1:nv,1:nv))
      call svbksb (R0(1:nn,1:nv), w(1:nv), v(1:nv,1:nv), bb(1:nn), x(1:nv))

      forall (j = 1:nv_use) momnt%k2l(use_to_i(j)) = x(j)
!dlr 4/12/03      forall (j = 1:nv_use) momnt_positron%k2l(use_to_i(j)) = x(j)
!dlr 4/12/03      forall (j = 1:nv_use) momnt_electron%k2l(use_to_i(j)) = x(j)

    endif ! k2^4 used?

! update new results

    do i = 1, n_var
      if (.not. momnt%use(i) .and. .not. con%compensate_natural) &
                                                          momnt%k2l(i)=0.
      if (momnt%use(i)) then
        ix = momnt%ix(i)
        if (con%compensate_natural) then
          lat%ele(ix)%value(k2$) = momnt%k2l(i)/lat%ele(ix)%value(l$)
        else
          lat%ele(ix)%value(k2$) = momnt%k2l(i)/lat%ele(ix)%value(l$)  &
                        + lat%ele(ix)%value(k2$)
        endif

        if (eastwest_sym) then
          if (west_to_east(i) == 0) then
            print *, 'ERROR IN MINIMIZE_MOMENTS: CANNOT FIND CORRESPONDING'
            print *, '      EAST SEXTUPOLE FOR: ', momnt%name(i)
            call err_exit
          endif
          momnt%k2l(west_to_east(i)) = momnt%k2l(i)
          ix = west_to_east(i)
          if (con%compensate_natural) then
            lat%ele(ix)%value(k2$) = momnt%k2l(i)/lat%ele(ix)%value(l$)
          else
            lat%ele(ix)%value(k2$) = momnt%k2l(i)/lat%ele(ix)%value(l$)  &
                          + lat%ele(ix)%value(k2$)
          endif
        endif

!        lat_pos%ele(ix)%value(k2$) = lat%ele(ix)%value(k2$)
!        lat_ele%ele(ix)%value(k2$) = lat%ele(ix)%value(k2$)
      endif
    end do

  endif  ! n_loops >= 0

!----------------------------------------------------------------------------
! write to file

  open(unit = 41, file = 'opt_sextupoles.dat', status='REPLACE')

  fmt1 = "(a7,'[K2] = ',e12.6,'; ',a7,'[K2] = ',e12.6)"
  fmt2 = "(1x,a7,'/',e10.4,',',a7,'/'e10.4,', &')"

  do i=1,momnt%n_sext
    if (momnt%name(i)(7:7) == 'W') then
      if (west_to_east(i) /= 0) then
        l1 = momnt%k2l(i) / lat%ele(momnt%ix(i))%value(l$)
        j = west_to_east(i)
        l2 = momnt%k2l(j) / lat%ele(momnt%ix(j))%value(l$)
        write(41, fmt1) momnt%name(i), l1, momnt%name(j), l2
        if(con%compensate_natural)then
            l1 = l1 - lat_old%ele(momnt%ix(i))%value(k2$)
            l2 = l2 - lat_old%ele(momnt%ix(j))%value(k2$)
        endif
        write(42, fmt2) momnt%name(i), l1, momnt%name(j), l2
      endif
    endif
  end do

  close (unit=41)

! calculate new chi2

  c(1:n_row) = matmul(row(1:n_row, 1:n_var), momnt%k2l(1:n_var))
  c = (c - b) * weight
  momnt%chi2 = dot_product(c, c)

  forall (i = 1:n_row) con%c(ix_sp(i))%contribution = 0
  if (con%use_sextupole_fom) then
    do i = 1, n_row
      j = ix_sp(i)
      con%c(j)%contribution = con%c(j)%contribution + c(i)**2
    enddo
    if (k2_4_row /= 0) con%c(ix_sp(k2_4_row))%contribution = &
                                       weight(k2_4_row)**2 * sum(momnt%k2l**4)
  endif

  do i = 1, n_row
    j = ix_sp(i)
    if (con%c(j)%weight /= 0) con%c(j)%actual_value = &
                      sqrt(con%c(j)%contribution/con%c(j)%weight)
  enddo

  return

!------------------------------------------------------------------------

contains

!------------------------------------------------------------------------

subroutine row_stuffit (arow, mom_row)

    implicit none

    real(rp) arow(*), mom_row(*)

    do i = 1, n_sext
      if (eastwest_sym) then
        if (west_to_east(i) == 0) exit
        i2 = momnt%ix(west_to_east(i))
        arow(i) = mom_row(i) + mom_row(i2)
      else
        arow(i) = mom_row(i)
      endif
    enddo

  end subroutine

!------------------------------------------------------------------------

subroutine next_row (j_sp)

    implicit none

    integer j_sp

    n_row = n_row + 1
    weight(n_row) = con%c(j_sp)%weight
    ix_sp(n_row) = j_sp

  end subroutine

end subroutine

!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------

subroutine eval_sex_moments (x, a, yfit, dyda)

  use moments_com_mod
  implicit none
  real(rp) :: x(:), a(:), yfit(:), dyda(:, :)
  integer na, nx
!
  na = ubound(a, 1)
  nx = ubound(x, 1)
  yfit = matmul(row(1:nx,1:na), a)
  yfit(k2_4_row) = sum(a**4)
  dyda = row(1:nx, 1:na)
  dyda(k2_4_row,:) = 4 * a**3

end subroutine