C SUBROUTINE FILE 'SRFC.FOR' FOR SPECIFICATION AND INTERPOLATION OF
C SMOOTH SURFACES IN THE MODEL IN RECTANGULAR GRIDS.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF THE FOLLOWING SUBROUTINES:
C     SRFC1...SUBROUTINE READING THE INPUT DATA FOR SMOOTH SURFACES.
C     SRFC2...SUBROUTINE EVALUATING THE FUNCTION VALUES AND THEIR FIRST
C             AND SECOND DERIVATIVES.  OUTSIDE THE SPECIFIED RECTANGULAR
C             GRID, THE FUNCTIONS ARE CONTINUED SMOOTHLY.
C SUBROUTINES SRFC1 AND SRFC2 SUPPORTING THE COMPLETE RAY TRACING
C ALGORITHM ONLY MEDIATE THE WORK OF SUBROUTINES VAL1 AND VAL2 WHICH
C MUST BE APPENDED.  IN ADDITION, SUBROUTINES CURVN1 (OR ITS ALTERNATIVE
C CURVB1), CURV2D (OR ITS ALTERNATIVE CURVBD), SURFB1, SURFBD, VAL3B1,
C VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ, INTRVL FROM THE
C SUBROUTINE PACKAGE 'FITPACK' BY ALAN KAYLOR CLINE, DEPARTMENT OF
C COMPUTER SCIENCES, UNIVERSITY OF TEXAS AT AUSTIN, ARE USED.  IN THE
C COMPLETE RAY TRACING, THIS SOFTWARE FILE 'SRFC.FOR' MAY BE  REPLACED
C BY ANY USER-DEFINED PACKAGE CONTAINING SUBROUTINES SRFC1 AND SRFC2
C WITH THE SAME NUMBER, TYPE AND MEANING OF THEIR PARAMETERS AS IN THIS
C FILE.
C
C IF MODEL VARIATIONS ARE TAKEN INTO ACCOUNT:
C     MODEL VARIATIONS ARE ASSUMED TO BE STORED WHILE EVALUATING THE
C     FUNCTION DESCRIBING A GIVEN SURFACE DURING THE INVOCATION OF
C     SUBROUTINE VAL OF FILE 'VAL.FOR' AND SUBSEQUENT ROUTINES OF FILE
C     'FIT.FOR'.  THE VARIATIONS ARE ASSUMED TO BE STORED IN REGISTER 1
C     OF THE SYSTEM VAR*.
C
C INPUT DATA (READ IN BY SUBROUTINE SRFC1):
C     THESE INPUT DATA DEFINE THE SURFACES.  THEY ARE READ IN BY
C     SUBROUTINE SRFC1.  THE NUMBER NSRFC OF THE SURFACES TO BE DEFINED
C     IS AN INPUT ARGUMENT OF SUBROUTINE SRFC1.  THE DATA ARE READ IN BY
C     THE LIST DIRECTED INPUT (FREE FORMAT).
C (1) NSRFC-TIMES (I.E. ONCE FOR EACH SURFACE) INPUT DATA (1A)+(1B):
C (1A) TEXTG,ISRFC
C     IDENTIFICATION OF THE SURFACE.
C     TEXTG...ANY STRING. ITS FIRST 3 CHARACTERS MUST DIFFER FROM 'END'.
C     ISRFC...INDEX OF THE SURFACE.
C (1B) 'INPUT DATA FOR ONE SURFACE', SEE BELOW.
C (2) TEXTE,AUX
C     END OF DATA.
C     TEXTE...STRING, THE FIRST 3 CHARACTERS OF WHICH MUST BE UPPER-CASE
C             'END'.
C     AUX...  ANY NUMBER OR A SLASH.
C FOR AN EXAMPLE REFER TO THE SAMPLE INPUT DATA FOR THE MODEL.
C
C INPUT DATA FOR ONE SURFACE:
C     THE DATA ARE READ IN BY THE LIST DIRECTED INPUT (FREE FORMAT).  IN
C     THE LIST OF INPUT DATA BELOW, EACH NUMBERED PARAGRAPH INDICATES
C     THE BEGINNING OF A NEW INPUT OPERATION (NEW READ STATEMENT).  IF
C     THE FIRST LETTER OF THE SYMBOLIC NAME OF THE INPUT VARIABLE IS
C     I-N, THE CORRESPONDING VALUE IN INPUT DATA MUST BE OF THE TYPE
C     INTEGER.  OTHERWISE, THE INPUT PARAMETER IS OF THE TYPE REAL.
C (1) IVAR1,IVAR2,IVAR3,SIGMA
C     THE FORM OF THE FUNCTION.
C     IVAR1,IVAR2,IVAR3... DENOTE THE FORM OF THE FUNCTION. THE FUNCTION
C             MUST BE OF THE FORM
C               F(X1,X2,X3) = W(A1,A2,A3)-B1-B2-B3   .
C             X1, X2, X3 ARE THE GENERAL COORDINATES. EACH OF A1, A2,
C             A3, B1, B2, B3 MUST BE EITHER:  (A) ONE OF GENERAL
C             COORDINATES X1, X2, X3,  OR (B) MUST BE LEFT OUT.  AT MOST
C             3 OF PARAMETERS A1-B3 MAY BE OF KIND (A).  NOTE THAT IVAR1
C             CONTROLS THE TYPE OF A1 AND B1, IVAR2 CONTROLS THE TYPE OF
C             A2 AND B2, IVAR3 CONTROLS THE TYPE OF A3 AND B3.
C             FOR IVAR1.EQ.0: A1, B1 ARE EMPTY (LEFT OUT),
C             FOR IVAR1.EQ.1: A1=X1, B1 IS EMPTY,
C             FOR IVAR1.EQ.2: A1=X2, B1 IS EMPTY,
C             FOR IVAR1.EQ.3: A1=X3, B1 IS EMPTY,
C             FOR IVAR1.EQ.-1: B1=X1, A1 IS EMPTY,
C             FOR IVAR1.EQ.-2: B1=X2, A1 IS EMPTY,
C             FOR IVAR1.EQ.-3: B1=X3, A1 IS EMPTY,
C             THE MEANING OF THE PARAMETERS IVAR2, IVAR3 IS SIMILAR.
C             EXAMPLES:
C             IVAR1: IVAR2: IVAR3: THE FORM OF THE FUNCTION:
C               1      2      3     F(X1,X2,X3)=W(X1,X2,X3)
C               3      1      2     F(X1,X2,X3)=W(X3,X1,X2)
C               1      2      0     F(X1,X2,X3)=W(X1,X2)
C               1      2     -3     F(X1,X2,X3)=W(X1,X2)-X3
C               1     -3      2     F(X1,X2,X3)=W(X1,X2)-X3
C             FUNCTION W IS INTERPOLATED BY MEANS OF SPLINES UNDER
C             TENSION.
C     SIGMA...IS THE TENSION FACTOR (ITS SIGN IS IGNORED).  THIS VALUE
C             INDICATES THE CURVINESS DESIRED.  IF ABS(SIGMA) IS NEARLY
C             ZERO (E.G. 0.001), THE RESULTING SURFACE IS APPROXIMATELY
C             THE TENSOR PRODUCT OF CUBIC SPLINES.  IF ABS(SIGMA) IS
C             LARGE (E.G. 50.), THE RESULTING SURFACE IS APPROXIMATELY
C             TRI-LINEAR.  IF SIGMA EQUALS ZERO, TENSOR PRODUCTS OF
C             CUBIC SPLINES RESULT.  A RECOMMENDED VALUE FOR SIGMA IS
C             APPROXIMATELY 1. IN ABSOLUTE VALUE.
C (2) NX(1),...,NX(NVAR)
C     THE NUMBERS OF GRID COORDINATES FOR THE INTERPOLATION.
C     THIS INPUT IS PERFORMED IF AT LEAST ONE OF IVAR1, IVAR2, IVAR3 IS
C     POSITIVE.
C     EACH OF NX(1),...,NX(NVAR) CORRESPONDS TO ONE POSITIVE VALUE OF
C     IVAR1, IVAR2, IVAR3 AND SPECIFIES THE NUMBER OF GRID COORDINATES
C     CORRESPONDING TO THAT INDEPENDENT VARIABLE OF FUNCTION W, SEE (1).
C     THE SIGN OF NX(1),...,NX(NVAR) IS IGNORED. NVAR (.LE.3) IS THE
C     NUMBER OF POSITIVE VALUES OF THE ABOVE QUANTITIES IVAR1, IVAR2,
C     IVAR3, I.E. THE NUMBER OF INDEPENDENT VARIABLES OF FUNCTION W,
C     SEE (1).
C (3) X1(1),...,X1(NX(1))
C     THE GRID COORDINATES CORRESPONDING TO THE FIRST INDEPENDENT
C     VARIABLE OF FUNCTION W, SEE (1).
C     THIS INPUT IS PERFORMED IF NX(1) IS SPECIFIED, SEE (2), AND IS NOT
C     ZERO. THE GRID COORDINATES MAY BE SPECIFIED IN ANY ORDER.
C (4) X2(1),...,X2(NX(2))
C     THE GRID COORDINATES CORRESPONDING TO THE SECOND INDEPENDENT
C     VARIABLE OF FUNCTION W, SEE (1).
C     THIS INPUT IS PERFORMED IF NX(2) IS SPECIFIED, SEE (2), AND IS NOT
C     ZERO. THE GRID COORDINATES MAY BE SPECIFIED IN ANY ORDER.
C (5) X3(1),...,X3(NX(3))
C     THE GRID COORDINATES CORRESPONDING TO THE THIRD INDEPENDENT
C     VARIABLE OF FUNCTION W, SEE (1).
C     THIS INPUT IS PERFORMED IF NX(3) IS SPECIFIED, SEE (2), AND IS NOT
C     ZERO. THE GRID COORDINATES MAY BE SPECIFIED IN ANY ORDER.
C (6) (((W(I,J,K),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1)),K=1,MAX(NX(3),1))
C     THE VALUES OF FUNCTION W AT GRID POINTS. FUNCTION VALUE W(I,J,K)
C     CORRESPONDS TO POINT (X1(I),X2(J),X3(K)).
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
      SUBROUTINE SRFC1(LUN,NSRFC)
      INTEGER LUN,NSRFC
