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