Program UPSAMPLE
C     Was TD_SNB_STD5  !
C     Was TD_SNB_STD4
C     Was IMU_1
	IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER(NPT=4800000,NPT2=10000,NUMCOL=6)
	DIMENSION BB(3,4),DISMAX(NPT2),VELMAX(NPT2),ACCMAX(NPT2)
     &,TDMAX(NPT2),TVMAX(NPT2),TAMAX(NPT2),DDXT(NPT),DDX(NPT)
     &,DXT(NPT),DX(NPT),XT(NPT),X(NPT)
     &,DUM1(NUMCOL),DUM2(NUMCOL) !NUMCOL= Acc Columns in input file
	Character*64 cfname(10)
	Character*32 crname
	CHARACTER*4 CDAMP(4)
C
	COMMON/CON/ONEPI,TWOPI,D2R
C ***
C     Set up constants
C ***
	ONEPI=4.0D0*DATAN(1.0D0)
	TWOPI=2.0D0*ONEPI
	TWOPI2=TWOPI**2
	D2R=TWOPI/180.0D0  ! Convert Degrees to radians
	G=386.0D0          ! Acceleration due to gravity in./sec^2
	FCON=SQRT(G)/TWOPI ! =3.12848....
C ***
C              ### Define File Name to ID Run ###
C ***
      crname='FM2_RateA'
C ***
C     Check Format of Input File !!! - See Read Statements - 2 sets
C              ### Input File Name ###
C ***
      cfname(6)='FM2_RateA.dat'
	OPEN(UNIT=12,FILE=cfname(6),status='old') ! Input
C***
C    Input can be defined analytically in VAL subroutine
C     Define 1/2 sine input pulse parameters
c	TONE=1.25D-4      ! Pulse time 1/2 sine input (sec) 4000 Hz
c	WIN=ONEPI/TONE
c	UDD=7.5D3*G		 ! 7,500 Gs peak acceleration
c	VDD=0.636D0*UDD*TONE ! Velocity change
	cfname(1) =TRIM(crname)//'AT.out'
	cfname(2) =TRIM(crname)//'DT.out'
	cfname(3) =TRIM(crname)//'TT.csv'
	cfname(4) =TRIM(crname)//'AS.out'
	cfname(5) =TRIM(crname)//'DS.out'
	cfname(8) =TRIM(crname)//'VT.out'
	cfname(9) =TRIM(crname)//'VS.out'
	cfname(10)=TRIM(crname)//'TU.csv'
	OPEN(UNIT=7, FILE=cfname(1),status='unknown') ! AT
	WRITE(7,'(25X,A)') cfname(1),cfname(6)
	OPEN(UNIT=8, FILE=cfname(2),status='unknown') ! DT
	WRITE(8,'(25X,A)') cfname(2),cfname(6)
c
	OPEN(UNIT=9, FILE=cfname(3),status='unknown') ! TT
	WRITE(9,'(28X,A)') cfname(3),cfname(6)
	OPEN(UNIT=13,FILE=cfname(10),status='unknown') ! TU
	WRITE(13,'(28X,A)') cfname(10),cfname(6)
c
	OPEN(UNIT=10,FILE=cfname(4),status='unknown') ! AS
	WRITE(10,'(25X,A)') cfname(4),cfname(6)
	OPEN(UNIT=11,FILE=cfname(5),status='unknown') ! DS
	WRITE(11,'(25X,A)') cfname(5),cfname(6)
c
	OPEN(UNIT=14,FILE=cfname(8),status='unknown') ! VT
	WRITE(14,'(25X,A)') cfname(8),cfname(6)
	OPEN(UNIT=15,FILE=cfname(9),status='unknown') ! VS
	WRITE(15,'(25X,A)') cfname(9),cfname(6)
