Commit ca58ed3a authored by agebhard's avatar agebhard
Browse files

deleted

parent 4f872db6
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
C PROGRAM TTIDBS(OUTPUT,TAPE6=OUTPUT) ID000070
C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO- ID000080
C GRAM PACKAGE. ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS ID000090
C ARE EXPECTED TO BE ZERO. ID000100
C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS ID000110
C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, ID000120
C THEREFORE, SYSTEM DEPENDENT. ID000130
C DECLARATION STATEMENTS ID000140
DIMENSION XD(30),YD(30),ZD(30), ID000150
1 XI(6),YI(5),ZI(6,5), ID000160
2 ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5), ID000170
3 IWK(1030),WK(240) ID000180
DATA NCP/4/ ID000190
DATA NDP/30/ ID000200
DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6), ID000210
1 XD(7), XD(8), XD(9), XD(10),XD(11),XD(12), ID000220
2 XD(13),XD(14),XD(15),XD(16),XD(17),XD(18), ID000230
3 XD(19),XD(20),XD(21),XD(22),XD(23),XD(24), ID000240
4 XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/ ID000250
5 11.16, 24.20, 19.85, 10.35, 19.72, 0.00, ID000260
6 20.87, 19.99, 10.28, 4.51, 0.00, 16.70, ID000270
7 6.08, 25.00, 14.90, 0.00, 9.66, 5.22, ID000280
8 11.77, 15.10, 25.00, 25.00, 14.59, 15.20, ID000290
9 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/ ID000300
DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6), ID000310
1 YD(7), YD(8), YD(9), YD(10),YD(11),YD(12), ID000320
2 YD(13),YD(14),YD(15),YD(16),YD(17),YD(18), ID000330
3 YD(19),YD(20),YD(21),YD(22),YD(23),YD(24), ID000340
4 YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/ ID000350
5 1.24, 16.23, 10.72, 4.11, 1.39, 20.00, ID000360
6 20.00, 4.62, 15.16, 20.00, 4.48, 19.65, ID000370
7 4.58, 11.87, 3.12, 0.00, 20.00, 14.66, ID000380
8 10.47, 17.19, 3.87, 0.00, 8.71, 0.00, ID000390
9 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/ ID000400
DATA ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6), ID000410
1 ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12), ID000420
2 ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18), ID000430
3 ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24), ID000440
4 ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/ ID000450
5 22.15, 2.83, 7.97, 22.33, 16.83, 34.60, ID000460
6 5.74, 14.72, 21.59, 15.61, 61.77, 6.31, ID000470
7 35.74, 4.40, 21.70, 58.20, 4.73, 40.36, ID000480
8 13.62, 12.57, 8.74, 12.00, 14.81, 21.60, ID000490
9 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/ ID000500
DATA NXI/6/, NYI/5/ ID000510
DATA XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/ ID000520
1 0.00, 5.00, 10.00, 15.00, 20.00, 25.00/ ID000530
DATA YI(1), YI(2), YI(3), YI(4), YI(5)/ ID000540
1 0.00, 5.00, 10.00, 15.00, 20.00/ ID000550
DATA ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1), ID000560
1 ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2), ID000570
2 ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3), ID000580
3 ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4), ID000590
4 ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/ ID000600
5 58.20, 39.55, 26.90, 21.71, 17.68, 12.00, ID000610
6 61.58, 39.39, 22.04, 21.29, 14.36, 8.04, ID000620
7 59.18, 27.39, 16.78, 13.25, 8.59, 5.36, ID000630
8 52.82, 40.27, 22.76, 16.61, 7.40, 2.88, ID000640
9 34.60, 14.05, 4.12, 3.17, 6.31, 0.60/ ID000650
DATA LUN/6/ ID000660
C CALCULATION ID000670
10 MD=1 ID000680
DO 12 IYI=1,NYI ID000690
DO 11 IXI=1,NXI ID000700
IF(IXI.NE.1.OR.IYI.NE.1) MD=2 ID000710
CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI), ID000720
1 ZI1(IXI,IYI),IWK,WK) ID000730
11 CONTINUE ID000740
12 CONTINUE ID000750
15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK) ID000760
DO 17 IYI=1,NYI ID000770
DO 16 IXI=1,NXI ID000780
DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI)) ID000790
DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI)) ID000800
16 CONTINUE ID000810
17 CONTINUE ID000820
C PRINTING OF INPUT DATA ID000830
20 WRITE (LUN,2020) NDP ID000840
DO 23 IDP=1,NDP ID000850
IF(MOD(IDP,5).EQ.1) WRITE (LUN,2021) ID000860
WRITE (LUN,2022) IDP,XD(IDP),YD(IDP),ZD(IDP) ID000870
23 CONTINUE ID000880
C PRINTING OF OUTPUT RESULTS ID000890
30 WRITE (LUN,2030) ID000900
WRITE (LUN,2031) YI ID000910
DO 33 IXI=1,NXI ID000920
WRITE (LUN,2032) XI(IXI),(ZI1(IXI,IYI),IYI=1,NYI) ID000930
33 CONTINUE ID000940
40 WRITE (LUN,2040) ID000950
WRITE (LUN,2031) YI ID000960
DO 43 IXI=1,NXI ID000970
WRITE (LUN,2032) XI(IXI),(DZI1(IXI,IYI),IYI=1,NYI) ID000980
43 CONTINUE ID000990
50 WRITE (LUN,2050) ID001000
WRITE (LUN,2031) YI ID001010
DO 53 IXI=1,NXI ID001020
WRITE (LUN,2032) XI(IXI),(ZI2(IXI,IYI),IYI=1,NYI) ID001030
53 CONTINUE ID001040
60 WRITE (LUN,2060) ID001050
WRITE (LUN,2031) YI ID001060
DO 63 IXI=1,NXI ID001070
WRITE (LUN,2032) XI(IXI),(DZI2(IXI,IYI),IYI=1,NYI) ID001080
63 CONTINUE ID001090
STOP ID001100
C FORMAT STATEMENTS ID001110
2020 FORMAT(1H1,6HTTIDBS/////3X,10HINPUT DATA,8X,5HNDP =,I3/// ID001120
1 30H I XD YD ZD /) ID001130
2021 FORMAT(1X) ID001140
2022 FORMAT(5X,I2,2X,3F7.2) ID001150
2030 FORMAT(1H1,6HTTIDBS/////3X,17HIDBVIP SUBROUTINE/// ID001160
1 26X,10HZI1(XI,YI)) ID001170
2031 FORMAT(7X,2HXI,4X,3HYI=/12X,5F7.2/) ID001180
2032 FORMAT(1X/1X,F9.2,2X,5F7.2) ID001190
2040 FORMAT(1X/////3X,10HDIFFERENCE/// ID001200
1 25X,11HDZI1(XI,YI)) ID001210
2050 FORMAT(1H1,6HTTIDBS/////3X,17HIDSFFT SUBROUTINE/// ID001220
1 26X,10HZI2(XI,YI)) ID001230
2060 FORMAT(1X/////3X,10HDIFFERENCE/// ID001240
1 25X,11HDZI2(XI,YI)) ID001250
END ID001260
SUBROUTINE IDBVIP(MD,NCP,NDP,XD,YD,ZD,NIP,XI,YI,ZI, ID001340
1 IWK,WK)
C THIS SUBROUTINE PERFORMS BIVARIATE INTERPOLATION WHEN THE PRO-
C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY
C DISTRIBUTED IN THE PLANE.
C THE INPUT PARAMETERS ARE
C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3),
C = 1 FOR NEW NCP AND/OR NEW XD-YD,
C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI,
C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT
C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP),
C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER),
C XD = ARRAY OF DIMENSION NDP CONTAINING THE X
C COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y
C COORDINATES OF THE DATA POINTS,
C ZD = ARRAY OF DIMENSION NDP CONTAINING THE Z
C COORDINATES OF THE DATA POINTS,
C NIP = NUMBER OF OUTPUT POINTS AT WHICH INTERPOLATION
C IS TO BE PERFORMED (MUST BE 1 OR GREATER),
C XI = ARRAY OF DIMENSION NIP CONTAINING THE X
C COORDINATES OF THE OUTPUT POINTS,
C YI = ARRAY OF DIMENSION NIP CONTAINING THE Y
C COORDINATES OF THE OUTPUT POINTS.
C THE OUTPUT PARAMETER IS
C ZI = ARRAY OF DIMENSION NIP WHERE INTERPOLATED Z
C VALUES ARE TO BE STORED.
C THE OTHER PARAMETERS ARE
C IWK = INTEGER ARRAY OF DIMENSION
C MAX0(31,27+NCP)*NDP+NIP
C USED INTERNALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A
C WORK AREA.
C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW
C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND
C YD ARRAYS MUST BE MADE WITH MD=1. THE CALL WITH MD=2 MUST BE
C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND
C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS. THE CALL WITH
C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP,
C AND NIP VALUES AND WITH THE SAME CONTENTS OF THE XD, YD, XI,
C AND YI ARRAYS. BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS
C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED.
C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM-
C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C THIS SUBROUTINE CALLS THE IDCLDP, IDLCTN, IDPDRV, IDPTIP, AND
C IDTANG SUBROUTINES.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),XI(1000),YI(1000),
1 ZI(1000),IWK(4100),WK(800)
COMMON/IDLC/NIT
COMMON/IDPI/ITPV
DATA LUN/6/
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES.
C (FOR MD=1,2,3)
10 MD0=MD
NCP0=NCP
NDP0=NDP
NIP0=NIP
C ERROR CHECK. (FOR MD=1,2,3)
20 IF(MD0.LT.1.OR.MD0.GT.3) GO TO 90
IF(NCP0.LT.2.OR.NCP0.GE.NDP0) GO TO 90
IF(NDP0.LT.4) GO TO 90
IF(NIP0.LT.1) GO TO 90
IF(MD0.GE.2) GO TO 21
IWK(1)=NCP0
IWK(2)=NDP0
GO TO 22
21 NCPPV=IWK(1)
NDPPV=IWK(2)
IF(NCP0.NE.NCPPV) GO TO 90
IF(NDP0.NE.NDPPV) GO TO 90
22 IF(MD0.GE.3) GO TO 23
IWK(3)=NIP
GO TO 30
23 NIPPV=IWK(3)
IF(NIP0.NE.NIPPV) GO TO 90
C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY. (FOR MD=1,2,3)
30 JWIPT=16
JWIWL=6*NDP0+1
JWIWK=JWIWL
JWIPL=24*NDP0+1
JWIWP=30*NDP0+1
JWIPC=27*NDP0+1
JWIT0=MAX0(31,27+NCP0)*NDP0
C TRIANGULATES THE X-Y PLANE. (FOR MD=1)
40 IF(MD0.GT.1) GO TO 50
CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
1 IWK(JWIWL),IWK(JWIWP),WK)
IWK(5)=NT
IWK(6)=NL
IF(NT.EQ.0) RETURN
C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1)
50 IF(MD0.GT.1) GO TO 60
CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC))
IF(IWK(JWIPC).EQ.0) RETURN
C LOCATES ALL POINTS AT WHICH INTERPOLATION IS TO BE PERFORMED.
C (FOR MD=1,2)
60 IF(MD0.EQ.3) GO TO 70
NIT=0
JWIT=JWIT0
DO 61 IIP=1,NIP0
JWIT=JWIT+1
CALL IDLCTN(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),
1 XI(IIP),YI(IIP),IWK(JWIT),IWK(JWIWK),WK)
61 CONTINUE
C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS.
C (FOR MD=1,2,3)
70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK)
C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3)
80 ITPV=0
JWIT=JWIT0
DO 81 IIP=1,NIP0
JWIT=JWIT+1
CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK,
1 IWK(JWIT),XI(IIP),YI(IIP),ZI(IIP))
81 CONTINUE
RETURN
C ERROR EXIT
90 WRITE (LUN,2090) MD0,NCP0,NDP0,NIP0
RETURN
C FORMAT STATEMENT FOR ERROR MESSAGE
2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S)./
1 7H MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6,
2 10X,5HNIP =,I6/
3 35H ERROR DETECTED IN ROUTINE IDBVIP/)
END
SUBROUTINE IDCLDP(NDP,XD,YD,NCP,IPC) ID002720
C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST
C TO EACH OF THE DATA POINT.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C COORDINATES OF THE DATA POINTS,
C NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA
C POINTS.
C THE OUTPUT PARAMETER IS
C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE
C POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C EACH OF THE NDP DATA POINTS ARE TO BE STORED.
C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST
C NOT EXCEED 25.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),IPC(400)
DIMENSION DSQ0(25),IPC0(25)
DATA NCPMX/25/, LUN/6/
C STATEMENT FUNCTION
DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
C PRELIMINARY PROCESSING
10 NDP0=NDP
NCP0=NCP
IF(NDP0.LT.2) GO TO 90
IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0) GO TO 90
C CALCULATION
20 DO 59 IP1=1,NDP0
C - SELECTS NCP POINTS.
X1=XD(IP1)
Y1=YD(IP1)
J1=0
DSQMX=0.0
DO 22 IP2=1,NDP0
IF(IP2.EQ.IP1) GO TO 22
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))