c as copied from EDC - 1996 Jan 5 c SUBROUTINE DISTRIB(RES,M,DMEAN,RMS,NUSE,SKEW,RKURT,QMEAN,NQ, 1 PEAK,FWHM,NPEAK,LEHM,X,Y,XFIT,YFIT,IERR) C C TO EXAMINE DISTRIBUTION OF SLR DATA, AND DETERMINE MEAN, PEAK and C LEADING EDGE HALF MAX. Also determines RMS, full-width-half-max, C skewness and kurtosis of the distribution. C C WRITTEN BY A.T. SINCLAIR (RGO). C C C FOR A SATELLITE THE DATA ARE RESIDUALS OF THE MEASURED QUANTITY C FROM A SMOOTHING FUNCTION, SUCH THAT ALL TRENDS ARE REMOVED. C FOR A CALIBRATION TARGET THE RAW RANGE MEASURES CAN BE USED. C C Same reference point must be used for calibration data as is used C for a satellite (ie. mean, peak or leading-edge-half-max). C C A BAND OF DATA WIDER THAN +/- 3.0*RMS AROUND THE APPROX PEAK C OF THE DATA DISTRIBUTION SHOULD BE INPUT TO THE SUBROUTINE, SO C THAT NO CLIPPING OCCURS BEFORE ENTRY TO THIS ROUTINE. C C The functions of the routine are: C 1. Determine the mean and RMS-from-mean of the data, using an C iterated rejection level of n*RMS. The value of n is hard coded C at the start of the routine (parameter AN). A value of n=3.0 C should normally be used, but for systems using first-photon C detection it may be considered that a tighter value might be C more suitable (but not less than 2.5). The data within this band C of n*RMS are used for calculating the skewness and kurtosis. C C 2. Determine the mean using a very tight rejection of 1.0*RMS, where C RMS is now held fixed at the value determined in (1). This gives C an estimate of the peak of the data distribution. C C 3. Use a "smoothing" method to derive a distribution function for C the data. This method regards each data point as being the centre C of a Gaussian distribution with sigma=60 ps, and sums the C contributions of all these Gaussians at 201 points spread over a C width of +/- 1000 ps (15 cm 1-way) about the mean. The peak of C this distribution function is located. C C 4. Fit a quartic polynomial to the peak of the data distribution C (actually to the part greater than 40% of the maximum). The reason C is that sometimes if the input data are sparse the smoothing method C may not give a single peak. Increasing the smoothing sigma will C usually fix this, but will also move the estimated location of the C peak towards the direction of skewness of the data. Instead the C peak of the fitted polynomial is used. C C So the routine returns two estimates of the peak, from the 1.0*RMS C and from the smoothing/polynomial fit. These should be in close C agreement (except that if the data distribution is poor the 1.0*RMS C iteration can fail to converge). C C UNITS: THE UNITS OF THE INPUT DATA ARE TWO-WAY LIGHT TIME IN PICOSEC. C THE UNITS OF OUTPUT QUANTITIES, RMS, MEAN, PEAK AND FWHM C ARE ALSO TWO-WAY LIGHT TIME IN PICOSEC. C C ALL ARGUMENTS OF SUBROUTINE ARE DOUBLE PRECISION OR INTEGER. C C INPUT: C RES(5000) ARRAY CONTAINING THE DATA VALUES (UP TO 5000) C M NUMBER OF DATA VALUES C C OUTPUT: C DMEAN ARITHMETIC MEAN VALUE OF DATA WITHIN n*RMS. C RMS RMS OF DATA FROM THE MEAN USING ITERATED n*RMS REJECTION. C NUSE NUMBER OF POINTS WITHIN n*RMS. C SKEW SKEWNESS OF DATA WITHIN n*RMS (dimensionless). C (NEGATIVE IF SKEWED TOWARDS SHORT RANGES) C RKURT KURTOSIS OF DATA WITHIN n*RMS (dimensionless) C C (value of n is hard-coded at start of the subroutine =AN) C C QMEAN ARITHMETIC MEAN USING 1.0*RMS REJECTION. C NQ NUMBER OF POINTS RETAINED IN FORMING ITERATED C MEAN USING 1.0*RMS REJECTION (SAME RMS AS ABOVE). C PEAK PEAK VALUE OF DISTRIBUTION FUNCTION OF DATA C USING SMOOTHING METHOD. C FWHM FULL-WIDTH HALF-MAXIMUM OF DISTRIBUTION C (FOR GAUSSIAN DISTRIBUTION, FWHM = 2.355*SIGMA) C NPEAK NUMBER OF PEAKS WITHIN FWHM. C IF VALUE IS GREATER THAN 1 THEN THE PEAK OF THE C FITTED POLYNOMIAL IS USED INSTEAD OF THE PEAK OF C THE DISTRIBUTION FUNCTION. C LEHM (DOUBLE PRECISION). LOCATION OF LEADING EDGE HALF C MAXIMUM OF THE DISTRIBUTION. Note that the leading C edge has shorter range measurement values than the C trailing edge. C X(201) ARRAY OF ORDINATE VALUES OF DISTRIBUTION, USING C SMOOTHING METHOD. C Y(201) ARRAY OF ABCISSA VALUES OF DISTRIBUTION, WITH PEAK C VALUE SET ARBITRARILY TO 100. C (PLOT Y AGAINST X TO SEE DISTRIBUTION OF DATA) C XFIT(201) ARRAY OR ORDINATE VALUES OF POLYNOMIAL FITTED TO THE C PEAK OF THE DATA DISTRIBUTION. C YFIT(201) ARRAY OF ABCISSA VALUES OF THE FITTED POLYNOMIAL. C IERR ERROR RETURN: C IERR=0 : OK. C IERR=1 : TOO MANY DATA POINTS INPUT. C IERR=2 : DETERMINATION OF MEAN USING REJ=n*RMS C NOT CONVERGED. C IERR=3 : DETERMINATION OF MEAN USING REJ=1*RMS C NOT CONVERGED. C IERR=4 : FAILURE TO INVERT MATRIX IN POLYNOMIAL FIT. C IERR=5 : LESS THAN 5 POINTS REMAINING IN POLYNOMIAL FIT. C C Note. For errors IERR=1 or 2 the subroutine will return immediately C with no further computation. C For IERR=3 it will continue with calculation of other quantities, C but the value of the 1*RMS mean will be unreliable. C For IERR=4 or IERR=5 the polynomial fit to the distribution C function is not carried out, and the value of PEAK is not C replaced by the peak of the polynomial if there are more than C one peak to the distribution function. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION RES(5000),DMEAN,RMS,SKEW,RKURT,QMEAN, 1 PEAK,FWHM,LEHM,X(201),Y(201),XFIT(201),YFIT(201),C(5) INTEGER*4 M,NUSE,NPEAK,IERR REAL SX,F,FUNC C IF(M.GT.5000) THEN IERR=1 RETURN ENDIF C C SET REJECTION LEVEL n FOR SCREENING OF DATA. C SHOULD BE 3.0 FOR SYSTEMS OPERATING MULTI=PHOTON DETECTION DEVICES C (EG MCPs), AND POSSIBLY 2.5 FOR SYSTEMS OPERATING SINGLE PHOTON DETECTION C DEVICES( EG PHOTO-DIODES). C AN=2.5 C C SET OTHER PARAMETER VALUES (IN 2-WAY PICOSEC) : C SIG : THE SIGMA OF THE GAUSSIAN SMOOTHING FUNCTIONS. C WIDE : THE FULL WIDTH OVER WHICH THE DISTRIBUTION IS EVALUATED. C STEP : THE STEP AT WHICH THE DISTRIBUTION IS EVALUATED: C SIG=60.D0 WIDE=2000.D0 STEP=10.D0 NDIS=IDNINT(WIDE/STEP) +1 IF(NDIS.GT.201) NDIS=201 C IERR=0 TOL=STEP*1.D-2 C DMEAN=0.D0 PMEAN=1.D20 REJ=1.D20 C2=0.D0 C3=0.D0 C4=0.D0 C C FORM ITERATED MEAN WITH n*RMS REJECTION. C ON LAST ITERATION, FORM MOMENTS ABOUT MEAN. C ITER=20 IEND=0 DO 100 I=1,ITER IF(I.EQ.ITER) IEND=2 SR=0.D0 SRR=0.D0 NPT=0 C DO 110 J=1,M RJ=RES(J) DR=RJ -DMEAN IF(DABS(DR).GT.REJ) GO TO 108 NPT=NPT +1 SR=SR +DR SRR=SRR +DR*DR IF(IEND.EQ.0) GO TO 108 C2=C2 +DR*DR C3=C3 +DR*DR*DR C4=C4 +DR*DR*DR*DR 108 CONTINUE 110 CONTINUE C C DELTA IS CORRECTION TO PREVIOUS ESTIMATE OF THE MEAN. C RMS IS RMS FROM NEW ESTIMATE OF THE MEAN. C THEN, REPLACE DMEAN BY THE NEW ESTIMATE OF THE MEAN. C C IF(IEND.NE.0) GO TO 120 IF(NPT.EQ.0) STOP ' ALL POINTS REJECTED FOR n*RMS MEAN' ANPT=DBLE(NPT) DELTA=SR/ANPT RMS=DSQRT(SRR/ANPT -DELTA*DELTA) DMEAN=DMEAN +DELTA REJ=AN*RMS IF(DABS(DMEAN-PMEAN).LT.TOL) IEND=1 PMEAN=DMEAN 100 CONTINUE C C 120 CONTINUE IF(IEND.EQ.2) THEN IERR=2 RETURN ENDIF C2=C2/ANPT C3=C3/ANPT C4=C4/ANPT SKEW=C3/(C2**1.5D0) RKURT=C4/(C2*C2) NUSE=NPT C C FORM ITERATED MEAN WITH 1*RMS REJECTION. C NOTE. RMS IS HELD FIXED AT VALUE OBTAINED ABOVE. C PMEAN=1.D20 QMEAN=DMEAN IEND=0 REJ=1.D0*RMS C DO 121 I=1,ITER IF(I.EQ.ITER) IEND=3 SR=0.D0 NPT=0 C DO 122 J=1,M RJ=RES(J) DR=RJ -QMEAN IF(DABS(DR).GT.REJ) GO TO 122 NPT=NPT +1 SR=SR +DR 122 CONTINUE C IF(NPT.EQ.0) STOP ' ALL POINTS REJECTED FOR 1.0*RMS MEAN' ANPT=DBLE(NPT) DELTA=SR/ANPT QMEAN=QMEAN +DELTA IF(DABS(QMEAN-PMEAN).LT.TOL) GO TO 123 PMEAN=QMEAN 121 CONTINUE C 123 CONTINUE IF(IEND.EQ.3) IERR=3 NQ=NPT C C EVALUATE DISTRIBUTION FUNCTION AT REQUIRED STEPS C AND MAKE COARSE LOCATION OF PEAK. C L=(NDIS -1)/2 XI=DMEAN -DBLE(L+1)*STEP YMAX=0.D0 C DO 130 I=1,NDIS XI=XI +STEP X(I)=XI SUM=0.D0 C DO 140 J=1,M DX=(XI -RES(J))/SIG SX=SNGL(DX) F=FUNC(SX) SUM=SUM +DBLE(F) 140 CONTINUE C Y(I)=SUM IF(SUM.GT.YMAX) THEN IMAX=I YMAX=SUM ENDIF 130 CONTINUE C C MAKE FINE LOCATION OF PEAK C FSTEP=STEP/10.D0 XI=X(IMAX-1) YMAX=0.D0 C DO 150 I=1,19 XI=XI +FSTEP SUM=0.D0 DO 155 J=1,M DX=(XI -RES(J))/SIG SX=SNGL(DX) F=FUNC(SX) SUM=SUM +DBLE(F) 155 CONTINUE IF(SUM.GT.YMAX) THEN XMAX=XI YMAX=SUM ENDIF 150 CONTINUE C C REPLACE COARSE DETERMINATION OF PEAK WITH FINE DETERMINATION C X(IMAX)=XMAX PEAK=XMAX Y(IMAX)=YMAX C C NORMALISE DISTRIBUTION TO MAXIMUM VALUE 100. C SELECT POINTS WITH Y>40.0 FOR POLYNOMIAL FITTING. C NFIT=0 FAC=100.D0/YMAX DO 160 I=1,NDIS YC=Y(I)*FAC Y(I)=YC IF(YC.LT.40.D0) GO TO 160 NFIT=NFIT +1 XFIT(NFIT)=X(I) YFIT(NFIT)=YC 160 CONTINUE C C C DETERMINE FWHM, AND COUNT NUMBER OF PEAKS WITHIN FWHM. C ALSO DETERMINE LEHM (leading edge half maximum) C NPEAK=1 ISLOPE=+1 DO 170 I=IMAX-1,1,-1 I1=I I2=I1 +1 Y1=Y(I1) Y2=Y(I2) IF(Y1.GT.Y2) THEN IF(ISLOPE.EQ.1) NPEAK=NPEAK +1 ISLOPE=-1 ELSE ISLOPE=+1 ENDIF IF(Y1.LT.50.D0) GO TO 171 170 CONTINUE 171 CONTINUE P=(50.0 -Y1)/(Y2 -Y1) XLO=X(I1) +STEP*P C ISLOPE=-1 DO 180 I=IMAX+1,NDIS I1=I I2=I1 -1 Y1=Y(I1) Y2=Y(I2) IF(Y1.GT.Y2) THEN IF(ISLOPE.EQ.-1) NPEAK=NPEAK +1 ISLOPE=+1 ELSE ISLOPE=-1 ENDIF IF(Y1.LT.50.D0) GO TO 181 180 CONTINUE 181 CONTINUE P=(Y2 -50.0)/(Y2 -Y1) XHI=X(I2) +STEP*P C FWHM=XHI -XLO LEHM=XLO C C C FIT A POLYNOMIAL TO THE PEAK C XZERO=XMAX XSPAN=(XFIT(NFIT) -XFIT(1))*0.5D0 CALL POLY4(NFIT,XFIT,YFIT,XZERO,XSPAN,C,JERR) IF(JERR.NE.0) THEN IERR=JERR +3 RETURN ENDIF C C EVALUATE THE FITTED CURVE, AND LOCATE MAXIMUM C YMAX=0.D0 XMAX=0.D0 X1=XFIT(1) X2=XFIT(NFIT) DX=(X2 -X1)/200.D0 XC=X1 -DX C DO 350 I=1,201 XC=XC +DX XQ=(XC -XZERO)/XSPAN X2=XQ*XQ X3=X2*XQ X4=X3*XQ C1=C(1) C2=C(2) C3=C(3) C4=C(4) C5=C(5) YC=C1 +C2*XQ +C3*X2 +C4*X3 +C5*X4 XFIT(I)=XC YFIT(I)=YC IF(YC.GT.YMAX) THEN YMAX=YC XMAX=XC ENDIF 350 CONTINUE C C C IMPROVE DETERMINATION OF THE MAX OF THE FITTED POLYNOMIAL BY FINDING C ZERO OF DERIVATIVE (NEWTON-RAPHSON ITERATION, 4 ITERATIONS) C XQ=(XMAX -XZERO)/XSPAN DO 230 I=1,4 X2=XQ*XQ X3=X2*XQ X4=X3*XQ F =C2 +2.0*C3*XQ +3.0*C4*X2 +4.0*C5*X3 FD=2.0*C3 +6.0*C4*XQ +12.0*C5*X2 XQ=XQ -F/FD 230 CONTINUE XMAX=XQ*XSPAN +XZERO C C IF DISTRIBUTION FUNCTION HAS MORE THAN ONE PEAK THEN USE THE PEAK C OF THE FITTED POLYNOMIAL AS THE SELECTED ESTIMATE INSTEAD. C NOTE: C This will not be done if IERR=4 or 5, as the subroutine will C already have returned before this point. C IF(NPEAK.GT.1) PEAK=XMAX C RETURN END C C ************************************************************* C REAL FUNCTION FUNC(X) REAL X,X2 X2=X*X FUNC=EXP(-0.5*X2) RETURN END C C ******************** C SUBROUTINE POLY4(NP,X,Y,XZERO,XSPAN,C,IERR) C C FITS POLYNOMIAL OF ORDER 4. C C INPUT: C NP - NUMBER OF POINTS C X - ARRAY (DIM 201) CONTAINING X COORDS OF POINTS C Y - ARRAY (DIM 201) CONTAINING Y COORDS OF POINTS C XZERO- ZERO OF X COORDS FOR THE POLYNOMIAL. C XSPAN- SPAN OF THE IMPUT X VALUES. C C OUTPUT: C C - ARRAY (DIM 5) CONTAINING POLYNOMIAL COEFFICIENTS C IERR - ERROR INDICATOR C IERR=0 : OK C IERR=1 : FAILURE TO INVERT MATRIX C IERR=2 : LESS THAN 5 POINTS DURING ITERATION C C THE OUTPUT POLYNOMIAL IS:- C XX = (X(I) - XZERO)/XSPAN C Y(I) = C(1) + C(2)*XX + C(3)*XX**2 + C(4)*XX**3 + C(5)*XX**4 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION X(201),Y(201),XZERO,XSPAN,C(5) DIMENSION A(5,5),B(5),XP(5),DC(5),SE(5) C IORD=4 ITER=10 RMSB=0.D0 RMSA=0.D0 IC=IORD +1 REJ=1.D20 DO 201 I=1,5 SE(I)=0.D0 201 C(I)=0.D0 C C DO 300 L=1,ITER C C Set anything that is a running total to be zero C FF = 0.D0 SXX = 0.D0 M=0 DO 200 I=1,IC B(I)=0.D0 DC(I)=0.D0 DO 200 J=1,IC 200 A(J,I)=0.D0 SIGMA=0.D0 C DO 100 I=1,NP XX=(X(I)-XZERO)/XSPAN XP(1)=1.D0 R=C(1) DO 101 J=1,IORD XP(J+1)=XP(J)*XX R=R +C(J+1)*XP(J+1) 101 CONTINUE F=Y(I) -R C IF(DABS(F).GT.REJ) GO TO 100 M=M +1 FF=FF +F*F SXX=SXX +XX*XX DO 102 J=1,IC B(J)=B(J) +F*XP(J) DO 102 K=1,IC A(K,J)=A(K,J) +XP(K)*XP(J) 102 CONTINUE 100 CONTINUE C C C IF(M.EQ.0)THEN IERR=2 RETURN ENDIF AM=DBLE(FLOAT(M)) RMSB=DSQRT(FF/AM) IF(M.LE.IC)THEN IERR=2 RETURN ENDIF C CALL DCHOLS(A,IC,IERR) IF(IERR.EQ.1) RETURN C VV=FF DO 103 J=1,IC DO 104 K=1,IC 104 DC(J)=DC(J) +A(J,K)*B(K) C(J)=C(J) +DC(J) VV=VV -DC(J)*B(J) 103 CONTINUE C C IF(VV.LT.0.D0) VV=0.D0 RMSA=DSQRT(VV/AM) REJ=3.0D0*RMSB FAC = DSQRT(VV/DBLE(M - IC)) DO 105 J=1,IC SE(J)=DSQRT(A(J,J))*FAC 105 CONTINUE 300 CONTINUE C C RETURN END C C ********************************* C SUBROUTINE DCHOLS(A,M,IERR) C C DOUBLE PRECISION CHOLESKI MATRIX INVERSION C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(5,5),B(5,5),X(5),S(5) IERR=0 C C SCALE MATRIX C DO 300 I=1,M SUM=A(I,I) IF(SUM.LE.0.D0) GO TO 400 300 S(I)=1.D0/DSQRT(A(I,I)) DO 301 I=1,M DO 301 J=1,M A(I,J)=A(I,J)*S(I)*S(J) 301 CONTINUE C DO 100 I=1,M SUM=A(I,I) IF(I.EQ.1) GO TO 10 DO 101 K=1,I-1 SUM=SUM -B(K,I)**2 101 CONTINUE 10 IF(SUM.LE.0.D0) GO TO 400 B(I,I)=DSQRT(SUM) IF(I.EQ.M) GO TO 100 DO 102 J=I+1,M SUM=A(J,I) IF(I.EQ.1) GO TO 11 DO 103 K=1,I-1 SUM=SUM -B(K,I)*B(K,J) 103 CONTINUE 11 B(I,J)=SUM/B(I,I) 102 CONTINUE 100 CONTINUE C C DO 200 I=1,M DO 206 J=1,M 206 X(J)=0.D0 X(I)=1.D0 DO 202 J=1,M SUM=X(J) IF(J.EQ.1) GO TO 20 DO 203 K=1,J-1 SUM=SUM -B(K,J)*X(K) 203 CONTINUE 20 X(J)=SUM/B(J,J) 202 CONTINUE DO 204 J=1,M M1=M -J +1 SUM=X(M1) IF(J.EQ.1) GO TO 21 DO 205 K=1,J-1 M2=M -K +1 SUM=SUM -B(M1,M2)*X(M2) 205 CONTINUE 21 CONTINUE X(M1)=SUM/B(M1,M1) 204 CONTINUE DO 22 J=1,M 22 A(J,I)=X(J) 200 CONTINUE C C UNSCALE MATRIX C DO 302 I=1,M DO 302 J=1,M A(I,J)=A(I,J)*S(I)*S(J) 302 CONTINUE RETURN 400 CONTINUE IERR=1 RETURN END