Skip to content

equations

Functions:

  • eq1_sol –

    Solve the first equilibrium ODE.

  • eq2_sol –

    Solve the second equilibrium ODE.

  • ic_eq_1 –

    Initial conditions for the first equilibrium ODE.

  • ic_eq_2 –

    Initial conditions for the second equilibrium ODE.

  • rhs_eq_1 –

    Compute the right-hand side of the first equilibrium ODE.

  • rhs_eq_2 –

    Compute the right-hand side of the second equilibrium ODE.

  • sigma –

    Calculate the sigma value based on the given parameters.

eq1_sol 

eq1_sol(alpha: float, theta_0: float, delta_h: float, r0: float = RMIN / 2.0) -> OdeSolution

Solve the first equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
127
128
129
130
131
132
133
134
135
136
def eq1_sol(alpha: float, theta_0: float, delta_h: float, r0: float = RMIN / 2.0) -> OdeSolution:
    """
    Solve the first equilibrium ODE.
    """
    y0 = ic_eq_1(alpha, theta_0, delta_h, r0)
    a = RA - np.abs(delta_h)

    sol = solve_ivp(rhs_eq_1, (r0, a), y0, args=(alpha, theta_0, delta_h), method="LSODA", dense_output=True, rtol=RTOL, atol=ATOL)

    return sol.sol

eq2_sol 

eq2_sol(alpha: float, theta_0: float, delta_h: float, delta_a: float, r0: float = RMIN / 2.0) -> Tuple[OdeSolution, bool]

Solve the second equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def eq2_sol(alpha: float, theta_0: float, delta_h: float, delta_a: float, r0: float = RMIN / 2.0) -> Tuple[OdeSolution, bool]:
    """
    Solve the second equilibrium ODE.
    """
    y0 = ic_eq_2(alpha, theta_0, delta_h, delta_a, r0)

    sol = solve_ivp(rhs_eq_2, (r0, RMAX), y0, args=(alpha, theta_0, delta_h), method="LSODA", dense_output=True, rtol=RTOL, atol=ATOL)

    # compute q from the solution
    r = sol.t
    bph, bth, delta, ddelta, bph2, bth2 = sol.y
    q = (bph + bph2) / (bth + bth2)

    a = RA - np.abs(delta_h)
    rfp_flag = np.min(q[r < a]) < 0.0

    return sol.sol, rfp_flag

ic_eq_1 

ic_eq_1(alpha: float, theta_0: float, delta_h: float, r0: float | ndarray = RMIN / 2.0) -> List[float] | ndarray

Initial conditions for the first equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def ic_eq_1(alpha: float, theta_0: float, delta_h: float, r0: float | np.ndarray = RMIN / 2.0) -> List[float] | np.ndarray:
    """
    Initial conditions for the first equilibrium ODE.
    """
    a = RA - np.abs(delta_h)

    bph_0 = r0
    bth_0 = 0.5 * 2.0 * theta_0 / a * R0 * r0
    delta_0 = -(r0**2) / (8.0 * R0)
    ddelta_0 = -r0 / (4.0 * R0)

    if isinstance(delta_h, np.ndarray) or isinstance(r0, np.ndarray):
        return np.concatenate([np.atleast_1d(bph_0), np.atleast_1d(bth_0), np.atleast_1d(delta_0), np.atleast_1d(ddelta_0)], axis=-1)
    else:
        return [bph_0, bth_0, delta_0, ddelta_0]

ic_eq_2 

ic_eq_2(alpha: float, theta_0: float, delta_h: float, delta_a: float, r0: float | ndarray = RMIN / 2.0) -> List[float] | ndarray

Initial conditions for the second equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def ic_eq_2(alpha: float, theta_0: float, delta_h: float, delta_a: float, r0: float | np.ndarray = RMIN / 2.0) -> List[float] | np.ndarray:
    """
    Initial conditions for the second equilibrium ODE.
    """
    a = RA - np.abs(delta_h)

    bph_0 = r0
    bth_0 = 0.5 * 2.0 * theta_0 / a * R0 * r0
    delta_0 = -(r0**2) / (8.0 * R0) - delta_a + delta_h
    ddelta_0 = -r0 / (4.0 * R0)
    bph2_0 = r0 * (3.0 / (2.0 * theta_0 / a * R0) ** 2 - delta_0 / R0)
    bth2_0 = 3.0 * r0 / (2.0 * 2.0 * theta_0 / a * R0)

    if isinstance(delta_h, np.ndarray) or isinstance(r0, np.ndarray):
        return np.concatenate([np.atleast_1d(bph_0), np.atleast_1d(bth_0), np.atleast_1d(delta_0), np.atleast_1d(ddelta_0), np.atleast_1d(bph2_0), np.atleast_1d(bth2_0)], axis=-1)
    else:
        return [bph_0, bth_0, delta_0, ddelta_0, bph2_0, bth2_0]

