c ------------------------------------------------------------------ c c driver for interpolation program c c anselmo@thermsa.eng.sunysb.edu c c modified to produce plotmtv-compatible data c program test parameter(nmax=3000) parameter(nmax_out=200) dimension xd(nmax),yd(nmax),zd(nmax) dimension xi(nmax_out),yi(nmax_out) dimension zi(nmax_out,nmax_out) dimension iwk(31*nmax+nmax*nmax) dimension wk(6*nmax) character *80 infile character *80 outfile write(*,*)'conversion program - random to uniform grid for CONTOUR' write(*,*)'input file:' read(*,999) infile write(*,*)'output file:' read(*,999) outfile 999 format(a80) c open files open(unit=8,file=infile) open(unit=9,file=outfile) c read data, calculate maximum and minimum i=1 xmax=-1e10 xmin=1e10 ymax=xmax ymin=xmin 100 continue read(8,*,end=200) xd(i),yd(i),zd(i) if(xd(i).gt.xmax) xmax=xd(i) if(xd(i).lt.xmin) xmin=xd(i) if(yd(i).gt.ymax) ymax=yd(i) if(yd(i).lt.ymin) ymin=yd(i) i=i+1 goto 100 200 continue close(unit=8) npts=i-1 c spit out data c do 300 i=1,npts c write(*,*) i,xd(i),yd(i),zd(i) c300 continue write(*,*) ' ' write(*,*) 'npts = ',npts write(*,*) 'xmax = ',xmax write(*,*) 'xmin = ',xmin write(*,*) 'ymax = ',ymax write(*,*) 'ymin = ',ymin dx=xmax-xmin dy=ymax-ymin write(*,*)'number of uniform grids in the x direction:' read(*,*) nxi write(*,*)'number of uniform grids in the y direction:' read(*,*) nyi c number of uniform grids in the x and y c nxi=10 c nyi=15 c create the grid values do 400 i=1,nxi frac=float(i-1)/float(nxi-1) xi(i)=xmin+frac*dx 400 continue do 500 i=1,nyi frac=float(i-1)/float(nyi-1) yi(i)=ymin+frac*dy 500 continue call IDSFFT (1,NPTS,XD,YD,ZD,NXI,NYI,NMAX_OUT,XI,YI,ZI,IWK,WK) c for simple x-y-z dump c do 600 i=1,nxi c do 600 j=1,nyi c write(9,*) xi(i),yi(j),zi(i,j) c600 continue c for CONTOUR-MTV format write(9,*) '$ DATA=CONTOUR' write(9,*) '% xmin=',xmin,' xmax=',xmax write(9,*) '% ymin=',ymin,' ymax=',ymax write(9,*) '% nx=',nxi,' ny=',nyi do 600 j=1,nyi write(9,*) (zi(i,j),i=1,nxi) 600 continue write(9,*) '$END' close(unit=9) stop end c ------------------------------------------------------------------ C C PACKAGE BIVAR NOTE: DOCUMENTATION FOR INDIVIDUAL ROUTINES C FOLLOWS THE GENERAL PACKAGE INFORMATION. C C LATEST REVISION JANUARY, 1985 C C PURPOSE TO PROVIDE BIVARIATE INTERPOLATION AND SMOOTH C SURFACE FITTING FOR VALUES GIVEN AT IRREGULARLY C DISTRIBUTED POINTS. C C THE RESULTING INTERPOLATING FUNCTION AND C ITS FIRST-ORDER PARTIAL DERIVATIVES ARE C CONTINUOUS. C C THE METHOD EMPLOYED IS LOCAL, I.E. A CHANGE C IN THE DATA IN ONE AREA OF THE PLANE DOES NOT C AFFECT THE INTERPOLATING FUNCTION EXCEPT IN C THAT LOCAL AREA. THIS IS ADVANTAGEOUS OVER C GLOBAL INTERPOLATION METHODS. C C ALSO, THE METHOD GIVES EXACT RESULTS WHEN ALL C POINTS LIE IN A PLANE. THIS IS ADVANTAGEOUS C OVER OTHER METHODS SUCH AS TWO-DIMENSIONAL C FOURIER SERIES INTERPOLATION. C C USAGE THIS PACKAGE CONTAINS TWO USER ENTRIES, C IDBVIP AND IDSFFT, BOTH REQUIRING INPUT C DATA TO BE GIVEN AT POINTS C ( X(I) , Y(I) ), I=1,...,N. C C IF THE USER DESIRES THE INTERPOLATED DATA C TO BE OUTPUT AT GRID POINTS, I.E. AT POINTS C ( XI(I) , YI(J) ), I=1,...,NX, J=1,...,NYI, C ROUTINE IDSFFT SHOULD BE USED. THIS IS C USEFUL FOR GENERATING AN INTERPOLATING C SURFACE. C C THE OTHER USER ENTRY POINT, IDBVIP, WILL C PRODUCE INTERPOLATED VALUES AT POINTS C ( XI(I) , YI(I) ), I=1,...,NIP. THIS IS C USEFUL FOR FILLING IN MISSING DATA POINTS C ON A GRID. C C I/O NONE, EXCEPT ERROR MESSAGES PRINTED VIA C ROUTINE ULIBER. C C PRECISION SINGLE C C REQUIRED LIBRARY ULIBER ON ULIB, WHICH IS AUTOMATICALLY C FILES LOADED ON NCAR'S CRAY MACHINES. C C LANGUAGE FORTRAN C C HISTORY THE ORIGINAL VERSION OF BIVAR WAS WRITTEN BY C HIROSHI AKIMA IN AUGUST 1975 AND REWRITTEN BY C HIM IN LATE 1976. IT WAS INCORPORATED INTO C NCAR'S PUBLIC SOFTWARE LIBRARIES IN JANUARY C 1977. IN AUGUST 1984 A NEW VERSION OF BIVAR, C INCORPORATING CHANGES DESCRIBED IN THE ROCKY C MOUNTAIN JOURNAL OF MATHEMATICS ARTICLE CITED C BELOW, WAS OBTAINED FROM DR. AKIMA BY MICHAEL C PERNICE OF NCAR'S SCIENTIFIC COMPUTING C DIVISION, WHO EVALUATED IT AND MADE IT C AVAILABLE IN FEBRUARY, 1985. C C PORTABILITY FORTRAN 66 C C ACCURACY ACCURATE TO MACHINE PRECISION ON THE INPUT C DATA POINTS. ACCURACY AT OTHER POINTS C GREATLY DEPENDS ON THE INPUT DATA. C C REFERENCES THE ORIGINAL METHOD IS DESCRIBED IN C C AKIMA, HIROSHI, 1978: A METHOD OF BIVARIATE C INTERPOLATION AND SMOOTH SURFACE FITTING C FOR VALUES GIVEN AT IRREGULARLY C DISTRIBUTED POINTS, C ACM-TOMS, VOL. 4, NO. 2, JUNE 1978. C C AN IMPROVEMENT TO THE ORIGINAL METHOD C WAS PRESENTED IN C C AKIMA, HIROSHI, 1984: ON ESTIMATING PARTIAL C DERIVATIVES FOR BIVARIATE INTERPOLATION C OF SCATTERED DATA, C ROCKY MOUNTAIN JOURNAL OF MATHEMATICS, C VOL. 14, NO. 1, WINTER 1984. C C METHOD THE X-Y PLANE IS DIVIDED INTO TRIANGULAR C CELLS, EACH CELL HAVING PROJECTIONS OF THREE C DATA POINTS IN THE PLANE AS ITS VERTICES, AND C A BIVARIATE QUINTIC POLYNOMIAL IN X AND Y IS C FITTED TO EACH TRIANGULAR CELL. C C THE COEFFICIENTS IN THE FITTED QUINTIC C POLYNOMIALS ARE DETERMINED BY CONTINUITY C REQUIREMENTS AND BY ESTIMATES OF PARTIAL C DERIVATIVES AT THE VERTICES AND ALONG THE C EDGES OF THE TRIANGLES. THE METHOD DESCRIBED C IN THE ROCKY MOUNTAIN JOURNAL REFERENCE C GUARANTEES THAT THE GENERATED SURFACE C DEPENDS CONTINUOUSLY ON THE TRIANGULATION. C C THE RESULTING INTERPOLATING FUNCTION IS C INVARIANT UNDER THE FOLLOWING TYPES OF LINEAR C COORDINATE TRANSFORMATIONS: C 1) A ROTATION OF X-Y COORDINATE SYSTEM C 2) LINEAR SCALE TRANSFORMATION OF Z-AXIS C 3) TILTING OF THE X-Y PLANE, I.E. NEW C COORDINATES (U,V,W) GIVEN BY C U = X C V = Y C W = Z + A*X + B*Y C WHERE A, B ARE ARBITRARY CONSTANTS. C C COMPLETE DETAILS OF THE METHOD ARE GIVEN C IN THE REFERENCE PUBLICATIONS. C C ******************************************************************** C C INDIVIDUAL USER ENTRY POINT DOCUMENTATION FOLLOWS. C C ******************************************************************** C C SUBROUTINE IDBVIP (MD,NDP,XD,YD,ZD,NIP,XI,YI,ZI,IWK,WK) C C DIMENSION OF XD(NDP), YD(NDP), ZD(NDP), XI(NIP), YI(NIP) C ARGUMENTS ZI(NIP), IWK(31*NDP+NIP), WK(8*NDP) C C PURPOSE TO PERFORM BIVARIATE INTERPOLATION WHEN THE C PROJECTIONS OF THE DATA POINTS IN THE X-Y C PLANE ARE IRREGULARLY DISTRIBUTED. C C USAGE CALL IDBVIP (MD,NDP,XD,YD,ZD,NIP,XI,YI,ZI, C IWK,WK) C C ARGUMENTS C C ON INPUT MD C MODE OF COMPUTATION (MUST BE 1, 2, OR 3, C ELSE AN ERROR RETURN OCCURS.) C = 1 IF THIS IS THE FIRST CALL TO THIS C SUBROUTINE, OR IF THE VALUE OF NDP C HAS BEEN CHANGED FROM THE PREVIOUS C CALL, OR IF THE CONTENTS OF THE XD C OR YD ARRAYS HAVE BEEN CHANGED FROM C THE PREVIOUS CALL. C = 2 IF THE VALUES OF NDP AND THE XD AND C YD ARRAYS ARE UNCHANGED FROM THE C PREVIOUS CALL, BUT NEW VALUES FOR C XI, YI ARE BEING USED. IF MD = 2 C AND NDP HAS BEEN CHANGED SINCE THE C PREVIOUS CALL TO IDBVIP, AN ERROR C RETURN OCCURS. C = 3 IF THE VALUES OF NDP, NIP, XD, C YD, XI, YI ARE UNCHANGED FROM THE C PREVIOUS CALL, I.E. IF THE ONLY C CHANGE ON INPUT TO IDBVIP IS IN THE C ZD ARRAY. IF MD=3 AND NDP OR NIP HAS C BEEN CHANGED SINCE THE PREVIOUS CALL C TO IDBVIP, AN ERROR RETURN OCCURS. C C BETWEEN THE CALL WITH MD=2 OR MD=3 AND C THE PRECEDING CALL, THE IWK AND WK WORK C ARRAYS SHOULD NOT BE DISTURBED. C C NDP C NUMBER OF DATA POINTS (MUST BE 4 OR C GREATER, ELSE AN ERROR RETURN OCCURS). C C XD C ARRAY OF DIMENSION NDP CONTAINING THE C X COORDINATES OF THE DATA POINTS. C C YD C ARRAY OF DIMENSION NDP CONTAINING THE C Y COORDINATES OF THE DATA POINTS. C C ZD C ARRAY OF DIMENSION NDP CONTAINING THE C Z COORDINATES OF THE DATA POINTS. C C NIP C THE NUMBER OF OUTPUT POINTS AT WHICH C INTERPOLATION IS TO BE PERFORMED (MUST BE C 1 OR GREATER, ELSE AN ERROR RETURN OCCURS). C C XI C ARRAY OF DIMENSION NIP CONTAINING THE X C COORDINATES OF THE OUTPUT POINTS. C C YI C ARRAY OF DIMENSION NIP CONTAINING THE Y C COORDINATES OF THE OUTPUT POINTS. C C IWK C INTEGER WORK ARRAY OF DIMENSION AT LEAST C 31*NDP + NIP C C WK C REAL WORK ARRAY OF DIMENSION AT LEAST 8*NDP C C ON OUTPUT ZI C ARRAY OF DIMENSION NIP WHERE INTERPOLATED C Z VALUES ARE TO BE STORED. C C SPECIAL CONDITIONS INADEQUATE WORK SPACE IWK AND WK MAY C MAY CAUSE INCORRECT RESULTS. C C THE DATA POINTS MUST BE DISTINCT AND THEIR C PROJECTIONS IN THE X-Y PLANE MUST NOT BE C COLLINEAR, OTHERWISE AN ERROR RETURN OCCURS. C ******************************************************************** C C SUBROUTINE IDSFFT (MD,NDP,XD,YD,ZD,NXI,NYI,NZI,XI,YI,ZI,IWK,WK) C C DIMENSION OF XD(NDP), YD(NDP), ZD(NDP), XI(NXI), C ARGUMENTS YI(NYI), ZI(NZI,NYI), WK(6*NDP), C IWK(31*NDP + NXI*NYI) C C PURPOSE THIS SUBROUTINE PERFORMS SMOOTH SURFACE C FITTING WHEN THE PROJECTIONS OF THE DATA C POINTS IN THE X-Y PLANE ARE IRREGULARLY C DISTRIBUTED IN THE PLANE. C C USAGE CALL IDSFFT (MD,NDP,XD,YD,ZD,NXI,NYI,NZI, C XI,YI,ZI,IWK,WK) C C ARGUMENTS C C ON INPUT MD C MODE OF COMPUTATION (MUST BE 1, 2, OR 3, C ELSE AN ERROR RETURN WILL OCCUR). C = 1 IF THIS IS THE FIRST CALL TO THIS C SUBROUTINE, OR IF THE VALUE OF NDP C HAS BEEN CHANGED FROM THE PREVIOUS C CALL, OR IF THE CONTENTS OF THE XD C OR YD ARRAYS HAVE BEEN CHANGED FROM C THE PREVIOUS CALL. C = 2 IF THE VALUES OF NDP AND THE XD, C YD ARRAYS ARE UNCHANGED FROM THE C PREVIOUS CALL, BUT NEW VALUES FOR C XI, YI ARE BEING USED. IF MD = 2 C AND NDP HAS BEEN CHANGED SINCE THE C PREVIOUS CALL TO IDSFFT, AN ERROR C RETURN OCCURS. C = 3 IF THE VALUES OF NDP, NXI, NYI, XD, C YD, XI, YI ARE UNCHANGED FROM THE C PREVIOUS CALL, I.E. IF THE ONLY CHANGE C ON INPUT TO IDSFFT IS IN THE ZD ARRAY. C IF MD = 3 AND NDP, NXI OR NYI HAS BEEN C CHANGED SINCE THE PREVIOUS CALL TO C IDSFFT, AN ERROR RETURN OCCURS. C C BETWEEN THE CALL WITH MD=2 OR MD=3 AND C THE PRECEDING CALL, THE IWK AND WK WORK C ARRAYS SHOULD NOT BE DISTURBED. C C NDP C NUMBER OF DATA POINTS (MUST BE 4 OR C GREATER, ELSE AN ERROR RETURN WILL OCCUR). C C XD C ARRAY OF DIMENSION NDP CONTAINING THE X C COORDINATES OF THE DATA POINTS. C C YD C ARRAY OF DIMENSION NDP CONTAINING THE Y C COORDINATES OF THE DATA POINTS. C C ZD C ARRAY OF DIMENSION NDP CONTAINING THE Z C COORDINATES OF THE DATA POINTS. C C NXI C NUMBER OF OUTPUT GRID POINTS IN THE X- C DIRECTION (MUST BE 1 OR GREATER, ELSE C AN ERROR RETURN WILL OCCUR). C C NYI C NUMBER OF OUTPUT GRID POINTS IN THE Y- C DIRECTION (MUST BE 1 OR GREATER, ELSE C AN ERROR RETURN WILL OCCUR). C C NZI C FIRST DIMENSION OF ZI AS DECLARED IN THE C CALLING PROGRAM. NZI MUST BE GREATER THAN C OR EQUAL TO NXI, ELSE AN ERROR RETURN WILL C OCCUR. C C XI C ARRAY OF DIMENSION NXI CONTAINING THE C X COORDINATES OF THE OUTPUT GRID POINTS. C C YI C ARRAY OF DIMENSION NYI CONTAINING THE C Y COORDINATES OF THE OUTPUT GRID POINTS. C C IWK C INTEGER WORK ARRAY OF DIMENSION AT C LEAST 31*NDP + NXI*NYI C C WK C REAL WORK ARRAY OF DIMENSION AT LEAST 6*NDP C C ON OUTPUT ZI C REAL, TWO-DIMENSIONAL ARRAY OF DIMENSION C (NZI,NYI), STORING THE INTERPOLATED Z C VALUES AT THE OUTPUT GRID POINTS. C C SPECIAL CONDITIONS INADEQUATE WORK SPACE IWK AND WK MAY C MAY CAUSE INCORRECT RESULTS. C C THE DATA POINTS MUST BE DISTINCT AND THEIR C PROJECTIONS IN THE X-Y PLANE MUST NOT BE C COLLINEAR, OTHERWISE AN ERROR RETURN OCCURS. C ******************************************************************** SUBROUTINE IDBVIP(MD,NDP,XD,YD,ZD,NIP,XI,YI,ZI, 1 IWK,WK) C THIS SUBROUTINE CALLS THE IDLCTN, IDPDRV, IDPTIP, AND IDTANG C SUBROUTINES. C DECLARATION STATEMENTS DIMENSION XD(NDP), YD(NDP), ZD(NDP), XI(NIP), 1 YI(NIP), ZI(NIP), IWK(31*NDP + NIP), WK(8*NDP) COMMON/IDLC/ITIPV,DMMY1(13) COMMON/IDPT/ITPV,DMMY(27) C C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR C C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES. C (FOR MD=1,2,3) 10 MD0=MD NDP0=NDP NIP0=NIP C ERROR CHECK. (FOR MD=1,2,3) 20 IF (MD0.LT.1.OR.MD0.GT.3) THEN CALL ULIBER (32, 1' IDBVIP (BIVAR) - INPUT PARAMETER MD OUT OF RANGE',49) STOP 'ULIBER32' ENDIF IF (NDP0.LT.4) THEN CALL ULIBER (33, 1' IDBVIP (BIVAR) - INPUT PARAMETER NDP OUT OF RANGE',50) STOP 'ULIBER33' ENDIF IF (NIP0.LT.1) THEN CALL ULIBER (34, 1' IDBVIP (BIVAR) - INPUT PARAMETER NIP OUT OF RANGE',50) STOP 'ULIBER34' ENDIF IF(MD0.GT.1) GO TO 21 IWK(1)=NDP0 GO TO 22 21 NDPPV=IWK(1) IF (NDP0.NE.NDPPV) THEN CALL ULIBER (50, 1' IDBVIP (BIVAR) - MD=2 OR 3 BUT NDP WAS CHANGED SINCE LAST CALL', 2 63) STOP 'ULIBER50' ENDIF 22 IF(MD0.GT.2) GO TO 23 IWK(3)=NIP GO TO 30 23 NIPPV=IWK(3) IF (NIP0.LT.NIPPV) THEN CALL ULIBER (51, 1' IDBVIP (BIVAR) - MD=3 BUT NIP WAS CHANGED SINCE LAST CALL', 2 58) STOP 'ULIBER51' ENDIF 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 JWIT0=31*NDP0 JWWPD=5*NDP0+1 C TRIANGULATES THE X-Y PLANE. (FOR MD=1) 40 IF(MD0.GT.1) GO TO 41 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 GO TO 50 41 NT=IWK(5) NL=IWK(6) C LOCATES ALL POINTS AT WHICH INTERPOLATION IS TO BE PERFORMED. C (FOR MD=1,2) 50 IF(MD0.GT.2) GO TO 60 ITIPV=0 JWIT=JWIT0 DO 51 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) 51 CONTINUE C ESTIMATES PARTIAL DERIVATIVES AT ALL DATA POINTS. C (FOR MD=1,2,3) 60 CALL IDPDRV(NDP0,XD,YD,ZD,NT,IWK(JWIPT),WK,WK(JWWPD)) C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3) 70 ITPV=0 JWIT=JWIT0 DO 71 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)) 71 CONTINUE RETURN END SUBROUTINE IDGRID(XD,YD,NT,IPT,NL,IPL,NXI,NYI,XI,YI, 1 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(1), YD(1), IPT(1), IPL(1), XI(1), YI(1), NGP(1), 1 IGP(1) C STATEMENT FUNCTIONS SPDT(U1,V1,U2,V2,U3,V3)=(U1-U2)*(U3-U2)+(V1-V2)*(V3-V2) VPDT(U1,V1,U2,V2,U3,V3)=(U1-U3)*(V2-V3)-(V1-V3)*(U2-U3) C PRELIMINARY PROCESSING 10 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. 20 JNGP0=0 JNGP1=2*(NT0+2*NL0)+1 JIGP0=0 JIGP1=NXINYI+1 DO 39 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 22 IXI=1,NXI0 IF(XI(IXI).GE.XMN.AND.XI(IXI).LE.XMX) GO TO 21 IF(INSD.EQ.0) GO TO 22 IXIMX=IXI-1 GO TO 23 21 IF(INSD.EQ.1) GO TO 22 INSD=1 IXIMN=IXI 22 CONTINUE IF(INSD.EQ.0) GO TO 38 IXIMX=NXI0 23 DO 37 IYI=1,NYI0 YII=YI(IYI) IF(YII.LT.YMN.OR.YII.GT.YMX) GO TO 37 DO 36 IXI=IXIMN,IXIMX XII=XI(IXI) L=0 IF(VPDT(X1,Y1,X2,Y2,XII,YII)) 36,25,26 25 L=1 26 IF(VPDT(X2,Y2,X3,Y3,XII,YII)) 36,27,28 27 L=1 28 IF(VPDT(X3,Y3,X1,Y1,XII,YII)) 36,29,30 29 L=1 30 IZI=NXI0*(IYI-1)+IXI IF(L.EQ.1) GO TO 31 NGP0=NGP0+1 JIGP0=JIGP0+1 IGP(JIGP0)=IZI GO TO 36 31 IF(JIGP1.GT.NXINYI) GO TO 33 DO 32 JIGP1I=JIGP1,NXINYI IF(IZI.EQ.IGP(JIGP1I)) GO TO 36 32 CONTINUE 33 NGP1=NGP1+1 JIGP1=JIGP1-1 IGP(JIGP1)=IZI 36 CONTINUE 37 CONTINUE 38 JNGP0=JNGP0+1 NGP(JNGP0)=NGP0 JNGP1=JNGP1-1 NGP(JNGP1)=NGP1 39 CONTINUE C DETERMINES GRID POINTS OUTSIDE THE DATA AREA. C - IN SEMI-INFINITE RECTANGULAR AREA. 40 DO 79 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 42 IXI=1,NXI0 IF(XI(IXI).GE.XMN.AND.XI(IXI).LE.XMX) GO TO 41 IF(INSD.EQ.0) GO TO 42 IXIMX=IXI-1 GO TO 43 41 IF(INSD.EQ.1) GO TO 42 INSD=1 IXIMN=IXI 42 CONTINUE IF(INSD.EQ.0) GO TO 58 IXIMX=NXI0 43 DO 57 IYI=1,NYI0 YII=YI(IYI) IF(YII.LT.YMN.OR.YII.GT.YMX) GO TO 57 DO 56 IXI=IXIMN,IXIMX XII=XI(IXI) L=0 IF(VPDT(X1,Y1,X2,Y2,XII,YII)) 46,45,56 45 L=1 46 IF(SPDT(X2,Y2,X1,Y1,XII,YII)) 56,47,48 47 L=1 48 IF(SPDT(X1,Y1,X2,Y2,XII,YII)) 56,49,50 49 L=1 50 IZI=NXI0*(IYI-1)+IXI IF(L.EQ.1) GO TO 51 NGP0=NGP0+1 JIGP0=JIGP0+1 IGP(JIGP0)=IZI GO TO 56 51 IF(JIGP1.GT.NXINYI) GO TO 53 DO 52 JIGP1I=JIGP1,NXINYI IF(IZI.EQ.IGP(JIGP1I)) GO TO 56 52 CONTINUE 53 NGP1=NGP1+1 JIGP1=JIGP1-1 IGP(JIGP1)=IZI 56 CONTINUE 57 CONTINUE 58 JNGP0=JNGP0+1 NGP(JNGP0)=NGP0 JNGP1=JNGP1-1 NGP(JNGP1)=NGP1 C - IN SEMI-INFINITE TRIANGULAR AREA. 60 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 62 IXI=1,NXI0 IF(XI(IXI).GE.XMN.AND.XI(IXI).LE.XMX) GO TO 61 IF(INSD.EQ.0) GO TO 62 IXIMX=IXI-1 GO TO 63 61 IF(INSD.EQ.1) GO TO 62 INSD=1 IXIMN=IXI 62 CONTINUE IF(INSD.EQ.0) GO TO 78 IXIMX=NXI0 63 DO 77 IYI=1,NYI0 YII=YI(IYI) IF(YII.LT.YMN.OR.YII.GT.YMX) GO TO 77 DO 76 IXI=IXIMN,IXIMX XII=XI(IXI) L=0 IF(SPDT(X1,Y1,X2,Y2,XII,YII)) 66,65,76 65 L=1 66 IF(SPDT(X3,Y3,X2,Y2,XII,YII)) 70,67,76 67 L=1 70 IZI=NXI0*(IYI-1)+IXI IF(L.EQ.1) GO TO 71 NGP0=NGP0+1 JIGP0=JIGP0+1 IGP(JIGP0)=IZI GO TO 76 71 IF(JIGP1.GT.NXINYI) GO TO 73 DO 72 JIGP1I=JIGP1,NXINYI IF(IZI.EQ.IGP(JIGP1I)) GO TO 76 72 CONTINUE 73 NGP1=NGP1+1 JIGP1=JIGP1-1 IGP(JIGP1)=IZI 76 CONTINUE 77 CONTINUE 78 JNGP0=JNGP0+1 NGP(JNGP0)=NGP0 JNGP1=JNGP1-1 NGP(JNGP1)=NGP1 79 CONTINUE RETURN END SUBROUTINE IDLCTN(NDP,XD,YD,NT,IPT,NL,IPL,XII,YII,ITI, 1 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(NDP), YD(NDP), IPT(3*NT), IPL(3*NL), IWK(18*NDP), 1 WK(8*NDP) DIMENSION IDSC(9) COMMON/IDLC/ITIPV,XS1,XS2,YS1,YS2,NTSC(9) C STATEMENT FUNCTIONS SPDT(U1,V1,U2,V2,U3,V3)=(U1-U2)*(U3-U2)+(V1-V2)*(V3-V2) VPDT(U1,V1,U2,V2,U3,V3)=(U1-U3)*(V2-V3)-(V1-V3)*(U2-U3) C PRELIMINARY PROCESSING 10 NDP0=NDP NT0=NT NL0=NL NTL=NT0+NL0 X0=XII Y0=YII C PROCESSING FOR A NEW SET OF DATA POINTS 20 IF(ITIPV.NE.0) GO TO 30 C - DIVIDES THE X-Y PLANE INTO NINE RECTANGULAR SECTIONS. XMN=XD(1) XMX=XMN YMN=YD(1) YMX=YMN DO 21 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) 21 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 22 ISC=1,9 NTSC(ISC)=0 IDSC(ISC)=0 22 CONTINUE IT0T3=0 JWK=0 DO 27 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 23 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 23 IF(YMX.LT.YS1.OR.YMN.GT.YS2) GO TO 24 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 24 IF(YMX.LT.YS2) GO TO 25 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 25 DO 26 ISC=1,9 IF(IDSC(ISC).EQ.0) GO TO 26 JIWK=9*NTSC(ISC)+ISC IWK(JIWK)=IT0 NTSC(ISC)=NTSC(ISC)+1 IDSC(ISC)=0 26 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 27 CONTINUE GO TO 60 C CHECKS IF IN THE SAME TRIANGLE AS PREVIOUS. 30 IT0=ITIPV IF(IT0.GT.NT0) GO TO 40 IT0T3=IT0*3 IP1=IPT(IT0T3-2) X1=XD(IP1) Y1=YD(IP1) IP2=IPT(IT0T3-1) X2=XD(IP2) Y2=YD(IP2) IF(VPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 60 IP3=IPT(IT0T3) X3=XD(IP3) Y3=YD(IP3) IF(VPDT(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 60 IF(VPDT(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 60 GO TO 80 C CHECKS IF ON THE SAME BORDER LINE SEGMENT. 40 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 50 IF(SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 60 IF(SPDT(X2,Y2,X1,Y1,X0,Y0).LT.0.0) GO TO 60 IF(VPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 60 GO TO 80 C CHECKS IF BETWEEN THE SAME TWO BORDER LINE SEGMENTS. 50 IF(SPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 60 IP3=IPL(3*IL2-1) X3=XD(IP3) Y3=YD(IP3) IF(SPDT(X3,Y3,X2,Y2,X0,Y0).LE.0.0) GO TO 80 C LOCATES INSIDE THE DATA AREA. C - DETERMINES THE SECTION IN WHICH THE POINT IN QUESTION LIES. 60 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 70 JIWK=-9+ISC DO 61 ITSC=1,NTSCI JIWK=JIWK+9 IT0=IWK(JIWK) JWK=IT0*4 IF(X0.LT.WK(JWK-3)) GO TO 61 IF(X0.GT.WK(JWK-2)) GO TO 61 IF(Y0.LT.WK(JWK-1)) GO TO 61 IF(Y0.GT.WK(JWK)) GO TO 61 IT0T3=IT0*3 IP1=IPT(IT0T3-2) X1=XD(IP1) Y1=YD(IP1) IP2=IPT(IT0T3-1) X2=XD(IP2) Y2=YD(IP2) IF(VPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 61 IP3=IPT(IT0T3) X3=XD(IP3) Y3=YD(IP3) IF(VPDT(X2,Y2,X3,Y3,X0,Y0).LT.0.0) GO TO 61 IF(VPDT(X3,Y3,X1,Y1,X0,Y0).LT.0.0) GO TO 61 GO TO 80 61 CONTINUE C LOCATES OUTSIDE THE DATA AREA. 70 DO 72 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 72 IF(SPDT(X1,Y1,X2,Y2,X0,Y0).LT.0.0) GO TO 71 IF(VPDT(X1,Y1,X2,Y2,X0,Y0).GT.0.0) GO TO 72 IL2=IL1 GO TO 75 71 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 75 72 CONTINUE IT0=1 GO TO 80 75 IT0=IL1*NTL+IL2 C NORMAL EXIT 80 ITI=IT0 ITIPV=IT0 RETURN END SUBROUTINE IDPDRV(NDP,XD,YD,ZD,NT,IPT,PD,WK) 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 NT = NUMBER OF TRIANGLES, C IPT = INTEGER ARRAY OF DIMENSION 3*NT CONTAINING THE C POINT NUMBERS OF THE VERTEXES OF THE TRIANGLES. 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 ITH C DATA POINT ARE TO BE STORED AS THE (5*I-4)TH, C (5*I-3)RD, (5*I-2)ND, (5*I-1)ST AND (5*I)TH C ELEMENTS, RESPECTIVELY, WHERE I = 1, 2, ..., C NDP. C THE OTHER PARAMETER IS C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A C WORK AREA. C DECLARATION STATEMENTS DIMENSION XD(NDP), YD(NDP), ZD(NDP), IPT(3*NT), PD(5*NDP), WK(NDP) DIMENSION IPTI(3),XV(3),YV(3),ZV(3),ZXV(3),ZYV(3), 1 W1(3),W2(3) DATA EPSLN/1.0E-6/ C PRELIMINARY PROCESSING 10 NDP0=NDP NT0=NT C CLEARS THE PD ARRAY. 20 JPDMX=5*NDP0 DO 21 JPD=1,JPDMX PD(JPD)=0.0 21 CONTINUE DO 22 IDP=1,NDP WK(IDP)=0.0 22 CONTINUE C ESTIMATES ZX AND ZY. 30 DO 34 IT=1,NT0 JPT0=3*(IT-1) DO 31 IV=1,3 JPT=JPT0+IV IDP=IPT(JPT) IPTI(IV)=IDP XV(IV)=XD(IDP) YV(IV)=YD(IDP) ZV(IV)=ZD(IDP) 31 CONTINUE DX1=XV(2)-XV(1) DY1=YV(2)-YV(1) DZ1=ZV(2)-ZV(1) DX2=XV(3)-XV(1) DY2=YV(3)-YV(1) DZ2=ZV(3)-ZV(1) VPX=DY1*DZ2-DZ1*DY2 VPY=DZ1*DX2-DX1*DZ2 VPZ=DX1*DY2-DY1*DX2 VPZMN=ABS(DX1*DX2+DY1*DY2)*EPSLN IF(ABS(VPZ).LE.VPZMN) GO TO 34 D12=SQRT((XV(2)-XV(1))**2+(YV(2)-YV(1))**2) D23=SQRT((XV(3)-XV(2))**2+(YV(3)-YV(2))**2) D31=SQRT((XV(1)-XV(3))**2+(YV(1)-YV(3))**2) W1(1)=1.0/(D31*D12) W1(2)=1.0/(D12*D23) W1(3)=1.0/(D23*D31) W2(1)=VPZ*W1(1) W2(2)=VPZ*W1(2) W2(3)=VPZ*W1(3) 32 DO 33 IV=1,3 IDP=IPTI(IV) JPD0=5*(IDP-1) WI=(W1(IV)**2)*W2(IV) PD(JPD0+1)=PD(JPD0+1)+VPX*WI PD(JPD0+2)=PD(JPD0+2)+VPY*WI WK(IDP)=WK(IDP)+VPZ*WI 33 CONTINUE 34 CONTINUE DO 36 IDP=1,NDP0 JPD0=5*(IDP-1) PD(JPD0+1)=-PD(JPD0+1)/WK(IDP) PD(JPD0+2)=-PD(JPD0+2)/WK(IDP) 36 CONTINUE C ESTIMATES ZXX, ZXY, AND ZYY. 40 DO 44 IT=1,NT0 JPT0=3*(IT-1) DO 41 IV=1,3 JPT=JPT0+IV IDP=IPT(JPT) IPTI(IV)=IDP XV(IV)=XD(IDP) YV(IV)=YD(IDP) JPD0=5*(IDP-1) ZXV(IV)=PD(JPD0+1) ZYV(IV)=PD(JPD0+2) 41 CONTINUE DX1=XV(2)-XV(1) DY1=YV(2)-YV(1) DZX1=ZXV(2)-ZXV(1) DZY1=ZYV(2)-ZYV(1) DX2=XV(3)-XV(1) DY2=YV(3)-YV(1) DZX2=ZXV(3)-ZXV(1) DZY2=ZYV(3)-ZYV(1) VPXX=DY1*DZX2-DZX1*DY2 VPXY=DZX1*DX2-DX1*DZX2 VPYX=DY1*DZY2-DZY1*DY2 VPYY=DZY1*DX2-DX1*DZY2 VPZ=DX1*DY2-DY1*DX2 VPZMN=ABS(DX1*DX2+DY1*DY2)*EPSLN IF(ABS(VPZ).LE.VPZMN) GO TO 44 D12=SQRT((XV(2)-XV(1))**2+(YV(2)-YV(1))**2) D23=SQRT((XV(3)-XV(2))**2+(YV(3)-YV(2))**2) D31=SQRT((XV(1)-XV(3))**2+(YV(1)-YV(3))**2) W1(1)=1.0/(D31*D12) W1(2)=1.0/(D12*D23) W1(3)=1.0/(D23*D31) W2(1)=VPZ*W1(1) W2(2)=VPZ*W1(2) W2(3)=VPZ*W1(3) 42 DO 43 IV=1,3 IDP=IPTI(IV) JPD0=5*(IDP-1) WI=(W1(IV)**2)*W2(IV) PD(JPD0+3)=PD(JPD0+3)+VPXX*WI PD(JPD0+4)=PD(JPD0+4)+(VPXY+VPYX)*WI PD(JPD0+5)=PD(JPD0+5)+VPYY*WI 43 CONTINUE 44 CONTINUE DO 46 IDP=1,NDP0 JPD0=5*(IDP-1) PD(JPD0+3)=-PD(JPD0+3)/WK(IDP) PD(JPD0+4)=-PD(JPD0+4)/(2.0*WK(IDP)) PD(JPD0+5)=-PD(JPD0+5)/WK(IDP) 46 CONTINUE RETURN END SUBROUTINE IDPTIP(XD,YD,ZD,NT,IPT,NL,IPL,PDD,ITI,XII,YII, 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(1), YD(1), ZD(1), IPT(1), IPL(1), PDD(1) COMMON/IDPT/ITPV,X0,Y0,AP,BP,CP,DP, 1 P00,P10,P20,P30,P40,P50,P01,P11,P21,P31,P41, 2 P02,P12,P22,P32,P03,P13,P23,P04,P14,P05 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) X0=XD(IDP) Y0=YD(IDP) Z0=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=Z0 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-X0 V=YII-Y0 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,NDP,XD,YD,ZD,NXI,NYI,NZI,XI,YI,ZI, 1 IWK,WK) C THIS SUBROUTINE CALLS THE IDGRID, IDPDRV, IDPTIP, AND IDTANG C SUBROUTINES. C DECLARATION STATEMENTS DIMENSION XD(NDP), YD(NDP), ZD(NDP), XI(NXI), 1 YI(NYI), ZI(NZI,NYI), IWK(31*NDP + NXI*NYI), WK(6*NDP) COMMON/IDPT/ITPV,DMMY(27) C C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR C C SETTING OF SOME INPUT PARAMETERS TO LOCAL VARIABLES. C (FOR MD=1,2,3) 10 MD0=MD NDP0=NDP NXI0=NXI NYI0=NYI C ERROR CHECK. (FOR MD=1,2,3) 20 IF (MD0.LT.1.OR.MD0.GT.3) THEN CALL ULIBER (39, 1' IDSFFT (BIVAR) - INPUT PARAMETER MD OUT OF RANGE',49) STOP 'ULIBER39' ENDIF IF (NDP0.LT.4) THEN CALL ULIBER (40, 1' IDSFFT (BIVAR) - INPUT PARAMETER NDP OUT OF RANGE',50) STOP 'ULIBER40' ENDIF IF (NXI0.LT.1.OR.NYI0.LT.1) THEN CALL ULIBER (41, 1' IDSFFT (BIVAR) - INPUT PARAMETER NXI OR NYI OUT OF RANGE',57) STOP 'ULIBER41' ENDIF IF (NXI0.GT.NZI) THEN CALL ULIBER (42, 1' IDSFFT (BIVAR) - INPUT PARAMETER NZI IS LESS THAN NXI',54) STOP 'ULIBER42' ENDIF IF(MD0.GT.1) GO TO 21 IWK(1)=NDP0 GO TO 22 21 NDPPV=IWK(1) IF (NDP0.NE.NDPPV) THEN CALL ULIBER (43, 1' IDSFFT (BIVAR) - MD=2 OR 3 BUT NDP WAS CHANGED SINCE LAST CALL', 2 63) STOP 'ULIBER43' ENDIF 22 IF(MD0.GT.2) GO TO 23 IWK(3)=NXI0 IWK(4)=NYI0 GO TO 30 23 NXIPV=IWK(3) NYIPV=IWK(4) IF (NXI0.NE.NXIPV) THEN CALL ULIBER (45, 1' IDSFFT (BIVAR) - MD=3 BUT NXI WAS CHANGED SINCE LAST CALL', 2 58) STOP 'ULIBER45' ENDIF IF (NYI0.NE.NYIPV) THEN CALL ULIBER (46, 1' IDSFFT (BIVAR) - MD=3 BUT NYI WAS CHANGED SINCE LAST CALL', 2 58) STOP 'ULIBER46' ENDIF 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 JWIGP0=31*NDP0 JWWPD=5*NDP0+1 C TRIANGULATES THE X-Y PLANE. (FOR MD=1) 40 IF(MD0.GT.1) GO TO 41 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 GO TO 50 41 NT=IWK(5) NL=IWK(6) C SORTS OUTPUT GRID POINTS IN ASCENDING ORDER OF THE TRIANGLE C NUMBER AND THE BORDER LINE SEGMENT NUMBER. (FOR MD=1,2) 50 IF(MD0.GT.2) GO TO 60 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) 60 CALL IDPDRV(NDP0,XD,YD,ZD,NT,IWK(JWIPT),WK,WK(JWWPD)) C INTERPOLATES THE ZI VALUES. (FOR MD=1,2,3) 70 ITPV=0 JIG0MX=0 JIG1MN=NXI0*NYI0+1 NNGP=NT+2*NL DO 79 JNGP=1,NNGP ITI=JNGP IF(JNGP.LE.NT) GO TO 71 IL1=(JNGP-NT+1)/2 IL2=(JNGP-NT+2)/2 IF(IL2.GT.NL) IL2=1 ITI=IL1*(NT+NL)+IL2 71 JWNGP=JWNGP0+JNGP NGP0=IWK(JWNGP) IF(NGP0.EQ.0) GO TO 76 JIG0MN=JIG0MX+1 JIG0MX=JIG0MX+NGP0 DO 72 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(IXI,IYI)) 72 CONTINUE 76 JWNGP=JWNGP0+2*NNGP+1-JNGP NGP1=IWK(JWNGP) IF(NGP1.EQ.0) GO TO 79 JIG1MX=JIG1MN-1 JIG1MN=JIG1MN-NGP1 DO 77 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(IXI,IYI)) 77 CONTINUE 79 CONTINUE RETURN END SUBROUTINE IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK) 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 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(NDP), YD(NDP), IPT(6*NDP - 15), IPL(6*NDP), 1 IWL(18*NDP), IWP(NDP), WK(NDP) DIMENSION ITF(2) DATA EPSLN/1.0E-6/, NREP/100/ C STATEMENT FUNCTIONS DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2 SPDT(U1,V1,U2,V2,U3,V3)=(U2-U1)*(U3-U1)+(V2-V1)*(V3-V1) VPDT(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) THEN CALL ULIBER (47, 1' IDTANG (BIVAR) - INPUT PARAMETER NDP OUT OF RANGE',50) STOP 'ULIBER47' ENDIF 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.) THEN CALL ULIBER (48, 1' IDTANG (BIVAR) - TWO OF THE INPUT DATA POINTS ARE IDENTICAL',60) STOP 'ULIBER48' ENDIF IF(DSQI.GE.DSQMN) GO TO 21 DSQMN=DSQI IPMN1=IP1 IPMN2=IP2 21 CONTINUE 22 CONTINUE 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 X1=XD(IPMN1) Y1=YD(IPMN1) X2=XD(IPMN2) Y2=YD(IPMN2) DO 36 JP=3,NDP0 IP=IWP(JP) SP=SPDT(XD(IP),YD(IP),X1,Y1,X2,Y2) VP=VPDT(XD(IP),YD(IP),X1,Y1,X2,Y2) IF(ABS(VP).GT.(ABS(SP)*EPSLN)) GO TO 37 36 CONTINUE CALL ULIBER (49,' IDTANG (BIVAR) - ALL COLLINEAR DATA POINTS',43) STOP 'ULIBER49' 37 IF(JP.EQ.3) GO TO 40 JPMX=JP DO 38 JPC=4,JPMX JP=JPMX+4-JPC 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(VPDT(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 FIRST INVISIBLE AND VISIBLE BORDER LINE SEG- C - MENTS, ILIV AND ILVS. DO 53 IL=1,NL0 IP2=IPL(3*IL-2) IP3=IPL(3*IL-1) X2=XD(IP2) Y2=YD(IP2) X3=XD(IP3) Y3=YD(IP3) SP=SPDT(X1,Y1,X2,Y2,X3,Y3) VP=VPDT(X1,Y1,X2,Y2,X3,Y3) IF(IL.NE.1) GO TO 51 IXVS=0 IF(VP.LE.(ABS(SP)*(-EPSLN))) IXVS=1 ILIV=1 ILVS=1 GO TO 53 51 IXVSPV=IXVS IF(VP.GT.(ABS(SP)*(-EPSLN))) GO TO 52 IXVS=1 IF(IXVSPV.EQ.1) GO TO 53 ILVS=IL IF(ILIV.NE.1) GO TO 54 GO TO 53 52 IXVS=0 IF(IXVSPV.EQ.0) GO TO 53 ILIV=IL IF(ILVS.NE.1) GO TO 54 53 CONTINUE IF(ILIV.EQ.1.AND.ILVS.EQ.1) ILVS=NL0 54 IF(ILVS.LT.ILIV) ILVS=ILVS+NL0 C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY. 55 IF(ILIV.EQ.1) GO TO 60 NLSH=ILIV-1 NLSHT3=NLSH*3 DO 56 JL1=1,NLSHT3 JL2=JL1+NLT3 IPL(JL2)=IPL(JL1) 56 CONTINUE DO 57 JL1=1,NLT3 JL2=JL1+NLSHT3 IPL(JL1)=IPL(JL2) 57 CONTINUE ILVS=ILVS-NLSH 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 IL=ILVS,NL0 ILT3=IL*3 IPL1=IPL(ILT3-2) IPL2=IPL(ILT3-1) IT =IPL(ILT3) 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(IL.NE.ILVS) GO TO 61 IPL(ILT3-1)=IP1 IPL(ILT3) =NT0 61 IF(IL.NE.NL0) GO TO 62 NLN=ILVS+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(IL.EQ.ILVS) IPL(ILT3)=IT IF(IL.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 IPL1=IWL(2*ILF-1) IPL2=IWL(2*ILF) 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. JWL1MN=2*NLFC+1 NLFT2=NLF*2 DO 77 JWL1=JWL1MN,NLFT2 JWL=JWL1+1-JWL1MN 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(VPDT(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 END FUNCTION IDXCHG(X,Y,I1,I2,I3,I4) 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(1), Y(1) EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ), 1 (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ) DATA EPSLN/1.0E-6/ 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(S3SQ,S4SQ)-AMIN1(S1SQ,S2SQ)).GT.EPSLN) 1 IDX=1 30 IDXCHG=IDX RETURN END SUBROUTINE ULIBER (IERR,MESS,LMESS) C SUBROUTINE ULIBER (IERR,MESS,LMESS) C C PURPOSE TO PRINT AN ERROR NUMBER AND AN ERROR MESSAGE C OR JUST AN ERROR MESSAGE. C C USAGE CALL ULIBER (IERR,MESS,LMESS) C C ARGUMENTS C ON INPUT IERR C THE ERROR NUMBER (PRINTED ONLY IF NON-ZERO). C C MESS C MESSAGE TO BE PRINTED. C C LMESS C NUMBER OF CHARACTERS IN MESS (.LE. 130). C C ARGUMENTS C ON OUTPUT NONE C C I/O THE MESSAGE IS WRITEN TO UNIT 101. C ****************************************************************** C REAL MESS(1) C IF (IERR.NE.0) WRITE (101,1001) IERR NWORDS=(LMESS+7)/8 WRITE (101,1002) (MESS(I),I=1,NWORDS) RETURN C 1001 FORMAT (6H0IERR=,I5) 1002 FORMAT (16A8,A2) C END