! 2011.07.06 - j. shanks ! program to simulate ac_eta, phase, and coupling data program sim_ac_and_cbar use bmad use bsim_interface use random_mod use dr_misalign_mod use correct_ring_mod2 use radiation_mod use cesr_basic_mod use sim_utils use nonlin_bpm_mod use sim_bpm_mod use mode3_mod use tune_tracker_mod use sim_ac_and_cbar_mod implicit none type (lat_struct) ring type (coord_struct), allocatable :: orb(:) logical err, everything_ok integer i, j, ix, jx, kx, dummy real(rp) harvest(7) integer :: seed = 0. real(rp) :: current = 1000000, skq02w_kick = 0. integer :: n_bpms = 0, nIter = 1 character(40) :: bpm_mask = "^DET\_[0-9]{2}[ewEW]$" ! default to using CESR naming convention ! For TBT tracking: real(rp) :: init_vec(6)=0., init_vec_orig(6) = 0. ! for starting oscillations for phase 'measurement' logical :: track6x6 = .false., damping = .false. integer :: n_turns = 1024, n_damping = 0 !if 6x6 tracking, n_damping = # of turns to wait before saving data ! Related to phase/coupling simulation real(rp), allocatable :: phi(:,:,:), delta_phi_ax(:), delta_phi_by(:) real(rp) :: dphi_ax_avg = 0, dphi_by_avg = 0 real(rp), allocatable :: cbar_meas(:,:,:), cbar_bmad(:,:,:) ! (:,2,2) = across all BPMs real(rp), allocatable :: A(:,:,:) ! = (n_bpms,2,2) = (i_bpm, (Aa,Ab),(Ba,Bb)) real(rp) :: rms_cbar(2,2)=0, rms_phia = 0, rms_phib = 0, rms_ac_x = 0, rms_ac_y = 0 ! Related to AC dispersion simulation: type(cesr_all_data_struct) ac_eta logical :: use_new_method = .false. ! parameters for misalignment; to be read in and loaded into bpm_error_sigmas and ma structures type(det_error_struct) :: bpm_error_sigmas type(det_struct), allocatable :: bpm(:) type(ma_struct) ma(20) logical :: misalign_magnets = .false. real(rp) :: shear_x = 0., gain_sigma = 0., timing_sigma = 0. real(rp) :: bpm_offset = 0., bpm_rotation = 0., bpm_noise = 0. ! parameters for tune tracker TYPE(tt_param_struct) :: tt_params(max_tt) !max_tt=3 set on tune_tracker module INTEGER :: n_TTs = 0 ! number of tune trackers INTEGER id(max_tt) real(rp) :: target_tunes(3) = 0. ! target tunes for q_tuning real(rp), allocatable :: dk1(:) ! for q_tuning ! File handling: integer ios, version, arg_num, iargc, readstatus integer :: lun_list(20) ! for recording available luns from lunget character(20) :: output='cesrv' character(100) lat_file, file_name, lat_file_name, path, phase_file_name, ac_file_name, rms_file_name ! intelligent bookkeeping (default to old): logical :: auto_bookkeeper = .true. namelist /parameters/ lat_file, seed, bpm_offset, bpm_rotation, bpm_noise, track6x6, n_damping, & shear_x, current, gain_sigma, timing_sigma, nIter, init_vec, n_turns, output, skq02w_kick,& n_TTs, tt_params, misalign_magnets, ma, use_new_method, target_tunes, bpm_mask, auto_bookkeeper arg_num=iargc() if(arg_num==0) then file_name='sim_ac_and_cbar.in' else call getarg(1, file_name) end if call string_trim (file_name, file_name, ix) open (unit= 1, file = file_name, status = 'old', iostat = ios) if(ios.ne.0)then write(*,*) write(*,*) 'ERROR: CANNOT OPEN FILE: ', trim(file_name) endif ! read in the parameters rewind(1) read(1, nml = parameters,iostat=readstatus) if(readstatus > 0) then print *,"CAN NOT READ FROM ",file_name, " IOSTAT = ", readstatus stop end if bmad_com%auto_bookkeeper = auto_bookkeeper !initialize randomizer's seed (harvest) if (seed < 1) then seed = 0 ! seed randomizer with CPU time end if call ran_seed_put(seed) output = trim(output) ! load parameters into bpm_error_sigmas, tt_params structure: bpm_error_sigmas%bpm_offset = bpm_offset bpm_error_sigmas%bpm_rotation = bpm_rotation bpm_error_sigmas%shear_x = shear_x bpm_error_sigmas%gain_sigma = gain_sigma bpm_error_sigmas%timing_sigma = timing_sigma init_vec_orig(:) = init_vec(:) write(*,*) "n_TTs = ", n_TTs ! init lattice if(lat_file(1:8) == 'digested')then call read_digested_bmad_file (lat_file, ring, version) else call fullfilename(lat_file,lat_file) call bmad_parser (lat_file, ring) endif dummy = splitfilename(lat_file,path,lat_file_name) lat_file_name = lat_file_name(1:len_trim(lat_file_name)-4) ! turn off all kicks / separators ! ring%ele(:)%value(hkick$) = 0 ! ring%ele(:)%value(vkick$) = 0 ! Turn on sk_q02w @ user-specified kick: !ring%ele(38)%value(k1$) = skq02w_kick !call make_mat6(ring%ele(38),ring%param) bmad_com%radiation_damping_on = .true. bmad_com%radiation_fluctuations_on = .false. ! tune trackers will dominate well over these effects call set_on_off(rfcavity$, ring, on$) call twiss3_at_start(ring,err) call twiss3_propagate_all(ring) call closed_orbit_calc(ring, orb, 6) call twiss_at_start(ring) call twiss_propagate_all(ring) ! Locate BPMs we wish to study: call find_bpms(ring, bpm_mask, bpm) n_bpms = ubound(bpm,1) write(*,*) "n_bpms = ", n_bpms ! Find resolution in terms of button signal error: call resolution_to_button_error(bpm(1), current, bpm_noise) bpm(:)%butn_res_sigma = bpm(1)%butn_res_sigma if ((n_TTs .gt. 0) .and. (track6x6 .eqv. .true.)) then ! Check if the element for each tune tracker's kicker has kick attributes. call check_tt_kickers(ring,tt_params,n_TTs) call tt_init(ring, orb, tt_params, n_TTs, id) endif call string_trim (file_name, file_name, ix) file_name = file_name(:ix-3) write(rms_file_name,'(a4,a,a4)') 'rms_', trim(file_name), '.out' lun_list(1) = lunget() open(unit=lun_list(1),file=trim(rms_file_name)) write(lun_list(1),'(a)') '!kx cbar22 cbar12 cbar11 phi_a phi_b' !===================================================================== !===== Begin averaging loop! ========================================= !===================================================================== allocate(cbar_bmad(n_bpms,2,2), delta_phi_ax(n_bpms), delta_phi_by(n_bpms)) do kx = 1,nIter ! number of iterations to average over !=== Generate BPM errors ============================================= ! Note that putting values here have no effect on tracking. They're ! only for storage lun_list(2) = lunget() open(unit=lun_list(2), file='bpms.ma', status='replace') write(lun_list(2), '(a7,7a14)') "!index", "x-offset", "y-offset", "tilt", "g1", "g2", "g3", "g4" write(lun_list(2), '(a7, 7a14)') "" ! call dr_misalign_...(blah) ! directly from ring_ma: dr_misalign_params%sigma_cutoff = 3. dr_misalign_params%alignment_multiplier = 1 if (misalign_magnets .eqv. .true.) then call dr_misalign(ring, ma) endif do i=1,n_bpms call bpm_errors(bpm(i), bpm_error_sigmas, ring%ele(bpm(i)%ix_quad)) write(lun_list(2), '(i7,7e14.4)') bpm(i)%ix_db, bpm(i)%x_offset, & bpm(i)%y_offset, bpm(i)%tilt, bpm(i)%gain(:) enddo close(unit=lun_list(2)) call closed_orbit_calc(ring, orb, 6) call lat_make_mat6(ring, -1, orb) call twiss_at_start(ring) call twiss_propagate_all(ring) call twiss3_at_start(ring, err) call twiss3_propagate_all(ring) ! q_tune to user-specified tunes: global_com%exit_on_error = .false. if (any(target_tunes .gt. 1.e-12)) everything_ok = set_tune3(ring, target_tunes) ! write misalignments in lattice format: call write_misalignments(ring, ma, -1) ! end block from ring_ma if (.not. any(init_vec_orig(:) > 1.e-12)) then ! floating point comparison; assume < 1.e-12 is negligible init_vec(:) = orb(0)%vec(:) else init_vec(:) = init_vec_orig(:) endif ! initialize each tune tracker do ix=1, n_TTs call reset_dTT(id(ix), tt_params(ix)) enddo call sim_ac_and_cbar_data(ring, bpm, n_turns, n_TTs, tt_params, id, n_damping, init_vec, use_new_method, ac_eta, phi, cbar_meas, A) ! populate cbar_bmad structure: do ix = 1, n_bpms call c_to_cbar(ring%ele(bpm(ix)%ix_lat),cbar_bmad(ix,:,:)) enddo if (n_TTs .gt. 0) then do ix = 1, n_TTs call dest_dTT(id(ix)) enddo endif lun_list(6) = lunget() open(unit=lun_list(6), file='amplitude.out', status='replace') write(lun_list(6),'(a7,8a18)') 'bpm_ix', 'Aa', 'Ab', 'Ba', 'Bb', 'Aa/sqrt(beta.a)', 'Ab/sqrt(beta.b)', & 'Ba/sqrt(beta.a)', 'Bb/sqrt(beta.b)' do ix=1,n_bpms write(lun_list(6),'(i7,8es18.4)') bpm(ix)%ix_db, A(ix,1,1), A(ix,1,2), A(ix,2,1), A(ix,2,2), & A(ix,1,1) / sqrt(ring%ele(bpm(ix)%ix_lat)%a%beta), & A(ix,1,2) / sqrt(ring%ele(bpm(ix)%ix_lat)%b%beta), & A(ix,2,1) / sqrt(ring%ele(bpm(ix)%ix_lat)%a%beta), & A(ix,2,2) / sqrt(ring%ele(bpm(ix)%ix_lat)%b%beta) enddo close(lun_list(6)) enddo ! matches [do kx=1, nIter] close(lun_list(1)) !=============================================================== lun_list(7) = lunget() open(unit=lun_list(7),file='cbar.out',status='replace') write(lun_list(7),'(a7,6a18)') '!bpm_ix', 'cbar12 (meas)', 'cbar12 (bmad)', 'cbar22 (meas)', & 'cbar22 (bmad)', 'cbar11 (meas)', 'cbar11 (bmad)' do ix=1,n_bpms write(lun_list(7),'(i7,6es18.4)') bpm(ix)%ix_db, cbar_meas(ix,1,2), cbar_bmad(ix,1,2), & cbar_meas(ix,2,2), cbar_bmad(ix,2,2), cbar_meas(ix,1,1), cbar_bmad(ix,1,1) enddo close(lun_list(7)) lun_list(8) = lunget() open(unit=lun_list(8), file='ac_eta.out',status='replace') lun_list(9) = lunget() open(unit=lun_list(9), file='phase.out', status='replace') write(lun_list(9),'(a7,6a18)')"bpm_ix","phi_ax","phi_by","phi_ax (des)","phi_by (des)", "phi_ax (meas-des)", "phi_by (meas-des)" write(lun_list(8),'(a7,8a18)')"bpm_ix","eta_x_meas","etap_x_meas","eta_y_meas","etap_y_meas",& "eta_x_calc","etap_x_calc","eta_y_calc","etap_y_calc" do i=1,n_bpms write(lun_list(9),'(i7,6es18.4)') bpm(i)%ix_db, phi(i,1,1)+dphi_ax_avg, phi(i,2,2)+dphi_by_avg, & ring%ele(bpm(i)%ix_lat)%a%phi, ring%ele(bpm(i)%ix_lat)%b%phi, & phi(i,1,1)+dphi_ax_avg-ring%ele(bpm(i)%ix_lat)%a%phi, & phi(i,2,2)+dphi_by_avg-ring%ele(bpm(i)%ix_lat)%b%phi write(lun_list(8),'(i7,8es18.4)') bpm(i)%ix_db, ac_eta%ac_eta_x(i)%value, ac_eta%ac_etap_x(i)%value, & ac_eta%ac_eta_y(i)%value, ac_eta%ac_etap_y(i)%value, & ring%ele(bpm(i)%ix_lat)%mode3%x%eta, ring%ele(bpm(i)%ix_lat)%mode3%x%etap,& ring%ele(bpm(i)%ix_lat)%mode3%y%eta, ring%ele(bpm(i)%ix_lat)%mode3%y%etap enddo close(lun_list(8)) close(lun_list(9)) write(phase_file_name,'(a,a13)') trim(file_name),'.all_data.dat' lun_list(10) = lunget() open(unit=lun_list(10),file=phase_file_name, status='replace') if (output == 'cesrv') then write(lun_list(10),'(a)') "&DATA_PARAMETERS ! Fortran90 namelist." write(lun_list(10),'(a)') " file_type = 'ALL DATA' ! Identifies file as an phase/cbar file." write(lun_list(10),'(a)') " comment = 'Simulated data'" write(lun_list(10),'(3a)') " lattice = '", trim(upcase(lat_file_name)), "'" write(lun_list(10),'(a)') "/" write(lun_list(10),'(a)') '' write(lun_list(10),'(a)') "&PHASE_PARAMETERS " write(lun_list(10),'(a)') " species = 1" write(lun_list(10),'(a,es18.4)') " horiz_freq = ", ring%a%tune/(2*3.1415926) write(lun_list(10),'(a,es18.4)') " vert_freq = ", ring%b%tune/(2*3.1415926) ! write(41,'(a)') " data_date = '1998-JAN-06 18:05:36' ! This is optional." ! write(41,'(a)') " save_set = 51220 ! This is optional." write(lun_list(10),'(a)') "/" write(lun_list(10),'(a)') "" write(lun_list(10),'(a)') " &ALL_DATA" write(lun_list(10),'(a10,3a18)') "!", "x", 'y', 'Valid?' do i=1,n_bpms write(lun_list(10),'(a,i0,a,2es18.4,a16)') 'phase(',bpm(i)%ix_db, ') = ', phi(i,1,1)+dphi_ax_avg, phi(i,2,2)+dphi_by_avg, 'T' write(lun_list(10),'(a,i0,a,es18.4,a34)') 'cbar11(',bpm(i)%ix_db, ') = ', cbar_meas(i,1,1), 'T' write(lun_list(10),'(a,i0,a,es18.4,a34)') 'cbar12(',bpm(i)%ix_db, ') = ', cbar_meas(i,1,2), 'T' write(lun_list(10),'(a,i0,a,es18.4,a34)') 'cbar22(',bpm(i)%ix_db, ') = ', cbar_meas(i,2,2), 'T' write(lun_list(10),'(a,i0,a,2es18.4,a16)') 'ac_eta(',bpm(i)%ix_db, ') = ', ac_eta%ac_eta_x(i)%value, ac_eta%ac_eta_y(i)%value, 'T' write(lun_list(10),'(a)') '' enddo write(lun_list(10),'(a)') '/' elseif (output == 'tao') then write(lun_list(10),'(a)') "set global lattice_calc_on = F" write(lun_list(10),'(a)') "set global plot_on = F" do i=1, n_bpms write(lun_list(10), '(a22,i2,a9,es14.4)') "set data bpm_phase.a[",bpm(i)%ix_db,"]|meas = ", phi(i,1,1) +dphi_ax_avg write(lun_list(10), '(a22,i2,a9,es14.4)') "set data bpm_phase.b[",bpm(i)%ix_db,"]|meas = ", phi(i,2,2) +dphi_by_avg write(lun_list(10), '(a22,i2,a9,es14.4)') "set data bpm_cbar.11b[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,1,1) write(lun_list(10), '(a22,i2,a9,es14.4)') "set data bpm_cbar.12a[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,1,2) write(lun_list(10), '(a22,i2,a9,es14.4)') "set data bpm_cbar.22a[",bpm(i)%ix_db,"]|meas = ", cbar_meas(i,2,2) enddo write(lun_list(10),'(a)') "set global lattice_calc_on = T" ! write(41,'(a)') "set global plot_on = T" ! write(41,'(a)') "place top phase" ! write(41,'(a)') "place bottom cbar" ! write(41,'(a)') "plot top meas-design" ! write(41,'(a)') "plot bottom meas-design" else write(*,*) "Not accepted output format: ", output write(*,*) "No program-specific analysis output generated." endif close(lun_list(10)) end program sim_ac_and_cbar