C***
C     Read input acceleration from unit 8
C       IRWSKP = No. of lines to skip
C       ICREAD = Column to get acceleration data from
C       XIFACT = Factor to multiply input (e.g. 386 if data is gs)
C       XOFACT = Factor to divide output (e.g. 386 if data is gs)
C       Upsample  by factor of L
C	  NPRINT=nn Print Step, Usually = L
C***
      L=8
	NPRINT=1 !L ! Print Step, Usually = L
C***  Next two lines use for data created by this program
c      IRWSKP=2 ! Lines to skip in data file header
c      ICREAD=1 ! Acc Column in Dum(X) Arrary - SET Parameter NUMCOL
      IRWSKP=0 ! Lines to skip in data file header
      ICREAD=6 ! Acc Column in Dum(X) Arrary - SET Parameter NUMCOL
      XIFACT= 1.0D0 !+386.0D0
      XOFACT= 1.0D0 !+386.0D0
c     IREAD =0 std input, time and 'NUMCOL' columns of data
c           =1 Packed data - upsampled output from this program
c           =2 Excalibur format input (no time column)
c           =3 acc only, single column, need to put in dt
      IREAD=1
	NAVE=0 !100 !IDINT(.001D0/DTI) No. of points to calc bias
C***
C     Read regular data, time and NUMCOLs of data
C     ### Set NUMCOL as a parameter !
C***
	IF(IREAD.EQ.0) THEN ! READ COLUMN "ICREAD"
       DO I=1,IRWSKP
	  READ(12,*) ! Skip title, Line I
	 ENDDO
	 READ(12,*) T1,DUM1
	 READ(12,*) T2,DUM2
       DDX(1)=DUM1(ICREAD)
       DDX(2)=DUM2(ICREAD)
	 DTI=T2-T1 ! ASSMUME THIS IS CORRECT
	 WRITE(*,19) T1,DDX(1),T2,DDX(2),DTI
	 WRITE(9,19) T1,DDX(1),T2,DDX(2),DTI
	 I=2
   20  CONTINUE
       I=I+1
	 READ(12,*,END=30,ERR=30) T,DUM1
	 DDX(I)=DUM1(ICREAD)
	 GO TO 20
   30  CONTINUE
       N1=I-1
	 WRITE(9,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(9,'(5X,A,I11)')    'No. of      Steps    ',N1
	 WRITE(*,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(*,'(5X,A,I11)')    'No. of      Steps    ',N1
C***
C	 Packed data - output from this program
C***
	ELSEIF(IREAD.EQ.1) THEN ! READ PACKED DATA
	 READ(12,*)
	 READ(12,*)
	 READ(12,*) DTI
	 T1=0.0D0
	 T2=T1+DTI
	 READ(12,*) DDX(1)
	 READ(12,*) DDX(2)
	 WRITE(*,19) T1,DDX(1),T2,DDX(2),DTI
	 WRITE(9,19) T1,DDX(1),T2,DDX(2),DTI
	 I=2
   21  CONTINUE
       I=I+1
	 READ(12,*,END=31,ERR=31) DDX(I)
	 GO TO 21
   31  CONTINUE
       N1=I-1
	 WRITE(9,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(9,'(5X,A,I11)')    'No. of      Steps    ',N1
	 WRITE(*,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(*,'(5X,A,I11)')    'No. of      Steps    ',N1
C***
C	 Excalibur input, dt=1.0D-6
C***
	ELSEIF(IREAD.EQ.2) THEN ! READ COLUMN "ICREAD", NO TIME COLUMN
       DO I=1,IRWSKP
	  READ(12,*) ! Skip title, Line I
	 ENDDO
	 READ(12,*) DUM1
	 READ(12,*) DUM2
       DDX(1)=DUM1(ICREAD)
       DDX(2)=DUM2(ICREAD)
C***   SPECIAL FOR EXCALIBUR IFV1
	 DTI=1.0D-6
	 T1=0.0D0
	 T2=T1+DTI
	 WRITE(*,19) T1,DDX(1),T2,DDX(2),DTI
	 WRITE(9,19) T1,DDX(1),T2,DDX(2),DTI
	 I=2
   22  CONTINUE
       I=I+1
	 READ(12,*,END=32,ERR=32) DUM1
	 DDX(I)=DUM1(ICREAD)
	 GO TO 22
   32  CONTINUE
       N1=I-1
	 WRITE(9,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(9,'(5X,A,I11)')    'No. of      Steps    ',N1
	 WRITE(*,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(*,'(5X,A,I11)')    'No. of      Steps    ',N1
C***
C	 Single column acc input
C***
	ELSEIF(IREAD.EQ.3) THEN ! READ SINGLE COLUMN
	 DTI=1.0D0/2.0D5 ! 200k Hz sample rate
	 T1=0.0D0
	 T2=T1+DTI
	 READ(12,*) DDX(1)
	 READ(12,*) DDX(2)
	 WRITE(*,19) T1,DDX(1),T2,DDX(2),DTI
	 WRITE(9,19) T1,DDX(1),T2,DDX(2),DTI
	 I=2
   23  CONTINUE
       I=I+1
	 READ(12,*,END=33,ERR=33) DDX(I)
	 GO TO 23
   33  CONTINUE
       N1=I-1
	 WRITE(9,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(9,'(5X,A,I11)')    'No. of      Steps    ',N1
	 WRITE(*,'(5X,A,1PE11.3)')'Start       Time     ',T1
	 WRITE(*,'(5X,A,I11)')    'No. of      Steps    ',N1
	ELSE
	 WRITE(*,*) 'Bad IREAD NO. ', IREAD
	 STOP
	ENDIF
C ***
C     Add some zeros to see damping 0.05 sec
c       This messes up Upsampling!
C ***
c	DO I=N1+1,N1+5*NAVE
c	  DDX(I)=0.0D0
c	ENDDO
c	N1=N1+5*NAVE
C ***
C     Subtract out bias
C ***
	DDXAVE=0.0D0
	DO I=1,NAVE
	 DDXAVE=DDXAVE+DDX(I)
	ENDDO
	IF(NAVE.GT.0) THEN
	  DDXAVE=DDXAVE/NAVE
	ELSE
	  DDXAVE=0.0D0
	ENDIF
	WRITE(9,'(5X,A,1PE11.3,A,6X,I11,4X,A)')
     & 'Bias is = ',DDXAVE,'  Gs ',NAVE,'  Points'
	WRITE(*,'(5X,A,1PE11.3,A,6X,I11,4X,A)')
     & 'Bias is = ',DDXAVE,'  Gs ',NAVE,'  Points'
	DO I=1,N1
	  DDX(I)=DDX(I)-DDXAVE
	ENDDO
C***  ########################
C     Upsample  by factor of L
C       DDXT is original  time history
C       DDX  is upsampled time history
C***  ########################
	NT=N1
      IF(L.GT.1) THEN
	  DO I=1,N1
	    DDXT(I)=DDX(I)
	  ENDDO
	  WRITE(9,'(5X,A,2X,3I11)')'Upsample'
     &    ,L,N1,L*N1-L-1
	  WRITE(*,'(5X,A,2X,3I11)')'Upsample'
     &    ,L,N1,L*N1-L+1
	  WRITE(*,*)'Start Upsamle'
	  DTI=DTI/FLOAT(L)
	  CALL UPSAMP(DDX,DDXT,N1,L)
C***
C     Echo upsampled data, use "compressed" input format
C***
  	  cfname(7)=TRIM(crname)//'.ups' ! Used only if L>1
	  OPEN(UNIT=14,FILE=cfname(7),status='unknown') ! UPsample Echo
	  WRITE(14,'(2X,2A)') cfname(6),cfname(7)
	  WRITE(14,'(5X,A,2X,3I11)')'Upsample'
     &    ,L,N1,L*N1-L-1
	   WRITE(14,'(1PE16.8)')DTI
	  DO I=1,N1
	   WRITE(14,'(1PE16.8)')DDX(I)
	  ENDDO
      ENDIF
C ***
C      Integrate the Original acceleration to velocity and displacement
C ***
	XT(1)=0.0D0
	DXT(1)=0.0D0
	DDXT(1)=DDX(1)*XIFACT
	DO I=2,N1
	  DDXT(I)=DDXT(I)*XIFACT ! Scale the input
	ENDDO
	DO I=2,N1
	  DXT(I)=DXT(I-1)+(DDXT(I)+DDXT(I-1))/2.0D0*DTI
	ENDDO
	DO I=2,N1
	  XT(I)=XT(I-1)+(DXT(I)+DXT(I-1))/2.0D0*DTI
	ENDDO
C
C     Min/Max - Zeros search
C
	DO I=1,NPT2
	  DISMAX(I)=0.0D0
	  VELMAX(I)=0.0D0
	  ACCMAX(I)=0.0D0
	ENDDO
	DO I=1,NPT2
        TAMAX(I)=0.0D0
        TVMAX(I)=0.0D0
        TDMAX(I)=0.0D0
	ENDDO
	IDISMAX=0
	IACCMAX=0
	IVELMAX=0
c
	DO JJ=1,4
	  BB(1,JJ)=0.0D0
	  BB(2,JJ)=0.0D0
	  BB(3,JJ)=0.0D0
	ENDDO
	LL=N1 ! Need to check this later for finding min/max
c
c	Echo input data and calcs
c
	WRITE(*,*)'Finish Upsamle'
	WRITE(*,*) 'No. of time steps = ',LL
C	Unit 7 is AT
	WRITE(7,*) 'No. of time steps = ',LL
C	Unit 8 is DT
	WRITE(8,*) 'No. of time steps = ',LL
C	Unit 9 is TT
	WRITE(9,*) 'No. of time steps = ',LL
C	Unit 10 is AS
	WRITE(10,*) 'No. of time steps = ',LL
C	Unit 11 is DS
	WRITE(11,*) 'No. of time steps = ',LL
C	Unit 14 is VT
	WRITE(14,*) 'No. of time steps = ',LL
C	Unit 15 is VS
	WRITE(15,*) 'No. of time steps = ',LL
C
	IF (LL.GT.NPT .OR. LL.LT.1) THEN
		WRITE(*,*) 'Number of time steps too large, zero,
     &                or negative, set to ', NPT
	    LL = NPT
	ENDIF
C***
C     Start out the counters with the first time step
C       Assume start at zero
C***
      IDISMAX=1
      TDMAX(IDISMAX)=BB(1,4)
	DISMAX(IDISMAX)=BB(1,2)
c
      IVELMAX=1
      TVMAX(IVELMAX)=BB(1,4)
	VELMAX(IVELMAX)=BB(1,2)
c
      IACCMAX=1
      TAMAX(IACCMAX)=BB(1,4)
	ACCMAX(IACCMAX)=BB(1,3)
C***
C     Loop over entire time history
C       The two are different lengths
C***
	DO 10 M=1,NT ! Do Original time history first
C***
C     Store output info into BB array for MIN/MAX Processing
C***
	DO J=1,4
        BB(3,J)=BB(2,J)
        BB(2,J)=BB(1,J)
      ENDDO
C***
C     Original time history
C***
	BB(1,1)=XT(M)
	BB(1,2)=DXT(M)
	BB(1,3)=DDXT(M)
C
	BB(1,4)=DTI*(M-1) ! Time
C
	IF(NPSTEP.EQ.NPRINT) THEN
        WRITE(9,2) BB(1,4),BB(1,1),BB(1,2),BB(1,3)
	  NPSTEP=0
	ENDIF
C
	IF(M.LT.3) GOTO 10  ! Need at least 3 time steps in BB
C***
C     Pick out zeros, max and min of DIS  = BB(J,2)
C       J=1 IS CURRENT TIME STEP
C       J=2 IS PREVIOUS TIME STEP
C       J=3 IS PREVIOUS-PREVIOUS TIME STEP
C     Note: a zero can also be a local min or max
C***
      IF(BB(1,2).EQ.0D0) THEN  ! CHECK FOR Exact ZERO
        WRITE(8,2) BB(1,4),BB(1,1),BB(1,2),BB(1,3)
        IDISMAX=IDISMAX+1
        TDMAX(IDISMAX)=BB(1,2)
	  DISMAX(IDISMAX)=BB(1,2)
c	  If sign change then went through zero
      ELSE
	  ZTEST=ABS(BB(1,2))*BB(2,2)/BB(1,2)/ABS(BB(2,2))
	  IF(ZTEST.LT.0D0) THEN  ! CHECK FOR ZERO=CHANGE SIGN
          WRITE(8,2) BB(1,4),BB(1,1),BB(1,2),BB(1,3)
          IDISMAX=IDISMAX+1
          TDMAX(IDISMAX)=BB(2,4)
	    DISMAX(IDISMAX)=BB(2,2)
	  ENDIF
	ENDIF
C***
C     Check for LOCAL min and max
C       Check for slope change
C***
	IF(BB(1,2).LT.BB(2,2) .AND. BB(3,2).LT.BB(2,2)) THEN
        WRITE(8,2) BB(1,4),BB(1,1),BB(1,2),BB(1,3)
          IDISMAX=IDISMAX+1
          TDMAX(IDISMAX)=BB(2,4)
	    DISMAX(IDISMAX)=BB(2,2)
	ELSEIF(BB(1,2).GT.BB(2,2) .AND. BB(3,2).GT.BB(2,2)) THEN
        WRITE(8,2) BB(1,4),BB(1,1),BB(1,2),BB(1,3)
          IDISMAX=IDISMAX+1
          TDMAX(IDISMAX)=BB(2,4)
	    DISMAX(IDISMAX)=BB(2,2)
	ENDIF
C
	NPSTEP=NPSTEP+1
C
  10	CONTINUE   ! End of M-dt loop
C
C     Write out summary info
C
      WRITE(11,*) '    Vaules based on Displacement - 1 Cycle'
      WRITE(11,*) '   ID    TMAX     FREQ     DAMPING     DIS1     DIS2'
      DO ID=6,IDISMAX,2
        FREQ=1.0D0/(TDMAX(ID)-TDMAX(ID-4))
	  DMAX2=DABS(DISMAX(ID))
	  DMAX1=DABS(DISMAX(ID-4))
	  DENOM=DSQRT(1.0D0+(TWOPI/LOG(DMAX2/DMAX1))**2)
	  DAMP=1.0D0/DENOM
	  WRITE(11,'(3X,I3,F12.6,4F10.4)')
     &    ID,TDMAX(ID),FREQ,DAMP,DISMAX(ID),DISMAX(ID-4)
	ENDDO  
C
      WRITE(10,*) '      Vaules based on Acceleration - 1 Cycle'
      WRITE(10,*) '   ID     TMAX      FREQ      DAMPING'
     &,'     ACC1     ACC2'
      DO ID=5,IACCMAX,2
        FREQ=1.0D0/(TAMAX(ID)-TAMAX(ID-4))
	  AMAX2=DABS(ACCMAX(ID))
	  AMAX1=DABS(ACCMAX(ID-4))
	  DENOM=DSQRT(1.0D0+(TWOPI/LOG(AMAX2/AMAX1))**2)
	  DAMP=1.0D0/DENOM
	  WRITE(10,'(3X,I3,F12.6,4F10.4)')
     &    ID,TAMAX(ID),FREQ,DAMP,ACCMAX(ID),ACCMAX(ID-4)
	ENDDO  
C
      STOP
C
C	WRITE(9,2)TIME,RX1,DX1,AX1/G,RX2,DX2,AX2/G,UL,DUL
   1  FORMAT(/,5X,'Time',10X,'Disp_O',6X,'Vel_O',5X,'Acc_O'
     &                   ,8X,'Disp_U',6X,'Vel_U',5X,'Acc_U'
     &,/,5X,"(Sec)", 9X,"(In.)",4X,"(In./Sec)",3X,"(Gs)"
     &             ,10X,"(In.)",4X,"(In./Sec)",3X,"(Gs)")
   2	FORMAT(1PE14.6,1x,3E11.3,2(2X,3E11.3))

   3	FORMAT(1X)

C	WRITE(*,4)
C	WRITE(*,5)TONE,TMAX,UDD,VDD,FN1,FN2,TAU1,TAU2,DT
   4	FORMAT(6X,'Pulse',7X,'Time',8X,'Acc',7X,'Delta'
     &,/5X,'Length',7X,'Max',8X,'Peak',5X,'Velocity'
     &,8X,'F1c',6X,' F2c ',6X,' TAU1',8X,' TAU2',6X,'  DT')
   5	FORMAT(3X,1PE11.3,1x,3E11.3,2(2X,3E11.3),4(2X,2E11.3))

C	WRITE(*,8)XM1*G,XM2*G,XK1,XK2,FN1U,FN2U,XZI1,XZI2
   8	FORMAT(/4X,'     W1     ','     W2     ','    K1    '
     &          ,'     K2     ','    F1u     ','   F2u   '
     &          ,'    ETA1    ','     ETA2   '
     &,/3x,1PE11.3,1X,3E11.3,2X,3E11.3,2X,3E11.3)

  18	FORMAT(/ 4X,'   GAP_POS  ','  GAP_NEG  '
     &           ,'    K_POS   ','   K_NEG   '
     &           ,'   FN_POS   ','  FN_NEG   '
     &           ,'    C_POS   ','   C_NEG   '
     &           ,'    Q_POS   ','   Q_NEG   '
     &           ,'   XZ_POS   ','  XZ_NEG   '
     &,/3X,1PE11.3,1x,3E11.3,2(2X,3E11.3),5(2X,2E11.3))
  19	FORMAT(/4X,'    T1      ','   DDX(1)   '
     &          ,'    T2      ','  DDX(2)     '
     &          ,'    DTI'
     &,/3X,1PE11.3,1x,3E11.3,2(2X,3E11.3),4(2X,2E11.3))
C
	END
C *****************************************************************************
	SUBROUTINE UPSAMP(DX,DXT,N1,L)
C *****************************************************************************
C     Upsample sequence DXT[N1] by a factor of L to form sequence DX[L*N1]
C ***
	IMPLICIT REAL*8 (A-H,O-Z)
	DIMENSION DX(*),DXT(*)
	COMMON /CON/ONEPI,TWOPI,D2R
C
      NXTT = N1
      NXT = NXTT * L
C      ReDim DX(NXT)
      DX(1) = DXT(1)
      DO N = 2 , NXT
       DX(N) = 0.0D0
       DO K = 1 , NXTT
        If (N .EQ. K * L) Then ! Avoid singularity
          DX(N) = DXT(K)
        Else
          DX(N) = DX(N) + DXT(K) *
     &    Sin(ONEPI * (N - K * L) / L) / (ONEPI * (N - K * L) / L)
        End If
       ENDDO ! Next K
      ENDDO ! Next N
!    Fix up phase shift and delete last L-1 terms, keep first term
      DO N = 2 , NXT - L + 1
       DX(N) = DX(N + L - 1)
      ENDDO ! Next N
      NXT = NXT - L + 1
	N1=NXT
C
	RETURN
	END