Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/interp.for
There are 2 other files named interp.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C INTERP.FOR (FILENAME ON LIBRARY DECTAPE)
C INTERP, 2.14.1 (CALLING NAME, SUBLST NO.)
C INTERPOLATION AND CURVE FITTING
C ADAPTED BY B. GRANET FROM (A) CACM 15,10 (OCT. 1972) 914-918
C (B) CACM 17,1 (JAN. 1974) "BIVARIATE INTERPOLATION AND
C SMOOTH SURFACE FITTING BASED ON LOCAL PROCEDURES.
C REPRINTING PRIVILEGES WERE GRANTED BY PERMISSION OF THE
C ASSOCIATION FOR COMPUTING MACHINERY BUT NOT FOR PROFIT.
C LIBRARY DECTAPE PROGS. USED: USAGE.MAC
C FORWMU PROGS. USED: TTYPTY, ALLCOR, DEVCHR, EXISTS,
C EXIST, GES, GETPPN, JOBNUM, PRINTS, RUNUUO
C APLB10 PROGS. USED: IO, GETFOR
C INTERNAL SUBR. USED: INTRPL, SFCFIT, CRVFIT, ITPLBV
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
1,IRSP,IDUMMY,JDUMMY
COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
INTEGER ARRAY(8)
DIMENSION A(1)
DATA ARRAY/'UN1CF','UN2CF','BIVCF','BIVIN','UNINT',
1'HEAD','FORM','HELP'/
IPAGCT=-1
IDLG=-1
IRSP=-4
NDEVI=4
NDEVO=6
NCOL=5
FMT=0
WRITE(IDLG,1)
C CALL USAGE('INTERP')
C---------------1, 0, NDEVO, NDEVI ARE INPUT. OTHER ARGS. ARE
C--------------- RETURNED. 1 MEANS OUTPUT? PRINTS, 0 MEANS INPUT? PRINTS
CALL IO(1,NDEVO,DEVNAM,IDVO,IFLNMO,IPJ,IPG,IBNK)
33 CALL IO(0,NDEVI,DEVNAM,IDVI,IFLNMI,IPJ,IPG,IBNK)
C---------------ICODE RETURNED, =0 MEANS TERMINAL JOB, =1 MEANS BATCH
CALL TTYPTY(ICODE)
46 WRITE(IDLG,3)
11 WRITE(IDLG,4)
READ(IRSP,5)JTYPE
DO 6 II=1,8
IF(ARRAY(II).EQ.JTYPE)
1GO TO(7,7,8,8,9,97,98,112),II
6 CONTINUE
GO TO 22
97 WRITE(IDLG,99)
54 READ(IRSP,100)(IDENT(I),I=1,4)
GO TO 46
98 FMT=1.0
GO TO 46
7 WRITE(IDLG,102)
103 READ(IRSP,29,ERR=106)L,M
N=(L-1)*M+1
IF(II.EQ.1)MD=1
IF(II.EQ.2)MD=2
MAX=2*(L+N)
CALL ALLCOR(MAX,IERR,I1,A)
IF(IERR.EQ.0)GO TO 114
WRITE(IDLG,115)
GO TO 7
114 I2=I1+L
I3=I2+L
I4=I3+N
CALL CRVFIT(NDEVO,MD,L,M,N,A(I1),A(I2),A(I3),A(I4))
GO TO 33
8 IF(II.EQ.3)WRITE(IDLG,104)
IF(II.EQ.4)WRITE(IDLG,109)
60 READ(IRSP,29,ERR=105)LX,LY,MX,MY
IF(II.EQ.4)GO TO 107
NU=(LX-1)*MX+1
NV=(LY-1)*MY+1
MAX=LX+LY+LX*LY+NU+NV+NU*NV
CALL ALLCOR(MAX,IERR,I1,A)
IF(IERR.EQ.0)GO TO 116
WRITE(IDLG,115)
GO TO 8
116 I2=I1+LX
I3=I2+LY
I4=I3+LX*LY
I5=I4+NU
I6=I5+NV
CALL SFCFIT(NDEVO,LX,LY,MX,MY,NU,NV,A(I1),A(I2),A(I3),
1A(I4),A(I5),A(I6))
GO TO 33
107 WRITE(IDLG,108)
111 READ(IRSP,29,ERR=63)NINTRP
MAX=LX+LY+LX*LY+3*NINTRP
CALL ALLCOR(MAX,IERR,I1,A)
IF(IERR.EQ.0)GO TO 117
WRITE(IDLG,115)
GO TO 8
117 I2=I1+LX
I3=I2+LY
I4=I3+LX*LY
I5=I4+NINTRP
I6=I5+NINTRP
CALL ITPLBV(NDEVO,LX,LY,NINTRP,A(I1),A(I2),A(I3),A(I4),
1A(I5),A(I6))
GO TO 33
9 WRITE(IDLG,110)
55 READ(IRSP,29,ERR=71)L,NINTRP
MAX=2*(L+NINTRP)
CALL ALLCOR(MAX,IERR,I1,A)
IF(IERR.EQ.0)GO TO 118
WRITE(IDLG,115)
GO TO 9
118 I2=I1+L
I3=I2+L
I4=I3+NINTRP
CALL INTRPL(L,NINTRP,A(I1),A(I2),A(I3),A(I4))
GO TO 33
112 WRITE(IDLG,113)
113 FORMAT(1X,'THE AVAILABLE OPTIONS ARE: '/
11X,'UN1CF-UNIVARIATE,SINGLE VALUED FUNCTION,CURVE FITTING'/
21X,'UN2CF-UNIVARIATE MULTIVALUED FUNCTION,CURVE FITTING'/
31X,'BIVCF-BIVARIATE CURVE FITTING'/
41X,'BIVIN- BIVARIATE INTERPOLATION'/
51X,'UNINT-UNIVARIATE,SINGLE VALUED FUNCTION,INTERPOLATION'/
61X,'HEAD-ENTER IDENTIFICATION FOR OUTPUT.'/
71X,'FORM-ENTER USER SPECIFIED FORMAT.'/)
GO TO 46
105 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,53)
GO TO 60
106 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,53)
GO TO 103
63 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,53)
GO TO 111
71 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,53)
GO TO 55
22 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,47)JTYPE
GO TO 11
1 FORMAT(1X,'WMU INTERPOLATION'/)
3 FORMAT(1X,'ENTER OPTION.'/)
4 FORMAT(/' *',$)
5 FORMAT(A5)
29 FORMAT(4I)
47 FORMAT(1X,A5,3X,' IS NOT VALID. TRY AGAIN.'/)
53 FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
99 FORMAT(1X,'ENTER IDENTIFICATION.'/)
100 FORMAT(16A5)
102 FORMAT(1X,'ENTER NO. OF INPUT PTS. AND SUBINTERVALS.'/)
104 FORMAT(1X,'ENTER NO. PTS. IN X AND Y COORDINATES AND NO. '/
11X,'OF SUBINTERVALS IN X AND Y COORD. SEPARATED BY COMMAS.'/)
108 FORMAT(1X,'ENTER NO. OF PTS. OF INTERP.'/)
109 FORMAT(1X,'ENTER NO. OF PTS. IN X AND Y COORD.',
1' SEPARATED BY COMMA.'/)
110 FORMAT(1X,'ENTER NO. INPUT PTS. AND NO. PTS. OF INTERP.'/)
115 FORMAT(1X,'NOT ABLE TO ALLOCATE CORE'/)
END
C---------------L, N ARE INPUT. OTHER ARGS. ARE RETURNED.
SUBROUTINE INTRPL(L,N,X,Y,U,V)
C INTERPOLATION OF A SINGLED-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
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB. IDENT,
C--------------- FMT, NDEVI, NDEVO, ICODE, IDVI ARE INPUT THRU COMMON
C--------------- /ALL/. IFMT IS RETURNED THRU COMMON /ALL/
COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
1,IRSP,IDUMMY,JDUMMY
COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
DIMENSION X(1),Y(1),U(1),V(1)
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
109 FORMAT(1X,'DATA BEING PROCESSED.'/)
270 FORMAT(1X,'ENTER PTS. OF INTERP.'/)
300 FORMAT(1X,' I',T10,'X(I)',T23,'Y(I)'/)
320 FORMAT(1X,I3,3X,F10.3,3X,F10.3)
490 FORMAT(1X,'ENTER I,XI,YI IN THIS ORDER.'/)
530 FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
810 FORMAT(F)
830 FORMAT(/)
1000 FORMAT(16A5)
J=0
K=0
ISTD=1
IF(FMT.NE.1)GO TO 101
C---------------IFMT, ISTD ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 96=NO. OF FMT. WORDS FOR OBJ. TIME FORMAT (6 LINES).
C--------------- 4 MEANS UNRESTRICTED FORMAT.
CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
101 IF(ISTD.EQ.1)IFMT(1)='(I,2F'
IF(ISTD.EQ.1)IFMT(2)=')'
WRITE(IDLG,270)
105 K=K+1
104 READ(IRSP,810,ERR=710)U(K)
IF(K.NE.N)GO TO 105
IF(IDVI.EQ.'TTY')GO TO 108
WRITE(IDLG,109)
GO TO 103
108 WRITE(IDLG,490)
100 J=J+1
103 READ(NDEVI,IFMT,ERR=520)I1,XI,YI
J=I1
X(J)=XI
Y(J)=YI
IF(J.NE.L)GO TO 100
GO TO 10
520 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,530)
GO TO 103
710 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,530)
GO TO 104
C---------------COME HERE FROM ST. 80+1 (NEAR END OF SUBR.)
110 WRITE(NDEVO,1000)(IDENT(I),I=1,4)
WRITE(NDEVO,300)
DO 310 I=1,N
310 WRITE(NDEVO,320)I,U(I),V(I)
WRITE(NDEVO,830)
RETURN
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 NUMERCIAL 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))
GO TO 110
C ERROR EXIT
90 WRITE (NDEVO,2090)
GO TO 99
91 WRITE (NDEVO,2091)
GO TO 99
95 WRITE (IU,2095)
GO TO 97
96 WRITE (NDEVO,2096)
97 WRITE (NDEVO,2097) I,X(I)
99 WRITE (NDEVO,2099) L0,N0
RETURN
C FORMAT SATATEMENTS
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,3H( =,I7/
1 36H ERROR DETECTED IN ROUTINE INTRPL)
END
C---------------X,Y,U,V ARE RETURNED. OTHER ARGS. ARE INPUT.
SUBROUTINE CRVFIT(IU,MD,L,M,N,X,Y,U,V)
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 (FOR ERROR
C MESSAGES
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
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/. FMT, NDEVI,
C--------------- NDEVO, IDENT, ICODE, IDVI ARE INPUT THRU COMMON
C--------------- /ALL/. IFMT IS RETURNED THRU COMMON /ALL/.
COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
1,IRSP,IDUMMY,JDUMMY
COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
DIMENSION X(1),Y(1),U(1),V(1)
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
109 FORMAT(1X,'DATA BEING PROCESSED.'/)
300 FORMAT(1X,' I',T10,'X(I)',T23,'Y(I)'/)
320 FORMAT(1X,I3,3X,F10.3,3X,F10.3)
490 FORMAT(1X,'ENTER I,XI,YI IN THIS ORDER.'/)
530 FORMAT(1X,'ERROR IN INPUT ,TRY AGAIN.'/)
830 FORMAT(/)
1000 FORMAT(16A5)
I=0
K=0
ISTD=1
IF(FMT.NE.1.0)GO TO 105
CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
105 IF(ISTD.EQ.1)IFMT(1)='(I,2F'
IF(ISTD.EQ.1)IFMT(2)=')'
IF(IDVI.EQ.'TTY')GO TO 108
WRITE(IDLG,109)
GO TO 100
108 WRITE(IDLG,490)
100 I=I+1
103 READ(NDEVI,IFMT,ERR=520)I1,XI,YI
I=I1
X(I)=XI
Y(I)=YI
IF(I.NE.L)GO TO 100
GO TO 10
520 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,530)
GO TO 103
110 WRITE(NDEVO,1000)(IDENT(I),I=1,4)
WRITE(NDEVO,300)
DO 310 I=1,N
310 WRITE(NDEVO,320)I,U(I),V(I)
WRITE(NDEVO,830)
RETURN
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
GO TO 110
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---------------IU,LX, LY ARE INPUT. OTHER ARGS. ARE RETURNED.
SUBROUTINE ITPLBV(IU,LX,LY,N,X,Y,Z,U,V,W)
C BIVARIATE INTERPOLATION
C THIS SUBROUTINE INTERPOLATES, FROM VALUES OF THE FUNCTION
C GIVEN AT INPUT GRID POINTS IN AN X-Y PLANE AND FOR A GIVEN
C SET OF POINTS IN THE PLANE, THE VALUES OF A SINGLE-VALUED
C BIVARIATE FUNCTION Z = Z(X,Y).
C THE METHOD IS BASED ON A PIECE-WISE FUNCTION COMPOSED OF
C A SET OF BICUBIC POLYNOMIALS IN X AND Y. EACH POLYNOMIAL
C IS APPLICABLE TO A RECTANGLE OF THE INPUT GRID IN THE X-Y
C PLANE. EACH POLYNOMIAL IS DETERMINED LOCALLY.
C THE INPUT PARAMETERS ARE
C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT
C LX = NUMBER OF INPUT GRID POINTS IN THE X COORDINATE
C (MUST BE 2 OR GREATER)
C LY = NUMBER OF INPUT GRID POINTS IN THE Y COORDINATE
C (MUST BE 2 OR GREATER)
C X = ARRAY OF DIMENSION LX STORING THE X COORDINATES
C OF INPUT GRID POINTS (IN ASCENDING ORDER)
C Y = ARRAY OF DIMENSION LY STORING THE Y COORDINATES
C OF INPUT GRID POINTS (IN ASENDING ORDER)
C Z = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (LX,LY)
C STORING THE VALUES OF THE FUNCTION (Z VALUES)
C AT INPUT GRID POINTS
C N = NUMBER OF POINTS AT WHICH INTERPOLATION OF THE
C Z VALUE IS DESIRED (MUST BE 1 OR GREATER)
C U = ARRAY OF DIMENSION N STORING THE X COORDINATES
C OF DESIRED POINTS
C V = ARRAY OF DIMENSION N STORING THE Y COORDINATES
C OF DESIRED POINTS
C THE OUTPUT PARAMETER IS
C W = ARRAY OF DIMENSION N WHERE THE INTERPOLATED Z
C VALUES AT DESIRED POINTS ARE TO BE DISPLAYED
C SOME VARIABLES INTERNALLY USED ARE
C ZA = DIVIDED DIFFERENCE OF Z WITH RESPECT TO X
C ZB = DIVIDED DIFFERENCE OF Z WITH RESPECT TO Y
C ZAB = SECOND ORDER DIVIDED DIFFERENCE OF Z WITH
C RESPECT TO X AND Y
C ZX = PARTIAL DERIVATIVE OF Z WITH RESPECT TO X
C ZY = PARTIAL DERIVATIVE OF Z WITH RESPECT TO Y
C ZXY = SECOND ORDER PARTIAL DERIVATIVE OF Z WITH
C RESPECT TO X AND Y
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/. IDENT,
C--------------- FMT, NDEVI, NDEVO, ICODE, IDVI ARE INPUT THRU COMMON
C--------------- /ALL/. IFMT IS RETURNED THRU COMMON /ALL/.
COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
1,IRSP,IDUMMY,JDUMMY
COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
DIMENSION X(1), Y(1), Z(1), U(1), V(1), W(1)
DIMENSION ZA(5,2), ZB(2,5), ZAB(3,3), ZX(4,4), ZY(4,4),
* ZXY(4,4)
EQUIVALENCE (Z3A1,ZA(1)), (Z3A2,ZA(2)), (Z3A3,ZA(3)),
* (Z3A4,ZA(4)), (Z3A5,ZA(5)), (Z4A1,ZA(6)), (Z4A2,ZA(7)),
* (Z4A3,ZA(8)), (Z4A4,ZA(9)), (Z4A5,ZA(10)), (Z3B1,ZB(1)),
* (Z3B2,ZB(3)), (Z3B3,ZB(5)), (Z3B4,ZB(7)), (Z3B5,ZB(9)),
* (Z4B1,ZB(2)), (Z4B2,ZB(4)), (Z4B3,ZB(6)), (Z4B4,ZB(8)),
* (Z4B5,ZB(10)), (ZA2B2,ZAB(1)), (ZA3B2,ZAB(2)),
* (ZA4B2,ZAB(3)), (ZA2B3,ZAB(4)), (ZA3B3,ZAB(5)),
* (ZA4B3,ZAB(6)), (ZA2B4,ZAB(7)), (ZA3B4,ZAB(8)),
* (ZA4B4,ZAB(9)), (ZX33,ZX(6)), (ZX43,ZX(7)),
* (ZX34,ZX(10)), (ZX44,ZX(11)), (ZY33,ZY(6)),
* (ZY43,ZY(7)), (ZY34,ZY(10)), (ZY44,ZY(11)),
* (ZXY33,ZXY(6)), (ZXY43,ZXY(7)), (ZXY34,ZXY(10)),
* (ZXY44,ZXY(11)), (P00,Z33), (P01,ZY33), (P10,ZX33),
* (P11,ZXY33)
EQUIVALENCE (LX0,ZX(1)), (LXM1,ZX(4)), (LXM2,ZX(13)),
* (LXP1,ZX(16)), (LY0,ZY(1)), (LYM1,ZY(4)), (LYM2,ZY(13)),
* (LYP1,ZY(16)), (IX,ZXY(1)), (IY,ZXY(4)), (IXPV,ZXY(13)),
* (IYPV,ZXY(16)), (IMN,JX), (IMX,JY), (JXM2,JX1),
* (JYM2,JY1), (UK,DX), (VK,DY), (A1,A5,B1,B5,ZX(2),A,Q0),
* (A2,ZX(5),B,Q1), (A4,ZX(8),C,Q2), (B2,ZY(2),D,Q3),
* (B4,ZY(14),E), (X2,ZX(3),A3SQ), (X4,ZX(9)), (X5,ZX(12)),
* (Y2,ZX(14)), (Y4,ZY(3),B3SQ), (Y5,ZX(15),P02),
* (Z23,ZY(5),P03), (Z24,ZY(8),P12), (Z32,ZY(9),P13),
* (Z34,ZY(12),P20), (Z35,ZY(15),P21), (Z42,ZXY(2),P22),
* (Z43,ZXY(5),P23), (Z44,ZXY(3),P30), (Z45,ZXY(8),P31),
* (Z53,ZXY(9),P32), (Z54,ZXY(12),P33), (W2,WY2,W4),
* (W3,WY3,W1,W5), (WX2,ZXY(14)), (WX3,ZXY(15))
C PRELIMINARY PROCESSING
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES
917 FORMAT(1X,'ENTER I,J,XI,YJ,ZIJ SEPARATED BY COMMAS.'/)
916 FORMAT(/)
900 FORMAT(1X,'ENTER PTS. OF INTERP. X,Y.'/)
901 FORMAT(2F)
904 FORMAT(1X,'ERROR IN INPUT ,TRY AGAIN.'/)
906 FORMAT(1X,'DATA BEING PROCESSED.'/)
912 FORMAT(16A5)
913 FORMAT(1X,' I',T10,'X(I)',T23,'Y(J)',T36,'Z(I,J)'/)
915 FORMAT(1X,I3,3X,F10.3,3X,F10.3,3X,F10.3)
DO 918 I=1,4
DO 918 J=1,4
ZX(I,J)=0
ZY(I,J)=0
ZXY(I,J)=0
918 CONTINUE
DO 919 I=1,5
DO 919 J=1,2
919 ZA(I,J)=0
DO 920 I=1,2
DO 920 J=1,5
920 ZB(I,J)=0
DO 921 I=1,3
DO 921 J=1,3
921 ZAB(I,J)=0
K=0
NOPTS=LX*LY
ISTD=1
IF(FMT.NE.1)GO TO 905
CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
905 IF(ISTD.EQ.1)IFMT(1)='(2I,3'
IF(ISTD.EQ.1)IFMT(2)='F)'
WRITE(IDLG,900)
903 K=K+1
READ(IRSP,901,ERR=902)U(K),V(K)
IF(K.NE.N)GO TO 903
K=0
IF(IDVI.EQ.'TTY')GO TO 907
WRITE(IDLG,906)
GO TO 908
907 WRITE(IDLG,917)
908 K=K+1
909 READ(NDEVI,IFMT,ERR=910)I1,J1,XI,YJ,ZIJ
I=I1
J=J1
X(I)=XI
Y(J)=YJ
IJ=(I-1)*LY+J
Z(IJ)=ZIJ
IF(K.NE.NOPTS)GO TO 908
GO TO 1
902 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,904)
GO TO 905
910 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,904)
GO TO 909
911 WRITE(NDEVO,912)(IDENT(I),I=1,4)
WRITE(NDEVO,913)
DO 914 I=1,N
914 WRITE(NDEVO,915)I,U(I),V(I),W(I)
WRITE(NDEVO,916)
RETURN
1 IU0 = IU
LX0 = LX
LXM1 = LX0 - 1
LXM2 = LXM1 - 1
LXP1 = LX0 + 1
LY0 = LY
LYM1 = LY0 - 1
LYM2 = LYM1 - 1
LYP1 = LY0 + 1
N0 = N
C ERROR CHECK
IF (LXM2.LT.0) GO TO 710
IF (LYM2.LT.0) GO TO 720
IF (N0.LT.1) GO TO 730
DO 10 IX=2,LX0
IF (X(IX-1)-X(IX)) 10, 740, 750
10 CONTINUE
DO 20 IY=2,LY0
IF (Y(IY-1)-Y(IY)) 20, 770, 780
20 CONTINUE
C INITIAL SETTING OF PREVIOUS VALUES OF IX AND IY
IXPV = 0
IYPV = 0
C MAIN DO-LOOP
DO 700 K=1,N0
UK = U(K)
VK = V(K)
C ROUTINE TO LOCATE THE DESIRED POINT
C TO FIND OUT THE IX VALUE FOR WHICH
C (U(K).GE.X(IX-1).AND.(U(K).LT.X(IX))
IF (LXM2.EQ.0) GO TO 80
IF (UK.GE.X(LX0)) GO TO 70
IF (UK.LT.X(1)) GO TO 60
IMN = 2
IMX = LX0
30 IX = (IMN+IMX)/2
IF (UK.GE.X(IX)) GO TO 40
IMX = IX
GO TO 50
40 IMN = IX + 1
50 IF (IMX.GT.IMN) GO TO 30
IX = IMX
GO TO 90
60 IX = 1
GO TO 90
70 IX = LXP1
GO TO 90
80 IX = 2
C TO FIND OUT THE IY VALUE FOR WHICH
C (V(K).GE.Y(IY-1)).AND.(V(K).LT.Y(IY))
90 IF (LYM2.EQ.0) GO TO 150
IF (VK.GE.Y(LY0)) GO TO 140
IF (VK.LT.Y(1)) GO TO 130
IMN = 2
IMX = LY0
100 IY = (IMN+IMX)/2
IF (VK.GE.Y(IY)) GO TO 110
IMX = IY
GO TO 120
110 IMN = IY + 1
120 IF (IMX.GT.IMN) GO TO 100
IY = IMX
GO TO 160
130 IY = 1
GO TO 160
140 IY = LYP1
GO TO 160
150 IY = 2
C TO CHECK IF THE DESIRED POINT IS IN THE SAME RECTANGLE
C AS THE PREVIOUS POINT. IF YES, SKIP TO THE COMPUTATION
C OF THE POLYNOMIAL
160 IF (IX.EQ.IXPV .AND. IY.EQ.IYPV) GO TO 690
IXPV = IX
IYPV = IY
C ROUTINES TO PICK UP NECESSARY X, Y, AND Z VALUES, TO
C COMPUTE THE ZA, ZB, AND ZAB VALUES, AND TO ESTIMATE THEM
C WHEN NECESSARY
JX = IX
IF (JX.EQ.1) JX = 2
IF (JX.EQ.LXP1) JX = LX0
JY = IY
IF (JY.EQ.1) JY = 2
IF (JY.EQ.LYP1) JY = LY0
JXM2 = JX - 2
JXML = JX - LX0
JYM2 = JY - 2
JYML = JY - LY0
C IN THE CORE AREA, I.E., IN THE RECTANGLE THAT CONTAINS
C THE DESIRED POINT
X3 = X(JX-1)
X4 = X(JX)
A3 = 1.0/(X4-X3)
Y3 = Y(JY-1)
Y4 = Y(JY)
B3 = 1.0/(Y4-Y3)
JJ1=(JX-2)*LY+JY-1
JJ2=JJ1+LY
JJ3=JJ1+1
JJ4=JJ3+LY
Z33=Z(JJ1)
Z43=Z(JJ2)
Z34=Z(JJ3)
Z44=Z(JJ4)
Z3A3 = (Z43-Z33)*A3
Z4A3 = (Z44-Z34)*A3
Z3B3 = (Z34-Z33)*B3
Z4B3 = (Z44-Z43)*B3
ZA3B3 = (Z4B3-Z3B3)*A3
C IN THE X DIRECTION
IF (LXM2.EQ.0) GO TO 230
IF (JXM2.EQ.0) GO TO 170
X2 = X(JX-2)
A2 = 1.0/(X3-X2)
JZ23=(JX-3)*LY+JY-1
JZ24=JZ23+1
Z23=Z(JZ23)
Z24=Z(JZ24)
Z3A2 = (Z33-Z23)*A2
Z4A2 = (Z34-Z24)*A2
IF (JXML.EQ.0) GO TO 180
170 X5 = X(JX+1)
A4 = 1.0/(X5-X4)
JZ53=JX*LY+JY-1
JZ54=JZ53+1
Z53=Z(JZ53)
Z54=Z(JZ54)
Z3A4 = (Z53-Z43)*A4
Z4A4 = (Z54-Z44)*A4
IF (JXM2.NE.0) GO TO 190
Z3A2 = Z3A3 + Z3A3 - Z3A4
Z4A2 = Z4A3 + Z4A3 - Z4A4
GO TO 190
180 Z3A4 = Z3A3 + Z3A3 - Z3A2
Z4A4 = Z4A3 + Z4A3 - Z4A2
190 ZA2B3 = (Z4A2-Z3A2)*B3
ZA4B3 = (Z4A4-Z3A4)*B3
IF (JX.LE.3) GO TO 200
A1 = 1.0/(X2-X(JX-3))
J3A1=(JX-4)*LY+JY-1
J4A1=J3A1+1
Z3A1=(Z23-Z(J3A1))*A1
Z4A1=(Z24-Z(J4A1))*A1
GO TO 210
200 Z3A1 = Z3A2 + Z3A2 - Z3A3
Z4A1 = Z4A2 + Z4A2 - Z4A3
210 IF (JX.GE.LXM1) GO TO 220
A5 = 1.0/(X(JX+2)-X5)
J3A5=(JX+1)*LY+JY-1
J4A5=J3A5+1
Z3A5=(Z(J3A5)-Z53)*A5
Z4A5=(Z(J4A5)-Z54)*A5
GO TO 240
220 Z3A5 = Z3A4 + Z3A4 - Z3A3
Z4A5 = Z4A4 + Z4A4 - Z4A3
GO TO 240
230 Z3A2 = Z3A3
Z4A2 = Z4A3
GO TO 180
C IN THE Y DIRECTION
240 IF (LYM2.EQ.0) GO TO 310
IF (JYM2.EQ.0) GO TO 250
Y2 = Y(JY-2)
B2 = 1.0/(Y3-Y2)
JZ32=(JX-2)*LY+JY-2
JZ42=(JX-1)*LY+JY-2
Z32=Z(JZ32)
Z42=Z(JZ42)
Z3B2 = (Z33-Z32)*B2
Z4B2 = (Z43-Z42)*B2
IF (JYML.EQ.0) GO TO 260
250 Y5 = Y(JY+1)
B4 = 1.0/(Y5-Y4)
JZ35=(JX-2)*LY+JY+1
JZ45=(JX-1)*LY+JY+1
Z35=Z(JZ35)
Z45=Z(JZ45)
Z3B4=(Z35-Z34)*B4
Z4B4 = (Z45-Z44)*B4
IF (JYM2.NE.0) GO TO 270
Z3B2 = Z3B3 + Z3B3 - Z3B4
Z4B2 = Z4B3 + Z4B3 - Z4B4
GO TO 270
260 Z3B4 = Z3B3 + Z3B3 - Z3B2
Z4B4 = Z4B3 + Z4B3 - Z4B2
270 ZA3B2 = (Z4B2-Z3B2)*A3
ZA3B4 = (Z4B4-Z3B4)*A3
IF (JY.LE.3) GO TO 280
B1 = 1.0/(Y2-Y(JY-3))
J3B1=(JX-2)*LY+JY-3
J4B1=(JX-1)*LY+JY-3
Z3B1=(Z32-Z(J3B1))*B1
Z4B1=(Z42-Z(J4B1))*B1
GO TO 290
280 Z3B1 = Z3B2 + Z3B2 - Z3B3
Z4B1 = Z4B2 + Z4B2 - Z4B3
290 IF (JY.GE.LYM1) GO TO 300
B5 = 1.0/(Y(JY+2)-Y5)
J3B5=(JX-2)*LY+JY+2
J4B5=(JX-1)*LY+JY+2
Z3B5=(Z(J3B5)-Z35)*B5
Z4B5=(Z(J4B5)-Z45)*B5
GO TO 320
300 Z3B5 = Z3B4 + Z3B4 - Z3B3
Z4B5 = Z4B4 + Z4B4 - Z4B3
GO TO 320
310 Z3B2 = Z3B3
Z4B2 = Z4B3
GO TO 260
C IN THE DIAGONAL DIRECTIONS
320 IF (LXM2.EQ.0) GO TO 400
IF (LYM2.EQ.0) GO TO 410
IF (JXML.EQ.0) GO TO 350
IF (JYM2.EQ.0) GO TO 330
J4B2=JZ53-1
ZA4B2=((Z53-Z(J4B2))*B2-Z4B2)*A4
IF (JYML.EQ.0) GO TO 340
330 J4B4=J4B2+3
ZA4B4=((Z(J4B4)-Z54)*B4-Z4B4)*A4
IF (JYM2.NE.0) GO TO 380
ZA4B2 = ZA4B3 + ZA4B3 - ZA4B4
GO TO 380
340 ZA4B4 = ZA4B3 + ZA4B3 - ZA4B2
GO TO 380
350 IF (JYM2.EQ.0) GO TO 360
J2B2=(JX-3)*LY+JY-2
ZA2B2=(Z3B2-(Z23-Z(J2B2))*B2)*A2
IF (JYML.EQ.0) GO TO 370
360 J2B4=J2B2+3
ZA2B4=(Z3B4-(Z(J2B4)-Z24)*B4)*A2
IF (JYM2.NE.0) GO TO 390
ZA2B2 = ZA2B3 + ZA2B3 - ZA2B4
GO TO 390
370 ZA2B4 = ZA2B3 + ZA2B3 - ZA2B2
GO TO 390
380 IF (JXM2.NE.0) GO TO 350
ZA2B2 = ZA3B2 + ZA3B2 - ZA4B2
ZA2B4=ZA3B4+ZA3B4-ZA4B4
GO TO 420
390 IF (JXML.NE.0) GO TO 420
ZA4B2 = ZA3B2 + ZA3B2 - ZA2B2
ZA4B4 = ZA3B4 + ZA3B4 - ZA2B4
GO TO 420
400 ZA2B2 = ZA3B2
ZA4B2 = ZA3B2
ZA2B4 = ZA3B4
ZA4B4 = ZA3B4
GO TO 420
410 ZA2B2 = ZA2B3
ZA2B4 = ZA2B3
ZA4B2 = ZA4B3
ZA4B4 = ZA4B3
C NUMERICAL DIFFERENTIATION --- TO DETERMINE PARTIAL
C DERIVATIVES ZX, ZY, AND ZXY AS WEIGHTED MEANS OF DIVIDED
C DIFFERENCES ZA, ZB, AND ZAB, RESPECTIVELY
420 DO 480 JY=2,3
DO 470 JX=2,3
W2 = ABS(ZA(JX+2,JY-1)-ZA(JX+1,JY-1))
W3 = ABS(ZA(JX,JY-1)-ZA(JX-1,JY-1))
SW = W2 + W3
IF (SW.EQ.0.0) GO TO 430
WX2 = W2/SW
WX3 = W3/SW
GO TO 440
430 WX2 = 0.5
WX3 = 0.5
440 ZX(JX,JY) = WX2*ZA(JX,JY-1) + WX3*ZA(JX+1,JY-1)
W2 = ABS(ZB(JX-1,JY+2)-ZB(JX-1,JY+1))
W3 = ABS(ZB(JX-1,JY)-ZB(JX-1,JY-1))
SW = W2 + W3
IF (SW.EQ.0.0) GO TO 450
WY2 = W2/SW
WY3 = W3/SW
GO TO 460
450 WY2 = 0.5
WY3 = 0.5
460 ZY(JX,JY) = WY2*ZB(JX-1,JY) + WY3*ZB(JX-1,JY+1)
ZXY(JX,JY) =
* WY2*(WX2*ZAB(JX-1,JY-1)+WX3*ZAB(JX,JY-1)) +
* WY3*(WX2*ZAB(JX-1,JY)+WX3*ZAB(JX,JY))
470 CONTINUE
480 CONTINUE
IF (IX.EQ.LXP1) GO TO 530
IF (IX.NE.1) GO TO 590
W2 = A4*(3.0*A3+A4)
W1 = 2.0*A3*(A3-A4)+W2
DO 500 JY = 2,3
ZX(1,JY)=(W1*ZA(1,JY-1)+W2*ZA(2,JY-1))/(W1+W2)
ZY(1,JY) = ZY(2,JY) + ZY(2,JY) - ZY(3,JY)
ZXY(1,JY) = ZXY(2,JY) + ZXY(2,JY) - ZXY(3,JY)
DO 490 JX1=2,3
JX = 5 - JX1
ZX(JX,JY) = ZX(JX-1,JY)
ZY(JX,JY) = ZY(JX-1,JY)
ZXY(JX,JY) = ZXY(JX-1,JY)
490 CONTINUE
500 CONTINUE
X3 = X3 - 1.0/A4
Z33 = Z33 - Z3A2/A4
DO 510 JY=1,5
ZB(2,JY) = ZB(1,JY)
510 CONTINUE
DO 520 JY=2,4
ZB(1,JY) = ZB(1,JY) - ZAB(1,JY-1)/A4
520 CONTINUE
A3 = A4
JX = 1
GO TO 570
530 W4 = A2*(3.0*A3+A2)
W5 = 2.0*A3*(A3-A2) + W4
DO 550 JY=2,3
ZX(4,JY) = (W4*ZA(4,JY-1)+W5*ZA(5,JY-1))/(W4+W5)
ZY(4,JY) = ZY(3,JY) + ZY(3,JY) - ZY(2,JY)
ZXY(4,JY) = ZXY(3,JY) + ZXY(3,JY) - ZXY(2,JY)
DO 540 JX=2,3
ZX(JX,JY) = ZX(JX+1,JY)
ZY(JX,JY) = ZY(JX+1,JY)
ZXY(JX,JY) = ZXY(JX+1,JY)
540 CONTINUE
550 CONTINUE
X3 = X4
Z33 = Z43
DO 560 JY=1,5
ZB(1,JY) = ZB(2,JY)
560 CONTINUE
A3 = A2
JX = 3
570 ZA(3,1) = ZA(JX+1,1)
DO 580 JY=1,3
ZAB(2,JY) = ZAB(JX,JY)
580 CONTINUE
C WHEN (V(K).LT.Y(1)).OR.(V(K).GT.Y(LY))
590 IF (IY.EQ.LYP1) GO TO 630
IF (IY.NE.1) GO TO 680
W2 = B4*(3.0*B3+B4)
W1 = 2.0*B3*(B3-B4) + W2
DO 620 JX=2,3
IF (JX.EQ.3 .AND. IX.EQ.LXP1) GO TO 600
IF (JX.EQ.2 .AND. IX.EQ.1) GO TO 600
ZY(JX,1) = (W1*ZB(JX-1,1)+W2*ZB(JX-1,2))/(W1+W2)
ZX(JX,1) = ZX(JX,2) + ZX(JX,2) - ZX(JX,3)
ZXY(JX,1) = ZXY(JX,2) + ZXY(JX,2) - ZXY(JX,3)
600 DO 610 JY1=2,3
JY = 5 - JY1
ZY(JX,JY) = ZY(JX,JY-1)
ZX(JX,JY) = ZX(JX,JY-1)
ZXY(JX,JY) = ZXY(JX,JY-1)
610 CONTINUE
620 CONTINUE
Y3 = Y3 - 1.0/B4
Z33 = Z33 - Z3B2/B4
Z3A3 = Z3A3 - ZA3B2/B4
Z3B3 = Z3B2
ZA3B3 = ZA3B2
B3 = B4
GO TO 670
630 W4 = B2*(3.0*B3+B2)
W5 = 2.0*B3*(B3-B2) + W4
DO 660 JX=2,3
IF (JX.EQ.3 .AND. IX.EQ.LXP1) GO TO 640
IF (JX.EQ.2 .AND. IX.EQ.1) GO TO 640
ZY(JX,4) = (W4*ZB(JX-1,4)+W5*ZB(JX-1,5))/(W4+W5)
ZX(JX,4) = ZX(JX,3) + ZX(JX,3) - ZX(JX,2)
ZXY(JX,4) = ZXY(JX,3) + ZXY(JX,3) - ZXY(JX,2)
640 DO 650 JY=2,3
ZY(JX,JY) = ZY(JX,JY+1)
ZX(JX,JY) = ZX(JX,JY+1)
ZXY(JX,JY) = ZXY(JX,JY+1)
650 CONTINUE
660 CONTINUE
Y3 = Y4
Z33 = Z33 + Z3B3/B3
Z3A3 = Z3A3 + ZA3B3/B3
Z3B3 = Z3B4
ZA3B3 = ZA3B4
B3 = B2
670 IF (IX.NE.1 .AND. IX.NE.LXP1) GO TO 680
JX = IX/LXP1 + 2
JX1 = 5 - JX
JY = IY/LYP1 + 2
JY1 = 5 - JY
ZX(JX,JY) = ZX(JX1,JY) + ZX(JX,JY1) - ZX(JX1,JY1)
ZY(JX,JY) = ZY(JX1,JY) + ZY(JX,JY1) - ZY(JX1,JY1)
ZXY(JX,JY) = ZXY(JX1,JY) + ZXY(JX,JY1) - ZXY(JX1,JY1)
C DETERMINATION OF THE COEFFICIENTS OF THE POLYNOMIAL
680 ZX3B3 = (ZX34-ZX33)*B3
ZX4B3 = (ZX44-ZX43)*B3
ZY3A3 = (ZY43-ZY33)*A3
ZY4A3 = (ZY44-ZY34)*A3
A = ZA3B3 - ZX3B3 - ZY3A3 + ZXY33
B = ZX4B3 - ZX3B3 - ZXY43 + ZXY33
C = ZY4A3 - ZY3A3 - ZXY34 + ZXY33
D = ZXY44 - ZXY43 - ZXY34 + ZXY33
E = A + A - B - C
A3SQ = A3*A3
B3SQ = B3*B3
P02=(2.0*(Z3B3-ZY33)+Z3B3-ZY34)*B3
P03=(-2.0*Z3B3+ZY34+ZY33)*B3SQ
P12=(2.0*(ZX3B3-ZXY33)+ZX3B3-ZXY34)*B3
P13=(-2.0*ZX3B3+ZXY34+ZXY33)*B3SQ
P20=(2.0*(Z3A3-ZX33)+Z3A3-ZX43)*A3
P21=(2.0*(ZY3A3-ZXY33)+ZY3A3-ZXY43)*A3
P22=(3.0*(A+E)+D)*A3*B3
P23=(-3.0*E-B-D)*A3*B3SQ
P30 = (-2.0*Z3A3+ZX43+ZX33)*A3SQ
P31 = (-2.0*ZY3A3+ZXY43+ZXY33)*A3SQ
P32 = (-3.0*E-C-D)*B3*A3SQ
P33 = (D+E+E)*A3SQ*B3SQ
C COMPUTATION OF THE POLYNOMIAL
690 DY = VK - Y3
Q0 = P00 + DY*(P01+DY*(P02+DY*P03))
Q1 = P10 + DY*(P11+DY*(P12+DY*P13))
Q2 = P20 + DY*(P21+DY*(P22+DY*P23))
Q3 = P30 + DY*(P31+DY*(P32+DY*P33))
DX = UK - X3
W(K) = Q0 + DX*(Q1+DX*(Q2+DX*Q3))
700 CONTINUE
C NORMAL EXIT
GO TO 911
C ERROR EXIT
710 WRITE (IU0,99999)
GO TO 800
720 WRITE (IU0,99998)
GO TO 800
730 WRITE (IU0,99997)
GO TO 800
740 WRITE (IU0,99996)
GO TO 760
750 WRITE (IU0,99995)
760 WRITE (IU0,99994) IX, X(IX)
GO TO 800
770 WRITE (IU0,99993)
GO TO 790
780 WRITE (IU0,99992)
790 WRITE (IU0,99991) IY, Y(IY)
800 WRITE (IU0,99990) LX0, LY0, N0
RETURN
C FORMAT STATEMENTS
99999 FORMAT(1X/23H *** LX = 1 OR LESS./)
99998 FORMAT(1X/23H *** LY = 1 OR LESS./)
99997 FORMAT(1X/22H *** N = 0 OR LESS./)
99996 FORMAT(1X/27H *** IDENTICAL X VALUES./)
99995 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./)
99994 FORMAT(7H IX =, I6, 10X, 7HX(IX) =, E12.3)
99993 FORMAT(1X/27H *** IDENTICAL Y VALUES./)
99992 FORMAT(1X/33H *** Y VALUES OUT OF SEQUENCE./)
99991 FORMAT(7H IY =, I6, 10X, 7HY(IY) =, E12.3)
99990 FORMAT(7H LX =, I6, 10X, 4HLY =, I6, 10X, 3HN =, I7/
*36H ERROR DETECTED IN ROUTINE ITPLBV)
END
C---------------IU, LX, LY, MX, MY, NU, NV ARE INPUT.
C--------------- OTHER ARGS. ARE RETURNED.
SUBROUTINE SFCFIT(IU,LX,LY,MX,MY,NU,NV,X,Y,Z,U,V,W)
C SMOOTH SURFACE FITTING
C THIS SUBROUTINE FITS A SMOOTH SURFACE OF A SINGLE-VALUED
C BIVARIATE FUNTION Z = Z(X,Y) TO A SET OF INPUT DATA
C POINTS GIVEN AT INPUT GRID POINTS IN AN X-Y PLANE. IT
C GENERATES A SET OF OUTPUT GRID POINTS BY EQUALLY DIVIDING
C THE X AND Y COORDINATES IN EACH INTERVAL BETWEEN A PAIR
C OF INPUT GRID POINTS, INTERPOLATES THE Z VALUE FOR THE
C X AND Y VALUES OF EACH OUTPUT GRID POINT, AND GENERATES
C A SET OF OUTPUT POINTS CONSISTING OF INPUT DATA POINTS
C AND THE INTERPOLATED POINTS.
C THE METHOD IS BASED ON A PIECE-WISE FUNTION COMPOSED OF
C A SET OF BICUBIC POLYNOMIALS IN X AND Y. EACH POLYNOMIAL
C IS APPLICABLE TO A RECTANGLE OF THE INPUT GRID IN THE X-Y
C PLANE. EACH POLYNOMIAL IS DETERMINED LOCALLY.
C THE INPUT PARAMETERS ARE
C IU = LOGICAL UNIT NUMBER OF STANDARD OUTPUT UNIT
C LX = NUMBER OF INPUT GRID POINTS IN THE X COORDINATE
C (MUST BE 2 OR GREATER)
C LY = NUMBER OF INPUT GRID POINTS IN THE Y COORDINATE
C (MUST BE 2 OR GREATER)
C X = ARRAY OF DIMENSION LX STORING THE X COORDINATES
C OF INPUT GRID POINTS (IN ASCENDING OR DESCENDING
C ORDER)
C Y = ARRAY OF DIMENSION LY STORING THE Y COORDINATES
C OF INPUT GRID POINTS (IN ASCENDING OR DESCENDING
C ORDER)
C Z = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (LX,LY)
C STORING THE VALUES OF THE FUNCTION AT INPUT
C GRID POINTS
C MX = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF
C INPUT GRID POINTS IN THE X COORDINATE
C (MUST BE 2 OR GREATER)
C MY = NUMBER OF SUBINTERVALS BETWEEN EACH PAIR OF
C INPUT GRID POINTS IN THE Y COORDINATE
C (MUST BE 2 OR GREATER)
C NU NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE
C = (LX-1)*MX+1
C NV = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE
C = (LY-1)*MY+1
C THE OUTPUT PARAMETERS ARE
C U = ARRAY OF DIMENSION NU WHERE THE X COORDINATES OF
C OUTPUT POINTS ARE TO BE DISPLAYED
C V = ARRAY OF DIMENSION NV WHERE THE Y COORINATES OF
C OUTPUT POINTS ARE TO BE DISPLAYED
C WHERE THE Z COORIDANETS OF OUTPUT POINTS ARE TO
C BE DISPLAYED
C SOME VARIABLES INTERNALLY USED ARE
C ZA = DIVIDED DIFFERENCE OF Z WITH RESPECT TO X
C ZB = DIVIDED DIFFERENCE OF Z WITH RESPECT TO Y
C ZAB = SECOND ORDER DIVIDED DIFFERENCE OF Z WITH
C RESPECT TO X AND Y
C ZX = PARTIAL DERIVATIVE OF Z WITH RESPECT TO X
C ZY = PARTIAL DERIVATIVE OF Z WITH RESPECT TO Y
C ZXY = SECOND ORDER PARTIAL DERIVATIVE OF Z WITH
C RESPECT TO X AND Y
C DECLARATION STATEMENTS
C---------------IDLG, IRSP ARE INPUT THRU COMMON /IOB/. IDENT, FMT, NDEVI,
C--------------- NDEVO, ICODE, IDVI ARE INPUT THRU COMMON /ALL/. IFMT
C--------------- IS RETURNED THRU COMMON /ALL/.
COMMON/IOB/LEFBK,IRTBK,IART,MAXPAG,IPAGE,IPAGCT,IDLG
1,IRSP,IDUMMY,JDUMMY
COMMON/ALL/II,IDENT(4),FMT,NDEVI,NDEVO,ICODE,IFMT(96),IDVI
DIMENSION X(1),Y(1),Z(1),U(1),V(1),W(1)
DIMENSION ZA(4,2), ZB(5), ZAB(2,3), ZX(2), ZY(2), ZXY(2)
EQUIVALENCE (Z3A2,ZA(1)), (Z3A3,ZA(2)), (Z3A4,ZA(3)),
* (Z3A5,ZA(4)), (Z4A2,ZA(5)), (Z4A3,ZA(6)), (Z4A4,ZA(7)),
* (Z4A5,ZA(8)), (Z4B1,ZB(1)), (Z4B2,ZB(2)), (Z4B3,ZB(3)),
* (Z4B4,ZB(4)), (Z4B5,ZB(5)), (ZA3B2,ZAB(1)),
* (ZA4B2,ZAB(2)), (ZA3B3,ZAB(3)), (ZA4B3,ZAB(4)),
* (ZA3B4,ZAB(5)), (ZA4B4,ZAB(6)), (ZX43,ZX(1)),
* (ZX44,ZX(2)), (ZY43,ZY(1)), (ZY44,ZY(2)),
* (ZXY43,ZXY(1)), (ZXY44,ZXY(2)), (P00,Z33), (P01,ZY33),
* (P10,ZX33), (P11,ZXY33)
EQUIVALENCE (IXM1,JX), (IXML,JY), (DU,DV,DX,DY),
* (FMX,RMX,FMY,RMY,SW,E), (W2,WY2,A,Q0), (W3,WY3,B,Q1),
* (WX2,C,Q2), (WX3,D,Q3), (Z3A2,P02), (Z4A2,P03),
* (Z4B1,P12), (Z4B2,P13), (Z4B4,P20), (Z4B5,P21),
* (ZA3B2,P22), (ZA3B4,P23)
C PRELIMINARY PROCESSING
C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES
919 FORMAT(1X,T14,1HI,5(F6.1,4X))
920 FORMAT(1X,60(1H-)/)
922 FORMAT(1X,I3,4X,F5.1,1HI,5(F7.2,3X))
924 FORMAT(//1X,T14,1HI,48(1H-)/1X,T14,1HI,T15,' J ='
1,/1X,T14,1HI,T15,I3,T25,I3,T35,I3,T45,I3,
2T55,I3)
926 FORMAT(1X,55(1H-)/1X,'I',5X,'X(I)',T14,1HI,T15,'Y(J) =')
906 FORMAT(1X,'DATA BEING PROCESSED.'/)
916 FORMAT(/)
912 FORMAT(16A5)
913 FORMAT(1X,T14,1HI,T32,'W(I,J)'/1X,T14,1HI,48(1H-)/1X,
1T14,1HI,T15,' J= ',/1X,T14,1HI,T15,'1',T25,'2',T35,'3',
2T45,'4',T55,'5',/1X,55(1H-),/1X,'I',5X,'X(I)',T14,
31HI,T15,'Y(J) =')
917 FORMAT(1X,'ENTER I,J,XI,YJ,ZIJ SEPARATED BY COMMAS.'/)
904 FORMAT(1X,'ERROR IN INPUT,TRY AGAIN.'/)
NCOL=5
DO 928 I=1,2
ZX(I)=0
ZY(I)=0
ZXY(I)=0
DO 928 J=1,3
ZAB(I,J)=0
928 CONTINUE
DO 929 I=1,5
929 ZB(I)=0
DO 930 I=1,4
DO 930 J=1,2
930 ZA(I,J)=0
K=0
NOPTS=LX*LY
ISTD=1
IF(FMT.NE.1)GO TO 905
CALL GETFOR(IRSP,IDLG,IFMT,ISTD,96,4)
905 IF(ISTD.EQ.1)IFMT(1)='(2I,3'
IF(ISTD.EQ.1)IFMT(2)='F)'
IF(IDVI.EQ.'TTY')GO TO 907
WRITE(IDLG,906)
GO TO 908
907 WRITE(IDLG,917)
908 K=K+1
909 READ(NDEVI,IFMT,ERR=910)I1,J1,XI,YJ,ZIJ
I=I1
J=J1
X(I)=XI
Y(J)=YJ
IJ=(I-1)*LY+J
Z(IJ)=ZIJ
IF(K.NE.NOPTS)GO TO 908
GO TO 1
911 WRITE(NDEVO,912)(IDENT(I),I=1,4)
WRITE(NDEVO,913)
NFIRST=1
NLAST=NCOL
IF(NLAST.GT.NV)NLAST=NV
925 IF(NFIRST.LE.NCOL)GO TO 918
WRITE(NDEVO,924)(I,I=NFIRST,NLAST)
WRITE(NDEVO,926)
918 WRITE(NDEVO,919)(V(I),I=NFIRST,NLAST)
WRITE(NDEVO,920)
DO 921 I=1,NU
IJ=(I-1)*NV
WRITE(NDEVO,922)I,U(I),(W(IJ+J),J=NFIRST,NLAST)
921 CONTINUE
IF(NLAST.GE.NV)GO TO 923
NFIRST=NFIRST+NCOL
NLAST=NLAST+NCOL
IF(NLAST.GT.NV)NLAST=NV
GO TO 925
923 WRITE(NDEVO,916)
RETURN
910 IF(ICODE.EQ.-1)CALL EXIT
WRITE(IDLG,904)
GO TO 909
1 IU0 = IU
LX0 = LX
LXM1 = LX0 - 1
LXM2 = LXM1 - 1
LY0 = LY
LYM1 = LY0 - 1
LYM2 = LYM1 - 1
MX0 = MX
MXP1 = MX0 + 1
MXM1 = MX0 - 1
MY0 = MY
MYP1 = MY0 + 1
MYM1 = MY0 - 1
NU0 = NU
NV0 = NV
C ERROR CHECK
IF (LXM2.LT.0) GO TO 400
IF (LYM2.LT.0) GO TO 410
IF (MXM1.LE.0) GO TO 420
IF (MYM1.LE.0) GO TO 430
IF (NU0.NE.LXM1*MX0+1) GO TO 440
IF (NV0.NE.LYM1*MY0+1) GO TO 450
IX = 2
IF (X(1)-X(2)) 10, 460, 30
10 DO 20 IX=3,LX0
IF (X(IX-1)-X(IX)) 20, 460, 470
20 CONTINUE
GO TO 50
30 DO 40 IX=3,LX0
IF (X(IX-1)-X(IX)) 470, 460, 40
40 CONTINUE
50 IY = 2
IF (Y(1)-Y(2)) 60, 490, 80
60 DO 70 IY=3,LY0
IF (Y(IY-1)-Y(IY)) 70, 490, 500
70 CONTINUE
GO TO 100
80 DO 90 IY=3,LY0
IF (Y(IY-1)-Y(IY)) 500, 490, 90
90 CONTINUE
C COMPUTATION OF THE U ARRAY
100 FMX = MX0
RMX = 1.0/FMX
KU = 1
X4 = X(1)
U(1) = X4
DO 120 IX=2,LX0
X3 = X4
X4 = X(IX)
DU = (X4-X3)*RMX
DO 110 JX=1,MXM1
KU = KU + 1
U(KU) = U(KU-1) + DU
110 CONTINUE
KU = KU + 1
U(KU) = X4
120 CONTINUE
C COMPUTATION OF THE V ARRAY
FMY = MY0
RMY = 1.0/FMY
KV = 1
Y4 = Y(1)
V(1) = Y4
DO 140 IY=2,LY0
Y3 = Y4
Y4 = Y(IY)
DV = (Y4-Y3)*RMY
DO 130 JY=1,MYM1
KV = KV + 1
V(KV) = V(KV-1) + DV
130 CONTINUE
KV = KV + 1
V(KV) = Y4
140 CONTINUE
C MAIN DO-LOOPS
JYMX = MY0
KV0 = 0
DO 390 IY=2,LY0
IYM2 = IY - 2
IYM3 = IYM2 - 1
IYML = IY - LY0
IYML1 = IYML + 1
IX6 = 0
IF (IYML.EQ.0) JYMX = MYP1
JXMX = MX0
KU0 = 0
DO 380 IX=1,LX0
IXM1 = IX - 1
IXML = IX - LX0
IF (IXML.EQ.0) JXMX = MXP1
C ROUTINES TO PICK UP NECESSARY X,Y, AND Z VALUES, TO
C COMPUTE THE ZA, ZB, AND ZAB VLAUES, AND TO ESTIMATE THEM
C WHEN NECESSARY
C PRELIMINARY WHEN IX.EQ.1
IF (IXM1.NE.0) GO TO 150
Y3 = Y(IY-1)
Y4 = Y(IY)
B3 = 1.0/(Y4-Y3)
B3SQ = B3*B3
IF (IYM2.GT.0) B2 = 1.0/(Y3-Y(IY-2))
IF (IYM3.GT.0) B1 = 1.0/(Y(IY-2)-Y(IY-3))
IF (IYML.LT.0) B4 = 1.0/(Y(IY+1)-Y4)
IF (IYML1.LT.0) B5 = 1.0/(Y(IY+2)-Y(IY+1))
GO TO 180
C TO SAVE THE OLD VALUES
150 Z3A2 = Z3A3
Z4A2 = Z4A3
X3 = X4
Z33 = Z43
Z3B3 = Z4B3
A3 = A4
A3SQ = A3*A3
Z3A3 = Z3A4
Z4A3 = Z4A4
ZA3B2 = ZA4B2
ZA3B3 = ZA4B3
ZA3B4 = ZA4B4
160 X4 = X5
Z43 = Z53
Z4B1 = Z5B1
Z4B2 = Z5B2
Z4B3 = Z5B3
Z4B4 = Z5B4
Z4B5 = Z5B5
A4 = A5
Z3A4 = Z3A5
Z4A4 = Z4A5
ZA4B2 = ZA5B2
ZA4B3 = ZA5B3
ZA4B4 = ZA5B4
170 X5 = X6
Z53 = Z63
Z54 = Z64
Z5B1 = Z6B1
Z5B2 = Z6B2
Z5B3 = Z6B3
Z5B4 = Z6B4
Z5B5 = Z6B5
C TO COMPUTE THE ZA,ZB, AND ZAB VALUES AND
C TO ESTIMATE THE ZB VALUES
C WHEN (IY.LE.3).OR.(IY.GE.LY-1)
180 IX6 = IX6 + 1
IF (IX6.GT.LX0) GO TO 260
X6 = X(IX6)
JZ63=(IX6-1)*LY+IY-1
Z63=Z(JZ63)
JZ64=JZ63+1
Z64=Z(JZ64)
Z6B3 = (Z64-Z63)*B3
IF (LYM2.EQ.0) GO TO 200
IF (IYM2.EQ.0) GO TO 190
JZ62=JZ64-2
Z62=Z(JZ62)
Z6B2 = (Z63-Z62)*B2
IF (IYML.NE.0) GO TO 190
Z6B4 = Z6B3 + Z6B3 - Z6B2
GO TO 210
190 JZ65=JZ64+1
Z65=Z(JZ65)
Z6B4 = (Z65-Z64)*B4
IF (IYM2.NE.0) GO TO 210
Z6B2 = Z6B3 + Z6B3 - Z6B4
GO TO 210
200 Z6B2 = Z6B3
Z6B4 = Z6B3
210 IF (IYM3.LE.0) GO TO 220
J6B1=(IX6-1)*LY+IY-3
Z6B1=(Z62-Z(J6B1))*B1
GO TO 230
220 Z6B1 = Z6B2 + Z6B2 - Z6B3
230 IF (IYML1.GE.0) GO TO 240
J6B5=(IX6-1)*LY+IY+2
Z6B5=(Z(J6B5)-Z65)*B5
GO TO 250
240 Z6B5 = Z6B4 + Z6B4 - Z6B3
250 IF (IX6.EQ.1) GO TO 170
A5 = 1.0/(X6-X5)
Z3A5 = (Z63-Z53)*A5
Z4A5 = (Z64-Z54)*A5
ZA5B2 = (Z6B2-Z5B2)*A5
ZA5B3 = (Z6B3-Z5B3)*A5
ZA5B4 = (Z6B4-Z5B4)*A5
IF (IX6.EQ.2) GO TO 160
GO TO 280
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN (IX.GE.LX-1).AND.(LX.GT.2)
260 IF (LXM2.EQ.0) GO TO 270
Z3A5 = Z3A4 + Z3A4 - Z3A3
Z4A5 = Z4A4 + Z4A4 -Z4A3
IF (IXML.EQ.0) GO TO 290
ZA5B2 = ZA4B2 + ZA4B2 - ZA3B2
ZA5B3 = ZA4B3 + ZA4B3 - ZA3B3
ZA5B4 = ZA4B4 + ZA4B4 - ZA3B4
GO TO 290
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN (IX.GE.LX-1).AND.(LX.EQ.2)
270 Z3A5 = Z3A4
Z4A5 = Z4A4
IF (IXML.EQ.0) GO TO 290
ZA5B2 = ZA4B2
ZA5B3 = ZA4B3
ZA5B4 = ZA4B4
C TO ESTIMATE THE ZA AND ZAB VALUES
C WHEN IX.EQ.1
280 IF (IXM1.NE.0) GO TO 290
Z3A3 = Z3A4 + Z3A4 - Z3A5
Z3A2 = Z3A3 + Z3A3 - Z3A4
Z4A3 = Z4A4 + Z4A4 - Z4A5
Z4A2 = Z4A3 + Z4A3 - Z4A4
ZA3B2 = ZA4B2 + ZA4B2 - ZA5B2
ZA3B3 = ZA4B3 + ZA4B3 - ZA5B3
ZA3B4 = ZA4B4 + ZA4B4 - ZA5B4
GO TO 300
C NUMERICAL DIFFERENTIATION --- TO DETERMINE PARTIAL
C DERIVATIVES ZX, ZY, AND ZXY AS WEIGHTED MEANS OF DIVIDED
C DIFFERENCES ZA, ZB, AND ZAB, RESPECTIVELY
C TO SAVE THE OLD VALUES WHEN IX.NE.1
290 ZX33 = ZX43
ZX34 = ZX44
ZY33 = ZY43
ZY34 = ZY44
ZXY33 = ZXY43
ZXY34 = ZXY44
C NEW COMPUTATION
300 DO 350 JY=1,2
W2 = ABS(ZA(4,JY)-ZA(3,JY))
W3 = ABS(ZA(2,JY)-ZA(1,JY))
SW = W2 + W3
IF (SW.EQ.0.0) GO TO 310
WX2 = W2/SW
WX3 = W3/SW
GO TO 320
310 WX2 = 0.5
WX3 = 0.5
320 ZX(JY) = WX2*ZA(2,JY) + WX3*ZA(3,JY)
W2 = ABS(ZB(JY+3)-ZB(JY+2))
W3 = ABS(ZB(JY+1)-ZB(JY))
SW = W2 + W3
IF (SW.EQ.0.0) GO TO 330
WY2 = W2/SW
WY3 = W3/SW
GO TO 340
330 WY2 = 0.5
WY3 = 0.5
340 ZY(JY) = WY2*ZB(JY+1) + WY3*ZB(JY+2)
ZXY(JY) = WY2*(WX2*ZAB(1,JY)+WX3*ZAB(2,JY)) +
* WY3*(WX2*ZAB(1,JY+1)+WX3*ZAB(2,JY+1))
350 CONTINUE
IF (IXM1.EQ.0) GO TO 380
C DETERMINATION FO THE COEFFICIENTS OF THE POLYNOMIAL
ZX3B3 = (ZX34-ZX33)*B3
ZX4B3 = (ZX44-ZX43)*B3
ZY3A3 = (ZY43-ZY33)*A3
ZY4A3 = (ZY44-ZY34)*A3
A = ZA3B3 - ZX3B3 - ZY3A3 + ZXY33
B = ZX4B3 - ZX3B3 - ZXY43 + ZXY33
C = ZY4A3 - ZY3A3 - ZXY34 + ZXY33
D = ZXY44 - ZXY43 - ZXY34 + ZXY33
E = A + A - B - C
P02 = (2.0*(Z3B3-ZY33)+Z3B3-ZY34)*B3
P03 = (-2.0*Z3B3+ZY34+ZY33)*B3SQ
P12 = (2.0*(ZX3B3-ZXY33)+ZX3B3-ZXY34)*B3
P13 = (-2.0*ZX3B3+ZXY34+ZXY33)*B3SQ
P20 = (2.0*(Z3A3-ZX33)+Z3A3-ZX43)*A3
P21 = (2.0*(ZY3A3-ZXY33)+ZY3A3-ZXY43)*A3
P22 = (3.0*(A+E)+D)*A3*B3
P23 = (-3.0*E-B-D)*A3*B3SQ
P30=(-2.0*Z3A3+ZX43+ZX33)*A3SQ
P31 = (-2.0*ZY3A3+ZXY43+ZXY33)*A3SQ
P32 = (-3.0*E-C-D)*B3*A3SQ
P33 = (D+E+E)*A3SQ*B3SQ
C COMPUTATION OF THE POLYNOMIAL
DO 370 JY=1,JYMX
KV = KV0 + JY
DY = V(KV) - Y3
Q0 = P00 + DY*(P01+DY*(P02+DY*P03))
Q1 = P10 + DY*(P11+DY*(P12+DY*P13))
Q2 = P20 + DY*(P21+DY*(P22+DY*P23))
Q3 = P30 + DY*(P31+DY*(P32+DY*P33))
DO 360 JX=1,JXMX
KU = KU0 + JX
DX = U(KU) - X3
JKUKV=(KU-1)*NV+KV
W(JKUKV) = Q0 + DX*(Q1+DX*(Q2+DX*Q3))
360 CONTINUE
370 CONTINUE
KU0 = KU0 + MX0
380 CONTINUE
KV0 = KV0 + MY0
390 CONTINUE
C NORMAL EXIT
GO TO 911
C ERROR EXIT
400 WRITE (IU0,99999)
GO TO 520
410 WRITE (IU0,99998)
GO TO 520
420 WRITE (IU0,99997)
GO TO 520
430 WRITE (IU0,99996)
GO TO 520
440 WRITE (IU0,99995)
GO TO 520
450 WRITE (IU0,99994)
GO TO 520
460 WRITE (IU0,99993)
GO TO 480
470 WRITE (IU0,99992)
480 WRITE (IU0,99991) IX, X(IX)
GO TO 520
490 WRITE (IU0,99990)
GO TO 510
500 WRITE (IU0,99989)
510 WRITE (IU0,99988) IY, Y(IY)
520 WRITE (IU0,99987) LX0, MX0, NU0, LY0, MY0, NV0
RETURN
C FORMAT STATEMENTS
99999 FORMAT(1X/23H *** LX = 1 OR LESS./)
99998 FORMAT(1X/23H *** LY = 1 OR LESS./)
99997 FORMAT(1X/22H *** N = 0 OR LESS./)
99996 FORMAT(1X/27H *** IDENTICAL X VALUES./)
99995 FORMAT(1X/33H *** X VALUES OUT OF SEQUENCE./)
99994 FORMAT(7H IX =, I6, 10X, 7HX(IX) =, E12.3)
99993 FORMAT(1X/27H *** IDENTICAL Y VALUES./)
99992 FORMAT(1X/33H *** Y VALUES OUT OF SEQUENCE./)
99991 FORMAT(7H IY =, I6, 10X, 7HY(IY) =, E12.3)
99990 FORMAT(1X/27H *** IDENTICAL Y VALUES./)
99989 FORMAT(1X/33H *** Y VALUES OUT OF SEQUENCE./)
99988 FORMAT(7H IY =, I6, 10X, 7HY(IY) =, E12.3)
99987 FORMAT(7H LX =, I6, 10X 4HMX =, I6, 10X, 4HNU =, I6/
* 7H LY =, I6, 10X, 4HMY =, I6, 10X, 4HNV =, I6/6H ERROR,
* 30H DETECTED IN ROUTINE SFCFIT)
END