Hurst-exponent/SaturnRingsHurst.m

82 lines
1.9 KiB
Matlab

% Exponente de Hurst en anillos de Saturno
% Programa para estudiar exponente de Hurst y dimensión fractal de la estructura de anillos
% Trabajo de Tesis - UNLP
% Horacio Daniel Salomone
% Start
%Image reading
AS1=imread('PIA22418.jpg');
%Correct to a circle
image(AS1), axis image;
%Shows the figure
subplot(2,2,1);
imshow(AS1),title('Imágen Original-Cassini - Anillos de Saturno 5');
%Cuts the figure and shows the cut
AS2=imcrop;
subplot(2,2,2);
imshow(AS2),title('Imagen Recortada-Anillos de Saturno 5');
%Loop to get the grayscale profiles
%image size AS2
sz=size(AS2);
n=sz(1,1);
%Intensity profile
Perfil='Vectores Perfil Intensidad.txt';
fid=fopen(Perfil,'a+');
formatSpec='%9.5f\r\n';
%Int vectors
for i=1:n
Int=improfile(AS2,[1 1018],[1 (0+i)]);
%Recorta ceros del vector Int
Int(isnan(Int))=0;
%Int(Int==0)=[];
l=1;
while Int(l)==0
l=l+1;
end
L=length(Int);
r=1;
while Int(L-r)==0
r=r+1;
end
Int(L-r+1:L)=[];
Int(1:l)=[];
end
%Hurst exponent calculation
for i=1:120
m=floor(L/(2*i));
for j=1:i
r=Int(1+(j-1)*m:j*m);
M=mean(r);
x=(r-M);
V=cumsum(x);
R(j)=max(V)-min(V);
S(j)=std(r);
end
tau(i)=m;
RS(i)=mean(R./S);
end
%Linear regression to obtain Hurst
subplot (2,2,3);
plot(log10(tau),log10(RS),'+')
xlabel('log(\tau)','FontSize',12)
ylabel('log(R/S)','FontSize',12)
hold on
q=polyfit(log10(tau),log10(RS),1);
t=0.5:0.1:3;
y=q(1)*t+(q(2));
plot(t,y,'r','LineWidth',2)
%text(1,2,['y =' num2str(q(1)),' x ' num2str(q(2))],'FontSize',12)
hold off