How to calculate the mean integrated curvature of an ellipsoid?

11 visualizaciones (últimos 30 días)
I am using the following code to compute the surface area, mean integrated curvature, and Euler characteristic of the ellipsoid parameterized by r(x,y). Here 'x' and 'y' represent the angles $\theta$ and $\phi$ in spherical coordinates. The code is
clc; clear;
%ellipsoid radii
a=value1; b=value2; c=value2;
syms x y
%parametric vector equation of the ellipsoid
r=[a*cos(y)*sin(x), b*sin(x)*sin(y), c*cos(x)];
%partial derivatives
r_x=diff(r,x);
r_xx=diff(r_x,x);
r_y=diff(r,y);
r_yy=diff(r_y,y);
r_xy=diff(r_x,y);
%normal vector to the surface at each point
n=cross(r_x,r_y)/norm(cross(r_x,r_y));
%coefficients of the second fundamental form
E=dot(r_x,r_x);
F=dot(r_x,r_y);
G=dot(r_y,r_y);
L=dot(r_xx,n);
M=dot(r_xy,n);
N=dot(r_yy,n);
A=E*G-F^2;
K=((E*N)+(G*L)-(2*F*M))/A;
H=(L*N-M^2)/A;
dA=sqrt(A);
dC=K*dA;
dX=H*dA;
%converting syms functions to matlab functions
dA_=matlabFunction(dA); dC_=matlabFunction(dC); dX_=matlabFunction(dX);
% numerical integration
S = integral2(dA_, 0, 2*pi, 0, pi); %total surface area
C = (1/2)*integral2(dC_, 0, 2*pi, 0, pi); % medium integrated curvature
X = (1/(2*pi))*integral2(dX_, 0, 2*pi, 0, pi); % Euler characteristics
S and X give me very well, but the problem is with C, which always gives me very small, of the order of 1e-16 to 1e-14, which is absurd since C must be of the order of 4*pi*r (for the case of a sphere a=b=c=r). Could someone tell me what is wrong?

Respuesta aceptada

Milan Bansal
Milan Bansal el 11 de Sept. de 2024
Editada: Milan Bansal el 13 de Sept. de 2024
Hi Johan Zúñiga
To calculate the mean integrated curvature of an ellipsoid, you need to evaluate the integral of the mean curvature HHH over the surface of the ellipsoid. For an ellipsoid, the mean curvature H at any point is the average of the two principal curvatures and
For a general surface, the mean curvature H is given by:
In the case of an ellipsoid, the mean integrated curvature can be computed by:
For a sphere with radius rrr, the mean integrated curvature is well-known to be:
For an ellipsoid with semi-axes a, b, and c, the exact analytical solution is not as straightforward as for the sphere. However, you can compute it numerically.
  • Mean Curvature H: Using the elements of the first and second fundamental forms, calculate the mean curvature H at each point on the ellipsoid surface using:
  • Numerically Integrate: Integrate the mean curvature H over the surface of the ellipsoid using integral2.
Here is how you can modify your code to calculate the curvature correctly:
clc; clear;
% Ellipsoid radii
a = 2; b = 5; c = 5;
syms x y
% Parametric vector equation of the ellipsoid
r = [a*cos(x)*sin(y), b*sin(x)*sin(y), c*cos(y)];
% Partial derivatives with respect to theta and phi
r_x = diff(r, x);
r_xx = diff(r_x, x);
r_y = diff(r, y);
r_yy = diff(r_y, y);
r_xy = diff(r_x, y);
% Normal vector to the surface
n = cross(r_x, r_y) / norm(cross(r_x, r_y));
% First fundamental form coefficients
E = dot(r_x, r_x);
F = dot(r_x, r_y);
G = dot(r_y, r_y);
% Second fundamental form coefficients
L = dot(r_xx, n);
M = dot(r_xy, n);
N = dot(r_yy, n);
% Determinant of the first fundamental form (area element factor)
A = E*G - F^2;
K=((E*N) + (G*L) - (2*F*M)) / A;
X=(L*N - M^2) / A;
% Mean curvature
H = (L*G - 2*M*F + N*E) / (2*A);
% Area element (sqrt(A))
dA = sqrt(A);
dC = K*dA;
dX = X*dA;
% Mean curvature element to integrate
dH = H * dA;
%converting syms functions to matlab functions
dA_ = matlabFunction(dA);
dC_ = matlabFunction(dC);
dX_ = matlabFunction(dX);
dH_ = matlabFunction(dH, 'Vars', [x, y]);
% Numerical integration over the parameter domain (theta = 0 to 2*pi, phi = 0 to pi)
S = integral2(dA_, 0, 2*pi, 0, pi) %total surface area
S = 200.0445
C_H = integral2(dH_, 0, 2*pi, 0, pi)
C_H = 52.3037
X = (1/(2*pi))*integral2(dX_, 0, 2*pi, 0, pi) % Euler characteristics
X = 2.0000
Hope this helps!

Más respuestas (0)

Categorías

Más información sobre Surface and Mesh Plots en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by