* SPLINES file SUBROUTINE USPC(Y,A,N) * * Calculation of B-spline coefficients from equidistant function * values. * - - * CALL USPC(Y,A,N) * - * The B-spline coefficients A(1)...A(N) are calculated from N given * equidistant function values Y(1)...Y(N). The not-a-knot condition * is assumed at both sides (i.e. no discontinuity in third * derivative at first and last inner knot). * Arrays Y(*) and A(*) may be identical. * For N > 100 the user has to include a common /SPCOM/ with at * least 2*N+6 words. * REAL Y(N),A(N),H(2,102) * ... DO I=2,N H(1,I)=Y(I)-Y(I-1) END DO H(2,1)=(5.0*H(1,2)+H(1,3))/2.0 H(2,2)=-H(2,1)+3.0*(H(1,3)+H(1,2)) H(1,2)=2.0 DO I=3,N-1 H(2,I)=-H(2,I-1)/H(1,I-1)+3.0*(H(1,I+1)+H(1,I)) H(1,I)=4.0-1.0/H(1,I-1) END DO H(2,N)=(5.0*H(1,N)+Y(N-1)-Y(N-2))/2.0 H(1,N)=1.0-2.0/H(1,N-1) * backward ... H(2,N)=(-2.0*H(2,N-1)/H(1,N-1)+H(2,N))/H(1,N) DO I=N-1,2,-1 H(2,I)=(H(2,I)-H(2,I+1))/H(1,I) END DO H(2,1)=(H(2,1)-2.0*H(2,2)) * ... and forward solution ALN=Y(N)+Y(N)-Y(N-1)-(2.0*H(2,N)+H(2,N-1))/3.0 DO I=1,N-1 A(I)=Y(I)+Y(I)-Y(I+1)+(2.0*H(2,I)+H(2,I+1))/3.0 END DO A(N)=ALN END FUNCTION USPD(X,A,NSP,XA,XB,NDER) * * calculate function value, derivative or integral at x * * * FD= USPD(X,A,N,XA,XB,ND) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * = function value at x for NDER =0 * = ND'th derivative for NDER =1, 2 or 3 * = integral between XA and x for NDER =-1 *) * * The result of the call is for all cases a linear combination of * B-spline coefficients. After the call the weights are stored in * the COMMON/SPCOM/I0,N0,S(200) * For ND=0...3 the returned value is: * sum(i=1...4) S(i)*A(i0+i) * for ND= -1 the returned value is: * sum(i=1...n0) S(i)*A(i) * * *) if, for nd=-1, n is larger than 200, the user has to include a * COMMON/SPCOM/I0,N0,S(N). * REAL*8 XX,HX,DX,U,V,DNSP,FAC,SUM REAL A(NSP) COMMON/SPCOM/I0,N0,S(200) * ... USPF=0.0 NDER=ND I0=-1 IF(NDER+1.EQ.0) THEN * integral USPF=USPI(XA,X,A,NSP,XA,XB) GOTO 100 ELSE IF(NDER.GT.3.OR.NDER+1.LT.0) THEN GOTO 100 END IF XX=AMIN1(AMAX1(X,XA),XB) DNSP=NSP-1 HX=(XB-XA)/DNSP DX=1.0D0+(XX-XA)/HX I0=DX IF(I0.GE.NSP) I0=NSP-1 V=DX-DFLOAT(I0) U=1.0D0-V I0=I0-2 IF(NDER.EQ.0) THEN * function value FAC=1.0D0/6.0D0 S(1)= (U*U*U) *FAC S(2)= (1.0+3.0*(1.0+U*V)*U) *FAC S(3)= (1.0+3.0*(1.0+U*V)*V) *FAC S(4)= (V*V*V) *FAC ELSE IF(NDER.EQ.1) THEN FAC=0.5D0/HX * first derivative S(1)=-(U*U) *FAC S(2)=-(1.0+(2.0*V-U)*U) *FAC S(3)= (1.0+(2.0*U-V)*V) *FAC S(4)= (V*V) *FAC ELSE IF(NDER.EQ.2) THEN FAC=1.0D0/(HX*HX) * second derivative S(1)= (U) *FAC S(2)= (V-2.0*U) *FAC S(3)= (U-2.0*V) *FAC S(4)= (V) *FAC ELSE IF(NDER.EQ.3) THEN * third derivative FAC=1.0D0/(HX*HX*HX) S(1)=-FAC S(2)= FAC+FAC+FAC S(3)=-S(2) S(4)=-S(1) END IF * apply not-a-knot condition at endpoints IF(I0.LT.0) THEN SS=S(1) S(1)=S(2)+4.0*SS S(2)=S(3)-6.0*SS S(3)=S(4)+4.0*SS S(4)=-SS I0=0 ELSE IF(I0.EQ.NSP-3) THEN SS=S(4) S(4)=S(3)+4.0*SS S(3)=S(2)-6.0*SS S(2)=S(1)+4.0*SS S(1)=-SS I0=NSP-4 END IF * final sum SUM =A(I0+1)*S(1)+A(I0+2)*S(2)+A(I0+3)*S(3)+A(I0+4)*S(4) USPD=SUM 100 RETURN END FUNCTION USPI(X1,X2,A,N,XA,XB) * * Integral of spline function between limits X1 and X2 * * FX= USPI(X1,X2,A,N,XA,XB) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * The integral over the spline function is calculated between the * limits X1 and X2. A limit below XA (above XB) is replaced by XA * (XB). The result is a linear combination of B-spline coefficients. * integral = sum(i=1...n0) S(i)*A(i) * After the call the weights are stored in the * COMMON/SPCOM/I0,N0,S(200). * If N is larger than 200, the user has to include a common * COMMON/SPCOM/I0,N0,S(N). * REAL A(N) COMMON/SPCOM/I0,N0,S(200) REAL BT(5) /16.0,17.0,28.0,23.0,24.0/,T(4) REAL TB(4,5)/15.0, 5.0, 5.0,-1.0, 1 16.0,16.0, 0.0, 0.0, 2 27.0,11.0, 1.0, 0.0, 3 22.0,12.0, 1.0, 0.0, 4 23.0,12.0, 1.0, 0.0/ DOUBLE PRECISION SUM * ... * clear array of weights DO I=1,N S(I)=0.0 END DO N0=0 HX=(XB-XA)/FLOAT(N-1) * lower limit of integration XJ=X1 DO 30 J=1,2 * transform to coordinate v relative to knots XJ=1.0+(XJ-XA)/HX IF(XJ.LT.1.0) GOTO 30 I0=XJ IF(I0.LE.N-1) THEN V=XJ-FLOAT(I0) ELSE I0=N-1 V =1.0 END IF U=1.0-V I0=I0-2 * calculate contribution from last knot interval T(1)=1.0-U*U*U*U T(2)=11.0-(((-3.0*U+4.0)*U+6.0)*U+4.0)*U T(3)=(((-3.0*V+4.0)*V+6.0)*V+4.0)*V T(4)=V*V*V*V * correct first and last interval and add contributions from * previous knot intervals IF(I0.LT.0) THEN TT=T(1) T(1)=T(2)+4.0*TT T(2)=T(3)-6.0*TT T(3)=T(4)+4.0*TT T(4)=-TT I0=0 ELSE IF(I0.LT.N-3) THEN T(1)=T(1)+TB(1,MIN0(I0+1,5)) T(2)=T(2)+TB(2,MIN0(I0+1,5)) T(3)=T(3)+TB(3,MIN0(I0+1,5)) T(4)=T(4)+TB(4,MIN0(I0+1,5)) ELSE IF(I0.EQ.N-3) THEN TT=T(4) T(4)=T(3)+4.0*TT+ 1.0 T(3)=T(2)-6.0*TT+12.0 T(2)=T(1)+4.0*TT+23.0 T(1)=-TT +24.0 I0=N-4 END IF * Subtract or add to array S(i) DO I=1,I0+4 IF(I.LE.I0) THEN B=BT(MIN0(I,5)) ELSE B=T(I-I0) END IF IF(J.EQ.1) THEN S(I)=-B ELSE S(I)=S(I)+B END IF N0=MAX0(N0,I0+4) END DO * Upper limit of integration 30 XJ=X2 * Integral by sum SUM=0.0 DO I=1,N IF(S(I).NE.0.0) THEN SUM=SUM+A(I)*S(I) S(I)=S(I)*HX/24.0 END IF END DO USPI=SUM*HX/24.0 100 RETURN END SUBROUTINE USPR(A,N,XA,XB) * * Print table of values of spline function, derivative, integral. * * CALL USPR(A,N,XA,XB) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * If N is larger than 200, the user has to include a common * COMMON/SPCOM/I0,N0,S(N) * REAL A(N),F(0:4) * ... WRITE(*,101) DO I=1,2*N-1 X=(FLOAT(2*N-1-I)*XA+FLOAT(I-1)*XB)/FLOAT(2*N-2) F( 0)=USPF(X,A,N,XA,XB) DO ND=1,3 F(ND)=USPD(X,A,N,XA,XB,ND) END DO F( 4)=USPI(XA,X,A,N,XA,XB) IF(MOD(I,2).EQ.1) THEN J=1+I/2 WRITE(*,102) J,A(J),X,F ELSE WRITE(*,103) X,F END IF END DO 100 RETURN 101 FORMAT('0',3X,'I',8X,'Coeff',10X,'X',10X,'F',10X, 1 2HF',10X,3HF'',10X,4HF''',10X,'Int F'/) 102 FORMAT(1X,I4,2X,3G15.5,5X,3G13.5,5X,G15.5) 103 FORMAT(22X,2G15.5,5X,3G13.5,5X,G15.5) END FUNCTION USPX(XL,FZ,A,N,XA,XB,NDER) * * Calculate next x-value with given function value or derivative * * -- -- - - -- -- ---- * X = USPX(XL,FZ,A,N,XA,XB,NDER) NDER = -1, 0, 1, 2 * * = next abscissa value with value = fd with x ge xa. * X > XB returned, if no (further) value can be found. * * The subprogram is not optimized for speed. * REAL A(N) * ... HX=(XB-XA)/FLOAT(N-1) EPS=0.1*HX USPX=XB+HX IF(XL.GT.XB) GOTO 100 IF(XL.LE.XA) THEN XX=XA ELSE XX=XL+EPS END IF * loop to find values with opposite sign X1=XX F1=USPD(X1,A,N,XA,XB,NDER)-FZ 10 X2=AMIN1(XB,X1+3.0*EPS) F2=USPD(X2,A,N,XA,XB,NDER)-FZ IF(F1*F2.LE.0.0) GOTO 20 IF(X2.EQ.XB) GOTO 100 X1=X2 F1=F2 GOTO 10 * enclosed 20 IF(F1.EQ.F2) GOTO 100 XX=(X1*F2-X2*F1)/(F2-F1) * Newton iteration to improve value IT=0 30 IF(IT.GT.10) GOTO 100 FN=USPD(XX,A,N,XA,XB,NDER)-FZ FD=USPD(XX,A,N,XA,XB,NDER+1) DX=-FN/FD XX=XX+DX IT=IT+1 IF(ABS(DX).GT.1.0E-3*EPS) GOTO 30 USPX=XX 100 RETURN END FUNCTION USPF(X,E,N,XA,XB) * * Compute function value or derivative for B-spline function * * F = USPF(X,A,N,XA,XB) function value at X * * X = abscissa value * A(1)...A(N) = B-spline coefficients * N = number of coefficients * XA,XB = range of B-spline function * * Extrapolation is done for X outside the range XA,XB. * REAL E(N) * ... S=FLOAT(N-1)/(XB-XA) R=S*(X-XA) I=MAX(1,MIN(INT(R),N-3)) Z=R-I A= ( E(I)+4.0*E(I+1)+ E(I+2) )/6.0 B=0.5*(-E(I) + E(I+2) ) C= ( E(I)-2.0*E(I+1)+ E(I+2) ) D= (-E(I)+3.0*E(I+1)-3.0*E(I+2)+E(I+3)) USPF=A+Z*(B+0.5*Z*(C+Z*D/3.0)) 100 RETURN END FUNCTION USPG(X,E,N,XA,XB,NDER) * * Compute function value or derivative for B-spline function * * D = USPG(...,NDER) NDER'th derivative at value X * * X = abscissa value * A(1)...A(N) = B-spline coefficients * N = number of coefficients * XA,XB = range of B-spline function * * Extrapolation is done for X outside the range XA,XB. * REAL E(N) S=FLOAT(N-1)/(XB-XA) R=S*(X-XA) I=MAX(1,MIN(INT(R),N-3)) Z=R-I A= ( E(I)+4.0*E(I+1)+ E(I+2) )/6.0 B=0.5*(-E(I) + E(I+2) ) C= ( E(I)-2.0*E(I+1)+ E(I+2) ) D= (-E(I)+3.0*E(I+1)-3.0*E(I+2)+E(I+3)) IF(NDER.EQ.0) THEN USPG=A+Z*(B+0.5*Z*(C+Z*D/3.0)) ELSE IF(NDER.EQ.1) THEN USPG=(B+Z*(C+0.5*Z*D))*S ELSE IF(NDER.EQ.2) THEN USPG=(C+Z*D)*S*S ELSE IF(NDER.EQ.3) THEN USPG=D*S*S*S ELSE USPG=0.0 END IF 100 RETURN END FUNCTION USPZ(A,NSP,XA,XB) * * calculate curvature (square of second der) for whole curve * REAL A(NSP) HX=FLOAT(NSP-1)/(XB-XA) SUM=0.0 DL=2.0*A(1)-5.0*A(2)+4.0*A(3)-A(4) DO 10 I=2,NSP-1 DN=A(I-1)-2.0*A(I)+A(I+1) SUM=SUM+DN*DN+DN*DL+DL*DL XC=FLOAT(I)/HX CC=DN*DN+DN*DL+DL*DL 10 DL=DN DN=-A(NSP-3)+4.0*A(NSP-2)-5.0*A(NSP-1)+2.0*A(NSP) USPZ=(SUM+DN*DN+DN*DL+DL*DL)*HX*HX*HX/(XB-XA) 100 RETURN END SUBROUTINE USPU(X,NSP,XA,XB,INULL,S,T) * * Calculate function values and second derivatives for the * individual B-splines * * S(INULL+I) = function values * T(INULL+I) = second derivatives * REAL S(4),T(4) * ... XX=MIN(MAX(X,XA),XB) HX=FLOAT(NSP-1)/(XB-XA) DX=1.0+(XX-XA)*HX I0=DX IF(I0.GE.NSP) I0=NSP-1 V=DX-FLOAT(I0) U=1.0-V INULL=I0-2 * * function value * S(1)= (U*U*U) /6.0 S(2)= (1.0+3.0*(1.0+U*V)*U) /6.0 S(3)= (1.0+3.0*(1.0+U*V)*V) /6.0 S(4)= (V*V*V) /6.0 * * second derivative * T(1)= (U) *HX*HX T(2)= (V-2.0*U) *HX*HX T(3)= (U-2.0*V) *HX*HX T(4)= (V) *HX*HX * * apply not-a-knot condition at endpoints * IF(INULL.LT.0) THEN SS=S(1) S(1)=S(2)+4.0*SS S(2)=S(3)-6.0*SS S(3)=S(4)+4.0*SS S(4)=-SS TT=T(1) T(1)=T(2)+4.0*TT T(2)=T(3)-6.0*TT T(3)=T(4)+4.0*TT T(4)=-TT INULL=0 ELSE IF(INULL.EQ.NSP-3) THEN SS=S(4) S(4)=S(3)+4.0*SS S(3)=S(2)-6.0*SS S(2)=S(1)+4.0*SS S(1)=-SS TT=T(4) T(4)=T(3)+4.0*TT T(3)=T(2)-6.0*TT T(2)=T(1)+4.0*TT T(1)=-TT INULL=NSP-4 END IF 100 RETURN END SUBROUTINE ESPB(X,N,XA,XB,INULL,S) * * returns S values (inside XA ... XB) * REAL S(4) * ... HX=(XB-XA)/FLOAT(N-3) XI=(X-XA)/HX IF(XI.LT.0.0) THEN * outside low INULL=0 V=0.0 ELSE INULL=XI IF(INULL.GT.N-4) THEN * outside high INULL=N-4 V=1.0 ELSE * inside V=XI-FLOAT(INULL) END IF END IF U=1.0-V * function values of the four splines S(1)= U*U*U/6.0 S(2)= (1.0+3.0*(1.0+U*V)*U)/6.0 S(3)= (1.0+3.0*(1.0+U*V)*V)/6.0 S(4)= V*V*V/6.0 * 100 RETURN END SUBROUTINE ESPC(Y,A,N) * * Calculation of coefficients from Y(.) values * * - - * CALL ESPC(Y,A,N) * - * Given n function values Y(1)...Y(n-2) at equidistant X-values, the * B-spline coeffcients A(1)...A(n) are calculated. Arrays Y and A * may be identical. The so called not-a-knot condition is used on * both sides (no discontinuity in third derivative at first and last * inner knot). * * For N > 100 the user has to include a common/spcom/ with at least * 2*n+6 words. * REAL Y(N),A(N) * ... DO I=N-1,2,-1 A(I)=Y(I-1) END DO * determination of spline coefficients with knot-a-knot condition CALL USPC(A(2),A(2),N-2) * ... with determination of first/last coefficient A(1)=4.0*A(2 )-6.0*A(3 )+4.0*A(4 )-A(5 ) A(N)=4.0*A(N-1)-6.0*A(N-2)+4.0*A(N-3)-A(N-4) RETURN END FUNCTION ESPF(X,A,N,XA,XB) * * calculate function value at x * * FX= ESPF(X,A,N,N,XA,XB) A(1)...A(N) = B-spline coefficients * XA, XB = Lower and upper limit * * = function value at X for XA < X < XB * For X outside the limits the x value is assumed at the limit * * The result of the call is a linear combination of B-spline * coefficients. After the call the weights are stored in * COMMON/SPCOL/I0,N0,S(200) * The returned value is: * SUM(I=1...4) S(I)*A(I0+I) * REAL A(N) COMMON/SPCOM/I0,N0,S(200) * ... ND=0 HX=(XB-XA)/FLOAT(N-3) XI=(X-XA)/HX IF(XI.LT.0.0) THEN * outside low I0=0 V=0.0 ELSE I0=XI IF(I0.GT.N-4) THEN * outside high I0=N-4 V=1.0 ELSE * inside V=XI-FLOAT(I0) END IF END IF U=1.0-V * function value FAC=1.0D0/6.0D0 S(1)= FAC* U*U*U S(2)= FAC*(1.0+3.0*(1.0+U*V)*U) S(3)= FAC*(1.0+3.0*(1.0+U*V)*V) S(4)= FAC* V*V*V * final sum ESPF=A(I0+1)*S(1)+A(I0+2)*S(2)+A(I0+3)*S(3)+A(I0+4)*S(4) 100 RETURN END FUNCTION ESPD(X,A,N,XA,XB,NDER) * * FD= ESPD(X,A,N,N,XA,XB,ND) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * FD= function value at X for ND=0 * = ND'th derivative for ND=1, 2 OR 3 * = integral between XA and X for ND=-1 *) * * For X outside the limits the X value is assumed at the limit. * * The result of the call is for all cases a linear combination of * B-spline coefficients. After the call the weights are stored in * the COMMON/SPCOM/I0,N0,S(200) * For ND=0...3 the returned value is: * sum(i=1...4) A(i)*A(i0+i) * For ND= -1 the returned value is: * sum(i=1...n0) S(i)*A(i) * * *) If, for ND=-1, N is larger than 200, the user has to include a * COMMON/SPCOM/I0,N0,S(N). * * REAL A(N) COMMON/SPCOM/I0,N0,S(200) * ... IF(NDER.GT.3.0.OR.NDER.LT.0) THEN I0=-1 ESPD=0.0 IF(NDER+1.EQ.0) ESPD=ESPI(XA,X,A,N,XA,XB) GOTO 100 END IF ND=NDER HX=(XB-XA)/FLOAT(N-3) XI=(X-XA)/HX IF(XI.LT.0.0) THEN * outside low I0=0 V=0.0 ELSE I0=XI IF(I0.GT.N-4) THEN * outside high I0=N-4 V=1.0 ELSE * inside V=XI-FLOAT(I0) END IF END IF U=1.0-V IF(ND.EQ.0) THEN * function value FAC=1.0D0/6.0D0 S(1)= FAC* U*U*U S(2)= FAC*(1.0+3.0*(1.0+U*V)*U) S(3)= FAC*(1.0+3.0*(1.0+U*V)*V) S(4)= FAC* V*V*V ELSE IF(ND.EQ.1) THEN * first derivative FAC=0.5D0/HX S(1)=-FAC* U*U S(2)=-FAC*(1.0+(2.0*V-U)*U) S(3)= FAC*(1.0+(2.0*U-V)*V) S(4)= FAC* V*V ELSE IF(ND.EQ.2) THEN * second derivative FAC=1.0D0/(HX*HX) S(1)= FAC* U S(2)= FAC*(V-2.0*U) S(3)= FAC*(U-2.0*V) S(4)= FAC* V ELSE IF(ND.EQ.3) THEN * third derivative FAC=1.0D0/(HX*HX*HX) S(1)=-FAC S(2)= FAC *3.0 S(3)=-S(2) S(4)=-S(1) END IF * final sum ESPD=A(I0+1)*S(1)+A(I0+2)*S(2)+A(I0+3)*S(3)+A(I0+4)*S(4) 100 RETURN END FUNCTION ESPI(X1,X2,A,N,XA,XB) * * Integral of spline function between limits X1 and X2 * * FX= ESPI(X1,X2,A,N,XA,XB) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * The integral over the spline function is calculated between the * limits X1 and X2. A limit below XA (above XB) is replaced by XA * (XB). The result is a linear combination of B-spline coefficients. * Integral = sum(i=1...n0) S(i)*A(i) * After the call the weights are stored in the * COMMON/SPCOM/I0,N0,S(200). * If N is larger than 200, the user has to include a common * COMMON/SPCOM/I0,N0,S(N). * REAL A(N) COMMON/SPCOM/I0,N0,S(200) REAL T(4),BT(4)/1.0,12.0,23.0,24.0/ REAL TB(3,3)/11.0,11.0, 1.0,22.0,12.0, 1.0,23.0,12.0,1.0/ DOUBLE PRECISION SUM * ... * clear array of weights DO I=1,N S(I)=0.0 END DO N0=0 HX=(XB-XA)/FLOAT(N-3) * lower limit of integration XJ=X1 DO J=1,2 * transform to coordinate v relative to knots XJ=(XJ-XA)/HX IF(XJ.GE.0.0) THEN I0=XJ IF(I0.LE.N-4) THEN V=XJ-FLOAT(I0) ELSE I0=N-4 V =1.0 END IF U=1.0-V * calculate contribution from last knot interval T(1)=1.0-U*U*U*U T(2)=11.0-(((-3.0*U+4.0)*U+6.0)*U+4.0)*U T(3)=(((-3.0*V+4.0)*V+6.0)*V+4.0)*V T(4)=V*V*V*V * add contributions from previous knot intervals IF(I0.GE.1) THEN T(1)=T(1)+TB(1,MIN0(I0,3)) T(2)=T(2)+TB(2,MIN0(I0,3)) T(3)=T(3)+TB(3,MIN0(I0,3)) END IF * subtract or add to array s(i) DO I=1,I0+4 IF(I.LE.I0) THEN B=BT(MIN0(I,4)) ELSE B=T(I-I0) END IF IF(J.EQ.1) THEN S(I)=-B ELSE S(I)=S(I)+B END IF N0=MAX0(N0,I0+4) END DO END IF * upper limit of integration XJ=X2 END DO * integral by sum SUM=0.0 DO I=1,N IF(S(I).GT.1.0E-10) THEN SUM=SUM+A(I)*S(I) S(I)=S(I)*HX/24.0 ELSE S(I)=0.0 END IF END DO ESPI=SUM*HX/24.0 100 RETURN END SUBROUTINE ESPR(A,N,XA,XB) * * Print table of values of spline function, derivative, integral. * * CALL ESPR(A,N,XA,XB) A(1)...A(N) = B-spline coefficients * XA, XB = lower and upper limit * * If N is larger than 200, the user has to include a common * COMMON/SPCOM/I0,N0,S(N). * REAL A(N),F(0:4) * ... WRITE(*,101) M=1 WRITE(*,102) M,A(1) DO I=1,2*N-5 X=(FLOAT(2*N-5-I)*XA+FLOAT(I-1)*XB)/FLOAT(2*N-6) F( 0)=ESPF(X,A,N,XA,XB) DO ND=1,3 F(ND)=ESPD(X,A,N,XA,XB,ND) END DO F( 4)=ESPI(XA,X,A,N,XA,XB) IF(MOD(I,2).EQ.1) THEN J=2+I/2 WRITE(*,102) J,A(J),X,F ELSE WRITE(*,103) X,F END IF END DO WRITE(*,102) N,A(N) 100 RETURN 101 FORMAT('0',3X,'I',8X,'Coeff',10X,'X',10X,'F',10X, + 2HF',10X,3HF'',10X,4HF''',10X,'INT F'/) 102 FORMAT(1X,I4,2X,3G15.5,5X,3G13.5,5X,G15.5) 103 FORMAT(22X,2G15.5,5X,3G13.5,5X,G15.5) END FUNCTION ESPX(XL,FZ,A,N,XA,XB,NDER) * * Calculate next x-value with given function value or derivative * * -- -- - - -- -- ---- * X = ESPX(XL,FZ,A,N,XA,XB,NDER) NDER = -1, 0, 1, 2 * * = next abscissa value with value = FD with X GE XA * X > XB returned, if no (further) value can be found. * * Subprogram is not optimized for speed. * REAL A(N) * ... * knot distance HX HX=(XB-XA)/FLOAT(N-3) EPS=0.1*HX ESPX=XB+HX IF(XL.GT.XB) GOTO 100 IF(XL.LE.XA) THEN XX=XA ELSE XX=XL+EPS END IF * * loop to find values with opposite sign * X1=XX F1=ESPD(X1,A,N,XA,XB,NDER)-FZ 10 X2=AMIN1(XB,X1+3.0*EPS) F2=ESPD(X2,A,N,XA,XB,NDER)-FZ IF(F1*F2.GT.0.0) THEN * still same signs of two function values IF(X2.EQ.XB) GOTO 100 X1=X2 F1=F2 GOTO 10 END IF * * enclosed * IF(F1.EQ.F2) GOTO 100 * approximate solution ... XX=(X1*F2-X2*F1)/(F2-F1) * * ... refined by Newton method * IT=0 30 IF(IT.LE.10) THEN * function value ... FN=ESPD(XX,A,N,XA,XB,NDER)-FZ * ... and derivative FD=ESPD(XX,A,N,XA,XB,NDER+1) IF(FD.NE.0.0) THEN * correction DX=-FN/FD XX=XX+DX IT=IT+1 IF(IT.LT.2.OR.ABS(DX).GT.1.0E-3*EPS) GOTO 30 ESPX=XX END IF END IF 100 RETURN END SUBROUTINE ESPQ(Y,Y2,NH,XA,XB,NKN,F,A) * * least square fit of spline to histogram * Y = histogram content, Y2 = histogram content variance * NH=4*(NKN-3) is recommended and corresponds to 4 bins per knot * REAL Y(*),Y2(*),F(*),A(*) COMMON/SPCOM/I0,N0,S(200) REAL Q(5,63) * ... * reset matrices DO J=1,NKN A(J)=0.0 DO I=1,5 Q(I,J)=0.0 END DO END DO * find minimum value of variance (not zero) VAR=0.0 DO K=1,NH IF(Y2(K).NE.0.0) THEN IF(VAR.EQ.0.0) VAR=Y2(K) VAR=MIN(VAR,Y2(K)) END IF END DO IF(VAR.EQ.0.0) VAR=1.0 * loop all bins XR=XA DO K=1,NH * left and right bin limits XL=XR XR=XA+(XB-XA)*FLOAT(K)/FLOAT(NH) C IF(Y(K).NE.0.0) THEN *C nonzero bin content WT=1.0/VAR IF(Y2(K).NE.0.0) WT=1.0/Y2(K) * integration over bin XL ... XR SUM= ESPI(XL,XR,A,NKN,XA,XB) * add to normal equations DO I=1,N0-I0 A(I0+I)=A(I0+I)+WT*S(I0+I)*Y(K) DO J=I,N0-I0 Q(J-I+1,I0+I)=Q(J-I+1,I0+I)+WT*S(I0+I)*S(I0+J) END DO END DO C END IF END DO * solve system with band matrix CALL BMDEQU(Q,5,NKN,A,DELTA) * loop to calculate fitted bin content XR=XA DO K=1,NH * left and right bin limits XL=XR XR=XA+(XB-XA)*FLOAT(K)/FLOAT(NH) * integration over bin XL...XR F(K)=ESPI(XL,XR,A,NKN,XA,XB) END DO * renormalize spline coefficients DO J=1,NKN A(J)=A(J)*(XB-XA)/FLOAT(NH) END DO END