C PROGRAM TTIDBS(OUTPUT,TAPE6=OUTPUT) ID000070 C THIS PROGRAM IS A TEST PROGRAM FOR THE IDBVIP/IDSFFT SUBPRO- ID000080 C GRAM PACKAGE. ALL ELEMENTS OF RESULTING DZI1 AND DZI2 ARRAYS ID000090 C ARE EXPECTED TO BE ZERO. ID000100 C THE LUN CONSTANT IN THE LAST DATA INITIALIZATION STATEMENT IS ID000110 C THE LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, ID000120 C THEREFORE, SYSTEM DEPENDENT. ID000130 C DECLARATION STATEMENTS ID000140 DIMENSION XD(30),YD(30),ZD(30), ID000150 1 XI(6),YI(5),ZI(6,5), ID000160 2 ZI1(6,5),ZI2(6,5),DZI1(6,5),DZI2(6,5), ID000170 3 IWK(1030),WK(240) ID000180 DATA NCP/4/ ID000190 DATA NDP/30/ ID000200 DATA XD(1), XD(2), XD(3), XD(4), XD(5), XD(6), ID000210 1 XD(7), XD(8), XD(9), XD(10),XD(11),XD(12), ID000220 2 XD(13),XD(14),XD(15),XD(16),XD(17),XD(18), ID000230 3 XD(19),XD(20),XD(21),XD(22),XD(23),XD(24), ID000240 4 XD(25),XD(26),XD(27),XD(28),XD(29),XD(30)/ ID000250 5 11.16, 24.20, 19.85, 10.35, 19.72, 0.00, ID000260 6 20.87, 19.99, 10.28, 4.51, 0.00, 16.70, ID000270 7 6.08, 25.00, 14.90, 0.00, 9.66, 5.22, ID000280 8 11.77, 15.10, 25.00, 25.00, 14.59, 15.20, ID000290 9 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/ ID000300 DATA YD(1), YD(2), YD(3), YD(4), YD(5), YD(6), ID000310 1 YD(7), YD(8), YD(9), YD(10),YD(11),YD(12), ID000320 2 YD(13),YD(14),YD(15),YD(16),YD(17),YD(18), ID000330 3 YD(19),YD(20),YD(21),YD(22),YD(23),YD(24), ID000340 4 YD(25),YD(26),YD(27),YD(28),YD(29),YD(30)/ ID000350 5 1.24, 16.23, 10.72, 4.11, 1.39, 20.00, ID000360 6 20.00, 4.62, 15.16, 20.00, 4.48, 19.65, ID000370 7 4.58, 11.87, 3.12, 0.00, 20.00, 14.66, ID000380 8 10.47, 17.19, 3.87, 0.00, 8.71, 0.00, ID000390 9 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/ ID000400 DATA ZD(1), ZD(2), ZD(3), ZD(4), ZD(5), ZD(6), ID000410 1 ZD(7), ZD(8), ZD(9), ZD(10),ZD(11),ZD(12), ID000420 2 ZD(13),ZD(14),ZD(15),ZD(16),ZD(17),ZD(18), ID000430 3 ZD(19),ZD(20),ZD(21),ZD(22),ZD(23),ZD(24), ID000440 4 ZD(25),ZD(26),ZD(27),ZD(28),ZD(29),ZD(30)/ ID000450 5 22.15, 2.83, 7.97, 22.33, 16.83, 34.60, ID000460 6 5.74, 14.72, 21.59, 15.61, 61.77, 6.31, ID000470 7 35.74, 4.40, 21.70, 58.20, 4.73, 40.36, ID000480 8 13.62, 12.57, 8.74, 12.00, 14.81, 21.60, ID000490 9 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/ ID000500 DATA NXI/6/, NYI/5/ ID000510 DATA XI(1), XI(2), XI(3), XI(4), XI(5), XI(6)/ ID000520 1 0.00, 5.00, 10.00, 15.00, 20.00, 25.00/ ID000530 DATA YI(1), YI(2), YI(3), YI(4), YI(5)/ ID000540 1 0.00, 5.00, 10.00, 15.00, 20.00/ ID000550 DATA ZI(1,1),ZI(2,1),ZI(3,1),ZI(4,1),ZI(5,1),ZI(6,1), ID000560 1 ZI(1,2),ZI(2,2),ZI(3,2),ZI(4,2),ZI(5,2),ZI(6,2), ID000570 2 ZI(1,3),ZI(2,3),ZI(3,3),ZI(4,3),ZI(5,3),ZI(6,3), ID000580 3 ZI(1,4),ZI(2,4),ZI(3,4),ZI(4,4),ZI(5,4),ZI(6,4), ID000590 4 ZI(1,5),ZI(2,5),ZI(3,5),ZI(4,5),ZI(5,5),ZI(6,5)/ ID000600 5 58.20, 39.55, 26.90, 21.71, 17.68, 12.00, ID000610 6 61.58, 39.39, 22.04, 21.29, 14.36, 8.04, ID000620 7 59.18, 27.39, 16.78, 13.25, 8.59, 5.36, ID000630 8 52.82, 40.27, 22.76, 16.61, 7.40, 2.88, ID000640 9 34.60, 14.05, 4.12, 3.17, 6.31, 0.60/ ID000650 DATA LUN/6/ ID000660 C CALCULATION ID000670 10 MD=1 ID000680 DO 12 IYI=1,NYI ID000690 DO 11 IXI=1,NXI ID000700 IF(IXI.NE.1.OR.IYI.NE.1) MD=2 ID000710 CALL IDBVIP(MD,NCP,NDP,XD,YD,ZD,1,XI(IXI),YI(IYI), ID000720 1 ZI1(IXI,IYI),IWK,WK) ID000730 11 CONTINUE ID000740 12 CONTINUE ID000750 15 CALL IDSFFT(1,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI2,IWK,WK) ID000760 DO 17 IYI=1,NYI ID000770 DO 16 IXI=1,NXI ID000780 DZI1(IXI,IYI)=ABS(ZI1(IXI,IYI)-ZI(IXI,IYI)) ID000790 DZI2(IXI,IYI)=ABS(ZI2(IXI,IYI)-ZI(IXI,IYI)) ID000800 16 CONTINUE ID000810 17 CONTINUE ID000820 C PRINTING OF INPUT DATA ID000830 20 WRITE (LUN,2020) NDP ID000840 DO 23 IDP=1,NDP ID000850 IF(MOD(IDP,5).EQ.1) WRITE (LUN,2021) ID000860 WRITE (LUN,2022) IDP,XD(IDP),YD(IDP),ZD(IDP) ID000870 23 CONTINUE ID000880 C PRINTING OF OUTPUT RESULTS ID000890 30 WRITE (LUN,2030) ID000900 WRITE (LUN,2031) YI ID000910 DO 33 IXI=1,NXI ID000920 WRITE (LUN,2032) XI(IXI),(ZI1(IXI,IYI),IYI=1,NYI) ID000930 33 CONTINUE ID000940 40 WRITE (LUN,2040) ID000950 WRITE (LUN,2031) YI ID000960 DO 43 IXI=1,NXI ID000970 WRITE (LUN,2032) XI(IXI),(DZI1(IXI,IYI),IYI=1,NYI) ID000980 43 CONTINUE ID000990 50 WRITE (LUN,2050) ID001000 WRITE (LUN,2031) YI ID001010 DO 53 IXI=1,NXI ID001020 WRITE (LUN,2032) XI(IXI),(ZI2(IXI,IYI),IYI=1,NYI) ID001030 53 CONTINUE ID001040 60 WRITE (LUN,2060) ID001050 WRITE (LUN,2031) YI ID001060 DO 63 IXI=1,NXI ID001070 WRITE (LUN,2032) XI(IXI),(DZI2(IXI,IYI),IYI=1,NYI) ID001080 63 CONTINUE ID001090 STOP ID001100 C FORMAT STATEMENTS ID001110 2020 FORMAT(1H1,6HTTIDBS/////3X,10HINPUT DATA,8X,5HNDP =,I3/// ID001120 1 30H I XD YD ZD /) ID001130 2021 FORMAT(1X) ID001140 2022 FORMAT(5X,I2,2X,3F7.2) ID001150 2030 FORMAT(1H1,6HTTIDBS/////3X,17HIDBVIP SUBROUTINE/// ID001160 1 26X,10HZI1(XI,YI)) ID001170 2031 FORMAT(7X,2HXI,4X,3HYI=/12X,5F7.2/) ID001180 2032 FORMAT(1X/1X,F9.2,2X,5F7.2) ID001190 2040 FORMAT(1X/////3X,10HDIFFERENCE/// ID001200 1 25X,11HDZI1(XI,YI)) ID001210 2050 FORMAT(1H1,6HTTIDBS/////3X,17HIDSFFT SUBROUTINE/// ID001220 1 26X,10HZI2(XI,YI)) ID001230 2060 FORMAT(1X/////3X,10HDIFFERENCE/// ID001240 1 25X,11HDZI2(XI,YI)) ID001250 END ID001260 SUBROUTINE IDBVIP(MD,NCP,NDP,XD,YD,ZD,NIP,XI,YI,ZI, ID001340 1 IWK,WK) C THIS SUBROUTINE PERFORMS BIVARIATE INTERPOLATION WHEN THE PRO- C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY C DISTRIBUTED IN THE PLANE. C THE INPUT PARAMETERS ARE C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3), C = 1 FOR NEW NCP AND/OR NEW XD-YD, C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI, C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI, C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI- C MATING PARTIAL DERIVATIVES AT EACH DATA POINT C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP), C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER), C XD = ARRAY OF DIMENSION NDP CONTAINING THE X C COORDINATES OF THE DATA POINTS, C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y C COORDINATES OF THE DATA POINTS, C ZD = ARRAY OF DIMENSION NDP CONTAINING THE Z C COORDINATES OF THE DATA POINTS, C NIP = NUMBER OF OUTPUT POINTS AT WHICH INTERPOLATION C IS TO BE PERFORMED (MUST BE 1 OR GREATER), C XI = ARRAY OF DIMENSION NIP CONTAINING THE X C COORDINATES OF THE OUTPUT POINTS, C YI = ARRAY OF DIMENSION NIP CONTAINING THE Y C COORDINATES OF THE OUTPUT POINTS. C THE OUTPUT PARAMETER IS C ZI = ARRAY OF DIMENSION NIP WHERE INTERPOLATED Z C VALUES ARE TO BE STORED. C THE OTHER PARAMETERS ARE C IWK = INTEGER ARRAY OF DIMENSION C MAX0(31,27+NCP)*NDP+NIP C USED INTERNALLY AS A WORK AREA, C WK = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A C WORK AREA. C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND C YD ARRAYS MUST BE MADE WITH MD=1. THE CALL WITH MD=2 MUST BE C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS. THE CALL WITH C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP, C AND NIP VALUES AND WITH THE SAME CONTENTS OF THE XD, YD, XI, C AND YI ARRAYS. BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED. C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM- C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, C THEREFORE, SYSTEM DEPENDENT. C THIS SUBROUTINE CALLS THE IDCLDP, IDLCTN, IDPDRV, IDPTIP, AND C IDTANG SUBROUTINES. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),ZD(100),XI(1000),YI(1000), 1 ZI(1000),IWK(4100),WK(800) COMMON/IDLC/NIT COMMON/IDPI/ITPV DATA LUN/6/ C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES. C (FOR MD=1,2,3) 10 MD0=MD NCP0=NCP NDP0=NDP NIP0=NIP C ERROR CHECK. (FOR MD=1,2,3) 20 IF(MD0.LT.1.OR.MD0.GT.3) GO TO 90 IF(NCP0.LT.2.OR.NCP0.GE.NDP0) GO TO 90 IF(NDP0.LT.4) GO TO 90 IF(NIP0.LT.1) GO TO 90 IF(MD0.GE.2) GO TO 21 IWK(1)=NCP0 IWK(2)=NDP0 GO TO 22 21 NCPPV=IWK(1) NDPPV=IWK(2) IF(NCP0.NE.NCPPV) GO TO 90 IF(NDP0.NE.NDPPV) GO TO 90 22 IF(MD0.GE.3) GO TO 23 IWK(3)=NIP GO TO 30 23 NIPPV=IWK(3) IF(NIP0.NE.NIPPV) GO TO 90 C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY. (FOR MD=1,2,3) 30 JWIPT=16 JWIWL=6*NDP0+1 JWIWK=JWIWL JWIPL=24*NDP0+1 JWIWP=30*NDP0+1 JWIPC=27*NDP0+1 JWIT0=MAX0(31,27+NCP0)*NDP0 C TRIANGULATES THE X-Y PLANE. (FOR MD=1) 40 IF(MD0.GT.1) GO TO 50 CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL), 1 IWK(JWIWL),IWK(JWIWP),WK) IWK(5)=NT IWK(6)=NL IF(NT.EQ.0) RETURN C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1) 50 IF(MD0.GT.1) GO TO 60 CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC)) IF(IWK(JWIPC).EQ.0) RETURN C LOCATES ALL POINTS AT WHICH INTERPOLATION IS TO BE PERFORMED. C (FOR MD=1,2) 60 IF(MD0.EQ.3) GO TO 70 NIT=0 JWIT=JWIT0 DO 61 IIP=1,NIP0 JWIT=JWIT+1 CALL IDLCTN(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL), 1 XI(IIP),YI(IIP),IWK(JWIT),IWK(JWIWK),WK) 61 CONTINUE C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS. C (FOR MD=1,2,3) 70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK) C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3) 80 ITPV=0 JWIT=JWIT0 DO 81 IIP=1,NIP0 JWIT=JWIT+1 CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK, 1 IWK(JWIT),XI(IIP),YI(IIP),ZI(IIP)) 81 CONTINUE RETURN C ERROR EXIT 90 WRITE (LUN,2090) MD0,NCP0,NDP0,NIP0 RETURN C FORMAT STATEMENT FOR ERROR MESSAGE 2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S)./ 1 7H MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6, 2 10X,5HNIP =,I6/ 3 35H ERROR DETECTED IN ROUTINE IDBVIP/) END SUBROUTINE IDCLDP(NDP,XD,YD,NCP,IPC) ID002720 C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST C TO EACH OF THE DATA POINT. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y C COORDINATES OF THE DATA POINTS, C NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA C POINTS. C THE OUTPUT PARAMETER IS C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE C POINT NUMBERS OF NCP DATA POINTS CLOSEST TO C EACH OF THE NDP DATA POINTS ARE TO BE STORED. C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST C NOT EXCEED 25. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, C THEREFORE, SYSTEM DEPENDENT. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),IPC(400) DIMENSION DSQ0(25),IPC0(25) DATA NCPMX/25/, LUN/6/ C STATEMENT FUNCTION DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2 C PRELIMINARY PROCESSING 10 NDP0=NDP NCP0=NCP IF(NDP0.LT.2) GO TO 90 IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0) GO TO 90 C CALCULATION 20 DO 59 IP1=1,NDP0 C - SELECTS NCP POINTS. X1=XD(IP1) Y1=YD(IP1) J1=0 DSQMX=0.0 DO 22 IP2=1,NDP0 IF(IP2.EQ.IP1) GO TO 22 DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2)) J1=J1+1 DSQ0(J1)=DSQI IPC0(J1)=IP2 IF(DSQI.LE.DSQMX) GO TO 21 DSQMX=DSQI JMX=J1 21 IF(J1.GE.NCP0) GO TO 23 22 CONTINUE 23 IP2MN=IP2+1 IF(IP2MN.GT.NDP0) GO TO 30 DO 25 IP2=IP2MN,NDP0 IF(IP2.EQ.IP1) GO TO 25 DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2)) IF(DSQI.GE.DSQMX) GO TO 25 DSQ0(JMX)=DSQI IPC0(JMX)=IP2 DSQMX=0.0 DO 24 J1=1,NCP0 IF(DSQ0(J1).LE.DSQMX) GO TO 24 DSQMX=DSQ0(J1) JMX=J1 24 CONTINUE 25 CONTINUE C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR. 30 IP2=IPC0(1) DX12=XD(IP2)-X1 DY12=YD(IP2)-Y1 DO 31 J3=2,NCP0 IP3=IPC0(J3) DX13=XD(IP3)-X1 DY13=YD(IP3)-Y1 IF((DY13*DX12-DX13*DY12).NE.0.0) GO TO 50 31 CONTINUE C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT. 40 NCLPT=0 DO 43 IP3=1,NDP0 IF(IP3.EQ.IP1) GO TO 43 DO 41 J4=1,NCP0 IF(IP3.EQ.IPC0(J4)) GO TO 43 41 CONTINUE DX13=XD(IP3)-X1 DY13=YD(IP3)-Y1 IF((DY13*DX12-DX13*DY12).EQ.0.0) GO TO 43 DSQI=DSQF(X1,Y1,XD(IP3),YD(IP3)) IF(NCLPT.EQ.0) GO TO 42 IF(DSQI.GE.DSQMN) GO TO 43 42 NCLPT=1 DSQMN=DSQI IP3MN=IP3 43 CONTINUE IF(NCLPT.EQ.0) GO TO 91 DSQMX=DSQMN IPC0(JMX)=IP3MN C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY. 50 J1=(IP1-1)*NCP0 DO 51 J2=1,NCP0 J1=J1+1 IPC(J1)=IPC0(J2) 51 CONTINUE 59 CONTINUE RETURN C ERROR EXIT 90 WRITE (LUN,2090) GO TO 92 91 WRITE (LUN,2091) 92 WRITE (LUN,2092) NDP0,NCP0 IPC(1)=0 RETURN C FORMAT STATEMENTS FOR ERROR MESSAGES 2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S).) 2091 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS.) 2092 FORMAT(8H NDP =,I5,5X,5HNCP =,I5/ 1 35H ERROR DETECTED IN ROUTINE IDCLDP/) END SUBROUTINE IDGRID(XD, YD, NT, IPT, NL, IPL, NXI, NYI, XI, YI, IDG 10 * NGP, IGP) C THIS SUBROUTINE ORGANIZES GRID POINTS FOR SURFACE FITTING BY C SORTING THEM IN ASCENDING ORDER OF TRIANGLE NUMBERS AND OF THE C BORDER LINE SEGMENT NUMBER. C THE INPUT PARAMETERS ARE C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y C COORDINATES OF THE DATA POINTS, WHERE NDP IS THE C NUMBER OF THE DATA POINTS, C NT = NUMBER OF TRIANGLES, C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES, C NL = NUMBER OF BORDER LINE SEGMENTS, C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE C POINT NUMBERS OF THE END POINTS OF THE BORDER C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE C NUMBERS, C NXI = NUMBER OF GRID POINTS IN THE X COORDINATE, C NYI = NUMBER OF GRID POINTS IN THE Y COORDINATE, C XI,YI = ARRAYS OF DIMENSION NXI AND NYI CONTAINING C THE X AND Y COORDINATES OF THE GRID POINTS, C RESPECTIVELY. C THE OUTPUT PARAMETERS ARE C NGP = INTEGER ARRAY OF DIMENSION 2*(NT+2*NL) WHERE THE C NUMBER OF GRID POINTS THAT BELONG TO EACH OF THE C TRIANGLES OR OF THE BORDER LINE SEGMENTS ARE TO C BE STORED, C IGP = INTEGER ARRAY OF DIMENSION NXI*NYI WHERE THE C GRID POINT NUMBERS ARE TO BE STORED IN ASCENDING C ORDER OF THE TRIANGLE NUMBER AND THE BORDER LINE C SEGMENT NUMBER. C DECLARATION STATEMENTS DIMENSION XD(100), YD(100), IPT(585), IPL(300), XI(101), * YI(101), NGP(800), IGP(10201) C STATEMENT FUNCTIONS SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3) SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2) C PRELIMINARY PROCESSING NT0 = NT NL0 = NL NXI0 = NXI NYI0 = NYI NXINYI = NXI0*NYI0 XIMN = AMIN1(XI(1),XI(NXI0)) XIMX = AMAX1(XI(1),XI(NXI0)) YIMN = AMIN1(YI(1),YI(NYI0)) YIMX = AMAX1(YI(1),YI(NYI0)) C DETERMINES GRID POINTS INSIDE THE DATA AREA. JNGP0 = 0 JNGP1 = 2*(NT0+2*NL0) + 1 JIGP0 = 0 JIGP1 = NXINYI + 1 DO 160 IT0=1,NT0 NGP0 = 0 NGP1 = 0 IT0T3 = IT0*3 IP1 = IPT(IT0T3-2) IP2 = IPT(IT0T3-1) IP3 = IPT(IT0T3) X1 = XD(IP1) Y1 = YD(IP1) X2 = XD(IP2) Y2 = YD(IP2) X3 = XD(IP3) Y3 = YD(IP3) XMN = AMIN1(X1,X2,X3) XMX = AMAX1(X1,X2,X3) YMN = AMIN1(Y1,Y2,Y3) YMX = AMAX1(Y1,Y2,Y3) INSD = 0 DO 20 IXI=1,NXI0 IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 10 IF (INSD.EQ.0) GO TO 20 IXIMX = IXI - 1 GO TO 30 10 IF (INSD.EQ.1) GO TO 20 INSD = 1 IXIMN = IXI 20 CONTINUE IF (INSD.EQ.0) GO TO 150 IXIMX = NXI0 30 DO 140 IYI=1,NYI0 YII = YI(IYI) IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 140 DO 130 IXI=IXIMN,IXIMX XII = XI(IXI) L = 0 IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 130, 40, 50 40 L = 1 50 IF (SIDE(X2,Y2,X3,Y3,XII,YII)) 130, 60, 70 60 L = 1 70 IF (SIDE(X3,Y3,X1,Y1,XII,YII)) 130, 80, 90 80 L = 1 90 IZI = NXI0*(IYI-1) + IXI IF (L.EQ.1) GO TO 100 NGP0 = NGP0 + 1 JIGP0 = JIGP0 + 1 IGP(JIGP0) = IZI GO TO 130 100 IF (JIGP1.GT.NXINYI) GO TO 120 DO 110 JIGP1I=JIGP1,NXINYI IF (IZI.EQ.IGP(JIGP1I)) GO TO 130 110 CONTINUE 120 NGP1 = NGP1 + 1 JIGP1 = JIGP1 - 1 IGP(JIGP1) = IZI 130 CONTINUE 140 CONTINUE 150 JNGP0 = JNGP0 + 1 NGP(JNGP0) = NGP0 JNGP1 = JNGP1 - 1 NGP(JNGP1) = NGP1 160 CONTINUE C DETERMINES GRID POINTS OUTSIDE THE DATA AREA. C - IN SEMI-INFINITE RECTANGULAR AREA. DO 450 IL0=1,NL0 NGP0 = 0 NGP1 = 0 IL0T3 = IL0*3 IP1 = IPL(IL0T3-2) IP2 = IPL(IL0T3-1) X1 = XD(IP1) Y1 = YD(IP1) X2 = XD(IP2) Y2 = YD(IP2) XMN = XIMN XMX = XIMX YMN = YIMN YMX = YIMX IF (Y2.GE.Y1) XMN = AMIN1(X1,X2) IF (Y2.LE.Y1) XMX = AMAX1(X1,X2) IF (X2.LE.X1) YMN = AMIN1(Y1,Y2) IF (X2.GE.X1) YMX = AMAX1(Y1,Y2) INSD = 0 DO 180 IXI=1,NXI0 IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 170 IF (INSD.EQ.0) GO TO 180 IXIMX = IXI - 1 GO TO 190 170 IF (INSD.EQ.1) GO TO 180 INSD = 1 IXIMN = IXI 180 CONTINUE IF (INSD.EQ.0) GO TO 310 IXIMX = NXI0 190 DO 300 IYI=1,NYI0 YII = YI(IYI) IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 300 DO 290 IXI=IXIMN,IXIMX XII = XI(IXI) L = 0 IF (SIDE(X1,Y1,X2,Y2,XII,YII)) 210, 200, 290 200 L = 1 210 IF (SPDT(X2,Y2,X1,Y1,XII,YII)) 290, 220, 230 220 L = 1 230 IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 290, 240, 250 240 L = 1 250 IZI = NXI0*(IYI-1) + IXI IF (L.EQ.1) GO TO 260 NGP0 = NGP0 + 1 JIGP0 = JIGP0 + 1 IGP(JIGP0) = IZI GO TO 290 260 IF (JIGP1.GT.NXINYI) GO TO 280 DO 270 JIGP1I=JIGP1,NXINYI IF (IZI.EQ.IGP(JIGP1I)) GO TO 290 270 CONTINUE 280 NGP1 = NGP1 + 1 JIGP1 = JIGP1 - 1 IGP(JIGP1) = IZI 290 CONTINUE 300 CONTINUE 310 JNGP0 = JNGP0 + 1 NGP(JNGP0) = NGP0 JNGP1 = JNGP1 - 1 NGP(JNGP1) = NGP1 C - IN SEMI-INFINITE TRIANGULAR AREA. NGP0 = 0 NGP1 = 0 ILP1 = MOD(IL0,NL0) + 1 ILP1T3 = ILP1*3 IP3 = IPL(ILP1T3-1) X3 = XD(IP3) Y3 = YD(IP3) XMN = XIMN XMX = XIMX YMN = YIMN YMX = YIMX IF (Y3.GE.Y2 .AND. Y2.GE.Y1) XMN = X2 IF (Y3.LE.Y2 .AND. Y2.LE.Y1) XMX = X2 IF (X3.LE.X2 .AND. X2.LE.X1) YMN = Y2 IF (X3.GE.X2 .AND. X2.GE.X1) YMX = Y2 INSD = 0 DO 330 IXI=1,NXI0 IF (XI(IXI).GE.XMN .AND. XI(IXI).LE.XMX) GO TO 320 IF (INSD.EQ.0) GO TO 330 IXIMX = IXI - 1 GO TO 340 320 IF (INSD.EQ.1) GO TO 330 INSD = 1 IXIMN = IXI 330 CONTINUE IF (INSD.EQ.0) GO TO 440 IXIMX = NXI0 340 DO 430 IYI=1,NYI0 YII = YI(IYI) IF (YII.LT.YMN .OR. YII.GT.YMX) GO TO 430 DO 420 IXI=IXIMN,IXIMX XII = XI(IXI) L = 0 IF (SPDT(X1,Y1,X2,Y2,XII,YII)) 360, 350, 420 350 L = 1 360 IF (SPDT(X3,Y3,X2,Y2,XII,YII)) 380, 370, 420 370 L = 1 380 IZI = NXI0*(IYI-1) + IXI IF (L.EQ.1) GO TO 390 NGP0 = NGP0 + 1 JIGP0 = JIGP0 + 1 IGP(JIGP0) = IZI GO TO 420 390 IF (JIGP1.GT.NXINYI) GO TO 410 DO 400 JIGP1I=JIGP1,NXINYI IF (IZI.EQ.IGP(JIGP1I)) GO TO 420 400 CONTINUE 410 NGP1 = NGP1 + 1 JIGP1 = JIGP1 - 1 IGP(JIGP1) = IZI 420 CONTINUE 430 CONTINUE 440 JNGP0 = JNGP0 + 1 NGP(JNGP0) = NGP0 JNGP1 = JNGP1 - 1 NGP(JNGP1) = NGP1 450 CONTINUE RETURN END SUBROUTINE IDLCTN(NDP, XD, YD, NT, IPT, NL, IPL, XII, YII, ITI, IDL 10 * IWK, WK) C THIS SUBROUTINE LOCATES A POINT, I.E., DETERMINES TO WHAT TRI- C ANGLE A GIVEN POINT (XII,YII) BELONGS. WHEN THE GIVEN POINT C DOES NOT LIE INSIDE THE DATA AREA, THIS SUBROUTINE DETERMINES C THE BORDER LINE SEGMENT WHEN THE POINT LIES IN AN OUTSIDE C RECTANGULAR AREA, AND TWO BORDER LINE SEGMENTS WHEN THE POINT C LIES IN AN OUTSIDE TRIANGULAR AREA. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y C COORDINATES OF THE DATA POINTS, C NT = NUMBER OF TRIANGLES, C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES, C NL = NUMBER OF BORDER LINE SEGMENTS, C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE C POINT NUMBERS OF THE END POINTS OF THE BORDER C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE C NUMBERS, C XII,YII = X AND Y COORDINATES OF THE POINT TO BE C LOCATED. C THE OUTPUT PARAMETER IS C ITI = TRIANGLE NUMBER, WHEN THE POINT IS INSIDE THE C DATA AREA, OR C TWO BORDER LINE SEGMENT NUMBERS, IL1 AND IL2, C CODED TO IL1*(NT+NL)+IL2, WHEN THE POINT IS C OUTSIDE THE DATA AREA. C THE OTHER PARAMETERS ARE C IWK = INTEGER ARRAY OF DIMENSION 18*NDP USED INTER- C NALLY AS A WORK AREA, C WK = ARRAY OF DIMENSION 8*NDP USED INTERNALLY AS A C WORK AREA. C DECLARATION STATEMENTS DIMENSION XD(100), YD(100), IPT(585), IPL(300), IWK(1800), * WK(800) DIMENSION NTSC(9), IDSC(9) COMMON /IDLC/ NIT C STATEMENT FUNCTIONS SIDE(U1,V1,U2,V2,U3,V3) = (U1-U3)*(V2-V3) - (V1-V3)*(U2-U3) SPDT(U1,V1,U2,V2,U3,V3) = (U1-U2)*(U3-U2) + (V1-V2)*(V3-V2) C PRELIMINARY PROCESSING NDP0 = NDP NT0 = NT NL0 = NL NTL = NT0 + NL0 X0 = XII Y0 = YII C PROCESSING FOR A NEW SET OF DATA POINTS IF (NIT.NE.0) GO TO 80 NIT = 1 C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS. XMN = XD(1) XMX = XMN YMN = YD(1) YMX = YMN DO 10 IDP=2,NDP0 XI = XD(IDP) YI = YD(IDP) XMN = AMIN1(XI,XMN) XMX = AMAX1(XI,XMX) YMN = AMIN1(YI,YMN) YMX = AMAX1(YI,YMX) 10 CONTINUE XS1 = (XMN+XMN+XMX)/3.0 XS2 = (XMN+XMX+XMX)/3.0 YS1 = (YMN+YMN+YMX)/3.0 YS2 = (YMN+YMX+YMX)/3.0 C - DETERMINES AND STORES IN THE IWK ARRAY TRIANGLE NUMBERS OF C - THE TRIANGLES ASSOCIATED WITH EACH OF THE NINE SECTIONS. DO 20 ISC=1,9 NTSC(ISC) = 0 IDSC(ISC) = 0 20 CONTINUE IT0T3 = 0 JWK = 0 DO 70 IT0=1,NT0 IT0T3 = IT0T3 + 3 I1 = IPT(IT0T3-2) I2 = IPT(IT0T3-1) I3 = IPT(IT0T3) XMN = AMIN1(XD(I1),XD(I2),XD(I3)) XMX = AMAX1(XD(I1),XD(I2),XD(I3)) YMN = AMIN1(YD(I1),YD(I2),YD(I3)) YMX = AMAX1(YD(I1),YD(I2),YD(I3)) IF (YMN.GT.YS1) GO TO 30 IF (XMN.LE.XS1) IDSC(1) = 1 IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(2) = 1 IF (XMX.GE.XS2) IDSC(3) = 1 30 IF (YMX.LT.YS1 .OR. YMN.GT.YS2) GO TO 40 IF (XMN.LE.XS1) IDSC(4) = 1 IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(5) = 1 IF (XMX.GE.XS2) IDSC(6) = 1 40 IF (YMX.LT.YS2) GO TO 50 IF (XMN.LE.XS1) IDSC(7) = 1 IF (XMX.GE.XS1 .AND. XMN.LE.XS2) IDSC(8) = 1 IF (XMX.GE.XS2) IDSC(9) = 1 50 DO 60 ISC=1,9 IF (IDSC(ISC).EQ.0) GO TO 60 JIWK = 9*NTSC(ISC) + ISC IWK(JIWK) = IT0 NTSC(ISC) = NTSC(ISC) + 1 IDSC(ISC) = 0 60 CONTINUE C - STORES IN THE WK ARRAY THE MINIMUM AND MAXIMUM OF THE X AND C - Y COORDINATE VALUES FOR EACH OF THE TRIANGLE. JWK = JWK + 4 WK(JWK-3) = XMN WK(JWK-2) = XMX WK(JWK-1) = YMN WK(JWK) = YMX 70 CONTINUE GO TO 110 C CHECKS IF IN THE SAME TRIANGLE AS PREVIOUS. 80 IT0 = ITIPV IF (IT0.GT.NT0) GO TO 90 IT0T3 = IT0*3 IP1 = IPT(IT0T3-2) X1 = XD(IP1) Y1 = YD(IP1) IP2 = IPT(IT0T3-1) X2 = XD(IP2) Y2 = YD(IP2) IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110 IP3 = IPT(IT0T3) X3 = XD(IP3) Y3 = YD(IP3) IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 110 IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 110 GO TO 170 C CHECKS IF ON THE SAME BORDER LINE SEGMENT. 90 IL1 = IT0/NTL IL2 = IT0 - IL1*NTL IL1T3 = IL1*3 IP1 = IPL(IL1T3-2) X1 = XD(IP1) Y1 = YD(IP1) IP2 = IPL(IL1T3-1) X2 = XD(IP2) Y2 = YD(IP2) IF (IL2.NE.IL1) GO TO 100 IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 110 IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 110 IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110 GO TO 170 C CHECKS IF BETWEEN THE SAME TWO BORDER LINE SEGMENTS. 100 IF (SPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 110 IP3 = IPL(3*IL2-1) X3 = XD(IP3) Y3 = YD(IP3) IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 170 C LOCATES INSIDE THE DATA AREA. C - DETERMINES THE SECTION IN WHICH THE POINT IN QUESTION LIES. 110 ISC = 1 IF (X0.GE.XS1) ISC = ISC + 1 IF (X0.GE.XS2) ISC = ISC + 1 IF (Y0.GE.YS1) ISC = ISC + 3 IF (Y0.GE.YS2) ISC = ISC + 3 C - SEARCHES THROUGH THE TRIANGLES ASSOCIATED WITH THE SECTION. NTSCI = NTSC(ISC) IF (NTSCI.LE.0) GO TO 130 JIWK = -9 + ISC DO 120 ITSC=1,NTSCI JIWK = JIWK + 9 IT0 = IWK(JIWK) JWK = IT0*4 IF (X0.LT.WK(JWK-3)) GO TO 120 IF (X0.GT.WK(JWK-2)) GO TO 120 IF (Y0.LT.WK(JWK-1)) GO TO 120 IF (Y0.GT.WK(JWK)) GO TO 120 IT0T3 = IT0*3 IP1 = IPT(IT0T3-2) X1 = XD(IP1) Y1 = YD(IP1) IP2 = IPT(IT0T3-1) X2 = XD(IP2) Y2 = YD(IP2) IF (SIDE(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 120 IP3 = IPT(IT0T3) X3 = XD(IP3) Y3 = YD(IP3) IF (SIDE(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 120 IF (SIDE(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 120 GO TO 170 120 CONTINUE C LOCATES OUTSIDE THE DATA AREA. 130 DO 150 IL1=1,NL0 IL1T3 = IL1*3 IP1 = IPL(IL1T3-2) X1 = XD(IP1) Y1 = YD(IP1) IP2 = IPL(IL1T3-1) X2 = XD(IP2) Y2 = YD(IP2) IF (SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 150 IF (SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 140 IF (SIDE(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 150 IL2 = IL1 GO TO 160 140 IL2 = MOD(IL1,NL0) + 1 IP3 = IPL(3*IL2-1) X3 = XD(IP3) Y3 = YD(IP3) IF (SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 160 150 CONTINUE IT0 = 1 GO TO 170 160 IT0 = IL1*NTL + IL2 C NORMAL EXIT 170 ITI = IT0 ITIPV = IT0 RETURN END SUBROUTINE IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD) ID008940 C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND C SECOND ORDER AT THE DATA POINTS. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X, C Y, AND Z COORDINATES OF THE DATA POINTS, C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI- C MATING PARTIAL DERIVATIVES AT EACH DATA POINT, C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING C THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO C EACH OF THE NDP DATA POINTS. C THE OUTPUT PARAMETER IS C PD = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED C ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA C POINTS ARE TO BE STORED. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),ZD(100),IPC(400),PD(500) REAL NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY C PRELIMINARY PROCESSING 10 NDP0=NDP NCP0=NCP NCPM1=NCP0-1 C ESTIMATION OF ZX AND ZY 20 DO 24 IP0=1,NDP0 X0=XD(IP0) Y0=YD(IP0) Z0=ZD(IP0) NMX=0.0 NMY=0.0 NMZ=0.0 JIPC0=NCP0*(IP0-1) DO 23 IC1=1,NCPM1 JIPC=JIPC0+IC1 IPI=IPC(JIPC) DX1=XD(IPI)-X0 DY1=YD(IPI)-Y0 DZ1=ZD(IPI)-Z0 IC2MN=IC1+1 DO 22 IC2=IC2MN,NCP0 JIPC=JIPC0+IC2 IPI=IPC(JIPC) DX2=XD(IPI)-X0 DY2=YD(IPI)-Y0 DNMZ=DX1*DY2-DY1*DX2 IF(DNMZ.EQ.0.0) GO TO 22 DZ2=ZD(IPI)-Z0 DNMX=DY1*DZ2-DZ1*DY2 DNMY=DZ1*DX2-DX1*DZ2 IF(DNMZ.GE.0.0) GO TO 21 DNMX=-DNMX DNMY=-DNMY DNMZ=-DNMZ 21 NMX=NMX+DNMX NMY=NMY+DNMY NMZ=NMZ+DNMZ 22 CONTINUE 23 CONTINUE JPD0=5*IP0 PD(JPD0-4)=-NMX/NMZ PD(JPD0-3)=-NMY/NMZ 24 CONTINUE C ESTIMATION OF ZXX, ZXY, AND ZYY 30 DO 34 IP0=1,NDP0 JPD0=JPD0+5 X0=XD(IP0) JPD0=5*IP0 Y0=YD(IP0) ZX0=PD(JPD0-4) ZY0=PD(JPD0-3) NMXX=0.0 NMXY=0.0 NMYX=0.0 NMYY=0.0 NMZ =0.0 JIPC0=NCP0*(IP0-1) DO 33 IC1=1,NCPM1 JIPC=JIPC0+IC1 IPI=IPC(JIPC) DX1=XD(IPI)-X0 DY1=YD(IPI)-Y0 JPD=5*IPI DZX1=PD(JPD-4)-ZX0 DZY1=PD(JPD-3)-ZY0 IC2MN=IC1+1 DO 32 IC2=IC2MN,NCP0 JIPC=JIPC0+IC2 IPI=IPC(JIPC) DX2=XD(IPI)-X0 DY2=YD(IPI)-Y0 DNMZ =DX1*DY2 -DY1*DX2 IF(DNMZ.EQ.0.0) GO TO 32 JPD=5*IPI DZX2=PD(JPD-4)-ZX0 DZY2=PD(JPD-3)-ZY0 DNMXX=DY1*DZX2-DZX1*DY2 DNMXY=DZX1*DX2-DX1*DZX2 DNMYX=DY1*DZY2-DZY1*DY2 DNMYY=DZY1*DX2-DX1*DZY2 IF(DNMZ.GE.0.0) GO TO 31 DNMXX=-DNMXX DNMXY=-DNMXY DNMYX=-DNMYX DNMYY=-DNMYY DNMZ =-DNMZ 31 NMXX=NMXX+DNMXX NMXY=NMXY+DNMXY NMYX=NMYX+DNMYX NMYY=NMYY+DNMYY NMZ =NMZ +DNMZ 32 CONTINUE 33 CONTINUE PD(JPD0-2)=-NMXX/NMZ PD(JPD0-1)=-(NMXY+NMYX)/(2.0*NMZ) PD(JPD0) =-NMYY/NMZ 34 CONTINUE RETURN END SUBROUTINE IDPTIP(XD,YD,ZD,NT,IPT,NL,IPL,PDD,ITI,XII,YII, ID010190 1 ZII) C THIS SUBROUTINE PERFORMS PUNCTUAL INTERPOLATION OR EXTRAPOLA- C TION, I.E., DETERMINES THE Z VALUE AT A POINT. C THE INPUT PARAMETERS ARE C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X, C Y, AND Z COORDINATES OF THE DATA POINTS, WHERE C NDP IS THE NUMBER OF THE DATA POINTS, C NT = NUMBER OF TRIANGLES, C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES, C NL = NUMBER OF BORDER LINE SEGMENTS, C IPL = INTEGER ARRAY OF DIMENSION 3*NL CONTAINING THE C POINT NUMBERS OF THE END POINTS OF THE BORDER C LINE SEGMENTS AND THEIR RESPECTIVE TRIANGLE C NUMBERS, C PDD = ARRAY OF DIMENSION 5*NDP CONTAINING THE PARTIAL C DERIVATIVES AT THE DATA POINTS, C ITI = TRIANGLE NUMBER OF THE TRIANGLE IN WHICH LIES C THE POINT FOR WHICH INTERPOLATION IS TO BE C PERFORMED, C XII,YII = X AND Y COORDINATES OF THE POINT FOR WHICH C INTERPOLATION IS TO BE PERFORMED. C THE OUTPUT PARAMETER IS C ZII = INTERPOLATED Z VALUE. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),ZD(100),IPT(585),IPL(300), 1 PDD(500) COMMON/IDPI/ITPV DIMENSION X(3),Y(3),Z(3),PD(15), 1 ZU(3),ZV(3),ZUU(3),ZUV(3),ZVV(3) REAL LU,LV EQUIVALENCE (P5,P50) C PRELIMINARY PROCESSING 10 IT0=ITI NTL=NT+NL IF(IT0.LE.NTL) GO TO 20 IL1=IT0/NTL IL2=IT0-IL1*NTL IF(IL1.EQ.IL2) GO TO 40 GO TO 60 C CALCULATION OF ZII BY INTERPOLATION. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. 20 IF(IT0.EQ.ITPV) GO TO 30 C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE C VERTEXES. 21 JIPT=3*(IT0-1) JPD=0 DO 23 I=1,3 JIPT=JIPT+1 IDP=IPT(JIPT) X(I)=XD(IDP) Y(I)=YD(IDP) Z(I)=ZD(IDP) JPDD=5*(IDP-1) DO 22 KPD=1,5 JPD=JPD+1 JPDD=JPDD+1 PD(JPD)=PDD(JPDD) 22 CONTINUE 23 CONTINUE C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM C AND VICE VERSA. 24 X0=X(1) Y0=Y(1) A=X(2)-X0 B=X(3)-X0 C=Y(2)-Y0 D=Y(3)-Y0 AD=A*D BC=B*C DLT=AD-BC AP= D/DLT BP=-B/DLT CP=-C/DLT DP= A/DLT C CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE C TRIANGLE FOR THE U-V COORDINATE SYSTEM. 25 AA=A*A ACT2=2.0*A*C CC=C*C AB=A*B ADBC=AD+BC CD=C*D BB=B*B BDT2=2.0*B*D DD=D*D DO 26 I=1,3 JPD=5*I ZU(I)=A*PD(JPD-4)+C*PD(JPD-3) ZV(I)=B*PD(JPD-4)+D*PD(JPD-3) ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD) ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD) ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD) 26 CONTINUE C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. 27 P00=Z(1) P10=ZU(1) P01=ZV(1) P20=0.5*ZUU(1) P11=ZUV(1) P02=0.5*ZVV(1) H1=Z(2)-P00-P10-P20 H2=ZU(2)-P10-ZUU(1) H3=ZUU(2)-ZUU(1) P30= 10.0*H1-4.0*H2+0.5*H3 P40=-15.0*H1+7.0*H2 -H3 P50= 6.0*H1-3.0*H2+0.5*H3 H1=Z(3)-P00-P01-P02 H2=ZV(3)-P01-ZVV(1) H3=ZVV(3)-ZVV(1) P03= 10.0*H1-4.0*H2+0.5*H3 P04=-15.0*H1+7.0*H2 -H3 P05= 6.0*H1-3.0*H2+0.5*H3 LU=SQRT(AA+CC) LV=SQRT(BB+DD) THXU=ATAN2(C,A) THUV=ATAN2(D,B)-THXU CSUV=COS(THUV) P41=5.0*LV*CSUV/LU*P50 P14=5.0*LU*CSUV/LV*P05 H1=ZV(2)-P01-P11-P41 H2=ZUV(2)-P11-4.0*P41 P21= 3.0*H1-H2 P31=-2.0*H1+H2 H1=ZU(3)-P10-P11-P14 H2=ZUV(3)-P11-4.0*P14 P12= 3.0*H1-H2 P13=-2.0*H1+H2 THUS=ATAN2(D-C,B-A)-THXU THSV=THUV-THUS AA= SIN(THSV)/LU BB=-COS(THSV)/LU CC= SIN(THUS)/LV DD= COS(THUS)/LV AC=AA*CC AD=AA*DD BC=BB*CC G1=AA*AC*(3.0*BC+2.0*AD) G2=CC*AC*(3.0*AD+2.0*BC) H1=-AA*AA*AA*(5.0*AA*BB*P50+(4.0*BC+AD)*P41) 1 -CC*CC*CC*(5.0*CC*DD*P05+(4.0*AD+BC)*P14) H2=0.5*ZVV(2)-P02-P12 H3=0.5*ZUU(3)-P20-P21 P22=(G1*H2+G2*H3-H1)/(G1+G2) P32=H2-P22 P23=H3-P22 ITPV=IT0 C CONVERTS XII AND YII TO U-V SYSTEM. 30 DX=XII-X0 DY=YII-Y0 U=AP*DX+BP*DY V=CP*DX+DP*DY C EVALUATES THE POLYNOMIAL. 31 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05)))) P1=P10+V*(P11+V*(P12+V*(P13+V*P14))) P2=P20+V*(P21+V*(P22+V*P23)) P3=P30+V*(P31+V*P32) P4=P40+V*P41 ZII=P0+U*(P1+U*(P2+U*(P3+U*(P4+U*P5)))) RETURN C CALCULATION OF ZII BY EXTRAPOLATION IN THE RECTANGLE. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. 40 IF(IT0.EQ.ITPV) GO TO 50 C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE END C POINTS OF THE BORDER LINE SEGMENT. 41 JIPL=3*(IL1-1) JPD=0 DO 43 I=1,2 JIPL=JIPL+1 IDP=IPL(JIPL) X(I)=XD(IDP) Y(I)=YD(IDP) Z(I)=ZD(IDP) JPDD=5*(IDP-1) DO 42 KPD=1,5 JPD=JPD+1 JPDD=JPDD+1 PD(JPD)=PDD(JPDD) 42 CONTINUE 43 CONTINUE C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM C AND VICE VERSA. 44 X0=X(1) Y0=Y(1) A=Y(2)-Y(1) B=X(2)-X(1) C=-B D=A AD=A*D BC=B*C DLT=AD-BC AP= D/DLT BP=-B/DLT CP=-BP DP= AP C CONVERTS THE PARTIAL DERIVATIVES AT THE END POINTS OF THE C BORDER LINE SEGMENT FOR THE U-V COORDINATE SYSTEM. 45 AA=A*A ACT2=2.0*A*C CC=C*C AB=A*B ADBC=AD+BC CD=C*D BB=B*B BDT2=2.0*B*D DD=D*D DO 46 I=1,2 JPD=5*I ZU(I)=A*PD(JPD-4)+C*PD(JPD-3) ZV(I)=B*PD(JPD-4)+D*PD(JPD-3) ZUU(I)=AA*PD(JPD-2)+ACT2*PD(JPD-1)+CC*PD(JPD) ZUV(I)=AB*PD(JPD-2)+ADBC*PD(JPD-1)+CD*PD(JPD) ZVV(I)=BB*PD(JPD-2)+BDT2*PD(JPD-1)+DD*PD(JPD) 46 CONTINUE C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. 47 P00=Z(1) P10=ZU(1) P01=ZV(1) P20=0.5*ZUU(1) P11=ZUV(1) P02=0.5*ZVV(1) H1=Z(2)-P00-P01-P02 H2=ZV(2)-P01-ZVV(1) H3=ZVV(2)-ZVV(1) P03= 10.0*H1-4.0*H2+0.5*H3 P04=-15.0*H1+7.0*H2 -H3 P05= 6.0*H1-3.0*H2+0.5*H3 H1=ZU(2)-P10-P11 H2=ZUV(2)-P11 P12= 3.0*H1-H2 P13=-2.0*H1+H2 P21=0.0 P23=-ZUU(2)+ZUU(1) P22=-1.5*P23 ITPV=IT0 C CONVERTS XII AND YII TO U-V SYSTEM. 50 DX=XII-X0 DY=YII-Y0 U=AP*DX+BP*DY V=CP*DX+DP*DY C EVALUATES THE POLYNOMIAL. 51 P0=P00+V*(P01+V*(P02+V*(P03+V*(P04+V*P05)))) P1=P10+V*(P11+V*(P12+V*P13)) P2=P20+V*(P21+V*(P22+V*P23)) ZII=P0+U*(P1+U*P2) RETURN C CALCULATION OF ZII BY EXTRAPOLATION IN THE TRIANGLE. C CHECKS IF THE NECESSARY COEFFICIENTS HAVE BEEN CALCULATED. 60 IF(IT0.EQ.ITPV) GO TO 70 C LOADS COORDINATE AND PARTIAL DERIVATIVE VALUES AT THE VERTEX C OF THE TRIANGLE. 61 JIPL=3*IL2-2 IDP=IPL(JIPL) X(1)=XD(IDP) Y(1)=YD(IDP) Z(1)=ZD(IDP) JPDD=5*(IDP-1) DO 62 KPD=1,5 JPDD=JPDD+1 PD(KPD)=PDD(JPDD) 62 CONTINUE C CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL. 67 P00=Z(1) P10=PD(1) P01=PD(2) P20=0.5*PD(3) P11=PD(4) P02=0.5*PD(5) ITPV=IT0 C CONVERTS XII AND YII TO U-V SYSTEM. 70 U=XII-X(1) V=YII-Y(1) C EVALUATES THE POLYNOMIAL. 71 P0=P00+V*(P01+V*P02) P1=P10+V*P11 ZII=P0+U*(P1+U*P20) RETURN END SUBROUTINE IDSFFT(MD,NCP,NDP,XD,YD,ZD,NXI,NYI,XI,YI,ZI, ID013070 1 IWK,WK) C THIS SUBROUTINE PERFORMS SMOOTH SURFACE FITTING WHEN THE PRO- C JECTIONS OF THE DATA POINTS IN THE X-Y PLANE ARE IRREGULARLY C DISTRIBUTED IN THE PLANE. C THE INPUT PARAMETERS ARE C MD = MODE OF COMPUTATION (MUST BE 1, 2, OR 3), C = 1 FOR NEW NCP AND/OR NEW XD-YD, C = 2 FOR OLD NCP, OLD XD-YD, NEW XI-YI, C = 3 FOR OLD NCP, OLD XD-YD, OLD XI-YI, C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI- C MATING PARTIAL DERIVATIVES AT EACH DATA POINT C (MUST BE 2 OR GREATER, BUT SMALLER THAN NDP), C NDP = NUMBER OF DATA POINTS (MUST BE 4 OR GREATER), C XD = ARRAY OF DIMENSION NDP CONTAINING THE X C COORDINATES OF THE DATA POINTS, C YD = ARRAY OF DIMENSION NDP CONTAINING THE Y C COORDINATES OF THE DATA POINTS, C ZD = ARRAY OF DIMENSION NDP CONTAINING THE Z C COORDINATES OF THE DATA POINTS, C NXI = NUMBER OF OUTPUT GRID POINTS IN THE X COORDINATE C (MUST BE 1 OR GREATER), C NYI = NUMBER OF OUTPUT GRID POINTS IN THE Y COORDINATE C (MUST BE 1 OR GREATER), C XI = ARRAY OF DIMENSION NXI CONTAINING THE X C COORDINATES OF THE OUTPUT GRID POINTS, C YI = ARRAY OF DIMENSION NYI CONTAINING THE Y C COORDINATES OF THE OUTPUT GRID POINTS. C THE OUTPUT PARAMETER IS C ZI = DOUBLY-DIMENSIONED ARRAY OF DIMENSION (NXI,NYI), C WHERE THE INTERPOLATED Z VALUES AT THE OUTPUT C GRID POINTS ARE TO BE STORED. C THE OTHER PARAMETERS ARE C IWK = INTEGER ARRAY OF DIMENSION C MAX0(31,27+NCP)*NDP+NXI*NYI C USED INTERNALLY AS A WORK AREA, C WK = ARRAY OF DIMENSION 5*NDP USED INTERNALLY AS A C WORK AREA. C THE VERY FIRST CALL TO THIS SUBROUTINE AND THE CALL WITH A NEW C NCP VALUE, A NEW NDP VALUE, AND/OR NEW CONTENTS OF THE XD AND C YD ARRAYS MUST BE MADE WITH MD=1. THE CALL WITH MD=2 MUST BE C PRECEDED BY ANOTHER CALL WITH THE SAME NCP AND NDP VALUES AND C WITH THE SAME CONTENTS OF THE XD AND YD ARRAYS. THE CALL WITH C MD=3 MUST BE PRECEDED BY ANOTHER CALL WITH THE SAME NCP, NDP, C NXI, AND NYI VALUES AND WITH THE SAME CONTENTS OF THE XD, YD, C XI, AND YI ARRAYS. BETWEEN THE CALL WITH MD=2 OR MD=3 AND ITS C PRECEDING CALL, THE IWK AND WK ARRAYS MUST NOT BE DISTURBED. C USE OF A VALUE BETWEEN 3 AND 5 (INCLUSIVE) FOR NCP IS RECOM- C MENDED UNLESS THERE ARE EVIDENCES THAT DICTATE OTHERWISE. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, C THEREFORE, SYSTEM DEPENDENT. C THIS SUBROUTINE CALLS THE IDCLDP, IDGRID, IDPDRV, IDPTIP, AND C IDTANG SUBROUTINES. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),ZD(100),XI(101),YI(101), 1 ZI(10201),IWK(13301),WK(500) COMMON/IDPI/ITPV DATA LUN/6/ C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES. C (FOR MD=1,2,3) 10 MD0=MD NCP0=NCP NDP0=NDP NXI0=NXI NYI0=NYI C ERROR CHECK. (FOR MD=1,2,3) 20 IF(MD0.LT.1.OR.MD0.GT.3) GO TO 90 IF(NCP0.LT.2.OR.NCP0.GE.NDP0) GO TO 90 IF(NDP0.LT.4) GO TO 90 IF(NXI0.LT.1.OR.NYI0.LT.1) GO TO 90 IF(MD0.GE.2) GO TO 21 IWK(1)=NCP0 IWK(2)=NDP0 GO TO 22 21 NCPPV=IWK(1) NDPPV=IWK(2) IF(NCP0.NE.NCPPV) GO TO 90 IF(NDP0.NE.NDPPV) GO TO 90 22 IF(MD0.GE.3) GO TO 23 IWK(3)=NXI0 IWK(4)=NYI0 GO TO 30 23 NXIPV=IWK(3) NYIPV=IWK(4) IF(NXI0.NE.NXIPV) GO TO 90 IF(NYI0.NE.NYIPV) GO TO 90 C ALLOCATION OF STORAGE AREAS IN THE IWK ARRAY. (FOR MD=1,2,3) 30 JWIPT=16 JWIWL=6*NDP0+1 JWNGP0=JWIWL-1 JWIPL=24*NDP0+1 JWIWP=30*NDP0+1 JWIPC=27*NDP0+1 JWIGP0=MAX0(31,27+NCP0)*NDP0 C TRIANGULATES THE X-Y PLANE. (FOR MD=1) 40 IF(MD0.GT.1) GO TO 50 CALL IDTANG(NDP0,XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL), 1 IWK(JWIWL),IWK(JWIWP),WK) IWK(5)=NT IWK(6)=NL IF(NT.EQ.0) RETURN C DETERMINES NCP POINTS CLOSEST TO EACH DATA POINT. (FOR MD=1) 50 IF(MD0.GT.1) GO TO 60 CALL IDCLDP(NDP0,XD,YD,NCP0,IWK(JWIPC)) IF(IWK(JWIPC).EQ.0) RETURN C SORTS OUTPUT GRID POINTS IN ASCENDING ORDER OF THE TRIANGLE C NUMBER AND THE BORDER LINE SEGMENT NUMBER. (FOR MD=1,2) 60 IF(MD0.EQ.3) GO TO 70 CALL IDGRID(XD,YD,NT,IWK(JWIPT),NL,IWK(JWIPL),NXI0,NYI0, 1 XI,YI,IWK(JWNGP0+1),IWK(JWIGP0+1)) C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS. C (FOR MD=1,2,3) 70 CALL IDPDRV(NDP0,XD,YD,ZD,NCP0,IWK(JWIPC),WK) C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3) 80 ITPV=0 JIG0MX=0 JIG1MN=NXI0*NYI0+1 NNGP=NT+2*NL DO 89 JNGP=1,NNGP ITI=JNGP IF(JNGP.LE.NT) GO TO 81 IL1=(JNGP-NT+1)/2 IL2=(JNGP-NT+2)/2 IF(IL2.GT.NL) IL2=1 ITI=IL1*(NT+NL)+IL2 81 JWNGP=JWNGP0+JNGP NGP0=IWK(JWNGP) IF(NGP0.EQ.0) GO TO 86 JIG0MN=JIG0MX+1 JIG0MX=JIG0MX+NGP0 DO 82 JIGP=JIG0MN,JIG0MX JWIGP=JWIGP0+JIGP IZI=IWK(JWIGP) IYI=(IZI-1)/NXI0+1 IXI=IZI-NXI0*(IYI-1) CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK, 1 ITI,XI(IXI),YI(IYI),ZI(IZI)) 82 CONTINUE 86 JWNGP=JWNGP0+2*NNGP+1-JNGP NGP1=IWK(JWNGP) IF(NGP1.EQ.0) GO TO 89 JIG1MX=JIG1MN-1 JIG1MN=JIG1MN-NGP1 DO 87 JIGP=JIG1MN,JIG1MX JWIGP=JWIGP0+JIGP IZI=IWK(JWIGP) IYI=(IZI-1)/NXI0+1 IXI=IZI-NXI0*(IYI-1) CALL IDPTIP(XD,YD,ZD,NT,IWK(JWIPT),NL,IWK(JWIPL),WK, 1 ITI,XI(IXI),YI(IYI),ZI(IZI)) 87 CONTINUE 89 CONTINUE RETURN C ERROR EXIT 90 WRITE (LUN,2090) MD0,NCP0,NDP0,NXI0,NYI0 RETURN C FORMAT STATEMENT FOR ERROR MESSAGE 2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S)./ 1 7H MD =,I4,10X,5HNCP =,I6,10X,5HNDP =,I6, 2 10X,5HNXI =,I6,10X,5HNYI =,I6/ 3 35H ERROR DETECTED IN ROUTINE IDSFFT/) END SUBROUTINE IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK) ID014770 C THIS SUBROUTINE PERFORMS TRIANGULATION. IT DIVIDES THE X-Y C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS C CORRESPONDING TO THE BORDER LINE SEGMENTS. C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE C ARE LISTED COUNTER-CLOCKWISE. POINT NUMBERS OF THE END POINTS C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE, C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE. C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS, C THEREFORE, SYSTEM DEPENDENT. C THIS SUBROUTINE CALLS THE IDXCHG FUNCTION. C THE INPUT PARAMETERS ARE C NDP = NUMBER OF DATA POINTS, C XD = ARRAY OF DIMENSION NDP CONTAINING THE C X COORDINATES OF THE DATA POINTS, C YD = ARRAY OF DIMENSION NDP CONTAINING THE C Y COORDINATES OF THE DATA POINTS. C THE OUTPUT PARAMETERS ARE C NT = NUMBER OF TRIANGLES, C IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE C POINT NUMBERS OF THE VERTEXES OF THE (IT)TH C TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND, C (3*IT-1)ST, AND (3*IT)TH ELEMENTS, C IT=1,2,...,NT, C NL = NUMBER OF BORDER LINE SEGMENTS, C IPL = INTEGER ARRAY OF DIMENSION 6*NDP, WHERE THE C POINT NUMBERS OF THE END POINTS OF THE (IL)TH C BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE C NUMBER ARE TO BE STORED AS THE (3*IL-2)ND, C (3*IL-1)ST, AND (3*IL)TH ELEMENTS, C IL=1,2,..., NL. C THE OTHER PARAMETERS ARE C IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED C INTERNALLY AS A WORK AREA, C IWP = INTEGER ARRAY OF DIMENSION NDP USED C INTERNALLY AS A WORK AREA, C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A C WORK AREA. C DECLARATION STATEMENTS DIMENSION XD(100),YD(100),IPT(585),IPL(600), 1 IWL(1800),IWP(100),WK(100) DIMENSION ITF(2) DATA RATIO/1.0E-6/, NREP/100/, LUN/6/ C STATEMENT FUNCTIONS DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2 SIDE(U1,V1,U2,V2,U3,V3)=(V3-V1)*(U2-U1)-(U3-U1)*(V2-V1) C PRELIMINARY PROCESSING 10 NDP0=NDP NDPM1=NDP0-1 IF(NDP0.LT.4) GO TO 90 C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT. 20 DSQMN=DSQF(XD(1),YD(1),XD(2),YD(2)) IPMN1=1 IPMN2=2 DO 22 IP1=1,NDPM1 X1=XD(IP1) Y1=YD(IP1) IP1P1=IP1+1 DO 21 IP2=IP1P1,NDP0 DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2)) IF(DSQI.EQ.0.0) GO TO 91 IF(DSQI.GE.DSQMN) GO TO 21 DSQMN=DSQI IPMN1=IP1 IPMN2=IP2 21 CONTINUE 22 CONTINUE DSQ12=DSQMN XDMP=(XD(IPMN1)+XD(IPMN2))/2.0 YDMP=(YD(IPMN1)+YD(IPMN2))/2.0 C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT C NUMBERS IN THE IWP ARRAY. 30 JP1=2 DO 31 IP1=1,NDP0 IF(IP1.EQ.IPMN1.OR.IP1.EQ.IPMN2) GO TO 31 JP1=JP1+1 IWP(JP1)=IP1 WK(JP1)=DSQF(XDMP,YDMP,XD(IP1),YD(IP1)) 31 CONTINUE DO 33 JP1=3,NDPM1 DSQMN=WK(JP1) JPMN=JP1 DO 32 JP2=JP1,NDP0 IF(WK(JP2).GE.DSQMN) GO TO 32 DSQMN=WK(JP2) JPMN=JP2 32 CONTINUE ITS=IWP(JP1) IWP(JP1)=IWP(JPMN) IWP(JPMN)=ITS WK(JPMN)=WK(JP1) 33 CONTINUE C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE C FIRST THREE DATA POINTS ARE NOT COLLINEAR. 35 AR=DSQ12*RATIO X1=XD(IPMN1) Y1=YD(IPMN1) DX21=XD(IPMN2)-X1 DY21=YD(IPMN2)-Y1 DO 36 JP=3,NDP0 IP=IWP(JP) IF(ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21).GT.AR) 1 GO TO 37 36 CONTINUE GO TO 92 37 IF(JP.EQ.3) GO TO 40 JPMX=JP JP=JPMX+1 DO 38 JPC=4,JPMX JP=JP-1 IWP(JP)=IWP(JP-1) 38 CONTINUE IWP(3)=IP C FORMS THE FIRST TRIANGLE. STORES POINT NUMBERS OF THE VER- C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM- C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN C THE IPL ARRAY. 40 IP1=IPMN1 IP2=IPMN2 IP3=IWP(3) IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3)) 1 .GE.0.0) GO TO 41 IP1=IPMN2 IP2=IPMN1 41 NT0=1 NTT3=3 IPT(1)=IP1 IPT(2)=IP2 IPT(3)=IP3 NL0=3 NLT3=9 IPL(1)=IP1 IPL(2)=IP2 IPL(3)=1 IPL(4)=IP2 IPL(5)=IP3 IPL(6)=1 IPL(7)=IP3 IPL(8)=IP1 IPL(9)=1 C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE. 50 DO 79 JP1=4,NDP0 IP1=IWP(JP1) X1=XD(IP1) Y1=YD(IP1) C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS. IP2=IPL(1) JPMN=1 DXMN=XD(IP2)-X1 DYMN=YD(IP2)-Y1 DSQMN=DXMN**2+DYMN**2 ARMN=DSQMN*RATIO JPMX=1 DXMX=DXMN DYMX=DYMN DSQMX=DSQMN ARMX=ARMN DO 52 JP2=2,NL0 IP2=IPL(3*JP2-2) DX=XD(IP2)-X1 DY=YD(IP2)-Y1 AR=DY*DXMN-DX*DYMN IF(AR.GT.ARMN) GO TO 51 DSQI=DX**2+DY**2 IF(AR.GE.(-ARMN).AND.DSQI.GE.DSQMN) GO TO 51 JPMN=JP2 DXMN=DX DYMN=DY DSQMN=DSQI ARMN=DSQMN*RATIO 51 AR=DY*DXMX-DX*DYMX IF(AR.LT.(-ARMX)) GO TO 52 DSQI=DX**2+DY**2 IF(AR.LE.ARMX.AND.DSQI.GE.DSQMX) GO TO 52 JPMX=JP2 DXMX=DX DYMX=DY DSQMX=DSQI ARMX=DSQMX*RATIO 52 CONTINUE IF(JPMX.LT.JPMN) JPMX=JPMX+NL0 NSH=JPMN-1 IF(NSH.LE.0) GO TO 60 C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY. NSHT3=NSH*3 DO 53 JP2T3=3,NSHT3,3 JP3T3=JP2T3+NLT3 IPL(JP3T3-2)=IPL(JP2T3-2) IPL(JP3T3-1)=IPL(JP2T3-1) IPL(JP3T3) =IPL(JP2T3) 53 CONTINUE DO 54 JP2T3=3,NLT3,3 JP3T3=JP2T3+NSHT3 IPL(JP2T3-2)=IPL(JP3T3-2) IPL(JP2T3-1)=IPL(JP3T3-1) IPL(JP2T3) =IPL(JP3T3) 54 CONTINUE JPMX=JPMX-NSH C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY. 60 JWL=0 DO 64 JP2=JPMX,NL0 JP2T3=JP2*3 IPL1=IPL(JP2T3-2) IPL2=IPL(JP2T3-1) IT =IPL(JP2T3) C - - ADDS A TRIANGLE TO THE IPT ARRAY. NT0=NT0+1 NTT3=NTT3+3 IPT(NTT3-2)=IPL2 IPT(NTT3-1)=IPL1 IPT(NTT3) =IP1 C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY. IF(JP2.NE.JPMX) GO TO 61 IPL(JP2T3-1)=IP1 IPL(JP2T3) =NT0 61 IF(JP2.NE.NL0) GO TO 62 NLN=JPMX+1 NLNT3=NLN*3 IPL(NLNT3-2)=IP1 IPL(NLNT3-1)=IPL(1) IPL(NLNT3) =NT0 C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER C - - LINE SEGMENTS. 62 ITT3=IT*3 IPTI=IPT(ITT3-2) IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63 IPTI=IPT(ITT3-1) IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63 IPTI=IPT(ITT3) C - - CHECKS IF THE EXCHANGE IS NECESSARY. 63 IF(IDXCHG(XD,YD,IP1,IPTI,IPL1,IPL2).EQ.0) GO TO 64 C - - MODIFIES THE IPT ARRAY WHEN NECESSARY. IPT(ITT3-2)=IPTI IPT(ITT3-1)=IPL1 IPT(ITT3) =IP1 IPT(NTT3-1)=IPTI IF(JP2.EQ.JPMX) IPL(JP2T3)=IT IF(JP2.EQ.NL0.AND.IPL(3).EQ.IT) IPL(3)=NT0 C - - SETS FLAGS IN THE IWL ARRAY. JWL=JWL+4 IWL(JWL-3)=IPL1 IWL(JWL-2)=IPTI IWL(JWL-1)=IPTI IWL(JWL) =IPL2 64 CONTINUE NL0=NLN NLT3=NLNT3 NLF=JWL/2 IF(NLF.EQ.0) GO TO 79 C - IMPROVES TRIANGULATION. 70 NTT3P3=NTT3+3 DO 78 IREP=1,NREP DO 76 ILF=1,NLF ILFT2=ILF*2 IPL1=IWL(ILFT2-1) IPL2=IWL(ILFT2) C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF C - - THE FLAGGED LINE SEGMENT. NTF=0 DO 71 ITT3R=3,NTT3,3 ITT3=NTT3P3-ITT3R IPT1=IPT(ITT3-2) IPT2=IPT(ITT3-1) IPT3=IPT(ITT3) IF(IPL1.NE.IPT1.AND.IPL1.NE.IPT2.AND. 1 IPL1.NE.IPT3) GO TO 71 IF(IPL2.NE.IPT1.AND.IPL2.NE.IPT2.AND. 1 IPL2.NE.IPT3) GO TO 71 NTF=NTF+1 ITF(NTF)=ITT3/3 IF(NTF.EQ.2) GO TO 72 71 CONTINUE IF(NTF.LT.2) GO TO 76 C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE C - - ON THE LINE SEGMENT. 72 IT1T3=ITF(1)*3 IPTI1=IPT(IT1T3-2) IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73 IPTI1=IPT(IT1T3-1) IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73 IPTI1=IPT(IT1T3) 73 IT2T3=ITF(2)*3 IPTI2=IPT(IT2T3-2) IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74 IPTI2=IPT(IT2T3-1) IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74 IPTI2=IPT(IT2T3) C - - CHECKS IF THE EXCHANGE IS NECESSARY. 74 IF(IDXCHG(XD,YD,IPTI1,IPTI2,IPL1,IPL2).EQ.0) 1 GO TO 76 C - - MODIFIES THE IPT ARRAY WHEN NECESSARY. IPT(IT1T3-2)=IPTI1 IPT(IT1T3-1)=IPTI2 IPT(IT1T3) =IPL1 IPT(IT2T3-2)=IPTI2 IPT(IT2T3-1)=IPTI1 IPT(IT2T3) =IPL2 C - - SETS NEW FLAGS. JWL=JWL+8 IWL(JWL-7)=IPL1 IWL(JWL-6)=IPTI1 IWL(JWL-5)=IPTI1 IWL(JWL-4)=IPL2 IWL(JWL-3)=IPL2 IWL(JWL-2)=IPTI2 IWL(JWL-1)=IPTI2 IWL(JWL) =IPL1 DO 75 JLT3=3,NLT3,3 IPLJ1=IPL(JLT3-2) IPLJ2=IPL(JLT3-1) IF((IPLJ1.EQ.IPL1.AND.IPLJ2.EQ.IPTI2).OR. 1 (IPLJ2.EQ.IPL1.AND.IPLJ1.EQ.IPTI2)) 2 IPL(JLT3)=ITF(1) IF((IPLJ1.EQ.IPL2.AND.IPLJ2.EQ.IPTI1).OR. 1 (IPLJ2.EQ.IPL2.AND.IPLJ1.EQ.IPTI1)) 2 IPL(JLT3)=ITF(2) 75 CONTINUE 76 CONTINUE NLFC=NLF NLF=JWL/2 IF(NLF.EQ.NLFC) GO TO 79 C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND. JWL=0 JWL1MN=(NLFC+1)*2 NLFT2=NLF*2 DO 77 JWL1=JWL1MN,NLFT2,2 JWL=JWL+2 IWL(JWL-1)=IWL(JWL1-1) IWL(JWL) =IWL(JWL1) 77 CONTINUE NLF=JWL/2 78 CONTINUE 79 CONTINUE C REARRANGES THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE C ARE LISTED COUNTER-CLOCKWISE. 80 DO 81 ITT3=3,NTT3,3 IP1=IPT(ITT3-2) IP2=IPT(ITT3-1) IP3=IPT(ITT3) IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3)) 1 .GE.0.0) GO TO 81 IPT(ITT3-2)=IP2 IPT(ITT3-1)=IP1 81 CONTINUE NT=NT0 NL=NL0 RETURN C ERROR EXIT 90 WRITE (LUN,2090) NDP0 GO TO 93 91 WRITE (LUN,2091) NDP0,IP1,IP2,X1,Y1 GO TO 93 92 WRITE (LUN,2092) NDP0 93 WRITE (LUN,2093) NT=0 RETURN C FORMAT STATEMENTS 2090 FORMAT(1X/23H *** NDP LESS THAN 4./8H NDP =,I5) 2091 FORMAT(1X/29H *** IDENTICAL DATA POINTS./ 1 8H NDP =,I5,5X,5HIP1 =,I5,5X,5HIP2 =,I5, 2 5X,4HXD =,E12.4,5X,4HYD =,E12.4) 2092 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS./ 1 8H NDP =,I5) 2093 FORMAT(35H ERROR DETECTED IN ROUTINE IDTANG/) END FUNCTION IDXCHG(X,Y,I1,I2,I3,I4) ID018560 C THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO C TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION C BY C. L. LAWSON. C THE INPUT PARAMETERS ARE C X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA C POINTS, C I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1, P2, C P3, AND P4 THAT FORM A QUADRILATERAL WITH P3 C AND P4 CONNECTED DIAGONALLY. C THIS FUNCTION RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EX- C CHANGE IS NECESSARY, AND 0 (ZERO) OTHERWISE. C DECLARATION STATEMENTS DIMENSION X(100),Y(100) EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ), 1 (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ) C PRELIMINARY PROCESSING 10 X1=X(I1) Y1=Y(I1) X2=X(I2) Y2=Y(I2) X3=X(I3) Y3=Y(I3) X4=X(I4) Y4=Y(I4) C CALCULATION 20 IDX=0 U3=(Y2-Y3)*(X1-X3)-(X2-X3)*(Y1-Y3) U4=(Y1-Y4)*(X2-X4)-(X1-X4)*(Y2-Y4) IF(U3*U4.LE.0.0) GO TO 30 U1=(Y3-Y1)*(X4-X1)-(X3-X1)*(Y4-Y1) U2=(Y4-Y2)*(X3-X2)-(X4-X2)*(Y3-Y2) A1SQ=(X1-X3)**2+(Y1-Y3)**2 B1SQ=(X4-X1)**2+(Y4-Y1)**2 C1SQ=(X3-X4)**2+(Y3-Y4)**2 A2SQ=(X2-X4)**2+(Y2-Y4)**2 B2SQ=(X3-X2)**2+(Y3-Y2)**2 C3SQ=(X2-X1)**2+(Y2-Y1)**2 S1SQ=U1*U1/(C1SQ*AMAX1(A1SQ,B1SQ)) S2SQ=U2*U2/(C2SQ*AMAX1(A2SQ,B2SQ)) S3SQ=U3*U3/(C3SQ*AMAX1(A3SQ,B3SQ)) S4SQ=U4*U4/(C4SQ*AMAX1(A4SQ,B4SQ)) IF(AMIN1(S1SQ,S2SQ).LT.AMIN1(S3SQ,S4SQ)) IDX=1 30 IDXCHG=IDX RETURN END