module correct_ring_mod2 use super_recipes_mod use bmad use sim_bpm_mod use dr_misalign_mod ! use ele_noise_mod use tune_tracker_mod use sim_ac_and_cbar_mod use sim_eta_mod use bsim_cesr_interface implicit none !========================================================================== ! Define some logicals and structures to be used here and elsewhere ! Measurement parameters integer, parameter :: orbit_x$ = 1, orbit_y$ = 2, eta_xx$ = 3, eta_yy$ = 4 integer, parameter :: phi_aa$ = 5, phi_bb$ = 6, cbar22$ = 7, cbar12$ = 8, cbar11$ = 9 character(40), parameter :: allowed_meas_names(9) = (/'orbit_x', 'orbit_y', 'eta_xx ', 'eta_yy ', 'phi_aa ', 'phi_bb ', & 'cbar22 ', 'cbar12 ', 'cbar11 ' /) type det_grp character(40) mask character(40) :: param_name = "" character(40) :: key_name integer :: key integer :: param = 0 real(rp) wt end type det_grp type cor_grp character(40) mask character(40) :: attrib_name = "" character(40) :: key_name integer :: key real(rp) wt end type cor_grp integer, parameter :: n_detcor_groups = 10 type correct_struct character(40) :: bpm_mask = '' type(det_grp) :: det(n_detcor_groups) type(cor_grp) :: cor(n_detcor_groups) end type correct_struct type correct_ring_params_struct real(rp) :: det_abs_res, det_diff_res, det_rot_res, n_lm_iterations, & eta_delta_e_e, sigma_cutoff logical :: write_elements = .true., skip_dup_ele end type correct_ring_params_struct type(correct_ring_params_struct) :: correct_ring_params type meas_struct character(40) :: name = "" character(40) :: param_name = "" integer :: loc = 0, param = 0 real(rp) :: wt = 0 !, rotation = 0, x_offset = 0, y_offset = 0 logical :: is_on = .true. end type meas_struct type parameter_ele_struct character(40) :: name = "" character(40) :: attrib_name = "" integer :: loc = 0 type(all_pointer_struct) :: attrib_ptr_model, attrib_ptr_meas real(rp) :: orig_val_model = 0. real(rp) :: wt = 0 end type parameter_ele_struct type measurement_struct real(rp) :: orbit_x=0, orbit_y=0, eta_x=0, eta_y=0, phi_a=0, phi_b=0 real(rp) :: cbar22=0, cbar12=0, cbar11=0 end type measurement_struct integer, parameter :: n_det_max = 1000 type(lat_struct) :: cr_model_ring ! defined / allocated in ring_ma2.f90, prior to calling this module type(coord_struct), allocatable :: cr_model_co(:) type(parameter_ele_struct) :: p_ele(n_det_max) contains !========================================================================== ! Routine to correct a misaligned ring subroutine correct_ring(ring, tt_params, id, bpm, correct, rms_param, n_turns_phase, n_turns_orbit, n_damping_phase, n_damping_orbit, init_vec, eta_method, verbose, meas_type, meas, seed, write_all_lats, status) implicit none type(lat_struct) :: ring type(correct_struct), target :: correct type(det_grp), pointer :: det_grp_ptr type(cor_grp), pointer :: cor_grp_ptr type(all_pointer_struct) :: attrib_ptr_model, attrib_ptr_meas real(rp), optional, intent(out) :: rms_param type(det_struct), target :: bpm(:) type(det_data_struct), allocatable :: bpm_save(:,:) type(tt_param_struct), target :: tt_params(:) integer, intent(in) :: n_turns_phase, n_turns_orbit, n_damping_phase, n_damping_orbit character(6), intent(in) :: eta_method INTEGER :: id(max_tt), seed real(rp) :: init_vec(:), orig_val_model logical err, err1, err2, err3 logical verbose, write_all_lats ! type(ele_noise_struct), optional :: ele_noise(:) type (super_mrqmin_storage_struct) storage type(all_pointer_struct) :: calib_ptr type (meas_struct) det(n_det_max) type (measurement_struct), allocatable :: meas(:) integer n_p_ele, n_det, n_det_use, loc, ix, jx type(coord_struct), allocatable :: co(:), eta(:) type(ele_struct), pointer :: ele, ele_model integer i_ele, i_det_grp, i_cor_grp, i_det, i, i_ele2 integer size_y, iy, iter, ip, rad_cache integer status, tt_stat real(rp) harvest(7), cbar(2,2), chisq, alamda, old_chisq, d_chisq real(rp), allocatable, dimension(:) :: x, y, wgt, a logical, allocatable, dimension(:) :: maska character(40) attrib_name, lat_type ! added 2011.07.06 -JSh INTEGER :: n_TTs = 2 ! number of tune trackers integer :: bpm_ix = 0, dat_ix = 0, bpm_ix_old = 0 real(rp) :: current = 1000000 integer :: n_bpms = 0 ! parameters for ac_eta calculation type (cesr_all_data_struct) ac_eta real(rp) :: rms_eta_x = 0, rms_eta_y = 0. logical :: use_new_method = .false. character(20) :: meas_type real(rp), allocatable :: phi(:,:,:), on_off_save(:) real(rp), allocatable :: cbar_meas(:,:,:), cbar_bmad(:,:,:) ! (:,2,2) = across all BPMs real(rp), allocatable :: Amp(:,:,:) ! = (n_bpms,2,2) = (i_bpm, (Aa,Ab),(Ba,Bb)) type(normal_modes_struct) modes integer ix_db, ix_lat ! debug pointers; workaround for a Totalview bug: type(tt_param_struct), pointer :: tt_ptr(:) type(det_struct), pointer :: bpm_ptr(:) if (.not. allocated(meas)) allocate(meas(n_det_max)) if (trim(eta_method) .eq. 'ac_old') then use_new_method = .false. elseif (trim(eta_method) .eq. 'ac_new') then use_new_method = .true. elseif (trim(eta_method) .ne. 'dc') then write(*,*) "eta_method ", trim(eta_method), " not a valid option. Must be 'ac_old', 'ac_new', or 'dc'." stop end if ! debug pointers; a workaround for Totalview bug: tt_ptr => tt_params bpm_ptr => bpm !-------------------------------------------- ! Locate detectors n_det = 0 n_det_use = 0 do i_det_grp = 1, n_detcor_groups det_grp_ptr => correct%det(i_det_grp) ! automatically cast to uppercase call upcase_string(det_grp_ptr%param_name) call upcase_string(det_grp_ptr%key_name) call upcase_string(det_grp_ptr%mask) if (len_trim(det_grp_ptr%param_name) == 0) cycle ! no name provided ! look up index of desired parameter call match_word(det_grp_ptr%param_name, allowed_meas_names, ix) det_grp_ptr%param = ix if (det_grp_ptr%param .eq. 0) then write(*,*) "Parameter ", trim(det_grp_ptr%param_name), & " not found for detectors. Cycling." cycle endif ele_loop: do i_det = 1, size(bpm) ele => ring%ele(bpm(i_det)%ix_lat) n_det = n_det + 1 if (ring%ele(bpm(i_det)%ix_lat)%is_on) n_det_use = n_det_use + 1 ! can't say det(i_det)%blah because n_det > size(bpm) for multiple correction iterations ! note: det_grp_ptr does not point to det, but rather correct%det(...) det(n_det)%name = ele%name det(n_det)%loc = bpm(i_det)%ix_lat det(n_det)%param = det_grp_ptr%param det(n_det)%param_name = det_grp_ptr%param_name det(n_det)%wt = det_grp_ptr%wt det(n_det)%is_on = ring%ele(bpm(i_det)%ix_lat)%is_on end do ele_loop end do !------------------------------------------------------- ! Locate elements to be varied for a particular correction n_p_ele = 0 do i_cor_grp = 1, n_detcor_groups cor_grp_ptr => correct%cor(i_cor_grp) call upcase_string(cor_grp_ptr%attrib_name) call upcase_string(cor_grp_ptr%key_name) call upcase_string(cor_grp_ptr%mask) if (len(trim(cor_grp_ptr%attrib_name)) == 0) cycle ! no ele specified; cycle cor_grp_ptr%key = key_name_to_key_index(cor_grp_ptr%key_name, .true.) ele_loop2: do i_ele = 1, ring%n_ele_max ! steerings are overlays, not tracked elements ele => ring%ele(i_ele) ele_model => cr_model_ring%ele(i_ele) ! check a few things first: ! regex match name: if (.not. (match_reg(trim(ele%name), trim(cor_grp_ptr%mask)))) cycle ele_loop2 ! match element key to desired key: if (.not. ((ele%key .eq. cor_grp_ptr%key) .or. (cor_grp_ptr%key .eq. -1))) cycle ele_loop2 ! is not a slave: if (ele%slave_status .eq. super_slave$) cycle ele_loop2 ! is free to vary: if (.not. attribute_free(ele, cor_grp_ptr%attrib_name, err_print_flag = .false.)) cycle ele_loop2 ! get pointers to attributes for each individual element; note we need to do this for both 'ring' and 'cr_model_ring' call pointer_to_attribute(ele, cor_grp_ptr%attrib_name, .true., attrib_ptr_meas, err) call pointer_to_attribute(ele_model, cor_grp_ptr%attrib_name, .true., attrib_ptr_model, err) n_p_ele = n_p_ele + 1 p_ele(n_p_ele)%name = ele%name p_ele(n_p_ele)%loc = i_ele p_ele(n_p_ele)%attrib_ptr_meas = attrib_ptr_meas p_ele(n_p_ele)%attrib_ptr_model = attrib_ptr_model p_ele(n_p_ele)%attrib_name = cor_grp_ptr%attrib_name p_ele(n_p_ele)%wt = cor_grp_ptr%wt p_ele(n_p_ele)%orig_val_model = value_of_attribute(ele_model, cor_grp_ptr%attrib_name) ! NEW form; %param is now a series of pointers "attrib_ptr_*" end do ele_loop2 ! Check for too many parameter elements if (n_p_ele == size(p_ele)) then write(*,*) "Too many parameter elements. Increase hard-coded size(p_ele)." call err_exit end if end do !========================================================================== ! Write out list of detectors and parameter elements if (correct_ring_params%write_elements) then write(*,*) " " write(*,'(a,i0,a,i0,a)') "Found ", n_det, " detectors. Using ", n_det_use, '.' if (verbose .eqv. .true.) then do i_ele = 1, n_det write(*,'(a10,a,a10,g10.2)', advance='no') & trim(det(i_ele)%name), ' ', & det(i_ele)%param_name, det(i_ele)%wt if (mod(i_ele,3) == 0) then write(*,*) else write(*,'(a3)', advance='no') " | " end if end do write(*,*) endif ! if verbose write(*,'(a,i0,a)') "Found ", n_p_ele, " corrector elements." if (verbose .eqv. .true.) then do i_ele = 1, n_p_ele write(*,'(a10,a,a10,g10.2)', advance='no') trim(p_ele(i_ele)%name), ' ',& p_ele(i_ele)%attrib_name, p_ele(i_ele)%wt if (mod(i_ele,3) == 0) then write(*,*) else write(*,'(a3)', advance='no') " | " end if end do write(*,*) else ! if NOT verbose write(*,*) "*************Verbose mode disabled!*****************" write(*,*) "To see full list of correctors, set verbose = .true." endif write(*,*) " " end if !========================================================================== ! Store simulated measurements with BPM errors ! Rewritten 7/6/11 to populate data structures with simulated measurements. -JSh !========================================================================== ! recall that BPM errors are populated at the same level as dr_misalign n_bpms = ubound(bpm,1) allocate(bpm_save(n_turns_orbit, n_bpms)) call set_on_off(rfcavity$, ring, save_state$, saved_values = on_off_save) call set_on_off(rfcavity$, ring, off$) ! call twiss_and_track(ring, co, status) call closed_orbit_calc(ring, co, 4, err_flag = err1) call lat_make_mat6(ring, -1, co, err_flag = err2) call twiss_at_start(ring, status = status) call twiss_propagate_all(ring, err_flag = err3) call set_on_off(rfcavity$, ring, restore_state$, saved_values = on_off_save) if (trim(meas_type) .eq. 'tracking') then call twiss3_at_start(ring, err) call twiss3_propagate_all(ring) ! init_vec has already been checked in main ring_ma2 routine ! reset each tune tracker now, after computing closed orbit ! for this specific correction on this misaligned lattice. tt_params(1)%offset = co(tt_params(1)%bpm_loc)%vec(1) tt_params(2)%offset = co(tt_params(2)%bpm_loc)%vec(3) do ix=1, n_TTs call reset_dTT(id(ix), tt_params(ix)) enddo ! simulate the data: write(*,*) "Simulating orbit measurement..." ! if (present(ele_noise)) then ! call sim_tbt_data(ring, bpm, n_turns_orbit, .true., 0, tt_params, id, n_damping_orbit, init_vec, bpm_save, ele_noise = ele_noise) ! disable tunetrackers for orbit measurement ! else call sim_tbt_data(ring, bpm, n_turns_orbit, .true., 0, tt_params, id, n_damping_orbit, init_vec, bpm_save) ! end if write(*,*) " " if (trim(eta_method) .ne. 'dc') then write(*,*) "Simulating phase/cbar and ac_eta measurements..." ! if (present(ele_noise)) then ! call sim_ac_and_cbar_data(ring, bpm, n_turns_phase, 3, tt_params, id, n_damping_phase, & ! init_vec, use_new_method, ac_eta, phi, cbar_meas, Amp, ele_noise) ! else call sim_ac_and_cbar_data(ring, bpm, n_turns_phase, 3, tt_params, id, n_damping_phase, & init_vec, use_new_method, ac_eta, phi, cbar_meas, Amp) ! end if else ! use DC dispersion write(*,*) "Simulating phase/cbar measurements..." ! if (present(ele_noise)) then ! call sim_cbar_data(ring, bpm, n_turns_phase, .true., 2, tt_params, id, n_damping_phase, & ! init_vec, phi, cbar_meas, Amp, ele_noise) ! else call sim_cbar_data(ring, bpm, n_turns_phase, .true., 2, tt_params, id, n_damping_phase, & init_vec, phi, cbar_meas, Amp) ! end if write(*,*) "" write(*,*) "Simulating DC dispersion..." ! if (present(ele_noise)) then ! call sim_dc_eta_data(ring, bpm, n_turns_orbit, & ! correct_ring_params%eta_delta_e_e, eta, ele_noise) ! else call sim_dc_eta_data(ring, bpm, n_turns_orbit, & correct_ring_params%eta_delta_e_e, eta) ! end if endif else ! must be meas_type = 'basic'; already filtered out bad options at startup call sim_eta_data(ring, bpm, correct_ring_params%eta_delta_e_e, eta) if (.not. allocated(cbar_bmad)) allocate(cbar_bmad(1:n_bpms,2,2), cbar_meas(1:n_bpms,2,2)) ! orbit and cbar: bpm(:)%vec(1) = co(bpm(:)%ix_lat)%vec(1) bpm(:)%vec(3) = co(bpm(:)%ix_lat)%vec(3) do jx = 1, n_bpms ele => ring%ele(bpm(jx)%ix_lat) call apply_bpm_errors(bpm(jx), current, lat_type) ! note: also transfers vec(6)-->amp(4) button signals ! transfer back to co for now: co(bpm(jx)%ix_lat)%vec(:) = bpm(jx)%vec ! note: this also populates the cbar_meas struct-- don't comment out! ! note: rotate_cbar presently just populates cbar_meas; don't trust calculation --2016.05.24 call rotate_cbar(ele, bpm(jx)%tilt, cbar_meas(jx,:,:)) enddo endif ! allocate and populate the 'bmad' (read: design) values if (.not. allocated(cbar_bmad)) allocate(cbar_bmad(n_bpms,2,2)) do ix = 1, n_bpms call c_to_cbar(ring%ele(bpm(ix)%ix_lat),cbar_bmad(ix,:,:)) enddo ! load data into 'meas' data structure: do dat_ix = 1, n_det loc = det(dat_ix)%loc ele => ring%ele(loc) ! there might be an easier way to align bpm(:) and det(:)... do ix = 1, size(bpm) if (bpm(ix)%ix_lat .eq. loc) then bpm_ix = ix exit endif enddo if (trim(meas_type) .eq. 'tracking') then ! new values: meas(dat_ix)%orbit_x = sum(bpm_save(:,bpm_ix)%vec(1))/size(bpm_save(:,bpm_ix)) ! average orbit over several turns meas(dat_ix)%orbit_y = sum(bpm_save(:,bpm_ix)%vec(3))/size(bpm_save(:,bpm_ix)) if (trim(eta_method) .ne. 'dc') then meas(dat_ix)%eta_x = ac_eta%ac_eta_x(bpm_ix)%value meas(dat_ix)%eta_y = ac_eta%ac_eta_y(bpm_ix)%value else meas(dat_ix)%eta_x = eta(loc)%vec(1) meas(dat_ix)%eta_y = eta(loc)%vec(3) endif meas(dat_ix)%phi_a = phi(bpm_ix,1,1) ! already phase-offset-corrected meas(dat_ix)%phi_b = phi(bpm_ix,2,2) meas(dat_ix)%cbar22 = cbar_meas(bpm_ix,2,2) meas(dat_ix)%cbar12 = cbar_meas(bpm_ix,1,2) meas(dat_ix)%cbar11 = cbar_meas(bpm_ix,1,1) else ! must be 'basic': meas(dat_ix)%orbit_x = co(loc)%vec(1) meas(dat_ix)%orbit_y = co(loc)%vec(3) meas(dat_ix)%eta_x = eta(loc)%vec(1) meas(dat_ix)%eta_y = eta(loc)%vec(3) meas(dat_ix)%phi_a = ele%a%phi meas(dat_ix)%phi_b = ele%b%phi meas(dat_ix)%cbar22 = cbar_meas(bpm_ix,2,2) meas(dat_ix)%cbar12 = cbar_meas(bpm_ix,1,2) meas(dat_ix)%cbar11 = cbar_meas(bpm_ix,1,1) endif enddo !========================================================================== ! Assemble the "target" values for the fitter. These are just the parameter ! element values followed by the simulated measurements. size_y = n_det_use + n_p_ele allocate(x(size_y), y(size_y), wgt(size_y)) allocate(a(n_p_ele), maska(n_p_ele)) iy = 1 y = 0. a = 0. ! Load element values do i_ele = 1, n_p_ele x(iy) = iy if (p_ele(i_ele)%wt .lt. 1.e-12) then ! floating-point comparison wgt(iy) = 0. else wgt(iy) = p_ele(i_ele)%wt endif maska(iy) = .true. iy = iy + 1 end do ! Load measurements do i_ele = 1, n_det if (.not. det(i_ele)%is_on) cycle x(iy) = i_ele if (det(i_ele)%wt .lt. 1.e-12) then ! floating-point comparison wgt(iy) = 0. else wgt(iy) = det(i_ele)%wt endif select case (det(i_ele)%param) case(orbit_x$) y(iy) = meas(i_ele)%orbit_x case(orbit_y$) y(iy) = meas(i_ele)%orbit_y case(eta_xx$) y(iy) = meas(i_ele)%eta_x case(eta_yy$) y(iy) = meas(i_ele)%eta_y case(phi_aa$) y(iy) = meas(i_ele)%phi_a case(phi_bb$) y(iy) = meas(i_ele)%phi_b case(cbar22$) y(iy) = meas(i_ele)%cbar22 case(cbar12$) y(iy) = meas(i_ele)%cbar12 case(cbar11$) y(iy) = meas(i_ele)%cbar11 case default write(*,*) "Invalid detector key" call err_exit end select iy = iy + 1 end do !========================================================================== ! Make specified number of calls to MRQMIN, then cleanup. Run at least ! n_lm_iterations times, then stop as soon as chisq doesn't change. write(*,*) "Calculating corrections..." alamda = -1. write (*,'(A8,3A12)') "Iter", "Chisq", "D_Chisq", "A_lambda" chisq = dot_product(y**2, wgt) old_chisq = chisq ! print starting chisq: write(*,'(i8,es12.4)') 0, chisq do iter =1, correct_ring_params%n_lm_iterations+1 if (iter == correct_ring_params%n_lm_iterations+1) alamda = 0 ! tell mrqmin we're done call super_mrqmin(y,wgt,a,chisq,lmfunc,storage,alamda, status, maska) if (status /= 0) then print*, "*** Singular matrix encountered--nothing more to do with this seed." goto 999 end if d_chisq = (chisq-old_chisq) write(*,'(i8,3es12.4)') iter, chisq, d_chisq, alamda old_chisq = chisq end do !call super_mrqmin(y,wgt,a,chisq,lmfunc,storage,alamda,status,maska) !========================================================================== ! Apply correction to misaligned ring, i.e., put in the negative of the ! parameter values calculated by the minimization. write(*,*) "Loading correction..." do i_ele = 1, n_p_ele ele => ring%ele(p_ele(i_ele)%loc) ! measured ring p_ele(i_ele)%attrib_ptr_meas%r = p_ele(i_ele)%attrib_ptr_meas%r - a(i_ele)*value_of_attribute(ele, 'CALIB', err) call set_flags_for_changed_attribute(ele, p_ele(i_ele)%attrib_ptr_meas%r) end do call lattice_bookkeeper(ring) call set_on_off(rfcavity$, ring, save_state$, saved_values = on_off_save) call set_on_off(rfcavity$, ring, off$) ! call twiss_and_track(ring, co, status) call closed_orbit_calc(ring, co, 4, err_flag = err1) call lat_make_mat6(ring, -1, co, err_flag = err2) call twiss_at_start(ring, status = status) call twiss_propagate_all(ring, err_flag = err3) call set_on_off(rfcavity$, ring, restore_state$, saved_values = on_off_save) ! write corrections to file now: if (write_all_lats) then call write_corrections(ring, p_ele, seed) else call write_corrections(ring, p_ele, -1) endif write(*,*) "Correction loaded." write(*,*) " " ! Clean up deallocate(x, y, wgt, a, maska, co) status = 0 999 return contains !========================================================================== ! Routine to provide the values for the chi^2 that MRQMIN will attempt to ! minimize. This is essentially just a wrapper for passing values to ! RING_TO_Y. subroutine lmfunc(a,y,dyda, status) implicit none real(rp), dimension(:), intent(in) :: a real(rp), dimension(:), intent(out) :: y real(rp), dimension(:,:), intent(out) :: dyda real(rp), dimension(size(a)) :: a2 real(rp), dimension(size(y)) :: y1, y2 real(rp), parameter :: delta = 1.e-7 integer status, bmad_ok integer ia call ring_to_y(a, y) if (any(y==999)) then dyda = 999 return end if a2 = a do ia = 1, size(a) a2(ia) = a(ia) + delta call ring_to_y(a2, y2) ! Central difference ! a2(ia) = a(ia) - delta ! call ring_to_y(a2, y1) ! dyda(:,ia) = (y2-y1)/(2.*delta) ! Forward difference dyda(:,ia) = (y2-y)/delta a2(ia) = a(ia) end do end subroutine lmfunc !========================================================================== ! Routine used in minimization for calculating how the relevant output ! parameters (Y) change as a function of the input parameters (A). subroutine ring_to_y(a, y) implicit none real(rp), dimension(:), intent(in) :: a real(rp), dimension(:), intent(out) :: y type(ele_struct), pointer :: ele integer :: ia, iy, ip, i_det, this_ia, bmad_ok, status real(rp), allocatable :: cbar_model(:,:,:), orb_model(:,:), eta_model(:,:), on_off_save(:) ! allocate and initialize orb_model, eta_model if (.not. allocated(orb_model)) allocate(orb_model(n_det,6), eta_model(n_det,2), cbar_model(n_det,2,2)) orb_model = 0 eta_model = 0 cbar_model = 0 ! Put the values A into the right places in the ring do ia = 1, size(a) ele => cr_model_ring%ele(p_ele(ia)%loc) p_ele(ia)%attrib_ptr_model%r = p_ele(ia)%orig_val_model + a(ia) ! note the omission of 'calib'; fitter should not know about the calibration error call set_flags_for_changed_attribute(ele, p_ele(ia)%attrib_ptr_model%r) end do call lattice_bookkeeper(cr_model_ring) ! Recalculate the ring parameters if (allocated(cr_model_co)) cr_model_co(0)%vec = 0. call set_on_off(rfcavity$, cr_model_ring, save_state$, saved_values = on_off_save) call set_on_off(rfcavity$, cr_model_ring, off$) ! call twiss_and_track(cr_model_ring, cr_model_co, status) call closed_orbit_calc(cr_model_ring, cr_model_co, 4, err_flag = err1) call lat_make_mat6(cr_model_ring, -1, cr_model_co, err_flag = err2) call twiss_at_start(cr_model_ring, status = status) call twiss_propagate_all(cr_model_ring, err_flag = err3) call set_on_off(rfcavity$, cr_model_ring, restore_state$, saved_values = on_off_save) if ((err1 .eqv. .false.) .and. (err2 .eqv. .false.) .and. (err3 .eqv. .false.) .and. (status == ok$)) then status = 1 else y = 999 return endif ! populate the orbit and dispersion measurement structures ! in order to allow for fitted BPM tilts to be applied do i_det = 1, n_det !if (.not. det(i_det)%is_on) cycle loc = det(i_det)%loc ele => cr_model_ring%ele(loc) orb_model(i_det, :) = cr_model_co(loc)%vec(:) eta_model(i_det, 1) = ele%x%eta eta_model(i_det, 2) = ele%y%eta call c_to_cbar(ele, cbar_model(i_det,:,:)) enddo ! ! special case: if BPM misalignments are being used as fit parameters, ! ! apply them here. Start by cycling over all detectors ! do i_det = 1, n_det ! ! ele => cr_model_ring%ele(det(i_det)%loc) ! ! ! align whether a marker is used both for fitting and measurement: ! this_ia = -1 ! do ia = 1, size(a) ! if (p_ele(ia)%loc .eq. ele%ix_ele) then ! matched element ! this_ia = ia ! exit ! endif ! enddo ! if (this_ia .eq. -1) cycle ! didn't match this detector to a varied parameter in p_ele ! ! ! if this element (which is also a detector) was fit using a tilt, rotate it: ! if (p_ele(ia)%attrib_name .eq. 'TILT') then ! call rotate_meas(orb_model(i_det,1),orb_model(i_det,3),a(ia)) ! orbit ! call rotate_meas(eta_model(i_det,1),eta_model(i_det,2),a(ia)) ! dispersion ! call rotate_cbar(ele, a(ia), cbar_model(i_det,:,:)) ! coupling (cbar) ! endif ! ! enddo if (status /= ok$) then y = 999 return end if ! Stick everything back into Y iy = 1 ! Load element values do ip = 1, n_p_ele y(iy) = a(ip) iy = iy + 1 end do ! Load measurement information do i_det = 1, n_det if (.not. det(i_det)%is_on) cycle loc = det(i_det)%loc ele => cr_model_ring%ele(loc) select case (det(i_det)%param) case(orbit_x$) y(iy) = orb_model(i_det, 1) !y(iy) = cr_model_co(loc)%vec(1) case(orbit_y$) y(iy) = orb_model(i_det, 3) !y(iy) = cr_model_co(loc)%vec(3) case(eta_xx$) y(iy) = eta_model(i_det, 1) !y(iy) = ele%x%eta case(eta_yy$) y(iy) = eta_model(i_det, 2) !y(iy) = ele%y%eta case(phi_aa$) y(iy) = ele%a%phi case(phi_bb$) y(iy) = ele%b%phi case(cbar22$) y(iy) = cbar_model(i_det,2,2) case(cbar12$) y(iy) = cbar_model(i_det,1,2) case(cbar11$) y(iy) = cbar_model(i_det,1,1) case default write(*,*) "Invalid detector key" call err_exit end select iy = iy + 1 end do ! Send this back if it's wanted if (present(rms_param)) rms_param = sqrt(dot_product(a, a)/size(a)) end subroutine ring_to_y end subroutine correct_ring !========================================================================== ! UPDATED routine to write misaligned lattice to a file subroutine write_misalignments(ma_ring, ma_params, seed) implicit none type(lat_struct), intent(in), target :: ma_ring type(ma_struct), intent(in), target :: ma_params(:) type(ma_struct), pointer :: ma type(ele_struct), pointer :: ele type (all_pointer_struct) a_ptr character(40) :: write_name = '', filename = '' integer i_param, i_ele, i_cor, ix, jx integer, intent(in) :: seed logical err if (seed > 0) then write(filename,'(a,i9.9,a)') 'ma_lat.',seed,'.bmad' else filename = 'ma_lat.bmad' endif open(unit=23,file=filename,status='replace') write(23,'(a)') 'expand_lattice' write(23,'(a)') '! Misalignments and errors: ' do i_param = 1, size(ma_params) if (ma_params(i_param)%amp == 0) cycle ma => ma_params(i_param) do i_ele = 1, ma_ring%n_ele_max ! must include lord elements above n_ele_track ele => ma_ring%ele(i_ele) if (ele%key /= ma%key) cycle if (.not. (ma%mask == "" .or. match_reg(ele%name, ma%mask))) cycle ! make sure element was free to vary: if (ele%slave_status .eq. super_slave$) cycle if (.not. attribute_free(ele, ma%attrib_name, err_print_flag = .false.)) cycle ! remove 'by_index' vs. 'by_name', in favor of using '##n' indexing: if (ele%ix_pointer .eq. 0) then write_name = ele%name else write(write_name,'(2a,i0)') trim(ele%name), '##', ele%ix_pointer endif write(23,'(4a,es18.9)') trim(write_name), '[', trim(ma%attrib_name), & '] := ', value_of_attribute(ele, ma%attrib_name, err) end do end do write(23,'(a)') '' write(23,'(a)') '' write(23,'(a)') '! Corrections:' write(23,'(a)') '' write(23,'(a)') '' close(23) end subroutine write_misalignments ! subroutine to write corrections to ma_lat.bmad subroutine write_corrections(ma_ring, p_ele, seed) implicit none type(lat_struct), intent(in), target :: ma_ring character(40) :: write_name = '', filename = '' type(parameter_ele_struct) :: p_ele(:) type(ele_struct), pointer :: ele integer i_param, i_ele, i_cor, ix, jx integer, intent(in) :: seed if (seed .gt. 0) then write(filename,'(a,i9.9,a)') 'ma_lat.', seed, '.bmad' else filename = 'ma_lat.bmad' endif open(unit=23,file=filename,access='append', status='old') do ix = 1, size(p_ele) if (len(trim(p_ele(ix)%name)) .eq. 0) cycle ele => ma_ring%ele(p_ele(ix)%loc) ! remove 'by_index' vs. 'by_name', in favor of using '##n' indexing: if (ele%ix_pointer .eq. 0) then write_name = ele%name else write(write_name,'(2a,i0)') trim(ele%name), '##', ele%ix_pointer endif write(23,'(4a,es18.9)') trim(write_name), '[', trim(p_ele(ix)%attrib_name), & '] := ', p_ele(ix)%attrib_ptr_meas%r end do close(23) end subroutine write_corrections !=========================================================================== !=========================================================================== !=========================================================================== subroutine init_corrector_calib(ring, correct) implicit none type(lat_struct), intent(inout), target :: ring type(correct_struct), target :: correct type (all_pointer_struct) a_ptr logical :: err = .false. type(ele_struct), pointer :: ele type(cor_grp), pointer :: cor_grp_ptr integer :: i_cor_grp, i_ele do i_cor_grp = 1, n_detcor_groups cor_grp_ptr => correct%cor(i_cor_grp) call upcase_string(cor_grp_ptr%attrib_name) call upcase_string(cor_grp_ptr%key_name) call upcase_string(cor_grp_ptr%mask) if (len(trim(cor_grp_ptr%attrib_name)) == 0) cycle ! no ele specified; cycle cor_grp_ptr%key = key_name_to_key_index(cor_grp_ptr%key_name, .true.) calib_ele_loop: do i_ele = 1, ring%n_ele_max ele => ring%ele(i_ele) ! check a few things first, to make sure it's a valid element to vary: ! regex match name: if (.not. (match_reg(trim(ele%name), cor_grp_ptr%mask))) cycle calib_ele_loop ! match element key to desired key: if (.not. ((ele%key .eq. cor_grp_ptr%key) .or. (cor_grp_ptr%key .eq. -1))) cycle calib_ele_loop ! is not a slave: if (ele%slave_status .eq. super_slave$) cycle calib_ele_loop ! is free to vary: if (.not. attribute_free(ele, cor_grp_ptr%attrib_name, err_print_flag = .false.)) cycle calib_ele_loop ! initialize all corrector calibrations to unity: call pointer_to_attribute(ele, 'CALIB', .true., a_ptr, err) a_ptr%r = 1.0 call set_flags_for_changed_attribute(ele, a_ptr%r) enddo calib_ele_loop enddo call lattice_bookkeeper(ring) end subroutine init_corrector_calib end module correct_ring_mod2