|
1 | 1 | from ocelot.cpbd.io import save_particle_array |
2 | | -from ocelot.common.globals import h_eV_s, m_e_eV, m_e_GeV, ro_e, speed_of_light |
3 | | -from ocelot.cpbd.beam import Twiss, beam_matching |
| 2 | +from ocelot.common.globals import h_eV_s, m_e_eV, m_e_GeV, ro_e, speed_of_light, q_e |
| 3 | +from ocelot.cpbd.beam import Twiss, beam_matching, global_slice_analysis, s_to_cur, get_envelope |
4 | 4 | from ocelot.utils.acc_utils import slice_bunching |
5 | 5 | from ocelot.common.ocelog import * |
6 | 6 | from ocelot.cpbd.beam import ParticleArray |
@@ -638,3 +638,118 @@ def apply(self, p_array, dz=0): |
638 | 638 | p_array.E = self.Eref |
639 | 639 | p_array.p()[:] = p_new[:] |
640 | 640 |
|
| 641 | + |
| 642 | +class IBS(PhysProc): |
| 643 | + """ |
| 644 | + Intrabeam Scattering (IBS) Physics Process. |
| 645 | +
|
| 646 | + This class models the intrabeam scattering process in particle beams. Two methods are implemented based on different formulations: |
| 647 | +
|
| 648 | + 1. **Huang Method:** Based on Z. Huang's work: *Intrabeam Scattering in an X-ray FEL Driver* (LCLS-TN-02-8, 2002). |
| 649 | + URL: https://www-ssrl.slac.stanford.edu/lcls/technotes/LCLS-TN-02-8.pdf |
| 650 | + 2. **Nagaitsev Method:** Based on S. Nagaitsev's work: *Intrabeam scattering formulas for fast numerical evaluation* |
| 651 | + (PRAB 8, 064403, 2005). URL: https://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.8.064403 |
| 652 | +
|
| 653 | + The methods can be selected using `self.method`, which accepts "Huang" or "Nagaitsev". |
| 654 | +
|
| 655 | + Key assumptions and features: |
| 656 | + - A **round beam approximation** is used, where `sigma_xy = (sigma_x + sigma_y) / 2`. |
| 657 | + - Beam parameters are calculated within a slice of width ±0.5 `sigma_tau` relative to the slice with maximum current |
| 658 | + (`self.slice = "Imax"` by default). |
| 659 | + - The Coulomb logarithm (`self.Clog`) is a configurable parameter, defaulting to 8 (based on Z. Huang's work). |
| 660 | + - A factor of 2 is applied to the Nagaitsev formula to align high-energy results with the Huang formula. |
| 661 | +
|
| 662 | + Attributes: |
| 663 | + - `method` (str): Selected method for IBS calculation ("Huang" or "Nagaitsev"). |
| 664 | + - `Clog` (float): Coulomb logarithm (default: 8) is used constant if update_Clog = False |
| 665 | + - `update_Clog` (bool): recalculate Clog on each step using Huang's formula without cut off. |
| 666 | + - `bounds` (list): Range of bounds for slice analysis in units of `sigma_tau` (default: `[-0.5, 0.5]`). |
| 667 | + - `slice` (str): Reference slice ("Imax" for maximum current by default). |
| 668 | +
|
| 669 | + Methods: |
| 670 | + - `get_beam_params(p_array)`: Computes beam parameters such as `sigma_xy`, `sigma_z`, and normalized emittance. |
| 671 | + - `apply(p_array, dz)`: Applies the IBS process to a particle array over a specified distance. |
| 672 | +
|
| 673 | + Raises: |
| 674 | + - `ValueError`: If an invalid method is specified. |
| 675 | + """ |
| 676 | + |
| 677 | + def __init__(self, method="Huang"): |
| 678 | + PhysProc.__init__(self) |
| 679 | + self.method = method |
| 680 | + self.Clog = 8 |
| 681 | + self.update_Clog = True |
| 682 | + self.bounds = [-0.5, 0.5] |
| 683 | + self.slice = "Imax" |
| 684 | + |
| 685 | + self.emit_n = None |
| 686 | + self.sigma_xy = None |
| 687 | + self.sigma_z = None |
| 688 | + self.sigma_dgamma0 = None |
| 689 | + |
| 690 | + def get_beam_params(self, p_array): |
| 691 | + """ |
| 692 | + Compute beam parameters such as sigma_xy, sigma_z, and normalized emittance. |
| 693 | +
|
| 694 | + :param p_array: ParticleArray |
| 695 | + Input particle array containing particle properties. |
| 696 | + """ |
| 697 | + tws = get_envelope(p_array, bounds=self.bounds, slice=self.slice) |
| 698 | + self.sigma_z = np.std(p_array.tau()) |
| 699 | + self.sigma_xy = (np.sqrt(tws.xx) + np.sqrt(tws.yy)) / 2. |
| 700 | + self.emit_n = (tws.emit_xn + tws.emit_yn) / 2. |
| 701 | + pc = np.sqrt(p_array.E ** 2 - m_e_GeV ** 2) |
| 702 | + self.sigma_dgamma0 = np.sqrt(tws.pp)*pc / m_e_GeV |
| 703 | + |
| 704 | + def estimate_Clog(self): |
| 705 | + """ |
| 706 | + Used formula from Hunag's paper without cut offs |
| 707 | +
|
| 708 | + :return: float - Coulomb logarithm |
| 709 | + """ |
| 710 | + Clog = np.log(self.emit_n * self.emit_n / (ro_e * self.sigma_xy)) |
| 711 | + return Clog |
| 712 | + |
| 713 | + def apply(self, p_array, dz): |
| 714 | + """ |
| 715 | + Apply the IBS process to a particle array over a given path length. |
| 716 | +
|
| 717 | + :param p_array: ParticleArray |
| 718 | + The particle array to be processed. |
| 719 | + :param dz: float |
| 720 | + The path length over which the IBS process is applied. |
| 721 | +
|
| 722 | + Raises: |
| 723 | + - ValueError: If an invalid method is specified. |
| 724 | + """ |
| 725 | + # Number of particles |
| 726 | + Nb = np.sum(p_array.q_array) / q_e |
| 727 | + |
| 728 | + # Update beam parameters |
| 729 | + self.get_beam_params(p_array) |
| 730 | + |
| 731 | + # estimate C log |
| 732 | + if self.update_Clog: |
| 733 | + Clog = self.estimate_Clog() |
| 734 | + else: |
| 735 | + Clog = self.Clog |
| 736 | + |
| 737 | + # Particle energy and momentum |
| 738 | + gamma = p_array.E / m_e_GeV |
| 739 | + pc = np.sqrt(p_array.E**2 - m_e_GeV**2) |
| 740 | + |
| 741 | + # Calculate sigma_gamma based on the selected method |
| 742 | + if self.method == "Huang": |
| 743 | + sigma_gamma = np.sqrt(Clog * ro_e**2 * Nb / (self.sigma_xy * self.emit_n * self.sigma_z) * dz / 4) |
| 744 | + |
| 745 | + elif self.method == "Nagaitsev": |
| 746 | + xi = (self.sigma_dgamma0 * self.sigma_xy / (gamma * self.emit_n))**2 |
| 747 | + F = (1 - xi**0.25) * np.log(xi + 1) / xi |
| 748 | + sigma_gamma = np.sqrt(Clog / 4 * ro_e**2 * Nb / (self.sigma_xy * self.emit_n * self.sigma_z) * F * dz) |
| 749 | + else: |
| 750 | + raise ValueError(f"Invalid method '{self.method}'. Choose 'Huang' or 'Nagaitsev'.") |
| 751 | + |
| 752 | + # Energy spread and update particle momenta |
| 753 | + sigma_e = sigma_gamma * m_e_GeV / pc |
| 754 | + p_array.p()[:] += np.random.randn(p_array.n) * sigma_e |
| 755 | + |
0 commit comments