SUBROUTINE RAT2D(U,X,Y,NXP,NYP,NX,NY,PP,A,IERR) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C Computes for given functional values U(I,J) at points X(I) and C Y(J) the coefficients A(I,J,K,L) of a birational spline C interpolation. The boundary points are determined through a C simple difference scheme. C C INPUT C U Data matrix (NXP,NYP). C X,Y Coordinates of the data. C NXP,NYP Physical length of the above. C NX,NY Length of the above to be processed. C PP Parameter of the birational spline interpolation. C PP > -1. C PP = 0: Bicubic Spline interpolation. C PP --> + Infinity: Towards bilinear interpolation. C IERR Flag for matrix dimensions C = 1: wrong matrix dimension C C OUTPUT C A Coefficient matrix (NX,NY,4,4) of the birational C spline interpolation. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (EPS=1.E-10,NNX=401,NNY=401,NB=4) C REAL U(NXP,NYP),X(NXP),Y(NYP),A(NXP,NYP,4,4) C ! REAL P(NNX,NNY),Q(NNX,NNY),R(NNX,NNY) ! REAL AX(NNX),BX(NNX),CX(NNX),RX(NNX) ! REAL AY(NNY),BY(NNY),CY(NNY),RY(NNY) ! REAL DX(NNX),DY(NNY) REAL P(NXP,NYP),Q(NXP,NYP),R(NXP,NYP) REAL AX(NXP),BX(NXP),CX(NXP),RX(NXP) REAL AY(NYP),BY(NYP),CY(NYP),RY(NYP) REAL DX(NXP),DY(NYP) C REAL B(NB,NB),C(NB,NB),D(NB,NB),E(NB,NB) C C...PARAMETERS AND CHECK OF ARRAY BOUNDS C ! MMX=NNX ! MMY=NNY MMX=NXP MMY=NYP C IF ((NX.LT.2).OR.(NY.LT.2)) THEN PRINT*,' RAT2D: NX, NY SMALLER THAN TWO, FULL STOP' PRINT*,' NX, NY: ',NX,NY IERR=1 STOP ENDIF C C...PARAMETERS AND NONEQUIDISTANT GRID C NX1=NX-1 NX2=NX-2 NY1=NY-1 NY2=NY-2 II=MOD(NX1,4) DO 10 I=1,II DX(I)=1./(X(I+1)-X(I)) 10 CONTINUE DO 11 I=1+II,NX1,4 DX(I)=1./(X(I+1)-X(I)) DX(I+1)=1./(X(I+2)-X(I+1)) DX(I+2)=1./(X(I+3)-X(I+2)) DX(I+3)=1./(X(I+4)-X(I+3)) 11 CONTINUE JJ=MOD(NY1,4) DO 12 J=1,JJ DY(J)=1./(Y(J+1)-Y(J)) 12 CONTINUE DO 13 J=1+JJ,NY1,4 DY(J)=1./(Y(J+1)-Y(J)) DY(J+1)=1./(Y(J+2)-Y(J+1)) DY(J+2)=1./(Y(J+3)-Y(J+2)) DY(J+3)=1./(Y(J+4)-Y(J+3)) 13 CONTINUE P1=PP+1. P2=PP+2. P3=PP+3. DK=1./P1 B(1,1)=P2*DK B(2,1)=-DK B(4,3)=-DK B(4,1)=DK B(2,3)=B(1,1) B(1,3)=-DK B(3,1)=-DK B(3,3)=DK E(1,1)=B(1,1) E(2,1)=B(2,1) E(3,1)=B(3,1) E(4,1)=B(4,1) E(1,3)=B(1,3) E(2,3)=B(2,3) E(3,3)=B(3,3) E(4,3)=B(4,3) C C...COMPUTE PARTIAL DERIVATIVES AT THE BOUNDARIES C JJ=MOD(NY,4) DO 49 J=1,JJ P(1,J)=(U(2,J)-U(1,J))*DX(1) P(NX,J)=(U(NX,J)-U(NX1,J))*DX(NX1) 49 CONTINUE TEMP=DX(1) DO 50 J=1+JJ,NY,4 P(1,J)=(U(2,J)-U(1,J))*TEMP P(1,J+1)=(U(2,J+1)-U(1,J+1))*TEMP P(1,J+2)=(U(2,J+2)-U(1,J+2))*TEMP P(1,J+3)=(U(2,J+3)-U(1,J+3))*TEMP 50 CONTINUE TEMP=DX(NX1) DO 51 J=1+JJ,NY,4 P(NX,J)=(U(NX,J)-U(NX1,J))*TEMP P(NX,J+1)=(U(NX,J+1)-U(NX1,J+1))*TEMP P(NX,J+2)=(U(NX,J+2)-U(NX1,J+2))*TEMP P(NX,J+3)=(U(NX,J+3)-U(NX1,J+3))*TEMP 51 CONTINUE II=MOD(NX,4) DO 59 I=1,II Q(I,1)=(U(I,2)-U(I,1))*DY(1) Q(I,NY)=(U(I,NY)-U(I,NY1))*DY(NY1) 59 CONTINUE TEMP=DY(1) DO 60 I=1+II,NX,4 Q(I,1)=(U(I,2)-U(I,1))*TEMP Q(I+1,1)=(U(I+1,2)-U(I+1,1))*TEMP Q(I+2,1)=(U(I+2,2)-U(I+2,1))*TEMP Q(I+3,1)=(U(I+3,2)-U(I+3,1))*TEMP 60 CONTINUE TEMP=DY(NY1) DO 61 I=1+II,NX,4 Q(I,NY)=(U(I,NY)-U(I,NY1))*TEMP Q(I+1,NY)=(U(I+1,NY)-U(I+1,NY1))*TEMP Q(I+2,NY)=(U(I+2,NY)-U(I+2,NY1))*TEMP Q(I+3,NY)=(U(I+3,NY)-U(I+3,NY1))*TEMP 61 CONTINUE R(1,1)=.5*((P(1,2)-P(1,1))*DY(1)+ & (Q(2,1)-Q(1,1))*DX(1)) R(1,NY)=.5*((P(1,NY)-P(1,NY1))*DY(NY1)+ & (Q(2,NY)-Q(1,NY))*DX(1)) R(NX,1)=.5*((P(NX,2)-P(NX,1))*DY(1)+ & (Q(NX,1)-Q(NX1,1))*DX(NX1)) R(NX,NY)=.5*((P(NX,NY)-P(NX,NY1))*DY(NY1)+ & (Q(NX,NY)-Q(NX1,NY))*DX(NX1)) C C...PREPARE THE COEFFICIENT MATRIX FOR X C 70 IF ((NX.EQ.2).AND.(NY.EQ.2)) GOTO 130 IF (NX.EQ.2) GOTO 90 II=MOD(NX2-1,4) DO 79 I=1,II AX(I)=DX(I+1) BX(I)=P2*(DX(I+1)+DX(I)) IF (ABS(BX(I)).LT.EPS) GOTO 82 79 CONTINUE DO 80 I=1+II,NX2-1,4 AX(I)=DX(I+1) AX(I+1)=DX(I+2) AX(I+2)=DX(I+3) AX(I+3)=DX(I+4) 80 CONTINUE DO 81 I=1+II,NX2-1,4 BX(I)=P2*(DX(I+1)+DX(I)) BX(I+1)=P2*(DX(I+2)+DX(I+1)) BX(I+2)=P2*(DX(I+3)+DX(I+2)) BX(I+3)=P2*(DX(I+4)+DX(I+3)) IF ((ABS(BX(I)).LT.EPS).OR.(ABS(BX(I+1)).LT.EPS).OR. & (ABS(BX(I+2)).LT.EPS).OR.(ABS(BX(I+3)).LT.EPS)) GOTO 82 81 CONTINUE GOTO 83 82 PRINT*,' RAT2D: BX(',I,') TOO SMALL, FULL STOP' STOP 83 CONTINUE BX(NX2)=P2*(DX(NX2+1)+DX(NX2)) C C...LU- DECOMPOSITION FOR X C DO 84 K=2,NX2 KM1=K-1 H=AX(KM1)/BX(KM1) CX(KM1)=H BX(K)=BX(K)-H*AX(KM1) 84 CONTINUE IVJ=0 IPU=1 CALL RATPE(U,NXP,NYP,P,RX,MMX,MMY,IVJ,IPU,P3,AX,BX,CX,DX,NX,NY) 90 IF (NY.EQ.2) GOTO 110 C C...PREPARE THE COEFFICIENT MATRIX FOR Y C JJ=MOD(NY2-1,4) DO 99 J=1,JJ AY(J)=DY(J+1) BY(J)=P2*(DY(J+1)+DY(J)) IF (ABS(BY(J)).LT.EPS) GOTO 102 99 CONTINUE DO 100 J=1+JJ,NY2-1,4 AY(J)=DX(J+1) AY(J+1)=DY(J+2) AY(J+2)=DY(J+3) AY(J+3)=DY(J+4) 100 CONTINUE DO 101 J=1+JJ,NY2-1,4 BY(J)=P2*(DY(J+1)+DY(J)) BY(J+1)=P2*(DY(J+2)+DY(J+1)) BY(J+2)=P2*(DY(J+3)+DY(J+2)) BY(J+3)=P2*(DY(J+4)+DY(J+3)) IF ((ABS(BY(J)).LT.EPS).OR.(ABS(BY(J+1)).LT.EPS).OR. & (ABS(BY(J+2)).LT.EPS).OR.(ABS(BY(J+3)).LT.EPS)) GOTO 102 101 CONTINUE GOTO 103 102 PRINT*,' RAT2D: BY(',J,') TOO SMALL, FULL STOP' STOP 103 CONTINUE BY(NY2)=P2*(DY(NY2+1)+DY(NY2)) C C...LU- DECOMPOSITION FOR Y C DO 104 K=2,NY2 KM1=K-1 H=AY(KM1)/BY(KM1) CY(KM1)=H BY(K)=BY(K)-H*AY(KM1) 104 CONTINUE IVJ=1 IPU=1 CALL RATPE(U,NXP,NYP,Q,RY,MMX,MMY,IVJ,IPU,P3,AY,BY,CY,DY,NY,NX) 110 IF (NX.EQ.2) GOTO 120 IVJ=0 IPU=NY1 CALL RATPE(Q,MMX,MMY,R,RX,MMX,MMY,IVJ,IPU,P3,AX,BX,CX,DX,NX,NY) 120 IF (NY.EQ.2) GOTO 130 IVJ=1 IPU=1 CALL RATPE(P,MMX,MMY,R,RY,MMX,MMY,IVJ,IPU,P3,AY,BY,CY,DY,NY,NX) C C...COMPUTE THE COEFFICIENT MATRIX C 130 DO 210 I=1,NX1 I1=I+1 FA=(X(I1)-X(I))/(PP+3.) B(1,2)=B(1,1)*FA B(2,2)=B(2,1)*FA B(3,2)=-B(1,2) B(4,2)=-B(2,2) B(1,4)=B(4,2) B(2,4)=B(2,2)*(PP+2.) B(3,4)=B(2,2) B(4,4)=-B(2,4) DO 200 J=1,NY1 J1=J+1 C(1,1)=U(I,J) C(1,2)=Q(I,J) C(2,1)=P(I,J) C(2,2)=R(I,J) C(1,3)=U(I,J1) C(1,4)=Q(I,J1) C(2,3)=P(I,J1) C(2,4)=R(I,J1) C(3,1)=U(I1,J) C(3,2)=Q(I1,J) C(4,1)=P(I1,J) C(4,2)=R(I1,J) C(3,3)=U(I1,J1) C(3,4)=Q(I1,J1) C(4,3)=P(I1,J1) C(4,4)=R(I1,J1) DO 160 K1=1,NB DO 150 K2=1,NB SUM=0. DO 140 K=1,NB SUM=SUM+B(K1,K)*C(K,K2) 140 CONTINUE D(K1,K2)=SUM 150 CONTINUE 160 CONTINUE FB=(Y(J1)-Y(J))/(PP+3.) E(1,2)=E(1,1)*FB E(2,2)=E(2,1)*FB E(3,2)=-E(1,2) E(4,2)=-E(2,2) E(1,4)=E(4,2) E(2,4)=E(2,2)*(PP+2.) E(3,4)=E(2,2) E(4,4)=-E(2,4) DO 190 K1=1,NB DO 180 K2=1,NB SUM=0. DO 170 K=1,NB SUM=SUM+D(K1,K)*E(K2,K) 170 CONTINUE A(I,J,K1,K2)=SUM 180 CONTINUE 190 CONTINUE 200 CONTINUE 210 CONTINUE RETURN END