3D rigid body energy level - from Matlab to Julia

General Usage
1

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

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

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
Screenshot from 2025-01-17 16-25-22624×566 8.96 KB

which I didn’t get with Julia.

2

You want the following (see the comments for explanation)

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
2 Likes
3

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:

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
image623×501 77.6 KB

3 Likes