C 
C Subroutine file 'srfc.for' for specification and interpolation of
C smooth surfaces in the model in rectangular grids.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     SRFC1...Subroutine reading the input data for smooth surfaces.
C             SRFC1
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             SRFC2
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.......................................................................
C
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,POWERW,/
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     POWERW..Given grid values (6) correspond to the POWERW-th power of
C             interpolated function W.  The given grid values (6) are
C             thus raised to the (1/POWERW)-th power immediately after
C             reading and then interpolated.  POWERW=1 is recommended.
C     /...    Obligatory slash at the end of line for future extensions.
C     Default: IVAR1=0, IVAR2=0, IVAR3=0, SIGMA=0, POWERW=1.
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=======================================================================
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 successively, 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
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: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage location:
      REAL POWER(1)
C
      CALL VAL2(1,IABS(ISRFC),1,COOR,F,POWER)
      RETURN
      END
C
C=======================================================================
C