disp(' ');
disp(' rainflow.m  ver 1.1  February 24, 2010 ');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' ASTM E 1049-85 (2005) Rainflow Counting Method ');
disp(' ');
%
clear a;
clear cc;
clear THM;
clear C;
clear y;
%
disp(' The input file must be a time history. ')
disp(' Select format: ');
disp('   1=amplitude ');
disp('   2=time & amplitude ');
ic=input(' ');
%
if(ic==1)
%
    disp(' ')
    disp(' Select file input method ');
    disp('   1=external ASCII file ');
    disp('   2=file preloaded into Matlab ');
    disp('   3=Excel file');
    file_choice = input('');
%
    if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g',[1 inf]);
        fclose(fid);
        THM=THM';
    end
    if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
    end
    if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
    end
    y=double(THM(:,1));
    figure(1);
    plot(y);
    grid on;
%
else	
%
    disp(' ')
    disp(' Select file input method ');
    disp('   1=external ASCII file ');
    disp('   2=file preloaded into Matlab ');
    disp('   3=Excel file');
    file_choice = input('');
%
    if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        fclose(fid);
        THM=THM';
    end
    if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
    end
    if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
    end
    y=double(THM(:,2)); 
    plot(y);
    grid on;
%
end
%
disp(' ');
disp(' The bin range is pre-set to 1 ');
disp(' ');
disp(' Multiply data by scale factor?');
disp(' 1=yes  2=no');
iscale=input(' ');
if(iscale==1)
    disp(' Enter scale factor ');
    scale=input(' ');
    y=y*scale;
    figure(1);
    plot(y);   
    grid on;
end
%
% 
% test case points for ASTM E 1049-85 example.
%
% a(1)=-2;  % A
% a(2)=1;   % B
% a(3)=-3;  % C
% a(4)=5;   % D
% a(5)=-1;  % E
% a(6)=3;   % F
% a(7)=-4;  % G
% a(8)=4;   % H
% a(9)=-2;  % I
%
k=1;
a(1)=y(1);
k=2;
m=max(size(y))-1;
for(i=2:m)
    slope1=(  y(i)-y(i-1));
    slope2=(y(i+1)-y(i));
    if((slope1*slope2)<=0)
        a(k)=y(i);
        k=k+1;
    end
end
a(k)=y(m+1);
cc=a;
%
num=round(max(a)-min(a))+1;
C=zeros(num,1);
%
n=1;
i=1;
j=2;
%
while(1)
    Y=round(abs(a(i)-a(i+1)));
    X=round(abs(a(j)-a(j+1)));
    out1=sprintf(' ai=%g aj=%g Y=%g X=%g',a(i),a(j),Y,X);
    disp(out1);
%
    if(X>=Y & Y>0)
        if(i==1)
           C(Y)=C(Y)+0.5;
           a(1)=[];
        else
           C(Y)=C(Y)+1;         
           a(i+1)=[];
           a(i)=[];
        end
        C
        i=1;
        j=2;    
    else
        i=i+1;
        j=j+1;
    end
%
    if((j+1)>max(size(a)))
        break;
    end
end
%
N=max(size(a));
disp(' ');
for(i=1:N-1) 
    Y=round(abs(a(i)-a(i+1)));
    if(Y>0)
        C(Y)=C(Y)+0.5;
    end
end
%
N=max(size(C));
disp(' ');
disp('  Range    Cycle ');
disp(' (units)   Counts ');
for(i=1:N)
    j=N+1-i;
    out1=sprintf('\t %d  %6.3g ',j,C(j));
    disp(out1);
end
%
figure(2);
h=bar(C);
grid on;
title('Rainflow');
ylabel('Cycle Counts');
xlabel('Range');
%