SUBROUTINE MODEL(MTEXT,LU) C C APPROXIMATION OF INTERFACES AND VELOCITY DISTRIBUTION C IN INDIVIDUAL LAYERS (B-SPLINE APPROXIMATION) C CHARACTER*80 MTEXT INTEGER NX,NY,NZ,IERR REAL X1(20),Y1(20),Z1(20),WG1(20),WG2(20,20),W(20,20,20),C(1000), * VX(5,20),VY(5,20),VZ(5,20),TEMP(100),SIGMA DIMENSION AUX(12),DEP(6),Y(2),IPR(101) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /AUXX/ MMX(20),MMY(20),MMXY(20) COMMON /EPAR/ aij(21,100),aijx(5,21,100),aijy(5,21,100), * aijz(5,21,100),XX(200),YY(200),ZZ(200),SSIGMA(20), * NXX(20),NYY(20),NZZ(20),NXXYZ(20) COMMON /DENS/ RHO(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NNX(20),NNY(20),BRD(6), 1 NINT,XINTA COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY C NW1=20 NW2=20 MPRINT=0 NINT=2 N1=10 N2=10 READ(LU,*) MPRINT,NINT,N1,N2 WRITE(LOU,101) MPRINT,NINT,N1,N2 IF(LUGRD.NE.0)THEN CALL FACETS(N1,N2,NINT) WRITE(LUGRD,108) WRITE(LUGRD,110) END IF NLAY=NINT-1 C C INPUT FOR 3D INTERFACES C MX2=0 MY2=0 MXY2=0 DO 2 I=1,NINT MX1=MX2+1 MY1=MY2+1 MXY1=MXY2+1 READ(LU,*) MX,MY WRITE(LOU,101) MX,MY NNX(I)=MX NNY(I)=MY MX2=MX1+MX-1 MY2=MY1+MY-1 MXY2=MXY1+MX*MY-1 READ(LU,*)(SX(J),J=MX1,MX2) READ(LU,*)(SY(J),J=MY1,MY2) WRITE(LOU,102)(SX(J),J=MX1,MX2) WRITE(LOU,102)(SY(J),J=MY1,MY2) M1=MXY1 DO 1 L=1,MX M2=M1+MY-1 READ(LU,*)(Z(J),J=M1,M2) WRITE(LOU,102)(Z(J),J=M1,M2) 1 M1=M2+1 C C DETERMINATION OF COEFFICIENTS OF BICUBIC SPLINES C APPROXIMATING INTERFACES C CALL BIAP(MX1,MX,MY1,MY,MXY1) MMX(I)=MX1 MMY(I)=MY1 MMXY(I)=MXY1 2 CONTINUE NB1=MMX(1) NB2=MMX(2)-1 BRD(1)=SX(NB1) BRD(2)=SX(NB2) NB1=MMY(1) NB2=MMY(2)-1 BRD(3)=SY(NB1) BRD(4)=SY(NB2) C C INPUT DATA FOR PRINTER PLOT OF 3-D INTERFACES C DO 6 INTR=1,NINT IF(MPRINT.GE.1) WRITE(LOU,109) INTR READ(LU,*) ZMIN,ZMAX IF(INTR.EQ.1)BRD(5)=ZMIN IF(INTR.EQ.NINT)BRD(6)=ZMAX IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX C C NUMERICAL FORM OF 3-D INTERFACES C ZMM=ZMAX-ZMIN ZMMM=ZMM/10. IF(MPRINT.GE.1) WRITE(LOU,102) ZMIN,ZMAX,ZMMM IF(MPRINT.GE.1) WRITE(LOU,106) BRD(1),BRD(2),BRD(3),BRD(4) DY=(BRD(4)-BRD(3))/FLOAT(N2-1) DX=(BRD(2)-BRD(1))/FLOAT(N1-1) MAUX=0 NDER=1 Y(2)=BRD(3)-DY DO 5 K=1,N2 Y(2)=Y(2)+DY Y(1)=BRD(1)-DX DO 4 L=1,N1 Y(1)=Y(1)+DX CALL DISC(Y,DEP) IF(LUGRD.NE.0)THEN DD=SQRT(1.+DEP(2)*DEP(2)+DEP(3)*DEP(3)) UN1=DEP(2)/DD UN2=DEP(3)/DD UN3=-1./DD KL=L+(K-1)*N1 WRITE(LUGRD,112)KL,Y(1),Y(2),DEP(1),UN1,UN2,UN3,INTR END IF AUX1=10.*(DEP(1)-ZMIN)/ZMM IPR(L)=IFIX(AUX1) IF(AUX1.LT.0.0.OR.AUX1.GT.10) IPR(L)=11 4 CONTINUE C C PRINTER PLOT OF 3-D INTERFACES C IF(MPRINT.GE.1) WRITE(LOU,105)(IPR(L),L=1,N1) 5 CONTINUE C C END OF LOOP OVER ALL INTERFACES C 6 CONTINUE IF(LUGRD.NE.0)WRITE(LUGRD,110) IF(LUGRD.NE.0)WRITE(LUGRD,110) C C INPUT OF ELASTIC PARAMETERS C ISQRT=0 IRHO=0 READ(LU,*)ISQRT,IRHO IF(MPRINT.EQ.0)WRITE(LOU,101)ISQRT,IRHO IF(ISQRT.EQ.0.AND.MPRINT.GT.0) WRITE(LOU,114) IF(ISQRT.NE.0.AND.MPRINT.GT.0) WRITE(LOU,113) IF(IRHO.NE.0)THEN DO 7 L=1,NLAY RHO(L)=1. 7 CONTINUE READ(LU,*)(RHO(L),L=1,NLAY) IF(MPRINT.GE.0)WRITE(LOU,102)(RHO(L),L=1,NLAY) END IF C C INPUT OF MATRIX OF ELASTIC PARAMETERS IN COMPRESSED FORM. C ELEMENTS OF THE MATRIX HAVE TO BE GIVEN IN (KM**2/SEC**2) C WHICH CORRESPONDS TO ELASTIC PARAMETERS DIVIDED BY DENSITY C NXXU=0 NYYU=0 NZZU=0 NXYZU=0 DO 19 L=1,NLAY IANI(L)=1 SIGMA=0. READ(LU,*)IANI(L),SIGMA WRITE(LOU,'(I10,f10.5)')IANI(L),SIGMA SSIGMA(L)=SIGMA NX=1 NY=1 NZ=1 X1(1)=0. Y1(1)=0. Z1(1)=0. READ(LU,*)NX,NY,NZ WRITE(LOU,100)NX,NY,NZ READ(LU,*)(X1(I),I=1,NX) WRITE(LOU,102)(X1(I),I=1,NX) READ(LU,*)(Y1(I),I=1,NY) WRITE(LOU,102)(Y1(I),I=1,NY) READ(LU,*)(Z1(I),I=1,NZ) WRITE(LOU,102)(Z1(I),I=1,NZ) NXYZ=NX*NY*NZ DO 8 I=1,NX II=NXXU+I XX(II)=X1(I) 8 CONTINUE DO 9 I=1,NY II=NYYU+I YY(II)=Y1(I) 9 CONTINUE DO 10 I=1,NZ II=NZZU+I ZZ(II)=Z1(I) 10 CONTINUE IF(IANI(L).NE.0) THEN IF(MPRINT.GE.1)WRITE(LOU,115) L III=1 ELSE IF(MPRINT.GE.1)WRITE(LOU,118) L III=15 END IF DO 18 II=1,21,III DO 12 I=1,NX DO 11 J=1,NY READ(LU,*)(W(I,J,K),K=1,NZ) WRITE(LOU,102)(W(I,J,K),K=1,NZ) IF(ISQRT.NE.0.AND.IANI(L).EQ.0)THEN DO 50 K=1,NZ W(I,J,K)=SQRT(W(I,J,K)) 50 CONTINUE END IF IF(NX.EQ.1.AND.NY.EQ.1.AND.NZ.EQ.1)GO TO 11 IF(NX.EQ.1)THEN IF(NY.EQ.1)THEN DO 51 K=1,NZ WG1(K)=W(1,1,K) 51 CONTINUE END IF IF(NZ.EQ.1)THEN WG1(J)=W(1,J,1) END IF IF(NY.NE.1.AND.NZ.NE.1)THEN DO 52 K=1,NZ WG2(J,K)=W(1,J,K) 52 CONTINUE END IF END IF IF(NY.EQ.1)THEN IF(NZ.EQ.1)THEN WG1(I)=W(I,1,1) END IF IF(NX.NE.1.AND.NZ.NE.1)THEN DO 53 K=1,NZ WG2(I,K)=W(I,1,K) 53 CONTINUE END IF END IF IF(NZ.EQ.1)THEN IF(NX.NE.1.AND.NY.NE.1)THEN WG2(I,J)=W(I,J,1) END IF END IF 11 CONTINUE 12 CONTINUE IF(NX.EQ.1.AND.NY.EQ.1.AND.NZ.EQ.1)THEN AIJ(II,NXYZU+1)=W(1,1,1) DO 54 M=1,5 AIJX(M,II,NXXU+1)=0. AIJY(M,II,NYYU+1)=0. AIJZ(M,II,NZZU+1)=0. 54 CONTINUE GO TO 18 END IF IF(NX.EQ.1)THEN IF(NY.EQ.1)THEN CALL CURVB1 (NZ,Z1,WG1,C,VZ,TEMP,SIGMA,IERR) END IF IF(NZ.EQ.1)THEN CALL CURVB1 (NY,Y1,WG1,C,VY,TEMP,SIGMA,IERR) END IF IF(NY.NE.1.AND.NZ.NE.1)THEN CALL SURFB1 (NY,NZ,Y1,Z1,WG2,NW1,C,VY,VZ,TEMP,SIGMA,IERR) END IF END IF IF(NY.EQ.1)THEN IF(NZ.EQ.1)THEN CALL CURVB1 (NX,X1,WG1,C,VX,TEMP,SIGMA,IERR) END IF IF(NX.NE.1.AND.NZ.NE.1)THEN CALL SURFB1 (NX,NZ,X1,Z1,WG2,NW1,C,VX,VZ,TEMP,SIGMA,IERR) END IF END IF IF(NZ.EQ.1)THEN IF(NX.NE.1.AND.NY.NE.1)THEN CALL SURFB1 (NX,NY,X1,Y1,WG2,NW1,C,VX,VY,TEMP,SIGMA,IERR) END IF END IF IF(NX.NE.1.AND.NY.NE.1.AND.NZ.NE.1)THEN CALL VAL3B1 (NX,NY,NZ,X1,Y1,Z1,W,NW1,NW2,C,VX,VY, * VZ,TEMP,SIGMA,IERR) END IF DO 13 I=1,NXYZ KK=NXYZU+I AIJ(II,KK)=C(I) 13 CONTINUE DO 17 M=1,5 DO 14 I=1,NX KK=NXXU+I AIJX(M,II,KK)=VX(M,I) IF(NX.EQ.1)AIJX(M,II,KK)=0. 14 CONTINUE DO 15 I=1,NY KK=NYYU+I AIJY(M,II,KK)=VY(M,I) IF(NY.EQ.1)AIJY(M,II,KK)=0. 15 CONTINUE DO 16 I=1,NZ KK=NZZU+I AIJZ(M,II,KK)=VZ(M,I) IF(NZ.EQ.1)AIJZ(M,II,KK)=0. 16 CONTINUE 17 CONTINUE 18 CONTINUE NXXU=NXXU+NX NYYU=NYYU+NY NZZU=NZZU+NZ NXYZU=NXYZU+NXYZ NXX(L)=NXXU NYY(L)=NYYU NZZ(L)=NZZU NXXYZ(L)=NXYZU 19 CONTINUE C WRITE(LOU,107)MTEXT C C FORMATS C 100 FORMAT(26I3) 101 FORMAT(14I5) 102 FORMAT(8F10.5) 103 FORMAT(5X,5(F7.3,13X),F7.3) 104 FORMAT(F7.3,1X,101I1) 105 FORMAT(8X,101I1) 106 FORMAT(1X,'(XMIN, XMAX)',2X,2F10.5,5X,'(YMIN, YMAX)',2F10.5) 107 FORMAT(////,3X,A///) 108 FORMAT(10H'VERTICES') 109 FORMAT(///,' INTERFACE NUMBER ',I5) 110 FORMAT('/') 111 FORMAT(6F10.5) 112 FORMAT(5H'VRTX,I3,1H',1X,6F10.5,1X,I2,1X,'/') 113 FORMAT(//' INTERPOLATION IN SQUARE ROOTS OF ELASTIC PARAMETERES'/) 114 FORMAT(//' INTERPOLATION IN VALUES OF ELASTIC PARAMETERES'/) 115 FORMAT(/' LAYER',I4,' IS ANISOTROPIC ',/,' MATRIX OF 1 ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2: 2 IT CONTAINS ELASTIC PARAMETERS/DENSITY'/) 118 FORMAT(/' LAYER',I4,' IS ISOTROPIC'/) RETURN END C C ********************************************************* C SUBROUTINE PARDIS(Y,IAY) C SAVE X1,Y1,Z1,SIGMA,NX,NY,NZ,NXYZ,NXXL,NYYL,NZZL,NXYZL,LLAY DIMENSION X1(20),Y1(20),Z1(20),C(1000), * VX(5,20),VY(5,20),VZ(5,20) DIMENSION Y(18),E(21),EX(21),EY(21),EZ(21) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /APROX1/ D(21),DX(21),DY(21),DZ(21),DXX(21), 1 DXY(21),DXZ(21),DYY(21),DYZ(21),DZZ(21) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /EPAR/ aij(21,100),aijx(5,21,100),aijy(5,21,100), * aijz(5,21,100),XX(200),YY(200),ZZ(200),SSIGMA(20), * NXX(20),NYY(20),NZZ(20),NXXYZ(20) COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) COMMON/DENS/RHO(20) C EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15) 1 ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24) 2 ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35) 2 ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55) 2 ,(E(20),A56),(E(21),A66) 1 ,(EX(1),DXA11),(EX(2),DXA12),(EX(3),DXA13),(EX(4),DXA14) 1 ,(EX(5),DXA15),(EX(6),DXA16),(EX(7),DXA22),(EX(8),DXA23) 1 ,(EX(9),DXA24),(EX(10),DXA25),(EX(11),DXA26),(EX(12),DXA33) 1 ,(EX(13),DXA34),(EX(14),DXA35),(EX(15),DXA36),(EX(16),DXA44) 1 ,(EX(17),DXA45),(EX(18),DXA46),(EX(19),DXA55),(EX(20),DXA56) 1 ,(EX(21),DXA66) EQUIVALENCE 1 (EY(1),DYA11),(EY(2),DYA12),(EY(3),DYA13),(EY(4),DYA14) 1 ,(EY(5),DYA15),(EY(6),DYA16),(EY(7),DYA22),(EY(8),DYA23) 1 ,(EY(9),DYA24),(EY(10),DYA25),(EY(11),DYA26),(EY(12),DYA33) 1 ,(EY(13),DYA34),(EY(14),DYA35),(EY(15),DYA36),(EY(16),DYA44) 1 ,(EY(17),DYA45),(EY(18),DYA46),(EY(19),DYA55),(EY(20),DYA56) 1 ,(EY(21),DYA66) 1 ,(EZ(1),DZA11),(EZ(2),DZA12),(EZ(3),DZA13),(EZ(4),DZA14) 1 ,(EZ(5),DZA15),(EZ(6),DZA16),(EZ(7),DZA22),(EZ(8),DZA23) 1 ,(EZ(9),DZA24),(EZ(10),DZA25),(EZ(11),DZA26),(EZ(12),DZA33) 1 ,(EZ(13),DZA34),(EZ(14),DZA35),(EZ(15),DZA36),(EZ(16),DZA44) 1 ,(EZ(17),DZA45),(EZ(18),DZA46),(EZ(19),DZA55),(EZ(20),DZA56) 1 ,(EZ(21),DZA66) C if(lay.eq.llay)go to 4 sigma=ssigma(lay) if(lay.gt.1)then nxxl=nxx(lay-1) nyyl=nyy(lay-1) nzzl=nzz(lay-1) nxyzl=nxxyz(lay-1) else nxxl=0 nyyl=0 nzzl=0 nxyzl=0 end if nxxu=nxx(lay) nyyu=nyy(lay) nzzu=nzz(lay) nxyzu=nxxyz(lay) NX=nxxu-nxxl do 1 i=1,nx ii=nxxl+i x1(i)=xx(ii) 1 continue NY=nyyu-nyyl do 2 i=1,ny ii=nyyl+i y1(i)=yy(ii) 2 continue NZ=nzzu-nzzl do 3 i=1,nz ii=nzzl+i z1(i)=zz(ii) 3 continue nxyz=nxyzu-nxyzl 4 continue IF(IANI(Lay).NE.0) THEN iii=1 else iii=15 end if do 10 ii=1,21,iii WX=0. WY=0. WZ=0. WXX=0. WXY=0. WXZ=0. WYY=0. WYZ=0. WZZ=0. do 5 i=1,nxyz kk=nxyzl+i c(i)=aij(ii,kk) 5 continue do 9 m=1,5 do 6 i=1,nx kk=nxxl+i vx(m,i)=aijx(m,ii,kk) 6 continue do 7 i=1,ny kk=nyyl+i vy(m,i)=aijy(m,ii,kk) 7 continue do 8 i=1,nz kk=nzzl+i vz(m,i)=aijz(m,ii,kk) 8 continue 9 continue if(nx.eq.1.and.ny.eq.1.and.nz.eq.1)then w=c(1) go to 20 end if if(nx.eq.1)then if(ny.eq.1)then call CURVBD (Y(3),W,WZ,WZZ,NZ,Z1,C,VZ,SIGMA) end if if(nz.eq.1)then call CURVBD (Y(2),W,WY,WYY,NY,Y1,C,VY,SIGMA) end if if(ny.ne.1.and.nz.ne.1)then call SURFBD (Y(2),Y(3),W,WY,WZ,WYY,WYZ,WZZ,NY,NZ,Y1, * Z1,C,VY,VZ,SIGMA) end if end if if(ny.eq.1)then if(nz.eq.1)then call CURVBD (Y(1),W,WX,WXX,NX,X1,C,VX,SIGMA) end if if(nx.ne.1.and.nz.ne.1)then call SURFBD (Y(1),Y(3),W,WX,WZ,WXX,WXZ,WZZ,NX,NZ,X1, * Z1,C,VX,VZ,SIGMA) end if end if if(nz.eq.1)then if(nx.ne.1.and.ny.ne.1)then call SURFBD (Y(1),Y(2),W,WX,WY,WXX,WXY,WYY,NX,NY,X1, * Y1,C,VX,VY,SIGMA) end if end if if(nx.ne.1.and.ny.ne.1.and.nz.ne.1)then call VAL3BD(Y(1),Y(2),Y(3),W,WX,WY,WZ,WXX,WXY,WYY,WYZ,WZZ,WXZ, * NX,NY,NZ,X1,Y1,Z1,C,VX,VY,VZ,SIGMA) end if 20 E(ii)=W EX(ii)=WX EY(ii)=WY EZ(ii)=WZ D(ii)=W DX(ii)=WX DY(ii)=WY DZ(ii)=WZ DXX(ii)=WXX DXY(ii)=WXY DXZ(ii)=WXZ DYY(ii)=WYY DYZ(ii)=WYZ DZZ(ii)=WZZ if(isqrt.ne.0.and.iani(lay).eq.0)then EE=2.*W E(ii)=W*W EX(ii)=WX*EE EY(ii)=WY*EE EZ(ii)=WZ*EE D(ii)=W*W DX(ii)=WX*EE DY(ii)=WY*EE DZ(ii)=WZ*EE DXX(ii)=WXX*EE+2.*WX*WX DXY(ii)=WXY*EE+2.*WX*WY DXZ(ii)=WXZ*EE+2.*WX*WZ DYY(ii)=WYY*EE+2.*WY*WY DYZ(ii)=WYZ*EE+2.*WY*WZ DZZ(ii)=WZZ*EE+2.*WZ*WZ end if 10 continue llay=lay IF(IANI(LAY).EQ.0) GOTO 12 A2546=A25+A46 A1266=A12+A66 A1355=A13+A55 A1456=A14+A56 A3645=A36+A45 A2344=A23+A44 IF(IAY.EQ.0)RETURN DO 11 I=1,21 11 AY(I+7,N)=E(I) RETURN 12 IF(IAY.EQ.0)RETURN AY(8,N)=A11 AY(9,N)=DXA11 AY(10,N)=DYA11 AY(11,N)=DZA11 AY(12,N)=A44 AY(13,N)=DXA44 AY(14,N)=DYA44 AY(15,N)=DZA44 RO=1.7+0.2*SQRT(A11) IF(IRHO.NE.0)RO=RHO(LAY) AY(16,N)=RO RETURN END C C ********************************************************* C C SUBROUTINES OF THE SOFTWARE PACKAGE 'FITPACK' BY A.K. CLINE C USED TO SPECIFY THE MODEL FOR THE COMPLETE RAY TRACING ALGORITHM. C C THIS FILE CONSISTS OF THE FOLLOWING PARTS: C (0) AUXILIARY SUBROUTINE C SNHCSH C COMMON TO ALL THE FOLLOWING PARTS. C (1) THE SUBROUTINES PREPARING THE PARAMETERS NECESSARY TO COMPUTE C AN INTERPOLATORY FUNCTION: C CURVB1 (B-SPLINE REPRESENTATION OF 1-D FUNCTION), C SURFB1 (B-SPLINE REPRESENTATION OF 2-D FUNCTION), C VAL3B1 (B-SPLINE REPRESENTATION OF 3-D FUNCTION), C VGEN (AUXILIARY SUBROUTINE), C TERMS (AUXILIARY SUBROUTINE), C TRIDEC (AUXILIARY SUBROUTINE), C TRISOL (AUXILIARY SUBROUTINE). C (2) THE SUBROUTINES EVALUATING THE VALUE, FIRST AND SECOND PARTIAL C DERIVATIVES OF THE INTERPOLATORY FUNCTION AT A GIVEN POINT: C CURVBD (B-SPLINE REPRESENTATION OF 1-D FUNCTION), C SURFBD (B-SPLINE REPRESENTATION OF 2-D FUNCTION), C VAL3BD (B-SPLINE REPRESENTATION OF 3-D FUNCTION), C DSPLNZ (AUXILIARY SUBROUTINE), C INTRVL (AUXILIARY EXTERNAL FUNCTION). C C TAKEN FROM: C FITPACK - A SOFTWARE PACKAGE FOR CURVE AND SURFACE FITTING C EMPLOYING SPLINES UNDER TENSION C BY ALAN KAYLOR CLINE, DEPARTMENT OF COMPUTER SCIENCES, C THE UNIVERSITY OF TEXAS AT AUSTIN, AUGUST 31, 1981. C NOTE 1: C TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS (1) C HAVE BEEN CHANGED TO (*), AND SUBROUTINE TRISOL HAS BEEN REVISED. C NOTE 2: C SUBROUTINES CURVB1 AND CURVBD DO NOT BELONG TO THE ORIGINAL C VERSION OF FITPACK. C NOTE 3: C TO GET THE ORIGINAL VERSIONS OF THE SUBROUTINES SURFBD AND VAL3BD, C THE STATEMENT WITH 'CALL VAR2' MUST BE REMOVED FROM EACH OF THEM. C THE STATEMENTS HAVE BEEN ADDED BY L.KLIMES FOR THE SAKE OF INVERSE C MODELLING TO THE SUBROUTINES CURVBD, SURFBD, AND VAL3BD. C THE THREE APPEARENCES OF THE STATEMENTS 'CALL VAR2' ARE DENOTED BY C 'V' IN THE FIRST COLUMN. THE THREE LINES SHOULD BE REMOVED OR C MODIFIED BEFORE COMPILATION. C C======================================================================= C C PART 0: C C======================================================================= C SUBROUTINE SNHCSH (SINHM,COSHM,X,ISW) C INTEGER ISW REAL SINHM,COSHM,X C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY A. K. CLINE AND R. J. RENKA C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE RETURNS APPROXIMATIONS TO C SINHM(X) = SINH(X)-X C COSHM(X) = COSH(X)-1 C AND C COSHMM(X) = COSH(X)-1-X*X/2 C WITH RELATIVE ERROR LESS THAN 6.16E-6 C C ON INPUT-- C C X CONTAINS THE VALUE OF THE INDEPENDENT VARIABLE. C C ISW INDICATES THE FUNCTION DESIRED C = -1 IF ONLY SINHM IS DESIRED, C = 0 IF BOTH SINHM AND COSHM ARE DESIRED, C = 1 IF ONLY COSHM IS DESIRED, C = 2 IF ONLY COSHMM IS DESIRED, C = 3 IF BOTH SINHM AND COSHMM ARE DESIRED. C C ON OUTPUT-- C C SINHM CONTAINS THE VALUE OF SINHM(X) IF ISW .LE. 0 OR C ISW .EQ. 3 (SINHM IS UNALTERED IF ISW .EQ.1 OR ISW .EQ. C 2). C C COSHM CONTAINS THE VALUE OF COSHM(X) IF ISW .EQ. 0 OR C ISW .EQ. 1 AND CONTAINS THE VALUE OF COSHMM(X) IF ISW C .GE. 2 (COSHM IS UNALTERED IF ISW .EQ. -1). C C AND C C X AND ISW ARE UNALTERED. C C----------------------------------------------------------- C DATA SP2/5.04850926418006E-04/, * SP1/3.62841692246321E-02/, * SQ1/-1.37157937097122E-02/ DATA CP2/1.31625490355985E-03/, * CP1/6.57866547762733E-02/, * CQ1/-1.75465241841312E-02/ DATA ZP2/1.40048186158693E-04/, * ZP1/1.67309141907440E-02/, * ZQ2/9.82154460147143E-05/, * ZQ1/-1.66024148976133E-02/ XX = X AX = ABS(XX) XS = XX*XX IF ((AX .GE. 2.20) .OR. (AX .GE. 1.17 .AND. * ISW .NE. 2)) EXPX = EXP(AX) C C SINHM APPROXIMATION C IF (ISW .EQ. 1 .OR. ISW .EQ. 2) GO TO 2 IF (AX .GE. 1.17) GO TO 1 SINHM = (((SP2*XS+SP1)*XS+1.)*XS*XX)/((SQ1*XS+1.)*6.) GO TO 2 1 SINHM = (EXPX-1./EXPX)/2.-AX IF (XX .LT. 0.) SINHM = -SINHM C C COSHM APPROXIMATION C 2 IF (ISW .NE. 0 .AND. ISW .NE. 1) GO TO 4 IF (AX .GE. 1.17) GO TO 3 COSHM = (((CP2*XS+CP1)*XS+1.)*XS)/((CQ1*XS+1.)*2.) GO TO 4 3 COSHM = (EXPX+1./EXPX)/2.-1. C C COSHMM APPROXIMATION C 4 IF (ISW .LE. 1) RETURN IF (AX .GE. 2.20) GO TO 5 COSHM = (((ZP2*XS+ZP1)*XS+1.)*XS*XS)/(((ZQ2*XS+ZQ1)*XS * +1.)*24.) RETURN 5 COSHM = (EXPX+1./EXPX)/2.-1.-XS/2. RETURN END C C======================================================================= C C PART 1: C C======================================================================= C SUBROUTINE CURVB1 (NX,X,W,C,VX,TEMP,SIGMA,IERR) C INTEGER NX,IERR REAL X(NX),W(NX),C(NX),VX(5,NX),TEMP(*),SIGMA C C COMPLEMENT TO FITPACK C BY ALAN KAYLOR CLINE C CODED -- OCTOBER 9, 1986 C BY LUDEK KLIMES C INST. GEOL. GEOTECHN. C CZECHOSL. ACAD. SCI., PRAGUE C C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO C COMPUTE AN INTERPOLATORY FUNCTION ON A ONE DIMENSIONAL C GRID. THE FUNCTION DETERMINED CAN BE C REPRESENTED BY SPLINES UNDER TENSION. FOR ACTUAL C MAPPING OF POINTS IT IS NECESSARY TO CALL THE SUBROUTINE C CURVBD, WHICH ALSO RETURNS FIRST AND SECOND DERIVATIVES. C C ON INPUT-- C C NX IS THE NUMBER OF GRID POINTS. C (NX SHOULD BE AT LEAST 2) C C X IS ARRAY OF THE NX COORDINATES OF THE GRID POINTS. C THESE SHOULD BE STRICTLY INCREASING. C C W IS AN ARRAY OF THE NX FUNCTIONAL VALUES AT THE C THE GRID POINTS, I. E. W(I) CONTAINS THE FUNCTIONAL C VALUE AT X(I) FOR I = 1,...,NX . C C C IS AN ARRAY OF AT LEAST NX LOCATIONS. THIS C PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS C DESTROYED ON OUTPUT. C C VX IS THE ARRAY OF AT LEAST 5 * NX LOCATIONS. C C TEMP IS AN ARRAY OF AT LEAST 3 * NX LOCATIONS C WHICH IS USED FOR SCRATCH STORAGE. C C SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATE C THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO C (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE C TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE C (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY C BI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF CUBIC C SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS C APPROXIMATELY 1. IN ABSOLUTE VALUE. C C ON OUTPUT-- C C C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE C INTERPOLATED FUNCTION IN A B-SPLINE FORM. C C VX CONTAINS B-SPLINE UNDER TENSION BASIS DATA. C C IERR CONTAINS AN ERROR FLAG. C = 0 FOR NORMAL RETURN, C = 1 IF NX IS LESS THAN 2, C = 2 IF THE X-ARRAY IS NOT STRICTLY C INCREASING. C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF C THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING C SEQUENCE). C C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS, C SNHCSH, TRIDEC, AND TRISOL. C C----------------------------------------------------------- C C COPY W INTO C C DO 1 I = 1,NX 1 C(I) = W(I) C C GENERATE BASIS FUNCTIONS ALONG X-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NX,X,SIGMA,VX,IERR) IF (IERR .NE. 0) RETURN DO 2 I = 2,NX 2 TEMP(I) = VX(5,I-1) NXPI = NX DO 3 I = 1,NX NXPI = NXPI+1 3 TEMP(NXPI) = 1. DO 4 I = 2,NX NXPI = NXPI+1 4 TEMP(NXPI) = VX(4,I) CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1), * TEMP(1),TEMP(NX+1),IERR) CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX, * 1,1) RETURN END C C======================================================================= C SUBROUTINE SURFB1 (NX,NY,X,Y,W,NW1,C,VX,VY,TEMP,SIGMA, * IERR) C INTEGER NX,NY,NW1,IERR REAL X(NX),Y(NY),W(NW1,NY),C(NX,NY),VX(5,NX),VY(5,NY), * TEMP(*),SIGMA C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO C COMPUTE AN INTERPOLATORY FUNCTION ON A TWO DIMENSIONAL C RECTANGULAR GRID. THE FUNCTION DETERMINED CAN BE C REPRESENTED AS A TENSOR PRODUCT OF SPLINES UNDER TENSION C FOR ACTUAL MAPPING OF POINTS IT IS NECESSARY TO CALL THE C SUBROUTINE SURFBD, WHICH ALSO RETURNS FIRST AND SECOND C PARTIAL DERIVATIVES. C C ON INPUT-- C C NX AND NY ARE THE NUMBER OF GRID LINES IN THE X- AND Y C DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR GRID. (NX C AND NY SHOULD BE AT LEAST 2.) C C X AND Y ARE ARRAYS OF THE NX AND NY COORDINATES OF THE C GRID LINES IN X- AND Y-DIRECTIONS, RESPECTIVELY. THESE C SHOULD BE STRICTLY INCREASING. C C W IS AN ARRAY OF THE NX * NY FUNCTIONAL VALUES AT THE C THE GRID POINTS, I. E. W(I,J) CONTAINS THE FUNCTIONAL C VALUE AT (X(I),Y(J)) FOR I = 1,...,NX, AND J = 1,...,NY. C C NW1 IS THE FIRST DIMENSION OF THE ARRAY W USED IN THE C CALLING PROGRAM (NW1 .GE. NX). C C C IS AN ARRAY OF AT LEAST NX * NY LOCATIONS. THIS C PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS C DESTROYED ON OUTPUT. C C VX AND VY ARE ARRAYS OF AT LEAST 5 * NX AND 5 * NY C LOCATIONS, RESPECTIVELY. C C TEMP IS AN ARRAY OF AT LEAST 3 * MAX(NX, NY) LOCATIONS C WHICH IS USED FOR SCRATCH STORAGE. C C AND C C SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATE C THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO C (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE C TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE C (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY C BI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF CUBIC C SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS C APPROXIMATELY 1. IN ABSOLUTE VALUE. C C ON OUTPUT-- C C C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE C INTERPOLATED FUNCTION IN A B-SPLINE TENSOR PRODUCTION C FORM. C C VX AND VY CONTAIN B-SPLINE UNDER TENSION BASIS DATA. C C IERR CONTAINS AN ERROR FLAG. C = 0 FOR NORMAL RETURN, C = 1 IF NX OR NY IS LESS THAN 2, C = 2 IF THE X- OR Y-ARRAYS ARE NOT STRICTLY C INCREASING. C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF C THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING C SEQUENCE). C C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS, C SNHCSH, TRIDEC, AND TRISOL. C C--------------------------------------------------------- - C C COPY W INTO C C DO 1 J = 1,NY DO 1 I = 1,NX 1 C(I,J) = W(I,J) C C GENERATE BASIS FUNCTIONS ALONG X-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NX,X,SIGMA,VX,IERR) IF (IERR .NE. 0) RETURN DO 2 I = 2,NX 2 TEMP(I) = VX(5,I-1) NXPI = NX DO 3 I = 1,NX NXPI = NXPI+1 3 TEMP(NXPI) = 1. DO 4 I = 2,NX NXPI = NXPI+1 4 TEMP(NXPI) = VX(4,I) CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1), * TEMP(1),TEMP(NX+1),IERR) CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX, * NY,1) C C GENERATE BASIS FUNCTIONS ALONG Y-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NY,Y,SIGMA,VY,IERR) IF (IERR .NE. 0) RETURN DO 5 J = 2,NY 5 TEMP(J) = VY(5,J-1) NYPJ = NY DO 6 J = 1,NY NYPJ = NYPJ+1 6 TEMP(NYPJ) = 1. DO 7 J = 2,NY NYPJ = NYPJ+1 7 TEMP(NYPJ) = VY(4,J) CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1), * TEMP(1),TEMP(NY+1),IERR) CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C,1, * NX,NX) RETURN END C C======================================================================= C SUBROUTINE VAL3B1 (NX,NY,NZ,X,Y,Z,W,NW1,NW2,C,VX,VY, * VZ,TEMP,SIGMA,IERR) C INTEGER NX,NY,NZ,NW1,NW2,IERR REAL X(NX),Y(NY),Z(NZ),W(NW1,NW2,NZ),C(NX,NY,NZ), * VX(5,NX),VY(5,NY),VZ(5,NZ),TEMP(*),SIGMA C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO C COMPUTE AN INTERPOLATORY FUNCTION ON A THREE DIMENSIONAL C RECTANGULAR GRID. THE FUNCTION DETERMINED CAN BE C REPRESENTED AS A TENSOR PRODUCT OF SPLINES UNDER TENSION. C FOR ACTUAL MAPPING OF POINTS IT IS NECESSARY TO CALL THE C SUBROUTINE VAL3BD, WHICH ALSO RETURNS FIRST AND SECOND C PARTIAL DERIVATIVES. C C ON INPUT-- C C NX, NY, AND NZ ARE THE NUMBER OF GRID LINES IN THE X-, C Y-, AND Z-DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR C GRID. (NX, NY, AND NZ SHOULD BE AT LEAST 2.) C C X, Y, AND Z ARE ARRAYS OF THE NX, NY, AND NZ COORDINATES C OF THE GRID LINES IN THE X-, Y-, AND Z-DIRECTIONS, C RESPECTIVELY. THESE SHOULD BE STRICTLY INCREASING. C C W IS AN ARRAY OF THE NX * NY * NZ FUNCTIONAL VALUES AT C THE GRID POINTS, I. E. W(I,J,K) CONTAINS THE FUNCTIONAL C VALUE AT (X(I),Y(J),Z(K)) FOR I = 1,...,NX, C J = 1,...,NY, AND K = 1,...,NZ. C C NW1 AND NW2 ARE THE FIRST TWO DIMENSIONS OF THE ARRAY W C USED IN THE CALLING PROGRAM (NW1 .GE. NX AND NW2 .GE. C NY). C C C IS AN ARRAY OF AT LEAST NX * NY * NZ LOCATIONS. THIS C PARAMETER MAY COINCIDE WITH W IN WHICH CASE W IS C DESTROYED ON OUTPUT. C C VX, VY, AND VZ ARE ARRAYS OF AT LEAST 5 * NX, 5 * NY, C AND 5 * NZ LOCATIONS, RESPECTIVELY. C C TEMP IS AN ARRAY OF AT LEAST 3 * MAX(NX, NY, NZ) C LOCATIONS WHICH IS USED FOR SCRATCH STORAGE. C C AND C C SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES C THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO C (E. G. .001) THE RESULTING SURFACE IS APPROXIMATELY THE C TENSOR PRODUCT OF CUBIC SPLINES. IF ABS(SIGMA) IS LARGE C (E. G. 50.) THE RESULTING SURFACE IS APPROXIMATELY C TRI-LINEAR. IF SIGMA EQUALS ZERO TENSOR PRODUCTS OF C CUBIC SPLINES RESULT. A STANDARD VALUE FOR SIGMA IS C APPROXIMATELY 1. IN ABSOLUTE VALUE. C C ON OUTPUT-- C C C CONTAINS THE COEFFICIENTS OF A REPRESENTATION OF THE C INTERPOLATED FUNCTION IN A B-SPLINE TENSOR PRODUCTION C FORM. C C VX, VY, AND VZ CONTAIN B-SPLINE UNDER TENSION BASIS C DATA. C C IERR CONTAINS AN ERROR FLAG. C = 0 FOR NORMAL RETURN, C = 1 IF NX, NY, OR NZ IS LESS THAN 2, C = 2 IF THE X-, Y-, OR Z-ARRAYS ARE NOT STRICTLY C INCREASING. C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT W IF C THIS PARAMETER AND C ARE IDENTICAL IN THE CALLING C SEQUENCE). C C THIS SUBROUTINE REFERENCES PACKAGE MODULES VGEN, TERMS, C SNHCSH, TRIDEC, AND TRISOL. C C----------------------------------------------------------- C C COPY W INTO C C DO 1 K = 1,NZ DO 1 J = 1,NY DO 1 I = 1,NX 1 C(I,J,K) = W(I,J,K) C C GENERATE BASIS FUNCTIONS ALONG X-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NX,X,SIGMA,VX,IERR) IF (IERR .NE. 0) RETURN DO 2 I = 2,NX 2 TEMP(I) = VX(5,I-1) NXPI = NX DO 3 I = 1,NX NXPI = NXPI+1 3 TEMP(NXPI) = 1. DO 4 I = 2,NX NXPI = NXPI+1 4 TEMP(NXPI) = VX(4,I) CALL TRIDEC (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1), * TEMP(1),TEMP(NX+1),IERR) CALL TRISOL (NX,TEMP(1),TEMP(NX+1),TEMP(2*NX+1),C,NX, * NY*NZ,1) C C GENERATE BASIS FUNCTIONS ALONG Y-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NY,Y,SIGMA,VY,IERR) IF (IERR .NE. 0) RETURN DO 5 J = 2,NY 5 TEMP(J) = VY(5,J-1) NYPJ = NY DO 6 J = 1,NY NYPJ = NYPJ+1 6 TEMP(NYPJ) = 1. DO 7 J = 2,NY NYPJ = NYPJ+1 7 TEMP(NYPJ) = VY(4,J) CALL TRIDEC (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1), * TEMP(1),TEMP(NY+1),IERR) DO 8 K = 1,NZ 8 CALL TRISOL (NY,TEMP(1),TEMP(NY+1),TEMP(2*NY+1),C(1,1,K), * 1,NX,NX) C C GENERATE BASIS FUNCTIONS ALONG Z-GRID C SET UP TRIDIAGONAL SYSTEM AND SOLVE C CALL VGEN (NZ,Z,SIGMA,VZ,IERR) IF (IERR .NE. 0) RETURN DO 9 K = 2,NZ 9 TEMP(K) = VZ(5,K-1) NZPK = NZ DO 10 K = 1,NZ NZPK = NZPK+1 10 TEMP(NZPK) = 1. DO 11 K = 2,NZ NZPK = NZPK+1 11 TEMP(NZPK) = VZ(4,K) CALL TRIDEC (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1), * TEMP(1),TEMP(NZ+1),IERR) CALL TRISOL (NZ,TEMP(1),TEMP(NZ+1),TEMP(2*NZ+1),C,1, * NX*NY,NX*NY) RETURN END C C======================================================================= C SUBROUTINE VGEN (N,X,SIGMA,V,IERR) C INTEGER N,IERR REAL X(N),SIGMA,V(5,N) C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE GENERATES AN ARRAY OF COEFFICIENTS USED BY C OTHER SUBROUTINES FOR THE DETERMINATION OF A B-SPLINE C UNDER TENSION BASIS. C C ON INPUT-- C C N IS THE NUMBER OF KNOTS DEFINING THE BASIS (N .GE. 2). C C X IS THE ARRAY OF THE N INCREASING KNOTS. ANY LINEAR C COMBINATION OF THE RESULTING BASIS WILL HAVE THIRD C DERIVATIVE DISCONTINUITIES ONLY AT THE INTERIOR KNOTS, C (I. E. X(2),...,X(N-1) ). C C SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES C THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO C (E. G. .001) THE BASIS FUNCTIONS ARE APPROXIMATELY CUBIC C SPLINES. IF ABS(SIGMA) IS LARGE (E. G. 50.) THE BASIS C FUNCTIONS ARE NEARLY PIECEWISE LINEAR. IF SIGMA EQUALS C ZERO A CUBIC SPLINE BASIS RESULTS. A STANDARD VALUE FOR C SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE. C C AND C C V IS AN ARRAY OF AT LEAST 5*N LOCATIONS. C C ON OUTPUT-- C C V CONTAINS CERTAIN COEFFICIENTS TO BE USED BY OTHER C SUBPROGRAMS FOR THE DETERMINATION OF THE B-SPLINE UNDER C TENSION BASIS. CONSIDERED AS A 5 BY N ARRAY, FOR I = 1, C ... , N, B-SPLINE BASIS FUNCTION I IS SPECIFIED BY-- C V(1,I) = SECOND DERIVATIVE AT X(I-1), FOR I .NE. 1, C V(2,I) = SECOND DERIVATIVE AT X(I), FOR ALL I, C V(3,I) = SECOND DERIVATIVE AT X(I+1), FOR I .NE. N, C V(4,I) = FUNCTION VALUE AT X(I-1), FOR I .NE. 1, C V(5,I) = FUNCTION VALUE AT X(I+1), FOR I .NE. N, C AND THE PROPERTIES THAT IT HAS-- C 1. FUNCTION VALUE 1 AT X(I), C 2. FUNCTION VALUE AND SECOND DERIVATIVE = 0 AT C X(1), ... , X(I-2), AND X(I+2), ... , X(N). C IN V(5,N) AND V(3,N) ARE CONTAINED FUNCTION VALUE AND C SECOND DERIVATIVE OF BASIS FUNCTION ZERO AT X(1), C RESPECTIVELY. IN V(4,1) AND V(1,1) ARE CONTAINED C FUNCTION VALUE AND SECOND DERIVATIVE OF BASIS FUNCTION C N+1 AT X(N), RESPECTIVELY. FUCTION VALUE AND SECOND C DERIVATIVE OF THESE TWO BASIS FUNCTIONS ARE ZERO AT ALL C OTHER KNOTS. ONLY BASIS FUNCTION ZERO HAS NON-ZERO C SECOND DERIVATIVE VALUE AT X(1) AND ONLY BASIS C FUNCTION N+1 HAS NON-ZERO SECOND DERIVATIVE AT X(N). C C IERR CONTAINS AN ERROR FLAG, C = 0 FOR NORMAL RETURN, C = 1 IF N IS LESS THAN 2, C = 2 IF X-VALUES ARE NOT STRICTLY INCREASING. C C AND C C N, X, AND SIGMA ARE UNALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULES TERMS AND C SNHCSH. C C----------------------------------------------------------- C NM1 = N-1 IERR = 0 IF (N .LE. 1) GO TO 3 IF (X(N) .LE. X(1)) GO TO 4 C C DENORMALIZE TENSION FACTOR C SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1)) C C GENERATE COEFFICIENTS FOR LEFT END BASIS FUNCTIONS C D3 = X(2)-X(1) IF (D3 .LE. 0.) GO TO 4 CALL TERMS (DIAG3,SDIAG3,SIGMAP,D3) D4 = D3 IF (N .GE. 3) D4 = X(3)-X(2) IF (D4 .LE. 0.) GO TO 4 CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4) A22 = DIAG3+SDIAG3 A23 = DIAG3+DIAG4+SDIAG3+SDIAG4 V(2,1) = 0. V(3,1) = 1./(D3*(DIAG3+DIAG4)+(D3+D4)*SDIAG4) V(5,1) = SDIAG4*D4*V(3,1) IF (N .EQ. 2) GO TO 2 A22 = 2.*A22 D1 = D3 D2 = D3 D3 = D4 DIAG1 = DIAG3 DIAG2 = DIAG3 DIAG3 = DIAG4 SDIAG1 = SDIAG3 SDIAG2 = SDIAG3 SDIAG3 = SDIAG4 C C GENERATE COEFFICIENTS FOR INTERIOR BASIS FUNCTIONS C DO 1 I = 2,NM1 IF (I .NE. NM1) D4 = X(I+2)-X(I+1) IF (D4 .LE. 0.) GO TO 4 IF (D4 .NE. D3) CALL TERMS (DIAG4,SDIAG4,SIGMAP,D4) A11 = DIAG1+DIAG2+SDIAG1*(1.+D1/D2) A12 = SDIAG2/A11 B1 = 1./(D2*A11) A33 = DIAG3+DIAG4+SDIAG4*(1.+D4/D3) A32 = SDIAG3/A33 B3 = 1./(D3*A33) A21 = A22 A22 = A23 A23 = DIAG3+DIAG4+SDIAG3+SDIAG4 V(2,I) = -(A21*B1+A23*B3)/(A22-A21*A12-A23*A32) V(1,I) = B1-A12*V(2,I) V(3,I) = B3-A32*V(2,I) V(4,I) = SDIAG1*D1*V(1,I) V(5,I) = SDIAG4*D4*V(3,I) C C SAVE CONSTANTS FOR NEXT ITERATION C D1 = D2 D2 = D3 D3 = D4 DIAG1 = DIAG2 DIAG2 = DIAG3 DIAG3 = DIAG4 SDIAG1 = SDIAG2 SDIAG2 = SDIAG3 1 SDIAG3 = SDIAG4 C C GENERATE COEFFICIENTS FOR RIGHT END BASIS FUNCTIONS C V(2,N) = 0. V(1,N) = 1./(D2*(DIAG1+DIAG2)+(D2+D1)*SDIAG1) V(4,N) = SDIAG1*D1*V(1,N) V(3,N) = V(1,3) V(5,N) = V(4,3) V(1,1) = V(3,N-2) V(4,1) = V(5,N-2) C C ADJUST BASIS FOR NATURAL END CONDITIONS C V(4,2) = V(4,2)-V(1,2)*V(5,N)/V(3,N) V(1,2) = 0. V(5,NM1) = V(5,NM1)-V(3,NM1)*V(4,1)/V(1,1) V(3,NM1) = 0. RETURN C C N EQUAL TO 2 C 2 V(4,1) = V(5,1) V(1,1) = V(3,1) V(3,1) = 0. V(5,1) = 0. V(1,2) = 0. V(2,2) = 0. V(3,2) = V(1,1) V(4,2) = 0. V(5,2) = V(4,1) RETURN C C TOO FEW KNOTS C 3 IERR = 1 RETURN C C X-VALUES NOT STRICTLY INCREASING C 4 IERR = 2 RETURN END C C======================================================================= C SUBROUTINE TERMS (DIAG,SDIAG,SIGMA,DEL) C REAL DIAG,SDIAG,SIGMA,DEL C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY A. K. CLINE AND R. J. RENKA C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE COMPUTES THE DIAGONAL AND SUPERDIAGONAL C TERMS OF THE TRIDIAGONAL LINEAR SYSTEM ASSOCIATED WITH C SPLINE UNDER TENSION INTERPOLATION. C C ON INPUT-- C C SIGMA CONTAINS THE TENSION FACTOR. C C AND C C DEL CONTAINS THE STEP SIZE. C C ON OUTPUT-- C C (SIGMA*DEL*COSH(SIGMA*DEL) - SINH(SIGMA*DEL) C DIAG = DEL*--------------------------------------------. C (SIGMA*DEL)**2 * SINH(SIGMA*DEL) C C SINH(SIGMA*DEL) - SIGMA*DEL C SDIAG = DEL*----------------------------------. C (SIGMA*DEL)**2 * SINH(SIGMA*DEL) C C AND C C SIGMA AND DEL ARE UNALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULE SNHCSH. C C----------------------------------------------------------- C IF (SIGMA .NE. 0.) GO TO 1 DIAG = DEL/3. SDIAG = DEL/6. RETURN 1 SIGDEL = SIGMA*DEL CALL SNHCSH (SINHM,COSHM,SIGDEL,0) DENOM = DEL/((SINHM+SIGDEL)*SIGDEL*SIGDEL) DIAG = DENOM*(SIGDEL*COSHM-SINHM) SDIAG = DENOM*SINHM RETURN END C C======================================================================= C SUBROUTINE TRIDEC (N,SUBDI,DIAGI,SUPD,SUBDO,DIAGO, * IERR) C INTEGER N,IERR REAL SUBDI(N),DIAGI(N),SUPD(N),SUBDO(N),DIAGO(N) C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE FACTORIZES A TRIDIAGONAL MATRIX IN ORDER C TO SOLVE SYSTEMS OF LINEAR EQUATIONS. THE FACTORIZATION C EMPLOYS GAUSSIAN ELIMINATION WITHOUT ANY INTERCHANGING OF C COLUMNS OR ROWS. THE SUBROUTINE TRISOL MAY BE CALLED TO C ACTUALLY SOLVE THE SYSTEM ONCE THE FACTORIZATION HAS BEEN C PERFORMED. C C ON INPUT-- C C N CONTAINS THE ORDER OF THE MATRIX (N .GE. 1). C C SUBDI IS AN ARRAY CONTAINING THE SUBDIAGONAL ELEMENTS OF C THE MATRIX IN POSITIONS 2, ... , N. C C DIAGI IS AN ARRAY CONTAINING THE DIAGONAL ELEMENTS OF C THE MATRIX. C C SUPD IS AN ARRAY CONTAINING THE SUPERDIAGONAL ELEMENTS C OF THE MATRIX IN POSITIONS 1, ... , N-1. C C AND C C SUBDO AND DIAGO ARE ARRAYS OF LENGTH N. (THE STORAGE C FOR THESE MAY COINCIDE WITH THAT FOR SUBDI AND DIAGI, C RESPECTIVELY, IN WHICH CASE THE ORIGINAL CONTENTS OF C SUBDI AND DIAGI WILL BE DESTROYED.) C C ON OUTPUT-- C C SUBDO AND DIAGO CONTAIN THE SUBDIAGONAL AND DIAGONAL OF C THE FACTORIZATION MATRIX. C C IERR CONTAINS AN ERROR FLAG, C = 0 FOR NORMAL RETURN, C = 1 IF N IS LESS THAN 1, C = 2 IF THE SYSTEM IS SINGULAR. C C AND C C N, SUBDI, DIAGI, AND SUPD ARE UNALTERED (UNLESS STORAGE C FOR SUBDI OR DIAGI COINCIDED WITH THAT FOR SUBDO C OR DIAGO, RESPECTIVELY). C C----------------------------------------------------------- C IF (N .LE. 0) GO TO 3 IERR = 2 DIAGO(1) = DIAGI(1) IF (N .EQ. 1) GO TO 2 C C FORWARD ELIMINATION C DO 1 I = 2,N IM1 = I-1 IF (DIAGO(IM1) .EQ. 0.) RETURN DIAGO(IM1) = 1./DIAGO(IM1) SUBDO(I) = SUBDI(I)*DIAGO(IM1) 1 DIAGO(I) = DIAGI(I)-SUBDO(I)*SUPD(IM1) 2 IF (DIAGO(N) .EQ. 0.) RETURN DIAGO(N) = 1./DIAGO(N) IERR = 0 RETURN C C N LESS THAN 1 C 3 IERR = 1 RETURN END C C======================================================================= C SUBROUTINE TRISOL (N,SUBD,DIAG,SUPD,RHS,MRHS,NUMRHS, * INCRHS) C INTEGER N,MRHS,NUMRHS,INCRHS REAL SUBD(N),DIAG(N),SUPD(N) REAL RHS(1+INCRHS*(N-1)+MRHS*(NUMRHS-1)) C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C REVISED -- DECEMBER 31, 1992 C BY LUDEK KLIMES C INSTITUTE OF GEOTECHNICS C CZECHOSL. ACAD. SCI., PRAGUE C C THIS SUBROUTINE SOLVES TRIDIAGONAL SYSTEMS OF LINEAR C EQUATIONS WITH MULTIPLE RIGHT HAND SIDES. THE RIGHT HAND C SIDES MAY BE STORED ROW-WISE OR COLUMN-WISE. THE C SUBROUTINE TRIDEC SHOULD BE CALLED EARLIER TO DETERMINE A C FACTORIZATION OF THE TRIDIAGONAL MATRIX. THE SOLUTION C VECTORS OVER-WRITE THE RIGHT HAND SIDES. C C ON INPUT-- C C N CONTAINS THE ORDER OF THE MATRIX (N .GE. 1). C C SUBD, DIAG, AND SUPD ARE ARRAYS OF LENGTH N CONTAINING C THE SUBDIAGONAL, DIAGONAL, AND SUPERDIAGONAL OF THE C FACTORIZATION, RESPECTIVELY. C C RHS IS AN ARRAY CONTAINING THE RIGHT HAND SIDES OF THE C TRIDIAGONAL SYSTEM. C C MRHS IS THE INCREMENT BETWEEN THE FIRST COMPONENTS OF C EACH OF THE RIGHT HAND SIDE VECTORS IN STORAGE (MRHS C .GE. 1). C C NUMRHS IS THE NUMBER OF RIGHT HAND SIDES TO BE SOLVED. C C AND C C INCRHS IS THE INCREMENT BETWEEN COMPONENTS WITHIN EACH C OF THE RIGHT HAND SIDE VECTORS IN STORAGE (INCRHS .GE. C 1). C C THE PARAMETERS N, SUBD, DIAG, AND SUPD MAY BE INPUT AS THE C PARAMETERS N, SUBDO, DIAGO, AND SUPD OUTPUT BY SUBROUTINE C TRIDEC, RESPECTIVELY. C C ON OUTPUT-- C C RHS CONTAINS THE SOLUTION VECTORS IN THE SAME STORAGE C STRUCTURE AS FOR THE RIGHT HAND SIDES. C C AND C C N, SUBD, DIAG, SUPD, MRHS, NUMRHS, AND INCRHS ARE C UNALTERED. C C----------------------------------------------------------- C NP1 = N+1 C C LOOP ON RIGHT HAND SIDES C DO 4 K = 1,NUMRHS C C FORWARD ELIMINATION C IRHS = 1+MRHS*(K-1) IF (N .EQ. 1) GO TO 2 DO 1 I = 2,N IM1RHS = IRHS IRHS = IRHS+INCRHS 1 RHS(IRHS) = RHS(IRHS)-SUBD(I)*RHS(IM1RHS) C C BACK SUBSTITUTION C 2 RHS(IRHS) = DIAG(N)*RHS(IRHS) IF (N .EQ. 1) GO TO 4 DO 3 IBAK = 2,N I = NP1-IBAK RHS(IM1RHS) = DIAG(I)*(RHS(IM1RHS)-SUPD(I) * *RHS(IRHS)) IRHS = IM1RHS 3 IM1RHS = IM1RHS-INCRHS 4 CONTINUE RETURN END C C======================================================================= C C PART 2: C C======================================================================= C SUBROUTINE CURVBD (XX,W,WX,WXX,NX,X,C,VX,SIGMA) C INTEGER NX REAL XX,W,WX,WXX,X(NX),VX(5,NX),C(NX),SIGMA C C COMPLEMENT TO FITPACK C BY ALAN KAYLOR CLINE C CODED -- OCTOBER 9, 1986 C BY LUDEK KLIMES C INST. GEOL. GEOTECHN. C CZECHOSL. ACAD. SCI., PRAGUE C C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE C FIRST PARTIAL DERIVATIVE, AND THE SECOND PARTIAL C DERIVATIVE OF A SPLINE UNDER TENSION IN ONE VARIABLE. C C ON INPUT-- C C XX CONTAINS THE X-COORDINATE OF THE POINT C AT WHICH THE INTERPOLATION IS TO BE PERFORMED C C NX IS THE NUMBER OF GRID POINTS C C X IS ARRAY CONTAINING THE X-GRID VALUES. C C C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN C TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE C EXPANSION OF THE FUNCTION, FOR I = 1,...,NX , C THE COEFFICIENT MULTIPLYING THE BASIS C FUNCTION I IS STORED IN C(I). C C VX IS THE ARRAY OF LENGTH 5*NX C CONTAINING THE B-SPLINE BASIS DATA C C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED). C C THE PARAMETERS NX, X, C, VX, AND SIGMA C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF CURVB1. C C ON OUTPUT-- C C W CONTAINS THE INTERPOLATED FUNCTION VALUE. C C WX CONTAINS THE FIRST DERIVATIVE . C C WXX CONTAINS THE SECOND DERIVATIVE . C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL, C AND SNHCSH. C C-------------------------------------------------------------- C REAL BX(3,4) C C EVALUATE BASIS FUNCTIONS AT XX C CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX) C C ACCUMULATE BASIS FUNCTIONS C SUM = 0. SUMX = 0. SUMXX = 0. DO 1 I = 1,4 II = ISTART+I-1 IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1 BX1I = BX(1,I) CI = C(II) SUM = SUM+CI*BX1I SUMX = SUMX+CI*BX(2,I) SUMXX = SUMXX+CI*BX(3,I) * CALL VAR2(II,BX1I,BX(2,I),0.,0.) 1 CONTINUE W = SUM WX = SUMX WXX = SUMXX RETURN END C C======================================================================= C SUBROUTINE SURFBD (XX,YY,W,WX,WY,WXX,WXY,WYY,NX,NY,X, * Y,C,VX,VY,SIGMA) C INTEGER NX,NY REAL XX,YY,W,WX,WY,WXX,WXY,WYY,X(NX),Y(NY),VX(5,NX), * VY(5,NY),C(NX,NY),SIGMA C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE TWO C FIRST PARTIAL DERIVATIVES, AND THE SIX SECOND PARTIAL C DERIVATIVES OF A TENSOR PRODUCT SPLINE UNDER TENSION IN C TWO VARIABLES. C C ON INPUT-- C C XX AND YY CONTAIN THE X- AND Y-COORDINES OF THE POINT C AT WHICH THE INTERPOLATION IS TO BE PERFORME .D. C C NX AND NY ARE THE NUMBER OF GRID LINES IN THE X- AND Y- C DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR GRID WHICH C SPECIFIED THE FUNCTION. C C X AND Y ARE ARRAYS CONTAINING THE X- AND Y-GRID VALUES, C RESPECTIVELY. C C C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN C TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE C EXPANSION OF THE FUNCTION, FOR I = 1,...,NX AND J = 1, C ...,NY, THE COEFFICIENT MULTIPLYING THE PRODUCT OF BASIS C FUNCTION I IN X AND BASIS FUNCTION J IN Y IS STORED IN C C(I,J). C C VX AND VY VZ ARE ARRAYS OF LENGTH 5*NX AND 5*NY, C RESPECTIVELY, CONTAINING THE B-SPLINE BASIS DATA FOR THE C X- AND Y-GRIDS. C C AND C C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED). C C THE PARAMETERS NX, NY, X, Y, Z, C, VX, VY, AND SIGMA C SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF SURFB1. C C ON OUTPUT-- C C W CONTAINS THE INTERPOLATED FUNCTION VALUE. C C WX AND WY CONTAIN THE X- AND Y-PARTIAL DERIVATIVES, C RESPECTIVELY. C C WXX, WXY, AND WYY CONTAIN THE XX-, XY-, AND YY-PARTIAL C DERIVATIVES, RESPECTIVELY. C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL, C AND SNHCSH. C C--------------------------------------------------------- ---- C REAL BX(3,4),BY(3,4) C C EVALUATE BASIS FUNCTIONS AT XX AND YY C CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX) CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY) C C ACCUMULATE TENSOR PRODUCTS C SUM = 0. SUMX = 0. SUMY = 0. SUMXX = 0. SUMXY = 0. SUMYY = 0. DO 2 J = 1,4 JJ = JSTART+J-1 IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2 BY1J = BY(1,J) BY2J = BY(2,J) BY3J = BY(3,J) DO 1 I = 1,4 II = ISTART+I-1 IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1 BX1I = BX(1,I) BX2I = BX(2,I) CIJ = C(II,JJ) SUM = SUM+CIJ*BX1I*BY1J SUMX = SUMX+CIJ*BX2I*BY1J SUMY = SUMY+CIJ*BX1I*BY2J SUMXX = SUMXX+CIJ*BX(3,I)*BY1J SUMXY = SUMXY+CIJ*BX2I*BY2J SUMYY = SUMYY+CIJ*BX1I*BY3J * CALL VAR2(II+NX*(JJ-1),BX1I*BY1J,BX2I*BY1J,BX1I*BY2J,0.) 1 CONTINUE 2 CONTINUE W = SUM WX = SUMX WY = SUMY WXX = SUMXX WXY = SUMXY WYY = SUMYY RETURN END C C======================================================================= C SUBROUTINE VAL3BD (XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY, * WYZ,WZZ,WXZ,NX,NY,NZ,X,Y,Z,C,VX,VY, * VZ,SIGMA) C INTEGER NX,NY,NZ REAL XX,YY,ZZ,W,WX,WY,WZ,WXX,WXY,WYY,WYZ,WZZ,WXZ, * X(NX),Y(NY),Z(NZ),VX(5,NX),VY(5,NY),VZ(5,NZ), * C(NX,NY,NZ),SIGMA C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE EVALUATES THE FUNCTION VALUE, THE THREE C FIRST PARTIAL DERIVATIVES, AND THE SIX SECOND PARTIAL C DERIVATIVES OF A TENSOR PRODUCT SPLINE UNDER TENSION IN C THREE VARIABLES. C C ON INPUT-- C C XX, YY, AND ZZ CONTAIN THE X-, Y-, AND Z-COORDINATES OF C THE POINT AT WHICH THE INTERPOLATION IS TO BE PERFORMED. C C NX, NY, AND NZ ARE THE NUMBER OF GRID LINES IN THE X-, C Y-, AND Z-DIRECTIONS, RESPECTIVELY, OF THE RECTANGULAR C GRID WHICH SPECIFIED THE FUNCTION. C C X, Y, AND Z ARE ARRAYS CONTAINING THE X-, Y-, AND Z-GRID C VALUES, RESPECTIVELY. C C C IS AN ARRAY OF COEFFICIENTS DESCRIBING THE FUNCTION IN C TERMS OF A B-SPLINE UNDER TENSION BASIS. IN THE C EXPANSION OF THE FUNCTION, FOR I = 1,...,NX, J = 1,..., C NY, AND K = 1,...,NZ, THE COEFFICIENT MULTIPLYING THE C PRODUCT OF BASIS FUNCTION I IN X, BASIS FUNCTION J IN Y, C AND BASIS FUNCTION K IN Z IS STORED IN C(I,J,K). C C VX, VY, AND VZ ARE ARRAYS OF LENGTH 5*NX, 5*NY, AND C 5*NZ, RESPECTIVELY, CONTAINING THE B-SPLINE BASIS DATA C FOR THE X-, Y-, AND Z-GRIDS. C C AND C C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED). C C THE PARAMETERS NX, NY, NZ, X, Y, Z, C, VX, VY, VZ, AND C SIGMA SHOULD BE INPUT UNALTERED FROM THE OUTPUT OF C VAL3B1. C C ON OUTPUT-- C C W CONTAINS THE INTERPOLATED FUNCTION VALUE. C C WX, WY, AND WZ CONTAIN THE X-, Y-, AND Z-PARTIAL C DERIVATIVES, RESPECTIVELY. C C WXX, WXY, WYY, WYZ, WZZ, AND WXZ CONTAIN THE XX-, XY- C YY-, YZ-, ZZ-, AND XZ-PARTIAL DERIVATIVES, RESPECTIVELY. C C AND C C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULES DSPLNZ, INTRVL, C AND SNHCSH. C C-------------------------------------------------------------- C REAL BX(3,4),BY(3,4),BZ(3,4) C C EVALUATE BASIS FUNCTIONS AT XX, YY, AND ZZ C CALL DSPLNZ (XX,NX,X,VX,SIGMA,ISTART,BX) CALL DSPLNZ (YY,NY,Y,VY,SIGMA,JSTART,BY) CALL DSPLNZ (ZZ,NZ,Z,VZ,SIGMA,KSTART,BZ) C C ACCUMULATE TENSOR PRODUCTS C SUM = 0. SUMX = 0. SUMY = 0. SUMZ = 0. SUMXX = 0. SUMXY = 0. SUMYY = 0. SUMYZ = 0. SUMZZ = 0. SUMXZ = 0. DO 3 K = 1,4 KK = KSTART+K-1 IF (KK .EQ. 0 .OR. KK .GT. NZ) GO TO 3 BZ1K = BZ(1,K) BZ2K = BZ(2,K) BZ3K = BZ(3,K) DO 2 J = 1,4 JJ = JSTART+J-1 IF (JJ .EQ. 0 .OR. JJ .GT. NY) GO TO 2 BY1J = BY(1,J) BY2J = BY(2,J) BY3J = BY(3,J) DO 1 I = 1,4 II = ISTART+I-1 IF (II .EQ. 0 .OR. II .GT. NX) GO TO 1 BX1I = BX(1,I) BX2I = BX(2,I) CIJK = C(II,JJ,KK) SUM = SUM+CIJK*BX1I*BY1J*BZ1K SUMX = SUMX+CIJK*BX2I*BY1J*BZ1K SUMY = SUMY+CIJK*BX1I*BY2J*BZ1K SUMZ = SUMZ+CIJK*BX1I*BY1J*BZ2K SUMXX = SUMXX+CIJK*BX(3,I)*BY1J*BZ1K SUMXY = SUMXY+CIJK*BX2I*BY2J*BZ1K SUMYY = SUMYY+CIJK*BX1I*BY3J*BZ1K SUMYZ = SUMYZ+CIJK*BX1I*BY2J*BZ2K SUMZZ = SUMZZ+CIJK*BX1I*BY1J*BZ3K SUMXZ = SUMXZ+CIJK*BX2I*BY1J*BZ2K * CALL VAR2(II+NX*(JJ-1+NY*(KK-1)),BX1I*BY1J*BZ1K, * * BX2I*BY1J*BZ1K,BX1I*BY2J*BZ1K,BX1I*BY1J*BZ2K) 1 CONTINUE 2 CONTINUE 3 CONTINUE W = SUM WX = SUMX WY = SUMY WZ = SUMZ WXX = SUMXX WXY = SUMXY WYY = SUMYY WYZ = SUMYZ WZZ = SUMZZ WXZ = SUMXZ RETURN END C C======================================================================= C SUBROUTINE DSPLNZ (T,N,X,V,SIGMA,ISTART,B) C INTEGER N,ISTART REAL T,X(N),V(5,N),SIGMA,B(3,4) C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY ALAN KAYLOR CLINE C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS SUBROUTINE EVALUATES AT A GIVEN POINT THE FOUR NON- C ZERO BASIS FUNCTIONS OF A B-SPLINE UNDER TENSION BASIS AND C THEIR FIRST AND SECOND DERIVATIVES. THE INDEX OF THE FIRST C NON-ZERO BASIS FUNCTION IS ALSO DETERMINED. (THE SENSE OF C THE WORD NON-ZERO IS EXTENDED TO INCLUDE THE SPECIAL CASE C WHERE THE GIVEN POINT COINCIDES WITH A KNOT IN WHICH CASE C THE LAST OF THE FOUR RETURNED FUNCTION VALUES MAY BE ZERO. C ) THE SUBROUTINE VGEN SHOULD BE CALLED EARLIER TO C DETERMINE CERTAIN NECESSARY COEFFICIENTS. C C ON INPUT-- C C T CONTAINS A REAL VALUE AT WHICH THE BASIS FUNCTIONS ARE C TO BE EVALUATED. C C N CONTAINS THE NUMBER OF KNOTS DEFINING THE BASIS. C C X CONTAINS THE ARRAY OF KNOTS. C C V CONTAINS THE ARRAY OF COEFFICIENTS DETERMINED BY VGEN C FOR CALCULATION OF BASIS FUNCTIONS. C C SIGMA CONTAINS THE TENSION FACTOR (ITS SIGN IS IGNORED). C C ISTART IS AN INTEGER VARIABLE. C C AND C C B IS A REAL ARRAY WITH 3 ROWS AND 4 COLUMNS. C C THE PARAMETERS N, X, V, AND SIGMA SHOULD BE INPUT C UNALTERED FROM THE OUTPUT OF VGEN. C C ON OUTPUT-- C C ISTART CONTAINS THE INDEX OF THE FIRST NON-ZERO BASIS C FUNCTION. THUS 0 .LE. ISTART .LE. N-2 AND THE N0N-ZERO C BASIS FUNCTIONS HAVE INDICES ISTART, ... , ISTART+3. C C B CONTAINS THE VALUES AT T OF BASIS FUNCTIONS ISTART, C ... , ISTART+3 IN B(1,1), ... , B(1,4), RESPECTIVELY. C FIRST AND SECOND DERIVATIVES OF THE CORRESPONDING C FUNCTIONS ARE CONTAINED IN B(2,1), ... , B(2,4), AND C B(3,1), ... , B(3,4), RESPECTIVELY. C C T, N, X, V, AND SIGMA ARE UNALTERED. C C THIS SUBROUTINE REFERENCES PACKAGE MODULES INTRVL AND C SNHCSH. C C----------------------------------------------------------- C C DENORMALIZE TENSION FACTOR C SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1)) C C DETERMINE INDEX OF FIRST NON-ZERO BASIS FUNCTION C I = INTRVL (T,X,N)-1 C C COMPUTE DISTANCES TO ADJACENT KNOTS AND LAGRANGIAN C WEIGHTS C DEL1 = T-X(I+1) DEL2 = X(I+2)-T DELS = X(I+2)-X(I+1) C10 = DEL2/DELS C20 = DEL1/DELS C11 = -1./DELS C21 = 1./DELS IF (SIGMAP .NE. 0.) GO TO 1 FAC = -DEL1*DEL2/(6.*DELS) CP10 = FAC*(DEL2+DELS) CP20 = FAC*(DEL1+DELS) CP11 = -(2.*DEL2*DEL2-DEL1*(DEL2+DELS))/(6.*DELS) CP21 = (2.*DEL1*DEL1-DEL2*(DEL1+DELS))/(6.*DELS) CP12 = C10 CP22 = C20 GO TO 2 1 DELP1 = SIGMAP*(DEL1+DELS)/2. DELP2 = SIGMAP*(DEL2+DELS)/2. CALL SNHCSH (SINHM1,COSHM1,SIGMAP*DEL1,0) CALL SNHCSH (SINHM2,COSHM2,SIGMAP*DEL2,0) CALL SNHCSH (SINHMS,DUMMY,SIGMAP*DELS,-1) CALL SNHCSH (SINHP1,DUMMY,SIGMAP*DEL1/2.,-1) CALL SNHCSH (SINHP2,DUMMY,SIGMAP*DEL2/2.,-1) CALL SNHCSH (DUMMY,COSHP1,DELP1,1) CALL SNHCSH (DUMMY,COSHP2,DELP2,1) SINHS = SINHMS+SIGMAP*DELS DENOM = SIGMAP*SIGMAP*DELS*SINHS CP10 = (SINHM2*DEL1-DEL2*(2.*(COSHP2+1.)*SINHP1 * +SIGMAP*COSHP2*DEL1))/DENOM CP20 = (SINHM1*DEL2-DEL1*(2.*(COSHP1+1.)*SINHP2 * +SIGMAP*COSHP1*DEL2))/DENOM CP11 = -(DELS*SIGMAP*COSHM2-SINHMS)/DENOM CP21 = (DELS*SIGMAP*COSHM1-SINHMS)/DENOM CP12 = (SINHM2+SIGMAP*DEL2)/SINHS CP22 = (SINHM1+SIGMAP*DEL1)/SINHS C C COMPUTE BASIS FUNCTION VALUES C 2 II = I IF (II .EQ. 0) II = N IIP1 = I+1 IIP2 = I+2 IIP3 = I+3 IF (IIP2 .EQ. N) IIP3 = 1 B(1,1) = C10*V(5,II)+CP10*V(3,II) B(1,2) = C10+C20*V(5,IIP1)+CP10*V(2,IIP1)+ * CP20*V(3,IIP1) B(1,3) = C10*V(4,IIP2)+C20+CP10*V(1,IIP2)+ * CP20*V(2,IIP2) B(1,4) = C20*V(4,IIP3)+CP20*V(1,IIP3) B(2,1) = C11*V(5,II)+CP11*V(3,II) B(2,2) = C11+C21*V(5,IIP1)+CP11*V(2,IIP1)+ * CP21*V(3,IIP1) B(2,3) = C11*V(4,IIP2)+C21+CP11*V(1,IIP2)+ * CP21*V(2,IIP2) B(2,4) = C21*V(4,IIP3)+CP21*V(1,IIP3) B(3,1) = CP12*V(3,II) B(3,2) = CP12*V(2,IIP1)+CP22*V(3,IIP1) B(3,3) = CP12*V(1,IIP2)+CP22*V(2,IIP2) B(3,4) = CP22*V(1,IIP3) ISTART = I RETURN END C C======================================================================= C FUNCTION INTRVL (T,X,N) C INTEGER N REAL T,X(N) C C FROM FITPACK -- AUGUST 31, 1981 C CODED BY A. K. CLINE AND R. J. RENKA C DEPARTMENT OF COMPUTER SCIENCES C UNIVERSITY OF TEXAS AT AUSTIN C C THIS FUNCTION DETERMINES THE INDEX OF THE INTERVAL C (DETERMINED BY A GIVEN INCREASING SEQUENCE) IN WHICH C A GIVEN VALUE LIES. C C ON INPUT-- C C T IS THE GIVEN VALUE. C C X IS A VECTOR OF STRICTLY INCREASING VALUES. C C AND C C N IS THE LENGTH OF X (N .GE. 2). C C ON OUTPUT-- C C INTRVL RETURNS AN INTEGER I SUCH THAT C C I = 1 IF T .LT. X(2) , C I = N-1 IF X(N-1) .LE. T , C OTHERWISE X(I) .LE. T .LT. X(I+1), C C NONE OF THE INPUT PARAMETERS ARE ALTERED. C C----------------------------------------------------------- C TT = T IF (TT .LT. X(2)) GO TO 4 IF (TT .GE. X(N-1)) GO TO 5 IL = 2 IH = N-1 C C LINEAR INTERPOLATION C 1 I = MIN0(IL+IFIX(FLOAT(IH-IL)*(TT-X(IL))/(X(IH)-X(IL))), * IH-1) IF (TT .LT. X(I)) GO TO 2 IF (TT .LT. X(I+1)) GO TO 3 C C TOO HIGH C IL = I+1 GO TO 1 C C TOO LOW C 2 IH = I GO TO 1 3 INTRVL = I RETURN C C LEFT END C 4 INTRVL = 1 RETURN C C RIGHT END C 5 INTRVL = N-1 RETURN END