module opt_mod use bmad use bmadz_interface use bmadz_mod use bmadz_utils use bunchcross_mod use zquad_lens_mod use pretz_mod use constraints_mod use nonlin_mod use twiss_max_mod use cbar_mod use calc_res_mod use mode3_mod use da_interface use da_mod contains !---------------------- subroutine bmad_calcg(lat, con, pc, quad, co, global, & twiss_max, pretz, cesr, nonlin, indep_var, G, moment, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (constraint_struct) con type (lat_struct) lat type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment real(rp), dimension(:), intent(out) :: G(indep_var_maxx) real(rp) fp, f0, xx real(rp), pointer :: a_ptr integer n, way, IGRAD data way/1/, IGRAD/18/ call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) f0 = con%figure_of_merit way = -way do n = 1, indep_var%n_var a_ptr => var_attribute(lat, indep_var%v(n)) xx = a_ptr a_ptr = xx + 0.5**IGRAD*way call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) fp = con%figure_of_merit a_ptr = xx G(n) = way*(fp-f0)*2.**IGRAD end do return end subroutine bmad_calcg ! !........................................................................ ! subroutine bmad_calc_dyda(lat, con, pc, quad, co, global, & twiss_max, pretz, cesr, nonlin, indep_var,dyda, moment, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (constraint_struct) con type (lat_struct) lat type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment real(rp), dimension(:,:), intent(out) :: dyda real(rp) fp, f0, xx real(rp), pointer :: a_ptr real(rp) dx real(rp), dimension(:), allocatable, save :: f_vec, G integer n, way, IGRAD integer ma data way/1/, IGRAD/18/ ma = con%n_constraint allocate(f_vec(ma), G(indep_var%n_var)) call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) f0 = con%figure_of_merit f_vec(1:ma) = con%c(1:ma)%contribution way = -way dx = 0.5**IGRAD*way do n = 1, indep_var%n_var a_ptr => var_attribute(lat, indep_var%v(n)) xx = a_ptr a_ptr = xx + dx call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) fp = con%figure_of_merit a_ptr = xx G(n) = way*(fp-f0)*2.**IGRAD dyda(1:ma,n) = (con%c(1:ma)%contribution - f_vec(1:ma))/dx end do deallocate(f_vec, G) return end subroutine bmad_calc_dyda !........................................................................ !........................................................................ subroutine bmad_minop(lat, indep_var, cesr, con, pc, quad, & co, global, twiss_max, pretz, nonlin, moment, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (indep_var_struct) comp_indep_var type (an_indep_var) thevar type (constraint_struct) con type (lat_struct) lat type (lat_struct), save :: lat_1 type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment real(rp) fig_start, fig_end real(rp) step, acc real*8 x8(indep_var_maxx) real *8 comp_vec(indep_var_maxx+3) character*60 new_bmad_file integer maxfun, idisp integer i logical init data init/.true./ ! maxfun = con%n_cycles idisp = 10 step = 0.5**13 acc = 2**(-9.5) call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) fig_start = con%figure_of_merit call showme_topten(con, lat, 0) lat_1 = lat call minop_b (lat_1, indep_var, cesr, con, pc, quad, co, & global, twiss_max, pretz, nonlin, maxfun, idisp, step, & acc, init, moment, sext_var) call calcf(lat_1, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat_1, cesr, con, pc, quad, co, global, twiss_max, pretz, nonlin, moment, indep_var, sext_var) fig_end = con%figure_of_merit if(fig_end <= fig_start) then lat = lat_1 call showme_topten(con, lat, 0) if(con%plot) call plotdo_bmadz('/xs', lat, co, pc) if(con%new_bmad_file<=' ') then new_bmad_file = 'new_bmad_input.' else new_bmad_file = con%new_bmad_file endif do i = 1, indep_var%n_var x8(i) = var_attribute(lat, indep_var%v(i)) end do if(con%sol_com >= 0 .and. con%sol_com <= 2) then comp_vec(1:indep_var%n_var)= x8(1:indep_var%n_var) call compensation_indep_var(con%sol_com, con%compensating_ele, lat, indep_var, comp_indep_var, comp_vec) call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, comp_vec(1:comp_indep_var%n_var)) else call opt_out(lat%input_file_name,new_bmad_file,indep_var, x8) endif ! call opt_out(lat%lattice, con%new_bmad_file, indep_var, x8) if (con%nonlinearity) call write_moments(lat, co, nonlin, moment) print *, 'fig_start = ', fig_start print *, 'fig_end = ', fig_end endif return end subroutine bmad_minop !........................................................................ !........................................................................ subroutine calcf(lat_1, con, pc, quad, co1, global, twiss_max, pretz, nonlin, indep_var, err_flag) implicit none type (lat_struct), target :: lat_1 type (lat_struct), save :: lat_2 type (coord_struct) orb0 type (coord_struct), save :: closed_start type (coord_struct), allocatable :: co1(:) type (coord_array_struct), save :: co2_arr(0:0) type (ele_struct), pointer :: ele0, ele2 type (ele_struct) ele1 type (pc_struct) pc type (zquad_struct) quad type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (indep_var_struct), optional :: indep_var type (constraint_struct) con type (nonlin_ele_struct) nonlin real(rp) horizontal_emit real(rp) sigma_x, sigma_y real(rp) unstable_factor integer i, status integer j integer ix character*10 ele_names(10) logical err logical :: debug_type = .false. logical :: use_taylor = .false. logical, save :: done = .false. logical, optional :: err_flag ! if (present(err_flag)) err_flag = .false. call reallocate_coord( co2_arr(0)%orbit, lat_1%n_ele_max ) global_com%exit_on_error = .false. do j = 1, 6 do i=0,lat_1%n_ele_track co2_arr(0)%orbit(i)%vec(j) = 0 end do enddo if( present (indep_var)) then do i = 1, indep_var%n_var ix = indep_var%v(i)%rindex call lat_make_mat6(lat_1,ix,co2_arr(0)%orbit, err_flag = err) end do else call lat_make_mat6(lat_1, -1, co2_arr(0)%orbit, err_flag = err) endif if (present(err_flag) .and. err) err_flag = .true. if (con%sol_com >= 0) call solenoid_compensation(lat_1, con%sol_com, con%compensating_ele ) if(.not. done)then closed_start%vec = 0 ele_names(1) = 'V_SEP*' ele_names(2) = 'H_SEP*' ele_names(3) = 'IP_L3*' call name_to_list (lat_1, ele_names(1:3)) done=.true. endif call make_hybrid_lat (lat_1, lat_2, use_taylor, co2_arr) call twiss_at_start (lat_2, status = status) if (present(err_flag) .and. status /= ok$) err_flag = .true. if (lat_2%param%stable .and. con%circular_machine) then lat_1%ele(0) = lat_2%ele(0) lat_1%a%tune = lat_2%a%tune lat_1%b%tune = lat_2%b%tune co1(0) = closed_start ! init guess call closed_orbit_calc( lat_1, co1, 4, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. closed_start = co1(0) ! save for init guess next time around call lat_make_mat6(lat_1, -1, co1, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. call twiss_at_start(lat_1, status = status) if (present(err_flag) .and. status /= ok$) err_flag = .true. else ! lat_1%ele(0) = lat_1%ele_init lat_1%ele(0)%mat6 = lat_2%ele(0)%mat6 ! co1(0) = co1(lat_1%n_ele_max) call track_all(lat_1, co1, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. call lat_make_mat6(lat_1, -1, co1, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. endif lat_1%param%stable = lat_2%param%stable lat_1%param%unstable_factor = lat_2%param%unstable_factor unstable_factor = lat_1%param%unstable_factor call twiss_propagate_all (lat_1, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. if (con%emitt_calc) call emittance_calc(lat_1, co1, con, global) if(global%wig%x_emit <0 )then global%wig%x_emit = abs(global%wig%x_emit) print *,' CALCF: antidamping, negative emittance' endif if(con%pretzel) then call twiss_to_pc (lat_1, global, pc) call track_to_pc (co1, lat_1%param, pc, lat_1) endif twiss_max%arc%a%beta=0. twiss_max%arc%b%beta=0. twiss_max%arc%a%eta=0. twiss_max%arc%b%eta=0. twiss_max%arc%a%sigma=0. twiss_max%arc%b%sigma=0. twiss_max%ir%a%beta=0. twiss_max%ir%b%beta=0. twiss_max%ir%a%eta=0. twiss_max%ir%b%eta=0. twiss_max%ir%a%sigma=0. twiss_max%ir%b%sigma=0. do j = 1, quad%n !get average beta and maximum beta in quads i=quad%lens(j)%ix call quad_beta_ave (lat_1%ele(i), quad%lens(j)%a%beta_ave, quad%lens(j)%b%beta_ave) ! propogate twiss to center of quad ele1 = lat_1%ele(i) ele1%value(k1$) = lat_1%ele(i)%value(k1$) ele1%value(l$) = 0.5 * lat_1%ele(i)%value(l$) if(ele1%key /= wiggler$)then call make_mat6(ele1, lat_1%param) call twiss_propagate1 ( lat_1%ele(i-1), ele1) endif ele0 => lat_1%ele(i-1) ! Twiss at the start of the quad ele2 => lat_1%ele(i) ! Twiss at the end of the quad quad%lens(j)%a%beta_max = max(ele0%a%beta, ele1%a%beta, ele2%a%beta) quad%lens(j)%b%beta_max = max(ele0%b%beta, ele1%b%beta, ele2%b%beta) ! quad%lens(j)%a%mobius_beta_max = max(ele0%a%mobius_beta, & ! ele1%a%mobius_beta, ele2%a%mobius_beta) ! quad%lens(j)%b%mobius_beta_max = max(ele0%b%mobius_beta, & ! ele1%b%mobius_beta, ele2%b%mobius_beta) ! quad%lens(j)%a%mobius_eta_max = max(ele0%a%mobius_eta, & ! ele1%a%mobius_eta, ele2%a%mobius_eta) ! quad%lens(j)%b%mobius_eta_max = max(ele0%b%mobius_eta, & ! ele1%b%mobius_eta, ele2%b%mobius_eta) quad%lens(j)%a%eta_max = max(ele0%a%eta, ele1%a%eta, ele2%a%eta) quad%lens(j)%b%eta_max = max(ele0%b%eta, ele1%b%eta, ele2%b%eta) ! get maximum beta and eta and sigma and location horizontal_emit = global%wig%x_emit if(con%fully_coupled)horizontal_emit = 0.5*global%wig%x_emit if (lat_1%ele(i)%type == 'ARC') then ! maximum beta and eta in arcs if(quad%lens(j)%a%beta_max>twiss_max%arc%a%beta)then twiss_max%arc%a%beta=quad%lens(j)%a%beta_max twiss_max%arc%a%ix_beta=quad%lens(j)%ix endif if(quad%lens(j)%b%beta_max>twiss_max%arc%b%beta)then twiss_max%arc%b%beta=quad%lens(j)%b%beta_max twiss_max%arc%b%ix_beta=quad%lens(j)%ix endif if(quad%lens(j)%a%eta_max>twiss_max%arc%a%eta)then twiss_max%arc%a%eta=quad%lens(j)%a%eta_max twiss_max%arc%a%ix_eta=quad%lens(j)%ix endif if(quad%lens(j)%b%eta_max>twiss_max%arc%b%eta)then twiss_max%arc%b%eta=quad%lens(j)%b%eta_max twiss_max%arc%b%ix_eta=quad%lens(j)%ix endif sigma_x=sqrt(max(0.0, quad%lens(j)%a%beta_max * horizontal_emit + & (global%wig%sige_e * quad%lens(j)%a%eta_max)**2)) lat_1%ele(i)%a%sigma=sigma_x twiss_max%arc%a%sigma=max(twiss_max%arc%a%sigma ,sigma_x) sigma_y=sqrt(max(0.0, quad%lens(j)%b%beta_max * global%wig%x_emit *.5)) lat_1%ele(i)%b%sigma=sigma_y twiss_max%arc%b%sigma=max(twiss_max%arc%b%sigma , sigma_y) else !maximum beta and eta in IR if(quad%lens(j)%a%beta_max>twiss_max%ir%a%beta)then twiss_max%ir%a%beta=quad%lens(j)%a%beta_max twiss_max%ir%a%ix_beta=quad%lens(j)%ix endif if(quad%lens(j)%b%beta_max>twiss_max%ir%b%beta)then twiss_max%ir%b%beta=quad%lens(j)%b%beta_max twiss_max%ir%b%ix_beta=quad%lens(j)%ix endif if(quad%lens(j)%a%eta_max>twiss_max%ir%a%eta)then twiss_max%ir%a%eta=quad%lens(j)%a%eta_max twiss_max%ir%a%ix_eta=quad%lens(j)%ix endif if(quad%lens(j)%b%eta_max>twiss_max%ir%b%eta)then twiss_max%ir%b%eta=quad%lens(j)%b%eta_max twiss_max%ir%b%ix_eta=quad%lens(j)%ix endif sigma_x=sqrt(max(0.0, quad%lens(j)%a%beta_max * horizontal_emit + & (global%wig%sige_e * quad%lens(j)%a%eta_max)**2)) lat_1%ele(i)%a%sigma= sigma_x twiss_max%ir%a%sigma=max(twiss_max%ir%a%sigma , sigma_x) sigma_y=sqrt(max(0.0, quad%lens(j)%b%beta_max * global%wig%x_emit *.5)) lat_1%ele(i)%b%sigma= sigma_y twiss_max%ir%b%sigma=max(twiss_max%ir%b%sigma, sigma_y) endif end do lat_1%param%unstable_factor = unstable_factor if (any(con%c(1:con%n_constraint)%variable == chrom$)) then call chromaticity_calc(lat_1, co1, global, err_flag = err) if (present(err_flag) .and. err) err_flag = .true. endif if (con%pretzel) then call pretzel(lat_1, con, global, pc, co1, quad, pretz) if(con%circular_machine .and. & (any(con%c(1:con%n_constraint)%variable == curly_d$))) & call curly_d(lat_1, co1, quad, global%wig%synch_int(2), pretz%curly_d) endif if (con%nonlinearity) then call nonlinearity(lat_1, co1, nonlin, con, err) if (present(err_flag) .and. err) err_flag = .true. endif ! if (any(con%c(1:con%n_constraint)%variable == d_beta_beam$).or. & ! any(con%c(1:con%n_constraint)%variable == d_alpha_beam$) .or. & ! any(con%c(1:con%n_constraint)%variable == d_xemit$) .or. & ! any(con%c(1:con%n_constraint)%variable == d_curly_d$) .or. & ! any(con%c(1:con%n_constraint)%variable == tonality$)) & ! call electron_positron (lat_1, lat_pos, lat_ele, con, global) ! call calc_z_tune(lat_1) end subroutine calcf !........................................................................ ! map constraint location to lat element ! ! elseif (variable == particle_tracking$) then ! call particle_tracking (lat, con, value) ! figure figure of merit !--------------------------------------------------------------------------- ! get ordered list of contrib(i) and map i(j) where i(1) is max contrib(i) etc. ! put everything together !........................................................................ subroutine fom (lat, cesr, CON, pc, quad, co, global, twiss_max, & pretz, nonlin, moment, indep_var, sext_var) implicit none type (lat_struct), target :: lat type (ele_struct), pointer :: slave, this_ele type (lat_struct), save :: lat_pretz_off, lat_pos, lat_ele type (coord_struct), allocatable :: co(:) type (coord_struct), allocatable, save :: vc_diff(:), vc7_diff(:) type (coord_struct), allocatable, save :: prz1_diff(:), prz13_diff(:) type (coord_struct), allocatable, save :: diff(:) type (coord_struct), allocatable, save :: co_positron(:), co_electron(:), co_pretz_off(:) type (db_struct) cesr type (zquad_struct) quad type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (constraint_struct) con type (indep_var_struct) indep_var, sext_var type (nonlin_ele_struct) nonlin type (moment_struct) moment type (ele_struct) ele_vert, ele_horz, ele_det, ele type cbar_struct real(rp) cbars(2,2) end type cbar_struct type (cbar_struct), allocatable, save :: cbar_ele(:) type (osc_param_struct) osc_param integer i,j, jsep(4), k, n integer l integer variable, plane, i1, i2, constraint integer location, ix, jshx, jsh_count integer nbunches integer n_turn real(rp) value, volt, v, vmax, x0, x_pos, y_pos, x_vel, y_vel real(rp) contrib, delta, weight, figure, target real(rp) I_x, I_y, I_wt, I_max real(rp) aperture, val_max, k2l real(rp) :: fact = 1 real(rp) v1, target_value, cbar(2,2) real(rp) max_tune, min_tune real(rp) s_vert, s_horz, s_det real(rp) d_amp_x, d_amp_y real(rp) axx, axy, ayy real(rp), allocatable :: etax2(:), betax1(:), betay1(:) real(rp) res_cos, res_sin, res_amp, max_x real(rp) res_cos_ele, res_sin_ele, res_amp_ele real(rp) slopes(4) real(rp) xix, xiy, dnux_dp2, dnuy_dp2 real(rp) betax1_rms, betay1_rms, etax2_rms real(rp) coord(6), jvar(2), ptc_fps_param(3) complex*8 h_nonlin(10), h_nonlin2(11) logical sextupole_constraint, error/.false./ logical verbose/.false./ call reallocate_coord( vc_diff, lat%n_ele_max ) call reallocate_coord( vc7_diff, lat%n_ele_max ) call reallocate_coord( prz1_diff, lat%n_ele_max ) call reallocate_coord( prz13_diff, lat%n_ele_max ) call reallocate_coord( diff, lat%n_ele_max ) call reallocate_coord( co_positron, lat%n_ele_max ) call reallocate_coord( co_electron, lat%n_ele_max ) call reallocate_coord( co_pretz_off, lat%n_ele_max ) allocate(cbar_ele(0:lat%n_ele_max)) ! if(any(con%c(1:con%n_constraint)%variable == dx_vcros5$))call dx_dvcros5_calc (lat, vc_diff) if(any(con%c(1:con%n_constraint)%variable == dx_vcros7$))call dx_dvcros7_calc (lat, vc7_diff) if(any(con%c(1:con%n_constraint)%variable == dx_pretz1$))call dx_dpretz1_calc (lat, prz1_diff) if(any(con%c(1:con%n_constraint)%variable == dx_pretz13$))call dx_dpretz13_calc (lat, prz13_diff) ! write(11,*)' FOM' ! do i=1,lat%n_ele_max ! if(lat%ele(i)%s > 20. .and. lat%ele(i)%s < 748.)cycle ! if(lat%ele(i)%key /= quadrupole$)cycle ! write(11, '(i7,1x,a16,2e12.4)')i,lat%ele(i)%name, lat%ele(i)%value(k1$), & ! lat%ele(i)%value(tilt$) ! end do if(any(con%c(1:con%n_constraint)%variable == dcbar12$))call cbar_v_e(lat, con%compensating_ele, slopes) figure = 0.0 I_x = 99.0 I_y = 99.0 I_wt = 99.0 coord(:) = 0.0_rp coord(1) = 22.0E-3 coord(3) = 2.0E-3 jvar(:) = 0.0_rp ptc_fps_param(:) = 0.0_rp if (any(con%c(1:con%n_constraint)%variable == d_beta_beam$).or. & any(con%c(1:con%n_constraint)%variable == d_alpha_beam$) .or. & any(con%c(1:con%n_constraint)%variable == d_xemit$) .or. & any(con%c(1:con%n_constraint)%variable == d_curly_d$) .or. & any(con%c(1:con%n_constraint)%variable == tonality$) .or. & any(con%c(1:con%n_constraint)%variable == elec_emitt$) .or. & any(con%c(1:con%n_constraint)%variable == beta_electrons$) .or. & any(con%c(1:con%n_constraint)%variable == eta_electrons$) .or. & any(con%c(1:con%n_constraint)%variable == etap_electrons$) .or. & any(con%c(1:con%n_constraint)%variable == qs_2qx_ele$) .or. & any(con%c(1:con%n_constraint)%variable == diff_displace$) .or. & any(con%c(1:con%n_constraint)%variable == elec_co_angle$) .or. & any(con%c(1:con%n_constraint)%variable == elec_displace$ )) then co_positron(0)%vec = co(0)%vec co_electron(0)%vec = co(0)%vec co_electron(0)%vec(2)=-co(0)%vec(2) co_electron(0)%vec(4)=-co(0)%vec(4) call electron_positron (lat, lat_pos, lat_ele, co_positron, co_electron, con, global) do i = 0,lat_ele%n_ele_track nonlin%nonele_ele(i)%tx%beta = lat_ele%ele(i)%a%beta nonlin%nonele_ele(i)%ty%beta = lat_ele%ele(i)%b%beta nonlin%nonele_ele(i)%tx%eta = lat_ele%ele(i)%a%eta nonlin%nonele_ele(i)%co_off(2)%vec(1:6) = co_electron(i)%vec(1:6) end do else lat_ele = lat lat_pos = lat co_positron = co endif if(any(con%c(1:con%n_constraint)%variable == cbar11$) .or. & any(con%c(1:con%n_constraint)%variable == cbar12$) .or. & any(con%c(1:con%n_constraint)%variable == cbar21$) .or. & any(con%c(1:con%n_constraint)%variable == cbar22$))then do j=1, lat%n_ele_track cbar_ele(j)%cbars(1:2,1:2) = 0. call c_to_cbar(lat%ele(j), cbar_ele(j)%cbars) end do endif if(any(con%c(1:con%n_constraint)%variable == osc_emit_max$) .or. & any(con%c(1:con%n_constraint)%variable == osc_dp_max$) .or. & any(con%c(1:con%n_constraint)%variable == osc_delta_s$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m56$) .or. & any(con%c(1:con%n_constraint)%variable == osc_tilde_m56$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m51$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m52$)) osc_param%linear=.true. if(any(con%c(1:con%n_constraint)%variable == osc_m511$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m512$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m516$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m522$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m526$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m533$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m534$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m536$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m544$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m546$) .or. & any(con%c(1:con%n_constraint)%variable == osc_sec_ord_x$) .or. & any(con%c(1:con%n_constraint)%variable == osc_sec_ord_y$) .or. & any(con%c(1:con%n_constraint)%variable == osc_m566$)) osc_param%t2mat = .true. if(any(con%c(1:con%n_constraint)%variable == osc_int_theta_x2$) .or. & any(con%c(1:con%n_constraint)%variable == osc_int_theta_y2$)) osc_param%theta2 = .true. if(osc_param%linear .or. osc_param%t2mat .or. osc_param%theta2)then do i=1,con%n_constraint if(con%c(i)%variable == osc_wavelength$)then osc_param%wavelength = con%c(i)%target_value i1 = con%c(i)%where1_ix i2 = con%c(i)%where2_ix endif end do if(osc_param%wavelength == 0.) then print *,' No wavelength specified for OSC_PARAMETERS. Set wavelength using OSC_WAVELENGTH in constraint file ' stop endif call osc_parameters(lat, i1,i2, co, osc_param, verbose) endif if(any(con%c(1:con%n_constraint)%variable == V15$) .or. & any(con%c(1:con%n_constraint)%variable == V16$) .or. & any(con%c(1:con%n_constraint)%variable == V25$) .or. & any(con%c(1:con%n_constraint)%variable == V35$) .or. & any(con%c(1:con%n_constraint)%variable == V26$))then call set_on_off( rfcavity$, lat, on$) call twiss3_at_start(lat,error) if(error)then print *,' FOM: twiss3_at_start returns error' stop endif call twiss3_propagate_all(lat) call set_on_off( rfcavity$, lat, off$) endif lat_pretz_off = lat forall(j=1:lat%n_ele_track)co_pretz_off(j)%vec(1:6)=0 call set_on_off( elseparator$, lat_pretz_off, off$) lat_pretz_off%param%n_part = 0 if(con%circular_machine)then call twiss_at_start(lat_pretz_off) call closed_orbit_calc(lat_pretz_off, co_pretz_off, 4) endif call lat_make_mat6(lat_pretz_off, -1, co_pretz_off) if(con%circular_machine)then call twiss_at_start(lat_pretz_off) else lat_pretz_off%ele(0) = lat%ele_init endif call twiss_propagate_all(lat_pretz_off) call set_on_off( elseparator$, lat_pretz_off, on$ ) if (con%calculate_moments .or. & any(con%c(1:con%n_constraint)%variable == alpha_xx$) .or. & any(con%c(1:con%n_constraint)%variable == alpha_xy$) .or. & any(con%c(1:con%n_constraint)%variable == alpha_yy$)) then call sextupole_moments (lat_pretz_off, co_positron, quad, con, moment) ! call calc_h(lat_pretz_off, h_nonlin ) call calc_h1_h2(lat_pretz_off, h_nonlin, h_nonlin2, 2) axx = real(h_nonlin2(1)) axy = real(h_nonlin2(2)) ayy = real(h_nonlin2(3)) endif if(any(con%c(1:con%n_constraint)%variable == dnux_dp2$) .or. & any(con%c(1:con%n_constraint)%variable == dnuy_dp2$))then call calc_xi(lat_pretz_off, xix, xiy) call calc_xi2(lat_pretz_off,dnux_dp2, dnuy_dp2) endif if (con%minimize_moments)then call minimize_moments (lat_pretz_off, moment, sext_var, con) do i = 1 , lat_pos%n_ele_track if(lat_pos%ele(i)%key == sextupole$ .and.lat_pos%ele(i)%value(tilt$) == 0) & lat_pos%ele(i)%value(k2$) = lat_pretz_off%ele(i)%value(k2$) lat_ele%ele(i)%value(k2$) = lat_pos%ele(i)%value(k2$) lat%ele(i)%value(k2$) = lat_pos%ele(i)%value(k2$) end do endif if(con%nonlinearity) call write_moments(lat_pretz_off, co_positron, nonlin, moment) ! if (any(con%c(1:con%n_constraint)%variable == alpha_xx$) .or. & ! any(con%c(1:con%n_constraint)%variable == alpha_xy$) .or. & ! any(con%c(1:con%n_constraint)%variable == alpha_yy$)) call calc_dnu_dJ(lat_pos, axx,axy,ayy) ! added 2015.03.09 - if (any(con%c(1:con%n_constraint)%variable == eta_x2$) .or. & any(con%c(1:con%n_constraint)%variable == beta_x1$) .or. & any(con%c(1:con%n_constraint)%variable == beta_y1$)) call calc_dbeta_ddelta(lat_pos, etax2, betax1, betay1) if (any(con%c(1:con%n_constraint)%variable == Qs_2Qx$) .or. & any(con%c(1:con%n_constraint)%variable == Qs_2Qx_track$)) & call sext_energy(lat_pos, co_positron(0), con, res_cos, res_sin, res_amp, max_x) if (any(con%c(1:con%n_constraint)%variable == Qs_2Qx_ele$) .or. & any(con%c(1:con%n_constraint)%variable == Qs_2Qx_cos_ele$)) & call sext_energy(lat_ele, co_electron(0), con, res_cos_ele, res_sin_ele, res_amp_ele, max_x) ! jx jy variant from tracking if (any(con%c(1:con%n_constraint)%variable == jx_var$) .or. & any(con%c(1:con%n_constraint)%variable == jy_var$)) then if (con%track%n_turn .gt. 0 .and. con%track%n_turn .lt. 1E4) then call calc_jwvar(lat_pretz_off, co_pretz_off(0), con%track%n_turn, con%track%start%vec, jvar) else call calc_jwvar(lat_pretz_off, co_pretz_off(0), 100, coord, jvar) ! if read_tracking = F, no input vec end if end if !ptc fixed points param if (any(con%c(1:con%n_constraint)%variable == ptc_axx$) .or. & any(con%c(1:con%n_constraint)%variable == ptc_g_real$) .or. & any(con%c(1:con%n_constraint)%variable == ptc_g_imag$)) then call ptc_calc_fps_param (lat_pretz_off, 5, ptc_fps_param) end if ! map constraint location to lat element main_loop: do i = 1, con%n_constraint sextupole_constraint = .false. contrib = 0.0 location=-1. ! no location identified variable = con%c(i)%variable plane = con%c(i)%plane constraint = con%c(i)%constraint weight = con%c(i)%weight target_value = con%c(i)%target_value i1 = con%c(i)%where1_ix i2 = con%c(i)%where2_ix ! if (variable == beta$) then if (constraint == maxz$) then call max_or_min(CON, LAT, QUAD, i, value, location) else if (plane == x_plane$) then value = lat%ele(i1)%a%beta elseif (plane == y_plane$) then value = lat%ele(i1)%b%beta else call plane_error (con%c(i)) endif endif elseif (variable == eta$) then if (constraint == maxz$) then call max_or_min(CON, LAT, QUAD, i, value, location) else if (plane == x_plane$) then value = lat%ele(i1)%a%eta elseif (plane == y_plane$) then value = lat%ele(i1)%b%eta else call plane_error (con%c(i)) endif endif elseif (variable == beta_electrons$) then if (constraint == maxz$) then call max_or_min(CON, LAT_ele, QUAD, i, value, location) else if (plane == x_plane$) then value = lat_ele%ele(i1)%a%beta elseif (plane == y_plane$) then value = lat_ele%ele(i1)%b%beta else call plane_error (con%c(i)) endif endif elseif (variable == eta_electrons$) then if (constraint == maxz$) then call max_or_min(CON, LAT_ele, QUAD, i, value, location) else if (plane == x_plane$) then value = lat_ele%ele(i1)%a%eta elseif (plane == y_plane$) then value = lat_ele%ele(i1)%b%eta else call plane_error (con%c(i)) endif endif elseif (variable == etap_electrons$) then if (constraint == maxz$) then call max_or_min(CON, LAT_ele, QUAD, i, value, location) else if (plane == x_plane$) then value = lat_ele%ele(i1)%a%etap elseif (plane == y_plane$) then value = lat_ele%ele(i1)%b%etap else call plane_error (con%c(i)) endif endif elseif (variable == curly_H$) then if (constraint == maxz$) then call max_or_min(CON, LAT, QUAD, i, value, location) else ele = lat%ele(i1) if (plane == x_plane$) then value = ele%a%beta*ele%a%etap**2 + ele%a%alpha * ele%a%eta * ele%a%etap & + ele%a%gamma * ele%a%eta**2 elseif (plane == y_plane$) then value = ele%b%beta*ele%b%etap**2 + ele%b%alpha * ele%b%eta * ele%b%etap & + ele%b%gamma * ele%b%eta**2 else call plane_error (con%c(i)) endif endif elseif (variable == d_path_length$) then if (plane == x_plane$) then !jsep(1)=cesr%h_sep(h_sep_08e$)%ix_lat !jsep(2)=cesr%h_sep(h_sep_08w$)%ix_lat !jsep(3)=cesr%h_sep(h_sep_45e$)%ix_lat !jsep(4)=cesr%h_sep(h_sep_45w$)%ix_lat !value = 0.0 !do j=1,4 ! this_ele => lat%ele(jsep(j)) ! if (this_ele%lord_status == super_lord$) this_ele => pointer_to_slave(this_ele, 1) ! value = value + this_ele%a%eta*this_ele%value(hkick$) !enddo if(any(con%c(1:con%n_constraint)%variable == d_beta_beam$) .or. & any(con%c(1:con%n_constraint)%variable == d_alpha_beam$) .or. & any(con%c(1:con%n_constraint)%variable == d_eta_beam$) .or. & any(con%c(1:con%n_constraint)%variable == beta_electrons$)) then n=lat%n_ele_track value = co_positron(n)%vec(5) - & co_positron(0)%vec(5) - & co_electron(n)%vec(5) + & co_electron(0)%vec(5) endif else call plane_error (con%c(i)) endif elseif (variable == alpha$) then if (plane == x_plane$) then value = lat%ele(i1)%a%alpha elseif (plane == y_plane$) then value = lat%ele(i1)%b%alpha else call plane_error (con%c(i)) endif elseif (variable == betax_betay$) then value = lat%ele(i1)%a%beta - lat%ele(i1)%b%beta elseif (variable == alphax_alphay$) then value = lat%ele(i1)%a%alpha - lat%ele(i1)%b%alpha elseif (variable == phix_phiy$) then value = lat%ele(i2)%a%phi - lat%ele(i1)%a%phi - & (lat%ele(i2)%b%phi - lat%ele(i1)%b%phi) elseif (variable == etap$) then if (plane == x_plane$) then value = lat%ele(i1)%a%etap elseif (plane == y_plane$) then value = lat%ele(i1)%b%etap else call plane_error (con%c(i)) endif elseif (variable == phase$) then if (plane == x_plane$) then value = lat%ele(i2)%a%phi - lat%ele(i1)%a%phi elseif (plane == y_plane$) then value = lat%ele(i2)%b%phi - lat%ele(i1)%b%phi else call plane_error (con%c(i)) endif value = value / twopi elseif (variable == unstable_ring$) then value = lat%param%unstable_factor if (con%nonlinearity) value= max(nonlin%growth_rate(1), & nonlin%growth_rate(2),nonlin%growth_rate(3)) elseif (variable == q$) then if (plane == x_plane$) then value = lat%a%tune / twopi elseif (plane == y_plane$) then value = lat%b%tune / twopi elseif (plane == z_plane$) then value = lat%z%tune / twopi else call plane_error (con%c(i)) endif if (value < 0.0) value = value + 1.0 elseif (variable == q_ele$) then if (plane == x_plane$) then value = lat_ele%a%tune / twopi elseif (plane == y_plane$) then value = lat_ele%b%tune / twopi elseif (plane == z_plane$) then value = lat_ele%z%tune / twopi else call plane_error (con%c(i)) endif if (value < 0.0) value = value + 1.0 elseif (variable == q_tot$) then if (plane == x_plane$) then value = lat%ele(lat%n_ele_track)%a%phi / twopi elseif (plane == y_plane$) then value = lat%ele(lat%n_ele_track)%b%phi / twopi endif ! elseif (variable == particle_tracking$) then ! call particle_tracking (lat, con, value) elseif (variable == pretz_aperture$) then if (con%c(i)%where1 == 'ARC') then if (plane == x_plane$) then value = pretz%x%arc_ape location = pretz%x%arc_ape_ix elseif (plane == y_plane$) then value = pretz%y%arc_ape location = pretz%y%arc_ape_ix else call plane_error (con%c(i)) endif endif if (con%c(i)%where1 == 'IR') then if (plane == x_plane$) then value = pretz%x%ir_ape location = pretz%x%ir_ape_ix elseif (plane == y_plane$) then value = pretz%y%ir_ape location = pretz%y%ir_ape_ix else call plane_error (con%c(i)) endif endif if (i2 >= i1 .and. i1 /= 0)then value = 0. do j=i1,i2 if (plane == x_plane$)then aperture = abs((co(j)%vec(1)))+ & pretz%n_sigma * lat%ele(j)%a%sigma if (aperture > value)then value=aperture location = j endif endif end do endif elseif (variable == beam_aperture$) then if (con%c(i)%where1 == 'ARC') then if (plane == x_plane$) then value = twiss_max%arc%a%sigma elseif (plane == y_plane$) then value = twiss_max%arc%b%sigma else call plane_error (con%c(i)) endif elseif (con%c(i)%where1 == 'IR') then if (plane == x_plane$) then value = twiss_max%ir%a%sigma elseif (plane == y_plane$) then value = twiss_max%ir%b%sigma else call plane_error (con%c(i)) endif elseif (i1 == i2 .and. i1/=0) then if (plane == x_plane$)then value = lat%ele(i1)%a%sigma elseif (plane == y_plane$)then value = lat%ele(i1)%b%sigma else call plane_error (con%c(i)) endif else print *, 'ERROR IN FOM: WHERE FOR APERATURE CONSTRAINT NOT' print *, ' "ARC" NOR "IR" NOR "ELEMENT NAME"' call err_exit endif elseif (variable == emittance$) then value = global%nowig%x_emit elseif (variable == wig_emittance$) then value = global%wig%x_emit elseif (variable == energy_spread$) then value = global%wig%sige_e elseif (variable == cross_angle$) then value = co(0)%vec(2) elseif (variable == momentum_comp$) then value = global%wig%synch_int(1)/lat%ele(lat%n_ele_track)%s elseif (variable == chrom$) then if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif if (plane == x_plane$) then value = global%chromx if (con%nonlinearity) value = nonlin%chrom_x elseif (plane == y_plane$) then value = global%chromy if (con%nonlinearity) value = nonlin%chrom_y else call plane_error (con%c(i)) endif elseif (variable == dchrom$) then if (.not. con%nonlinearity)weight=0 if (plane == x_plane$) then value = nonlin%dchrom_x elseif (plane == y_plane$) then value = nonlin%dchrom_y else call plane_error (con%c(i)) endif elseif (variable == c_mat_chrom12$) then if (.not. con%nonlinearity)weight=0 value = global%c_mat_chrom(1,2) elseif (variable == nonlin_q$) then if (.not. con%nonlinearity)weight=0 if (plane == x_plane$) then value = nonlin%tune_x elseif (plane == y_plane$) then value = nonlin%tune_y else call plane_error (con%c(i)) endif elseif (variable == tonality$) then if (con%minimize_moments .or. .not. con%nonlinearity) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif if (plane == x_plane$) then ! v1 = lat%a%tune/twopi value = lat_pos%ele(lat%n_ele_track)%a%phi & - lat_ele%ele(lat%n_ele_track)%a%phi elseif (plane == y_plane$) then ! v1 = lat%b%tune / twopi value = lat_pos%ele(lat%n_ele_track)%b%phi & - lat_ele%ele(lat%n_ele_track)%b%phi else call plane_error (con%c(i)) endif ! if (v1 < 0.) v1=1.+v1 ! value = (value-v1)*390.1 value = value / twopi * 390.1 /2. !divide by 2 for delta f from pretz off to on elseif (variable == lrbbi_del_q$) then if (i1 <= i2 .and. i1 > 0.)then value=0. do j=1,pc%total_pc if (pc%cross(j)%ix >= i1 .and. pc%cross(j)%ix <= i2)then if (plane == x_plane$)then if (abs(pc%cross(j)%dnuh) > value)then value = pc%cross(j)%dnuh location = pc%cross(j)%ix if (value > 0.)I_x = target_value / value endif endif if (plane == y_plane$)then if (abs(pc%cross(j)%dnuv) > value)then value = pc%cross(j)%dnuv location = pc%cross(j)%ix if (value > 0) I_y = target_value / value endif endif endif end do endif if (con%c(i)%where1 == 'GLOBAL')then if (plane == x_plane$) then value = pretz%x%dnu location = pretz%x%dnu_ix if (value > 0) I_x = target_value / value elseif (plane == y_plane$) then value = pretz%y%dnu location = pretz%y%dnu_ix if (value > 0) I_y = target_value / value else call plane_error (con%c(i)) endif endif elseif (variable == lrbbi_del_q_x$) then if (con%c(i)%where1 == 'GLOBAL')then if (plane == x_plane$) then value = pretz%x%dnu_x location = pretz%x%dnu_ix_x if (value > 0) I_x = target_value / value elseif (plane == y_plane$) then value = pretz%y%dnu_x location = pretz%y%dnu_ix_x if (value > 0) I_y = target_value / value else call plane_error (con%c(i)) endif endif elseif (variable == lrbbi_del_q_y$) then if (con%c(i)%where1 == 'GLOBAL')then if (plane == x_plane$) then value = pretz%x%dnu_y location = pretz%x%dnu_ix_y elseif (plane == y_plane$) then value = pretz%y%dnu_y location = pretz%y%dnu_ix_y else call plane_error (con%c(i)) endif endif elseif (variable == lrbbi_max_spread$) then if (con%c(i)%where1 == 'GLOBAL')then min_tune=99. max_tune=-99. nbunches = con%n_trains*con%n_cars if (plane == x_plane$) then do j=1,nbunches min_tune = min(pretz%x%dnu_tot(j),min_tune) max_tune = max(pretz%x%dnu_tot(j),max_tune) end do value = max_tune - min_tune elseif (plane == y_plane$) then do j=1,nbunches min_tune = min(pretz%y%dnu_tot(j),min_tune) max_tune = max(pretz%y%dnu_tot(j),max_tune) end do value = max_tune - min_tune else call plane_error (con%c(i)) endif endif elseif (variable == w_t$) then value = pretz%wt if (value > 0) I_wt = target_value / value elseif (variable == b_param_i$) then if (pretz%b_param < 0.00001) then value = -999.0 else value = 1.0 / pretz%b_param endif elseif (variable == pretz_efficiency$) then value=pretz%efficiency elseif (variable == b_efficiency$) then value=pretz%b_efficiency elseif (variable == vert_efficiency$) then value=pretz%vert_efficiency elseif (variable == sigma_sep$) then value = pretz%sigma_sep ! elseif (variable == displace$) then elseif (variable == displacement_$) then if( i2 /= 0 .and. (constraint == maxz$ .or. constraint == minz$ )) then call max_or_min_disp(con, co, i, value, location) else x_pos = co(con%c(i)%where1_ix)%vec(1) y_pos = co(con%c(i)%where1_ix)%vec(3) if (plane == x_plane$) then value = x_pos elseif (plane == y_plane$) then value = y_pos else value = sqrt(x_pos**2 + y_pos**2) endif endif elseif (variable == elec_displace$) then if( i2 /= 0 .and. (constraint == maxz$ .or. constraint == minz$ )) then call max_or_min_disp(con, co_electron, i, value, location) else x_pos = co_electron(con%c(i)%where1_ix)%vec(1) y_pos = co_electron(con%c(i)%where1_ix)%vec(3) if (plane == x_plane$) then value = x_pos elseif (plane == y_plane$) then value = y_pos else value = sqrt(x_pos**2 + y_pos**2) endif endif elseif (variable == diff_displace$) then forall(l=1:lat%n_ele_track) diff(l)%vec(1:6) = co_electron(l)%vec(1:6) - co_positron(l)%vec(1:6) if( i2 /= 0 .and. (constraint == maxz$ .or. constraint == minz$ )) then call max_or_min_disp(con, diff, i, value, location) else x_pos = diff(con%c(i)%where1_ix)%vec(1) y_pos = diff(con%c(i)%where1_ix)%vec(3) if (plane == x_plane$) then value = x_pos elseif (plane == y_plane$) then value = y_pos else value = sqrt(x_pos**2 + y_pos**2) endif endif elseif (variable == co_angle$) then if( i2 /= 0 .and. (constraint == maxz$ .or. constraint == minz$ )) then call max_or_min_disp(con, co, i, value, location) else x_vel = co(con%c(i)%where1_ix)%vec(2) y_vel = co(con%c(i)%where1_ix)%vec(4) if (plane == x_plane$) then value = x_vel elseif (plane == y_plane$) then value = y_vel else value = sqrt(x_vel**2 + y_vel**2) endif endif elseif (variable == elec_co_angle$) then if( i2 /= 0 .and. (constraint == maxz$ .or. constraint == minz$ )) then call max_or_min_disp(con, co_electron, i, value, location) else x_vel = co_electron(con%c(i)%where1_ix)%vec(2) y_vel = co_electron(con%c(i)%where1_ix)%vec(4) if (plane == x_plane$) then value = x_vel elseif (plane == y_plane$) then value = y_vel else value = sqrt(x_vel**2 + y_vel**2) endif endif elseif (variable == curly_d$) then value = pretz%curly_d elseif (variable == c11$) then value = lat%ele(i1)%c_mat(1,1) if(i2 > i1)then value=maxval(abs(lat%ele(i1:i2)%c_mat(1,1))) location = maxloc(abs(lat%ele(i1:i2)%c_mat(1,1)),1) endif elseif (variable == c12$) then value = lat%ele(i1)%c_mat(1,2) if(i2 > i1)then value=maxval(abs(lat%ele(i1:i2)%c_mat(1,2))) location = maxloc(abs(lat%ele(i1:i2)%c_mat(1,2)),1) endif elseif (variable == c21$) then value = lat%ele(i1)%c_mat(2,1) if(i2 > i1)then value=maxval(abs(lat%ele(i1:i2)%c_mat(2,1))) location = maxloc(abs(lat%ele(i1:i2)%c_mat(2,1)),1) endif elseif (variable == c22$) then value = lat%ele(i1)%c_mat(2,2) if(i2 > i1)then value=maxval(abs(lat%ele(i1:i2)%c_mat(2,2))) location = maxloc(abs(lat%ele(i1:i2)%c_mat(2,2)),1) endif elseif (variable == cbar11$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,1) if(i2 > i1)then value=maxval(abs(cbar_ele(i1:i2)%cbars(1,1))) location = maxloc(abs(cbar_ele(i1:i2)%cbars(1,1)),1) endif elseif (variable == cbar12$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,2) if(i2 > i1)then value=maxval(abs(cbar_ele(i1:i2)%cbars(1,2))) location = maxloc(abs(cbar_ele(i1:i2)%cbars(1,2)),1) endif elseif (variable == cbar21$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(2,1) if(i2 > i1)then value=maxval(abs(cbar_ele(i1:i2)%cbars(2,1))) location = maxloc(abs(cbar_ele(i1:i2)%cbars(2,1)),1) endif elseif (variable == cbar22$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(2,2) if(i2 > i1)then value=maxval(abs(cbar_ele(i1:i2)%cbars(2,2))) location = maxloc(abs(cbar_ele(i1:i2)%cbars(2,2)),1) endif elseif (variable == sum_cbar11_cbar22$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,1)+cbar(2,2) elseif (variable == dif_cbar11_cbar22$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,1)-cbar(2,2) elseif (variable == sum_cbar12_cbar21$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,2)+cbar(2,1) elseif (variable == dif_cbar12_cbar21$) then call c_to_cbar (lat%ele(i1), cbar) value = cbar(1,2)-cbar(2,1) elseif (variable == dcbar11$) then value = slopes(1) elseif (variable == dcbar12$) then value = slopes(2) elseif (variable == dcbar21$) then value = slopes(3) elseif (variable == dcbar22$) then value = slopes(4) elseif (variable == V15$) then value = lat%ele(i1)%mode3%V(1,5) elseif (variable == V25$) then value = lat%ele(i1)%mode3%V(2,5) elseif (variable == V35$) then value = lat%ele(i1)%mode3%V(3,5) elseif (variable == V16$) then value = lat%ele(i1)%mode3%V(1,6) elseif (variable == V26$) then value = lat%ele(i1)%mode3%V(2,6) elseif (variable == osc_wavelength$)then value = con%c(i)%target_value ! elseif(variable == osc_emit_max$ .or. variable == osc_dp_max$ .or. variable == osc_delta_s$ & ! .or. variable == osc_m56$ .or. variable == osc_tilde_m56$ .or. variable == osc_m51$ & ! .or. variable == osc_m52$.or. variable == osc_m511$.or. variable == osc_m512$ & ! .or. variable == osc_m516$.or. variable == osc_m522$.or. variable == osc_m526$ & ! .or. variable == osc_m566$ .or. variable == osc_int_theta2$)then ! figure out which osc calculations are necessary elseif(variable == osc_emit_max$)then value = osc_param%osc_emit_max elseif(variable == osc_dp_max$)then value = osc_param%osc_dp_max elseif(variable == osc_m56$)then value = osc_param%Tmat(5,6) elseif(variable == osc_tilde_m56$)then value = osc_param%tilde_m56 elseif(variable == osc_delta_s$)then value = osc_param%delta_s elseif(variable == osc_m51$)then value = osc_param%Tmat(5,1) elseif(variable == osc_m52$)then value = osc_param%Tmat(5,2) elseif(variable == osc_m511$)then value = osc_param%T5mat(1,1) elseif(variable == osc_m512$)then value = osc_param%T5mat(1,2) elseif(variable == osc_m516$)then value = osc_param%T5mat(1,6) elseif(variable == osc_m522$)then value = osc_param%T5mat(2,2) elseif(variable == osc_m526$)then value = osc_param%T5mat(2,6) elseif(variable == osc_m533$)then value = osc_param%T5mat(3,3) elseif(variable == osc_m534$)then value = osc_param%T5mat(3,4) elseif(variable == osc_m536$)then value = osc_param%T5mat(3,6) elseif(variable == osc_m544$)then value = osc_param%T5mat(4,4) elseif(variable == osc_m546$)then value = osc_param%T5mat(4,6) elseif(variable == osc_m566$)then value = osc_param%T5mat(6,6) elseif(variable == osc_int_theta_x2$)then value = osc_param%int_theta_x2 elseif(variable == osc_int_theta_y2$)then value = osc_param%int_theta_y2 elseif(variable == osc_sec_ord_x$)then value = osc_param%second_order_x elseif(variable == osc_sec_ord_y$)then value = osc_param%second_order_y elseif (variable == dx_vcros5$) then if(plane == x_plane$) value = vc_diff(i1)%vec(1) if(plane == y_plane$) value = vc_diff(i1)%vec(3) elseif (variable == dx_vcros7$) then if(plane == x_plane$) value = vc7_diff(i1)%vec(1) if(plane == y_plane$) value = vc7_diff(i1)%vec(3) elseif (variable == dx_pretz1$) then if(plane == x_plane$) value = prz1_diff(i1)%vec(1) if(plane == y_plane$) value = prz1_diff(i1)%vec(3) elseif (variable == dx_pretz13$) then if(plane == x_plane$) value = prz13_diff(i1)%vec(1) if(plane == y_plane$) value = prz13_diff(i1)%vec(3) elseif (variable == d_beta_beam$) then if(plane == x_plane$) value = lat_pos%ele(i1)%a%beta - lat_ele%ele(i1)%a%beta if(plane == y_plane$) value = lat_pos%ele(i1)%b%beta - lat_ele%ele(i1)%b%beta elseif (variable == d_alpha_beam$) then if(plane == x_plane$) value = lat_pos%ele(i1)%a%alpha - lat_ele%ele(i1)%a%alpha if(plane == y_plane$) value = lat_pos%ele(i1)%b%alpha - lat_ele%ele(i1)%b%alpha elseif (variable == d_eta_beam$) then if(plane == x_plane$) value = lat_pos%ele(i1)%a%eta - lat_ele%ele(i1)%a%eta if(plane == y_plane$) value = lat_pos%ele(i1)%b%eta - lat_ele%ele(i1)%b%eta elseif (variable == avg_eta_beam$) then if(plane == x_plane$) value = (lat_pos%ele(i1)%a%eta + lat_ele%ele(i1)%a%eta)*0.5 if(plane == y_plane$) value = (lat_pos%ele(i1)%b%eta + lat_ele%ele(i1)%b%eta)*0.5 elseif (variable == d_xemit$) then value = global%mode_electron%a%emittance - global%wig%x_emit elseif (variable == elec_emitt$) then value = global%mode_electron%a%emittance elseif (variable == d_curly_D$) then value = global%mode_electron%a%j_damp - global%mode_positron%a%j_damp elseif (variable == pos_feedbk_phase$ .or. variable == ele_feedbk_phase$) then s_vert = 55.09 !location of vertical kicker s_horz = 56.55 !location of horizontal kicker s_det = 25.18 !location of detector 5W call twiss_and_track_at_s(lat,s_vert,ele_vert) call twiss_and_track_at_s(lat,s_horz,ele_horz) call twiss_and_track_at_s(lat,s_det,ele_det) ele_det%a%phi = lat%ele(i1)%a%phi !element i1 is detector ele_det%b%phi = lat%ele(i1)%b%phi n_turn = int(target_value)+1 if(variable == pos_feedbk_phase$)then if(plane == x_plane$) & value = (ele_horz%a%phi - ele_det%a%phi + n_turn * lat%a%tune)/twopi if(plane == y_plane$) & value = (ele_vert%b%phi - ele_det%b%phi + n_turn * lat%b%tune)/twopi endif if (variable == ele_feedbk_phase$) then if(plane == x_plane$) & value = (-ele_horz%a%phi + ele_det%a%phi + n_turn * lat%a%tune)/twopi if(plane == y_plane$) & value = (-ele_vert%b%phi + ele_det%b%phi + n_turn * lat%b%tune)/twopi endif value = value - int(value) elseif (variable == delta_beta_quad$) then call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == delta_beta$) then sextupole_constraint = .true. if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == delta_alpha$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == delta_phi$) then call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == det_2by2_ul$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == det_2by2_lr$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == det_4by4$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == trace$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) elseif (variable == dbeta_dpretz$) then sextupole_constraint = .true. if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif call max_or_min_delta_beta(con, lat_pretz_off, nonlin, moment, i, value, location, weight) elseif (variable == sex_k2_4$) then sextupole_constraint = .true. if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif val_max = 0 con%c(i)%contribution = 0 do j = 1, indep_var%n_var ix = indep_var%v(j)%rindex k2l = lat%ele(ix)%value(k2$) * lat%ele(ix)%value(l$) con%c(i)%contribution = con%c(i)%contribution + k2l**4 if (k2l >= val_max) then location = ix value = k2l endif enddo con%c(i)%contribution = con%c(i)%contribution * weight**con%fom_power elseif (variable == sext_moments$) then sextupole_constraint = .true. value = con%c(i)%actual_value if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == sext_resonances$) then sextupole_constraint = .true. value = con%c(i)%actual_value if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == sext_symmetry$) then sextupole_constraint = .true. value = con%c(i)%actual_value if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == sext_antisymmetry$) then sextupole_constraint = .true. value = con%c(i)%actual_value if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == chi2$) then sextupole_constraint = .true. value = sqrt(moment%chi2) elseif (variable == sign_symmetry$) then sextupole_constraint = .true. value = con%c(i)%actual_value if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == dbeta_dcos$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == dbeta_dsin$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == coupling_a_real$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == coupling_a_image$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == coupling_b_real$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == coupling_b_image$) then sextupole_constraint = .true. call max_or_min_delta_beta(con, lat, nonlin, moment, i, value, location, weight) if (con%minimize_moments) then if (constraint /= target$) call target_error (con%c(i)) constraint = int_garbage$ endif elseif (variable == sync_beta_path$ )then call sync_beta_path(lat, d_amp_x, d_amp_y) if(plane == x_plane$) value = d_amp_x if(plane == y_plane$) value = d_amp_y elseif (variable == sync_beta_volt$ )then call sync_beta_volt(lat, d_amp_x, d_amp_y) if(plane == x_plane$) value = d_amp_x if(plane == y_plane$) value = d_amp_y elseif (variable == alpha_xx$)then value = axx elseif (variable == alpha_xy$ )then value = axy elseif (variable == alpha_yy$ )then value = ayy elseif (variable == eta_x2$)then ! return rms for now etax2_rms = 0. jsh_count = 0 do jshx=1,ubound(etax2,1) if (etax2(jshx) .eq. 0.) cycle jsh_count = jsh_count + 1 ! filter out all elements where etax2 not computed etax2_rms = etax2_rms + (etax2(jshx))**2 !*value_of_attribute(lat_pos%ele(jshx),'L'))**2 enddo etax2_rms = sqrt(etax2_rms) /jsh_count !/ lat_pos%ele(lat_pos%n_ele_track)%s value = etax2_rms elseif (variable == beta_x1$)then betax1_rms = 0. jsh_count = 0 do jshx=1,ubound(betax1,1) if (betax1(jshx) .eq. 0.) cycle jsh_count = jsh_count + 1 ! filter out all elements where betax1 not computed betax1_rms = betax1_rms + (betax1(jshx))**2 !*value_of_attribute(lat_pos%ele(jshx),'L'))**2 enddo betax1_rms = sqrt(betax1_rms) /jsh_count !/ lat_pos%ele(lat_pos%n_ele_track)%s value = betax1_rms elseif (variable == beta_y1$)then betay1_rms = 0. jsh_count = 0 do jshx=1,ubound(betay1,1) if (betay1(jshx) .eq. 0.) cycle jsh_count = jsh_count + 1 ! filter out all elements where betay1 not computed betay1_rms = betay1_rms + (betay1(jshx))**2 !*value_of_attribute(lat_pos%ele(jshx),'L'))**2 enddo betay1_rms = sqrt(betay1_rms) /jsh_count !/ lat_pos%ele(lat_pos%n_ele_track)%s value = betay1_rms elseif (h11001r$ == variable )then value = real(h_nonlin(1)) elseif (h00111r$ == variable )then value = real(h_nonlin(2)) elseif (h20001r$ == variable )then value = real(h_nonlin(3)) elseif (h00201r$ == variable )then value = real(h_nonlin(4)) elseif (h10002r$ == variable )then value = real(h_nonlin(5)) elseif (h21000r$ == variable )then value = real(h_nonlin(6)) elseif (h30000r$ == variable )then value = real(h_nonlin(7)) elseif (h10110r$ == variable )then value = real(h_nonlin(8)) elseif (h10020r$ == variable )then value = real(h_nonlin(9)) elseif (h10200r$ == variable )then value = real(h_nonlin(10)) elseif (h11001i$ == variable )then value = aimag(h_nonlin(1)) elseif (h00111i$ == variable )then value = aimag(h_nonlin(2)) elseif (h20001i$ == variable )then value = aimag(h_nonlin(3)) elseif (h00201i$ == variable )then value = aimag(h_nonlin(4)) elseif (h10002i$ == variable )then value = aimag(h_nonlin(5)) elseif (h21000i$ == variable )then value = aimag(h_nonlin(6)) elseif (h30000i$ == variable )then value = aimag(h_nonlin(7)) elseif (h10110i$ == variable )then value = aimag(h_nonlin(8)) elseif (h10020i$ == variable )then value = aimag(h_nonlin(9)) elseif (h10200i$ == variable )then value = aimag(h_nonlin(10)) elseif (variable == h31000a$ )then value = abs(h_nonlin2(4)) elseif (variable == h11200a$ )then value = abs(h_nonlin2(5)) elseif (variable == h40000r$ )then value = real(h_nonlin2(6)) elseif (variable == h40000i$ )then value = aimag(h_nonlin2(6)) elseif (variable == h20020a$ )then value = abs(h_nonlin2(7)) elseif (variable == h20110a$ )then value = abs(h_nonlin2(8)) elseif (variable == h20200a$ )then value = abs(h_nonlin2(9)) elseif (variable == h00310a$ )then value = abs(h_nonlin2(10)) elseif (variable == h00400a$ )then value = abs(h_nonlin2(11)) elseif (variable == jx_var$ )then value = jvar(1) elseif (variable == jy_var$ )then value = jvar(2) elseif (variable == ptc_axx$ )then value = ptc_fps_param(1) elseif (variable == ptc_g_real$ )then value = ptc_fps_param(2) elseif (variable == ptc_g_imag$ )then value = ptc_fps_param(3) elseif (dnux_dp2$ == variable)then value = dnux_dp2 elseif (dnuy_dp2$ == variable)then value = dnuy_dp2 elseif (Qs_2Qx$ == variable )then value = res_amp elseif (Qs_2Qx_cos$ == variable )then value = res_cos elseif (Qs_2Qx_sin$ == variable )then value = res_sin elseif (Qs_2Qx_ele$ == variable )then value = res_amp_ele elseif (Qs_2Qx_cos_ele$ == variable )then value = res_cos_ele elseif (Qs_2Qx_sin_ele$ == variable )then value = res_sin_ele elseif (variable == Qs_2Qx_track$ )then value = max_x elseif (variable == variable_limit$) then vmax = 0 contrib = 0 value = 0 do k=1, indep_var%n_var x0 = var_attribute(lat, indep_var%v(k)) if (x0 < indep_var%v(k)%min_val) then v = indep_var%v(k)%min_val - x0 contrib = contrib + (v * weight)**2 elseif (x0 > indep_var%v(k)%max_val) then v = x0 - indep_var%v(k)%max_val contrib = contrib + (v * weight)**2 else v = 0 endif if (v > vmax) then vmax = v location = indep_var%v(k)%rindex value = x0 endif end do do k=1,sext_var%n_var x0 = var_attribute(lat, sext_var%v(k)) if (x0 < sext_var%v(k)%min_val) then v = sext_var%v(k)%min_val - x0 contrib = contrib + (v * weight)**2 elseif (x0 > sext_var%v(k)%max_val) then v = x0 - sext_var%v(k)%max_val contrib = contrib + (v * weight)**2 else v = 0 endif if (v > vmax) then vmax = v location = sext_var%v(k)%rindex value = x0 endif end do contrib = (max(0.0, vmax) * weight)**2 elseif (variable /= i_bunch$) then print *, 'ERROR IN FOM: UNRECOGNIZED VARIABLE: ', variable, i if (variable > 0) print *, var_name(variable) call err_exit endif ! figure figure of merit con%c(i)%actual_value = value if (variable /= variable_limit$) then if (con%c(i)%use_amplitude) value = abs(value) delta = value - target_value ! normalize by sigma if used if (con%use_sigma_norm) delta = delta / con%c(i)%target_sigma if (constraint == minz$) contrib = min(0.0, delta) if (constraint == maxz$) contrib = max(0.0, delta) if (constraint == target$) contrib = delta if (constraint == int_garbage$) contrib = con%c(i)%contribution if (con%c(i)%use_relative) contrib = contrib / con%c(i)%target_value endif if (variable == lrbbi_del_q$ .or. variable == w_t$ .or. & variable == i_bunch$) contrib = 0. ! con%c(i)%actual_value = value con%c(i)%location = location if (sextupole_constraint .and. .not. con%use_sextupole_fom) then con%c(i)%contribution = 0.0 else con%c(i)%contribution = contrib endif figure = figure + (contrib * weight)**con%fom_power enddo main_loop !--------------------------------------------------------------------------- ! get ordered list of contrib(i) and map i(j) where i(1) is max contrib(i) etc. do i = 1, con%n_constraint variable = con%c(i)%variable if (variable == lrbbi_del_q$ .or. variable == w_t$ .or. & variable == i_bunch$ .or. & variable == lrbbi_del_q_x$ .or. variable == lrbbi_del_q_y$) then constraint = con%c(i)%constraint if (variable == lrbbi_del_q$ .or. variable == w_t$ .or. & variable == lrbbi_del_q_x$ .or. variable == lrbbi_del_q_y$) then value = con%c(i)%actual_value I_max = min(I_x,I_y,I_wt) con%c(i)%actual_value = value*I_max elseif (variable == i_bunch$) then I_max = min(I_x,I_y,I_wt) value = I_max con%c(i)%actual_value = value endif delta = con%c(i)%actual_value - con%c(i)%target_value if (con%use_sigma_norm) delta = delta / con%c(i)%target_sigma weight = con%c(i)%weight if (constraint == minz$) contrib = min(0.0, delta) if (constraint == maxz$) contrib = max(0.0, delta) if (constraint == target$) contrib = delta con%c(i)%contribution = contrib figure = figure + (contrib * weight)**con%fom_power endif end do ! put everything together con%figure_of_merit = figure con%figure_of_merit_minop = figure*1.e-4*fact deallocate(cbar_ele) ! write(11,*)' figure of merit ',figure return end subroutine fom !........................................................................ !........................................................................ ! real(rp) b(indep_var_maxx,1) ! ! call showme_topten(sp1, lat, 0) ! b(1:n_par,1) = target_tem(1:n_par) - value(1:n_par) ! call gaussj(matrix(1:n_par,1:n_par), b(1:n_par,:)) subroutine linear_opt (lat, indep_var, cesr, sp1, pc, quad, & co, global, twiss_max, pretz, nonlin, moment1, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (indep_var_struct) comp_indep_var type (an_indep_var) thevar type (lat_struct) lat type (constraint_struct) sp1 type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment1 real(rp) vdel(indep_var_maxx), vec(indep_var_maxx) real(rp) fig_start, fig_end real(rp) value(indep_var_maxx), target(indep_var_maxx) ! real(rp) b(indep_var_maxx,1) real(rp) b(40), x(40) real(rp) Matrix(40, 40) real(rp) target_tem(indep_var_maxx), value_fin(indep_var_maxx) real(rp) comp_vec(indep_var_maxx+3) real(rp), pointer :: a_ptr integer n, n_par, i, j integer k character*80 fmt character*60 new_bmad_file ! if(sp1%n_constraint /= indep_var%n_var) then print *, ' Number of variables is not equal to the number of constraints ' print *, ' in LINEAR' call err_exit endif n_par = indep_var%n_var do n=1, n_par vec(n) = var_attribute(lat, indep_var%v(n)) enddo vdel(1:n_par) = indep_var%v(1:n_par)%delta call calcf (lat, sp1, pc, quad, co, global, twiss_max, pretz, nonlin, & indep_var) call fom (lat, cesr, sp1, pc, quad, co, global, twiss_max, pretz, & nonlin, moment1, indep_var, sext_var) fig_start = sp1%figure_of_merit ! call showme_topten(sp1, lat, 0) target(1:n_par) = sp1%c(1:n_par)%target_value value_fin(1:n_par) = sp1%c(1:n_par)%actual_value print '(1x, a, i3)', 'cycle = ', j do n = 1, n_par thevar = indep_var%v(n) print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', var_attribute(lat, indep_var%v(n)) enddo print *,' value target ' do n = 1, n_par print '(i6,2e12.4)', n, sp1%c(n)%actual_value, target(n) enddo do k = 1,sp1%n_loops value_fin(1:n_par) = sp1%c(1:n_par)%actual_value do j = 1, sp1%n_cycles value(1:n_par) = sp1%c(1:n_par)%actual_value do n=1, n_par vec(n) = var_attribute(lat, indep_var%v(n)) enddo target_tem(1:n_par) = value_fin(1:n_par) + & (target(1:n_par) - value_fin(1:n_par))/(sp1%n_cycles-j+1) if(k == sp1%n_loops)target_tem(1:n_par) = target(1:n_par) ! b(1:n_par,1) = target_tem(1:n_par) - value(1:n_par) b(1:n_par) = target_tem(1:n_par) - value(1:n_par) do n = 1, n_par a_ptr => var_attribute(lat, indep_var%v(n)) a_ptr = vec(n) + vdel(n) call calcf (lat, sp1, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom (lat, cesr, sp1, pc, quad, co, global, twiss_max, pretz, nonlin, moment1, indep_var, sext_var) do i = 1,n_par Matrix(n,i) = (sp1%c(i)%actual_value - value(i))/vdel(n) end do a_ptr = vec(n) end do call gjdet(matrix, b, x, n_par) ! call gaussj(matrix(1:n_par,1:n_par), b(1:n_par,:)) do n= 1, n_par a_ptr => var_attribute(lat, indep_var%v(n)) a_ptr = vec(n) + x(n) enddo call calcf(lat, sp1, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, sp1, pc, quad, co, global, twiss_max, pretz, nonlin, moment1, indep_var, sext_var) print '(1x, a, i3)', 'cycle = ', j do n = 1, n_par thevar = indep_var%v(n) print *, n, ' ', thevar%ele_name, thevar%val_name, ' ', var_attribute(lat, indep_var%v(n)) enddo print *,' value target ' do n = 1, n_par print '(i6,2e12.4)', n, sp1%c(n)%actual_value, target_tem(n) enddo end do end do if (sp1%plot) call plotdo_bmadz('/xs', lat, co, pc) do n= 1, n_par vec(n) = var_attribute(lat, indep_var%v(n)) enddo if(sp1%new_bmad_file <= ' ') then new_bmad_file = 'new_bmad_input.' else new_bmad_file = sp1%new_bmad_file endif if(sp1%sol_com > 0 .and. sp1%sol_com <= 2) then comp_vec(1:indep_var%n_var)= vec(1:indep_var%n_var) call compensation_indep_var(sp1%sol_com, sp1%compensating_ele, lat, indep_var, comp_indep_var, comp_vec) call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, comp_vec(1:comp_indep_var%n_var)) else call opt_out(lat%input_file_name,new_bmad_file,indep_var, vec) endif ! call opt_out(lat%input_file_name,new_bmad_file,indep_var, vec) call showme_topten(sp1, lat, 0) return end subroutine !----------------------------------------------------------------- !----------------------------------------------------------------- !----------------------------------------------------------------- subroutine manual_opt (lat, indep_var, cesr, con, pc, quad, & co, global, twiss_max, pretz, nonlin, moment, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (indep_var_struct) comp_indep_var type (an_indep_var) thevar type (lat_struct) lat type (constraint_struct) con type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment type (ele_struct) ele character*80 line, last_line character*4 plot_type character*60 new_bmad_file logical :: keep_trying = .true., hardcopy_flag = .false. real*8 vec(indep_var_maxx) real*8 comp_vec(indep_var_maxx+3) integer ix, n ! call calcf (lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, & indep_var) call fom (lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) call showme_topten(con, lat, 0) plot_type = ' ' if (con%plot) call plotdo_bmadz('/xs', lat, co, pc) do while (keep_trying) print '(a, $)', ' BMADZ> ' read(*,'(a)') line ix = index(line, '!') if (ix /= 0) line = line(:ix-1) ! strip off comments call str_upcase(line, line) call string_trim(line, line, ix) if (ix == 0) then ! nothing typed. do the same thing line = last_line endif last_line = line if(line(1:4) == 'EXIT')then return elseif(line(1:2) == 'WR')then !write new input do n = 1, indep_var%n_var vec(n) = var_attribute(lat, indep_var%v(n)) print *,' manual_opt:' ,n,vec(n) enddo if(con%new_bmad_file <= ' ') then new_bmad_file = 'new_bmad_input.' else new_bmad_file = con%new_bmad_file endif if(con%sol_com >= 0 .and. con%sol_com <= 2) then comp_vec(1:indep_var%n_var)= vec(1:indep_var%n_var) call compensation_indep_var(con%sol_com, con%compensating_ele, lat, indep_var, comp_indep_var, comp_vec) call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, comp_vec(1:comp_indep_var%n_var)) else call opt_out(lat%input_file_name,new_bmad_file,indep_var, vec) endif ! call opt_out(lat%input_file_name,new_bmad_file,indep_var, vec) print *,' New_bmad_input. created.' elseif(line(1:2) == 'SH')then call showme_topten(con, lat, 0) else call find_change( line, lat) call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) if (con%plot) call plotdo_bmadz('/xs', lat, co, pc) call showme_topten(con, lat, 0) endif end do return end subroutine manual_opt !----------------------------------------------------------------- !----------------------------------------------------------------- !----------------------------------------------------------------- subroutine find_change_backup(line, lat) implicit none type (lat_struct) lat type (ele_struct) ele real(rp) value character*16 attribute character line*(*) integer ix, n, ix_ele, ix_attr, change integer :: delta$ = 0, set_to$ = 1 integer ele_key ! ix = index(line, '!') if (ix /= 0) line = line(:ix-1) ! strip off comments call str_upcase(line, line) call string_trim(line, line, ix) if(line(1:2) == 'HE' .or. line(1:1) == '?')then print *,' Type element name (Q04W), <= or delta>,, and ' else ! first word is the lat element ix_ele = -1 do n=1, lat%n_ele_max if(lat%ele(n)%name == line(1:ix))then ix_ele = n ele_key = lat%ele(n)%key ele = lat%ele(n) endif end do if(ix_ele < 0)then print *,' Cannot find element ',line(1:ix), ' Try again' return endif call string_trim(line(ix+1:), line, ix) !second word is = or delta if(line(1:1) == '=')change=set_to$ if(line(1:1) == 'D')change=delta$ call string_trim(line(ix+1:), line, ix) !third word is attribute attribute = line(1:ix) ele%key = ele_key ix_attr = attribute_index(ele, attribute) if(ix_attr == 0)then call type_ele (ele, .true., 0, .false., 0) print *,' Attributes of ',ele%name,' shown above. Try again' return endif call string_trim(line(ix+1:), line, ix) !fourth word is value read(line(1:ix), *, err=10) value if(change == set_to$)lat%ele(ix_ele)%value(ix_attr) = value if(change == delta$)lat%ele(ix_ele)%value(ix_attr) = & lat%ele(ix_ele)%value(ix_attr) + value print '(1x, a, a8, a, a4, a, es20.10)', & 'Element = ', lat%ele(ix_ele)%name, ' Change ', attribute, & ' to ', lat%ele(ix_ele)%value(ix_attr) return 10 print *, ' Bad value ' endif return end subroutine !........................................................................ ! ! Subroutine : max_or_min(CON, LAT, QUAD, i_con, value, location) ! ! Description: If the constrained parameter is beta,alpha,eta, or etap this ! routine will find the maximum or minimum or target of the ! parameter in the range ! ! Arguments : ! INPUT: CON -- constraints structure ! QUAD ! i_con -- constraint number(integer) ! ! OUTPUT: CON.c(i_con).actual_value (real), value of the parameter ! of the i_con constraint ! CON.c(i_con).location (integer), LAT element number where ! minimum or maximum or target value of constrained ! parameter appears ! !........................................................................ subroutine max_or_min(CON, LAT, QUAD, i_con, value, location) implicit none type (lat_struct) lat type (zquad_struct) quad type (constraint_struct) con type (twiss_struct) t1, t2 type (beta_struct) xy integer i, i1, i2, j, i_con, location real(rp) param real(rp) value real(rp) curly_h logical first_min/.true./,first_max/.true./ ! i1 = con%c(i_con)%where1_ix i2 = con%c(i_con)%where2_ix if (i1 > i2) then print *, 'ERROR IN MAX_OR_MIN: LOCATION1 AFTER LOCATION2!' print *, 'Variable:', var_name(con%c(i_con)%variable) print *, 'Location/index 1: ', con%c(i_con)%where1, i1 print *, 'Location/index 2: ', con%c(i_con)%where2, i2 call err_exit endif if (i1 == 0) then print *, 'ERROR IN MAX_OR_MIN: LOCATION1 =0!' print *, 'Variable:', var_name(con%c(i_con)%variable) print *, 'Location/index 1: ', con%c(i_con)%where1, i1 print *, 'Location/index 2: ', con%c(i_con)%where2, i2 call err_exit endif ! location = -1 ! dummy do i = i1, i2 j = lat%ele(i)%ix_pointer if (j /= 0 .or. i == i1) then ! if (j == 0 .or.con%c(i_con)%variable == curly_h$ .or. & ! con%c(i_con)%variable == beta_electrons$ .or. con%c)then if (con%c(i_con)%plane == x_plane$) then t1 = lat%ele(i-1)%a t2 = lat%ele(i)%a else t1 = lat%ele(i-1)%b t2 = lat%ele(i)%b endif xy%beta_ave = (t1%beta + t2%beta) / 2 xy%beta_max = max(t1%beta, t2%beta) xy%eta_max = max(t1%eta, t2%eta) curly_h = (t1%beta*t1%etap**2 + t1%alpha * t1%eta * t1%etap + t1%gamma * t1%eta**2 + & t2%beta*t2%etap**2 + t2%alpha * t2%eta * t2%etap + t2%gamma * t2%eta**2) / 2 ! else ! if (con%c(i_con)%plane == x_plane$) then ! xy = quad%lens(j)%a ! elseif (con%c(i_con)%plane == y_plane$) then ! xy = quad%lens(j)%b ! else ! call plane_error (con%c(i)) ! endif ! endif if (con%c(i_con)%variable == beta$ .or. con%c(i_con)%variable == beta_electrons$) then param = xy%beta_max elseif (con%c(i_con)%variable == eta$ .or. con%c(i_con)%variable == eta_electrons$) then param = xy%eta_max elseif (con%c(i_con)%variable == curly_h$)then param = curly_h else print *, 'ERROR IN MAX_OR_MIN: I CANNOT HANDLE THIS CONSTRAINT!' print *, 'Variable:', var_name(con%c(i_con)%variable) call err_exit endif ! if (con%c(i_con)%constraint == minz$) then if(isnan(value)) value = 99. if (param < value .or. location < 0) then value = param location = i endif elseif (con%c(i_con)%constraint == maxz$) then if(isnan(value)) value=-99. if (param > value .or. location < 0) then value = param location = i endif else print *, 'ERROR IN MAX_OR_MIN: CONSTRAINT TYPE IS NOT MAX OR MIN!' print *, 'Variable:', var_name(con%c(i_con)%variable) print *, 'Constraint:', con%c(i_con)%constraint call err_exit endif endif ! j ne 0 end do return end subroutine max_or_min !+ ! If the constrained parameter is beta,alpha,eta, or etap this ! routine will find the maximum or minimum or target of the parameter ! in the range for electron beam ! ! INPUT: CON -- constraints structure ! QUAD ! i_con -- constraint number(integer) ! nonlin -- nonlin_struct ! ! OUTPUT: CON.c(i_con).actual_value (real), value of the parameter of the i_con constraint ! CON.c(i_con).location (integer), LAT element number where minimum or maximum ! or target value of constrained parameter ! appears !- subroutine max_or_min_electrons(CON, LAT, nonlin, QUAD, i_con, value, location) implicit none type (lat_struct) lat type (zquad_struct) quad type (constraint_struct) con type (twiss_struct) t1, t2 type (beta_struct) xy type (nonlin_ele_struct) nonlin integer i, i1, i2, j, i_con, location real(rp) value, param real(rp) etap_max ! i1 = con%c(i_con)%where1_ix i2 = con%c(i_con)%where2_ix if (i1 > i2) then print *, 'ERROR IN MAX_OR_MIN: LOCATION1 AFTER LOCATION2!' print *, 'Variable:', var_name(con%c(i_con)%variable) print *, 'Location/index 1: ', con%c(i_con)%where1, i1 print *, 'Location/index 2: ', con%c(i_con)%where2, i2 call err_exit endif ! location = -1 ! dummy do i = i1, i2 j = lat%ele(i)%ix_pointer if (j /= 0 .or. i == i1) then ! if (j == 0) then if (con%c(i_con)%plane == x_plane$) then t1 = nonlin%nonele_ele(i-1)%tx t2 = nonlin%nonele_ele(i)%tx else t1 = nonlin%nonele_ele(i-1)%ty t2 = nonlin%nonele_ele(i)%ty endif xy%beta_max = max(t1%beta, t2%beta) xy%eta_max = max(t1%eta, t2%eta) etap_max = max(t1%etap, t2%etap) ! else ! if (con%c(i_con)%plane == x_plane$) then ! xy = quad%lens(j)%a ! elseif (con%c(i_con)%plane == y_plane$) then ! xy = quad%lens(j)%b ! else ! call plane_error (con%c(i)) ! endif ! endif if (con%c(i_con)%variable == beta_electrons$) then param = xy%beta_max elseif (con%c(i_con)%variable == eta_electrons$) then param = xy%eta_max elseif (con%c(i_con)%variable == etap_electrons$) then param = etap_max else print *, 'ERROR IN MAX_OR_MIN: I CANNOT HANDLE THIS CONSTRAINT!' print *, 'Variable:', var_name(con%c(i_con)%variable) call err_exit endif ! if (con%c(i_con)%constraint == minz$) then if (param < value .or. location < 0) then value = param location = i endif elseif (con%c(i_con)%constraint == maxz$) then if (param > value .or. location < 0) then value = param location = i endif else print *, 'ERROR IN MAX_OR_MIN: CONSTRAINT TYPE IS NOT MAX OR MIN!' print *, 'Variable:', var_name(con%c(i_con)%variable) print *, 'Constraint:', con%c(i_con)%constraint call err_exit endif endif ! j ne 0 end do return end subroutine max_or_min_electrons !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine minop_b(lat, indep_var, cesr, con, pc, quad, co, & global, twiss_max, pretz, nonlin, maxfun, idisp, step, & acc, init, moment, sext_var) implicit none type (indep_var_struct) indep_var, sext_var type (indep_var_struct) comp_indep_var type (constraint_struct) con type (lat_struct) lat type (db_struct) cesr type (pc_struct) pc type (twiss_max_struct) twiss_max type (pretz_struct) pretz type (global_struct) global type (zquad_struct) quad type (coord_struct), allocatable :: co(:) type (nonlin_ele_struct) nonlin type (moment_struct) moment ! RECORD OF LARGEST AND NEXT LARGEST CONTRIB TO F real(rp) X(indep_var_maxx), G(INDEP_VAR_MAXX),XA(INDEP_VAR_MAXX) real(rp) x0(indep_var_maxx) real(rp) GA(INDEP_VAR_MAXX),H(INDEP_VAR_MAXX,INDEP_VAR_MAXX) real(rp) GG(INDEP_VAR_MAXX,INDEP_VAR_MAXX) real(rp) V(INDEP_VAR_MAXX),U(INDEP_VAR_MAXX), WA(INDEP_VAR_MAXX) real(rp) WB(INDEP_VAR_MAXX), WC(INDEP_VAR_MAXX), WT(INDEP_VAR_MAXX) real(rp) WS(INDEP_VAR_MAXX), Y(INDEP_VAR_MAXX), Z(INDEP_VAR_MAXX) real(rp) GZ(INDEP_VAR_MAXX), GV(INDEP_VAR_MAXX), PZ(INDEP_VAR_MAXX) real(rp) f,step,acc,fprev,dss,gsq, c, ggdiag, hdiag, epssq, wasq, gggg real(rp) hggg, d, ca, cb, theta, fa, dg, dga, dggd, wbsq, b0, dsq real(rp) flag, c2, c3, sum, c1, zgv, vgv, zwc, zgz, vwc, vsq, e real(rp) zv, uv, usq, sigma real(rp) a,b,cc,phi real(rp), pointer :: a_ptr real*8 x8(indep_var_maxx) real*8 comp_vec(indep_var_maxx+3) integer maxfun, nfcall,ngcall,iptest,i,j, iupd integer idisp integer lastout character*60 new_bmad_file logical init ! ****************************************************************** ! * INITIATLIZE THE FOLLOWING PARAMETERS * ! * NFCALL & NGCALL = NO. OF FUNCTION & GRADIENT CALLS RESP.. * ! ****************************************************************** lastout=0 NFCALL = 1 FPREV=0 !RECORD LAST F NGCALL = 1 IUPD = 1 EPSSQ = ACC * ACC IPTEST = 1 DSS = STEP*STEP do i = 1, indep_var%n_var x0(i) = 1.0 x(i) = var_attribute(lat, indep_var%v(i))/x0(i) x8(i)=x(i) end do call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) f = con%figure_of_merit CALL BMAD_CALCG (lat, con, pc, quad, co, global, & twiss_max, pretz, cesr, nonlin, indep_var, G, moment, sext_var) IF(.NOT.INIT) GO TO 8 !SAVE HESSIAN IF NOT FIRST RUN GSQ = 0 DO 3 I = 1,indep_var%n_var 3 GSQ = GSQ + G(I) * G(I) ! ! ****************************************************************** ! * PREPARE FOR THE FIRST ITERATION BY INITIALIZING MATRICES H&GG * ! ****************************************************************** ! C = -SQRT( DSS / GSQ ) GGDIAG = 0.01 * SQRT( GSQ) / STEP HDIAG = 1. / GGDIAG DO 6 I= 1,indep_var%n_var DO 5 J = 1,indep_var%n_var H(I,J) = 0.0D0 GG(I,J) = 0.0D0 5 CONTINUE H(I,I) = HDIAG GG(I,I) = GGDIAG 6 Z(I) = C * G(I) ! 8 GSQ = 0 DO 9 I = 1, indep_var%n_var 9 GSQ = GSQ + G(I) * G(I) IF(F>.00002.AND.1.05*FEPSSQ).AND.(F>.00006)) GO TO 18 RETURN ! RETURN ON GRADIENT, OR FUNCTION VALUE ! ****************************************************************** ! * TEST WHETHER MAXFUN CALLS OF CALCFG HAVE BEEN MADE * ! ****************************************************************** ! 18 FPREV=F IF(NFCALL==1)then call showme_topten(con, lat, 0) if(con%plot) call plotdo_bmadz('/xs', lat, co, pc) if(con%new_bmad_file<=' ') then new_bmad_file = 'new_bmad_input.' else new_bmad_file = con%new_bmad_file endif if(con%sol_com >= 0 .and. con%sol_com <= 2) then comp_vec(1:indep_var%n_var)= x8(1:indep_var%n_var) call compensation_indep_var(con%sol_com,con%compensating_ele, & lat, indep_var, comp_indep_var, comp_vec) call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, & comp_vec(1:comp_indep_var%n_var)) else call opt_out(lat%input_file_name,new_bmad_file,indep_var, x8) endif ! call opt_out(lat%input_file_name, new_bmad_file, indep_var, x8) endif IF(NFCALL>MAXFUN) RETURN !FUNC CALLS EXHAUSTED ! ****************************************************************** ! * CALCUALTE NEWTON POINT * ! ****************************************************************** 26 WASQ = 0.0D0 DO 28 I = 1,indep_var%n_var WA(I) = 0.0D0 DO 27 J= 1,indep_var%n_var 27 WA(I) = WA(I) - H(I,J) * G(J) 28 WASQ = WASQ + WA(I) * WA(I) IF ( WASQ <= DSS) GO TO 39 GGGG = 0 HGGG = 0.0D0 DO 30 I = 1, indep_var%n_var WB(I) = 0.0D0 DO 29 J = 1,indep_var%n_var 29 WB(I) = WB(I) + GG(I,J) * G(J) HGGG = HGGG - WA(I) * G(I) 30 GGGG = GGGG + WB(I) * G(I) ! D = HGGG * GGGG ! D = GSQ * GSQ / D D = (GSQ/HGGG)*(GSQ/GGGG) D = 0.2 + 0.8 * D IF ( D * D * WASQ > DSS ) GO TO 32 D = SQRT ( DSS / WASQ ) DO 31 I = 1,indep_var%n_var 31 WA(I) = D * WA(I) GO TO 41 32 CONTINUE ! IF ((GGGG *(GGGG * DSS))> GSQ ** 3 ) GO TO 34 IF (((GGGG/GSQ) *(GGGG/GSQ) * DSS)> GSQ ) GO TO 34 ! ! ****************************************************************** ! * THE CAUCHY POINT IS OUTSIDE CIRCLE, SO PICK THE POINT ON THE * ! * BOUNDARY * ! ****************************************************************** C = -SQRT ( DSS / GSQ) DO 33 I = 1,indep_var%n_var 33 WA(I) = C * G(I) GO TO 41 ! ! ****************************************************************** ! * THE CAUCHY POINT IS INSIDE CIRCLE, USE DOG LEG STEP. * ! ****************************************************************** 34 DO 35 I = 1,indep_var%n_var 35 WA(I) = D * WA(I) CA = 0.0 CB = 0 C = -GSQ / GGGG ! ****************************************************************** ! * SET THE STEEPEST DESCENT CORRECTION IN WB AND THE DIFFERENCE * ! * BETWEEN WA AND WB IN WC. * ! ****************************************************************** DO 36 I = 1, indep_var%n_var WB(I) = C * G(I) WC(I)= WA(I) - WB(I) CA = CA + WB(I) * WC(I) 36 CB = CB + WC(I) * WC(I) ! ! * FIND CORRECTION VECTOR AT THE INTERSECTION OF WA-WB AND CIRCLE * C = DSS -C * C * GSQ THETA = SIGN (C / (ABS(CA) + SQRT( CA * CA + C * CB)), CA) ! * TEST WHETHER TO USE THE GENERALIZED NEWTON STEP * IF (THETA-1.) 37, 39, 39 37 DO 38 I = 1, indep_var%n_var 38 WA(I) = WB(I) + THETA * WC(I) GO TO 41 39 CONTINUE ! 41 DO 42 I = 1,indep_var%n_var 42 XA(I) = X(I) + WA(I) NFCALL = NFCALL + 1 do i = 1, indep_var%n_var a_ptr => var_attribute(lat, indep_var%v(i)) a_ptr = xa(i)*x0(i) x8(i)=xa(i) end do call calcf(lat, con, pc, quad, co, global, twiss_max, pretz, nonlin, indep_var) call fom(lat, cesr, con, pc, quad, co, global, twiss_max, pretz, & nonlin, moment, indep_var, sext_var) fa = con%figure_of_merit 43 DSQ = 0.0D0 DO 44 I = 1,indep_var%n_var 44 DSQ = DSQ + WA(I) * WA(I) ! ****************************************************************** ! * TEST IF FUNCTION VALUE IS DECREASED * ! ****************************************************************** if(fa < f) then lastout=nfcall IF(((NFCALL-lastout)>=idisp)) then call showme_topten(con, lat, 0) if(con%plot) call plotdo_bmadz('/xs', lat, co, pc) endif if(con%new_bmad_file<=' ') then new_bmad_file = 'new_bmad_input.' else new_bmad_file = con%new_bmad_file endif if(con%sol_com >= 0 .and. con%sol_com <= 2) then comp_vec(1:indep_var%n_var)= x8(1:indep_var%n_var) call compensation_indep_var(con%sol_com, con%compensating_ele, & lat, indep_var, comp_indep_var, comp_vec) call opt_out(lat%input_file_name,new_bmad_file,comp_indep_var, & comp_vec(1:comp_indep_var%n_var)) else call opt_out(lat%input_file_name,new_bmad_file,indep_var, x8) endif ! call opt_out(lat%lattice, new_bmad_file, indep_var%n_var, x8) endif IF (FA < F) GO TO 46 IF(F<.00006) RETURN !MUST IMPROVE EACH ITER OR STOP DSS = 0.25 * DSQ GO TO 18 ! * CALCULATE SOME NUMBERS FOR REVISING STEP BOUND * 46 NGCALL = NGCALL + 1 do i = 1, indep_var%n_var a_ptr => var_attribute(lat, indep_var%v(i)) a_ptr = xa(i)*x0(i) end do CALL BMAD_CALCG (lat, con, pc, quad, co, global, twiss_max, & pretz, cesr, nonlin, indep_var, GA, moment, sext_var) DG = 0.0 DGA = 0.0 DGGD = 0.0 WBSQ = 0.0 DO 47 I = 1, indep_var%n_var WC(I) = 0.0D0 DO 47 J = 1,indep_var%n_var 47 WC(I) = WC(I) + GG(I,J) * WA(J) ! ! * SET (Y-GS) IN WB * DO 48 I = 1,indep_var%n_var U(I) = WC(I) WB(I) = GA(I) -G(I) - WC(I) DG = DG + G(I) * WA(I) DGA = DGA + GA(I) * WA(I) DGGD = DGGD + WC(I) * WA(I) 48 WBSQ = WBSQ + WB(I) *WB(I) ! ****************************************************************** ! * TEST WHETHER TO DECREASE THE STEP-BOUND * ! ****************************************************************** IF (FA-F-0.1 * DG-0.05 * DGGD) 51, 51, 50 50 DSS = 0.25 * DSQ IF(DSS>1E-15) GO TO 54 !GIVE UP, STEPS TOO SMALL print *, ' ESCAPED THE STEP SIZE TRAP....' RETURN ! ****************************************************************** ! * TEST WHETHER TO INCREASE THE STEP-BOUND * ! ****************************************************************** 51 DSS = DSQ IF (WBSQ - 0.25 * GSQ) 53, 53, 52 52 IF (DG -DGA -DGA) 54, 53, 53 53 DSS = 4.0 * DSQ ! * SET THE DIFFERENCE BETWEN GRADIENTS IN WC * 54 DO 55 I = 1,indep_var%n_var 55 WC(I) = GA(I) - G(I) F = FA B0 = 0.0 DO 56 I = 1,indep_var%n_var X(I) = XA(I) G(I) = GA(I) 56 B0 = B0 + WC(I) * WA(I) ! ****************************************************************** ! * TEST WHETHER WE NEED TO UPDATE H AND G * ! ****************************************************************** IF (B0 >= 1.0D-30) GO TO 58 ! WRITE(6,57) NFCALL ! 57 FORMAT('FUNCTION DECREASED ON',I4,'TH ITERATION, NO UPDATE IS MADE ! 1 ,SINCE Y*Z IS NEG. OR TOO SMALL') GO TO 8 58 IF (IUPD == 1) GO TO 63 IF (FLAG == 0.0) GO TO 60 DO 59 I = 1, indep_var%n_var 59 Z(I) = WA(I) GO TO 63 60 C2=0.0 C3 = 0.0 SUM = 0.0 DO 61 I = 1, indep_var%n_var C2 = C2 + WS(I) * WA(I) C3 = C3 + WT(I) * WA(I) 61 SUM = SUM +Y(I) * WA(I) DO 62 I = 1, indep_var%n_var 62 Z(I) = C1 *(C2 * Z(I) + C3 * V(I)) - (SUM/ZWC) * Z(I) 63 IUPD = IUPD+1 DO 64 I = 1, indep_var%n_var GV(I) = U(I) - WC(I) V(I) = WA(I) GZ(I) = 0.0 DO 64 J = 1, indep_var%n_var GZ(I) = GZ(I) + GG(I,J) * Z(J) 64 V(I) =V(I) - H(I,J) * WC(J) ZGV = 0.0 VGV = 0.0 ZWC = 0.0 ZGZ = 0.0 VWC = 0.0 DO 65 I = 1, indep_var%n_var ZWC = ZWC + Z(I) * WC(I) ZGZ = ZGZ + Z(I) * GZ(I) VWC = VWC + V(I) * WC(I) ZGV = ZGV + Z(I) * GV(I) VGV = VGV + V(I) * GV(I) 65 CONTINUE C1 = 1. / (ZGZ * VGV - ZGV * ZGV) C2 = ZWC * VGV - VWC * ZGV C3 = -ZWC * ZGV + VWC * ZGZ VSQ =0.0 ZV =0.0 A = 0.0 E = 0.0 DO 66 I = 1,indep_var%n_var Y(I) = C1 * (C2 * GZ(I) + C3 * GV(I)) WS(I) = VGV * GZ(I) - ZGV * GV(I) WT(I) = -ZGV * GZ(I) + ZGZ * GV(I) ZV = ZV + Z(I) * V(I) VSQ = VSQ + V(I) * V(I) A = A + WS(I) * WA(I) 66 E = E + WT(I) * WA(I) B = 0.0 DO 67 I = 1,indep_var%n_var PZ(I) = C1 * ( A * Z(I) + E * V(I)) 67 B = B +PZ(I) * Y(I) IF (B >= 1.0D-30) GO TO 70 B = B0 ! WRITE (6,68) 68 FORMAT(15X,'IMAGE PRODUCT NEG.OR TOO SMALL,ORIGINAL STEP USED') DO 69 I = 1, indep_var%n_var PZ(I) = WA(I) 69 GA(I) = WC(I) GO TO 72 70 DO 71 I = 1, indep_var%n_var 71 GA(I) = Y(I) 72 CC = 0.0 A = 0.0 DO 74 I = 1, indep_var%n_var WB(I) = 0.0 GZ(I) = 0.0 DO 73 J = 1, indep_var%n_var GZ(I) = GZ(I) + GG(I,J) *PZ(J) 73 WB(I) = WB(I) + H(I,J) * GA(J) A = A + GA(I) * WB(I) 74 CC = CC + PZ(I) * GZ(I) ! ****************************************************************** ! * OPTIMAL CONDITIONING * ! ****************************************************************** IF ( B <= 2.0D0 * A * CC / (A + CC)) GO TO 75 PHI = B / (B - A) GO TO 76 75 PHI = B * (CC - B)/(A * CC - B * B) 76 USQ = 0.0 UV = 0.0 ! ****************************************************************** ! * TEST WHETHER Z & V ARE LINEARLY INDEPENDENT * ! ****************************************************************** DO 77 I = 1, indep_var%n_var U(I) = Z(I) - ( ZV / VSQ ) * V(I) USQ = USQ + U(I) * U(I) 77 UV = UV + U(I) * V(I) FLAG = 0.0 IF (1.0D6 * UV * UV >= VSQ * USQ ) FLAG =1.0 ! ****************************************************************** ! * IF FLAG = 1 THEN Z & V ARE LINEARLY DEPENDENT * ! ****************************************************************** SIGMA = 1.+ PHI * (A * CC/ (B * B) - 1.) IF ( SIGMA < 1.0D-30) WRITE (6,78) 78 FORMAT(20X,' G IS ALMOST SINGULAR') DO 79 I = 1, indep_var%n_var WA(I) = (1./B) * PZ(I) - (1./A) * WB(I) WC(I) = (1./B) * GZ(I) -CC/(B * B) * GA(I) 79 U(I) = GA(I)- GZ(I) ! ****************************************************************** ! * START TO UPDATE MATRICES H & GG * ! ****************************************************************** DO 80 I = 1, indep_var%n_var DO 80 J = 1, indep_var%n_var GG(I,J) = GG(I,J) + (1./B) * (U(I) * GA(J) + GA(I) * U(J)) & -(( B- CC ) / (B * B)) * GA(I) * GA(J) & -(PHI* A / SIGMA) * WC(I) * WC(J) H(I,J) = H(I,J) - (1./A) * WB(I) * WB(J) + (1./B) * PZ(I) * PZ(J) & + PHI * A * WA(I) * WA(J) 80 CONTINUE GO TO 8 END SUBROUTINE ! subroutine plane_error (c) implicit none type (con_struct) c print *, 'ERROR IN FOM: PLANE IS NOT "X" NOR "Y"' print *, ' FOR CONSTRAINT: ', var_name(c%variable) print *, ' AT LOCATION: ', c%where1 call err_exit end subroutine plane_error !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine target_error (c) implicit none type (con_struct) c print *, 'ERROR IN FOM: FOR MINIMIZE_MOMENTS YOU MUST USE A "TARGET"' print *, ' CONSTRAINT FOR: ', var_name(c%variable) call err_exit end subroutine target_error end module opt_mod