I am looking at plotting the the slowness profiles of Titinium (Inverser of velocity) by solving the christoffell equation.
I have managed to get a polar plot of this but would like to visualise this in 3d. For my polar plot I have done.
c11 = 162000000000;
c12=92000000000;
c44=47000000000;
c13= 69000000000;
c33= 181000000000;
c66= 35000000000
rho = 4430;
C=[c11,c12,c13,0,0,0;
c12,c11,c13,0,0,0;
c13,c13,c33,0,0,0;
0,0,0,c44,0,0;
0,0,0,0,c44,0;
0,0,0,0,0,c66];
phase_vel3=zeros(1,181);
r=0;
for theta=-pi:pi/180:pi
r=r+1;
n=[cos(theta),0,sin(theta)];
L_11=(C(1,1)*n(1)^2+C(6,6)*n(2)^2+C(5,5)*n(3)^2+2*C(1,6)*n(1)*n(2)+2*C(1,5)*n(1)*n(3)+2*C(5,6)*n(2)*n(3));
L_12=(C(1,6)*n(1)^2+C(2,6)*n(2)^2+C(4,5)*n(3)^2+(C(1,2)+C(6,6))*n(1)*n(2)+(C(1,4)+C(5,6))*n(1)*n(3)+(C(4,6)+C(2,5))*n(2)*n(3));
L_13=(C(1,5)*n(1)^2+C(4,6)*n(2)^2+C(3,5)*n(3)^2+(C(1,4)+C(5,6))*n(1)*n(2)+(C(1,3)+C(5,5))*n(1)*n(3)+(C(3,6)+C(4,5))*n(2)*n(3));
L_22=(C(6,6)*n(1)^2+C(2,2)*n(2)^2+C(4,4)*n(3)^2+2*C(2,6)*n(1)*n(2)+2*C(4,6)*n(1)*n(3)+2*C(2,4)*n(2)*n(3));
L_23=(C(5,6)*n(1)^2+C(2,4)*n(2)^2+C(3,4)*n(3)^2+(C(4,6)+C(2,5))*n(1)*n(2)+(C(3,6)+C(4,5))*n(1)*n(3)+(C(2,3)+C(4,4))*n(2)*n(3));
L_33=(C(5,5)*n(1)^2+C(4,4)*n(2)^2+C(3,3)*n(3)^2+2*C(4,5)*n(1)*n(2)+2*C(3,5)*n(1)*n(3)+2*C(3,4)*n(2)*n(3));
Christoffel_mat=[L_11, L_12, L_13;
L_12,L_22,L_23;
L_13,L_23,L_33];
[ev,d]=eig(Christoffel_mat);
% eigvecs = polarisation vectors
pl=ev(:,3);
% eigvals ~ slowness (phase)
phase_vel3(r)=sqrt(rho./d(3,3)); % quasi-longitudinal
r=0;
for theta=-pi:pi/180:pi
r=r+1;
n=[cos(theta),0,sin(theta)];
L_11=(C(1,1)*n(1)^2+C(6,6)*n(2)^2+C(5,5)*n(3)^2+2*C(1,6)*n(1)*n(2)+2*C(1,5)*n(1)*n(3)+2*C(5,6)*n(2)*n(3));
L_12=(C(1,6)*n(1)^2+C(2,6)*n(2)^2+C(4,5)*n(3)^2+(C(1,2)+C(6,6))*n(1)*n(2)+(C(1,4)+C(5,6))*n(1)*n(3)+(C(4,6)+C(2,5))*n(2)*n(3));
L_13=(C(1,5)*n(1)^2+C(4,6)*n(2)^2+C(3,5)*n(3)^2+(C(1,4)+C(5,6))*n(1)*n(2)+(C(1,3)+C(5,5))*n(1)*n(3)+(C(3,6)+C(4,5))*n(2)*n(3));
L_22=(C(6,6)*n(1)^2+C(2,2)*n(2)^2+C(4,4)*n(3)^2+2*C(2,6)*n(1)*n(2)+2*C(4,6)*n(1)*n(3)+2*C(2,4)*n(2)*n(3));
L_23=(C(5,6)*n(1)^2+C(2,4)*n(2)^2+C(3,4)*n(3)^2+(C(4,6)+C(2,5))*n(1)*n(2)+(C(3,6)+C(4,5))*n(1)*n(3)+(C(2,3)+C(4,4))*n(2)*n(3));
L_33=(C(5,5)*n(1)^2+C(4,4)*n(2)^2+C(3,3)*n(3)^2+2*C(4,5)*n(1)*n(2)+2*C(3,5)*n(1)*n(3)+2*C(3,4)*n(2)*n(3));
Christoffel_mat=[L_11, L_12, L_13;
L_12,L_22,L_23;
L_13,L_23,L_33];
[ev,d]=eig(Christoffel_mat);
% eigvecs = polarisation vectors
pl=ev(:,3);
% eigvals ~ slowness (phase)
phase_vel3(r)=sqrt(rho./d(3,3)); % quasi-longitudinal
figure
polar(0:pi/180:2*pi,1./phase_vel3)
Any assistance how to visualise this in 3d would be appreciated.
Thanks
Dave :)

 Réponse acceptée

[moved this from Comment to Answer as I had originally intended]
Please post the simplest possible code fragment that clearly illustrates what you want to do. This is always a good idea, in order to raise the odss of getting a useful answer.
I cannot tell what you want to plot on the third dimension.
Matlab cannot directly plot in cylindrical or spherical coordinates, as far as I know. But it does have a built-in function to convert from cylindrical to Cartesian, which I use below.
z=0:200; theta=z*pi/25; r=1-z/400;
[x,y,z]=pol2cart(theta,r,z); %convert cylindrical to Cartesian coords
figure;
subplot(121); polar(theta,r); %polar plot
subplot(122); plot3(x,y,z); %3D plot
I hope that helps.

4 commentaires

Hi William,
Apologies for the clustered post. Maybe I am misunderstanding what exactly I am needing.
So, if I plot my slownes curve using
figure
polar(0:pi/180:2*pi,1./phase_vel3)
I get something that looks like this
WHat I am trying to achieve , but not sure how to go about it should look like the image below. Velocity varying in z-direction
Thanks
Dave :)
William Rose
William Rose le 10 Août 2022
I know you already know this, but anyway, let me say it anyway, to motivate the question I still have:
The polar plot you showed is a plot of r versus theta, where r=1/v. The second plot is a sphere that satisfies x^2+y^2+z^2=360, and the color of the sphere is assigned based on the z-coordinate of the surface.
I assume you want x and y to be assigned as the polar coordinates of theta and r, where r=1/v. What value do you want to assign to z, in your 3D plot?
David Harra
David Harra le 11 Août 2022
Hi William
I think that is my issue/ lack of understanding.
What you siad makes perfect sense, but I actually thought 1/v would be plotted in the z-direction as this shows the variation in velocity. I might need to revisit some of the maths etc.
Thanks for the input- its massively appreciated.
Thanks
Dave :)
William Rose
William Rose le 17 Août 2022
@David Harra, you're welcome.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur 2-D and 3-D Plots dans Centre d'aide et File Exchange

Produits

Version

R2021b

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by