<p>I tried to plot the paths of  particle of different sizes in the same polar coordinate using the code given below using matplotlib. somewhat like this:<a href="https://i.sstatic.net/jkfQm.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jkfQm.jpg" alt="path of different particle sizes" /></a></p>
<p>Instead I got this plot:
<a href="https://i.sstatic.net/vRuKj.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/vRuKj.png" alt="plot for only one value of particle size 1.0e-11m is visibe" /></a>while the path for first four particle sizes are invisible as the curves are overlapping.This means the curves are now following the same line. <a href="https://i.sstatic.net/Qy5N2.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Qy5N2.png" alt="same paths curve" /></a>. if the simulation is run separately for different particle sizes they follow different paths as shown: <a href="https://i.sstatic.net/Vnrtb.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Vnrtb.jpg" alt="different path for different particles" /></a></p>
<p>if the particle size 1.0 e-10m the path of the 2.3e-11m particle size alone is visible as below:<a href="https://i.sstatic.net/nWzVm.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/nWzVm.png" alt="plot showing the path of only 2.3e-11m particle" /></a>. The python code used in the simulations is as given below:</p>
<pre><code>from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

# Constants
μ = 379312077e8
R = 60268e3
g_10 = 21141e-9
Ω = 1.7e-4
ρ = 1e3
b_values = [0.95e-10, 0.94e-10, 0.26e-10, 0.25e-10, 0.23e-10, 0.1e-10]
V = 10
ε = 8.85e-12
j_2 = 1.629071e-2

# Function to define the ODEs
def odes(t, p,b):
    r, x, θ, y, ϕ, z = p

# Constants for the given b
    m = (4/3) * pi * (b**3) * ρ
    q = 4 * pi * ε * b * V
    β = q / m

# Other constants
    μ_0 = 4 * pi * 1e-7
    B_θ = (μ_0 / (4 * pi)) * (R / r)**3 * g_10 * sin(θ)
    B_r = (μ_0 / (4 * pi)) * 2 * (R / r)**3 * g_10 * cos(θ)

# Defining the ODEs
    drdt = x
    dxdt = r * (y**2 + (z + Ω)**2 * sin(θ)**2 - β * z * sin(θ) * B_θ) - (μ / r**2) * (1 - (3/2) * j_2 * (R / r)**2 * (3 * cos(θ)**2 - 1))
    dθdt = y
    dydt = (-2 * x * y + r * (z + Ω)**2 * sin(θ) * cos(θ) + (β * r * z * sin(θ) * B_r)) / r + (3 * μ / r**2) * j_2 * (R / r)**2 * sin(θ) * cos(θ)
    dϕdt = z
    dzdt = (-2 * (z + Ω) * (x * sin(θ) + r * y * cos(θ)) + (β * (x * B_θ - r * y * B_r)) / (r * sin(θ)))

    return np.array([drdt, dxdt, dθdt, dydt, dϕdt, dzdt])
 # Define the event function
def event_func(t, p,b):
    return p[0] - R

event_func.terminal = True
event_func.direction = -1

# Time window
t_span = (0, 1000)
t = np.linspace(t_span[0], t_span[1], 1001)

# Initial conditions
r0 = 1.12 * R
p0 = np.array([1.12 * R, 0.0, 90.0 * (pi / 180), 0.0, 0.0, 0.0203 * (pi / 180)])

# Create a figure with a polar projection
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})

# Define colors for plotting
colors = ['r', 'g', 'b', 'm', 'y', 'k']

# Iterate through different values of b
for b, color in zip(b_values, colors):
    
# Solve IVP for the current value of b
    sol = solve_ivp(odes, t_span, p0, args=(b,), t_eval=t,events=(event_func,), method=&quot;DOP853&quot;,
                dense_output=True, vectorized=False, atol=1e-4, rtol=1e-5)

    # Extract r and ϕ values
    r = sol.y[0, :]
    ϕ = sol.y[4, :]

    # Plot r against ϕ for the current value of b
    ax.plot(ϕ, r, label=f'b = {b:.2e}', color=color)

# Set radial limits
ax.set_rmax(1.5 * R)

# Set radial ticks
ax.set_rticks([R, 1.12 * R, 1.18 * R, 1.2 * R, 1.3 * R])

# Add a circle at the center
circle = plt.Circle((0, 0), R, transform=ax.transData._b, facecolor='yellow',
                fill=True, edgecolor='red', linewidth=0.5, zorder=2, alpha=0.5)
ax.add_artist(circle)

# Set plot title and legend
ax.set_title(&quot;Polar Plot of r vs Azimuthal Angle (ϕ)&quot;)
ax.legend()
plt.show()
</code></pre>
<p>**Why on running the simulation for different sizes changes the paths to cirular ones?**if anyone could help me to get the desired plot I would be immensely grateful. Thanks</p>

CC BY-SA 4.0
2023-11-01T12:21:03.787Z
Hi Lunthang Peter, from your chosen tags I realize that you're trying to plot orbits, but since it's purely a Python question, I don't think it belongs here. You should probably try StackOverflow instead.
pela
CC BY-SA 4.0
2023-11-02T18:13:25.560Z
@pela. Thank you. I shall do as you suggested.
Lunthang Peter
<astrophysics><python><planetary-formation><n-body-simulations>
How to plot different curves in the same polar cordinates in matplotlib?
Lunthang Peter
# How to plot different curves in the same polar cordinates in matplotlib?