C
C THIS SUBROUTINE READS THE INPUT DATA FOR THE SMOOTH SURFACES,
C DETERMINES THE PARAMETERS NECESSARY TO COMPUTE AN INTERPOLATORY
C FUNCTION ON A THREE DIMENSIONAL RECTANGULAR GRID, AND STORES THEM IN
C THE MEMORY.  THE FUNCTION DETERMINED CAN BE REPRESENTED AS A TENSOR
C PRODUCT OF SPLINES UNDER TENSION.  FOR ACTUAL MAPPING OF POINTS IT IS
C NECESSARY TO CALL THE SUBROUTINE SRFC2, WHICH ALSO RETURNS THE FIRST
C AND SECOND PARTIAL DERIVATIVES.  SUBROUTINE SRFC1 MAY BE CALLED
C SEVERAL TIMES.  THE SURFACES ARE INDEXED SUCCESIVELY, FOLLOWING THE
C SURFACES DEFINED DURING THE PREVIOUS INVOCATIONS.
C
C INPUT:
C     LUN...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE INPUT DATA.
C     NSRFC...NUMBER OF THE SURFACES FOR WHICH THE INPUT DATA ARE
C             SPECIFIED DURING THE CURRENT INVOCATION OF SRFC1.
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C NO OUTPUT.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL VAL1
C     VAL1, SORTV, READV... FILE 'VAL.FOR'.
C     CURVN1 OR CURVB1 (ALTERNATIVES), SURFB1, VAL3B1, SNHCSH, VGEN,
C              TERMS, TRIDEC, TRISOL... SUBROUTINE PACKAGE 'FITPACK'
C              (FILE 'FIT.FOR').
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      CHARACTER*3 TFUNCT(1)
      DATA TFUNCT/'   '/
