diff --git a/SaturnRingsHurst.m b/SaturnRingsHurst.m new file mode 100644 index 0000000..73b5bca --- /dev/null +++ b/SaturnRingsHurst.m @@ -0,0 +1,81 @@ +% 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 +