$STORAGE:2
C  REV    280694 IP  SPAD DATA VERSION / NEW FILTERING ALGORITHM
C  REV    060993 MC  !!! MUST BE LINKED WITH LIBRARY RADAR.LIB
C      NEW DATA FORMAT SINCE AUG.88
C      NEW DATA FORMAT SINCE SEP.94
C  REV   100890 MC  WRITES THE RMS IN PSEC TO VALUE IPRES FOR CATALOG
C                   PROGRAM MUST COPY THE 2. LINE
C
C     INPUT FILE:   RANGXXX.DAT
C
C     OUPUT FILE:   RANGXXX.DAT
C
C     PROGRAM FOR MEASURED DATA ANALYSIS BY MEANS OF
C     POLYNOMIAL APPROXIMATION OF MEASURED DATA
C     AUTOMATICAL VERSION
C     STORES THE RESIDUALS IN THE SAME FILE
C
      PROGRAM RDFP
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION NAMS(3),IPLOT(42),DELTA(2000)
      INTEGER HISTO(50),RLAB
      DOUBLE PRECISION X(2000),Y(2000),CP(20)
      CHARACTER*64 NAME
      DOUBLE PRECISION DINP(3)
      DIMENSION IW(2000),JBUF(7)
      DATA JBUF/2H--,2HRM,2HTG,2HN4,2HSC,2HMD,2HBR/
      DIFF=0.
      RMSF=1.D12
      RMS=1.D11
      JUMP=1
      MAXDEG=2
      ISCC=0
      SCC=1.D12
      OPEN(UNIT=1,FILE='CON',STATUS='UNKNOWN')
C
C     INPUT OF INITIAL VALUES
C
      WRITE(*,1)