C
      CALL VAL1(LUN,1,NSRFC,1,TFUNCT)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE SRFC2(ISRFC,COOR,F)
      INTEGER ISRFC
      REAL COOR(3),F(10)
C
C THIS SUBROUTINE EVALUATES THE FUNCTIONS DESCRIBING VARIOUS SMOOTH
C SURFACES IN THE MODEL AT A GIVEN POINT. THE THREE FIRST AND SIX SECOND
C PARTIAL DERIVATIVES ARE ALSO EVALUATED.  THE SPECIFIED FUNCTIONS ARE
C REPRESENTED AS A TENSOR PRODUCT OF SPLINES UNDER TENSION.  THE
C COEFFICIENTS OF THESE FUNCTIONS ARE PREPARED IN SUBROUTINE SRFC1, IN
C WHICH THE INPUT DATA CONCERNING THE FUNCTION OF EACH SURFACE ARE READ
C IN.
C
C INPUT:
C     ISRFC...INDEX OF A SURFACE.
C     COOR... ARRAY CONTAINING COORDINATES X1, X2, X3 OF THE GIVEN
C             POINT.
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     F...    THE VALUE AND THE FIRST AND SECOND PARTIAL DERIVATIVES F,
C             F1, F2, F3, F11, F12, F22, F13, F23, F33 OF THE FUNCTION
C             F(X1,X2,X3) DETERMINING THE SURFACE ISRFC AT THE GIVEN
C             POINT.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL VAL2
C     VAL2... FILE 'VAL.FOR'.
C     CURV2D OR CURVBD (ALTERNATIVES), SURFBD, VAL3BD, SNHCSH, DSPLNZ,
C              INTRVL... SUBROUTINE PACKAGE 'FITPACK' (FILE 'FIT.FOR').
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATION:
      REAL POWER(1)
C
      CALL VAL2(1,ISRFC,1,COOR,F,POWER)
      RETURN
      END
C
C=======================================================================
C
