[Scipy-svn] r4637 - in branches/Interpolate1D: . docs extensions tests

scipy-svn@scip... scipy-svn@scip...
Mon Aug 11 16:39:20 CDT 2008


Author: fcady
Date: 2008-08-11 16:39:17 -0500 (Mon, 11 Aug 2008)
New Revision: 4637

Added:
   branches/Interpolate1D/extensions/interp_526.f
   branches/Interpolate1D/extensions/interp_526a.pyf
Modified:
   branches/Interpolate1D/__init__.py
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/info.py
   branches/Interpolate1D/interpolate2d.py
   branches/Interpolate1D/setup.py
   branches/Interpolate1D/tests/regression_test.py
   branches/Interpolate1D/tests/test_interpolate2d.py
Log:
incorporation of TOMS algorithm 526

Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/__init__.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -15,3 +15,6 @@
 # wrapped by interpolate*.py files
 from fitpack_wrapper import Spline
 
+# wrapper around Fortran implementing Tom's algorithm 526
+from algorithm526_wrapper import algorithm526
+

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-08-11 21:39:17 UTC (rev 4637)
@@ -564,6 +564,8 @@
     'linear', 'cubic' and 'spline', however, work in both cases, and we try to give analogous
     methods the same name.  But some methods are particular to, or have only been written
     for, one praticular dimensionality.
+#. In particular, 2D supports the keywork '526', which implements TOMS algorithm 526.
+    See below for more information.
 
 As in 1D, linear interpolation is used by default, while out of bounds returns NaN.
 
@@ -639,6 +641,13 @@
 #) get_coeffs
     Same as Spline
     
+-------------------------------
+TOMS Algorithm 526
+-------------------------------
+TOMS (Transactions on Mathematical Software) algorithm 526 is an
+algorithm for interpolation of from scattered 2-dimensional data.  It is
+described in "ALGORITHM 526: Bivariate Interpolation and Smooth
+Fitting for Irregularly Distributed Data Points [El ]" by Hiroshi Akima. 
 
 
 ================================================

Added: branches/Interpolate1D/extensions/interp_526.f
===================================================================
--- branches/Interpolate1D/extensions/interp_526.f	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/extensions/interp_526.f	2008-08-11 21:39:17 UTC (rev 4637)
@@ -0,0 +1,1681 @@
+CF77FLAGS(gnu)=-fno-automatic -finit-local-zero
+      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
+C        FIXME: This is where the thing dies if we have more than 1 output value
+        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      WRITE(*,*) 'IDLCTN 0.5', NIT
+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      WRITE(*,*) 'IDLCTN 1'
+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
+C      WRITE(*,*) 'IDLCTN after if90'
+      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)
+C      WRITE(*,*) 'XXXXXXXXXXXXXXXXXXXXIDLCTN after 90', IP1
+      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

Added: branches/Interpolate1D/extensions/interp_526a.pyf
===================================================================
--- branches/Interpolate1D/extensions/interp_526a.pyf	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/extensions/interp_526a.pyf	2008-08-11 21:39:17 UTC (rev 4637)
@@ -0,0 +1,25 @@
+!    -*- f90 -*-
+! Note: the context of this file is case sensitive.
+
+python module interp_526 ! in 
+    ! blubber foo bar
+    interface  ! in :interp_526
+        subroutine idbvip(md,ncp,ndp,xd,yd,zd,nip,xi,yi,zi,iwk,wk) ! in :interp_526:interp_526a.f
+            integer optional,check((md>=1)&&(md<=3)):: md = 1
+            integer optional,check((ncp>=2)&&(ncp<ndp)),depend(ndp):: ncp = (ndp-1<4?ndp-1:4)
+            integer depend(xd),intent(hide),check(ndp>=4) :: ndp=len(xd)
+            real dimension(ndp) :: xd
+            real dimension(ndp),depend(ndp),check(len(yd)==ndp) :: yd
+            real dimension(ndp),depend(ndp) :: zd
+            integer depend(xi),intent(hide) :: nip=len(xi)
+            real dimension(nip) :: xi
+            real dimension(nip),depend(nip),check(len(yi)==nip) :: yi
+            real dimension(nip),depend(nip),intent(out) :: zi
+            integer dimension((MAX(31,27+ncp))*ndp+nip),depend(ncp,ndp,nip),intent(hide,cache) :: iwk
+            real dimension(8 * ndp),depend(ndp),intent(hide,cache) :: wk
+        end subroutine idbvip
+    end interface 
+end python module interp_526
+
+! This file was auto-generated with f2py (version:2_2527).
+! See http://cens.ioc.ee/projects/f2py2e/

Modified: branches/Interpolate1D/info.py
===================================================================
--- branches/Interpolate1D/info.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/info.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -56,6 +56,7 @@
         logarithmic :  logarithmic interpolation
         block : block interpolation
         block_average_above : block average above interpolation
