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