LOTPS
SUBROUTINE LOTPS (MODE,NPPR,NPI,XI,YI,FI,NXO,XO,NYO,YO,IWK,NIWK,
1 NIWKU,WK,NWK,NWKU,FO,KER)
C***START PROLOGUE LOTPS
C
C THIS VERSION IS DATED 03/04/82.
C
C RICHARD FRANKE
C DEPARTMENT OF MATHEMATICS
C NAVAL POSTGRADUATE SCHOOL
C MONTEREY, CALIFORNIA 93940
C (408)646-2758 / 2206
C
C
C
C REFERENCE
C SMOOTH INTERPOLATION OF SCATTERED DATA BY LOCAL THIN
C PLATE SPLINES, COMPUTERS AND MATHEMATICS WITH
C APPLICATIONS 8(1982)???-???+8
C OR
C NAVAL POSTGRADUATE SCHOOL TR#NPS-53-81-002, 1981
C (AVAILABLE FROM NTIS, AD-A098 232/2)
C
C ABSTRACT
C SUBROUTINE LOTPS SERVES AS THE USER INTERFACE FOR A SET OF
C SUBROUTINES WHICH SOLVE THE SCATTERED DATA INTERPOLATION
C PROBLEM. A SMOOTH FUNCTION PASSING THROUGH THE GIVEN POINTS
C (XI(K),YI(K),FI(K)),K=1,...,NPI IS CONSTRUCTED.
C THE RESULT RETURNED IS AN ARRAY OF VALUES, FO(I,J), OF THE INT-
C ERPOLATION FUNCTION AT GRID POINTS, (XO(I),YO(J)),I=1,...,NXO,
C J=1,...,NYO.
C THE METHOD USED INVOLVES CONSTRUCTION OF LOCALLY DEFINED 'THIN
C PLATE SPLINES', WHICH ARE THEN BLENDED TOGETHER SMOOTHLY
C THROUGH THE USE OF A PARTITION OF UNITY DEFINED ON A
C RECTANGULAR GRID ON THE PLANE. THE FUNCTIONS IN THE PARTITION
C OF UNITY ARE UNIVARIATE PIECEWISE HERMITE CUBIC POLYNOMIALS.
C
C CAUTIONS
C THE USER SHOULD BE AWARE THAT FOR SOME DATA THE INTERPOLATION
C FUNCTION MAY BE ILL-BEHAVED. SOME INVESTIGATION OF ITS
C BEHAVIOR FOR THE TYPE OF DATA TO BE INPUT SHOULD BE UNDERTAKEN
C BEFORE IMBEDDING ANY SCHEME FOR SCATTERED DATA INTERPOLATION
C INTO ANOTHER PROGRAM.
C
C DESCRIPTION OF ARGUMENTS
C
C MODE - INPUT. INDICATES THE STATUS OF THE CALCULATION.
C = 1, SET UP THE PROBLEM. COMPUTE THE COEFFICIENTS
C FOR THE LOCAL APPROXIMATIONS BY THIN PLATE
C SPLINES, AND RETURN THE GRID OF INTERPOLATED
C FUNCTION VALUES INDICATED BY NXO, XO, NYO, YO
C IN THE ARRAY FO.
C = 2, THIS MODE VALUE IS A CONVENIENCE FOR USERS WHO
C WISH TO CALL THE ROUTINE TO EVALUATE THE
C SURFACE REPEATEDLY ON DIFFERENT GRIDS OF
C POINTS. A CALL TO LOTPS WITH MODE = 1 HAS
C BEEN MADE PREVIOUSLY, NOW CALCULATE
C THE GRID OF INTERPOLATED POINTS INDICATED
C BY NXO, XO, NYO, YO IN IN THE ARRAY FO. THE
C PROGRAM ASSUMES THAT THE ARRAYS XI, YI, IWK,
C AND WK ARE UNCHANGED FROM THE PREVIOUS CALL.
C NPPR - INPUT. DESIRED AVERAGE NUMBER OF POINTS PER REGION.
C THE SUGGESTED VALUE FOR THE NOVICE USER IS TEN,
C WHICH USUALLY GIVES GOOD RESULTS. THIS PAR-
C AMETER HAS TO DO WITH THE LOCAL PROPERTY OF THE
C SURFACE. THE INFLUENCE REGION OF A POINT HAS
C AREA WHICH IS ROUGHLY PROPORTIONAL TO NPPR.
C UNDER CERTAIN CONDITIONS, SUCH AS TO PRESERVE
C ROTATIONAL INVARIANCE, OR TO FORCE CERTAIN
C SETS OF POINTS TO BELONG TO THE SAME REGION,
C THE USER MAY SPECIFY HIS OWN GRID LINES.
C IF THE USER WISHES TO SPECIFY HIS OWN GRID LINES
C X TILDA AND Y TILDA, HE MAY DO SO BY SETTING
C NPPR = 0 AND SETTING NECESSARY VALUES IN THE
C ARRAYS IWK AND WK, AS NOTED BELOW. DATA WHICH
C HAS A POOR DISTRIBUTION OVER THE REGION OF INT-
C EREST SHOULD PROBABLY HAVE THE GRID SPECIFIED.
C THIS IS ALSO ADVISABLE IF THE X-Y POINTS OCCUR
C ALONG LINES. SEE THE REFERENCE FOR ADDITIONAL
C DETAILS.
C NPI - INPUT. NUMBER OF INPUT DATA POINTS.
C XI - \
C YI - INPUT ARRAYS. THE DATA POINTS (XI,YI,FI), I=1,...,NPI.
C FI - /
C NXO - INPUT. THE NUMBER OF XO VALUES AT WHICH THE INTERP-
C OLATION FUNCTION IS TO BE CALCULATED.
C XO - INPUT ARRAY. THE VALUES OF X AT WHICH THE INTERPOLATION
C FUNCTION IS TO BE CALCULATED. THESE SHOULD
C BE IN INCREASING ORDER FOR MOST EFFICIENT
C EVALUATION, HOWEVER, THEY ONLY NEED TO BE
C MONOTONIC.
C NYO - INPUT. THE NUMBER OF YO VALUES AT WHICH THE INTERP-
C OLATION FUNCTION IS TO BE CALCULATED.
C YO - INPUT ARRAY. THE VALUES OF Y AT WHICH THE INTERPOLATION
C FUNCTION IS TO BE CALCULATED. THESE SHOULD
C BE IN INCREASING ORDER FOR MOST EFFICIENT
C EVALUATION, HOWEVER, THEY ONLY NEED TO BE
C MONOTONIC.
C IWK - INPUT/OUTPUT ARRAY. THIS ARRAY IS OUTPUT WHEN MODE = 1
C AND IS INPUT WHEN MODE = 2. THIS MUST BE
C AN ARRAY DIMENSIONED APPROXIMATELY 7*NPI. THE
C EXACT DIMENSION IS NOT KNOWN A PRIORI, BUT
C WILL BE RETURNED AS THE VALUE OF NIWKU.
C WHEN NPPR IS INPUT AS ZERO THE USER MUST
C SPECIFY THE NUMBER OF VERTICAL GRID LINES (THE
C NUMBER OF X TILDA VALUES) IN IWK(1) AND THE
C NUMBER OF HORIZONTAL GRID LINES (THE NUMBER OF
C Y TILDA VALUES) IN IWK(2).
C NIWK - INPUT. ON ENTRY WITH MODE = 1 THIS MUST BE SET TO THE
C DIMENSION OF THE ARRAY IWK IN THE CALLING
C PROGRAM.
C NIWKU- OUTPUT. THE ACTUAL NUMBER OF LOCATIONS NEEDED IN THE
C ARRAY IWK.
C WK - INPUT/OUTPUT ARRAY. THIS ARRAY IS OUTPUT WHEN MODE = 1
C AND IS INPUT WHEN MODE = 2. THIS MUST BE AN
C ARRAY DIMENSIONED APPROXIMATELY 7*NPI PLUS
C THE NUMBER NEEDED TO SET UP AND SOLVE THE SYSTEM
C OF EQUATIONS FOR THE LOCAL APPROXIMATIONS. FOR
C NPPR NONZERO THIS WILL BE ABOUT 2.5*NPPR*NPPR
C PLUS 11*NPPR. THE EXACT DIMENSION IS NOT KNOWN
C A PRIORI, BUT WILL BE RETURNED AS THE VALUE OF
C NWKU.
C WHEN NPPR IS INPUT AS ZERO THE USER MUST SPECIFY
C THE VALUES OF X TILDA AND Y TILDA AS FOLLOWS.
C WK(2), ... , WK(NXG+1) ARE THE NXG (= IWK(1))
C X GRID VALUES, X(I) TILDA, IN INCREASING ORDER.
C TYPICALLY WK(1) = MIN X(I), ALTHOUGH IT NEED
C NOT BE. WK(1) MUST BE LESS THAN OR EQUAL TO
C WK(2), AND SHOULD BE LESS THAN OR EQUAL TO
C MIN X(I). WK(NXG+2) IS USUALLY MAX X(I), AL-
C THOUGH IT NEED NOT BE. WK(NXG+2) MUST BE
C GREATER THAN WK(NXG+1), AND SHOULD BE GREATER
C THAN OR EQUAL TO MAX X(I).
C THE VALUES OF WK(NXG+3), ... , WK(NXG+NYG+4)
C ARE THE Y GRID VALUES, Y(I) TILDA, AND MUST
C SATISFY DUAL CONDITIONS.
C NWK - INPUT. ON ENTRY WITH MODE = 1 THIS MUST BE SET TO THE
C DIMENSION OF THE ARRAY WK IN THE CALLING
C PROGRAM.
C NWKU - OUTPUT. THE ACTUAL NUMBER OF LOCATIONS NEEDED IN THE
C ARRAY WK.
C FO - OUTPUT ARRAY. VALUES OF THE INTERPOLATION FUNCTION AT
C THE GRID OF POINTS INDICATED BY NXO, XO, NYO, YO
C FO IS ASSUMED TO BE DIMENSIONED (NXO,NYO) IN THE
C CALLING PROGRAM.
C KER - OUTPUT. RETURN INDICATOR.
C = 0, NORMAL RETURN.
C = NONZERO, ERROR CONDITION ENCOUNTERED.
C
C ERROR MESSAGES
C NO. 1 FATAL SINGULAR MATRIX IN THE CALCULATION OF
C LOCAL THIN PLATE SPLINES. TRY LARGER
C VALUE FOR NPPR AND/OR MINPTS. (MINPTS
C IS IN SUBROUTINE LOLIP.)
C NO. 2 RECOVERABLE FIRST CALL TO LOTPS MUST BE WITH MODE=1
C NO. 3 FATAL PREVIOUS ERROR RETURN FROM SUBROUTINE
C LOCAL NOT CORRECTED.
C NO. 4 FATAL ARRAY IWK AND/OR WK NOT DIMENSIONED LARGE
C ENOUGH. REDIMENSION AS GIVEN BY NIWKU
C AND NWKU.
C NO. 5 RECOVERABLE MODE IS OUT OF RANGE.
C
C SUBROUTINES USED
C
C THIS PACKAGE: LOGRD, LOLIP, LOCAL, LOEVL.
C LINPACK: SGECO, SGESL
C SLATEC: SSORT, XERROR
C
C***END PROLOGUE