subroutine wave_cbar (phase, tune, wave) use cesrv_struct use cesrv_interface implicit none type (d2_data_struct) phase type (d2_data_struct) tune type (wave_struct), target :: wave type (p1_plot_struct), pointer :: plot1 real(rp) phi_s_ring, phi_r_ring, phi_0s, phi_0r, phi_r_ave real(rp) cos_s(250), sin_s(250), cos_r(250), sin_r(250) real(rp) phi_s(250), phi_r(250) real(rp) coef_ba(4), coef_a(4), coef_b(4), rms_a(5), rms_b(5), rms_ba(5) real(rp) amp_sa, amp_sb, amp_ra, amp_rb, amp_sba, amp_rba, x integer i_a1, i_a2, i_b1, n_a, n_b integer m, m_min, m_max, nc, n_pts_min / 5 / integer i, ix, j, n_data_use character fmt*80 ! load VAR_ array from what is plotted and find typical value plot1 => wave%p2%plot1 call plotting_data_calc (plot1) phi_s_ring = tune%x%d(1)%model + tune%y%d(1)%model phi_r_ring = tune%x%d(1)%model - tune%y%d(1)%model n_data_use = plot1%n_use if (n_data_use == 0) then print *, 'ERROR: NO DATA FOR THE WAVE ANALYSIS.' return endif plot1%n_use = 0 do i = 1, n_data_use ix = wave%p2%ix(i) x = wave%p2%x(i) phi_s(i) = phase%x%d(ix)%model + phase%y%d(ix)%model phi_r(i) = phase%x%d(ix)%model - phase%y%d(ix)%model cos_s(i) = cos(phi_s(i)) sin_s(i) = sin(phi_s(i)) cos_r(i) = cos(phi_r(i)) sin_r(i) = sin(phi_r(i)) plot1%y(i+n_data_use) = plot1%y(i) wave%p2%ix(i+n_data_use) = wave%p2%ix(i) + 100 wave%p2%x(i+n_data_use) = wave%p2%x(i) + 100 phi_s(i+n_data_use) = phi_s(i) + phi_s_ring phi_r(i+n_data_use) = phi_r(i) + phi_r_ring cos_s(i+n_data_use) = cos(phi_s(i+n_data_use)) sin_s(i+n_data_use) = sin(phi_s(i+n_data_use)) cos_r(i+n_data_use) = cos(phi_r(i+n_data_use)) sin_r(i+n_data_use) = sin(phi_r(i+n_data_use)) if (x > 50 .and. plot1%n_use == 0) plot1%n_use = i - 1 + n_data_use enddo !----------------------------------------- ! Fit to regions... ! First: Find how many points are in the regions and where they start. ! Second: Find model coeficients COEF_A, and COEF_B. call wave_region_setup (wave, i_a1, i_b1, n_a, n_b) if (n_a < n_pts_min) then print *, 'ERROR: REGION A HAS LESS THAN', n_pts_min, ' POINTS' return elseif (n_b < n_pts_min .and. (wave%ix_b1 /= 0 .or. wave%ix_b2 /= 0)) then print *, 'ERROR: REGION B HAS LESS THAN', n_pts_min, ' POINTS' return endif call wave_anal (plot1%y(i_a1:), sin_s(i_a1:), cos_s(i_a1:), & sin_r(i_a1:), cos_r(i_a1:), 4, n_a, coef_a, rms_a) call wave_anal (plot1%y(i_b1:), sin_s(i_b1:), cos_s(i_b1:), & sin_r(i_b1:), cos_r(i_b1:), 4, n_b, coef_b, rms_b) wave%p2%plot2%n_use = plot1%n_use wave%p2%plot3%n_use = plot1%n_use do i = 1, plot1%n_use wave%p2%plot2%y(i) = plot1%y(i) - coef_a(1) * sin_s(i) - & coef_a(2) * cos_s(i) - coef_a(3) * sin_r(i) - coef_a(4) * cos_r(i) wave%p2%plot3%y(i) = plot1%y(i) - coef_b(1) * sin_s(i) - & coef_b(2) * cos_s(i) - coef_b(3) * sin_r(i) - coef_b(4) * cos_r(i) if (n_b < n_pts_min) wave%p2%plot3%y(i) = 0 enddo ! compute the kick if (n_b < n_pts_min) then wave%n_cross = 0 return endif do i = 1, 4 coef_ba(i) = coef_b(i) - coef_a(i) rms_ba(i) = sqrt(rms_a(i)**2 + rms_b(i)**2) enddo amp_sa = sqrt (coef_a(1)**2 + coef_a(2)**2) amp_ra = sqrt (coef_a(3)**2 + coef_a(4)**2) amp_sb = sqrt (coef_b(1)**2 + coef_b(2)**2) amp_rb = sqrt (coef_b(3)**2 + coef_b(4)**2) amp_sba = sqrt (coef_ba(1)**2 + coef_ba(2)**2) amp_rba = sqrt (coef_ba(3)**2 + coef_ba(4)**2) wave%amp_sba = amp_sba wave%amp_rba = amp_rba if (coef_ba(1) == 0 .and. coef_ba(2) == 0) then phi_0s = 0 else phi_0s = atan2 (coef_ba(1), coef_ba(2)) endif if (coef_ba(3) == 0 .and. coef_ba(4) == 0) then phi_0r = 0 else phi_0r = atan2 (coef_ba(3), coef_ba(4)) endif i_a2 = i_a1 + n_a - 1 ! end of region A m_min = int((phi_s(i_a2) - phi_0s)/pi) + 1 m_max = int((phi_s(i_b1) - phi_0s)/pi) nc = 0 do i = m_min, m_max nc = nc + 1 wave%phi_s(nc) = pi * i + phi_0s wave%kick(nc) = -2 * (coef_ba(1) * sin(wave%phi_s(nc)) + & coef_ba(2) * cos(wave%phi_s(nc))) do j = i_b1, i_a2, -1 if (phi_s(j) < wave%phi_s(nc)) then wave%ix_cross(nc) = wave%p2%x(j) exit endif enddo phi_r_ave = (phi_r(j) + phi_r(j+1)) / 2 m = nint((phi_r_ave - phi_0r) / pi) wave%phi_r(nc) = phi_0r + m * pi if (wave%phi_s(nc) > phi_s_ring) then wave%phi_s(nc) = wave%phi_s(nc) - phi_s_ring wave%phi_r(nc) = wave%phi_r(nc) - phi_r_ring endif enddo wave%n_cross = nc ! rms to amplitude ratios if (amp_ra == 0 .or. amp_sa == 0 .or. amp_rb == 0 .or. amp_sb == 0) then wave%chi_a = 0 wave%rms_ra = 0 wave%rms_sa = 0 wave%rms_sb = 0 wave%rms_rb = 0 wave%rms_ks = 0 wave%rms_kr = 0 wave%rms_sphi = 0 wave%rms_rphi = 0 else wave%chi_a = abs(amp_sba - amp_rba) / max(amp_sba + amp_rba, 1e-10) wave%rms_ra = sqrt(rms_a(3)**2 + rms_a(4)**2) / max(amp_ra, 1e-10) wave%rms_sa = sqrt(rms_a(1)**2 + rms_a(2)**2) / max(amp_sa, 1e-10) wave%rms_sb = sqrt(rms_b(1)**2 + rms_b(2)**2) / max(amp_sb, 1e-10) wave%rms_rb = sqrt(rms_b(3)**2 + rms_b(4)**2) / max(amp_rb, 1e-10) wave%rms_ks = sqrt(rms_ba(1)**2 * coef_ba(1)**2 + & rms_ba(2)**2 * coef_ba(2)**2) / max(amp_sba**2, 1e-10) wave%rms_kr = sqrt(rms_ba(3)**2 * coef_ba(3)**2 + & rms_ba(4)**2 * coef_ba(4)**2) / max(amp_rba**2, 1e-10) wave%rms_sphi = sqrt(rms_ba(2)**2 * coef_ba(1)**2 + & rms_ba(1)**2 * coef_ba(2)**2) / max(amp_sba**2, 1e-10) wave%rms_rphi = sqrt(rms_ba(4)**2 * coef_ba(3)**2 + & rms_ba(3)**2 * coef_ba(4)**2) / max(amp_rba**2, 1e-10) endif ! typeout info fmt = '(1x, a10, 6f9.4)' print *, 'NUMBER OF POSSIBLE KICKER LOCATIONS:', wave%n_cross print *, ' Lambda Rho Gamma Zeta', & ' Amp_s Amp_r' print fmt, 'A Region:', (coef_a(i), i = 1, 4), amp_sa, amp_ra print fmt, 'B Region:', (coef_b(i), i = 1, 4), amp_sb, amp_rb print fmt, 'B-A:', (coef_ba(i), i = 1, 4), amp_sba, amp_rba print * print *, ' rms_L rms_R rms_G rms_Z' print fmt, 'A Region:', (rms_a(i), i = 1, 4) print fmt, 'B Region:', (rms_b(i), i = 1, 4) print fmt, 'B-A:', (rms_ba(i), i = 1, 4) end subroutine