1     FORMAT(20X/' RANGING DATA FITTING PROGRAM '
     #    /      ' **************************** '//
     # ' ENTER THE NO. OF RANGING DATA FILE >>> '\)
      READ(*,*) INUM
      WRITE(NAME,2) INUM
2     FORMAT('RANG',I3.3,'.DAT')
C     WRITE(*,950)
C950  FORMAT(' First point: ',\)
C     READ(*,*) IFIRST
C     WRITE(*,951)
C951  FORMAT(' Last point : '\)
C     READ(*,*) ILAST
      IFIRST=1
      ILAST=32111
C
C  WHICH COLLUMN IS PROCESSING 1-RANGE, 2-ONLINE DIFF, 3-OFP DIFF
C
      MODE=3
C
      WRITE(*,952)
952   FORMAT(' Starting polynom degree: '\)
      READ(*,*) IPOLY
      IPOLY=IPOLY-1
      OPEN(UNIT=2,FILE=NAME,STATUS='OLD')
      OPEN(UNIT=100,FILE='RANGWR.DAT',STATUS='UNKNOWN')
C
      READ(2,786)IYR,MO,IDY,SNAME
786   FORMAT(3I3,2X,F8.0)
      READ(2,787)ITEMP,IHUM,IPRES,C1,CC,LAB2,RLAB,MLAB
787   FORMAT(I4,I3,I5,F10.2,F8.0,3I2)
      IF( MLAB .EQ. 1 ) THEN
          WRITE(*,*) ' The file has been processed!'
          WRITE(*,*) ' RMS will not be written to the file!'
      ENDIF
      IB0=0
      IXY=0
C
3     CONTINUE
      IB0=IB0+1
      IPO=IPO+1
      READ(2,4,END=10)NPOINT,IHR,MIN,SEC,DINP(1),DINP(2),DINP(3),
     #         POLYDF,NDELTA,IVAHA
C
4     FORMAT(I5,I3,I3,F11.7,F13.8,F7.3,F7.2,F6.2,I4,I2)
C
C     SELECT POINTS
C
      IF(NPOINT.LT.IFIRST) GO TO 3
      IF(IVAHA.NE.1.) GO TO 3
      IF(NPOINT.LT.0.) GO TO 3
      IF(NPOINT.EQ.0) GO TO 10
      IF(NPOINT.GT.ILAST) GOTO 10
201   IXY=IXY+1
      IF(IXY.GT.2000) STOP 2000
      X(IXY)=IHR*1.D0+MIN/60.D0+SEC/3600.D0
      Y(IXY)=DINP(MODE)
      IW(IXY)=NPOINT
      IF(IXY.EQ.1.) GO TO 3
      IF(X(IXY).LT.X(IXY-1)) WRITE(*,203) NPOINT
203   FORMAT(' DATA SEQUENCE ERROR ON PT#: ',I3)
      GO TO 3
10    CONTINUE
      TIMEST=X(1)
      DO 461 ISH=1,IXY
461   X(ISH)=X(ISH)-TIMEST
C
      N=IXY
      IF(N.GT.2.) GO TO 7
      WRITE(*,8) N
8     FORMAT(// 10X,' WARNING '/' TOTAL NO.OF RANGING IS :',I2 )
      STOP
C
C     POLYNOMIAL FIT
C
7     CONTINUE
      IPOLY=IPOLY+1
300   CONTINUE
      ISCC=0
      DO 310 I=3,N
      TG=X(I)-X(I-1)
      IF(TG.GT.TGAP) TGAP=TG
310   CONTINUE
      RMD=(X(N)-X(3))/TGAP
      MAXDEG=RMD
      M=IPOLY
      IFLAG=0
C
C     STOP THE IMPROVING
C
C     THE MAX.DEGREE DUE TO THE TIME GAP HAS BEEN COMENTED OUT 26.6.93 IP
C
      IF(RMS.GT.RMSF)  JUMP=2
C     IF(M.EQ.MAXDEG)  JUMP=3
      IF(M.EQ.N/4)     JUMP=4
      IF(SCC.LT.0.2)   JUMP=5
      IF(M.GT.9)       JUMP=6
      IF(IFBRK(IDUMMY).EQ.-1) JUMP=7
      RMSF=RMS
      IF(JUMP.GT.1) WRITE(*,3100) JBUF(JUMP)
3100  FORMAT(' SEQUENCE STOPPED BY CRITERION   ',1A2/)
      IF(JUMP.GT.1) THEN
       WRITE(*,3101)
3101   FORMAT(' Poly degree: ',\)
       READ(*,*) M
      ENDIF
      CALL DAPRX(X,Y,N,M,CP)
      GO TO (302,303,303,303,303,303,303),JUMP
C
C
303   CONTINUE
C
302   CONTINUE
      DO 12 I=1,N
      XX=X(I)
      CALL DPVAL(XX,CP,M,D)
      D=Y(I)-D
      D=D*1.D6
      RMS=RMS+D*D
C
      GO TO (304,3033,3033,3033,3033,3033,3033),JUMP
C
C
3033  CONTINUE
      XX=X(I)+TIMEST
      IHR=XX
      MIN=(XX-IHR)*60.D0
      SEC=(XX-IHR*1.D0-MIN/60.D0)*3600.D0
      ISEC=SEC+0.5D0
C
304   IF(MODE.EQ.2.) D=D*1.D-3
      IF(MODE.EQ.3.) D=D*1.D-6
      DELTA(I)=D
C
C     CALCULATE THE S.C.C.
C
      SIGD=D*DIFF
      IF(SIGD.GT.0.) ISCC=ISCC+1
      DIFF=D
C
12    CONTINUE
C
      FD=N-M
      RMS=DSQRT(RMS/FD)
      IPOLY=M-1
      IF(MODE.EQ.2.) RMS=RMS*0.001
      IF(MODE.EQ.3.) RMS=RMS*1.D-6
      SCC=ISCC
      SCC=SCC/N
      WRITE(1,306)N ,IPOLY,RMS,SCC
306   FORMAT(' # OF POINTS ',I4,' POLDEG ',I2,' RMS ',
     # F6.3,'   NS    SCC',F6.3)
C
C     REJECT POINTS FOR POLY.CONSTRUCTION
C
      GO TO (311,301,301,301,301,301,301),JUMP
311   RMS3=RMS*2.
      DO 200 I=1,N
      IF(DABS(DELTA(I)).LT.RMS3) GO TO 200
C      WRITE(*,2010) IW(I)
C2010  FORMAT(' DELETED POINT NO. ',I5)
      IREJ=IREJ+1
      DO 202 II=I,N-1
      IW(II)=IW(II+1)
      X(II)=X(II+1)
      Y(II)=Y(II+1)
202   DELTA(II)=DELTA(II+1)
      N=N-1
      IFLAG=1
200   CONTINUE
C
      GO TO (7,300),IFLAG+1
C     END OF THE IMPROVE LOOP
C
301   WRITE(1,37)
37    FORMAT(///)
C
36    RMSM=RMS*150.
      WRITE(1,16)N,IREJ,IPOLY,SCC,RMS,RMSM
16    FORMAT(
     #/,' POINTS ACCEPTED ',I6,' REJECTED ',I3,
     #/,' POLYNOM DEGREE  ',I6,
     #/,' SERIAL COR.COEF.',F6.3,
     #/,' RMS            : ',F6.3,' (ns)',F7.1,'(mm)')
C
C     STORE THE RESIDUALS ON THE FILE RANGWR.DAT
C
C
      REWIND 2
      READ(2,786) IYR,MO,IDY,SNAME
      WRITE(100,786) IYR,MO,IDY,SNAME
      READ(2,787) ITEMP,IHUM,IPRES,C1,CC,LAB2,RLAB,MLAB
      IF( MLAB .NE. 1 )  THEN
            IPRES = RMS*1000.0
            IF( IPRES .GT. 9999 )   IPRES = 9999
      ENDIF
      WRITE(100,787) ITEMP,IHUM,IPRES,C1,CC,LAB2,RLAB,MLAB
C
      IPT=0
C     DELETE CRITERIUM FOR WRITTING TO THE OUTPUT FILE
      RMS3=RMS*2.0
500   CONTINUE
      READ(2,4,END=1000) NPOINT,IHR,MIN,SEC,DINP,POLYDF,NDELTA,IVAHA
      IF(NPOINT.LT.IFIRST) GO TO 500
      IF(NPOINT.GT.ILAST) GO TO 1000
      XX=IHR*1.D0+MIN/60.D0+SEC/3600.D0
      XX=XX-TIMEST
      CALL DPVAL(XX,CP,M,D)
      D=DINP(3)-D
      POLYDF=D
      IVAHA=1
      IF(DABS(D).GT.RMS3) IVAHA=0
      IF(POLYDF.GT.99.) POLYDF=99.9
      IF(POLYDF.LT.-99.) POLYDF=-99.9
      WRITE(100,4) NPOINT,IHR,MIN,SEC,DINP,POLYDF,NDELTA,IVAHA
      IF( IVAHA.EQ.1 )  IPT=IPT+1
      GO TO 500
C
C     COPY THE SCRATCH FILE TO THE ORIGINAL
C
1000  WRITE(*,555)
555   FORMAT(' Store to file   1-YES: ',\)
      READ(*,*) KODE
      CLOSE(UNIT=2)
      CLOSE(UNIT=100)
      IF(KODE.EQ.1) THEN
C
C         DELETE RANGXX.DAT AND RENAME RANGWR.DAT TO RANGXX.DAT
          CALL DELETE(NAME)
          CALL RENAME('RANGWR.DAT',NAME)
          WRITE(*,1111) IPT,NAME
1111      FORMAT(/,2X,I5,' valid points was written to file ',A12,/)
      ELSE
          CALL DELETE('RANGWR.DAT')
      ENDIF
C
      STOP
      END
C
C
      SUBROUTINE DAPRX(X,Y,N,M,CP)
C     ****************
C
C     APROXIMATION OF THE MEASURED VALUES BY THE POLYNOM
C
C     PARAMETERS    X,Y   ... INPUT VALUES IN D.P. ARRAYS
C                   N     ... NO. OF POINTS
C                   M      ... DEGREE OF POLYNOM
C                   CP    ... OUTPUT ARRAY CONTAINING THE
C                              COEFFICIENTS
C
      DOUBLE PRECISION X(1),Y(1),C(20,11),CP(20),B(40),DPV
      DOUBLE PRECISION U,S,P
      IF((N-M).GT.1) GO TO 67
      WRITE(*,66)
66    FORMAT(' *** ERROR *** IN APPROXIMATION , M >N-1 !!! ',/)
      RETURN
67    CONTINUE
      M1=M
      M=M+1
      M2=2*M1
      DO 490 K=1,M2
      B(K)=0.D0
      DO 490 I=1,N
      B(K)=B(K)+X(I)**K
490   CONTINUE
C
      DO 570 I=1,M
      DO 570 J=1,M
      I2=I+J-2
      IF(I2.NE.0) GO TO 560
      C(1,1)=N
      GO TO 570
560   C(I,J)=B(I2)
570   CONTINUE
C
      C(1,M+1)=0.D0
      DO 630 I=1,N
630   C(1,M+1)=C(1,M+1)+Y(I)
      DO 690 K=2,M
      C(K,M+1)=0.D0
      DO 690 I=1,N
      C(K,M+1)=C(K,M+1)+Y(I)*(X(I)**(K-1))
690   CONTINUE
C
      N1=M+1
      J=1
      K=1
      I=1
740   I3=I+K
      DPV=C(I,J)*C(I3,J)
      IF(DPV.NE.0.D0) GO TO 790
      U=1.D0
      S=0.0D0
      GO TO 820
790   P=DSQRT(C(I,J)*C(I,J)+C(I3,J)*C(I3,J))
      U=C(I,J)/P
      S=-1.D0*C(I3,J)/P
820   DO 860 J1=J,N1
      P=U*C(I,J1)-S*C(I3,J1)
      C(I3,J1)=S*C(I,J1)+U*C(I3,J1)
      C(I,J1)=P
860   CONTINUE
      K=K+1
      KH=K+I-M
      IF(KH.LE.0) GO TO 740
      K=1
      I=I+1
      J=J+1
      IF((J-M).LT.0) GO TO 740
C
      I=M
950   C(I,N1)=C(I,N1)/C(I,I)
      K=I-1
970   IF(K.LE.0) GO TO 1010
      C(K,N1)=C(K,N1)-C(K,I)*C(I,N1)
      K=K-1
      GO TO 970
1010  I=I-1
      IF((I-1).GE.0) GO TO 950
C
      DO 1011 I=1,20
1011  CP(I)=C(I,M+1)
      RETURN
      END
C
      SUBROUTINE DPVAL(X,CP,M,Y)
C     ****************
C
C     ROUTINEWE FOR THE VALUE OF THE POLYNOM EVALUATION
C     THE HURWITZ SCHEME IS USED
C
      DOUBLE PRECISION X,CP(10) ,HX ,Y
      IF(M-1)2,3,4
4     Y=CP(1)
      HX=1.D0
      DO 1 I=1,M-1
      HX=HX*X
1     Y=Y+CP(I+1)*HX
      RETURN
3     Y=CP(1)
      RETURN
2     Y=0.D0
      RETURN
      END
C
C
C
      SUBROUTINE HISTP(IHISTO,N)
C     ****************
C
C     PLOT THE HISTOGRAM ON THE LU
C     INPUT     IHISTO ... LINEAR VECTOR CONTAINING THE HIST.HEIGHTS
C               N      ... DIMENSION OF IHISTO
C
      DIMENSION IHISTO(1),IPLOT(50)
      IHM=0
      DO 1 I=1,N
      IF(IHM-IHISTO(I)) 2,2,1
2     IHM=IHISTO(I)
1     CONTINUE
C
      DO 4 I=1,IHM
      DO 3 IS=1,N
3     IPLOT(IS)=32+32*256
      DO 7 J=1,N
      IF(IHISTO(J)-IHM+I)7,7,5
5     IPLOT(J)=2H**
7     CONTINUE
4     WRITE(1,6)(IPLOT(JJ),JJ=1,N)
6     FORMAT(34X,50A2)
      RETURN
      END
