!+ ! In the calling program define the structure by: ! RECORD / NONLIN_ELE_STRUCT / NONLIN ! ! Example values: ! ! NONLIN.CHROM_X ! Delta tune /delta E ! NONLIN.NON_ELE(i).X.DBETA ! (Delta beta_x / delta E) ! NONLIN.NON_ELE(i).X.DPHI ! (Delta phi_x / delta E) where phi is ! phase from start to element i ! NONLIN.NON_ELE(i).DET.FOUR_BY4(j) ! determinant of 4 X 4 matrix from ! ! start to element for jth amplitude ! NONLIN.NON_ELE(i).DET.TWO_BY2_ul(j) ! determinant of 2 X 2 upper left ! ! block of jacobian matrix from ! ! start to element i ! NONLIN.NON_ELE(i).DET.TWO_BY2_lr(j) ! determinant of 2 X 2 lower right ! ! block of jacobian matrix from start ! ! to element i ! NONLIN.NON_ELE(i).X.DBETA_DCOS ! fractional change in beta for ! orbit that propogates as cos(phi) ! from ip and closed orbit ! module nonlin_mod use bmad use bmadz_mod use zquad_lens_mod use constraints_mod type chrom_struct real(rp) dbeta, dalpha, dphi, eta, etap, dbeta_dcos, dbeta_dsin real(rp) beta, alpha, phi end type type determinant_struct real(rp) four_by4(10), two_by2_ul(10), two_by2_lr(10) end type type nonele_struct type (chrom_struct) a, b ! dTwiss/d delta at end of element type (twiss_struct) tx,ty,tz type (determinant_struct) det type (coord_struct) co_off(3) end type type nonlin_ele_struct type (coord_struct) delta_orb(10) type (nonele_struct), pointer :: non_ele(:), nonele_pos(:) type (nonele_struct), pointer :: nonele_ele(:) real(rp) amplitude_x(10), amplitude_y(10) real(rp) energy_offset(3) !energy_offset(1 to 3)= -dE/E, 0, +dE/E real(rp) chrom_x, chrom_y real(rp) dchrom_x, dchrom_y real(rp) tonality_x, tonality_y !kHz real(rp) growth_rate(3) real(rp) tune_x, tune_y, trace(10) integer n_det_calc, energies end type type dvar_struct real(rp) dvar0 real(rp) dvar(100) end type type dmat_dvar_struct type (dvar_struct), pointer :: ele(:) end type type moment_struct type (dmat_dvar_struct) dc_mat(2,2) real(rp) k2l(275) ! sextupole strength * sextupole length real(rp) m1(275), m2(275), m3(275), m4(275), m5(275) real(rp) m6(275), m7(275), m8(275), m9(275), m10(275) real(rp), pointer :: dbeta_dpretz_x(:, :) !, dbeta_dpretz_x_0(:) !delta_beta / pretz real(rp), pointer :: dbeta_dpretz_y(:, :) !, dbeta_dpretz_y_0(:) real(rp), pointer :: dbeta_dcos_x(:, :) !delta_beta / pretz + cos real(rp), pointer :: dbeta_dcos_y(:, :) real(rp), pointer :: dbeta_dsin_x(:, :) !delta_beta / pretz + sin real(rp), pointer :: dbeta_dsin_y(:, :) real(rp), pointer :: dbeta_x_0_norm(:) !delta_beta /delta_energy real(rp), pointer :: dbeta_y_0_norm(:) real(rp), pointer :: dbeta_x(:, :), dbeta_x_0(:) !delta_beta /delta_energy real(rp), pointer :: dbeta_y(:, :), dbeta_y_0(:) real(rp), pointer :: a_real(:, :), a_image(:, :) !coupling coefficient real(rp), pointer :: b_real(:, :), b_image(:, :) real(rp) chrom_x(275), chrom_y(275), chrom_x_0, chrom_y_0, chrom_x_ir, chrom_y_ir real(rp) tonality_x(275), tonality_y(275) real(rp), pointer :: mdbeta_x(:), mdbeta_y(:) real(rp), pointer :: mdbeta_dpretz_x(:) real(rp), pointer :: mdbeta_dpretz_y(:) real(rp), pointer :: mdbeta_dcos_x(:) real(rp), pointer :: mdbeta_dcos_y(:) real(rp), pointer :: mdbeta_dsin_x(:) real(rp), pointer :: mdbeta_dsin_y(:) real(rp), pointer :: ma_real(:), ma_image(:) real(rp), pointer :: mb_real(:), mb_image(:) real(rp) mchrom_x, mchrom_y, mtonality_x, mtonality_y, mtonality_x_west real(rp) mtonality_y_west real(rp) s_h, s_v real(rp) chi2 !chi2 of moment minimization integer n_sext !number of sextupoles (length of vector) integer ix(275) !lat element number of sextupole character*40 name(275) !element name of sextupole logical use(275) end type type det_ele_struct real(rp) four_by4(10), two_by2_ul(10), two_by2_lr(10), trace(10) end type type select_mat real(rp) select(6,6) end type type bmadz_track_struct type (coord_struct), allocatable :: orb(:) end type !------------------------------------------------------------------------ ! To calculate the change in action due to a single sextupole: ! Given the action-angle coords of a particle ! ! sqrt_j(3) = (/ sqrt(j_a), sqrt(j_b), sqrt(j_z) /) ! phi(3) = (/ phi_a, phi_b, phi_z /) ! ! Then the total change in the a mode action DJA is ! ! type (sextupole_res_struct) res ! call sextupole_resonance_calc (sextupole_ele, res) ! ! dja = 0 ! sqrt_j_vec = (/ sqrt(j_a), sqrt(j_b), sqrt(j_z) /) ! action for a particle ! phi_vec = (/ phi_a, phi_b, phi_z /) ! angle for a particle ! do i = 1, res%n_res ! theta_a = sum(res%a_pole(i)%iq * phi_vec) ! factor = sextupole%value(k2$) * sextupole%value(l$) * & ! product(sqrt_j_vec**res%a_pole(i)%ij) ! dja = dja + factor * (res%a_pole(i)%dj_cos * cos(theta_a) + & ! res%a_pole(i)%dj_sin * sin(theta_a)) ! enddo ! ! with a similar equation for the b mode. type single_res real(rp) dj_cos, dj_sin ! change in action components. real(rp) amp ! amplitude of resonance real(rp) cos, sin ! sin and cos resonance components integer iq(3) ! p,q,r array as in p*Q_a + q*Q_b + r*Q_s integer ij(3) ! power of resonance: J_a^ij(1) + J_b^ij(2) + J_z^ij(3) end type type sextupole_res_struct type (single_res) a(18), b(18) ! individual resonanes type (single_res) x(18), y(18) ! individual resonanes real(rp) x_weight(18), y_weight(18) ! weight in merit function integer n_res ! number of resonances integer a_map(-2:3, -2:3, -2:2) integer b_map(-2:3, -2:3, -2:2) integer x_map(-2:3, -2:3, -2:2) integer y_map(-2:3, -2:3, -2:2) end type !------------------------------------------------------------------------ contains subroutine allocate_nonlin_ele( nonlin, n_ele_maxx ) implicit none type (nonlin_ele_struct) :: nonlin integer n_ele_maxx type (determinant_struct) :: zero_det type (coord_struct):: zero_coord(3) logical alloc_it ! allocate nonlin pointers alloc_it = .false. if (associated (nonlin%non_ele)) then if (ubound(nonlin%non_ele, 1) /= n_ele_maxx) then call deallocate_nonlin_ele (nonlin) alloc_it = .true. endif else alloc_it = .true. endif if (alloc_it) then allocate( nonlin%non_ele(0:n_ele_maxx) ) allocate( nonlin%nonele_pos(0:n_ele_maxx) ) allocate( nonlin%nonele_ele(0:n_ele_maxx) ) endif ! Initialize all to zero zero_det%four_by4 = 0 zero_det%two_by2_ul = 0 zero_det%two_by2_lr = 0 zero_coord(1)%vec(:) = 0 zero_coord(2)%vec(:) = 0 zero_coord(3)%vec(:) = 0 nonlin%non_ele(:) = nonele_struct(chrom_struct(0,0,0,0,0,0,0,0,0,0), & chrom_struct(0,0,0,0,0,0,0,0,0,0),twiss_struct(0,0,0,0,0,0,0,0,0,0), & twiss_struct(0,0,0,0,0,0,0,0,0,0), twiss_struct(0,0,0,0,0,0,0,0,0,0), & zero_det, zero_coord) nonlin%nonele_pos(:) = nonlin%non_ele(1) nonlin%nonele_ele(:) = nonlin%non_ele(1) end subroutine allocate_nonlin_ele !------------------------------------------------------------------------ subroutine deallocate_nonlin_ele( nonlin ) implicit none type (nonlin_ele_struct) :: nonlin integer error_stat deallocate( nonlin%non_ele, stat=error_stat ) deallocate( nonlin%nonele_pos, stat=error_stat ) deallocate( nonlin%nonele_ele, stat=error_stat ) end subroutine deallocate_nonlin_ele !------------------------------------------------------------------------ subroutine allocate_moment( moment, n_ele_maxx ) implicit none type (moment_struct) :: moment integer n_ele_maxx ! if (associated (moment%a_real)) then if (ubound(moment%a_real, 1) == n_ele_maxx) return call deallocate_moment (moment) endif allocate( moment%dbeta_dpretz_x(0:n_ele_maxx, 275) ) allocate( moment%dbeta_dpretz_y(0:n_ele_maxx, 275) ) ! allocate( moment%dbeta_dpretz_x_0(0:n_ele_maxx) ) ! allocate( moment%dbeta_dpretz_y_0(0:n_ele_maxx) ) allocate( moment%dbeta_dcos_x(0:n_ele_maxx, 275) ) allocate( moment%dbeta_dcos_y(0:n_ele_maxx, 275) ) allocate( moment%dbeta_dsin_x(0:n_ele_maxx, 275) ) allocate( moment%dbeta_dsin_y(0:n_ele_maxx, 275) ) allocate( moment%dbeta_x_0_norm(0:n_ele_maxx) ) allocate( moment%dbeta_y_0_norm(0:n_ele_maxx) ) allocate( moment%dbeta_x(0:n_ele_maxx, 275) ) allocate( moment%dbeta_x_0(0:n_ele_maxx) ) allocate( moment%dbeta_y(0:n_ele_maxx, 275) ) allocate( moment%dbeta_y_0(0:n_ele_maxx) ) allocate( moment%a_real(0:n_ele_maxx, 275) ) allocate( moment%a_image(0:n_ele_maxx, 275) ) allocate( moment%b_real(0:n_ele_maxx, 275) ) allocate( moment%b_image(0:n_ele_maxx, 275) ) allocate( moment%mdbeta_x(0:n_ele_maxx) ) allocate( moment%mdbeta_y(0:n_ele_maxx) ) allocate( moment%mdbeta_dpretz_x(0:n_ele_maxx) ) allocate( moment%mdbeta_dpretz_y(0:n_ele_maxx) ) allocate( moment%mdbeta_dcos_x(0:n_ele_maxx) ) allocate( moment%mdbeta_dcos_y(0:n_ele_maxx) ) allocate( moment%mdbeta_dsin_x(0:n_ele_maxx) ) allocate( moment%mdbeta_dsin_y(0:n_ele_maxx) ) allocate( moment%ma_real(0:n_ele_maxx) ) allocate( moment%ma_image(0:n_ele_maxx) ) allocate( moment%mb_real(0:n_ele_maxx) ) allocate( moment%mb_image(0:n_ele_maxx) ) end subroutine allocate_moment !------------------------------------------------------------------------ subroutine deallocate_moment( moment ) implicit none type (moment_struct) :: moment integer error_stat deallocate( moment%dbeta_dpretz_x, stat=error_stat ) deallocate( moment%dbeta_dpretz_y, stat=error_stat ) ! deallocate( moment%dbeta_dpretz_x_0) ! deallocate( moment%dbeta_dpretz_y_0) deallocate( moment%dbeta_dcos_x, stat=error_stat ) deallocate( moment%dbeta_dcos_y, stat=error_stat ) deallocate( moment%dbeta_dsin_x, stat=error_stat ) deallocate( moment%dbeta_dsin_y, stat=error_stat ) deallocate( moment%dbeta_x_0_norm, stat=error_stat ) deallocate( moment%dbeta_y_0_norm, stat=error_stat ) deallocate( moment%dbeta_x, stat=error_stat ) deallocate( moment%dbeta_x_0, stat=error_stat ) deallocate( moment%dbeta_y, stat=error_stat ) deallocate( moment%dbeta_y_0, stat=error_stat ) deallocate( moment%a_real, stat=error_stat ) deallocate( moment%a_image, stat=error_stat ) deallocate( moment%b_real, stat=error_stat ) deallocate( moment%b_image, stat=error_stat ) deallocate( moment%mdbeta_x, stat=error_stat ) deallocate( moment%mdbeta_y, stat=error_stat ) deallocate( moment%mdbeta_dpretz_x, stat=error_stat ) deallocate( moment%mdbeta_dpretz_y, stat=error_stat ) deallocate( moment%mdbeta_dcos_x, stat=error_stat ) deallocate( moment%mdbeta_dcos_y, stat=error_stat ) deallocate( moment%mdbeta_dsin_x, stat=error_stat ) deallocate( moment%mdbeta_dsin_y, stat=error_stat ) deallocate( moment%ma_real, stat=error_stat ) deallocate( moment%ma_image, stat=error_stat ) deallocate( moment%mb_real, stat=error_stat ) deallocate( moment%mb_image, stat=error_stat ) end subroutine deallocate_moment !........................................................................ ! ! Subroutine : CHROMATICITY_CALC (LAT, CO, GLOBAL, err_flag) ! ! Description: Subroutine to calculate the chromaticities by computing the ! tune change when then energy is changed. ! NOTE: CURRENTLY LINEARITY IS ASSUMED. ! ! Arguments : ! Input: ! LAT -- lat_struct: Lat ! CO(:) -- Coord_struct: THIS IS NOT USED SINCE LINEARITY IS ASSUMED!!! ! GLOBAL -- Global_struct: Structure ! .DELTA_E_CHROM -- Real(rp): Delta energy used for the calculation. ! If 0 then default of 1.0e-4 is used. ! STATUS -- Common block status structure ! .TYPE_OUT -- Logical: If .true. then will type a message for ! non ok$ STATUS ! ! Output: ! GLOBAL.CHROMX -- Real(rp): Horizontal chromaticity. ! GLOBAL.CHROMY -- Real(rp): Vertical chromaticity. ! STATUS -- Common block status structure ! .STATUS -- Calculation status. See MAT_SYMP_DECOUPLE for info. !........................................................................ subroutine chromaticity_calc(lat, co, global, err_flag) implicit none type (lat_struct) lat type (coord_struct), allocatable :: co(:) type (global_struct) global logical err_flag ! co(0)%vec(5) = co(0)%vec(5) ! so no compiler complaints global%delta_e_chrom = 1.e-6 call chrom_calc_z (lat, global%delta_e_chrom, global%chromx, global%chromy, & global%c_mat_chrom, err_flag) return end subroutine chromaticity_calc !........................................................................ ! ! Subroutine : EMITTANCE_CALC (LAT, CO, CONS, GLOBAL) ! ! Description: Subroutine to calcuate the horizontal emittance, energy ! spread, and other global parameters. ! Note: ! Momentum compaction = .SYNCH_INT(1) / lat_length ! ! Arguments : ! Input: ! LAT -- lat_struct: Lat ! CO(:) -- Coord_struct: Closed orbit array ! CON -- Constraint_struct: ! ! Output: ! GLOBAL.NOWIG.X_EMIT -- Real(rp): "Horizontal" emittance. ! GLOBAL.NOWIG.SIGE_E -- Real(rp): Sigma_E/E energy spread. ! GLOBAL.WIG.X_EMIT -- Real(rp): "Horizontal" emittance ! including wigglers. ! GLOBAL.WIG.SIGE_E -- Real(rp): Sigma_E/E energy spread ! including wigglers. ! GLOBAL.NOWIG.E_LOSS -- Real(rp): Energy loss/turn (Mev) ! GLOBAL.WIG.E_LOSS -- Real(rp): Energy loss/turn including ! wigglers ! GLOBAL.NOWIG.SYNCH_INT(5) -- Real(rp): array of synchrotron ! integrals ! GLOBAL.WIG.SYNCH_INT(5) -- Real(rp): array of synchrotron ! integrals w/ wigs ! ! ! Mod/Commons: ! ! Calls : ! ! Author : ! ! Modified : ! !........................................................................ subroutine emittance_calc (lat, co, con, global) implicit none type (lat_struct) lat type (coord_struct), allocatable :: co(:) type (constraint_struct) con type (global_struct) global type (normal_modes_struct), save :: m_bend, m_wig, mode real(rp) i1b, i2b, i3b, i4b, i5b, i1w, i2w, i3w, i4w, i5w real(rp) emit_factor, sigE_E_factor real(rp) arg integer n integer, save :: ix_positrons = 0 logical :: wig_as_bends = .false., wig_as_map = .true. ! if(wig_as_bends)then call emit_calc (lat, bends$, m_bend) call emit_calc (lat, wigglers$, m_wig) endif if(wig_as_map)then ! if(lat%param%symmetry == ew_antisymmetry$) then ! print *,' ERROR: subroutine radiation_integrals cannot be used with ' ! print *,' east - west symmetry' ! stop ! endif call radiation_integrals (lat, co, mode, ix_positrons) global%wig%synch_int(1:3) = mode%synch_int(1:3) global%wig%synch_int(4:5) = mode%a%synch_int(4:5) global%wig%e_loss = mode%e_loss global%wig%x_emit = mode%a%emittance global%wig%sige_e = mode%sigE_E global%nowig%synch_int(1:3) = mode%synch_int(1:3) global%nowig%synch_int(4:5) = mode%a%synch_int(4:5) global%nowig%e_loss = mode%e_loss global%nowig%x_emit = mode%a%emittance global%nowig%sige_e = mode%sigE_E return endif ! if (lat%param%symmetry == ew_antisymmetry$) then ! n = 2 ! else n = 1 ! endif if(wig_as_bends)then i1b = m_bend%synch_int(1) i2b = m_bend%synch_int(2) i3b = m_bend%synch_int(3) i4b = m_bend%a%synch_int(4) i5b = m_bend%a%synch_int(5) i1w = m_wig%synch_int(1) i2w = m_wig%synch_int(2) i3w = m_wig%synch_int(3) i4w = m_wig%a%synch_int(4) i5w = m_wig%a%synch_int(5) emit_factor = m_bend%a%emittance / (i5b / (i2b - i4b)) sigE_E_factor = 0. arg = (i3b / (2*i2b + i4b)) if(arg > 0.)sigE_E_factor = m_bend%sigE_E / sqrt(i3b / (2*i2b + i4b)) global%nowig%synch_int(1) = n * i1b global%nowig%synch_int(2) = n * i2b global%nowig%synch_int(3) = n * i3b global%nowig%synch_int(4) = n * i4b global%nowig%synch_int(5) = n * i5b global%nowig%e_loss = n * m_bend%e_loss !ev global%nowig%x_emit = m_bend%a%emittance global%nowig%sige_e = m_bend%sigE_E global%wig%synch_int(1) = n * (i1b + i1w) global%wig%synch_int(2) = n * (i2b + i2w) global%wig%synch_int(3) = n * (i3b + i3w) global%wig%synch_int(4) = n * (i4b + i4w) global%wig%synch_int(5) = n * (i5b + i5w) global%wig%e_loss = global%nowig%e_loss + n * m_wig%e_loss global%wig%x_emit = emit_factor * (i5b + i5w) / (i2b + i2w - i4b - i4w) arg = ((i3b + i3w) / (2*i2b + 2*i2w + i4b + i4w)) if(arg > 0.)global%wig%sige_e = sigE_E_factor * & sqrt((i3b + i3w) / (2*i2b + 2*i2w + i4b + i4w)) global%ewig%synch_int(1) = & (global%nowig%synch_int(1) + global%wig%synch_int(1)) / 2 global%ewig%synch_int(2) = & (global%nowig%synch_int(2) + global%wig%synch_int(2)) / 2 global%ewig%synch_int(3) = & (global%nowig%synch_int(3) + global%wig%synch_int(3)) / 2 global%ewig%synch_int(4) = & (global%nowig%synch_int(4) + global%wig%synch_int(4)) / 2 global%ewig%synch_int(5) = & (global%nowig%synch_int(5) + global%wig%synch_int(5)) / 2 global%ewig%e_loss = (global%nowig%e_loss + global%wig%e_loss) / 2 global%ewig%x_emit = emit_factor * & (2*i5b + i5w) / (2*i2b + i2w - 2*i4b - i4w) arg = ((2*i3b + i3w) / (4*i2b + 2*i2w + 2*i4b + i4w)) if(arg > 0.)global%ewig%sige_e = sigE_E_factor * & sqrt((2*i3b + i3w) / (4*i2b + 2*i2w + 2*i4b + i4w)) endif return end subroutine emittance_calc !........................................................................ ! ! Subroutine : NONLINEAR_CLOSED_ORBIT (LAT, CO, LAT_OUT, CO_OUT, NONLIN, ! DET_ELE) ! ! Description: Subroutine to calcuate the closed orbit and jacobian ! (6X6 matrix) at every element including nonlinearities. ! ! Arguments : ! Input: ! LAT -- lat_struct: Lat ! CO -- Coord_struct: Starting point of closed orbit neglecting ! nonlinearities or zero if there are no synchrotron ! oscillations, the energy offset is co%vec(6) ! NONLIN -- Nonlin_struct: ! %N_DET_CALC -- The number of n_det_calc (default=3) ! %DELTA_ORB -- OFFSETS used for start of tracking to find Jacobian. ! ! Output: ! CO_OUT -- Coord_struct: Closed orbit at each element ! LAT_OUT -- lat_struct: Lat with jacobian ! DET_ELE -- Det_ele_struct: Det_ele determinants of transfer matrices ! from start to element for each iteration !........................................................................ subroutine nonlinear_closed_orbit (lat, co, lat_out, co_out, nonlin, det_ele, track_state) implicit none type (lat_struct) lat, lat_out type (coord_struct) co type (coord_struct), allocatable :: co_out(:) type (nonlin_ele_struct) nonlin type (det_ele_struct), allocatable :: det_ele(:) type (bmadz_track_struct) t(7) integer i, j, k, n_det_calc, ntrk, n, nk1, nt, track_state real(rp) vec_init(7,6), unit(6,6), z(6), rco(6), y(7,7), R(6,6), NegR(6,6) real(rp) R0(6,6), R_temp(6,6), R_ft(6,6), R_inverse(6,6) real(rp) diff_mat(6,6), inverse_diff_mat(6,6) real(rp) f,det real(rp) xemit logical debug real(rp) xemit_0 ! do i=1,7 allocate( t(i)%orb(0:lat%n_ele_max) ) end do debug = .false. xemit_0=1.e-10 f=0.2 lat%n_ele_track = lat%n_ele_track n = lat%n_ele_track ! ntrk = 5 nk1 = ntrk-1 ! if(lat%param%rf_on)ntrk = 7 call mat_make_unit(unit) co_out(n)%vec = co%vec n_det_calc = nonlin%n_det_calc if(n_det_calc == 0) n_det_calc = 1 if (debug) print '(a, 6e12.4)', 'Input CO:', co_out(n)%vec do i=1,n_det_calc if(nonlin%delta_orb(i)%vec(1) == 0.)then if (i == 1) then print *, 'ERROR IN NONLINEAR_CLOSED_ORBIT: NONLIN%DELTA_ORB NOT SET!' call err_exit else nonlin%delta_orb(i)%vec = nonlin%delta_orb(i-1)%vec * f**(i-1) endif endif call get_init_vec(lat, nonlin%delta_orb(i), ntrk, vec_init) xemit=nonlin%delta_orb(i)%vec(1) nonlin%amplitude_x(i) = nonlin%delta_orb(i)%vec(1) nonlin%amplitude_y(i) = nonlin%delta_orb(i)%vec(3) if(xemit > xemit_0)then do j=1,ntrk do k=1,6 t(j)%orb(0)%vec(k) = co_out(n)%vec(k) + vec_init(j,k) if(abs(t(j)%orb(0)%vec(k)) > 1.)then print *, ' NONLINEAR_CLOSED_ORBIT: orbit amp >1. ' print '(6e12.4)', t(j)%orb(0)%vec return endif end do call track_all(lat, t(j)%orb, 0, track_state) if (track_state /= moving_forward$) return end do call find_jacobian(t, ntrk, n, y ) R0(1:nk1,1:nk1) = y(1:nk1,1:nk1) ! call mat_transfer(y, R0, 1, 1, 1, 1,nk1, 7, 6) do j=1,nk1 z(j) = y(j,ntrk) end do diff_mat(1:nk1,1:nk1) = unit(1:nk1,1:nk1) - R0(1:nk1,1:nk1) call mat_inverse(diff_mat(1:nk1,1:nk1), inverse_diff_mat(1:nk1,1:nk1)) co_out(n)%vec(1:nk1) = matmul(inverse_diff_mat(1:nk1,1:nk1), z(1:nk1)) det = determinant (R0(1:nk1,1:nk1)) if(debug)then print *, ' iteration ', i, ' determinant=', det, ' xemit=', xemit print '(4e12.4)', ((R0(j,k),k=1,nk1),j=1,nk1) print *, ' closed orbit ' print '(4e12.4)', co_out(n)%vec endif det_ele(0)%four_by4(i) = det det_ele(0)%two_by2_ul(i) = R0(1,1)*R0(2,2)-R0(1,2)*R0(2,1) det_ele(0)%two_by2_lr(i) = R0(3,3)*R0(4,4)-R0(3,4)*R0(4,3) det_ele(0)%trace(i) = R0(1,1)+R0(2,2)+R0(3,3)+R0(4,4) endif do j=1,n call find_jacobian(t, ntrk, j, y) R(1:nk1,1:nk1) = y(1:nk1,1:nk1) det = determinant (R(1:nk1,1:nk1)) ! full turn at n = R*R0*R^{-1} nt = nk1 call mat_inverse(R(1:nt,1:nt), R_inverse(1:nt,1:nt)) R_ft(1:nk1,1:nk1) = matmul (R(1:nk1,1:nk1), & matmul(R0(1:nk1,1:nk1), R_inverse(1:nk1,1:nk1))) ! call mat_mult (R0, R_inverse, R_temp, nk1, 6) ! call mat_mult (R, R_temp, R_ft, nk1, 6) det_ele(j)%four_by4(i) = det !matrix from start to n det_ele(j)%two_by2_ul(i) = R_ft(1,1)*R_ft(2,2)-R_ft(1,2)*R_ft(2,1) det_ele(j)%two_by2_lr(i) = R_ft(3,3)*R_ft(4,4)-R_ft(3,4)*R_ft(4,3) if (i == n_det_calc)then !last iteration copy R to lat_out co_out(0)%vec = co_out(n)%vec ! call mat_transfer(R, lat_out%ele(j)%mat6, 1, 1, 1, 1, 6, 6, 6) lat_out%ele(j)%mat6(1:6,1:6) = R(1:6,1:6) ! Added the following two lines to keep twiss_propagate1 happy... - mjf, dcs lat_out%ele(j)%mat6(5,5) = 1.0 lat_out%ele(j)%mat6(6,6) = 1.0 do k=1,nk1 z(k) = y(k,ntrk) end do rco(1:nk1) = matmul (lat_out%ele(j)%mat6(1:nk1,1:nk1), & co_out(0)%vec(1:nk1)) ! call mat_vec_mult(lat_out%ele(j)%mat6(1:6,1:6), & ! co_out(0)%vec, rco, 0, nk1, 6) co_out(j)%vec(1:nk1) = z(1:nk1) + rco(1:nk1) endif end do end do !n_det_calc lat_out%ele(0)%mat6(1:nk1,1:nk1) = lat_out%ele(n)%mat6(1:nk1,1:nk1) do i=1,7 deallocate( t(i)%orb ) end do return end subroutine !........................................................................ ! ! Subroutine : NONLINEAR_TWISS(LAT, CO, LAT_OUT, CO_OUT, NONLIN, DET_ELE, track_state) ! ! Description: Subroutine to calcuate the twiss parameters based on jacobian ! computed about the closed orbit including nonlinearities ! ! Arguments : ! Input: ! LAT -- lat_struct: Lat ! CO -- Coord_struct: Starting point of closed orbit ! neglecting nonlinearities or zero ! If there are no synchrotron ! oscillations, the energy offset ! is co.vec(6) ! ! Output: ! ! LAT_OUT -- lat_struct: Transfer matrices about the closed ! orbit and twiss parameters for element i are ! lat_out.ele(i).mat6 and lat_out.ele(i).x.beta ! etc. ! CO_OUT -- Coord_struct: Closed orbit ! NONLIN -- Nonlin_ele_struct: Nonlin ! DET_ELE -- Det_ele_struct : Det_ele determinants at each ! element for each iteration ! ! Mod/Commons: ! ! Calls : ! ! Author : ! ! Modified : ! !........................................................................ subroutine nonlinear_twiss(lat, co, lat_out, co_out, nonlin, det_ele, err_flag, track_state) implicit none type (lat_struct) lat , lat_out type (coord_struct) co type (coord_struct), allocatable :: co_out(:) type (nonlin_ele_struct) nonlin type (det_ele_struct), allocatable :: det_ele(:) type (twiss_struct) twiss1, twiss2 real(rp) :: t0_4(4,4) real(rp) u(4,4), v(4,4), ubar(4,4), vbar(4,4), g(4,4) integer i, status, track_state real(rp) rate1, rate2 logical err_flag ! lat_out = lat call nonlinear_closed_orbit(lat,co, lat_out, co_out,nonlin, det_ele, track_state) if (track_state /= moving_forward$ ) return ! call mat_transfer (lat_out.ele(0).mat6, t0_4, 1, 1, 1, 1, 4, 6, 4) t0_4(1:4,1:4) = lat_out%ele(0)%mat6(1:4, 1:4) twiss1%beta = 0. twiss2%beta = 0. call mat_symp_decouple (t0_4, status, u, v, ubar, vbar, g, twiss1, twiss2, lat_out%ele(0)%gamma_c, .false.) err_flag = (status /= ok$) rate1 = sqrt(max(abs(u(1,1) + u(2,2)) - 2, 0.0)) rate2 = sqrt(max(abs(u(3,3) + u(4,4)) - 2, 0.0)) lat_out%param%unstable_factor = max(rate1, rate2) if(twiss1%beta /= 0. .and. twiss2%beta /= 0.)then lat_out%ele(0)%c_mat = v(1:2,3:4) lat_out%ele(0)%a = twiss1 lat_out%ele(0)%b = twiss2 lat_out%a%tune = twiss1%phi lat_out%b%tune = twiss2%phi if(relative_mode_flip(lat%ele(0),lat_out%ele(0)))then call do_mode_flip (lat_out%ele(0)) lat_out%ele(0)%mode_flip = .false. lat_out%a%tune = twiss2%phi lat_out%b%tune = twiss1%phi endif do i=1,lat%n_ele_track if(lat_out%ele(i)%key == marker$)then call twiss_propagate1 (lat_out%ele(i-1), lat_out%ele(i)) else call twiss_propagate1 (lat_out%ele(0), lat_out%ele(i)) endif end do else lat_out%ele(0)%a%beta = 0. lat_out%ele(0)%b%beta = 0. lat_out%ele(0)%a%alpha = 0. lat_out%ele(0)%b%alpha = 0. lat_out%ele(0)%a%gamma = 0. lat_out%ele(0)%b%gamma = 0. lat_out%a%tune = 0. lat_out%b%tune = 0. lat_out%param%unstable_factor = 1. endif return end subroutine nonlinear_twiss !........................................................................ ! ! Subroutine : NONLINEARITY(LAT, CO, NONLIN, CON) ! ! Description: Compute the nonlinear closed orbit for delta E/E = 0. and ! for an energy slightly above and slightly below 0. For each ! closed orbit compute jacobian and associated twiss parameters ! at the IP and from the IP to the end of each element and then ! construct chromatic function (1/Beta)(d Beta)/(d delta). The ! difference of the off energy closed orbits yields dispersion ! ! Arguments : ! INPUT: ! LAT -- lat_struct: Lat ! CO -- Coord_struct: Co is the starting point for calculating ! the nonlinear closed orbit ! NONLIN -- Nonlin_struct: ! .ENERGY_OFFSET(I) -- The energy offsets for closed orbit. ! .N_DET_CALC -- Number of n_det_calc to compute closed orbit. ! .DELTA_ORB -- Effective emittance for start of tracking to ! find Jacobian. Coord_struct ! CON -- Constraint_struct: ! ! OUTPUT: ! CO -- Coord_struct: Co is overwritten with the on energy ! nonlinear closed orbit ! NONLIN -- Nonlin_struct: Includes energy dependence of twiss parameters, ! determinants with max deviation from 1 ! for all energies ! ! Mod/Commons: ! ! Calls : ! ! Author : ! ! Modified : ! !........................................................................ subroutine nonlinearity(lat, co, nonlin, con, err_flag) implicit none type (lat_struct) lat type (lat_struct), save :: lat_off(3) type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (constraint_struct) con type (coord_struct), allocatable, save :: co_out(:) type (det_ele_struct), allocatable, save :: det_ele(:) type (det_ele_struct), allocatable, save :: det_max(:) real(rp) Qx1, Qx2, Qx3, Qy1, Qy2, Qy3 real(rp) delta integer i, j, k, track_state logical flip1, flip2, flip3, err_flag call reallocate_coord( co_out, lat%n_ele_max ) if ( .not. allocated(det_ele) ) allocate( det_ele(0:lat%n_ele_max) ) if ( .not. allocated(det_max) ) allocate( det_max(0:lat%n_ele_max) ) ! ! do k = 1,3 ! if ( .not. associated( lat_off(k)%ele ) ) then ! call init_lat( lat_off(k), lat%n_ele_max ) ! endif ! end do ! save_sym = lat%param%symmetry ! symmetry = no_symmetry$ ! call set_symmetry(symmetry, lat) call set_on_off (rfcavity$, lat, off$, co) do k = 0,lat%n_ele_track do j = 1,nonlin%n_det_calc det_max(k)%four_by4(j) =0. det_max(k)%two_by2_ul(j) =0. det_max(k)%two_by2_lr(j) =0. end do end do do i=1,nonlin%energies co(0)%vec(6) = nonlin%energy_offset(i) call nonlinear_twiss(lat, co(0), lat_off(i), co_out, nonlin, det_ele, err_flag, track_state) if(track_state /= moving_forward$)then nonlin%chrom_x = 999. nonlin%chrom_y = 999. return endif ! find max deviation of determinant from 1 for each energy and stash do k = 0,lat%n_ele_track nonlin%non_ele(k)%co_off(i) = co_out(k) do j = 1,nonlin%n_det_calc det_max(k)%four_by4(j) = & max(det_max(k)%four_by4(j), abs(1.-det_ele(k)%four_by4(j))) det_max(k)%two_by2_ul(j) = & max(det_max(k)%two_by2_ul(j), abs(1.-det_ele(k)%two_by2_ul(j))) det_max(k)%two_by2_lr(j) = & max(det_max(k)%two_by2_lr(j), abs(1.-det_ele(k)%two_by2_lr(j))) end do end do if(lat_off(i)%param%unstable_factor > 0.)then print '(a,1x,a,es12.4,a,4es12.4)',' NONLINEAR_TWISS finds unstable lattice ',& 'energy offset = ', nonlin%energy_offset(i), & ' disp vec = ', nonlin%delta_orb(nonlin%n_det_calc)%vec(1:2) , nonlin%delta_orb(nonlin%n_det_calc)%vec(3:4) nonlin%chrom_x = 99. nonlin%chrom_y = 99. return endif end do delta = nonlin%energy_offset(3)-nonlin%energy_offset(1) do i=0,lat%n_ele_track do k = 1,nonlin%n_det_calc nonlin%non_ele(i)%det%four_by4(k) = det_max(i)%four_by4(k) nonlin%non_ele(i)%det%two_by2_ul(k) = det_max(i)%two_by2_ul(k) nonlin%non_ele(i)%det%two_by2_lr(k) = det_max(i)%two_by2_lr(k) end do flip1 = lat_off(1)%ele(i)%mode_flip flip2 = lat_off(2)%ele(i)%mode_flip flip3 = lat_off(3)%ele(i)%mode_flip if ((flip1 .eqv. .not. flip2) .or. (flip2 .eqv. .not. flip3)) then nonlin%non_ele(i)%a%dbeta = 0 nonlin%non_ele(i)%b%dbeta = 0 nonlin%non_ele(i)%a%dalpha = 0 nonlin%non_ele(i)%b%dalpha = 0 nonlin%non_ele(i)%a%dphi = 0 nonlin%non_ele(i)%b%dphi = 0 nonlin%non_ele(i)%a%eta = 0 nonlin%non_ele(i)%b%eta = 0 nonlin%non_ele(i)%a%etap = 0 nonlin%non_ele(i)%b%etap = 0 else nonlin%non_ele(i)%a%dbeta =(lat_off(3)%ele(i)%a%beta & -lat_off(1)%ele(i)%a%beta)/delta /lat_off(2)%ele(i)%a%beta nonlin%non_ele(i)%b%dbeta =(lat_off(3)%ele(i)%b%beta & -lat_off(1)%ele(i)%b%beta)/delta /lat_off(2)%ele(i)%b%beta nonlin%non_ele(i)%a%dalpha =(lat_off(3)%ele(i)%a%alpha & -lat_off(1)%ele(i)%a%alpha)/delta nonlin%non_ele(i)%b%dalpha =(lat_off(3)%ele(i)%b%alpha & -lat_off(1)%ele(i)%b%alpha)/delta nonlin%non_ele(i)%a%dphi =(lat_off(3)%ele(i)%a%phi & -lat_off(1)%ele(i)%a%phi)/delta nonlin%non_ele(i)%b%dphi =(lat_off(3)%ele(i)%b%phi & -lat_off(1)%ele(i)%b%phi)/delta nonlin%non_ele(i)%a%eta = (nonlin%non_ele(i)%co_off(1)%vec(1) & - nonlin%non_ele(i)%co_off(3)%vec(1))/delta nonlin%non_ele(i)%a%etap = (nonlin%non_ele(i)%co_off(1)%vec(2) & - nonlin%non_ele(i)%co_off(3)%vec(2))/delta nonlin%non_ele(i)%b%eta = (nonlin%non_ele(i)%co_off(3)%vec(3) & - nonlin%non_ele(i)%co_off(1)%vec(3))/delta nonlin%non_ele(i)%b%etap = (nonlin%non_ele(i)%co_off(3)%vec(4) & - nonlin%non_ele(i)%co_off(1)%vec(4))/delta endif nonlin%non_ele(i)%a%beta = lat_off(2)%ele(i)%a%beta nonlin%non_ele(i)%b%beta = lat_off(2)%ele(i)%b%beta nonlin%non_ele(i)%a%alpha = lat_off(2)%ele(i)%a%alpha nonlin%non_ele(i)%b%alpha = lat_off(2)%ele(i)%b%alpha nonlin%non_ele(i)%a%phi = lat_off(2)%ele(i)%a%phi nonlin%non_ele(i)%b%phi = lat_off(2)%ele(i)%b%phi co(i) = nonlin%non_ele(i)%co_off(2) end do Qx3= lat_off(3)%a%tune/twopi if(Qx3<0.)Qx3=Qx3+1. Qx2= lat_off(2)%a%tune/twopi if(Qx2<0.)Qx2=Qx2+1. Qx1= lat_off(1)%a%tune/twopi if(Qx1<0.)Qx1=Qx1+1. Qy1= lat_off(1)%b%tune/twopi if(Qy1<0.)Qy1=Qy1+1. Qy2= lat_off(2)%b%tune/twopi if(Qy2<0.)Qy2=Qy2+1. Qy3= lat_off(3)%b%tune/twopi if(Qy3<0.)Qy3=Qy3+1. nonlin%tune_x = Qx2 nonlin%tune_y = Qy2 nonlin%chrom_x = (Qx3-Qx1)/delta nonlin%chrom_y = (Qy3-Qy1)/delta nonlin%dchrom_x = (Qx3-Qx2)/nonlin%energy_offset(3) - & (Qx1-Qx2)/nonlin%energy_offset(1) nonlin%dchrom_y = (Qy3-Qy2)/nonlin%energy_offset(3) - & (Qy1-Qy2)/nonlin%energy_offset(1) nonlin%trace(1:nonlin%n_det_calc) = det_ele(0)%trace(1:nonlin%n_det_calc) ! call set_symmetry(save_sym, lat) ! deallocate( det_ele ) return end subroutine !........................................................................ ! ! Subroutine : SEXTUPOLE_MOMENTS(LAT, CO, QUAD, CON, MOMENT) ! ! Description: Compute first order phase space distortion due to sextupoles. ! The scalar product of the moments vector (the length of the ! vector is just the number of sextupoles), and the sextupole ! strength vector gives the distortions. See Tom Pelaia's thesis ! page 21. There are ten such moments Also compute first order ! energy and pretzel dependence of beta and tunes as vectors ! that multiply the sextupole strength vector. All quantities ! are linear in the sextupole strengths. ! ! Arguments : ! INPUT: ! LAT -- lat_struct: Lat - This is presumably the lat structure ! containing twiss parameters and eta derived from ! matrix multiplication ! CO(0:*) -- Coord_struct: Co is the closed orbit ! QUAD -- zquad_struct: list of quads and average beta ! CON -- Constraint_struct: ! ! OUTPUT: ! MOMENT -- Moment_struct : Moment - Sextupole moments ! ! Mod/Commons: ! ! Calls : ! ! Author : ! ! Modified : ! !----------------------------------------------------------------- ! exact calculation ! find twiss parameters with an energy offset ! find dbeta/dsextupole ! extrapolate to zero sextupole strength ! !........................................................................ subroutine sextupole_moments(lat, co, quad, con, moment) implicit none type (lat_struct), target :: lat type (lat_struct), save :: lat1, lat2 type (coord_struct), allocatable :: co(:) type (coord_struct), allocatable, save :: orb1(:), orb2(:), co_opposite(:) type (zquad_struct) quad type (moment_struct) moment type (constraint_struct) con type (ele_struct), pointer :: e0, e1, e2 type (ele_struct) ave real(rp) pi_nu1, pi_nu2, sin_2pi_nu1, sin_2pi_nu2, sin_pi_nu1, sin_3pi_nu1 real(rp) denom_plus2, denom_minus2 real(rp) beta1(275), beta2(275), root_beta1, root_beta2 real(rp) eta1(275), phi1(275), phi2(275) real(rp) x(275), k2_new, k2_old real(rp) cos_stuff_1, phi0_1, cos_stuff_2, phi0_2, arg1, arg2 real(rp) four_pi real(rp) beta1_0, beta2_0, phiq1, phiq2, del_x, del_y, k2l real(rp) ks, rb1_3, rb1_b2, p1, p2 real(rp) ampx_cos, ampx_sin, ampy_cos, ampy_sin real(rp) :: amp = 1 real(rp) xi_plus, xi_minus, cos_xi_sum,sin_xi_sum, cos_xi_dif, sin_xi_dif real(rp) pi_delta_nu, pi_sigma_nu, sin_pi_delta, sin_pi_sigma real(rp) trM_N, del_k2l, del_k2 real(rp) root_betax, root_betay real(rp) A(2,2), A_inv(2,2), det_A real(rp) ks_x, ks_y real(rp) k_plus, k_minus integer k,i, ix, j integer n_twopi logical :: nominal = .true. ! call reallocate_coord( orb1, lat%n_ele_max ) call reallocate_coord( orb2, lat%n_ele_max ) call reallocate_coord( co_opposite, lat%n_ele_max ) four_pi = 2 * twopi pi_nu1 = lat%a%tune/2 pi_nu2 = lat%b%tune/2 sin_pi_nu1 = sin(pi_nu1) sin_3pi_nu1 = sin(3*pi_nu1) denom_plus2 = sin(pi_nu1 + 2*pi_nu2) denom_minus2 = sin(pi_nu1 - 2*pi_nu2) sin_2pi_nu1 = sin(2*pi_nu1) sin_2pi_nu2 = sin(2*pi_nu2) pi_delta_nu = pi_nu1 - pi_nu2 pi_sigma_nu = pi_nu1 + pi_nu2 sin_pi_delta = sin(pi_delta_nu) sin_pi_sigma = sin(pi_sigma_nu) trM_N = 2*(cos(2*pi_nu1) - cos(2*pi_nu2)) moment%chrom_x_0 = 0. moment%chrom_y_0 = 0. moment%chrom_x_ir = 0. moment%chrom_y_ir = 0. i=0 A(1,1)=0. A(1,2)=0. A(2,1)=0. A(2,2)=0. do k=0,lat%n_ele_track co_opposite(k)%vec(1:2) = -co(k)%vec(1:2) co_opposite(k)%vec(3:6) = co(k)%vec(3:6) end do do k = 1, lat%n_ele_track moment%dbeta_x_0(k) = 0. moment%dbeta_y_0(k) = 0. if((lat%ele(k)%key == sextupole$ .and. & lat%ele(k)%value(tilt$) == 0.) .or. lat%ele(k)%key == wiggler$ .or. & lat%ele(k)%key == ab_multipole$)then i=i+1 moment%name(i) = lat%ele(k)%name moment%ix(i) = k beta1(i) = 0.5*(lat%ele(k)%a%beta + lat%ele(k-1)%a%beta) root_beta1 = sqrt(beta1(i)) beta2(i) = 0.5*(lat%ele(k)%b%beta + lat%ele(k-1)%b%beta) root_beta2 = sqrt(beta2(i)) phi1(i) = 0.5*(lat%ele(k)%a%phi + lat%ele(k-1)%a%phi) phi2(i) = 0.5*(lat%ele(k)%b%phi + lat%ele(k-1)%b%phi) eta1(i) = 0.5*(lat%ele(k)%a%eta + lat%ele(k-1)%a%eta) x(i) = 0.5*(co(k)%vec(1) + co(k-1)%vec(1)) p1=phi1(i) p2=phi2(i) rb1_3 = root_beta1**3 rb1_b2 = root_beta1*beta2(i) if(lat%ele(k)%key == sextupole$)then moment%k2l(i) = lat%ele(k)%value(k2$) * lat%ele(k)%value(l$) else if(x(i) /= 0)then call lat_make_mat6(lat, k, co) k_plus = lat%ele(k)%mat6(2,1) call lat_make_mat6(lat, k, co_opposite) k_minus = lat%ele(k)%mat6(2,1) moment%k2l(i) = -(k_plus-k_minus)/2/x(i) else moment%k2l(i) = 0. endif endif moment%m1(i) = rb1_3 * sin(p1) / (4 * sin_pi_nu1) moment%m2(i) = rb1_3 * sin(3*p1) / (4 * sin_3pi_nu1) moment%m3(i) = rb1_3 * cos(p1) / (4 * sin_pi_nu1) moment%m4(i) = rb1_3 * cos(3*p1) / (4 * sin_3pi_nu1) moment%m5(i) = rb1_b2 * sin(p1) / (2 * sin_pi_nu1) moment%m6(i) = rb1_b2 * sin(p1 + 2*p2) / (4 * denom_plus2) moment%m7(i) = rb1_b2 * sin(p1 - 2*p2) / (4 * denom_minus2) moment%m8(i) = rb1_b2 * cos(p1) / (2 * sin_pi_nu1) moment%m9(i) = rb1_b2 * cos(p1 + 2*p2) / (4 * denom_plus2) moment%m10(i) = rb1_b2 * cos(p1 - 2*p2) / (4 * denom_minus2) moment%chrom_x(i) = beta1(i) * eta1(i) / four_pi moment%chrom_y(i) = -beta2(i) * eta1(i) / four_pi if(beta1(i) > beta2(i))then A(1,1)=A(1,1)+moment%chrom_x(i) A(2,1)=A(2,1)+moment%chrom_y(i) else A(1,2)=A(1,2)+moment%chrom_x(i) A(2,2)=A(2,2)+moment%chrom_y(i) endif moment%tonality_x(i) = beta1(i) * x(i)/four_pi moment%tonality_y(i) = -beta2(i) * x(i)/four_pi moment%n_sext = i endif !sextupole end do ! do j=1,quad%n !chromaticity i=quad%lens(j)%ix call twiss_at_element (lat%ele(i), average = ave) quad%lens(j)%a%beta_ave = ave%a%beta quad%lens(j)%b%beta_ave = ave%b%beta beta1_0 = quad%lens(j)%a%beta_ave beta2_0 = quad%lens(j)%b%beta_ave ks = lat%ele(i)%value(k1$) * lat%ele(i)%value(l$) if( lat%ele(i)%key == wiggler$ )then ks_x = -lat%ele(i)%mat6(2,1) moment%chrom_x_0 = moment%chrom_x_0 - beta1_0 * ks_x / four_pi ks_y = lat%ele(i)%mat6(4,3) moment%chrom_y_0 = moment%chrom_y_0 + beta2_0 * ks_y / four_pi else moment%chrom_x_0 = moment%chrom_x_0 - beta1_0 * ks / four_pi moment%chrom_y_0 = moment%chrom_y_0 + beta2_0 * ks / four_pi if(lat%ele(i)%s < lat%ele(lat%n_ele_track)%s/12. .or. & lat%ele(i)%s > 11./12.*(lat%ele(lat%n_ele_max)%s ))then moment%chrom_x_ir = moment%chrom_x_ir - beta1_0 * ks / four_pi moment%chrom_y_ir = moment%chrom_y_ir + beta2_0 * ks / four_pi endif endif end do root_betay = sqrt(lat%ele(0)%b%beta) root_betax = sqrt(lat%ele(0)%a%beta) det_A = A(1,1)*A(2,2)-A(1,2)*A(2,1) A_inv(1,1)=A(2,2) A_inv(2,2)=A(1,1) A_inv(1,2)=-A(1,2) A_inv(2,1)=-A(2,1) moment%s_h= (-A_inv(1,1)*moment%chrom_x_0 - & A_inv(1,2)*moment%chrom_y_0)/det_A moment%s_v= (-A_inv(2,1)*moment%chrom_x_0 - & A_inv(2,2)*moment%chrom_y_0)/det_A do k = 1, lat%n_ele_track phi0_1 = 0.5*(lat%ele(k)%a%phi + lat%ele(k-1)%a%phi) phi0_2 = 0.5*(lat%ele(k)%b%phi + lat%ele(k-1)%b%phi) do j=1,quad%n beta1_0 = quad%lens(j)%a%beta_ave beta2_0 = quad%lens(j)%b%beta_ave i=quad%lens(j)%ix ks = lat%ele(i)%value(k1$) * lat%ele(i)%value(l$) phiq1 = 0.5*(lat%ele(i)%a%phi + lat%ele(i-1)%a%phi) phiq2 = 0.5*(lat%ele(i)%b%phi + lat%ele(i-1)%b%phi) arg1 = phi0_1-phiq1 if(arg1 < 0.)arg1 = arg1 + 2*pi_nu1 cos_stuff_1 = cos(2*(arg1 - pi_nu1)) arg2 = phi0_2-phiq2 if(arg2 < 0.)arg2 = arg2 + 2*pi_nu2 cos_stuff_2 = cos(2*(arg2 - pi_nu2)) if( lat%ele(i)%key == wiggler$ )then ks_x = -lat%ele(i)%mat6(2,1) moment%dbeta_x_0(k) = moment%dbeta_x_0(k)+ & beta1_0 * ks_x*cos_stuff_1/2/sin_2pi_nu1 ks_y = lat%ele(i)%mat6(4,3) moment%dbeta_y_0(k) = moment%dbeta_y_0(k) & - beta2_0 * ks_y*cos_stuff_2/2/sin_2pi_nu2 else moment%dbeta_x_0(k) = moment%dbeta_x_0(k)+ & beta1_0 * ks*cos_stuff_1/2/sin_2pi_nu1 moment%dbeta_y_0(k) = moment%dbeta_y_0(k) & - beta2_0 * ks*cos_stuff_2/2/sin_2pi_nu2 endif end do moment%dbeta_x_0_norm(k) = moment%dbeta_x_0(k) * sin_2pi_nu1 moment%dbeta_y_0_norm(k) = moment%dbeta_y_0(k) * sin_2pi_nu2 do i=1,moment%n_sext arg1 = phi0_1-phi1(i) if(arg1 < 0.)arg1 = arg1 + 2*pi_nu1 cos_stuff_1 = cos(2*(arg1 - pi_nu1)) arg2 = phi0_2-phi2(i) if(arg2 < 0.)arg2 = arg2 + 2*pi_nu2 cos_stuff_2 = cos(2*(arg2 - pi_nu2)) xi_plus = arg1 + arg2 xi_minus = arg1 - arg2 cos_xi_sum = cos(xi_plus+pi_sigma_nu) sin_xi_sum = sin(xi_plus+pi_sigma_nu) cos_xi_dif = cos(xi_minus+pi_delta_nu) sin_xi_dif = sin(xi_minus+pi_delta_nu) ampx_cos = sqrt(beta1(i)) * amp * (cos(phi1(i))) ampx_sin = sqrt(beta1(i)) * amp * (sin(phi1(i))) ampy_cos = sqrt(beta2(i)) * amp * (cos(phi2(i))) ampy_sin = sqrt(beta2(i)) * amp * (sin(phi2(i))) moment%dbeta_x(k,i) = -beta1(i) * eta1(i)*cos_stuff_1/2/sin_2pi_nu1 moment%dbeta_y(k,i) = beta2(i) * eta1(i)*cos_stuff_2/2/sin_2pi_nu2 moment%dbeta_dpretz_x(k,i) = -beta1(i) * x(i) * cos_stuff_1 & /2/sin_2pi_nu1 moment%dbeta_dpretz_y(k,i) = beta2(i) * x(i) * cos_stuff_2 & /2/sin_2pi_nu2 moment%dbeta_dcos_x(k,i) = -beta1(i) * ampx_cos * cos_stuff_1 & /2/sin_2pi_nu1 moment%dbeta_dcos_y(k,i) = beta2(i) * ampx_cos * cos_stuff_2 & /2/sin_2pi_nu2 moment%dbeta_dsin_x(k,i) = -beta1(i) * ampx_sin * cos_stuff_1 & /2/sin_2pi_nu1 moment%dbeta_dsin_y(k,i) = beta2(i) * ampx_sin * cos_stuff_2 & /2/sin_2pi_nu2 moment%a_real(k,i) = sqrt(beta1(i)*beta2(i))* 2 * sin_pi_delta & * ampy_sin * cos_xi_sum /trM_N moment%a_image(k,i) = -sqrt(beta1(i)*beta2(i))* 2 * sin_pi_delta & * ampy_sin * sin_xi_sum /trM_N moment%b_real(k,i) = sqrt(beta1(i)*beta2(i))* 2 * sin_pi_sigma & * ampy_sin * cos_xi_dif /trM_N moment%b_image(k,i) = -sqrt(beta1(i)*beta2(i))* 2 * sin_pi_sigma & * ampy_sin * sin_xi_dif /trM_N end do end do !----------------------------------------------------------------- ! exact calculation if (con%dbeta_exact) then lat1 = lat where (lat1%ele%key == elseparator$) lat1%ele%value(hkick$) = 0 lat1%ele%value(vkick$) = 0 end where ! find twiss parameters with an energy offset call set_on_off (rfcavity$, lat1, off$, co) orb1(:)%vec(6) = -con%delta_e orb1(:)%vec(1) = -con%delta_e*lat1%ele(:)%a%eta orb1(:)%vec(2) = -con%delta_e*lat1%ele(:)%a%etap orb1(:)%vec(3) = -con%delta_e*lat1%ele(:)%b%eta orb1(:)%vec(4) = -con%delta_e*lat1%ele(:)%b%etap call lat_make_mat6(lat1, -1, orb1) call twiss_at_start (lat1) call twiss_propagate_all (lat1) lat2 = lat1 forall (i = 0:lat%n_ele_max) orb2(i)%vec = -orb1(i)%vec call lat_make_mat6(lat2, -1, orb2) call twiss_at_start (lat2) call twiss_propagate_all (lat2) do k = 1, lat%n_ele_track e0 => lat%ele(k) e1 => lat1%ele(k) e2 => lat2%ele(k) moment%dbeta_x_0(k) = & (e2%gamma_c**2 * e2%a%beta - e1%gamma_c**2 * e1%a%beta) / & (2 * e0%gamma_c**2 * e0%a%beta * con%delta_e) moment%dbeta_y_0(k) = & (e2%gamma_c**2 * e2%b%beta - e1%gamma_c**2 * e1%b%beta) / & (2 * e0%gamma_c**2 * e0%b%beta * con%delta_e) moment%dc_mat(1,1)%ele(k)%dvar0 = (e2%c_mat(1,1) - e1%c_mat(1,1)) / (2 * con%delta_e) moment%dc_mat(1,2)%ele(k)%dvar0 = (e2%c_mat(1,2) - e1%c_mat(1,2)) / (2 * con%delta_e) moment%dc_mat(2,1)%ele(k)%dvar0 = (e2%c_mat(2,1) - e1%c_mat(2,1)) / (2 * con%delta_e) moment%dc_mat(2,2)%ele(k)%dvar0 = (e2%c_mat(2,2) - e1%c_mat(2,2)) / (2 * con%delta_e) enddo moment%chrom_x_0 = (lat2%a%tune - lat1%a%tune) / (2 * twopi * con%delta_e) moment%chrom_y_0 = (lat2%b%tune - lat1%b%tune) / (2 * twopi * con%delta_e) ! find dbeta/dsextupole do i = 1, moment%n_sext ix = moment%ix(i) if (lat%ele(ix)%value(k2$) < 0) then del_k2l = -0.1 else del_k2l = 0.1 endif del_k2 = del_k2l / lat%ele(ix)%value(l$) k2_old = lat%ele(ix)%value(k2$) k2_new = lat%ele(ix)%value(k2$) + del_k2 lat1%ele(ix)%value(k2$) = k2_new call lat_make_mat6(lat1, ix, orb1) call twiss_at_start (lat1) call twiss_propagate_all (lat1) lat2%ele(ix)%value(k2$) = k2_new call lat_make_mat6(lat2, ix, orb2) call twiss_at_start (lat2) call twiss_propagate_all (lat2) do k = 1, lat%n_ele_track e0 => lat%ele(k) e1 => lat1%ele(k) e2 => lat2%ele(k) moment%dbeta_x(k,i) = & (e2%gamma_c**2 * e2%a%beta - e1%gamma_c**2 * e1%a%beta) / & (2 * e0%gamma_c**2 * e0%a%beta * con%delta_e) moment%dbeta_y(k,i) = & (e2%gamma_c**2 * e2%b%beta - e1%gamma_c**2 * e1%b%beta) / & (2 * e0%gamma_c**2 * e0%b%beta * con%delta_e) moment%dc_mat(1,1)%ele(k)%dvar(i) = (e2%c_mat(1,1) - e1%c_mat(1,1)) / (2 * con%delta_e) moment%dc_mat(1,2)%ele(k)%dvar(i) = (e2%c_mat(1,2) - e1%c_mat(1,2)) / (2 * con%delta_e) moment%dc_mat(2,1)%ele(k)%dvar(i) = (e2%c_mat(2,1) - e1%c_mat(2,1)) / (2 * con%delta_e) moment%dc_mat(2,2)%ele(k)%dvar(i) = (e2%c_mat(2,2) - e1%c_mat(2,2)) / (2 * con%delta_e) moment%dc_mat(1,1)%ele(k)%dvar(i) = (moment%dc_mat(1,1)%ele(k)%dvar(i) - & moment%dc_mat(1,1)%ele(k)%dvar0) / del_k2l moment%dc_mat(1,2)%ele(k)%dvar(i) = (moment%dc_mat(1,2)%ele(k)%dvar(i) - & moment%dc_mat(1,2)%ele(k)%dvar0) / del_k2l moment%dc_mat(2,1)%ele(k)%dvar(i) = (moment%dc_mat(2,1)%ele(k)%dvar(i) - & moment%dc_mat(2,1)%ele(k)%dvar0) / del_k2l moment%dc_mat(2,2)%ele(k)%dvar(i) = (moment%dc_mat(2,2)%ele(k)%dvar(i) - & moment%dc_mat(2,2)%ele(k)%dvar0) / del_k2l enddo del_x = (lat2%a%tune - lat1%a%tune) / (2 * twopi * con%delta_e) del_y = (lat2%b%tune - lat1%b%tune) / (2 * twopi * con%delta_e) moment%chrom_x(i) = (del_x - moment%chrom_x_0) / del_k2l moment%chrom_y(i) = (del_y - moment%chrom_y_0) / del_k2l moment%dbeta_x(:,i) = (moment%dbeta_x(:,i) - moment%dbeta_x_0) / del_k2l moment%dbeta_y(:,i) = (moment%dbeta_y(:,i) - moment%dbeta_y_0) / del_k2l lat1%ele(ix)%value(k2$) = k2_old lat2%ele(ix)%value(k2$) = k2_old call lat_make_mat6(lat1, ix, orb1) call lat_make_mat6(lat2, ix, orb2) enddo ! extrapolate to zero sextupole strength do i = 1, moment%n_sext ix = moment%ix(i) k2l = lat%ele(ix)%value(k2$) * lat%ele(ix)%value(l$) moment%dbeta_x_0(:) = moment%dbeta_x_0(:) - moment%dbeta_x(:,i) * k2l moment%dbeta_y_0(:) = moment%dbeta_y_0(:) - moment%dbeta_y(:,i) * k2l moment%chrom_x_0 = moment%chrom_x_0 - moment%chrom_x(i) * k2l moment%chrom_y_0 = moment%chrom_y_0 - moment%chrom_y(i) * k2l enddo endif return end subroutine sextupole_moments !........................................................................ ! ! Subroutine : WRITE_MOMENTS(LAT, CO, NONLIN, MOMENT) ! ! Description: Write coordinate of trajectory at the IP and scalar product of ! moments and sextupole strength vector to a file ! ! Arguments : ! INPUT: ! LAT -- lat_struct : Lat ! MOMENT -- Moment_struct: moment ! CO -- Coord_struct: Co is the closed orbit ! ! Mod/Commons: ! ! Calls : ! ! Author : ! ! Modified : ! ! ! call beta_v_amplitude (lat, nonlin) ! endif !----------------------------------------------------------------------- subroutine write_moments(lat, co, nonlin, moment) implicit none type (lat_struct) lat type (coord_struct), allocatable :: co(:) type (moment_struct) moment type (nonlin_ele_struct) nonlin real(rp) dby_pretz, dbx_pretz, sum1, sum2, sum3, sum4, sum5, sum6 real(rp) sum7, sum8, sum9, sum10 integer k, n integer m,mp logical :: debug = .true. ! n=moment%n_sext sum1 = dot_product(moment%m1 , moment%k2l) sum2 = dot_product(moment%m2 , moment%k2l) sum3 = dot_product(moment%m3 , moment%k2l) sum4 = dot_product(moment%m4 , moment%k2l) sum5 = dot_product(moment%m5 , moment%k2l) sum6 = dot_product(moment%m6 , moment%k2l) sum7 = dot_product(moment%m7 , moment%k2l) sum8 = dot_product(moment%m8 , moment%k2l) sum9 = dot_product(moment%m9 , moment%k2l) sum10 = dot_product(moment%m1, moment%k2l) moment%mchrom_x = dot_product(moment%chrom_x , moment%k2l) + & moment%chrom_x_0 moment%mchrom_y = dot_product(moment%chrom_y , moment%k2l) + & moment%chrom_y_0 moment%mtonality_x = dot_product(moment%tonality_x , moment%k2l) moment%mtonality_y = dot_product(moment%tonality_y , moment%k2l) moment%mtonality_x_west = & dot_product(moment%tonality_x(1:n/2), moment%k2l(1:n/2)) moment%mtonality_y_west = & dot_product(moment%tonality_y(1:n/2), moment%k2l(1:n/2)) nonlin%tonality_x=moment%mtonality_x * 390.1 nonlin%tonality_y=moment%mtonality_y * 390.1 m=lat%n_ele_track mp=lat%n_ele_track moment%mdbeta_x(1:m) = matmul(moment%dbeta_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_y(1:m) = matmul(moment%dbeta_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dpretz_x(1:m) = matmul(moment%dbeta_dpretz_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dpretz_y(1:m) = matmul(moment%dbeta_dpretz_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dcos_x(1:m) = matmul(moment%dbeta_dcos_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dcos_y(1:m) = matmul(moment%dbeta_dcos_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dsin_x(1:m) = matmul(moment%dbeta_dsin_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dsin_y(1:m) = matmul(moment%dbeta_dsin_y(1:m,1:n), moment%k2l(1:n)) moment%ma_real(1:m) = matmul(moment%a_real(1:m,1:n), moment%k2l(1:n)) moment%ma_image(1:m) = matmul(moment%a_image(1:m,1:n), moment%k2l(1:n)) moment%mb_real(1:m) = matmul(moment%b_real(1:m,1:n), moment%k2l(1:n)) moment%mb_image(1:m) = matmul(moment%b_image(1:m,1:n), moment%k2l(1:n)) if(debug)then open(unit=21, file = 'pretz_energy_moments.dat') write(21, 7)co(0)%vec 7 format('!',6f12.6) write(21, 3) open(unit=24, file = 'dbeta_dcos_moments.dat') write(24, 7)co(0)%vec write(24, 33) open(unit=25, file = 'coupling_moments.dat') write(25, 7)co(0)%vec write(25, 34) open(unit=26) ! call beta_v_amplitude (lat, nonlin) open(unit=27, file = 'dbeta_dcos_calc.dat') write(27, 7)co(0)%vec write(27, '(7a)')'!',' ELE# ',' ELEMENT NAME ',' dbeta_x/dcos ', & ' dbeta_y/dcos ',' dbeta_x/dsin ',' dbeta_y/dsin ' open (unit=23, file = 'pretz_energy_calc.dat' ) write(23, 7)co(0)%vec write(23, 8) 3 format('!',' ELE# ',' ELEMENT NAME ',' dbeta_x/ddelta ', & ' dbeta_y/ddelta ',' dbeta_x/dpretz ',' dbeta_y/dpretz') 33 format('!',' ELE# ',' ELEMENT NAME ',' dbeta_dcos_x ', & ' dbeta_dcos_y ',' dbeta_dsin_x ',' dbeta_dsin_y') 34 format('!',' ELE# ',' ELEMENT NAME ',' a_real ', & ' a_imaginary ',' b_real ',' b_imaginary ') 8 format('!',' ELE# ',' ELEMENT NAME ',' DBETA_X/DDELTA ', & ' DBETA_Y/DDELTA ',' DBETA_X/DPRETZ ',' DBETA_Y/DPRETZ') do k = 1, lat%n_ele_track dbx_pretz = (nonlin%non_ele(k)%a%beta - lat%ele(k)%a%beta)/ & lat%ele(k)%a%beta dby_pretz = (nonlin%non_ele(k)%b%beta - lat%ele(k)%b%beta)/ & lat%ele(k)%b%beta write(21, 1)k, lat%ele(k)%name, moment%mdbeta_x(k)+ & moment%dbeta_x_0(k), & moment%mdbeta_y(k)+moment%dbeta_y_0(k), & moment%mdbeta_dpretz_x(k), & moment%mdbeta_dpretz_y(k) write(24, 1)k, lat%ele(k)%name, moment%mdbeta_dcos_x(k), & moment%mdbeta_dcos_y(k), & moment%mdbeta_dsin_x(k), & moment%mdbeta_dsin_y(k) write(25, 1)k, lat%ele(k)%name, moment%ma_real(k), & moment%ma_image(k), & moment%mb_real(k), & moment%mb_image(k) write(23, 1)k, lat%ele(k)%name, & nonlin%non_ele(k)%a%dbeta, & nonlin%non_ele(k)%b%dbeta, & dbx_pretz, dby_pretz write(27, 1)k, lat%ele(k)%name, & nonlin%non_ele(k)%a%dbeta_dcos, & nonlin%non_ele(k)%b%dbeta_dcos, & nonlin%non_ele(k)%a%dbeta_dsin, & nonlin%non_ele(k)%b%dbeta_dsin 1 format(1x,i4,1x,a16, 4(f12.5,4x)) write(26,'(a, 4f12.4)')lat%ele(k)%name, moment%dbeta_x_0(k), & nonlin%non_ele(k)%a%dbeta, moment%dbeta_y_0(k), & nonlin%non_ele(k)%b%dbeta ! endif end do close(unit=21) close(unit=23) close(unit=24) close(unit=25) close(unit=26) close(unit=27) write(22, 7)co(0)%vec write(22, 2) moment%mchrom_x, moment%mchrom_y 2 format(' chrom_x = ',f8.3,5x,'chrom_y = ',f8.3) write(22,4) moment%mtonality_x*390.1, moment%mtonality_y*390.1 4 format(' tonality_x(kHz) = ', f8.3,5x,'tonality_y(kHz) = ',f8.3) write(22,9) moment%mtonality_x_west*390.1, & moment%mtonality_y_west * 390.1 9 format(' tonality_x_west(kHz) = ', f8.3,5x,'tonality_y_west(kHz) = ',f8.3) write(22,5)sum1,sum2,sum3,sum4, sum5 5 format(' M1 = ',e12.4,3x,'M2 = ',e12.4,3x,'M3 = ',e12.4,3x,'M4 = ',e12.4, & 3x,'M5 = ',e12.4) write(22,6)sum6,sum7,sum8,sum9, sum10 6 format(' M6 = ',e12.4,3x,'M7 = ',e12.4,3x,'M8 = ',e12.4,3x,'M9 = ',e12.4, & 3x,'M10 = ',e12.4) endif !(debug) end subroutine write_moments !------------------------------------------------------ !------------------------------------------------------- !------------------------------------------------------ subroutine write_moments_debug(lat, co, nonlin, moment) implicit none type (lat_struct) lat type (coord_struct), allocatable :: co(:) type (moment_struct) moment type (nonlin_ele_struct) nonlin real(rp) dby_pretz, dbx_pretz, sum1, sum2, sum3, sum4, sum5, sum6 real(rp) sum7, sum8, sum9, sum10 integer k, n integer m,mp logical :: debug = .true. ! n=moment%n_sext sum1 = dot_product(moment%m1 , moment%k2l) sum2 = dot_product(moment%m2 , moment%k2l) sum3 = dot_product(moment%m3 , moment%k2l) sum4 = dot_product(moment%m4 , moment%k2l) sum5 = dot_product(moment%m5 , moment%k2l) sum6 = dot_product(moment%m6 , moment%k2l) sum7 = dot_product(moment%m7 , moment%k2l) sum8 = dot_product(moment%m8 , moment%k2l) sum9 = dot_product(moment%m9 , moment%k2l) sum10 = dot_product(moment%m1, moment%k2l) moment%mchrom_x = dot_product(moment%chrom_x , moment%k2l) + & moment%chrom_x_0 moment%mchrom_y = dot_product(moment%chrom_y , moment%k2l) + & moment%chrom_y_0 moment%mtonality_x = dot_product(moment%tonality_x , moment%k2l) moment%mtonality_y = dot_product(moment%tonality_y , moment%k2l) moment%mtonality_x_west = & dot_product(moment%tonality_x(1:n/2), moment%k2l(1:n/2)) moment%mtonality_y_west = & dot_product(moment%tonality_y(1:n/2), moment%k2l(1:n/2)) nonlin%tonality_x=moment%mtonality_x * 390.1 nonlin%tonality_y=moment%mtonality_y * 390.1 m=lat%n_ele_track mp=lat%n_ele_track moment%mdbeta_x(1:m) = matmul(moment%dbeta_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_y(1:m) = matmul(moment%dbeta_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dpretz_x(1:m) = matmul(moment%dbeta_dpretz_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dpretz_y(1:m) = matmul(moment%dbeta_dpretz_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dcos_x(1:m) = matmul(moment%dbeta_dcos_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dcos_y(1:m) = matmul(moment%dbeta_dcos_y(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dsin_x(1:m) = matmul(moment%dbeta_dsin_x(1:m,1:n), moment%k2l(1:n)) moment%mdbeta_dsin_y(1:m) = matmul(moment%dbeta_dsin_y(1:m,1:n), moment%k2l(1:n)) moment%ma_real(1:m) = matmul(moment%a_real(1:m,1:n), moment%k2l(1:n)) moment%ma_image(1:m) = matmul(moment%a_image(1:m,1:n), moment%k2l(1:n)) moment%mb_real(1:m) = matmul(moment%b_real(1:m,1:n), moment%k2l(1:n)) moment%mb_image(1:m) = matmul(moment%b_image(1:m,1:n), moment%k2l(1:n)) if(debug)then open(unit=21, file = 'pretz_energy_moments.dat') write(21, 7)co(0)%vec 7 format('!',6f12.6) write(21, 3) open(unit=24, file = 'dbeta_dcos_moments.dat') write(24, 7)co(0)%vec write(24, 33) open(unit=25, file = 'coupling_moments.dat') write(25, 7)co(0)%vec write(25, 34) open(unit=26) ! call beta_v_amplitude (lat, nonlin) open(unit=27, file = 'dbeta_dcos_calc.dat') write(27, 7)co(0)%vec write(27, '(7a)')'!',' ELE# ',' ELEMENT NAME ',' dbeta_x/dcos ', & ' dbeta_y/dcos ',' dbeta_x/dsin ',' dbeta_y/dsin ' open (unit=23, file = 'pretz_energy_calc.dat' ) write(23, 7)co(0)%vec write(23, 8) 3 format('!',' ELE# ',' ELEMENT NAME ',' dbeta_x/ddelta ', & ' dbeta_y/ddelta ',' dbeta_x/dpretz ',' dbeta_y/dpretz') 33 format('!',' ELE# ',' ELEMENT NAME ',' dbeta_dcos_x ', & ' dbeta_dcos_y ',' dbeta_dsin_x ',' dbeta_dsin_y') 34 format('!',' ELE# ',' ELEMENT NAME ',' a_real ', & ' a_imaginary ',' b_real ',' b_imaginary ') 8 format('!',' ELE# ',' ELEMENT NAME ',' DBETA_X/DDELTA ', & ' DBETA_Y/DDELTA ',' DBETA_X/DPRETZ ',' DBETA_Y/DPRETZ') do k = 1, lat%n_ele_track dbx_pretz = (nonlin%non_ele(k)%a%beta - lat%ele(k)%a%beta)/ & lat%ele(k)%a%beta dby_pretz = (nonlin%non_ele(k)%b%beta - lat%ele(k)%b%beta)/ & lat%ele(k)%b%beta write(21, 1)k, lat%ele(k)%name, moment%mdbeta_x(k)+ & moment%dbeta_x_0(k), & moment%mdbeta_y(k)+moment%dbeta_y_0(k), & moment%mdbeta_dpretz_x(k), & moment%mdbeta_dpretz_y(k) write(24, 1)k, lat%ele(k)%name, moment%mdbeta_dcos_x(k), & moment%mdbeta_dcos_y(k), & moment%mdbeta_dsin_x(k), & moment%mdbeta_dsin_y(k) write(25, 1)k, lat%ele(k)%name, moment%ma_real(k), & moment%ma_image(k), & moment%mb_real(k), & moment%mb_image(k) write(23, 1)k, lat%ele(k)%name, & nonlin%non_ele(k)%a%dbeta, & nonlin%non_ele(k)%b%dbeta, & dbx_pretz, dby_pretz write(27, 1)k, lat%ele(k)%name, & nonlin%non_ele(k)%a%dbeta_dcos, & nonlin%non_ele(k)%b%dbeta_dcos, & nonlin%non_ele(k)%a%dbeta_dsin, & nonlin%non_ele(k)%b%dbeta_dsin 1 format(1x,i4,1x,a16, 4(f12.5,4x)) write(26,'(a, 4f12.4)')lat%ele(k)%name, moment%dbeta_x_0(k), & nonlin%non_ele(k)%a%dbeta, moment%dbeta_y_0(k), & nonlin%non_ele(k)%b%dbeta ! endif end do close(unit=21) close(unit=23) close(unit=24) close(unit=25) close(unit=26) close(unit=27) write(22, 7)co(0)%vec write(22, 2) moment%mchrom_x, moment%mchrom_y 2 format(' chrom_x = ',f8.3,5x,'chrom_y = ',f8.3) write(22,4) moment%mtonality_x*390.1, moment%mtonality_y*390.1 4 format(' tonality_x(kHz) = ', f8.3,5x,'tonality_y(kHz) = ',f8.3) write(22,9) moment%mtonality_x_west*390.1, & moment%mtonality_y_west * 390.1 9 format(' tonality_x_west(kHz) = ', f8.3,5x,'tonality_y_west(kHz) = ',f8.3) write(22,5)sum1,sum2,sum3,sum4, sum5 5 format(' M1 = ',e12.4,3x,'M2 = ',e12.4,3x,'M3 = ',e12.4,3x,'M4 = ',e12.4, & 3x,'M5 = ',e12.4) write(22,6)sum6,sum7,sum8,sum9, sum10 6 format(' M6 = ',e12.4,3x,'M7 = ',e12.4,3x,'M8 = ',e12.4,3x,'M9 = ',e12.4, & 3x,'M10 = ',e12.4) endif !(debug) end subroutine write_moments_debug !----------------------------------------------------------------------- !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine matmulvec(xxx,yyy,m,n,mp,np, zzz) implicit none integer m,mp,n,np real(rp) zzz(0:mp),xxx(0:mp,np),yyy(np) integer i,j do i=1,m zzz(i) = 0. do j = 1,n zzz(i) = zzz(i) + xxx(i,j)*yyy(j) end do end do return end subroutine matmulvec end module nonlin_mod