# 3D rigid body energy level - from Matlab to Julia

**URL:** <https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880>\
**Category:** General Usage\
**Tags:** plotting\
**Created:** [January 17, 2025, 3:28pm UTC](https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880 "2025-01-17T15:28:08Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [January 17, 2025, 3:28pm UTC](https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880/1 "2025-01-17T15:28:09Z")

</div>

I would like to plot the energy level of the 3D rigid body, with Matlab found the following code

```julia
clear all
figure
for H=[.01:0.05:0.46 .5 .54:0.05:1.5]
I1=2;I2=1;I3=2/3;
y1=[-1:.0001:0];
y2=sqrt(((2*H-1/I3)-(1/I1-1/I3)*y1.^2)/(1/I2-1/I3));
y3=sqrt(((2*H-1/I2)-(1/I1-1/I2)*y1.^2)/(1/I3-1/I2));

%extract the real solutions y2 and y3:
m=1;
for k=1:length(y1)
    if( imag(y2(k))==0 & imag(y3(k))==0 )
      y1_new(m) = y1(k);
      y2_new(m) = y2(k);
      y3_new(m) = y3(k);
      m=m+1;
    end
end

hold on
if(m>1)
 plot3(y1_new,y2_new,y3_new)
 plot3(y1_new,-y2_new,y3_new)
 plot3(y1_new,y2_new,-y3_new)
 plot3(y1_new,-y2_new,-y3_new)

 plot3(-y1_new,y2_new,y3_new)
 plot3(-y1_new,-y2_new,y3_new)
 plot3(-y1_new,y2_new,-y3_new)
 plot3(-y1_new,-y2_new,-y3_new)
end
end

```

However, with Julia I tried the following

```julia
using Plots
I1=2;I2=1;I3=2/3;
for H in vcat(.01:0.05:0.46, .5, .54:0.05:1.5)
y1=-1:.0001:0
y2=sqrt.(Complex.(((2 .*H-1/I3).-(1/I1-1/I3) *y1 .^2)/(1/I2-1/I3)));
y3=sqrt.(Complex.(((2 .*H-1/I2).-(1/I1-1/I2) *y1 .^2)/(1/I3-1/I2)));

m=1;
y1_new = [];
y2_new = [];
y3_new = [];
for k = 1:length(y1)
    if( imag(y2[k])==0 && imag(y3[k])==0 )
      push!(y1_new, y1[k]);
      push!(y2_new, y2[k]);
      push!(y3_new, y3[k]);
      m=m+1;
    end
end
if(m>1)
plot3d(y1_new,y2_new,y3_new)
plot3d(y1_new,-y2_new,y3_new)
plot3d(y1_new,y2_new,-y3_new)
plot3d(y1_new,-y2_new,-y3_new)

plot3d(-y1_new,y2_new,y3_new)
plot3d(-y1_new,-y2_new,y3_new)
plot3d(-y1_new,y2_new,-y3_new)
plot3d(-y1_new,-y2_new,-y3_new)
end
end

```

Using Matlab code, the plot looks like

 ![Screenshot from 2025-01-17 16-25-22](https://global.discourse-cdn.com/julialang/original/3X/b/b/bbf99201ecef1571b878c3d351d82b729df76d50.png)

which I didn’t get with Julia.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 17, 2025, 4:28pm UTC](https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880/2 "2025-01-17T16:28:16Z")

</div>

You want the following (see the comments for explanation)

```julia
using Plots
I1=2;I2=1;I3=2/3;
for H in vcat(.01:0.05:0.46, .5, .54:0.05:1.5)
y1=-1:.0001:0
y2=sqrt.(Complex.(((2 .*H-1/I3).-(1/I1-1/I3) *y1 .^2)/(1/I2-1/I3)));
y3=sqrt.(Complex.(((2 .*H-1/I2).-(1/I1-1/I2) *y1 .^2)/(1/I3-1/I2)));

m=1;
plot() # clear the plot
y1_new = Float64[]; #make these Float64 because otherwise y2 and y3 are complex, so we will be pushing complex numbers (with 0 imaginary part) rather than real numbers.
y2_new = Float64[];
y3_new = Float64[];
for k = 1:length(y1)
    if( imag(y2[k])==0 && imag(y3[k])==0 )
      push!(y1_new, y1[k]);
      push!(y2_new, y2[k]);
      push!(y3_new, y3[k]);
      m=m+1;
    end
end
if(m>1)
# use plot3d! rather than plot3d to update rather than replace the plot
plot3d!(y1_new,y2_new,y3_new)
plot3d!(y1_new,-y2_new,y3_new)
plot3d!(y1_new,y2_new,-y3_new)
plot3d!(y1_new,-y2_new,-y3_new)

plot3d!(-y1_new,y2_new,y3_new)
plot3d!(-y1_new,-y2_new,y3_new)
plot3d!(-y1_new,y2_new,-y3_new)
plot3d!(-y1_new,-y2_new,-y3_new)
end
end
plot!() # show the final plot

```

---

<div class="post-metadata">

**Author:** ![p\_f](https://avatars.discourse-cdn.com/v4/letter/p/45deac/32.png) [@p\_f](https://discourse.julialang.org/u/p_f)\
**Post date:** [January 17, 2025, 6:03pm UTC](https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880/3 "2025-01-17T18:03:29Z")

</div>

With Oscar’s solution, I think the first `plot()` should be moved outside the first for loop.

I came from MATLAB, and personally much prefer Makie.jl over Plots.jl.

Using GLMakie.jl, you get a nice window you can rotate around a bit like in MATLAB. And the layouting features are very good.

Code:

```julia
using GLMakie

I1=2;I2=1;I3=2/3;

fig = Figure() #create a figure
ax = Axis3(fig[1,1]) #put 3d axis into fig. 

for H in vcat(.01:0.05:0.46, .5, .54:0.05:1.5)
  y1=-1:.0001:0
  y2=sqrt.(Complex.(((2 .*H-1/I3).-(1/I1-1/I3) *y1 .^2)/(1/I2-1/I3)));
  y3=sqrt.(Complex.(((2 .*H-1/I2).-(1/I1-1/I2) *y1 .^2)/(1/I3-1/I2)));
  
  m=1;
  
  y1_new = Float64[]; 
  y2_new = Float64[];
  y3_new = Float64[];
  for k = 1:length(y1)
    if( imag(y2[k])==0 && imag(y3[k])==0 )
      push!(y1_new, y1[k]);
      push!(y2_new, y2[k]);
      push!(y3_new, y3[k]);
      m=m+1;
    end
  end
  if m>1
    lines!(ax,y1_new,y2_new,y3_new)
    lines!(ax,y1_new,-y2_new,y3_new)
    lines!(ax,y1_new,y2_new,-y3_new)
    lines!(ax,y1_new,-y2_new,-y3_new)
    
    lines!(ax,-y1_new,y2_new,y3_new)
    lines!(ax,-y1_new,-y2_new,y3_new)
    lines!(ax,-y1_new,y2_new,-y3_new)
    lines!(ax,-y1_new,-y2_new,-y3_new)
  end
end

fig

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/a/aa796f6e53bdaee7d8018ef2d503c2f45ed936e8.png)

---

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [January 20, 2025, 9:14am UTC](https://discourse.julialang.org/t/3d-rigid-body-energy-level-from-matlab-to-julia/124880/4 "2025-01-20T09:14:12Z")

</div>

Dear @Oscar_Smith , thanks for you update but as @p_f said, `plot()` should be moved outside the first for loop!! Otherwise, it give you an empty figure!