## Issue

I’m facing a mathematical/programming problem which I don’t really understand. In the following `α,β,θc`

are given parameters, `ϕ`

is the variable. I have to solve the equation `Cp=0`

, where:

and:

The solutions that I have computed are:

Coding the above solutions with Python/Numpy, I find that:

- If
`β=0`

, then both solutions are correct:`Cp(ϕ1)=Cp(ϕ2)=0`

. - However, if
`β!=0`

then`Cp(ϕ1)!=0`

(which is wrong) and`Cp(ϕ2)=0`

. Why? I checked again and again the signs and terms on my coded expressions, they seem fine to me…

Here is the code:

```
import numpy as np
theta_c = np.deg2rad(10)
alpha = np.deg2rad([1, 2, 4, 6, 8, 10, 12, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85])
beta = np.deg2rad(15)
lamb = np.cos(alpha) * np.cos(beta)
nu = -np.sin(beta)
omega = np.sin(alpha) * np.cos(beta)
common_1 = lamb * nu * np.tan(theta_c)
common_2 = omega * np.sqrt(omega**2 + nu**2 - lamb**2 * np.tan(theta_c)**2)
common_3 = omega**2 + nu**2
phi_1 = np.arcsin((-common_1 + common_2) / common_3)
phi_2 = np.arcsin((-common_1 - common_2) / common_3)
Cp = lambda phi: lamb * np.sin(theta_c) + nu * np.cos(theta_c) * np.sin(phi) - omega * np.cos(phi) * np.cos(theta_c)
print(Cp(phi_1))
# OUTPUT: this is wrong!
# [-0.02357418 -0.04432021 -0.0779439 -0.10278244 -0.12111675 -0.13494261
# -0.1457658 -0.15858164 -0.17562354 -0.19112922 -0.20703468 -0.22411243
# -0.24265822 -0.26274848 -0.28434621 -0.3073487 -0.33161092 -0.35695785
# -0.38319138 -0.41009444 -0.43743318 -0.4649581 ]
print(Cp(phi_2))
# OUTPUT: all values are very close to zero. Good solution.
# [-1.04083409e-17 -4.16333634e-17 0.00000000e+00 -1.38777878e-17
# -2.77555756e-17 -2.77555756e-17 0.00000000e+00 0.00000000e+00
# 5.55111512e-17 -1.11022302e-16 1.11022302e-16 5.55111512e-17
# -1.11022302e-16 0.00000000e+00 -3.88578059e-16 -1.11022302e-16
# -2.77555756e-16 1.66533454e-16 5.55111512e-17 2.22044605e-16
# 5.55111512e-17 -2.22044605e-16]
```

Alternatively, I could use `fsolve`

from Scipy optimize, which I have already verified it computes the correct solutions. However, this procedure will be executed many times by another iterative algorithm, so I’d rather use the analytical solutions.

## Solution

Seems like your analytical solution has a mistake.

I rewrote the target function in form

```
Cp(ϕ) = 2 * (A + B sin ϕ - C cos ϕ)^2
```

Solutions of `Cp = 0`

come from `A + B sin ϕ - C cos ϕ = 0`

equation, which I solved using Weierstrass substitution (WS). The solutions `ϕ₁`

, `ϕ₂`

I found are

```
# t is parameter in WS
-B ± sqrt(B^2 - (A^2 - C^2))
t₁ and t₂ = ----------------------------
A + C
ϕ₁ = 2atan(t₁)
ϕ₂ = 2atan(t₂)
```

Here is code snippet written in Julia

```
using Plots
θ = deg2rad(10)
α = deg2rad(10)
β = deg2rad(15)
λ = cos(α) * cos(β)
ν = -sin(β)
ω = sin(α) * cos(β)
# Cp = 2 * (A + B sin ϕ - C cos ϕ)^2
A = λ * sin(θ)
B = ν * cos(θ)
C = ω * cos(θ)
t₁ = (-B + sqrt(B^2 - (A^2 - C^2)))/(A + C)
t₂ = (-B - sqrt(B^2 - (A^2 - C^2)))/(A + C)
ϕ₁ = 2atan(t₁)
ϕ₂ = 2atan(t₂)
Cfunc(ϕ) = 2*(A + B * sin(ϕ) - C * cos(ϕ))^2
plot(Cfunc; xlims=(-10, 10), label="Cₚ")
scatter!([ϕ₁, ϕ₂], Cfunc.([ϕ₁, ϕ₂]); label="Solutions")
```

producing

I played with value of `β`

, seems like the code above works good for both zero and non-zero degrees. Also, I didn’t care about function domains, you should reproduce the derivation to ensure that all parameters and solution are well-defined.

Answered By – Stepan Zakharov

Answer Checked By – Pedro (BugsFixing Volunteer)