(Asked by: username_0 on 2023-10-30)

I tried to plot the paths of particle of different sizes in the same polar coordinate using the code given below using matplotlib. somewhat like this:[![path of different particle sizes](https://i.sstatic.net/jkfQm.jpg)](https://i.sstatic.net/jkfQm.jpg)

Instead I got this plot:
[![plot for only one value of particle size 1.0e-11m is visibe](https://i.sstatic.net/vRuKj.png)](https://i.sstatic.net/vRuKj.png)while the path for first four particle sizes are invisible as the curves are overlapping.This means the curves are now following the same line. [![same paths curve](https://i.sstatic.net/Qy5N2.png)](https://i.sstatic.net/Qy5N2.png). if the simulation is run separately for different particle sizes they follow different paths as shown: [![different path for different particles](https://i.sstatic.net/Vnrtb.jpg)](https://i.sstatic.net/Vnrtb.jpg)

if the particle size 1.0 e-10m the path of the 2.3e-11m particle size alone is visible as below:[![plot showing the path of only 2.3e-11m particle](https://i.sstatic.net/nWzVm.png)](https://i.sstatic.net/nWzVm.png). The python code used in the simulations is as given below:

```
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

# Constants
μ = 379312077e8
R = 60268e3
g_10 = 21141e-9
Ω = 1.7e-4
ρ = 1e3
b_values = [0.95e-10, 0.94e-10, 0.26e-10, 0.25e-10, 0.23e-10, 0.1e-10]
V = 10
ε = 8.85e-12
j_2 = 1.629071e-2

# Function to define the ODEs
def odes(t, p,b):
    r, x, θ, y, ϕ, z = p

# Constants for the given b
    m = (4/3) * pi * (b**3) * ρ
    q = 4 * pi * ε * b * V
    β = q / m

# Other constants
    μ_0 = 4 * pi * 1e-7
    B_θ = (μ_0 / (4 * pi)) * (R / r)**3 * g_10 * sin(θ)
    B_r = (μ_0 / (4 * pi)) * 2 * (R / r)**3 * g_10 * cos(θ)

# Defining the ODEs
    drdt = x
    dxdt = r * (y**2 + (z + Ω)**2 * sin(θ)**2 - β * z * sin(θ) * B_θ) - (μ / r**2) * (1 - (3/2) * j_2 * (R / r)**2 * (3 * cos(θ)**2 - 1))
    dθdt = y
    dydt = (-2 * x * y + r * (z + Ω)**2 * sin(θ) * cos(θ) + (β * r * z * sin(θ) * B_r)) / r + (3 * μ / r**2) * j_2 * (R / r)**2 * sin(θ) * cos(θ)
    dϕdt = z
    dzdt = (-2 * (z + Ω) * (x * sin(θ) + r * y * cos(θ)) + (β * (x * B_θ - r * y * B_r)) / (r * sin(θ)))

    return np.array([drdt, dxdt, dθdt, dydt, dϕdt, dzdt])
 # Define the event function
def event_func(t, p,b):
    return p[0] - R

event_func.terminal = True
event_func.direction = -1

# Time window
t_span = (0, 1000)
t = np.linspace(t_span[0], t_span[1], 1001)

# Initial conditions
r0 = 1.12 * R
p0 = np.array([1.12 * R, 0.0, 90.0 * (pi / 180), 0.0, 0.0, 0.0203 * (pi / 180)])

# Create a figure with a polar projection
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})

# Define colors for plotting
colors = ['r', 'g', 'b', 'm', 'y', 'k']

# Iterate through different values of b
for b, color in zip(b_values, colors):
    
# Solve IVP for the current value of b
    sol = solve_ivp(odes, t_span, p0, args=(b,), t_eval=t,events=(event_func,), method="DOP853",
                dense_output=True, vectorized=False, atol=1e-4, rtol=1e-5)

    # Extract r and ϕ values
    r = sol.y[0, :]
    ϕ = sol.y[4, :]

    # Plot r against ϕ for the current value of b
    ax.plot(ϕ, r, label=f'b = {b:.2e}', color=color)

# Set radial limits
ax.set_rmax(1.5 * R)

# Set radial ticks
ax.set_rticks([R, 1.12 * R, 1.18 * R, 1.2 * R, 1.3 * R])

# Add a circle at the center
circle = plt.Circle((0, 0), R, transform=ax.transData._b, facecolor='yellow',
                fill=True, edgecolor='red', linewidth=0.5, zorder=2, alpha=0.5)
ax.add_artist(circle)

# Set plot title and legend
ax.set_title("Polar Plot of r vs Azimuthal Angle (ϕ)")
ax.legend()
plt.show()

```

\*\*Why on running the simulation for different sizes changes the paths to cirular ones?\*\*if anyone could help me to get the desired plot I would be immensely grateful. Thanks

**Tags:** astrophysics, python, planetary-formation, n-body-simulations

> **Comments:**
> **username_1**: Hi username_0 username_0, from your chosen tags I realize that you're trying to plot orbits, but since it's purely a Python question, I don't think it belongs here. You should probably try StackOverflow instead.  
> **username_0**: @username_1. Thank you. I shall do as you suggested.  

---

## 0 answers

