PROGRAM RAINFLOW C C RAINFLOW.F VER 1.3 MARCH 20, 2014 C BY TOM IRVINE C EMAIL: TOM@VIBRATIONDATA.COM C C ASTM E 1049-85 (2005) RAINFLOW COUNTING METHOD C C23456789 C PARAMETER(MAX=10000000) C INTEGER, ALLOCATABLE :: T(:) INTEGER, ALLOCATABLE :: ST(:) C DOUBLE PRECISION, ALLOCATABLE :: Y(:),A(:),TEMP(:) DOUBLE PRECISION, ALLOCATABLE :: C(:) DOUBLE PRECISION, ALLOCATABLE :: CC(:) DOUBLE PRECISION, ALLOCATABLE :: AA(:,:) DOUBLE PRECISION, ALLOCATABLE :: A_MEAN(:,:),B(:,:) C DOUBLE PRECISION, ALLOCATABLE :: AVERAGEMEAN(:) DOUBLE PRECISION, ALLOCATABLE :: MAXMEAN(:) DOUBLE PRECISION, ALLOCATABLE :: MINMEAN(:) DOUBLE PRECISION, ALLOCATABLE :: MAXPEAK(:) DOUBLE PRECISION, ALLOCATABLE :: MINVALLEY(:) DOUBLE PRECISION, ALLOCATABLE :: MAXAMP(:) DOUBLE PRECISION, ALLOCATABLE :: AVERAGEAMP(:) C DOUBLE PRECISION AMP(MAX) DOUBLE PRECISION L(14) DOUBLE PRECISION AMAX,AAMAX DOUBLE PRECISION SLOPE1,SLOPE2 DOUBLE PRECISION BX,D C DOUBLE PRECISION P1,P2,TP1,TP2 DOUBLE PRECISION X,YY C INTEGER IFLAG,KFLAG INTEGER I,IP1,J,JP1,K,M,NUM INTEGER MSA,MSA_ORIG INTEGER TC INTEGER IOstatus INTEGER stat_alloc C CHARACTER :: B_file*20,C_file*20 C PRINT'(" ")' PRINT'(" ASTM E 1049-85 (2005) RAINFLOW COUNTING METHOD")' C CALL DATA_READ(AMP,NUM) C PRINT'(" ")' PRINT'(" Enter name of binned data output file: ")' READ(*,*)B_file C_file='amp_cycles.txt' OPEN(UNIT=11,FILE=B_file,ACTION='WRITE') OPEN(UNIT=12,FILE=C_file,ACTION='WRITE') C ALLOCATE(T(NUM)) ALLOCATE(Y(NUM)) ALLOCATE(A(NUM)) C DO I=1,NUM Y(I)=AMP(I) A(I)=0. T(I)=0 ENDDO C A(1)=Y(1) T(1)=1 K=2 C M=NUM WRITE(*,133)M 133 FORMAT(/,' Total number of input points = ',I10,/) C PRINT'(" Begin slope calculation ")' C C K = number of peaks C DO I=2,NUM-1 SLOPE1=( Y(I)-Y(I-1)) SLOPE2=(Y(I+1)-Y(I)) IF((SLOPE1*SLOPE2).LE.0)THEN A(K)=Y(I) T(K)=I K=K+1 ENDIF ENDDO C A(K)=Y(NUM) T(K)=T(K-1)+1 WRITE(*,512)K 512 FORMAT(/,' Number of Peaks = ',I10,/) C ALLOCATE(TEMP(K)) TEMP(1:K)=A(1:K) DEALLOCATE(A,stat=stat_alloc) ALLOCATE(A(K)) A=TEMP C TEMP(1:K)=T(1:K) DEALLOCATE(T,stat=stat_alloc) ALLOCATE(T(K)) T=TEMP C ALLOCATE(CC(K)) CC=A ALLOCATE(AA(K,2)) C DO I=1,K AA(I,1)=T(I) AA(I,2)=A(I) C WRITE(*,199)A(I) C199 FORMAT(G10.4) ENDDO DEALLOCATE(T,stat=stat_alloc) DEALLOCATE(A,stat=stat_alloc) C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C RULES FOR THIS METHOD ARE AS FOLLOWS: LET X DENOTE C RANGE UNDER CONSIDERATION Y, PREVIOUS RANGE ADJACENT TO X AND C S, STARTING POINT IN THE HISTORY. C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C AAMAX=0. ALLOCATE(B(NUM,4)) ALLOCATE(A_MEAN(NUM,4)) C A_MEAN=0. B=0. C C ST is the index status. C C ST(I)=0 array entry still active C ST(I)=1 array entry deleted C ALLOCATE(ST(K)) ST=0 C C Find AI C KV=1 I1=1 C I=1 J=2 JP1=3 IFLAG=1 IB=0 C prevAI=-999. prevAJ=-999. C AI=0. AJ=0. C DO WHILE (I.LE.K-3 .and. (prevAI.ne.AI.and.prevAJ.ne.AJ)) C prevAI=AI prevAJ=AJ C IP1=J C AI=AA(I,2) AIP1=AA(J,2) AJ=AA(J,2) AJP1=AA(JP1,2) C YY=(ABS(AIP1-AI)) X=(ABS(AJP1-AJ)) C C write(*,543)AI,AJ,AJP1 543 format(3(1x,g10.4)) IF(X.GE.YY .AND. YY.GT.0.)THEN IF(IFLAG.EQ.1)THEN B(KV,1)=YY B(KV,2)=0.5 B(KV,3)=AI B(KV,4)=AJ ST(I)=1 C WRITE(*,112)'p1a ',B(KV,1),B(KV,2),B(KV,3),B(KV,4) I=J ISTART=J+1 IEND=K-1 CALL FIND_NEXT(ISTART,IEND,ST,J,K) ISTART=J+1 IEND=K CALL FIND_NEXT(ISTART,IEND,ST,JP1,K) ELSE B(KV,1)=YY B(KV,2)=1 B(KV,3)=AI B(KV,4)=AJ ST(I)=1 ST(IP1)=1 C WRITE(*,112)'p1b ',B(KV,1),B(KV,2),B(KV,3),B(KV,4) ISTART=1 IEND=K-2 CALL FIND_NEXT(ISTART,IEND,ST,I,K) ISTART=I+1 IEND=K-1 CALL FIND_NEXT(ISTART,IEND,ST,J,K) ISTART=J+1 IEND=K CALL FIND_NEXT(ISTART,IEND,ST,JP1,K) ENDIF IF(YY.GT.AAMAX)THEN P1=AA(I,2) P2=AA(I+1,2) TP1=AA(I,1) TP2=AA(I+1,1) AAMAX=YY ENDIF KV=KV+1 C I1=I IFLAG=1 C ELSE C IFLAG=IFLAG+1 I=J ISTART=J+1 IEND=K-1 CALL FIND_NEXT(ISTART,IEND,ST,J,K) ISTART=J+1 IEND=K CALL FIND_NEXT(ISTART,IEND,ST,JP1,K) C ENDIF C ENDDO C write(*,*)' ref 890 ' C C23456789 C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C COUNT EACH RANGE THAT HAS NOT BEEN PREVIOUSLY COUNTED C AS ONE-HALF CYCLE. C PRINT'(" ")' C I1=1 ISTART=I1 IEND=K-1 CALL FIND_NEXT(ISTART,IEND,ST,I,K) C DO WHILE(I1.LE.K) ISTART=I+1 IEND=K CALL FIND_NEXT(ISTART,IEND,ST,J,K) C AI=AA(I,2) AJ=AA(J,2) C YY=ABS(AJ-AI) C ST(I)=1 IF(YY.GT.0.)THEN B(KV,1)=YY B(KV,2)=0.5 B(KV,3)=AI B(KV,4)=AJ C WRITE(*,112)'p2 ',B(KV,1),B(KV,2),B(KV,3),B(KV,4) C112 FORMAT(4(1X,G12.5)) ENDIF C IF(YY.GT.AAMAX)THEN P1=AI P2=AIP1 TP1=AA(I,1) TP2=AA(IP1,1) AAMAX=YY ENDIF KV=KV+1 I=J I1=I+1 ENDDO C write(*,*)' ref 900 ' C23456789 C DEALLOCATE(ST,stat=stat_alloc) C C AMAX=MAX(Y)-MIN(Y) C PRINT'(" BEGIN BIN SORTING ")' C AMAX=MAXVAL(B(1:M,1)) L(1)=0 L(2)=2.5 L(3)=5 L(4)=10 L(5)=15 L(6)=20 L(7)=30 L(8)=40 L(9)=50 L(10)=60 L(11)=70 L(12)=80 L(13)=90 L(14)=100 L=L*AMAX/100 C C write(*,*)' AMAX= ',AMAX C NUM=SIZE(L)-1 C ALLOCATE(AVERAGEMEAN(NUM)) ALLOCATE(MAXMEAN(NUM)) ALLOCATE(MINMEAN(NUM)) ALLOCATE(MAXPEAK(NUM)) ALLOCATE(MINVALLEY(NUM)) ALLOCATE(MAXAMP(NUM)) ALLOCATE(AVERAGEAMP(NUM)) ALLOCATE(C(NUM)) C DO I=1,NUM C(I)=0. AVERAGEMEAN(I)=0 MAXMEAN(I)=-1.0E+09 MINMEAN(I)= 1.0E+09 C MAXPEAK(I)=-1.0E+09 MINVALLEY(I)= 1.0E+09 C MAXAMP(I)=0 AVERAGEAMP(I)=0 ENDDO C C print'(" ref 5 ")' C DO I=1,KV-1 YY=B(I,1) Yh=YY/2. C WRITE(12,535)Yh,B(I,2) 535 FORMAT(G10.4,2X,F3.1) C DO IJK=1,NUM IF((YY.GE.L(IJK)).AND.(YY.LE.L(IJK+1)))THEN C(IJK)=C(IJK)+B(I,2) BM=(B(I,3)+B(I,4))/2 IF(B(I,3).GT.MAXPEAK(IJK))THEN MAXPEAK(IJK)=B(I,3) ENDIF IF(B(I,4).GT.MAXPEAK(IJK))THEN MAXPEAK(IJK)=B(I,4) ENDIF IF(B(I,3).LT.MINVALLEY(IJK))THEN MINVALLEY(IJK)=B(I,3) ENDIF IF(B(I,4).LT.MINVALLEY(IJK))THEN MINVALLEY(IJK)=B(I,4) ENDIF C AVERAGEAMP(IJK)=AVERAGEAMP(IJK)+B(I,1)*B(I,2) AVERAGEMEAN(IJK)=AVERAGEMEAN(IJK)+BM*B(I,2) C IF( BM .GT. MAXMEAN(IJK))THEN MAXMEAN(IJK)=BM ENDIF IF( BM .LT. MINMEAN(IJK))THEN MINMEAN(IJK)=BM ENDIF C IF(B(I,1).GT.MAXAMP(IJK))THEN MAXAMP(IJK)=B(I,1) ENDIF EXIT ENDIF ENDDO ENDDO C DO IJK=1,NUM IF( C(IJK).GT.0.)THEN AVERAGEAMP(IJK)=AVERAGEAMP(IJK)/C(IJK) AVERAGEMEAN(IJK)=AVERAGEMEAN(IJK)/C(IJK) ENDIF ENDDO C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C23456789 WRITE(11,555) WRITE(11,556) WRITE(11,557) 555 FORMAT(' AMPLITUDE = (PEAK-VALLEY)/2 ',/) 556 FORMAT(' RANGE LIMITS CYCLE AVERAGE MAX 1 MIN AVE MAX MIN MAX') 557 FORMAT(' (UNITS) COUNTS AMP AMP 1 MEAN MEAM MEAN VALLEY PEAK',/) C PRINT'(" ")' PRINT'(" AMPLITUDE = (PEAK-VALLEY)/2 ")' PRINT'(" ")' PRINT'(" RANGE LIMITS CYCLE AVERAGE MAX")' PRINT'(" (UNITS) COUNTS AMP AMP")' C MAXAMP=MAXAMP/2 AVERAGEAMP=AVERAGEAMP/2 C DO I=1,NUM J=NUM+1-I C IF(C(J).LT.0.5)THEN AVERAGEAMP(J)=0. MAXAMP(J)=0. MINMEAN(J)=0. AVERAGEMEAN(J)=0. MAXMEAN(J)=0. MINVALLEY(J)=0. MAXPEAK(J)=0. ENDIF C C23456789 WRITE(*,207)L(J),L(J+1),C(J),AVERAGEAMP(J),MAXAMP(J) 207 FORMAT(1X,G10.4,' to ',G10.4,1X,F8.1,2X,G10.4,2X,G10.4) WRITE(11,201)L(J),L(J+1),C(J),AVERAGEAMP(J),MAXAMP(J), 1 MINMEAN(J),AVERAGEMEAN(J),MAXMEAN(J),MINVALLEY(J),MAXPEAK(J) 201 FORMAT(1X,G10.4,' to ',G10.4,2X,F8.1,1X,7(1X,G10.4)) ENDDO CLOSE(11) CLOSE(12) C WRITE(*,301)AAMAX 301 FORMAT(/,' MAX RANGE = ',G10.4,/) C TC=SUM(C) WRITE(*,401)TC 401 FORMAT(/,' TOTAL CYCLES = ',I10,/) C PRINT'(" ")' PRINT '(" Calculate Relative Damage Index? 1=yes 2=no ")' READ(*,*)IDI C IF(IDI.EQ.1)THEN PRINT'(" ")' PRINT '(" Enter fatigue exponent ")' READ(*,*)BX C D=0. DO I=1,KV-1 D=D+B(I,2)*(0.5*B(I,1))**BX ENDDO C WRITE(*,717)D 717 FORMAT(/,' Damage Index = ',G10.4,/) C ENDIF C C DEALLOCATE(Y,stat=stat_alloc) DEALLOCATE(TEMP,stat=stat_alloc) DEALLOCATE(CC,stat=stat_alloc) DEALLOCATE(A_MEAN,stat=stat_alloc) DEALLOCATE(B,stat=stat_alloc) DEALLOCATE(C,stat=stat_alloc) C DEALLOCATE(AVERAGEMEAN,stat=stat_alloc) DEALLOCATE(MAXMEAN,stat=stat_alloc) DEALLOCATE(MINMEAN,stat=stat_alloc) DEALLOCATE(MAXPEAK,stat=stat_alloc) DEALLOCATE(MINVALLEY,stat=stat_alloc) DEALLOCATE(MAXAMP,stat=stat_alloc) DEALLOCATE(AVERAGEAMP,stat=stat_alloc) C WRITE(*,*) WRITE(*,*)' Output files: ' WRITE(*,*)' ',B_file,' - binned data ' WRITE(*,*)' ',C_file,' - amplitude & cycles ' C STOP END C SUBROUTINE FIND_NEXT(ISTART,IEND,ST,IRETURN,K) INTEGER JK,ISTART,IEND,IRETURN INTEGER ST(K) IF(ISTART.LE.0)ISTART=1 DO JK=ISTART,IEND IF(ST(JK).EQ.0)THEN IRETURN=JK EXIT ENDIF ENDDO RETURN END C SUBROUTINE DATA_READ(AMP,NUM) PARAMETER(MAX=10000000) DOUBLE PRECISION AMP(MAX) DOUBLE PRECISION TT,A INTEGER I,IC,NUM INTEGER ios,IOstatus,stat_alloc CHARACTER :: A_file*20 C PRINT'(" ")' PRINT'(" The input file must be a time history. ")' PRINT'(" Select format: ")' PRINT'(" 1=amplitude ")' PRINT'(" 2=time & amplitude ")' READ(*,*)IC C ios=1 DO WHILE(ios.NE.0) PRINT '( " Enter the input filename:")' READ(*,*)A_file C OPEN(UNIT=10,FILE=A_file,STATUS='OLD',ACTION='READ',IOSTAT=ios) C IF(ios.EQ.0)THEN WRITE(*,720) 720 FORMAT(/,'File opened.',/) ELSE WRITE(*,721) 721 FORMAT(/,' Error: File does not exist.',/) ENDIF C ENDDO C DO I=1,MAX IF(IC.EQ.1)THEN READ(10,*,IOSTAT=IOstatus) A ELSE READ(10,*,IOSTAT=IOstatus) TT,A ENDIF IF (IOstatus.EQ.0) THEN NUM=I ELSE EXIT ENDIF ENDDO C REWIND(10) C DO I=1,NUM IF(IC.EQ.1)THEN READ(10,*,IOSTAT=IOstatus) AMP(I) ELSE READ(10,*,IOSTAT=IOstatus) TT,AMP(I) ENDIF ENDDO C C WRITE(*,404)NUM C404 FORMAT(/,' NUM = ',I8,/) C CLOSE(10) RETURN END