+        algorithm526 : interpolation of 2D scattered data by TOMS Algorithm 526
 
 """
         

Modified: branches/Interpolate1D/interpolate2d.py
===================================================================
--- branches/Interpolate1D/interpolate2d.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/interpolate2d.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -2,6 +2,7 @@
 from numpy import NaN, array
 import numpy as np
 from fitpack_wrapper import Spline2d
+from algorithm526_wrapper import algorithm526
 
 def atleast_1d_and_contiguous(ary, dtype = np.float64):
     # FIXME : don't have in 2 places
@@ -19,7 +20,8 @@
                     'Quartic' : Spline2d(kx=4, ky=4), 'quartic' : Spline2d(kx=4, ky=4),
                     'Quar' : Spline2d(kx=4, ky=4), 'quar' : Spline2d(kx=4, ky=4),
                     'Quintic' : Spline2d(kx=5, ky=5), 'quintic' : Spline2d(kx=5, ky=5),
-                    'Quin' : Spline2d(kx=5, ky=5), 'quin' : Spline2d(kx=5, ky=5)
+                    'Quin' : Spline2d(kx=5, ky=5), 'quin' : Spline2d(kx=5, ky=5),
+                    '526' : algorithm526, 'algorithm526':algorithm526,
                 }
                 
 # dictionary of types for casting.  key = possible datatype, value = datatype it is cast to
@@ -133,11 +135,12 @@
         self._init_xyz(x, y, z, bad_data)
         
         self.kind = self._init_interp_method(kind)
+        
         self.out = self._init_interp_method(out)
         
     def _init_xyz(self, x, y, z, bad_data):
-        # all must be vectors.  Should allow meshgrid, though?  Perhaps.
-        # Mention that in documentation.
+
+        # FIXME : perhaps allow 2D input if it is inthe form of meshgrid
         
         if bad_data is not None:
             try: # check that bad_data contains only numerical values
@@ -248,9 +251,9 @@
                                 (min(self._y) <= newy) & (newy <= max(self._y))        
         
         # filling array of interpolated z-values
-        result = np.zeros(np.shape(newx))
-        if sum(in_range_mask) > 0: # if there are in-range values.  hack to deal 
-                                                # with behavior of np.vectorize on arrays of length 0
+        result = np.zeros(np.shape(newx), dtype = self._zdtype)
+        if sum(in_range_mask) > 0:  # if there are in-range values.  hack to deal
+                                               # with behavior of np.vectorize on arrays of length 0
             result[in_range_mask] = self.kind(newx[in_range_mask], newy[in_range_mask])        
         if sum(~in_range_mask) > 0:
             result[~in_range_mask] = self.out(newx[~in_range_mask], newy[~in_range_mask])

Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/setup.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -30,15 +30,21 @@
                         
     # ND Image routines for ND interpolation
     config.add_extension('_nd_image',
-                        sources=["extensions/ndimage/nd_image.c",
-                                    "extensions/ndimage/ni_filters.c",
-                                    "extensions/ndimage/ni_fourier.c",
-                                    "extensions/ndimage/ni_interpolation.c",
-                                    "extensions/ndimage/ni_measure.c",
-                                    "extensions/ndimage/ni_morphology.c",
-                                    "extensions/ndimage/ni_support.c"],
+                        sources = ["extensions/ndimage/nd_image.c",
+                                        "extensions/ndimage/ni_filters.c",
+                                        "extensions/ndimage/ni_fourier.c",
+                                        "extensions/ndimage/ni_interpolation.c",
+                                        "extensions/ndimage/ni_measure.c",
+                                        "extensions/ndimage/ni_morphology.c",
+                                        "extensions/ndimage/ni_support.c"],
                         include_dirs=['extensions/ndimage']+[get_include()],
                         )
+    
+    # implements algorithm 526 for 2D interpolation
+    config.add_extension('interp_526',
+                        sources = ['extensions/interp_526a.pyf',
+                                        'extensions/interp_526.f']
+                        )
                         
     config.add_data_dir('docs')
 

Modified: branches/Interpolate1D/tests/regression_test.py
===================================================================
--- branches/Interpolate1D/tests/regression_test.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/tests/regression_test.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -38,7 +38,8 @@
         Test1_dict[test_name] = t2-t1
 
     current_dict['Test1'] = Test1_dict
-    
+
+# 2-dimensional
 if True:
     test_list = [name for name in dir(Test2) if name.find('test_') == 0]
     Test2_dict = {}
@@ -52,6 +53,7 @@
 
     current_dict['Test2'] = Test2_dict
 
+# N-dmensional
 if True:
     test_list = [name for name in dir(TestN) if name.find('test_') == 0]
     TestN_dict = {}

Modified: branches/Interpolate1D/tests/test_interpolate2d.py
===================================================================
--- branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-11 16:17:36 UTC (rev 4636)
+++ branches/Interpolate1D/tests/test_interpolate2d.py	2008-08-11 21:39:17 UTC (rev 4637)
@@ -129,7 +129,7 @@
     def test_string_linear(self):
         """ make sure : string 'linear' works
         """
-        N = 7
+        N = 700
         X, Y = meshgrid(arange(N), arange(N))
         Z = X + Y
         x, y, z = map(ravel, [X, Y, Z] )
@@ -145,7 +145,7 @@
     def test_string_quadratic(self):
         """ make sure : string 'quadratic' works
         """
-        N = 7
+        N = 700
         X, Y = meshgrid(arange(N), arange(N))
         Z = X + Y
         x, y, z = map(ravel, [X, Y, Z] )
@@ -161,7 +161,7 @@
     def test_string_cubic(self):
         """make sure : string "cubic" works
         """
-        N = 7
+        N = 700
         X, Y = meshgrid(arange(N), arange(N))
         Z = X + Y
         x, y, z = map(ravel, [X, Y, Z] )
@@ -174,5 +174,22 @@
         
         self.assertAllclose(newz, newx+newy)
         
+    def test_string_526(self):
+        """ make sure : keyword '526' works
+            ie that TOMS algorithm 526 works
+        """
+        N = 700
+        X, Y = meshgrid(arange(N), arange(N))
+        Z = X + Y
+        x, y, z = map(ravel, [X, Y, Z] )
+        
+        newx = np.arange(1, N)-0.5
+        newy = newx
+        
+        interp_func = Interpolate2d(x, y, z, kind='526', out='526')
+        newz = interp_func(newx, newy)
+        
+        self.assertAllclose(newz, newx+newy)
+        
 if __name__ == '__main__':
     unittest.main()                 
\ No newline at end of file



More information about the Scipy-svn mailing list