82 lines
1.9 KiB
Mathematica
82 lines
1.9 KiB
Mathematica
|
% 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
|
||
|
|
|