Extract area from elipsoidal/cirle trajectory.

Dear all, how are you?

I created the following code to simulate an anticlockwise rotating trajectory forced by an external velocity. Now, I want to calculate the area each time the trajectory makes a closed circle, as described in the attached image. I would appreciate any help in finding the best way to identify the areas inside the circles for each loop/circle in the trajectory.

For circle detection, I used the polyshape function as described in the attached figure. Unfortunately, this requires me to cut the series each time a circular trajectory is identified, and then manually choose the area I need, which takes too much time

Thanks for your time!

clear all,close all;

Fs = 1/0.5;                     
duration_hours = 25 * 24;   
t = 0:1/Fs:duration_hours;  
period = 20;                
phase = -pi/2;

num_periods = floor(duration_hours / period);
amplitudes = 0.1/2 + 0.1 * rand(1, num_periods); % Random amplitudes between 0.1 and 0.3

y = zeros(size(t));
for i = 1:num_periods
    start_idx = (i-1) * period * Fs + 1;
    end_idx = i * period * Fs;
    uI(start_idx:end_idx) = amplitudes(i) * sin(2*pi*(t(start_idx:end_idx) - t(start_idx))/period);
end

amplitudes = 0.1/2 + 0.1 * rand(1, num_periods); % Random amplitudes between 0.1 and 0.3
for i = 1:num_periods
    start_idx = (i-1) * period * Fs + 1;
    end_idx = i * period * Fs;
    vI(start_idx:end_idx) = amplitudes(i) * sin(2*pi*(t(start_idx:end_idx) - t(start_idx))/period);
end

uI = uI(10:end); 
vI = vI(1:length(uI));

cv = complex(uI*100,vI*100);
dt=0.5/24;         
dt=dt*3600*24;      
cxI=cumsum(cv)*dt;   
cxI=cxI/100/1000;     

x = real(cxI + complex(cumsum(repmat(0.1,1191,1)),0)' );
y = imag(cxI + complex(cumsum(repmat(0.1,1191,1)),0)' );
clearvars -except x y;
plot(x,y),axis equal;

% for the first circle detected
pgonAll = polyshape(x(1:48),y(1:48),'Simplify',false);
pgonAll = simplify(pgonAll);
pgonEach = regions(pgonAll);

figure
plot(x(1:48),y(1:48))
hold on
plot(pgonEach(3))

Answer 1

[Matt J]

load data

ymin=min(y);
[xx,yy]=deal(x,y);

xx(end+1)=x(end); yy(end+1)=ymin;

p=regions(polyshape(xx,yy));

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

h=arrayfun(@(P)min(P.Vertices(:,2)), p );

p=p(h>0.001+ymin);

Areas=area(p)'

Areas = 1x29

1.8247 0.9686 0.0081 1.0051 0.5726 0.0032 0.1837 0.5910 2.8787 1.0814 0.0241 1.7894 2.5199 0.7346 1.0110 0.0168 0.3967 2.1162 1.7427 2.5049 0.7968 1.8758 1.2651 0.6767 1.3316 1.8829 0.2996 1.1132 0.0252

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

plot(x,y); hold on
plot(p,'FaceColor','r'); hold off
axis([50,100,-5,5]); daspect([1,1,1])