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