C Program INV1 to evaluate the derivatives of the travel time with C respect to the model coefficients. C C Just a preliminary demo version, illustrating the usage of the C routines designed to calculate the variations of travel times with C respect to model parameters (FORTRAN77 source code files 'var.for', C 'spsp.for', 'ap.for'). C C Program INV1 assumes all model parameters (coefficients) stored in the C common block /VALC/ as in the submitted versions of user-defined model C specification FORTRAN77 source code files 'srfc.for', 'parm.for', and C 'val.for'. Thus, unlike the other parts of the complete ray tracing, C the INV1 program cannot work with user's modifications of the C subroutines SRFC1, SRFC2, PARM1, and PARM2. C C Main input data file read from the interactive device (*): C (1) 'INV1' C 'INV1' name of the input file described below. C Default: 'INV1'='inv1.out'. C C Input data 'INV1': C This main data file contains the names of other input and output C files. The data are read in by the list directed input (free C format). In the list of input data below, each numbered paragraph C indicates the beginning of a new input operation (new READ C statement). All input variables are of the type CHARACTER. Only C the first 80 characters of the strings are significant. C (1) 'POINTS','FTT','MODEL','DATA','SOFT','INV1LOG',/ C 'POINTS'... Input data file containing the coordinates of shot and C receiver points. Its structure is described below. C Ignored (not opened) if the 'data' filename below is C blank. C 'FTT'... Input data file containing the travel times from the C field measurements. Its structure is described below. C Ignored (not opened) if the 'DATA' filename below is C blank. C 'MODEL'... Input data file containing the model parameters. C See the file 'model.for' of the package 'MODEL'. C 'DATA'..Output file containing the computed travel times and their C derivatives with respect to the model coefficients. C Its structure is described below. C If the filename is blank, no file with objective prior C information is generated. C 'SOFT'..Output file containing the subjective prior information. C Its structure is described below. C If the filename is blank, no file with subjective prior C information is generated. C 'INV1LOG'... Output log file. Just for your information. C Not opened and generated if the 'DATA' filename above is C blank. C /... An obligatory slash at the end of line. C Default: 'MODEL'='model.dat', 'DATA'='data.out', C 'SOFT'='soft.out', 'INV1LOG'='inv1log.out'. C (2) DIST,VPOWER,ERRINT,ERRVEL,/ C DIST... Maximum distance between the source or receiver point and C the initial or end point of a synthetic ray. C VPOWER... Velocity power for the computation of the travel-time C check sum. If the VPOWER-th velocity power is expressed, C in all blocks of the model, in terms of explicit functions C of model coordinates, linearly homogeneous in their C coefficients (e.g. B-splines), the travel time minus its C check sum (see the log output file) should be zero within C rounding errors. Otherwise, the check sum may have no C sense. C VPOWER=0: no check sum is evaluated and printed. C ERRINT..Weighting factor of the posterior interface B-spline C coefficient error. C ERRVEL..Weighting factor of the posterior velocity B-spline C coefficient error. C /... An obligatory slash at the end of line. C Default: VPOWER=0.0, ERRINT=1.0, ERRVEL=1.0. C (3) Any times the following data: C 'RAYS','ENDPOINTS','INITIALPOINTS' C 'RAYS'... File with the quantities stored along rays (see C C.R.T.5.5.1). C 'ENDPOINTS'... File with the quantities at the points of C intersection of rays with the specified surface at which C the receivers are situated for the case of two-point ray C tracing (see C.R.T.5.5.2). C If this filename is not blank, just the two-point rays C with minimum travel time at each receiver are considered. C If this filename is blank, all rays are taken into C account. C Attention: All rays taken into account must start in some C of the specified sources and terminate in some of the C specified receivers, see the input file 'FTT'. C 'INITIALPOINTS'... File with the quantities at the initial points C of rays, corresponding to the above file rays or endpoints C (see C.R.T.6.1). C (4) / (a slash). C (5) (NW1(I),NW2(I),NW3(I),I=1,NW),/ C List of partial derivatives included in the Sobolev scalar product C which is assumed to represent subjective prior information about C the model. C NW1,NW2,NW3... Orders of partial derivatives with respect to C X1,X2,X3 coordinates. C (6) ((WCS(I,J),I=1,J),J=1,NW) C Elements of the constant symmetric weighting matrix of the Sobolev C scalar product. C WCS(I,J)... Coefficient of the product of C (NW1(I),NW2(I),NW3(I))-th and (NW1(J),NW2(J),NW3(J))-th C partial derivatives of functions in the Sobolev scalar C product. C C Input data file 'POINTS': C This data file contains the coordinates of shot and receiver C points. The data are read in by the list directed input C (free format). In the list of input data below, each numbered C paragraph indicates the beginning of a new input operation (new C READ statement). The CHARACTER strings are explicitly mentioned C in this description. Otherwise, if the first letter of the C symbolic name of the input variable is I-N, the corresponding C value in input data must be of the type INTEGER. Otherwise, the C input parameter is of the type REAL. C (1) Several strings terminated by / (a slash). C (2) List of the sources and receivers: Any times the following data: C (2.1) POINT,X1,X2,X3 C POINT...CHARACTER*11 string containing the name of the source or C receiver point. C X1,X2,X3... Coordinates of the source or receiver point. C (3) / (a slash) or the end of file. C C Input data file 'FTT': C This data file contains the travel times from the field C measurements. The data are read in by the list directed input C (free format). In the list of input data below, each numbered C paragraph indicates the beginning of a new input operation (new C READ statement). The CHARACTER strings are explicitly mentioned C in this description. Otherwise, if the first letter of the C symbolic name of the input variable is I-N, the corresponding C value in input data must be of the type INTEGER. Otherwise, the C input parameter is of the type REAL. C (1) Several strings terminated by / (a slash). C (2) List of the travel times: Any times the following data (2.1): C (2.1) SOURCE,RECEIVER,TFIELD,TERR C SOURCE..CHARACTER*11 string containing the name of the source C point. C RECEIVER... CHARACTER*11 string containing the name of the C receiver point. The source and receiver points may be C mutually interchanged. C TFIELD..Travel time from a field measurement. C TERR... Error of the travel time from a field measurement. C (3) / (a slash) or the end of file. C C Output file 'DATA': C (1) ND,NM,(INDM(I),I=1,NM) C ND... The number of data (i.e. the number of equations). C NM... Number of unknown model coefficients. C INDM... Indices of the model coefficients influencing the travel C times. The indices correspond to the relative location in C the memory. C B-spline coefficients are listed in the same order as the C grid velocities in the file 'MODEL'. C (2) ND-times the following data (2.1): C (2.1) KD,RD,ED,WD,(GD(I),I=1,NM) C KD... Index of the field travel time within the file 'ftt'. C The field travel times are indexed consecutively 1,2,3,... C RD... Field travel time minus the computed synthetic travel C time. In the case of multiple two-point rays, the first C arrival of them is considered. C ED... Prior error of the above travel time difference. C Identical to TERR, file 'FTT', part (3). I.e. the C synthetic travel time is assumed to be sufficiently C accurate. C WD... Field travel time (stored for the purpose to assess a C posterior data misfit covariance weighting matrix). C GD... Derivatives of the synthetic travel time with respect to C the model coefficients INDM. C C Output file 'SOFT': C (1) (INDM(I),I=1,NM) C INDM... Indices of the model coefficients influencing the travel C times. C (2) (WM(I),I=1,NM) C RS... Parameters (coefficients) of the initial model. C (3) (WM(I),I=1,NM) C WM... Relative posterior errors of the model parameters. C WM=ERRVEL for material parameters (see (2) of the input C data 'INV1'), C WM=ERRINT for structural interfaces (see (2) of the input C data 'INV1'). C Default: ERRVEL=1.0, ERRINT=1.0. C (4) ((CS(I,J),I=1,J),J=1,NM): C CS... Relative inverse subjective prior information covariance C matrix. Note that subjective data are assumed to be minus C the above initial model parameters, and the matrix GS C projecting model parameters onto the subjective data is C assumed to be identity. C C Output log file 'LOG': C (1) For each considered ray: C (1.1) SOURCE,RECEIVER,TFIELD,TDIF,SDIST,RDIST,CHECK C SOURCE..Name of the source point. C RECEIVER... Name of the receiver point. C TFIELD..Travel time from a field measurement. C TDIF... Field travel time minus the minimum synthetic travel time. C SDIST...Distance between the source and the initial point of the C synthetic ray. C RDIST...Distance between the receiver and the end point of the C synthetic ray. C CHECK...Synthetic travel time minus the travel time resulting from C the derivatives of the theoretical travel time with C respect to the model coefficients. This quantity should C not exceed in order the numerical error of the synthetic C travel time. C In this version defined just for the models described in C the terms of velocity. C C Date: 1994, January 23 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Common block /VALC/: INCLUDE 'val.inc' C None of the storage locations of the common block are altered. C C Common block /POINTC/: INCLUDE 'pointc.inc' C C----------------------------------------------------------------------- C C Filenames: CHARACTER*80 FILE0,FILE1,FILE2,FILE3,FILE4,FILE5,FILE6 C C Logical unit numbers: INTEGER LU0,LU1,LU2,LU3,LU4,LU5,LU7,IU1,IU2 PARAMETER (LU0=10) PARAMETER (LU1=11) PARAMETER (LU2=12) PARAMETER (LU3=13) PARAMETER (LU4=14) PARAMETER (LU5=15) PARAMETER (LU7=17) C C Input data: INTEGER MP,MT PARAMETER (MP=1000) PARAMETER (MT=4000) CHARACTER*1 TEXT(10) CHARACTER*11 POINT(MP),SRC,REC INTEGER NP,NT,KS(MT),KR(MT) REAL DIST,DIST2,VPOWER,ERRINT,ERRVEL REAL COOR(3,MP),TFIELD(MT),TERR(MT) INTEGER MW PARAMETER (MW=12) INTEGER NW,NW1(MW+1),NW2(MW),NW3(MW) REAL WCS(MW*(MW+1)/2) C POINT...Names of the source and receiver points. C NP... Number of source and receiver points. C NT... Number of field travel times. C KS(I)...Index of the source point corresponding to the I-th field C travel time. C KR(I)...Index of the receiver point corresponding to the I-th C field travel time. C DIST... Maximum distance between the source or receiver point and C the initial or end point of a synthetic ray. C DIST2...DIST**2 C VPOWER... Velocity power for the computation of the travel-time C check SUM. C COOR... Coordinates of the source or receiver points. C TFIELD..Field travel times. C TERR... Field travel time errors. C NW... Number of specified partial derivatives. C NW1,NW2,NW3... Orders of partial derivatives. C WCS(I,J)... Coefficient matrix of the Sobolev scalar product. C C Output data: variations of the synthetic travel time: INTEGER NSUM,NM,INDM(NPAR) REAL SUM(NPAR) C NM... Number of the unknown model parameters. C INDM... Indices of the unknown model parameters. C C Output data: subjective prior information: INTEGER MM,MCS C MM... Maximum number of model parameters. PARAMETER (MM=200) PARAMETER (MCS=MM*(MM+1)/2) REAL CS(MCS) C C Auxiliary storage locations: INTEGER IS,IR,IT,ND,IRAYTT,I,J,K REAL TTMIN(MT) REAL TT,TTAUX,TDIF,XI1,XI2,XI3,XE1,XE2,XE3,PI(6),PE(6) REAL SI,SE,RI,RE,SDIST,RDIST,WEIGHT C IS... Index of the source point. C IR... Index of the receiver point. C IT... Index of the field travel time. C ND... Number of synthetic travel times corresponding to the C field travel times. C IRAYTT..Index of the last processed ray. C I,J,K...Temporary storage locations. C TTMIN...Minimum synthetic travel times corresponding to the C individual field travel times. C TT... Synthetic travel time. C TTAUX...Temporary storage location. C TDIF... Field travel time minus the minimum synthetic travel time. C XI1,XI2,XI3,XE1,XE2,XE3... Coordinates of the initial and end C points of the last processed ray. C PI,PE...Slowness vectors at the initial and end points of the C last processed ray. C SI,SE,RI,RE... Squares of the distances between the source or C receiver points and the initial or end points of the ray. C SDIST...Distance between the source and the initial point of the C synthetic ray. C RDIST...Distance between the receiver and the end point of the C synthetic ray. C WEIGHT..Weighting factor corresponding to the derivatives C (average quadratic value of the derivatives on input) C C More than 23*MP+16*MT+4*MM*(MM+1) Bytes of memory required. C C....................................................................... C C Opening data files and reading the input data: C C Main input data file read from the interactive device (*): WRITE(*,'(A)') ' Enter the name of the main input data file: ' FILE0='inv1.dat' READ(*,*) FILE0 WRITE(*,'(A)') '+ ' C C Input data 'INV1': OPEN(LU0,FILE=FILE0,STATUS='OLD') FILE2='model.dat' FILE4='data.out' FILE5='soft.out' FILE6='inv1log.out' READ(LU0,*) FILE0,FILE1,FILE2,FILE4,FILE5,FILE6 VPOWER=0. ERRINT=1. ERRVEL=1. READ(LU0,*) DIST,VPOWER,ERRINT,ERRVEL DIST2=DIST*DIST OPEN(LU4,FILE=FILE2,STATUS='OLD') CALL MODEL1(LU4) CLOSE(LU4) C C Number of unknown model parameters: CALL SOFT(2,0,0,0,0,0,0,0.,NM,INDM,CS) WRITE(*,'(A,I4,A)') '+',NM,' model parameters' C IF(FILE4.EQ.' ') THEN C No objective information file generated. 10 CONTINUE FILE1=' ' READ(LU0,*,END=70) FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') THEN C GO TO (subjective prior information) GO TO 80 END IF GO TO 10 END IF C C....................................................................... C C Reading source and receiver points: C OPEN(LU4,FILE=FILE0,STATUS='OLD') NP=-1 READ(LU4,*,END=2) TEXT 1 CONTINUE NP=NP+1 IF(NP.GE.MP) THEN PAUSE 'Error: Too many source and receiver points' STOP END IF POINT(NP+1)=' ' READ(LU4,*,END=2) POINT(NP+1),(COOR(I,NP+1),I=1,3) IF(POINT(NP+1).NE.' ') GO TO 1 2 CONTINUE CLOSE(LU4) C C Reading field travel times: C OPEN(LU4,FILE=FILE1,STATUS='OLD') NT=0 READ(LU4,*,END=2) TEXT 3 CONTINUE NT=NT+1 IF(NT.GT.MT) THEN PAUSE 'Error: Too many field travel times' STOP END IF SRC=' ' READ(LU4,*,END=9) SRC,REC,TFIELD(NT),TERR(NT) IF(SRC.EQ.' ') THEN GO TO 9 END IF DO 4 I=1,NP IF(SRC.EQ.POINT(I)) THEN KS(NT)=I GO TO 5 END IF 4 CONTINUE PAUSE 'Error: Source name not recognized' STOP 5 CONTINUE DO 6 I=1,NP IF(REC.EQ.POINT(I)) THEN KR(NT)=I GO TO 7 END IF 6 CONTINUE PAUSE 'Error: Receiver name not recognized' STOP 7 CONTINUE GO TO 3 9 CONTINUE NT=NT-1 CLOSE(LU4) C C....................................................................... C C Computing quantities describing objective prior information: C OPEN(LU5,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(LU7,FILE=FILE6) WRITE(*,*) C KS(NT+1)=NP+1 KR(NT+1)=NP+1 TFIELD(NT+1)=0. IRAY=0 IWAVE=0 NSUM=IPAR(IPAR(IPAR(2))) DO 12 I=1,NT TTMIN(I)=999999. 12 CONTINUE C C Loop for the files with computed rays 20 CONTINUE FILE1=' ' READ(LU0,*,END=70) FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') THEN GO TO 70 END IF I=INDEX(FILE1,' ') J=INDEX(FILE2,' ') K=INDEX(FILE3,' ') WRITE(*,'(''+Processing: '',3A)') FILE1(1:I),FILE2(1:J), * FILE3(1:K) WRITE(*,*) OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD') IF(FILE2.EQ.' ') THEN IU1=0 IU2=LU1 ELSE IU1=LU1 IU2=LU2 OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD') END IF OPEN(LU3,FILE=FILE3,FORM='UNFORMATTED',STATUS='OLD') IRAYTT=0 C C Loop for the endpoints of rays 30 CONTINUE C Reading the results of the complete ray tracing CALL AP00(IU1,IU2,LU3) IF(IPT.LE.1.AND.IRAYTT.NE.0)THEN C New ray - recording results for the last ray iraytt C loop for field travel times - searching for two-point ray DO 39 IT=1,NT IS=KS(IT) IR=KR(IT) SI=(COOR(1,IS)-XI1)**2+(COOR(2,IS)-XI2)**2 * +(COOR(3,IS)-XI3)**2 RI=(COOR(1,IR)-XI1)**2+(COOR(2,IR)-XI2)**2 * +(COOR(3,IR)-XI3)**2 SE=(COOR(1,IS)-XE1)**2+(COOR(2,IS)-XE2)**2 * +(COOR(3,IS)-XE3)**2 RE=(COOR(1,IR)-XE1)**2+(COOR(2,IR)-XE2)**2 * +(COOR(3,IR)-XE3)**2 IF(SE.LE.SI.AND.RI.LE.RE) THEN C interchanging source and receiver points SI=SE RE=RI IS=KR(IT) IR=KS(IT) END IF IF((SI.LE.DIST2.AND.RE.LE.DIST2)) THEN C Synthetic ray may correspond to the field travel time C check for ray directions near the source C (allowable angle deviation +-30 deg: cosine**2=0.750) C (allowable angle deviation +-15 deg: cosine**2=0.933) TTAUX=(COOR(1,IR)-COOR(1,IS))*(XE1-XI1) * +(COOR(2,IR)-COOR(2,IS))*(XE2-XI2) * +(COOR(3,IR)-COOR(3,IS))*(XE3-XI3) IF(TTAUX.GT.0..AND.TTAUX*TTAUX.GT. * 0.933*((COOR(1,IR)-COOR(1,IS))**2 * +(COOR(2,IR)-COOR(2,IS))**2 * +(COOR(3,IR)-COOR(3,IS))**2) * *((XE1-XI1)**2+(XE2-XI2)**2+(XE3-XI3)**2) ) THEN C Synthetic ray corresponds to the field travel time TTAUX=TT-PI(1)*(COOR(1,IS)-XI1) * -PI(2)*(COOR(2,IS)-XI2) * -PI(3)*(COOR(3,IS)-XI3) * +PE(1)*(COOR(1,IR)-XE1) * +PE(2)*(COOR(2,IR)-XE2) * +PE(3)*(COOR(3,IR)-XE3) IF(TTAUX.LT.TTMIN(IT)) THEN TTMIN(IT)=TTAUX C Possible minimum synthetic travel time SDIST=SQRT(SI) RDIST=SQRT(RE) WRITE(LU5) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM) END IF END IF END IF 39 CONTINUE IF(NT.EQ.0) THEN WRITE(LU5) NT+1,TT,TT,0.,0.,(SUM(I),I=1,NSUM) END IF IRAYTT=0 END IF IF(IWAVE.EQ.0) THEN GO TO 60 END IF C *** for future extensions (selection of two-point rays): C IF(IU1.NE.0) THEN C CALL AP30(IREC) C IF(IREC.EQ.0) THEN C IF(IPT.LE.1.) THEN C WRITE(*,'(''+WAVE:'',I3,'' RAY:'',I4,'' POINT:'', C * I4)') IWAVE,IRAY,IPT C END IF C GO TO 30 C END IF C END IF C *** IF(IPT.EQ.1.OR.MOD(IPT,10).EQ.0) THEN WRITE(*,'(''+WAVE:'',I3,'' RAY:'',I4,'' POINT:'',I4)') * IWAVE,IRAY,IPT END IF IRAYTT=IRAY XI1=YI(3) XI2=YI(4) XI3=YI(5) XE1=Y(3) XE2=Y(4) XE3=Y(5) CALL AP01(TT,TTAUX) CALL AP02(PI,PE) CALL AP29(NSUM,SUM) GO TO 30 C End of the loop for endpoints of rays C 60 CONTINUE CLOSE(LU1) CLOSE(LU2) GO TO 20 C C All minimum travel times and their derivatives are stored in the C scratch file LU5. C C....................................................................... C C Writing objective prior information: C 70 CONTINUE OPEN(LU4,FILE=FILE4) I=MAX0(INDEX(FILE4,' ')-1,11) WRITE(*,'(''+Generating the output: '',A)') FILE4(1:I) ND=0 DO 71 I=1,NT IF(TTMIN(I).LT.999999.) THEN ND=ND+1 END IF 71 CONTINUE C ND is the number of equations. C REWIND(LU5) WRITE(LU4,'(I9,5I13)') ND,NM WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM) WRITE(LU7,'(2A)') ' SOURCE RECEIVER TFIELD ', * 'TFIELD-TT SDIST RDIST TT-CHECKSUM' 73 CONTINUE READ(LU5,END=79) IT,TT,TTAUX,SDIST,RDIST,(SUM(I),I=1,NSUM) TTMIN(NT+1)=TTAUX IF(TTAUX.EQ.TTMIN(IT)) THEN C Minimum synthetic travel time C C System of linear equations: TDIF=TFIELD(IT)-TTAUX WRITE(LU4,'(I9,4X,6G13.6)') IT,TDIF,TERR(IT),TFIELD(IT) WRITE(LU4,'(6G13.6)') (SUM(INDM(I)),I=1,NM) C C Check sums and log output: IF(VPOWER.NE.0.) THEN TTAUX=0. DO 74 I=1,NM J=INDM(I) IF(IPAR(IPAR(IPAR(1))).LT.J) THEN IF(SUM(J).NE.0.) THEN TTAUX=TTAUX+RPAR(J)*SUM(J) END IF END IF 74 CONTINUE IS=KS(IT) IR=KR(IT) TTAUX=TT+VPOWER*TTAUX WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR), * TFIELD(IT),TDIF,SDIST,RDIST,TTAUX ELSE WRITE(LU7,'(2(1X,A),5F12.6)') POINT(IS),POINT(IR), * TFIELD(IT),TDIF,SDIST,RDIST END IF END IF GO TO 73 C 79 CONTINUE CLOSE(LU4) CLOSE(LU5) CLOSE(LU7) C C....................................................................... C C Subjective prior information: C 80 CONTINUE IF(FILE5.NE.' ') THEN C IF(NM*(NM+1)/2.GT.MCS) THEN PAUSE 'Error: Too many model parameters' STOP END IF C C Opening output file: OPEN(LU4,FILE=FILE5) I=MAX0(INDEX(FILE5,' ')-1,11) WRITE(*,'('' Generating the output: '',A)') FILE5(1:I) C C Reading prior subjective information coefficients: DO 81 I=1,MW+1 NW1(I)=-1 81 CONTINUE READ(LU0,*) (NW1(I),NW2(I),NW3(I),I=1,MW),NW1(MW+1) DO 82 I=1,MW IF(NW1(I).EQ.-1) THEN NW=I-1 GO TO 83 END IF 82 CONTINUE PAUSE 'Error: Too large energy coefficient matrix' STOP 83 CONTINUE READ(LU0,*) (WCS(I),I=1,NW*(NW+1)/2) C C Generating prior subjective information covariance matrix: DO 84 I=1,MCS CS(I)=0. 84 CONTINUE DO 89 J=1,NW DO 88 I=1,J WEIGHT=WCS(J*(J-1)/2+I) IF(WEIGHT.NE.0.) THEN CALL SOFT(2,NW1(I),NW2(I),NW3(I),NW1(J),NW2(J),NW3(J), * WEIGHT,NM,INDM,CS) END IF 88 CONTINUE 89 CONTINUE C C (1,2) Writing the initial model parameters: WRITE(LU4,'(I9,5I13)') (INDM(I),I=1,NM) WRITE(LU4,'(6G13.6)') (RPAR(INDM(I)),I=1,NM) C C (3) Writing relative posterior errors of the model parameters: DO 93 I=1,NM J=INDM(I) IF(IPAR(IPAR(IPAR(1))).LT.J) THEN SUM(I)=ERRVEL ELSE SUM(I)=ERRINT END IF 93 CONTINUE WRITE(LU4,'(24F3.0)') (SUM(I),I=1,NM) C C (4) Writing the prior subjective information covariance matrix: DO 94 J=1,NM WRITE(LU4,'(6G13.6)') (CS(I),I=J*(J-1)/2+1,J*(J+1)/2) 94 CONTINUE CLOSE(LU4) C END IF C STOP END C C======================================================================= C SUBROUTINE SOFT(NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,WEIGHT,NM,INDM,CS) C INTEGER NCLASS,NA1,NA2,NA3,NB1,NB2,NB3,NM,INDM(*) REAL WEIGHT,CS(*) C C This subroutine accumulates the prior subjective information C covariance matrix describing the smoothness of the functions C interpolated by means of subroutines of the file 'val.for'. C C Input: C NCLASS..Number of classes required. C NA1,NA2,NA3... Orders of X1,X2,X3-partial derivatives of the first C basis function in the product being integrated. C NB1,NB2,NB3... Orders of X1,X2,X3-partial derivatives of the C second basis function in the product being integrated. C WEIGHT..Weighting factor corresponding to the derivatives C (average quadratic value of the derivatives on input) C CS... Symmetric matrix to be regularized: C Contribution corresponding to the derivatives C is added to the matrix CS C None of the input parameters are altered. C C Output: C NM... Number of the unknown model parameters. C INDM... Indices of the unknown model parameters. C CS... Resulting symmetric subjective prior information C covariance matrix. C C Common block: INCLUDE 'val.inc' C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SPSP C C Error messages: C 359a.. Incorrect number of classes: C The number of classes required is zero, negative or C greater than the number of classes defined. C 359b.. Insufficient memory: C The dimension of the buffer B is not sufficient. C C Date: 1992, September 7 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C INTEGER MM,MB C MM... Maximum number of model parameters. PARAMETER (MM=200) PARAMETER (MB=MM*(MM+1)/2+3*(20*21/2)) REAL B(MB) INTEGER ICLASS,JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,NVAR INTEGER NX(3),NX1,NX2,NX3 EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3) INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7 EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3) EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6) EQUIVALENCE (JADR(7),JADR7) INTEGER I1,I2,I3,I,J,N REAL WDENS,SIGMA C C B... Auxiliary storage locations for partial covariance C matrices. C ICLASS..Index of the class. C JGROUP..Address of a current group. C LFUNCT,MFUNCT,JFUNCT... Addresses of the first, last and arbitrary C functions of the group. C LADR,MADR,IADR... Addresses of the first, last and arbitrary C parameters of the current function. C NVAR... Number of the independent variables A1, A2, A3 of the C interpolated function W. C NX=(NX1,NX2,NX3)... Numbers of grid lines. C JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... Addresses of C parameters describing the interpolated function (grid C coordinates, B-spline coefficients, B-spline basis C functions). C I1,I2,I3,I,J,N... Local auxiliary variables. C WDENS...Weighting factor divided by volume (area,length). C SIGMA...Tension factor. C C....................................................................... C NM=0 IF(NCLASS.LT.1.OR.IPAR(0).LT.NCLASS) THEN PAUSE 'Error 359a in SOFT: Incorrect number of classes' STOP END IF C Loop for classes: DO 92 ICLASS=1,NCLASS C Loop for groups of the current class: DO 91 JGROUP=IPAR(ICLASS-1)+1,IPAR(ICLASS) LFUNCT=IPAR(JGROUP-1)+1 MFUNCT=IPAR(JGROUP) MADR =IPAR(LFUNCT-1) C C Loop for functions of the current group: DO 90 JFUNCT=LFUNCT,MFUNCT C Starting and end addresses of the parameters describing the C function LADR=MADR+1 MADR=IPAR(JFUNCT) IF(LADR.LE.MADR) THEN WDENS=WEIGHT C Tension factor SIGMA=RPAR(LADR+5) C C The number, types and values of the independent variables C of the interpolated function: C Initial address IADR=LADR+6 C Initial number of the independent variables NVAR=0 NX1=1 NX2=1 NX3=1 JADR1=0 JADR2=0 JADR3=0 JADR4=0 C Loop for the possible independent variables: DO 20 I=LADR+2,LADR+4 C Type of the possible independent variable: J=IPAR(I) IF(J.GT.0) THEN N=IABS(IPAR(IADR)) IF(N.GE.2) THEN NVAR=NVAR+1 NX(NVAR)=N ELSE IF(N.EQ.1) THEN JADR(NVAR+1)=JADR(NVAR+1)+1 END IF IADR=IADR+1 END IF 20 CONTINUE JADR5=0 JADR6=0 JADR7=0 C C Interpolated function W: JADR1=IADR+JADR1 IF(NVAR.LE.0) THEN C No independent variable. JADR4=JADR1 ELSE JADR2=JADR1+NX1+JADR2 WDENS=WDENS/(RPAR(JADR2-1)-RPAR(JADR1)) IF(NVAR.EQ.1) THEN C One independent variable: JADR4=JADR2 JADR5=JADR2+NX1 ELSE JADR3=JADR2+NX2+JADR3 WDENS=WDENS/(RPAR(JADR3-1)-RPAR(JADR2)) IF(NVAR.EQ.2) THEN C Two independent variables: JADR4=JADR3 JADR5=JADR3+NX1*NX2 JADR6=JADR4+5*NX1 ELSE C Three independent variables: JADR4=JADR3+NX3+JADR4 JADR5=JADR4+NX1*NX2*NX3 JADR6=JADR5+5*NX1 JADR7=JADR6+5*NX2 WDENS=WDENS/(RPAR(JADR4-1)-RPAR(JADR3)) END IF END IF END IF C WDENS is weighting factor divided by volume (area,length). C JADR4=JADR4-1 N=NX1*NX2*NX3 IF(WEIGHT.EQ.0.) THEN DO 72 I2=1,N INDM(NM+I2)=JADR4+I2 72 CONTINUE NM=NM+N ELSE IF(NA1.EQ.NB1.AND.NA2.EQ.NB2.AND.NA3.EQ.NB3) THEN I=N*(N+1)/2 J=0 ELSE I=N*N J=1 END IF IF(NA1.EQ.NB1) THEN I1=NX1*(NX1+1)/2 ELSE I1=NX1*NX1 END IF IF(NA2.EQ.NB2) THEN I2=NX2*(NX2+1)/2 ELSE I2=NX2*NX2 END IF IF(NA3.EQ.NB3) THEN I3=NX3*(NX3+1)/2 ELSE I3=NX3*NX3 END IF I1=I1+I I2=I2+I1 I3=I3+I2 IF(I3.GT.MB) THEN PAUSE 'Error 359b in SOFT: Insufficient memory' STOP END IF CALL SPSP(NA1,NA2,NA3,NB1,NB2,NB3,NX1,NX2,NX3, * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3), * RPAR(JADR5),RPAR(JADR6),RPAR(JADR7), * SIGMA,WDENS,B(1),B(I+1),B(I1+1),B(I2+1)) I=NM*(NM+1)/2 DO 82 I2=1,N INDM(NM+I2)=JADR4+I2 I=I+NM DO 81 I1=1,I2 I=I+1 IF(J.EQ.0) THEN CS(I)=CS(I)+B(I2*(I2-1)/2+I1) ELSE CS(I)=CS(I)+B(N*(I2-1)+I1)+B(N*(I1-1)+I2) END IF 81 CONTINUE 82 CONTINUE NM=NM+N C Contribution corresponding to the derivatives C is added to the matrix CS. END IF C END IF 90 CONTINUE 91 CONTINUE 92 CONTINUE C End of loops for functions. C RETURN END C C======================================================================= C INCLUDE 'modelv.for' INCLUDE 'metric.for' INCLUDE 'srfc.for' INCLUDE 'parmv.for' INCLUDE 'valv.for' INCLUDE 'fitv.for' INCLUDE 'var.for' INCLUDE 'spsp.for' INCLUDE 'means.for' INCLUDE 'ap.for' INCLUDE 'apvar.for' C C======================================================================= C