rhs_eq_1 

rhs_eq_1(r: float | ndarray, y: List[float] | Tuple[float] | ndarray, alpha: float, theta_0: float, delta_h: float) -> Tuple[float | ndarray]

Compute the right-hand side of the first equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def rhs_eq_1(r: float | np.ndarray, y: List[float] | Tuple[float] | np.ndarray, alpha: float, theta_0: float, delta_h: float) -> Tuple[float | np.ndarray]:
    """
    Compute the right-hand side of the first equilibrium ODE.
    """
    # if y is a 2D array extract each variable keeping the dimensions
    if isinstance(y, np.ndarray) and y.ndim == 2:
        bph, bth, delta, ddelta = y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:4]
    else:
        bph, bth, delta, ddelta = y[0], y[1], y[2], y[3]

    sigma_r = sigma(r, alpha, theta_0, delta_h)

    dbph_dr = bph / r - sigma_r * r / R0 * bth
    dbth_dr = -bth / r + sigma_r * R0 / r * bph
    ddelta_dr = ddelta
    ddelta_dr2 = -ddelta / r * (1.0 + 2.0 * r * dbth_dr / bth) - 1.0 / R0

    if isinstance(y, np.ndarray) and y.ndim == 2:
        return dbph_dr, dbth_dr, ddelta_dr, ddelta_dr2
    else:
        return [dbph_dr, dbth_dr, ddelta_dr, ddelta_dr2]

rhs_eq_2 

rhs_eq_2(r: float | ndarray, y: List[float] | Tuple[float] | ndarray, alpha: float, theta_0: float, delta_h: float) -> Tuple[float | ndarray]

Compute the right-hand side of the second equilibrium ODE.

Source code in src/fpga_profile_reco/data/equations.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
def rhs_eq_2(r: float | np.ndarray, y: List[float] | Tuple[float] | np.ndarray, alpha: float, theta_0: float, delta_h: float) -> Tuple[float | np.ndarray]:
    """
    Compute the right-hand side of the second equilibrium ODE.
    """
    # if y is a 2D array extract each variable keeping the dimensions
    if isinstance(y, np.ndarray) and y.ndim == 2:
        bph, bth, delta, ddelta, bph2, bth2 = (y[:, 0:1], y[:, 1:2], y[:, 2:3], y[:, 3:4], y[:, 4:5], y[:, 5:6])
    else:
        bph, bth, delta, ddelta, bph2, bth2 = y[0], y[1], y[2], y[3], y[4], y[5]

    sigma_r = sigma(r, alpha, theta_0, delta_h)

    dbph_dr = bph / r - sigma_r * r / R0 * bth
    dbth_dr = -bth / r + sigma_r * R0 / r * bph
    ddelta_dr = ddelta
    ddelta_dr2 = -ddelta / r * (1.0 + 2.0 * r * dbth_dr / bth) - 1.0 / R0
    dbph2_dr = bph2 / r - r / R0 * ((delta / R0 + r / (2.0 * R0) * ddelta - r**2 / (2.0 * R0**2)) * (R0 / r * dbph_dr - R0 / r**2 * bph) + R0 / r * bph * (1.5 * ddelta / R0 + r * ddelta_dr2 / (2.0 * R0) - r / R0**2)) - sigma_r * r / R0 * bth2
    dbth2_dr = -bth2 / r - R0 / r * ((r**2 / (2.0 * R0**2) + ddelta**2 / 2.0 + r * ddelta / (2.0 * R0) - delta / R0) * (bth / R0 + r / R0 * dbth_dr)) - bth * (r / R0**2 + ddelta * ddelta_dr2 - ddelta / (2.0 * R0) + r / (2.0 * R0) * ddelta_dr2) + sigma_r * R0 / r * bph2

    if isinstance(y, np.ndarray) and y.ndim == 2:
        return dbph_dr, dbth_dr, ddelta_dr, ddelta_dr2, dbph2_dr, dbth2_dr
    else:
        return [dbph_dr, dbth_dr, ddelta_dr, ddelta_dr2, dbph2_dr, dbth2_dr]

sigma 

sigma(r: float | ndarray, alpha: float, theta_0: float, delta_h: float) -> float | ndarray

Calculate the sigma value based on the given parameters.

Source code in src/fpga_profile_reco/data/equations.py 
24
25
26
27
28
29
30
def sigma(r: float | np.ndarray, alpha: float, theta_0: float, delta_h: float) -> float | np.ndarray:
    """
    Calculate the sigma value based on the given parameters.
    """
    a = RA - np.abs(delta_h)
    sigma_r = np.where(r <= a, (2.0 * theta_0 / a) * (1.0 - (r / a) ** alpha), 0.0)
    return sigma_r