Principal axis + phase and inclination.

Hello, I would appreciate any help/assistance in calculating the following variables using the attached code.:

Semi-major axis (Mayor), Semi-minor axis (Minor), Inclination and Phase

of Horizontal velocities u and v are represented as complex numbers complex(u,v) with the purpose of plotting the ellipse as shown in the attached figure.

%horizontal velocity ellipse rotating clockwise.

% Example data for time series of u and v velocities
u = [ -3.4500   -2.3500    0.8500   -1.3500    0.2000    0.7500    3.8500    5.0000    4.9000    8.0000    6.2000   14.6500   21.8000   14.8500   13.0500];
v = [ 13.1000   12.3000    8.0000    3.8000    2.9000    1.8500    3.7500    1.0000    3.3000    4.8500   -3.9000    9.9000   11.4000   16.5500   20.3000];

% Combine u and v velocities into complex vector
wind_complex = complex(u, v);

% Compute covariance matrix
cov_matrix = cov(real(wind_complex), imag(wind_complex));

% Find eigenvalues and eigenvectors
[eigenvectors, eigenvalues] = eig(cov_matrix);

% Major and minor axes
major_axis_length = max(diag(eigenvalues));
minor_axis_length = min(diag(eigenvalues));

% Calculate inclination and phase
inclination = atan2(eigenvectors(2,1), eigenvectors(1,1)); % angle in radians
phase = atan2(mean(imag(wind_complex)), mean(real(wind_complex))); % angle in radians

% Convert inclination and phase to degrees
inclination_degrees = rad2deg(inclination);
phase_degrees = rad2deg(phase);

% Display results
disp(['Major Axis Length: ', num2str(major_axis_length)]);
disp(['Minor Axis Length: ', num2str(minor_axis_length)]);
disp(['Inclination (degrees): ', num2str(inclination_degrees)]);
disp(['Phase (degrees): ', num2str(phase_degrees)]);
Major Axis Length: 67.5505
Minor Axis Length: 29.7218
Inclination (degrees): -54.0886
Phase (degrees): 51.446
>> 

Answer 1

[William Rose]

@Claudio Iturra,

I assume you want to do a best-fit ellipse to the wind data. Here is a link to my notes on how to find the best-fit ellipse to a set of x,y points. The notes refer to center-of-presssure data , but it could be any set of x,y points.

https://www1.udel.edu/biology/rosewc/kaap686/reserve/cop/center%20of%20pressure.html

I apologize for the poor quality of the file. This old file online is the only copy of my notes that I can locate at the moment. Key points from the notes:

The following procedure finds the “best fit” ellipse.

• Start with N pairs (x1, y1), (x2, y2), …, (xN, yN).
• Compute the means (xm, ym).
• Compute the covariance matrix Cov= (cxx , cxy ; cyx , cyy) where
• cxx= [ ∑(xi-xm)² ] / (N-1)
• cxy = cyx = [ ∑(xi-xm)(yi-ym) ] / (N-1)
• cyy = [ ∑(yi-ym)² ] / (N-1)
• Find the eigenvalues (λ1, λ2) of the covariance matrix (λ1 > λ2).
• Find the unit eigenvectors (v11 ; v21) and (v12 ; v22) of the covariance matrix.
• The best fit ellipse has its center at (xm, ym).
• The semimajor axis of the 95% ellipse has direction (v11; v21) and length sqrt(5.99 λ1).
• The semiminor axis of the 95% ellipse has direction (v12; v22) and length sqrt(5.99 λ2).
• The area of the 95% ellipse is A = 5.99 π sqrt(λ1 λ2)
• The angle between the semimajor axis and the y-axis is θ = sin-1(v11).

Plot your wind data:

% Example data for time series of u and v velocities
u = [-3.45 -2.35 0.85 -1.35 0.20 0.75 3.85 5.00 4.90 8.00  6.20 14.65 21.80 14.85 13.05];
v = [13.10 12.30 8.00  3.80 2.90 1.85 3.75 1.00 3.30 4.85 -3.90  9.90 11.40 16.55 20.30];
figure; plot(u,v,'k*'); 
grid on; xlabel('u (m/s)'); ylabel('v (m/s)'); hold on

Calculate best fit ellipse. I notice that you do not take the square root of the eigenvalues, in your calculations. You should. And you should scale them by a suitable factor, such as sqrt(5.99) to catch 95% of the points, or sqrt(4.61) to catch 90% of the points. (The scaling factors come from the chi squared distribution, and are based on the assumption that the noise in the data is Gaussian and independent for u and v.)

[William Rose]

@Claudio Iturra,

Here is a script that computes the semimajor and semiminor axis length and the inclination of the majopr axis with respect to the +x axis. It makes the figure below.

The script also displays these results:

>> bestFitEllipse

Major Axis Length = 20.14.

Minor Axis Length = 13.36.

Inclination = -144.1 degrees.

I am not sure what you mean by the phase angle. If the data lay on or close to an ellipse, and if there were an associated time vector, then I suppose the phase could be defined as the angle of the vector from the ellipse center to the point at time zero, or to the first point in the data set. But this data do not appear to lie on or close to an ellipse.