C
C Program GRDTE to compute the values of a real or complex function, C described in terms of the Taylor expansions of its amplitude and C phase, on a given grid C C Version: 6.20 C Date: 2008, June 10 C C Coded by Karel Zacek C Department of Geophysics, Charles University Prague C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Names of the output files: C GRDTER='string' C GRDTEI='string'... Names of the output files with the calculated C grid values of the given complex-valued function. C File GRDTER will contain the real part and file GRDTEI C will contain the imaginary part of the function. C If the string is blank, the corresponding file is not C generated. C For general description of files with gridded data refer C to file forms.htm. C Defaults: GRDTER=' ', GRDTEI=' ' C Data specifying the dimensions of the grid for discretization: C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C O1=real, O2=real, O3=real ... Coordinates of the origin of the C grid. C Defaults: O1=0., O2=0., O3=0. C D1=real... Grid interval along the X1 axis. C Default: D1=1. C D2=real... Grid interval along the X2 axis. C Default: D2=1. C D3=real... Grid interval along the X3 axis. C Default: D3=1. C Data specifying the origins of individual discretized functions: C N10, N20, N30, O10, O20, O30, D10, D20 and D30... Specify the set C of origins X0j for individual shifted functions. Their C meaning is analogous to parameters N1, N2, N3, O1, O2, O3, C D1, D2 and D3. C Defaults: N10=1, N20=1, N30=1, O10=0., O20=0., O30=0., C D10=1., D20=1., D30=1. C Additional data specific to this program: C COEFFICIENT=real C where COEFFICIENT stands for any of the following 80 names C AR0, AI0, TR0, TI0, C AR1, AR2, AR3, AI1, AI2, AI3, C TR1, TR2, TR3, TI1, TI2, TI3, C AR11, AR12, AR22, AR13, AR23, AR33, C AI11, AI12, AI22, AI13, AI23, AI33, C TR11, TR12, TR22, TR13, TR23, TR33, C TI11, TI12, TI22, TI13, TI23, TI33, C AR111, AR112, AR122, AR222, AR113, C AR123, AR223, AR133, AR233, AR333, C AI111, AI112, AI122, AI222, AI113, C AI123, AI223, AI133, AI233, AI333, C TR111, TR112, TR122, TR222, TR113, C TR123, TR223, TR133, TR233, TR333, C TI111, TI112, TI122, TI222, TI113, C TI123, TI223, TI133, TI233, TI333. C The calculated complex-valued function has the form of C F(Xm) = A(Xm) exp(T(n)) , C where the complex-valued amplidtude is C A(Xm) = A0 + Aj*Yj + Ajk*Yj*Yk/2 + Ajkl*Yj*Yk*Yl/6 C and the complex-valued phase is C T(Xm) = T0 + Tj*Yj + Tjk*Yj*Yk/2 + Tjkl*Yj*Yk*Yl/6 , C with C Yj = Xj - X0j . C Here C A0 = AR0 + i*AI0 , C Aj = ARj + i*AIj , C Ajk = ARjk + i*AIjk , C Ajkl = ARjkl + i*AIjkl , C T0 = TR0 + i*TI0 , C Tj = TRj + i*TIj , C Tjk = TRjk + i*TIjk , C Tjkl = TRjkl + i*TIjkl . C The origins X0j of individual discretized functions are C specified by input parameters N10, N20, N30, O10, O20, C O30, D10, D20 and D30. The individual functions are C written as individual snapshots. It means that they may C be processed or displayed like N4=N10*N20*N30 spatial C grids. C Defaults: COEFFICIENT=0. C Optional parameters specifying the form of the real-valued quantities C written to the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C NUMLIN=positive integer ... Number of reals in one line of the C output file. See the description in file C forms.for. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C....................................................................... C INTEGER N1,N2,N3,N10,N20,N30,LU1,LU2,N1N2N3,PT,PTM REAL D1,D2,D3,O1,O2,O3,D10,D20,D30,O10,O20,O30 INTEGER IN1,IN2,IN3,IN10,IN20,IN30,I,J,K,L REAL X(3),X0(3) REAL FR,FI,AR,AI,TR,TI,AR0,AI0,TR0,TI0,DX,DXI,DXJ,DXK REAL AR1(3),AR2(3,3),AR3(3,3,3),TR1(3),TR2(3,3),TR3(3,3,3) REAL AI1(3),AI2(3,3),AI3(3,3,3),TI1(3),TI2(3,3),TI3(3,3,3) REAL E,C,S CHARACTER*80 FILE,OUTR,OUTI PARAMETER (LU1=1,LU2=2) C C----------------------------------------------------------------------- C C Reading main input data: WRITE(*,'(A)') '+GRDTE: Enter input filename: ' FILE=' ' READ(*,*) FILE IF(FILE.EQ.' ') THEN C GRDTE-01 CALL ERROR('GRDTE-01: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. END IF WRITE(*,'(A)') '+GRDTE: Reading, calculating... ' C C Reading grid dimensions: CALL RSEP1(LU1,FILE) CALL RSEP3T('GRDTER',OUTR,' ') CALL RSEP3T('GRDTEI',OUTI,' ') IF (OUTR.EQ.' ' .AND. OUTI.EQ.' ') THEN C GRDFFT-03 CALL ERROR('GRDFFT-03: No output files specified.') ENDIF CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3I('N10',N10,1) CALL RSEP3I('N20',N20,1) CALL RSEP3I('N30',N30,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3R('D10',D10,1.) CALL RSEP3R('D20',D20,1.) CALL RSEP3R('D30',D30,1.) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3R('O10',O10,0.) CALL RSEP3R('O20',O20,0.) CALL RSEP3R('O30',O30,0.) C N1N2N3=N1*N2*N3 IF(2*N1N2N3.GT.MRAM) THEN C GRDTE-02 CALL ERROR('GRDTE-02: Too small array RAM(MRAM)') C Too small array RAM(MRAM) to allocate both N1*N2*N3 real and C N1*N2*N3 imaginary output grid values. C If possible, increase dimension MRAM in include file C ram.inc. END IF IF(OUTR.NE.' ') THEN OPEN (LU1,FILE=OUTR,FORM='FORMATTED') END IF IF(OUTI.NE.' ') THEN OPEN (LU2,FILE=OUTI,FORM='FORMATTED') END IF PTM=N10*N20*N30 PT=0 C C Reading function parameters CALL RSEP3R('AR0',AR0,0.) CALL RSEP3R('AI0',AI0,0.) CALL RSEP3R('TR0',TR0,0.) CALL RSEP3R('TI0',TI0,0.) CALL RSEP3R('AR1',AR1(1),0.) CALL RSEP3R('AI1',AI1(1),0.) CALL RSEP3R('TR1',TR1(1),0.) CALL RSEP3R('TI1',TI1(1),0.) CALL RSEP3R('AR2',AR1(2),0.) CALL RSEP3R('AI2',AI1(2),0.) CALL RSEP3R('TR2',TR1(2),0.) CALL RSEP3R('TI2',TI1(2),0.) CALL RSEP3R('AR3',AR1(3),0.) CALL RSEP3R('AI3',AI1(3),0.) CALL RSEP3R('TR3',TR1(3),0.) CALL RSEP3R('TI3',TI1(3),0.) CALL RSEP3R('AR11',AR2(1,1),0.) CALL RSEP3R('AR12',AR2(1,2),0.) CALL RSEP3R('AR22',AR2(2,2),0.) CALL RSEP3R('AR13',AR2(1,3),0.) CALL RSEP3R('AR23',AR2(2,3),0.) CALL RSEP3R('AR33',AR2(3,3),0.) CALL RSEP3R('AI11',AI2(1,1),0.) CALL RSEP3R('AI12',AI2(1,2),0.) CALL RSEP3R('AI22',AI2(2,2),0.) CALL RSEP3R('AI13',AI2(1,3),0.) CALL RSEP3R('AI23',AI2(2,3),0.) CALL RSEP3R('AI33',AI2(3,3),0.) CALL RSEP3R('TR11',TR2(1,1),0.) CALL RSEP3R('TR12',TR2(1,2),0.) CALL RSEP3R('TR22',TR2(2,2),0.) CALL RSEP3R('TR13',TR2(1,3),0.) CALL RSEP3R('TR23',TR2(2,3),0.) CALL RSEP3R('TR33',TR2(3,3),0.) CALL RSEP3R('TI11',TI2(1,1),0.) CALL RSEP3R('TI12',TI2(1,2),0.) CALL RSEP3R('TI22',TI2(2,2),0.) CALL RSEP3R('TI13',TI2(1,3),0.) CALL RSEP3R('TI23',TI2(2,3),0.) CALL RSEP3R('TI33',TI2(3,3),0.) CALL RSEP3R('AR111',AR3(1,1,1),0.) CALL RSEP3R('AR112',AR3(1,1,2),0.) CALL RSEP3R('AR122',AR3(1,2,2),0.) CALL RSEP3R('AR222',AR3(2,2,2),0.) CALL RSEP3R('AR113',AR3(1,1,3),0.) CALL RSEP3R('AR123',AR3(1,2,3),0.) CALL RSEP3R('AR223',AR3(2,2,3),0.) CALL RSEP3R('AR133',AR3(1,3,3),0.) CALL RSEP3R('AR233',AR3(2,3,3),0.) CALL RSEP3R('AR333',AR3(3,3,3),0.) CALL RSEP3R('AI111',AI3(1,1,1),0.) CALL RSEP3R('AI112',AI3(1,1,2),0.) CALL RSEP3R('AI122',AI3(1,2,2),0.) CALL RSEP3R('AI222',AI3(2,2,2),0.) CALL RSEP3R('AI113',AI3(1,1,3),0.) CALL RSEP3R('AI123',AI3(1,2,3),0.) CALL RSEP3R('AI223',AI3(2,2,3),0.) CALL RSEP3R('AI133',AI3(1,3,3),0.) CALL RSEP3R('AI233',AI3(2,3,3),0.) CALL RSEP3R('AI333',AI3(3,3,3),0.) CALL RSEP3R('TR111',TR3(1,1,1),0.) CALL RSEP3R('TR112',TR3(1,1,2),0.) CALL RSEP3R('TR122',TR3(1,2,2),0.) CALL RSEP3R('TR222',TR3(2,2,2),0.) CALL RSEP3R('TR113',TR3(1,1,3),0.) CALL RSEP3R('TR123',TR3(1,2,3),0.) CALL RSEP3R('TR223',TR3(2,2,3),0.) CALL RSEP3R('TR133',TR3(1,3,3),0.) CALL RSEP3R('TR233',TR3(2,3,3),0.) CALL RSEP3R('TR333',TR3(3,3,3),0.) CALL RSEP3R('TI111',TI3(1,1,1),0.) CALL RSEP3R('TI112',TI3(1,1,2),0.) CALL RSEP3R('TI122',TI3(1,2,2),0.) CALL RSEP3R('TI222',TI3(2,2,2),0.) CALL RSEP3R('TI113',TI3(1,1,3),0.) CALL RSEP3R('TI123',TI3(1,2,3),0.) CALL RSEP3R('TI223',TI3(2,2,3),0.) CALL RSEP3R('TI133',TI3(1,3,3),0.) CALL RSEP3R('TI233',TI3(2,3,3),0.) CALL RSEP3R('TI333',TI3(3,3,3),0.) DO 2 J=1,3 DO 1 I=1,J IF(I.NE.J) THEN AR2(I,J)=2.*AR2(I,J) AI2(I,J)=2.*AI2(I,J) TR2(I,J)=2.*TR2(I,J) TI2(I,J)=2.*TI2(I,J) END IF 1 CONTINUE 2 CONTINUE DO 5 K=1,3 DO 4 J=1,K DO 3 I=1,J IF(I.NE.J.AND.J.NE.K.AND.K.NE.I) THEN AR3(I,J,K)=6.*AR3(I,J,K) AI3(I,J,K)=6.*AI3(I,J,K) TR3(I,J,K)=6.*TR3(I,J,K) TI3(I,J,K)=6.*TI3(I,J,K) ELSE IF(I.NE.J.OR.J.NE.K.OR.K.NE.I) THEN AR3(I,J,K)=2.*AR3(I,J,K) AI3(I,J,K)=2.*AI3(I,J,K) TR3(I,J,K)=2.*TR3(I,J,K) TI3(I,J,K)=2.*TI3(I,J,K) END IF 3 CONTINUE 4 CONTINUE 5 CONTINUE C C C Loop over Xi0: DO 120 IN30=1,N30 X0(3)=O30+D30*(IN30-1) DO 110 IN20=1,N20 X0(2)=O20+D20*(IN20-1) DO 100 IN10=1,N10 X0(1)=O10+D10*(IN10-1) L=0 C C Loop over Xi: DO 90 IN3=1,N3 X(3)=O3+D3*(IN3-1) DO 80 IN2=1,N2 X(2)=O2+D2*(IN2-1) DO 70 IN1=1,N1 X(1)=O1+D1*(IN1-1) AR=AR0 AI=AI0 TR=TR0 TI=TI0 C Calculating the functional value: DO 10 I=1,3 DX=X(I)-X0(I) AR=AR+AR1(I)*DX AI=AI+AI1(I)*DX TR=TR+TR1(I)*DX TI=TI+TI1(I)*DX 10 CONTINUE DO 30 J=1,3 DXJ=(X(J)-X0(J))/2. DO 20 I=1,J DXI=X(I)-X0(I) DX=DXI*DXJ AR=AR+AR2(I,J)*DX AI=AI+AI2(I,J)*DX TR=TR+TR2(I,J)*DX TI=TI+TI2(I,J)*DX 20 CONTINUE 30 CONTINUE DO 60 K=1,3 DXK=(X(K)-X0(K))/6. DO 50 J=1,K DXJ=X(J)-X0(J) DO 40 I=1,J DXI=X(I)-X0(I) DX=DXI*DXJ*DXK AR=AR+AR3(I,J,K)*DX AI=AI+AI3(I,J,K)*DX TR=TR+TR3(I,J,K)*DX TI=TI+TI3(I,J,K)*DX 40 CONTINUE 50 CONTINUE 60 CONTINUE E=EXP(TR) C=COS(TI) S=SIN(TI) FR=(AR*C-AI*S)*E FI=(AR*S+AI*C)*E L=L+1 RAM(L)=FR RAM(N1N2N3+L)=FI C End of calculation at one gridpoint C 70 CONTINUE 80 CONTINUE 90 CONTINUE C C Writing output grid values IF(OUTR.NE.' ') THEN CALL WARRAY(LU1,' ','FORMATTED',.FALSE.,0.,.FALSE.,0., * N1N2N3,RAM) END IF IF(OUTI.NE.' ') THEN CALL WARRAY(LU2,' ','FORMATTED',.FALSE.,0.,.FALSE.,0., * N1N2N3,RAM(N1N2N3+1)) END IF C C End of loop over Xi C PT=PT+1 WRITE(*,'(''+GRDTE: '',I16,'' loops over Xi0 of'',I9)') PT,PTM 100 CONTINUE 110 CONTINUE 120 CONTINUE C C End of loop over Xi0 C IF(OUTR.NE.' ') THEN CLOSE(LU1) END IF IF(OUTI.NE.' ') THEN CLOSE(LU2) END IF C WRITE(*,'(A)') '+GRDTE: Done.' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= C