C @(#)nagi.for 17.1.1.1 (ES0-DMD) 01/25/02 17:40:35 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C SUBROUTINE E02ADZ(MFIRST,MLAST,MTOT,KPLUS1,NROWS,KALL,NDV,X,Y,W, * XMIN,XMAX,INUP1,NU,WORK1,WORK2,A,S,SERR,EPS, * IFAIL) C MARK 7 RELEASE. NAG COPYRIGHT 1978. C MARK 8 REVISED. IER-228 (APR 1980). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C E02ADZ COMPUTES WEIGHTED LEAST-SQUARES POLYNOMIAL C APPROXIMATIONS TO AN ARBITRARY SET OF DATA POINTS, C WITH, IF REQUIRED, SEVERAL SETS OF VALUES OF THE C DEPENDENT VARIABLE. C C FORSYTHE-CLENSHAW METHOD WITH MODIFICATIONS DUE TO C REINSCH AND GENTLEMAN. C C STARTED - 1973. C COMPLETED - 1978. C AUTHORS - MGC AND GTA. C C WORK1 AND WORK2 ARE WORKSPACE AREAS. C WORK1(1, R) CONTAINS THE VALUE OF X(R) TRANSFORMED C TO THE RANGE -1 TO +1. C WORK1(2, R) CONTAINS THE WEIGHTED VALUE OF THE CURRENT C ORTHOGONAL POLYNOMIAL (OF DEGREE I) AT THE R TH C DATA POINT. C WORK2(1, J) CONTAINS THE COEFFICIENT OF THE CHEBYSHEV C POLYNOMIAL OF DEGREE J - 1 IN THE CHEBYSHEV-SERIES C REPRESENTATION OF THE CURRENT ORTHOGONAL POLYNOMIAL C (OF DEGREE I). C WORK2(2, J) CONTAINS THE COEFFICIENT OF THE CHEBYSHEV C POLYNOMIAL OF DEGREE J - 1 IN THE CHEBYSHEV-SERIES C REPRESENTATION OF THE PREVIOUS ORTHOGONAL POLYNOMIAL C (OF DEGREE I - 1). C C .. Scalar Arguments .. DOUBLE PRECISION XMAX, XMIN INTEGER IFAIL, INUP1, KALL, KPLUS1, MFIRST, MLAST, MTOT, * NDV, NROWS C .. Array Arguments .. DOUBLE PRECISION A(NDV,NROWS,KPLUS1), EPS(NDV,MLAST), NU(INUP1), * S(NDV,KPLUS1), SERR(KPLUS1), W(MLAST), * WORK1(2,MTOT), WORK2(2,KPLUS1), X(MLAST), * Y(NDV,MLAST) C .. Local Scalars .. DOUBLE PRECISION ALPIP1, BETAI, BJ, BJP1, BJP2, CIL, D, DF, DI, * DIM1, DJ, EPSLR, FACTOR, PIJ, SIGMAI, WR, WRPR, * WRPRSQ, X1, XCAPR, XM INTEGER I, IERROR, II, IM1, INU, IPLUS1, IPLUS2, J, * JPLUS1, JPLUS2, JREV, K, L, M, MDIST, MR, R LOGICAL WNZ C .. Local Arrays .. DOUBLE PRECISION CI(10) C .. Intrinsic Functions .. INTRINSIC SQRT C .. Executable Statements .. K = KPLUS1 - 1 INU = INUP1 - 1 C C TEST THE VALIDITY OF THE DATA. C C CHECK INPUT PARAMETERS. C M = MLAST - MFIRST + 1 I = KPLUS1 - INU IERROR = 5 IF (MFIRST.LT.1 .OR. INUP1.LT.1 .OR. KPLUS1.LT.INUP1 .OR. M.LT. * I .OR. NDV.LT.1 .OR. (KALL.NE.1 .AND. KALL.NE.0)) GO TO 600 C C CHECK THAT THE VALUES OF X(R) ARE NON-DECREASING AND C DETERMINE THE NUMBER (MDIST) OF DISTINCT VALUES OF X(R) C WITH NON-ZERO WEIGHT C IERROR = 2 MDIST = 1 IF (W(MFIRST).EQ.0.0D+0) MDIST = 0 L = MFIRST + 1 IF (L.GT.MLAST) GO TO 40 WNZ = W(MFIRST) .NE. 0.0D+0 DO 20 R = L, MLAST IF (X(R).LT.X(R-1)) GO TO 600 IF (X(R).GT.X(R-1)) WNZ = .FALSE. IF (W(R).EQ.0.0D+0 .OR. WNZ) GO TO 20 MDIST = MDIST + 1 WNZ = .TRUE. 20 CONTINUE C C CHECK THAT XMIN.LT.XMAX AND THAT XMIN AND XMAX SPAN THE DATA C X VALUES. C 40 IERROR = 1 IF (XMIN.GT.X(MFIRST) .OR. XMAX.LT.X(MLAST) .OR. XMIN.GE.XMAX) * GO TO 600 C C IF THE NUMBER OF DISTINCT VALUES OF X(R) WITH NON-ZERO C WEIGHT IS LESS THAN THE NUMBER OF INDEPENDENT COEFFICIENTS C IN THE FIT OF MAXIMUM DEGREE K THERE IS NO UNIQUE C POLYNOMIAL C APPROXIMATION OF THAT DEGREE. C L = K - INU IERROR = 3 IF (MDIST.LE.L) GO TO 600 C C CHECK THAT NROWS HAS BEEN SET SUFFICIENTLY LARGE. C IERROR = 5 IF (KALL.EQ.1 .AND. NROWS.LT.KPLUS1) GO TO 600 IF (INUP1.EQ.1) GO TO 80 C C NORMALIZE THE FORCING FACTOR SO THAT ITS LEADING COEFFICIENT C IS UNITY, CHECKING THAT THIS COEFFICIENT WAS NOT ZERO. C IERROR = 4 DI = NU(INUP1) IF (DI.EQ.0.0D0) GO TO 600 DO 60 I = 1, INUP1 WORK2(1,I) = NU(I)/DI WORK2(2,I) = 0.0D0 60 CONTINUE C 80 IERROR = 0 X1 = XMIN XM = XMAX D = XM - X1 C C THE INITIAL VALUES OF EPS(L,R) (L = 1,2,....NDV AND R = C MFIRST, MFIRST+1,....MLAST) OF THE WEIGHTED RESIDUALS AND C THE VALUES WORK1(1,R)(R=1,2...M) OF THE NORMALIZED C INDEPENDENT VARIABLE ARE COMPUTED. N.B. WORK1(1,R) IS C COMPUTED FROM THE EXPRESSION BELOW RATHER THAN THE MORE C NATURAL FORM (2.0*X(R) - X1 - XM)/D C SINCE THE FORMER GUARANTEES THE COMPUTED VALUE TO DIFFER FROM C THE TRUE VALUE BY AT MOST 4.0*MACHINE ACCURACY, WHEREAS THE C LATTER HAS NO SUCH GUARANTEE. C C MDIST IS NOW USED TO RECORD THE TOTAL NUMBER OF DATA POINTS C WITH NON-ZERO WEIGHT. C MDIST = 0 DO 120 R = MFIRST, MLAST WR = W(R) IF (WR.NE.0.0D0) MDIST = MDIST + 1 MR = R - MFIRST + 1 DO 100 L = 1, NDV EPS(L,R) = WR*Y(L,R) 100 CONTINUE WORK1(1,MR) = ((X(R)-X1)-(XM-X(R)))/D 120 CONTINUE IM1 = INU*KALL + 1 BETAI = 0.0D0 DO 160 JPLUS1 = 1, KPLUS1 SERR(JPLUS1) = 0.0D0 DO 140 L = 1, NDV A(L,IM1,JPLUS1) = 0.0D0 140 CONTINUE 160 CONTINUE DO 560 IPLUS1 = INUP1, KPLUS1 C C SET STARTING VALUES FOR DEGREE I. C II = (IPLUS1-1)*KALL + 1 IPLUS2 = IPLUS1 + 1 IF (IPLUS1.EQ.KPLUS1) GO TO 240 IF (KALL.EQ.0) GO TO 220 DO 200 JPLUS1 = IPLUS2, KPLUS1 DO 180 L = 1, NDV A(L,II,JPLUS1) = 0.0D0 180 CONTINUE 200 CONTINUE 220 WORK2(1,IPLUS2) = 0.0D0 WORK2(2,IPLUS2) = 0.0D0 240 ALPIP1 = 0.0D0 DI = 0.0D0 DO 260 L = 1, NDV CI(L) = 0.0D0 260 CONTINUE WORK2(1,IPLUS1) = 1.0D0 IF (KPLUS1.GT.1) WORK2(2,1) = WORK2(1,2) DO 440 R = MFIRST, MLAST IF (W(R).EQ.0.0D0) GO TO 440 MR = R - MFIRST + 1 XCAPR = WORK1(1,MR) C C THE WEIGHTED VALUE WORK1(2, R) OF THE ORTHOGONAL POLYNOMIAL C OF DEGREE I AT X = X(R) IS COMPUTED BY RECURRENCE FROM ITS C CHEBYSHEV-SERIES REPRESENTATION. C IF (IPLUS1.GT.1) GO TO 280 WRPR = W(R)*0.5D0*WORK2(1,1) WORK1(2,MR) = WRPR GO TO 400 280 J = IPLUS2 IF (XCAPR.GT.0.5D0) GO TO 360 IF (XCAPR.GE.-0.5D0) GO TO 320 C C GENTLEMAN*S MODIFIED RECURRENCE. C FACTOR = 2.0D0*(1.0D0+XCAPR) DJ = 0.0D0 BJ = 0.0D0 DO 300 JREV = 2, IPLUS1 J = J - 1 DJ = WORK2(1,J) - DJ + FACTOR*BJ BJ = DJ - BJ 300 CONTINUE WRPR = W(R)*(0.5D0*WORK2(1,1)-DJ+0.5D0*FACTOR*BJ) WORK1(2,MR) = WRPR GO TO 400 C C CLENSHAW*S ORIGINAL RECURRENCE. C 320 FACTOR = 2.0D0*XCAPR BJP1 = 0.0D0 BJ = 0.0D0 DO 340 JREV = 2, IPLUS1 J = J - 1 BJP2 = BJP1 BJP1 = BJ BJ = WORK2(1,J) - BJP2 + FACTOR*BJP1 340 CONTINUE WRPR = W(R)*(0.5D0*WORK2(1,1)-BJP1+0.5D0*FACTOR*BJ) WORK1(2,MR) = WRPR GO TO 400 C C REINSCH*S MODIFIED RECURRENCE. C 360 FACTOR = 2.0D0*(1.0D0-XCAPR) DJ = 0.0D0 BJ = 0.0D0 DO 380 JREV = 2, IPLUS1 J = J - 1 DJ = WORK2(1,J) + DJ - FACTOR*BJ BJ = BJ + DJ 380 CONTINUE WRPR = W(R)*(0.5D0*WORK2(1,1)+DJ-0.5D0*FACTOR*BJ) WORK1(2,MR) = WRPR C C THE COEFFICIENTS CI(L) OF THE I TH ORTHOGONAL POLYNOMIAL C L=1,2....NDV AND THE COEFFICIENTS ALPIP1 AND BETAI IN THE C THREE-TERM RECURRENCE RELATION FOR THE ORTHOGONAL C POLYNOMIALS ARE COMPUTED. C 400 WRPRSQ = WRPR**2 DI = DI + WRPRSQ DO 420 L = 1, NDV CI(L) = CI(L) + WRPR*EPS(L,R) 420 CONTINUE ALPIP1 = ALPIP1 + WRPRSQ*XCAPR 440 CONTINUE DO 460 L = 1, NDV CI(L) = CI(L)/DI 460 CONTINUE IF (IPLUS1.NE.INUP1) BETAI = DI/DIM1 ALPIP1 = 2.0D0*ALPIP1/DI C C THE WEIGHTED RESIDUALS EPS(L,R)(L=1,2....NDV AND R=MFIRST, C MFIRST+1....MLAST) FOR DEGREE I ARE COMPUTED, TOGETHER C WITH THEIR SUM OF SQUARES, SIGMAI C DF = MDIST - (IPLUS1-INU) DO 500 L = 1, NDV CIL = CI(L) SIGMAI = 0.0D0 DO 480 R = MFIRST, MLAST IF (W(R).EQ.0.0D0) GO TO 480 MR = R - MFIRST + 1 EPSLR = EPS(L,R) - CIL*WORK1(2,MR) EPS(L,R) = EPSLR SIGMAI = SIGMAI + EPSLR**2 480 CONTINUE C C THE ROOT MEAN SQUARE RESIDUAL S(L, I + 1) FOR DEGREE I C IS THEORETICALLY UNDEFINED IF M = I + 1 - INU (THE C CONDITION FOR THE POLYNOMIAL TO PASS EXACTLY THROUGH THE C DATA POINTS). SHOULD THIS CASE ARISE THE R.M.S. RESIDUAL C IS SET TO ZERO. C IF (DF.LE.0.0D0) S(L,IPLUS1) = 0.0D0 IF (DF.GT.0.0D0) S(L,IPLUS1) = SQRT(SIGMAI/DF) 500 CONTINUE C C THE CHEBYSHEV COEFFICIENTS A(L, I+1, 1), A(L, I+1, 2).... C A(L, I+1, I+1) IN THE POLYNOMIAL APPROXIMATION OF DEGREE I C TO EACH SET OF VALUES OF THE INDEPENDENT VARIABLE C (L=1,2,...,NDV) TOGETHER WITH THE COEFFICIENTS C WORK2(1, 1), WORK2(1, 2), ..., WORK2(1, I + 1), IN THE C CHEBYSHEV-SERIES REPRESENTATION OF THE (I + 1) TH C ORTHOGONAL POLYNOMIAL ARE COMPUTED. C DO 540 JPLUS1 = 1, IPLUS1 JPLUS2 = JPLUS1 + 1 PIJ = WORK2(1,JPLUS1) SERR(JPLUS1) = SERR(JPLUS1) + PIJ**2/DI DO 520 L = 1, NDV A(L,II,JPLUS1) = A(L,IM1,JPLUS1) + CI(L)*PIJ 520 CONTINUE IF (JPLUS1.EQ.KPLUS1) GO TO 560 WORK2(1,JPLUS1) = WORK2(1,JPLUS2) + WORK2(2,JPLUS1) - * ALPIP1*PIJ - BETAI*WORK2(2,JPLUS2) WORK2(2,JPLUS2) = PIJ 540 CONTINUE DIM1 = DI IM1 = II 560 CONTINUE DO 580 IPLUS1 = 1, KPLUS1 SERR(IPLUS1) = 1.0D0/SQRT(SERR(IPLUS1)) 580 CONTINUE 600 IFAIL = IERROR RETURN END SUBROUTINE E02AEF(NPLUS1,A,XCAP,P,IFAIL) C NAG LIBRARY SUBROUTINE E02AEF C C E02AEF EVALUATES A POLYNOMIAL FROM ITS CHEBYSHEV- C SERIES REPRESENTATION. C C CLENSHAW METHOD WITH MODIFICATIONS DUE TO REINSCH C AND GENTLEMAN. C C USES NAG LIBRARY ROUTINES P01AAF AND X02AAF. C USES INTRINSIC FUNCTION ABS. C C STARTED - 1973. C COMPLETED - 1976. C AUTHOR - MGC AND JGH. C C NAG COPYRIGHT 1975 C MARK 5 RELEASE C MARK 7 REVISED IER-140 (DEC 1978) C MARK 9 REVISED. IER-352 (SEP 1981) C MARK 11.5(F77) REVISED. (SEPT 1985.) C C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='E02AEF') C .. Scalar Arguments .. DOUBLE PRECISION P, XCAP INTEGER IFAIL, NPLUS1 C .. Array Arguments .. DOUBLE PRECISION A(NPLUS1) C .. Local Scalars .. DOUBLE PRECISION BK, BKP1, BKP2, DK, ETA, FACTOR INTEGER IERROR, K, KREV, N, NPLUS2 C .. Local Arrays .. CHARACTER*1 P01REC(1) C .. External Functions .. DOUBLE PRECISION X02AAF INTEGER P01ABF EXTERNAL X02AAF, P01ABF C .. Intrinsic Functions .. INTRINSIC ABS C .. Executable Statements .. IERROR = 0 ETA = X02AAF(ETA) C INSERT CALL TO X02AAF C C ETA IS THE SMALLEST POSITIVE NUMBER SUCH THAT C THE COMPUTED VALUE OF 1.0 + ETA EXCEEDS UNITY. C IF (NPLUS1.GE.1) GO TO 20 IERROR = 2 GO TO 180 20 IF (ABS(XCAP).LE.1.0D0+4.0D0*ETA) GO TO 40 IERROR = 1 P = 0.0D0 GO TO 180 40 IF (NPLUS1.GT.1) GO TO 60 P = 0.5D0*A(1) GO TO 180 60 N = NPLUS1 - 1 NPLUS2 = N + 2 K = NPLUS2 IF (XCAP.GT.0.5D0) GO TO 140 IF (XCAP.GE.-0.5D0) GO TO 100 C C GENTLEMAN*S MODIFIED RECURRENCE. C FACTOR = 2.0D0*(1.0D0+XCAP) DK = 0.0D0 BK = 0.0D0 DO 80 KREV = 1, N K = K - 1 DK = A(K) - DK + FACTOR*BK BK = DK - BK 80 CONTINUE P = 0.5D0*A(1) - DK + 0.5D0*FACTOR*BK GO TO 180 C C CLENSHAW*S ORIGINAL RECURRENCE. C 100 FACTOR = 2.0D0*XCAP BKP1 = 0.0D0 BK = 0.0D0 DO 120 KREV = 1, N K = K - 1 BKP2 = BKP1 BKP1 = BK BK = A(K) - BKP2 + FACTOR*BKP1 120 CONTINUE P = 0.5D0*A(1) - BKP1 + 0.5D0*FACTOR*BK GO TO 180 C C REINSCH*S MODIFIED RECURRENCE. C 140 FACTOR = 2.0D0*(1.0D0-XCAP) DK = 0.0D0 BK = 0.0D0 DO 160 KREV = 1, N K = K - 1 DK = A(K) + DK - FACTOR*BK BK = BK + DK 160 CONTINUE P = 0.5D0*A(1) + DK - 0.5D0*FACTOR*BK 180 IF (IERROR) 200, 220, 200 200 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC) RETURN 220 IFAIL = 0 RETURN END SUBROUTINE E02CAF(M,N,K,L,X,Y,F,W,NX,A,NA,XMIN,XMAX,NUX,INUXP1, * NUY,INUYP1,WORK,NWORK,IFAIL) C MARK 7 RELEASE. NAG COPYRIGHT 1978. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C E02CAF IS A CALLING ROUTINE WHICH SETS UP WORK SPACE ARRAYS C FOR THE SUBROUTINE E02CAZ AND THEN CALLS IT TO OBTAIN AN C APPROXIMATION TO THE LEAST SQUARES POLYNOMIAL SURFACE FIT C TO DATA ARBITRARILY DISTRIBUTED ON LINES PARALLEL TO ONE C INDEPENDENT COORDINATE AXIS. C C STARTED - 1978. C COMPLETED - 1978. C AUTHOR - GTA. C C C FIND MAXIMUM OF N AND THE ELEMENTS OF M, AND LARGER OF K AND C L, AND THE SUM OF THE ELEMENTS OF M C C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='E02CAF') C .. Scalar Arguments .. INTEGER IFAIL, INUXP1, INUYP1, K, L, N, NA, NWORK, NX C .. Array Arguments .. DOUBLE PRECISION A(NA), F(NX), NUX(INUXP1), NUY(INUYP1), W(NX), * WORK(NWORK), X(NX), XMAX(N), XMIN(N), Y(N) INTEGER M(N) C .. Local Scalars .. INTEGER I, IERROR, IW, KP1, MI, MJ, MJP1, MSUM, NCI, * NCII, NFI, NRESID, NS, NW, NWI, NWT, NYW C .. Local Arrays .. CHARACTER*1 P01REC(1) C .. External Functions .. INTEGER P01ABF EXTERNAL P01ABF C .. External Subroutines .. EXTERNAL E02CAZ C .. Executable Statements .. IERROR = 1 IF (INUXP1.LT.1 .OR. INUYP1.LT.1 .OR. INUXP1.GT.K+1 .OR. * INUYP1.GT.L+1) GO TO 80 IF (N.LT.L-INUYP1+2) GO TO 80 MSUM = 0 MI = N DO 20 I = 1, N MJ = M(I) IF (MJ.LT.K-INUXP1+2) GO TO 80 MSUM = MSUM + MJ IF (MJ.GT.MI) MI = MJ 20 CONTINUE MJ = L IF (K.GT.L) MJ = K MJP1 = MJ + 1 C C CHECK THAT THE INTEGER PARAMETERS HAVE REASONABLE VALUES C KP1 = K + 1 MJ = KP1*(L+1) NW = 2*MI + 2*N*KP1 + 5*MJP1 + 2*N + MSUM IF (NX.LT.MSUM .OR. NA.LT.MJ .OR. NWORK.LT.NW) GO TO 80 C C CHECK THAT THE Y(I) ARE STRICTLY INCREASING C IF (N.EQ.1) GO TO 60 IERROR = 3 DO 40 I = 2, N IF (Y(I).LE.Y(I-1)) GO TO 80 40 CONTINUE IERROR = 1 C C ALLOCATE WORKSPACE FOR E02CAZ C 60 NWI = 1 NCI = NWI + N NYW = NCI + N*KP1 NCII = NYW + N*KP1 NWT = NCII + MJP1 NS = NWT + MJP1 NFI = NS + MJP1 NRESID = NFI + N NW = NRESID + MSUM IW = NWORK - NW + 1 CALL E02CAZ(M,N,KP1,L,X,Y,F,W,NX,A,NA,XMIN,XMAX,NUX,INUXP1,NUY, * INUYP1,WORK(NWI),WORK(NCI),WORK(NYW),WORK(NCII) * ,WORK(NWT),WORK(NS),WORK(NFI),WORK(NRESID) * ,MSUM,MJP1,WORK(NW),IW,IERROR) 80 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC) RETURN END SUBROUTINE E02CAZ(M,N,KP1,L,X,Y,F,W,NX,A,NA,XMIN,XMAX,NUX,INUXP1, * NUY,INUYP1,WI,CI,YW,CII,WT,S,FI,RESID,MTOT,MJ1, * WORK,NWORK,IFAIL) C MARK 7 RELEASE. NAG COPYRIGHT 1978. C MARK 8C REVISED. IER-270 (OCT. 1980). C MARK 11.5(F77) REVISED. (SEPT 1985.) C C THIS SUBROUTINE OBTAINS A BIVARIATE CHEBYSHEV SERIES C APPROXIMATION OF DEGREE K IN X AND L IN Y TO A SET OF DATA C POINTS ARBITRARILY DISTRIBUTED ON LINES PARALLEL TO THE X C AXIS. C SUBROUTINE E02ADZ IS FIRST USED TO OBTAIN CHEBYSHEV SERIES C COEFFICIENTS OF THE WEIGHTED LEAST-SQUARES CURVE FIT OF C DEGREE K TO THE DATA ON EACH LINE IN TURN. THE SAME ROUTINE C IS THEN USED TO FIT THESE COEFFICIENTS AS POLYNOMIALS OF C DEGREE L IN Y, GIVING THE COEFFICIENTS OF THE BIVARIATE FIT. C THE FIT CAN BE CONSTRAINED TO CONTAIN GIVEN POLYNOMIAL C FACTORS IN X AND IN Y. C THE SUBROUTINE IS BASED ON THE METHOD GIVEN IN SECTIONS 5,6 C AND 8 OF CLENSHAW, C.W. AND HAYES, J.G., CURVE AND SURFACE C FITTING, J. INST. MATHS APPLICS, 1, PP.164-188, 1965. C C STARTED - 1978. C COMPLETED - 1978. C AUTHOR - GTA. C C .. Scalar Arguments .. INTEGER IFAIL, INUXP1, INUYP1, KP1, L, MJ1, MTOT, N, NA, * NWORK, NX C .. Array Arguments .. DOUBLE PRECISION A(NA), CI(N,KP1), CII(MJ1), F(NX), FI(N), * NUX(INUXP1), NUY(INUYP1), RESID(MTOT), S(MJ1), * W(NX), WI(N), WORK(NWORK), WT(MJ1), X(NX), * XMAX(N), XMIN(N), Y(N), YW(N,KP1) INTEGER M(N) C .. Local Scalars .. INTEGER I, IERROR, J, LP1, MFIRST, MLAST, T C .. External Subroutines .. EXTERNAL E02ADZ C .. Executable Statements .. LP1 = L + 1 MLAST = 0 IFAIL = 0 C C FIT THE DATA ON EACH OF THE N LINES Y = Y(I), I = 1,2,...N. C DO 60 I = 1, N IERROR = 1 MFIRST = MLAST + 1 MLAST = MLAST + M(I) CALL E02ADZ(MFIRST,MLAST,MTOT,KP1,1,0,1,X,F,W,XMIN(I),XMAX(I) * ,INUXP1,NUX,WORK(2*KP1+1),WORK(1) * ,CII,S,WT,RESID,IERROR) IF (IERROR.EQ.0) GO TO 20 IFAIL = IERROR + 1 IF (IFAIL.EQ.6) IFAIL = 1 GO TO 200 C C STORE THE COEFFICIENTS OF THE DEGREE K FIT TO DATA ON THE ITH C LINE C 20 DO 40 J = 1, KP1 CI(I,J) = CII(J) YW(I,J) = WT(J) 40 CONTINUE 60 CONTINUE C C FIT THE COEFFICIENTS A(J) ON EACH LINE AS A POLYNOMIAL IN Y. C IF (N.GT.1) GO TO 100 DO 80 J = 1, KP1 A(J) = 2.0D+0*CI(1,J) 80 CONTINUE RETURN 100 DO 180 J = 1, KP1 DO 120 I = 1, N FI(I) = CI(I,J) WI(I) = YW(I,J) 120 CONTINUE IERROR = 1 CALL E02ADZ(1,N,N,LP1,1,0,1,Y,FI,WI,Y(1),Y(N) * ,INUYP1,NUY,WORK(2*LP1+1),WORK(1) * ,CII,S,WT,RESID,IERROR) IF (IERROR.EQ.0) GO TO 140 IFAIL = IERROR + 1 IF (IFAIL.EQ.6) IFAIL = 1 GO TO 200 140 T = (J-1)*LP1 DO 160 I = 1, LP1 T = T + 1 A(T) = CII(I) 160 CONTINUE 180 CONTINUE 200 RETURN END SUBROUTINE E02CBF(MFIRST,MLAST,K,L,X,XMIN,XMAX,Y,YMIN,YMAX,FF,A, * NA,WORK,NWORK,IFAIL) C MARK 7 RELEASE. NAG COPYRIGHT 1978. C MARK 11.5(F77) REVISED. (SEPT 1985.) C C THIS SUBROUTINE EVALUATES A POLYNOMIAL OF DEGREE K AND L C RESPECTIVELY IN THE INDEPENDENT VARIABLES X AND Y. THE C POLYNOMIAL IS GIVEN IN DOUBLE CHEBYSHEV SERIES FORM C A(I,J) * TI(XCAP) * TJ(YCAP), C SUMMED OVER I = 0,1,...K AND J = 0,1,...L WITH THE CONVENTION C THAT TERMS WITH EITHER I OR J ZERO ARE HALVED AND THE TERM C WITH BOTH I AND J ZERO IS MULTIPLIED BY 0.25. HERE TI(XCAP) C IS THE CHEBYSHEV POLYNOMIAL OF THE FIRST KIND OF DEGREE I C WITH ARGUMENT XCAP=((X - XMIN) - (XMAX - X))/(XMAX - XMIN). C TJ(YCAP) IS DEFINED SIMILARLY. THE COEFFICIENT A(I,J) C SHOULD BE STORED IN ELEMENT (L + 1)*I + J + 1 OF THE SINGLE C DIMENSION ARRAY A. THE EVALUATION IS PERFORMED FOR A SINGLE C GIVEN VALUE OF Y WITH EACH X VALUE GIVEN IN X(R), FOR R = C MFIRST, MFIRST+1,....,MLAST. C C STARTED - 1978. C COMPLETED - 1978. C AUTHOR - GTA. C C .. Parameters .. CHARACTER*6 SRNAME PARAMETER (SRNAME='E02CBF') C .. Scalar Arguments .. DOUBLE PRECISION XMAX, XMIN, Y, YMAX, YMIN INTEGER IFAIL, K, L, MFIRST, MLAST, NA, NWORK C .. Array Arguments .. DOUBLE PRECISION A(NA), FF(MLAST), WORK(NWORK), X(MLAST) C .. Local Scalars .. DOUBLE PRECISION D, XCAP, YCAP INTEGER I, IERROR, KP1, LP1, M, R C .. Local Arrays .. CHARACTER*1 P01REC(1) C .. External Functions .. INTEGER P01ABF EXTERNAL P01ABF C .. External Subroutines .. EXTERNAL E02AEF C .. Executable Statements .. KP1 = K + 1 LP1 = L + 1 M = MLAST - MFIRST + 1 C C CHECK THAT THE INTEGER INPUT PARAMETERS HAVE REASONABLE C VALUES C IERROR = 1 IF (M.LE.0 .OR. K.LT.0 .OR. L.LT.0 .OR. NA.LT.KP1*LP1 .OR. * NWORK.LT.KP1) GO TO 80 C C CHECK THAT THE Y RANGE IS REASONABLE AND THAT THE GIVEN C VALUE OF Y IS NOT OUTSIDE IT C IERROR = 2 IF (YMIN.GE.YMAX .OR. Y.LT.YMIN .OR. Y.GT.YMAX) GO TO 80 D = XMAX - XMIN C C CHECK THAT THE X RANGE IS REASONABLE AND THAT NONE OF C THE GIVEN VALUES OF X IS OUTSIDE IT C IERROR = 3 IF (D.LE.0.0D+0) GO TO 80 DO 20 R = MFIRST, MLAST IF (X(R).LT.XMIN .OR. X(R).GT.XMAX) GO TO 80 20 CONTINUE C C CALCULATE YCAP, THE NORMALIZED VALUE OF Y C YCAP = ((Y-YMIN)-(YMAX-Y))/(YMAX-YMIN) IERROR = 1 R = -L C C EVALUATE THE COEFFICIENTS OF THE POLYNOMIAL FOR THE GIVEN Y C DO 40 I = 1, KP1 R = R + LP1 CALL E02AEF(LP1,A(R),YCAP,WORK(I),IERROR) IERROR = IERROR + 1 IF (IERROR.NE.1) GO TO 80 40 CONTINUE C C EVALUATE THE POLYNOMAL AT THE GIVEN X VALUES C DO 60 R = MFIRST, MLAST XCAP = ((X(R)-XMIN)-(XMAX-X(R)))/D IERROR = 1 CALL E02AEF(KP1,WORK,XCAP,FF(R),IERROR) IF (IERROR.EQ.0) GO TO 60 IERROR = 3 GO TO 80 60 CONTINUE 80 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC) RETURN END CX DOUBLE PRECISION FUNCTION G05DDF(A,B) C MARK 6 RELEASE NAG COPYRIGHT 1976 C MARK 11 REVISED. IER-441 (FEB 1984). C MARK 11.5(F77) REVISED. (SEPT 1985.) C WRITTEN BY N.M.MACLAREN C UNIVERSITY OF CAMBRIDGE COMPUTER LABORATORY C THIS RETURNS A REAL NUMBER, NORMALLY DISTRIBUTED (GAUSSIAN) C WITH MEAN A AND STANDARD DEVIATION B. C THE METHOD USED IS A MODIFICATION OF CACM ALGORITHM 488 C (BY R.P.BRENT). C THE FOLLOWING CAN BE EXTENDED IF A TRUNCATION PROBABILITY C OF 1.0E-12 IS REGARDED AS UNSATISFACTORY. C THE CODE ASSUMES THAT THE CONTENTS OF COMMON BLOCK /BG05CA/ C ARE SAVED BETWEEN CALLS OF THIS ROUTINE AND CALLS OF OTHER C G05 ROUTINES THAT REFERENCE IT. A SAVE STATEMENT C ENSURES THAT THIS IS SO. C .. Scalar Arguments .. CX DOUBLE PRECISION A, B C .. Scalars in Common .. CX DOUBLE PRECISION STORE1, STORE2 C .. Local Scalars .. CX DOUBLE PRECISION HALF, ONE, T, U, V, W CX INTEGER N C .. Local Arrays .. CX DOUBLE PRECISION D(41) C .. External Functions .. CX DOUBLE PRECISION G05CAF CX EXTERNAL G05CAF C .. Common blocks .. CX COMMON /BG05CA/STORE1, STORE2 C .. Save statement .. CX SAVE /BG05CA/ C .. Data statements .. CX DATA ONE/1.0D0/, HALF/0.5D0/ CX DATA D(1), D(2), D(3), D(4), D(5), CX * D(6), D(7), D(8), D(9), D(10), CX * D(11), D(12), D(13), D(14)/0.0D0, CX * 0.674489750196082D0, CX * 1.150349380376008D0, CX * 1.534120544352546D0, CX * 1.862731867421652D0, CX * 2.153874694061456D0, CX * 2.417559016236505D0, CX * 2.660067468617460D0, CX * 2.885634912426757D0, CX * 3.097269078198785D0, CX * 3.297193345691964D0, CX * 3.487104104114431D0, CX * 3.668329285121323D0, CX * 3.841930685501911D0/ CX DATA D(15), D(16), D(17), D(18), CX * D(19), D(20), D(21), D(22), CX * D(23), D(24), D(25), D(26), CX * D(27)/4.008772594168585D0, CX * 4.169569323349106D0, CX * 4.324919040826046D0, CX * 4.475328424654204D0, CX * 4.621231001499247D0, CX * 4.763001034267814D0, CX * 4.900964207963193D0, CX * 5.035405969463927D0, CX * 5.166578119728753D0, CX * 5.294704084854598D0, CX * 5.419983174916868D0, CX * 5.542594057802940D0, CX * 5.662697617459439D0/ CX DATA D(28), D(29), D(30), D(31), CX * D(32), D(33), D(34), D(35), CX * D(36), D(37), D(38), D(39), CX * D(40)/5.780439324478935D0, CX * 5.895951216739571D0, CX * 6.009353565530745D0, CX * 6.120756285971941D0, CX * 6.230260137989044D0, CX * 6.337957754553790D0, CX * 6.443934526538564D0, CX * 6.548269367831731D0, CX * 6.651035379893011D0, CX * 6.752300431407015D0, CX * 6.852127665896068D0, CX * 6.950575947916750D0, CX * 7.047700256664409D0/ CX DATA D(41)/7.143552034352190D0/ C .. Executable Statements .. C FIRST CALL TO G05CAF COMES BEFORE FIRST REFERENCE TO STORE1 C TO MAKE SURE THAT STORE1 IS INITIALIZED CX V = G05CAF(0.0D0) CX U = STORE1 CX DO 20 N = 1, 39 CX IF (U.GT.HALF) GO TO 40 CX U = U + U CX 20 CONTINUE CX N = 40 CX 40 T = D(N) CX U = V CX 60 W = (D(N+1)-T)*U CX V = W*(W*HALF+T) CX 80 U = G05CAF(0.0D0) CX IF (V.LE.U) GO TO 100 CX V = G05CAF(0.0D0) CX IF (U.GT.V) GO TO 80 CX U = (V-U)/(ONE-U) CX GO TO 60 CX 100 U = (U-V)/(ONE-V) CX IF (U.GT.HALF) GO TO 120 CX STORE1 = U + U CX G05DDF = A + B*(W+T) CX RETURN CX 120 STORE1 = U + U - ONE CX G05DDF = A - B*(W+T) CX RETURN CX END INTEGER FUNCTION P01AAF(IFAIL,ERROR,SRNAME) C MARK 1 RELEASE. NAG COPYRIGHT 1971 C MARK 3 REVISED C MARK 4A REVISED, IER-45 C MARK 4.5 REVISED C MARK 7 REVISED (DEC 1978) C MARK 11 REVISED (FEB 1984) C MARK 11.5(F77) REVISED. (SEPT 1985.) C RETURNS THE VALUE OF ERROR OR TERMINATES THE PROGRAM. C TEST IF NO ERROR DETECTED C .. Scalar Arguments .. DOUBLE PRECISION SRNAME INTEGER ERROR, IFAIL C .. Local Scalars .. INTEGER NOUT C .. External Subroutines .. EXTERNAL X04AAF C .. Intrinsic Functions .. INTRINSIC MOD C .. Executable Statements .. IF (ERROR.EQ.0) GO TO 40 C DETERMINE OUTPUT UNIT FOR MESSAGE CALL X04AAF(0,NOUT) C TEST FOR SOFT FAILURE IF (MOD(IFAIL,10).EQ.1) GO TO 20 C HARD FAILURE WRITE (NOUT,FMT=99999) SRNAME, ERROR C ******************** IMPLEMENTATION NOTE ******************** C THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN C IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR C ABORT THE PROGRAM. C ************************************************************* STOP C SOFT FAIL C TEST IF ERROR MESSAGES SUPPRESSED 20 IF (MOD(IFAIL/10,10).EQ.0) GO TO 40 WRITE (NOUT,FMT=99999) SRNAME, ERROR 40 P01AAF = ERROR RETURN C 99999 FORMAT (/' ERROR DETECTED BY NAG LIBRARY ROUTINE ',A8,' - IFAIL ', * '= ',I5,//) END INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C P01ABF is the error-handling routine for the NAG Library. C C P01ABF either returns the value of IERROR through the routine C name (soft failure), or terminates execution of the program C (hard failure). Diagnostic messages may be output. C C If IERROR = 0 (successful exit from the calling routine), C the value 0 is returned through the routine name, and no C message is output C C If IERROR is non-zero (abnormal exit from the calling routine), C the action taken depends on the value of IFAIL. C C IFAIL = 1: soft failure, silent exit (i.e. no messages are output) C IFAIL =-1: soft failure, noisy exit (i.e. messages are output) C IFAIL = 0: hard failure, noisy exit C C For compatibility with certain routines included before Mark 12 C P01ABF also allows an alternative specification of IFAIL in which C it is regarded as a decimal integer with least significant digits C cba. Then C C a = 0: hard failure a = 1: soft failure C b = 0: silent exit b = 1: noisy exit C C except that hard failure now always implies a noisy exit. C C S.Hammarling, M.P.Hooper and J.J.du Croz, NAG Central Office. C C .. Scalar Arguments .. INTEGER IERROR, IFAIL, NREC CHARACTER*(*) SRNAME C .. Array Arguments .. CHARACTER*(*) REC(*) C .. Local Scalars .. INTEGER I, NERR CHARACTER*72 MESS C .. External Subroutines .. EXTERNAL X04AAF, X04BAF C .. Intrinsic Functions .. INTRINSIC ABS, MOD C .. Executable Statements .. IF (IERROR.NE.0) THEN C Abnormal exit from calling routine IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. * (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN C Noisy exit CALL X04AAF(0,NERR) DO 20 I = 1, NREC CALL X04BAF(NERR,REC(I)) 20 CONTINUE WRITE (MESS,FMT=99999) SRNAME, IERROR CALL X04BAF(NERR,MESS) IF (ABS(MOD(IFAIL,10)).NE.1) THEN C Hard failure CALL X04BAF(NERR, * ' ** NAG hard failure - execution terminated' * ) STOP ELSE C Soft failure CALL X04BAF(NERR, * ' ** NAG soft failure - control returned') END IF END IF END IF P01ABF = IERROR RETURN C 99999 FORMAT (' ** ABNORMAL EXIT from NAG Library routine ',A,': IFAIL', * ' =',I6) END DOUBLE PRECISION FUNCTION X02AAF(X) C MARK 12 RE-ISSUE. NAG COPYRIGHT 1986. C C RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE. C RETURNS B**(1-P) OTHERWISE C I.E. RETURNS THE SAME VALUE AS X02AJF C C .. Scalar Arguments .. DOUBLE PRECISION X C .. External Functions .. DOUBLE PRECISION X02AJF EXTERNAL X02AJF C .. Executable Statements .. X02AAF = X02AJF() RETURN END DOUBLE PRECISION FUNCTION X02AJF() C MARK 12 RELEASE. NAG COPYRIGHT 1986. C C RETURNS (1/2)*B**(1-P) IF ROUNDS IS .TRUE. C RETURNS B**(1-P) OTHERWISE C C .. Executable Statements .. C IN THEORY THIS SHOULD BE 2.0**(-56) BUT 2.0**(-55) HAS BEEN FOUND C TO BE MORE PRACTICAL IN THE PAST. X02AJF = 2.0D0**(-55) RETURN END SUBROUTINE X04AAF(I,NERR) C MARK 7 RELEASE. NAG COPYRIGHT 1978 C MARK 7C REVISED IER-190 (MAY 1979) C MARK 11.5(F77) REVISED. (SEPT 1985.) C IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER C (STORED IN NERR1). C IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO C VALUE SPECIFIED BY NERR. C C .. Scalar Arguments .. INTEGER I, NERR C .. Local Scalars .. INTEGER NERR1 C .. Save statement .. SAVE NERR1 C .. Data statements .. DATA NERR1/6/ C .. Executable Statements .. IF (I.EQ.0) NERR = NERR1 IF (I.EQ.1) NERR1 = NERR RETURN END SUBROUTINE X04BAF(NOUT,REC) C MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986. C C X04BAF writes the contents of REC to the unit defined by NOUT. C C Trailing blanks are not output, except that if REC is entirely C blank, a single blank character is output. C If NOUT.lt.0, i.e. if NOUT is not a valid Fortran unit identifier, C then no output occurs. C C .. Scalar Arguments .. INTEGER NOUT CHARACTER*(*) REC C .. Local Scalars .. INTEGER I C .. Intrinsic Functions .. INTRINSIC LEN C .. Executable Statements .. IF (NOUT.GE.0) THEN C Remove trailing blanks DO 20 I = LEN(REC), 2, -1 IF (REC(I:I).NE.' ') GO TO 40 20 CONTINUE C Write record to external file 40 WRITE (NOUT,FMT=99999) REC(1:I) END IF RETURN C 99999 FORMAT (A) END