C ALGORITHM 433 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 10, C P. 914. SUBROUTINE INTRPL(IU,L,X,Y,N,U,V) INTR 10 C INTERPOLATION OF A SINGLE-VALUED FUNCTION C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION C GIVEN AS ORDINATES OF INPUT DATA POINTS IN AN X-Y PLANE C AND FOR A GIVEN SET OF X VALUES (ABSCISSAS), THE VALUES OF C A SINGLE-VALUED FUNCTION Y = Y(X). C THE INPUT PARAMETERS ARE C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT C L = NUMBER OF INPUT DATA POINTS C (MUST BE 2 OR GREATER) C X = ARRAY OF DIMENSION L STORING THE X VALUES C (ABSCISSAS) OF INPUT DATA POINTS C (IN ASCENDING ORDER) C Y = ARRAY OF DIMENSION L STORING THE Y VALUES C (ORDINATES) OF INPUT DATA POINTS C N = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE C Y VALUE (ORDINATE) IS DESIRED C (MUST BE 1 OR GREATER) C U = ARRAY OF DIMENSION N STORING THE X VALUES C (ABSCISSAS) OF DESIRED POINTS C THE OUTPUT PARAMETER IS C V = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Y C VALUES (ORDINATES) ARE TO BE DISPLAYED C DECLARATION STATEMENTS DIMENSION X(L),Y(L),U(N),V(N) EQUIVALENCE (P0,X3),(Q0,Y3),(Q1,T3) REAL M1,M2,M3,M4,M5 EQUIVALENCE (UK,DX),(IMN,X2,A1,M1),(IMX,X5,A5,M5), 1 (J,SW,SA),(Y2,W2,W4,Q2),(Y5,W3,Q3) C PRELIMINARY PROCESSING 10 L0=L LM1=L0-1 LM2=LM1-1 LP1=L0+1 N0=N IF(LM2.LT.0) GO TO 90 IF(N0.LE.0) GO TO 91 DO 11 I=2,L0 IF(X(I-1)-X(I)) 11,95,96 11 CONTINUE IPV=0 C MAIN DO-LOOP DO 80 K=1,N0 UK=U(K) C ROUTINE TO LOCATE THE DESIRED POINT 20 IF(LM2.EQ.0) GO TO 27 IF(UK.GE.X(L0)) GO TO 26 IF(UK.LT.X(1)) GO TO 25 IMN=2 IMX=L0 21 I=(IMN+IMX)/2 IF(UK.GE.X(I)) GO TO 23 22 IMX=I GO TO 24 23 IMN=I+1 24 IF(IMX.GT.IMN) GO TO 21 I=IMX GO TO 30 25 I=1 GO TO 30 26 I=LP1 GO TO 30 27 I=2 C CHECK IF I = IPV 30 IF(I.EQ.IPV) GO TO 70 IPV=I C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND C TO ESTIMATE THEM IF NECESSARY 40 J=I IF(J.EQ.1) J=2 IF(J.EQ.LP1) J=L0 X3=X(J-1) Y3=Y(J-1) X4=X(J) Y4=Y(J) A3=X4-X3 M3=(Y4-Y3)/A3 IF(LM2.EQ.0) GO TO 43 IF(J.EQ.2) GO TO 41 X2=X(J-2) Y2=Y(J-2) A2=X3-X2 M2=(Y3-Y2)/A2 IF(J.EQ.L0) GO TO 42 41 X5=X(J+1) Y5=Y(J+1) A4=X5-X4 M4=(Y5-Y4)/A4 IF(J.EQ.2) M2=M3+M3-M4 GO TO 45 42 M4=M3+M3-M2 GO TO 45 43 M2=M3 M4=M3 45 IF(J.LE.3) GO TO 46 A1=X2-X(J-3) M1=(Y2-Y(J-3))/A1 GO TO 47 46 M1=M2+M2-M3 47 IF(J.GE.LM1) GO TO 48 A5=X(J+2)-X5 M5=(Y(J+2)-Y5)/A5 GO TO 50 48 M5=M4+M4-M3 C NUMERICAL DIFFERENTIATION 50 IF(I.EQ.LP1) GO TO 52 W2=ABS(M4-M3) W3=ABS(M2-M1) SW=W2+W3 IF(SW.NE.0.0) GO TO 51 W2=0.5 W3=0.5 SW=1.0 51 T3=(W2*M2+W3*M3)/SW IF(I.EQ.1) GO TO 54 52 W3=ABS(M5-M4) W4=ABS(M3-M2) SW=W3+W4 IF(SW.NE.0.0) GO TO 53 W3=0.5 W4=0.5 SW=1.0 53 T4=(W3*M3+W4*M4)/SW IF(I.NE.LP1) GO TO 60 T3=T4 SA=A2+A3 T4=0.5*(M4+M5-A2*(A2-A3)*(M2-M3)/(SA*SA)) X3=X4 Y3=Y4 A3=A2 M3=M4 GO TO 60 54 T4=T3 SA=A3+A4 T3=0.5*(M1+M2-A4*(A3-A4)*(M3-M4)/(SA*SA)) X3=X3-A4 Y3=Y3-M2*A4 A3=A4 M3=M2 C DETERMINATION OF THE COEFFICIENTS 60 Q2=(2.0*(M3-T3)+M3-T4)/A3 Q3=(-M3-M3+T3+T4)/(A3*A3) C COMPUTATION OF THE POLYNOMIAL 70 DX=UK-P0 80 V(K)=Q0+DX*(Q1+DX*(Q2+DX*Q3)) RETURN C ERROR EXIT 90 WRITE (IU,2090) GO TO 99 91 WRITE (IU,2091) GO TO 99 95 WRITE (IU,2095) GO TO 97 96 WRITE (IU,2096) 97 WRITE (IU,2097) I,X(I) 99 WRITE (IU,2099) L0,N0 RETURN C FORMAT STATEMENTS 2090 FORMAT(1X/22H *** L = 1 OR LESS./) 2091 FORMAT(1X/22H *** N = 0 OR LESS./) 2095 FORMAT(1X/27H *** IDENTICAL X VALUES./) 2096 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./) 2097 FORMAT(6H I =,I7,10X,6HX(I) =,E12.3) 2099 FORMAT(6H L =,I7,10X,3HN =,I7/ 1 36H ERROR DETECTED IN ROUTINE INTRPL) END SUBROUTINE CRVFIT(IU,MD,L,X,Y,M,N,U,V) CRVF1670 C SMOOTH CURVE FITTING C THIS SUBROUTINE FITS A SMOOTH CURVE TO A GIVEN SET OF IN- C PUT DATA POINTS IN AN X-Y PLANE. IT INTERPOLATES POINTS C IN EACH INTERVAL BETWEEN A PAIR OF DATA POINTS AND GENER- C ATES A SET OF OUTPUT POINTS CONSISTING OF THE INPUT DATA C POINTS AND THE INTERPOLATED POINTS. IT CAN PROCESS EITHER C A SINGLE-VALUED FUNCTION OR A MULTIPLE-VALUED FUNCTION. C THE INPUT PARAMETERS ARE C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT C MD = MODE OF THE CURVE (MUST BE 1 OR 2) C = 1 FOR A SINGLE-VALUED FUNCTION C = 2 FOR A MULTIPLE-VALUED FUNCTION C L = NUMBER OF INPUT DATA POINTS C (MUST BE 2 OR GREATER) C X = ARRAY OF DIMENSION L STORING THE ABSCISSAS OF C INPUT DATA POINTS (IN ASCENDING OR DESCENDING C ORDER FOR MD = 1) C Y = ARRAY OF DIMENSION L STORING THE ORDINATES OF C INPUT DATA POINTS C M = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF C INPUT DATA POINTS (MUST BE 2 OR GREATER) C N = NUMBER OF OUTPUT POINTS C = (L-1)*M+1 C THE OUTPUT PARAMETERS ARE C U = ARRAY OF DIMENSION N WHERE THE ABSCISSAS OF C OUTPUT POINTS ARE TO BE DISPLAYED C V = ARRAY OF DIMENSION N WHERE THE ORDINATES OF C OUTPUT POINTS ARE TO BE DISPLAYED C DECLARATION STATEMENTS DIMENSION X(L),Y(L),U(N),V(N) EQUIVALENCE (M1,B1),(M2,B2),(M3,B3),(M4,B4), 1 (X2,P0),(Y2,Q0),(T2,Q1) REAL M1,M2,M3,M4 EQUIVALENCE (W2,Q2),(W3,Q3),(A1,P2),(B1,P3), 1 (A2,DZ),(SW,R,Z) C PRELIMINARY PROCESSING 10 MD0=MD MDM1=MD0-1 L0=L LM1=L0-1 M0=M MM1=M0-1 N0=N IF(MD0.LE.0) GO TO 90 IF(MD0.GE.3) GO TO 90 IF(LM1.LE.0) GO TO 91 IF(MM1.LE.0) GO TO 92 IF(N0.NE.LM1*M0+1) GO TO 93 GO TO (11,16), MD0 11 I=2 IF(X(1)-X(2)) 12,95,14 12 DO 13 I=3,L0 IF(X(I-1)-X(I)) 13,95,96 13 CONTINUE GO TO 18 14 DO 15 I=3,L0 IF(X(I-1)-X(I)) 96,95,15 15 CONTINUE GO TO 18 16 DO 17 I=2,L0 IF(X(I-1).NE.X(I)) GO TO 17 IF(Y(I-1).EQ.Y(I)) GO TO 97 17 CONTINUE 18 K=N0+M0 I=L0+1 DO 19 J=1,L0 K=K-M0 I=I-1 U(K)=X(I) 19 V(K)=Y(I) RM=M0 RM=1.0/RM C MAIN DO-LOOP 20 K5=M0+1 DO 80 I=1,L0 C ROUTINES TO PICK UP NECESSARY X AND Y VALUES AND C TO ESTIMATE THEM IF NECESSARY IF(I.GT.1) GO TO 40 30 X3=U(1) Y3=V(1) X4=U(M0+1) Y4=V(M0+1) A3=X4-X3 B3=Y4-Y3 IF(MDM1.EQ.0) M3=B3/A3 IF(L0.NE.2) GO TO 41 A4=A3 B4=B3 31 GO TO (33,32), MD0 32 A2=A3+A3-A4 A1=A2+A2-A3 33 B2=B3+B3-B4 B1=B2+B2-B3 GO TO (51,56), MD0 40 X2=X3 Y2=Y3 X3=X4 Y3=Y4 X4=X5 Y4=Y5 A1=A2 B1=B2 A2=A3 B2=B3 A3=A4 B3=B4 IF(I.GE.LM1) GO TO 42 41 K5=K5+M0 X5=U(K5) Y5=V(K5) A4=X5-X4 B4=Y5-Y4 IF(MDM1.EQ.0) M4=B4/A4 GO TO 43 42 IF(MDM1.NE.0) A4=A3+A3-A2 B4=B3+B3-B2 43 IF(I.EQ.1) GO TO 31 GO TO (50,55), MD0 C NUMERICAL DIFFERENTIATION 50 T2=T3 51 W2=ABS(M4-M3) W3=ABS(M2-M1) SW=W2+W3 IF(SW.NE.0.0) GO TO 52 W2=0.5 W3=0.5 SW=1.0 52 T3=(W2*M2+W3*M3)/SW IF(I-1) 80,80,60 55 COS2=COS3 SIN2=SIN3 56 W2=ABS(A3*B4-A4*B3) W3=ABS(A1*B2-A2*B1) IF(W2+W3.NE.0.0) GO TO 57 W2=SQRT(A3*A3+B3*B3) W3=SQRT(A2*A2+B2*B2) 57 COS3=W2*A2+W3*A3 SIN3=W2*B2+W3*B3 R=COS3*COS3+SIN3*SIN3 IF(R.EQ.0.0) GO TO 58 R=SQRT(R) COS3=COS3/R SIN3=SIN3/R 58 IF(I-1) 80,80,65 C DETERMINATION OF THE COEFFICIENTS 60 Q2=(2.0*(M2-T2)+M2-T3)/A2 Q3=(-M2-M2+T2+T3)/(A2*A2) GO TO 70 65 R=SQRT(A2*A2+B2*B2) P1=R*COS2 P2=3.0*A2-R*(COS2+COS2+COS3) P3=A2-P1-P2 Q1=R*SIN2 Q2=3.0*B2-R*(SIN2+SIN2+SIN3) Q3=B2-Q1-Q2 GO TO 75 C COMPUTATION OF THE POLYNOMIALS 70 DZ=A2*RM Z=0.0 DO 71 J=1,MM1 K=K+1 Z=Z+DZ U(K)=P0+Z 71 V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3)) GO TO 79 75 Z=0.0 DO 76 J=1,MM1 K=K+1 Z=Z+RM U(K)=P0+Z*(P1+Z*(P2+Z*P3)) 76 V(K)=Q0+Z*(Q1+Z*(Q2+Z*Q3)) 79 K=K+1 80 CONTINUE RETURN C ERROR EXIT 90 WRITE (IU,2090) GO TO 99 91 WRITE (IU,2091) GO TO 99 92 WRITE (IU,2092) GO TO 99 93 WRITE (IU,2093) GO TO 99 95 WRITE (IU,2095) GO TO 98 96 WRITE (IU,2096) GO TO 98 97 WRITE (IU,2097) 98 WRITE (IU,2098) I,X(I),Y(I) 99 WRITE (IU,2099) MD0,L0,M0,N0 RETURN C FORMAT STATEMENTS 2090 FORMAT(1X/31H *** MD OUT OF PROPER RANGE./) 2091 FORMAT(1X/22H *** L = 1 OR LESS./) 2092 FORMAT(1X/22H *** M = 1 OR LESS./) 2093 FORMAT(1X/25H *** IMPROPER N VALUE./) 2095 FORMAT(1X/27H *** IDENTICAL X VALUES./) 2096 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./) 2097 FORMAT(1X/33H *** IDENTICAL X AND Y VALUES./) 2098 FORMAT(7H I =,I4,10X,6HX(I) =,E12.3, 1 10X,6HY(I) =,E12.3) 2099 FORMAT(7H MD =,I4,8X,3HL =,I5,8X, 1 3HM =,I5,8X,3HN =,I5/ 2 36H ERROR DETECTED IN ROUTINE CRVFIT) END