From 8e0d3ba48d506d2c366febe6064ae7a37da3f767 Mon Sep 17 00:00:00 2001 From: okayamayu8 Date: Mon, 8 Jul 2024 13:37:05 +0900 Subject: [PATCH 1/2] change beam_cut file --- .../.ipynb_checkpoints/beam_cut-checkpoint.py | 191 ++++++++++++++++++ grasp2alm/beam_cut.py | 4 +- 2 files changed, 193 insertions(+), 2 deletions(-) create mode 100644 grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py diff --git a/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py b/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py new file mode 100644 index 0000000..5f0992b --- /dev/null +++ b/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py @@ -0,0 +1,191 @@ +# -*- encoding: utf-8 -*- + +from dataclasses import dataclass, field +import numpy as np +import matplotlib.pyplot as plt +from typing import Union +from .beam_polar import BeamPolar + +@dataclass +class BeamCut: + """Class to hold the data from a beam cut file of GRASP. + + Attributes: + header (str): Record with identification text. + vini (float): Initial value. + vinc (float): Increment. + vnum (int): Number of values in cut. + c (np.ndarray): Constant. + icomp (int): Polarization control parameter. + icut (int): Control parameter of cut. + ncomp (int): Number of field components. + ncut (int): Number of cuts. + amp (np.ndarray): Amplitude. + + Methods: + __init__(self, filepath): Initializes a BeamCut object. + __post_init__(self): Performs post-initialization tasks. + to_polar(self, copol_axis="x"): Converts the beam + to polar coordinates. + plot(self, pol='co', color_resol=20, figsize=6, cmap="inferno", + return_fields=False): Plots the beam. + + """ + header: str = "" + vini: float = 0.0 + vinc: float = 0.0 + vnum: int = 0 + c: np.ndarray = None # TODO check is it needed variable + + icomp: int = 0 + icut: int = 0 + ncomp: int = 0 + ncut: int = 0 + amp: np.ndarray = None + + def __init__(self, filepath): + super().__init__() + self.filepath = filepath + self.filename = filepath.split("/")[-1] + self.__post_init__() + + def __post_init__(self): + if not self.filepath.endswith(".cut"): + raise ValueError("Error in BeamCut.__post_init__: The file is not a GRASP cut file.") + with open(self.filepath, "r") as fi: + self.header = fi.readline().strip() + while True: + line = fi.readline() + if not line: + break + data = line.split() + if len(data) == 7: + self.vini, self.vinc, self.vnum, c, self.icomp, self.icut, self.ncomp = map(float, data) + self.vnum, self.icomp, self.icut, self.ncomp = map(int, (self.vnum, self.icomp, self.icut, self.ncomp)) + self.c = np.append(self.c, c) + self.ncut += 1 + if self.ncomp > 2: + raise ValueError("Three field components present. Beam package can only handle two field components.") + if self.vnum % 2 == 0: + raise ValueError("The number of pixels in a cut must be odd.") + + self.amp = np.zeros((self.ncomp, self.vnum, self.ncut), dtype=complex) + fi.seek(0) + cnt = 0 + while True: + line = fi.readline() + if not line: + break + data = line.split() + if len(data) == 7: + self.vini, self.vinc, self.vnum, self.c, self.icomp, self.icut, self.ncomp = map(float, data) + self.vnum, self.icomp, self.icut, self.ncomp = map(int, (self.vnum, self.icomp, self.icut, self.ncomp)) + for i in range(self.vnum): + line = fi.readline() + data = line.split() + tmp1, tmp2, tmp3, tmp4 = map(float, data) + self.amp[0, i, cnt] = complex(tmp1, tmp2) + self.amp[1, i, cnt] = complex(tmp3, tmp4) + cnt += 1 + + def to_polar(self, copol_axis="x"): + """Converts beam in "cut" format to Stokes parameters + on a polar grid. Assumes that cuts are evenly spaced + in theta. The value of copol specifies the alignment + of the co-polar basis ('x' or 'y') of the input GRASP file. + + Args: + copol_axis (str): The axis of copolarization. Must be either 'x' or 'y'. + + Returns: + BeamPolar: The beam in polar coordinates. + + Raises: + ValueError: If the beam is not in the expected format. + + """ + copol_axis = copol_axis.lower() + + if self.icomp != 3: + raise ValueError("Error in BeamCut.to_polar: beam is not in linear 'co' and 'cx' components") + if self.icut != 1: + raise ValueError("Error in BeamCut.to_polar: beam is not in phi cuts") + if self.ncomp != 2: + raise ValueError("Error in BeamCut.to_polar: beam has the wrong number of components") + assert copol_axis in ["x", "y"], "Error in BeamCut.to_polar: copol_axis must be 'x' or 'y'" + + nphi = int(2 * self.ncut) + ntheta = int(self.vnum // 2)+1 + theta_rad_min = 0.0 + theta_rad_max = np.deg2rad(np.abs(self.vini)) + beam_polar = BeamPolar(nphi, ntheta, theta_rad_min, theta_rad_max, self.filename) + amp_tmp = np.zeros((2, nphi, ntheta), dtype=complex) + + for i in range(self.ncut): + amp_tmp[:, i, :] = self.amp[:, ntheta-1:self.vnum+1, i] + amp_tmp[:, self.ncut + i, :] = self.amp[:, ntheta-1::-1, i] + + if copol_axis == "x": + sign = -1 + elif copol_axis == "y": + sign = 1 + + c = amp_tmp[0, :, :] + x = amp_tmp[1, :, :] + + modc2 = np.abs(c)**2 + modx2 = np.abs(x)**2 + acaxs = c * np.conj(x) + + beam_polar.stokes[0, :, :] = modc2 + modx2 + beam_polar.stokes[1, :, :] = sign * (modc2 - modx2) + beam_polar.stokes[2, :, :] = sign * 2.0 * np.real(acaxs) + beam_polar.stokes[3, :, :] = 2.0 * np.imag(acaxs) + return beam_polar + + def plot(self, pol='co', color_resol=20, figsize=6, cmap="inferno", return_fields=False): + """Plot the beam pattern. + + Args: + pol (str): The polarization to plot. Must be either 'co' or 'cx'. + color_resol (int): The number of color levels in the plot. + figsize (int): The size of the figure. + cmap (str): The colormap to use for the plot. + return_fields (bool): Whether to return the x, y, and z values. + + Returns: + If return_fields is False: + None + + If return_fields is True: + x (ndarray): The x values of the plot. + y (ndarray): The y values of the plot. + z (ndarray): The z values of the plot. + + """ + assert pol == 'co' or pol == 'cx', "Error in BeamCut.plot: pol must be 'co' or 'cx'" + + theta = np.deg2rad(np.linspace(self.vini, self.vini + self.vinc * (self.vnum - 1), self.vnum)) + phi = np.deg2rad(np.linspace(0.0, 180.0, self.ncut)) + + theta_grid, phi_grid = np.meshgrid(theta, phi) + if pol == 'co': + z = np.real(self.amp[0])**2 + np.imag(self.amp[0])**2 + dBz = 10*np.log10(z) + elif pol == 'cx': + z = np.real(self.amp[1])**2 + np.imag(self.amp[1])**2 + dBz = 10*np.log10(z) + x = np.rad2deg(theta_grid * np.cos(phi_grid)) + y = np.rad2deg(theta_grid * np.sin(phi_grid)) + levels = np.linspace(np.min(dBz), np.max(dBz), color_resol) + + if not return_fields: + plt.figure(figsize=(figsize,figsize)) + plt.axes().set_aspect("equal") + plt.title(f"{self.filename} : {pol}") + plt.xlabel("Degrees") + plt.ylabel("Degrees") + plt.contourf(x, y, dBz.T, levels=levels, cmap=cmap, extend='both') + plt.colorbar(orientation="vertical", label="dBi") + else: + return (x, y, z.T) diff --git a/grasp2alm/beam_cut.py b/grasp2alm/beam_cut.py index ce5da4b..5f0992b 100644 --- a/grasp2alm/beam_cut.py +++ b/grasp2alm/beam_cut.py @@ -115,14 +115,14 @@ def to_polar(self, copol_axis="x"): assert copol_axis in ["x", "y"], "Error in BeamCut.to_polar: copol_axis must be 'x' or 'y'" nphi = int(2 * self.ncut) - ntheta = int(self.vnum // 2) + ntheta = int(self.vnum // 2)+1 theta_rad_min = 0.0 theta_rad_max = np.deg2rad(np.abs(self.vini)) beam_polar = BeamPolar(nphi, ntheta, theta_rad_min, theta_rad_max, self.filename) amp_tmp = np.zeros((2, nphi, ntheta), dtype=complex) for i in range(self.ncut): - amp_tmp[:, i, :] = self.amp[:, ntheta:self.vnum-1, i] + amp_tmp[:, i, :] = self.amp[:, ntheta-1:self.vnum+1, i] amp_tmp[:, self.ncut + i, :] = self.amp[:, ntheta-1::-1, i] if copol_axis == "x": From 58f19e3953c2a0a5b2f8314d09bae01fbed8e527 Mon Sep 17 00:00:00 2001 From: okayamayu8 <49719415+okayamayu8@users.noreply.github.com> Date: Mon, 8 Jul 2024 16:46:17 +0900 Subject: [PATCH 2/2] Delete grasp2alm/.ipynb_checkpoints directory --- .../.ipynb_checkpoints/beam_cut-checkpoint.py | 191 ------------------ 1 file changed, 191 deletions(-) delete mode 100644 grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py diff --git a/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py b/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py deleted file mode 100644 index 5f0992b..0000000 --- a/grasp2alm/.ipynb_checkpoints/beam_cut-checkpoint.py +++ /dev/null @@ -1,191 +0,0 @@ -# -*- encoding: utf-8 -*- - -from dataclasses import dataclass, field -import numpy as np -import matplotlib.pyplot as plt -from typing import Union -from .beam_polar import BeamPolar - -@dataclass -class BeamCut: - """Class to hold the data from a beam cut file of GRASP. - - Attributes: - header (str): Record with identification text. - vini (float): Initial value. - vinc (float): Increment. - vnum (int): Number of values in cut. - c (np.ndarray): Constant. - icomp (int): Polarization control parameter. - icut (int): Control parameter of cut. - ncomp (int): Number of field components. - ncut (int): Number of cuts. - amp (np.ndarray): Amplitude. - - Methods: - __init__(self, filepath): Initializes a BeamCut object. - __post_init__(self): Performs post-initialization tasks. - to_polar(self, copol_axis="x"): Converts the beam - to polar coordinates. - plot(self, pol='co', color_resol=20, figsize=6, cmap="inferno", - return_fields=False): Plots the beam. - - """ - header: str = "" - vini: float = 0.0 - vinc: float = 0.0 - vnum: int = 0 - c: np.ndarray = None # TODO check is it needed variable - - icomp: int = 0 - icut: int = 0 - ncomp: int = 0 - ncut: int = 0 - amp: np.ndarray = None - - def __init__(self, filepath): - super().__init__() - self.filepath = filepath - self.filename = filepath.split("/")[-1] - self.__post_init__() - - def __post_init__(self): - if not self.filepath.endswith(".cut"): - raise ValueError("Error in BeamCut.__post_init__: The file is not a GRASP cut file.") - with open(self.filepath, "r") as fi: - self.header = fi.readline().strip() - while True: - line = fi.readline() - if not line: - break - data = line.split() - if len(data) == 7: - self.vini, self.vinc, self.vnum, c, self.icomp, self.icut, self.ncomp = map(float, data) - self.vnum, self.icomp, self.icut, self.ncomp = map(int, (self.vnum, self.icomp, self.icut, self.ncomp)) - self.c = np.append(self.c, c) - self.ncut += 1 - if self.ncomp > 2: - raise ValueError("Three field components present. Beam package can only handle two field components.") - if self.vnum % 2 == 0: - raise ValueError("The number of pixels in a cut must be odd.") - - self.amp = np.zeros((self.ncomp, self.vnum, self.ncut), dtype=complex) - fi.seek(0) - cnt = 0 - while True: - line = fi.readline() - if not line: - break - data = line.split() - if len(data) == 7: - self.vini, self.vinc, self.vnum, self.c, self.icomp, self.icut, self.ncomp = map(float, data) - self.vnum, self.icomp, self.icut, self.ncomp = map(int, (self.vnum, self.icomp, self.icut, self.ncomp)) - for i in range(self.vnum): - line = fi.readline() - data = line.split() - tmp1, tmp2, tmp3, tmp4 = map(float, data) - self.amp[0, i, cnt] = complex(tmp1, tmp2) - self.amp[1, i, cnt] = complex(tmp3, tmp4) - cnt += 1 - - def to_polar(self, copol_axis="x"): - """Converts beam in "cut" format to Stokes parameters - on a polar grid. Assumes that cuts are evenly spaced - in theta. The value of copol specifies the alignment - of the co-polar basis ('x' or 'y') of the input GRASP file. - - Args: - copol_axis (str): The axis of copolarization. Must be either 'x' or 'y'. - - Returns: - BeamPolar: The beam in polar coordinates. - - Raises: - ValueError: If the beam is not in the expected format. - - """ - copol_axis = copol_axis.lower() - - if self.icomp != 3: - raise ValueError("Error in BeamCut.to_polar: beam is not in linear 'co' and 'cx' components") - if self.icut != 1: - raise ValueError("Error in BeamCut.to_polar: beam is not in phi cuts") - if self.ncomp != 2: - raise ValueError("Error in BeamCut.to_polar: beam has the wrong number of components") - assert copol_axis in ["x", "y"], "Error in BeamCut.to_polar: copol_axis must be 'x' or 'y'" - - nphi = int(2 * self.ncut) - ntheta = int(self.vnum // 2)+1 - theta_rad_min = 0.0 - theta_rad_max = np.deg2rad(np.abs(self.vini)) - beam_polar = BeamPolar(nphi, ntheta, theta_rad_min, theta_rad_max, self.filename) - amp_tmp = np.zeros((2, nphi, ntheta), dtype=complex) - - for i in range(self.ncut): - amp_tmp[:, i, :] = self.amp[:, ntheta-1:self.vnum+1, i] - amp_tmp[:, self.ncut + i, :] = self.amp[:, ntheta-1::-1, i] - - if copol_axis == "x": - sign = -1 - elif copol_axis == "y": - sign = 1 - - c = amp_tmp[0, :, :] - x = amp_tmp[1, :, :] - - modc2 = np.abs(c)**2 - modx2 = np.abs(x)**2 - acaxs = c * np.conj(x) - - beam_polar.stokes[0, :, :] = modc2 + modx2 - beam_polar.stokes[1, :, :] = sign * (modc2 - modx2) - beam_polar.stokes[2, :, :] = sign * 2.0 * np.real(acaxs) - beam_polar.stokes[3, :, :] = 2.0 * np.imag(acaxs) - return beam_polar - - def plot(self, pol='co', color_resol=20, figsize=6, cmap="inferno", return_fields=False): - """Plot the beam pattern. - - Args: - pol (str): The polarization to plot. Must be either 'co' or 'cx'. - color_resol (int): The number of color levels in the plot. - figsize (int): The size of the figure. - cmap (str): The colormap to use for the plot. - return_fields (bool): Whether to return the x, y, and z values. - - Returns: - If return_fields is False: - None - - If return_fields is True: - x (ndarray): The x values of the plot. - y (ndarray): The y values of the plot. - z (ndarray): The z values of the plot. - - """ - assert pol == 'co' or pol == 'cx', "Error in BeamCut.plot: pol must be 'co' or 'cx'" - - theta = np.deg2rad(np.linspace(self.vini, self.vini + self.vinc * (self.vnum - 1), self.vnum)) - phi = np.deg2rad(np.linspace(0.0, 180.0, self.ncut)) - - theta_grid, phi_grid = np.meshgrid(theta, phi) - if pol == 'co': - z = np.real(self.amp[0])**2 + np.imag(self.amp[0])**2 - dBz = 10*np.log10(z) - elif pol == 'cx': - z = np.real(self.amp[1])**2 + np.imag(self.amp[1])**2 - dBz = 10*np.log10(z) - x = np.rad2deg(theta_grid * np.cos(phi_grid)) - y = np.rad2deg(theta_grid * np.sin(phi_grid)) - levels = np.linspace(np.min(dBz), np.max(dBz), color_resol) - - if not return_fields: - plt.figure(figsize=(figsize,figsize)) - plt.axes().set_aspect("equal") - plt.title(f"{self.filename} : {pol}") - plt.xlabel("Degrees") - plt.ylabel("Degrees") - plt.contourf(x, y, dBz.T, levels=levels, cmap=cmap, extend='both') - plt.colorbar(orientation="vertical", label="dBi") - else: - return (x, y, z.T)