disp(' ');
disp(' rainflow_main.m   ver 1.6   July 9, 2014');
disp(' ');
disp(' by Tom Irvine');
disp(' ');
disp(' ASTM E 1049-85 (2005) Rainflow Counting Method ');
disp(' ');
%
clear B;
clear a;
clear cc;
clear THM;
clear THM_hold;
clear t;
clear C;
clear y;
clear L;
clear max;
clear size;
clear m;
clear k;
clear iscale;
clear length;
clear slope1;
clear slope2;
clear BIG;
clear dat;
clear data_s;
%
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)); 
    THM_hold=THM;
    plot(y);
    grid on;
%
end
clear THM;
%
m=length(y)-1;
%
out1=sprintf(' total input points =%d ',m);
disp(out1);
%
dchoice=-1.;
%
exponent=1;
%
disp(' ');
disp(' Calculate Relative Damage Index?  1=yes 2=no ');
ic=input(' ');
if(ic==1)
    dchoice=1.;  % needs to be double 
%
    disp(' ');
    disp(' Enter fatigue exponent ');
    exponent=input(' ');
%
end
%
[L,C,AverageAmp,MaxAmp,AverageMean,MinValley,MaxPeak,D,ac1,ac2,cL]...
                                         =rainflow_ac_mex(y,dchoice,exponent);
%
sz=size(ac1);
if(sz(2)>sz(1))
    ac1=ac1';
    ac2=ac2';
end
%
ncL=int64(cL);
%
amp_cycles=[ ac1(1:ncL) ac2(1:ncL) ];
clear ac1;
clear ac2;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
L=flipud(L);
C=flipud(C);
AverageAmp=flipud(AverageAmp);
MaxAmp=flipud(MaxAmp);
AverageMean=flipud(AverageMean);
MinValley=flipud(MinValley);
MaxPeak=flipud(MaxPeak);
%
  disp(' ');
  disp(' Round the cycle and amplitude values to nearest integer ? ');
  disp(' 1=yes 2=no');
  rv=input(' ');
%
  clear BIG;
%
  N=max(size(C));
%
  ijk=1;
  for i=1:N-1
      if( abs(L(i))>1.0e-09 && abs(L(i+1))>1.0e-09)
          ijk=ijk+1;
      end
  end
  N=ijk;
%
  BIG=zeros(N,8);
  disp(' ');
  disp('  Amplitude = (peak-valley)/2 ');
  disp(' ');
  disp('        Range Limits         Cycle      Average     Max      Average   Min     Max ');
  disp('          (units)            Counts      Amp        Amp      Mean      Valley  Peak');
%
  MaxAmp=MaxAmp/2;
  AverageAmp=AverageAmp/2;
%
  for j=1:N
%
      if(C(j)==0)
         AverageAmp(j)=0.;
         MaxAmp(j)=0.; 
         AverageMean(j)=0.;
         MinValley(j)=0.;
         MaxPeak(j)=0.; 
      end
%
      if(rv==2)
          out1=sprintf('\t %7.4g to %7.4g \t %g \t %6.3g \t %6.3g \t %6.3g\t %6.3g\t %6.3g',L(j+1),L(j),C(j),AverageAmp(j),MaxAmp(j),AverageMean(j),MinValley(j),MaxPeak(j));
      else
          out1=sprintf('\t %7.4g to %7.4g \t %g \t %6.3g \t %6.3g \t %6.3g\t %6.3g\t %6.3g',L(j+1),L(j),round(C(j)),round(AverageAmp(j)),round(MaxAmp(j)),round(AverageMean(j)),round(MinValley(j)),round(MaxPeak(j)));
      end
      disp(out1);
      BIG(j,1)=L(j+1);
      BIG(j,2)=L(j);
      BIG(j,3)=C(j);
      BIG(j,4)=AverageAmp(j);
      BIG(j,5)=MaxAmp(j);
      BIG(j,6)=AverageMean(j);
      BIG(j,7)=MinValley(j); 
      BIG(j,8)=MaxPeak(j);   
  end  
%
TC=sum(C);
if(rv==2)
   out1=sprintf('\n Total Cycles =%g \n',TC);
else
   out1=sprintf('\n Total Cycles =%g \n',round(TC));   
end
disp(out1);
%
if(ic==1)
    disp(' ');
    out1=sprintf(' Relative Damage Index = %9.4g ',D(1));
    disp(out1);
end    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ');
disp(' Present data in Matlab uitable:  1=yes 2=no ')
ip=input(' ');
%
close all hidden;
%
if(ip==1)
%
%
    table1 = uitable;
    set(table1,'ColumnWidth',{26})
    f = figure('Position',[100 100 900 350]);
    dat =  BIG(:,1:8);
columnname =   {'Lower Range','Upper Range','Cycles','Ave Amp',...
       'Max Amp','Ave Mean','Min Valley','Max Peak' };
columnformat = {'numeric', 'numeric','numeric','numeric','numeric',...
       'numeric','numeric','numeric'};
columneditable = [false false false false false false false false];   
%
sz=size(BIG);
%
for i = 1:sz(1)
    for j=1:sz(2)
            if(j==3)
                tempStr = sprintf('%10.1f', dat(i,j));                
            else
                tempStr = sprintf('%8.4g', dat(i,j));
            end    
            data_s{i,j} = tempStr; 
    end    
end
%
table1 = uitable('Units','normalized','Position',...
            [0.1 0.1 0.9 0.9], 'Data', data_s,... 
            'ColumnName', columnname,...
            'ColumnFormat', columnformat,...
            'ColumnEditable', columneditable,...
            'RowName',[]);
%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
fig_num=3;
figure(fig_num);
fig_num=fig_num+1;  
if(ic==1)
    plot(y);
else
    plot(THM_hold(:,1),THM_hold(:,2));
end
grid on;
%
  figure(fig_num);
  fig_num=fig_num+1;
  h=bar(C);
  grid on;
  title('Rainflow');
  ylabel('Cycle Counts');
  xlabel('Range');
%
%
RM=BIG(:,1:3);
%
disp('Output array: ');
disp('   amp_cycles  (amplitude & cycles)');