C C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CTITLESADM3DX - AMPLITUDE DISTRIBUTION FOR STACK CA CA DESIGNER JAMES SUN CA AUTHOR JAMES SUN CA LANGUAGE FORTRAN 77 CA SYSTEM CRAY/IBM CA WRITTEN 12/02/87 C REVISED 08/12/88 CAD MODIFIED TO MEET EDP STANDARDS C REVISED 02/06/90 JCS MODIFIED TO USE BEASLEY'S WEIGHT C REVISED 08/06/90 JCS ADD USER INPUT PARAMETER 'MXAPER' C FOR MAXIMUM CROSSLINE APERTURE C REVISED 12/17/91 JJC MODIFIED TO MEET SPARC STANDARDS. C C C PURPOSE: : DISTRIBUTE THE ENERGY TO OTHER LOCATIONS ACCORDING C TO DMO IN CASE OF VELOCITY OPTION. CA CA CALLING PROCEDURE: CA SUBROUTINE SADM3DX(HX,HY,IX0,IY0,IFLV,A,TSQRT,FCUT, CA + T0MAX,STRCH,WGHT,P,Q,GAMMA,WK1,WK2,WK3,WK4,OTR, CA + NT,NX,NY,IABORT,WORKN,WORKQ) CA C CALLING ARGUMENTS CA CA INPUT: HX : HALF OF X OFFSET R4 CA INPUT: HY : HALF OF Y OFFSET R4 CA INPUT: IX0 : CDP INDEX OF INPUT TRACE I4 CA INPUT: IY0 : LINE INDEX OF INPUT TRACE I4 CA INPUT: IFLV : FIRST LIVE VALUE I4 CA IN/OUT A : R4 CA INPUT: TSQRT : WEIGHTING FUNCTION (ARRAY) R4 CA INPUT: FCUT : HIGH CUT FREQUENCH ARRAY (ARRAY) R4 CA INPUT: T0MAX : MAXIMUM TIME FOR DMO OPERATION (ARRAY) R4 CA INPUT: STRCH : DMO TIME STRETCH ARRAY R4 CA INPUT: WGHT : DMO WEIGHT ARRAY R4 CA INPUT: P : INPUT TRACE AFTER FFT CONVERSION R4 CA INPUT: Q : WORK ARRAY R4 CA INPUT: OTR : INPUT TRACE R4 CA INPUT: NT : NUMBER OF SAMPLES I4 CA INPUT: NX : NUMBER OF OUTPUT GRIDS IN X I4 CA INPUT: NY : NUMBER OF OUTPUT GRIDS IN Y I4 CA IN/OUT KPWRKS : WORK FILE POINTER FOR DIRECT ACCESS I4 CA INPUT: IABORT : ABORT FLAG I4 C SUBROUTINE SADM3DX(HX,HY,IX0,IY0,IFLV,A,TSQRT,FCUT, + T0MAX,STRCH,WGHT,P,Q,GAMMA,WK1,WK2,WK3,WK4,OTR, + NT,NX,NY,IABORT,WORKN,WORKQ) C C IMPLICIT INTEGER (A-Z) C C REAL A REAL AHI REAL AHX REAL AHY REAL AIT0 REAL AITN REAL BITN REAL BK REAL AK REAL AMAX1 REAL AMIN1 REAL AMX REAL AMY REAL ANF REAL ANT1 REAL AXDST REAL DF REAL DRSQ REAL DX REAL DY REAL FCUT REAL FIX0 REAL FIY0 REAL FLOAT REAL HSQ REAL HX REAL HY REAL OTR REAL P REAL PA REAL PB REAL PI REAL PID2 REAL Q REAL STRCH REAL T0MAX REAL TEMP REAL TSQRT REAL WGHT REAL GAMMA REAL WK1 REAL WK2 REAL WK3 REAL WK4 REAL WORKN REAL WORKQ C C COMMON /USER/ SLOCAL(50),ULOCAL(220) C C EQUIVALENCE (ANT1 , ULOCAL( 8)) EQUIVALENCE (DF , ULOCAL( 18)) EQUIVALENCE (DX , ULOCAL( 25)) EQUIVALENCE (DY , ULOCAL( 26)) EQUIVALENCE (ANF , ULOCAL( 29)) EQUIVALENCE (NF1 , ULOCAL( 51)) EQUIVALENCE (IPR , ULOCAL( 56)) EQUIVALENCE (NW2 , ULOCAL(119)) EQUIVALENCE (IVFLAG , ULOCAL(120)) EQUIVALENCE (MXAPER , ULOCAL(121)) EQUIVALENCE (MXSUM , ULOCAL(122)) EQUIVALENCE (MYSUM , ULOCAL(124)) EQUIVALENCE (NF , ULOCAL(126)) EQUIVALENCE (NO , ULOCAL(129)) C C EQUIVALENCE (NT , ULOCAL(132)) C EQUIVALENCE (NWD21 , ULOCAL(139)) C C EQUIVALENCE (NX , ULOCAL(140)) C EQUIVALENCE (NY , ULOCAL(144)) C EQUIVALENCE (PI , ULOCAL(147)) EQUIVALENCE (YES , ULOCAL(176)) EQUIVALENCE (FIX0 , ULOCAL(181)) EQUIVALENCE (FIY0 , ULOCAL(182)) C C DIMENSION REQUIREMENT: C TSQRT: NT C FCUT : 1001 C T0MAX: 1001 C STRCH: 1001 C WGHT : 1001 C P : (NT )*(NF+1) C Q : (NW+2) C OTR : NT C C C INPUT TRACE: Q C DIMENSION A(NT,NX,NY) C C DIMENSION TSQRT(1) DIMENSION FCUT(1) DIMENSION T0MAX(1) DIMENSION STRCH(1) DIMENSION WGHT(1) DIMENSION P(NT,1) DIMENSION Q(1) DIMENSION OTR(1) DIMENSION GAMMA(1) DIMENSION WK1(NT,1) DIMENSION WK2(1) DIMENSION WK3(1) DIMENSION WK4(1) DIMENSION WORKN(1) DIMENSION WORKQ(1) C C INTEGER ORTN,CDPN,CDPT,TICD,XDST,NS,SI,SSP,FN,THL,LNNO C C INTEGER YES,NO C C C C COMMON /HEAD/ ORTN,CDPN,CDPT,TICD,XDST,NS,SI,SSP,FN,THL,LNNO C C DATA PID2/1.570796327/ C C C C THE LOGIC USED TO DISTRIBUTE THE ENERGIES IS DESCRIBED GRAPHICALLY C AS FOLLOWS: C C C C C -A 0 +A C ____________________________________________ C * | * C * | * C * | * C * * * C | C * |T * C |I C |M C |E C * || * C || C |V C * | * C | C * | * C * C C C ABOVE GRAPH SHOWS A TRACE AND THE '*' SHOW THE THE RELATION C BETWEEN ENERGY AT 0 AND ENERGY AT ANY DISTANCE X AWAY. C C C WHEN THE SAME GRAPH WHEN VIEWED IN XY PLANE THREE POSSIBILITIES C MAY OCCUR. ONE WHEN MX = MY = 0 AND OTHER TWO WHEN MX > MY C AND MX <= MY. C C C WHEN MX > MY WHEN MX <= MY C C C B| . B| C | | . C 0 | . | C | | | C | | . | C M | | . C Y | . | C | | C | . 0 | C | | | . C X-------------------------- | | C -A 0 A | C <---------MX----------------> M | C -B Y | . C | C | C | C | . C | C | C | C |________________________ C -A X A C <------------MX------------> C -B C C C ABOVE FIGURES SHOW THE SAME DISTANCE FROM -A TO A ALONG X AXIS. C THE -A TO 0 IS SYMMETRIC AND THEREFORE NOT SHOWN HERE. C DISTANCE -A TO A IS MX AND DISTANCE -B TO B IS MY. C A TRACE COMES IN AT 0,0 IN ABOVE FIGURE WHICH IS IX0,IY0 IN THE C FOLLOWING PROGRAM. IN THE PROGRAM PT A AND -A ARE THE IXMX AND C AND IXMN AND PTS B AND -B ARE IYMX AND IYMN RESPECTIVELY. C THE LOCATION OF ANY POINT '.' ON THE GRAPH IS CALCULATED GIVEN C AX0,AY0, AND KNOWING THE SLOPE AND CONSTANT OF THE STRAIGHT LINE. C IX AND IY BECOME THE X AND Y COMPONENTS OF ANY POINT '.'. THE C ENERGY ALONG ANY POINT '.' IS DISTRIBUTED BASED ON WHERE THE C POINT LIES ON THE GRID AS SHOWN BELOW C C C |___|___|___|___|___DY3 THE ENERGY IS DISTRIBUTED C | | | | | BASED ON WHERE A POINT '.' C | | | | Z. LIES. FOR POINT X THE ENERGY C |___|___|___|___|___DY2 WILL DISTRIBUTED TO THE LEVEL C | | Y. | | OF (DX1,DY1). FOR PT Y TO THE C X. | | | | LEVEL OF (DX3,DY2) AND SO ON. C |___|___|___|___|___DY1 THE SAME DISTRIBUTION IS DONE C DX1 DX2 DX3 DX4 DX5 BUT ALONG HORIZONTAL FOR THE C SITUATION WHEN MY >= MX. C C C IFLAG=0 C C IF(DX.NE.0.) AHX=HX/DX IF(DY.NE.0.) AHY=HY/DY AMX=ABS(AHX) AMY=ABS(AHY) C MX=AMX C MY=AMY MX=NINT(0.6*AMX) MY=NINT(0.6*AMY) C C WRITE(IPR,7777) HX,DX,AHX,AMX,MX,NINT(AMX) C7777 FORMAT(' HX,DX,AHX,AMX =',4E15.8,' MX,NINT(AMX) =',2I5) C HSQ=HX*HX+HY*HY AXDST=SQRT(HSQ) C C AHI=1./AMAX1(AXDST,1.) C C WRITE(IPR,9875) LNNO, CDPN,NINT(2.*AXDST) C9875 FORMAT(' LNNO,CDPN,XDST =',3I5) C C CALL SCOPY(NT,OTR,1,Q(MOP),1) CALL SCOPY(NT,OTR,1,Q( 1),1) C C=================================================================== C IF(MX.EQ.0 .AND. MY.EQ.0) THEN C C=================================================================== C IF(IX0.LT.1 .OR. IX0.GT.NX) GO TO 120 IF(IY0.LT.1 .OR. IY0.GT.NY) GO TO 120 C C IF(IFLAG.EQ.0) THEN CALL SADM3DL(P,Q,NT,WORKN,WORKQ) IFLAG=1 ENDIF C C DO 100 IT=IFLV,NT A(IT,IX0,IY0)=A(IT,IX0,IY0)+P(IT,NF1) 100 CONTINUE C C IF(IXO.EQ.1 .AND. IYO.EQ.1) THEN C WRITE(IPR,9780) LNNO, CDPN,NINT(2.*AXDST),IFLV C9780 FORMAT(' *0* LNNO,CDPN,XDST =',3I5,' IFLV =',I5) C CALL PTST1R('A ',A(1,IX0,IY0),250,IPR) C ENDIF C 120 CONTINUE C C=================================================================== C ELSE IF(MX.GE.MY) THEN C C=================================================================== C IXMN=MAX0(IX0-MX, 1) IXMX=MIN0(IX0+MX, NX) C C DO 200 IX=IXMN,IXMX C C IY=(AHY/AHX)*FLOAT(IX-IX0)+0.5+FLOAT(IY0) C IY=NINT((AHY/AHX)*(FLOAT(IX)-FIX0)+FIY0) IF(IY.GT.NY .OR. IY.LT.1 .OR. ABS(IY-IY0).GT.MY) GO TO 200 IF(IABS(IY-IY0+1).GT.MXAPER) GO TO 200 C C IF(IFLAG.EQ.0) THEN C JT=MOP JT=1 DO 140 IT=1,NT Q(JT)=Q(JT)*ABS(DX)/(ABS(DX)+ABS(HX)*TSQRT(IT)) 140 JT=JT+1 C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C WRITE(IPR,9881) LNNO,CDPN,NINT(2.*AXDST),DX,HX C9881 FORMAT(' *X* LNNO,ICDPN,XDST =',3I5,' DX,HX =',2E12.5) C CALL PTST1R('OTR ',OTR ,NT ,IPR) C CALL PTST1R('TSQT',TSQRT ,NT ,IPR) C CALL PTST1R('Q ',Q(MOP) ,NT ,IPR) C ENDIF C ENDIF CALL SADM3DL(P,Q,NT,WORKN,WORKQ) IFLAG=1 ENDIF C C DRSQ=((FLOAT(IX)-FIX0)*DX)**2+((FLOAT(IY)-FIY0)*DY)**2 IKK=NINT(DRSQ/HSQ*1000.)+1 IF(IKK.GT.1001) GO TO 200 JFLV=MAX0(IFLV,2) C C TEMP=AMX*FCUT(IKK) C C IT0MX= MIN1(AXDST*T0MAX(IKK),ANT1/STRCH(IKK))+1 C C IF(IVFLAG.EQ.NO) THEN C C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C WRITE(IPR,9781) LNNO,CDPN,NINT(2.*AXDST),JFLV,IT0MX,IFLV C9781 FORMAT(' *X* LNNO,ICDPN,XDST =',3I5,' JFLV,IT0MX =',2I5, C + ' IFLV =',I5) C CALL PTST1R('Q ',Q ,250,IPR) C CALL PTST2R('P ',P,NT,NF1,250,1,IPR) C CALL PTST1R('A ',A(1,IX,IY),250,IPR) C ENDIF C ENDIF C DO 160 IT0=JFLV,IT0MX AIT0=IT0-1 AK=AMIN1(TEMP/AIT0+1.,ANF) KK=AK AK=AK-KK AITN=AIT0*STRCH(IKK)+1. ITN=AITN AITN=AITN-ITN PA=P(ITN ,KK)+AK*(P(ITN ,KK+1)-P(ITN ,KK)) PB=P(ITN+1,KK)+AK*(P(ITN+1,KK+1)-P(ITN+1,KK)) A(IT0,IX,IY)=A(IT0,IX,IY)+WGHT(IKK)*(PA+AITN*(PB-PA)) 160 CONTINUE C C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C CALL PTST1R('A ',A(1,IX,IY),250,IPR) C ENDIF C ENDIF C ELSE C C DO 180 IT0=JFLV,IT0MX AIT0=IT0-1 AK=AMIN1(TEMP/AIT0+1.,ANF) KK=AK AK=AK-KK ITN=WK1(IT0,IKK) AITN=WK1(IT0,IKK)-ITN PA=P(ITN ,KK)+AK*(P(ITN ,KK+1)-P(ITN ,KK)) PB=P(ITN+1,KK)+AK*(P(ITN+1,KK+1)-P(ITN+1,KK)) A(IT0,IX,IY)=A(IT0,IX,IY)+WGHT(IKK)*(PA+AITN*(PB-PA)) 180 CONTINUE C C ENDIF C C 200 CONTINUE C C C C=================================================================== C ELSE C C=================================================================== C IYMN=MAX0(IY0-MY, 1) IYMX=MIN0(IY0+MY, NY) C C WRITE(IPR,7971) LNNO,CDPN,IY0,MY,NY,IYMN,IYMX C7971 FORMAT(' LNNO,CDPN =',2I5,' IY0,MY,NY =',3I5,' IYMN,IYMX =',2I5) C DO 280 IY=IYMN,IYMX IF(IABS(IY-IY0+1).GT.MXAPER) GO TO 280 C C IX=(AHX/AHY)*FLOAT(IY-IY0)+0.5+FLOAT(IX0) C IX=NINT((AHX/AHY)*(FLOAT(IY)-FIY0)+FIX0) C WRITE(IPR,7972) IY,IX,IX0,MX,AHX,AHY,FIY0,FIX0 C7972 FORMAT(' IY,IX,IX0,MX =',4I5,' AHX,AHY,FIY0,FIX0 =',4E12.5) IF(IX.GT.NX .OR. IX.LT.1 .OR. ABS(IX-IX0).GT.MX) GO TO 280 C C IF(IFLAG.EQ.0) THEN C JT=MOP JT=1 DO 220 IT=1,NT Q(JT)=Q(JT)*ABS(DY)/(ABS(DY)+ABS(HY)*TSQRT(IT)) 220 JT=JT+1 C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C WRITE(IPR,9882) LNNO,CDPN,NINT(2.*AXDST),DY,HY C9882 FORMAT(' *Y* LNNO,ICDPN,XDST =',3I5,' DY,HY =',2E12.5) C CALL PTST1R('OTR ',OTR ,NT ,IPR) C CALL PTST1R('TSQT',TSQRT ,NT ,IPR) C CALL PTST1R('Q ',Q ,NT ,IPR) C ENDIF C ENDIF CALL SADM3DL(P,Q,NT,WORKN,WORKQ) IFLAG=1 ENDIF C C DRSQ=((FLOAT(IX)-FIX0)*DX)**2+((FLOAT(IY)-FIY0)*DY)**2 IKK=NINT(DRSQ/HSQ*1000.)+1 IF(IKK.GT.1001) GO TO 280 JFLV=MAX0(IFLV,2) C C C TEMP=AMX*FCUT(IKK) TEMP=AMY*FCUT(IKK) C C IT0MX= MIN1(AXDST*T0MAX(IKK),ANT1/STRCH(IKK))+1 C C IF(IVFLAG.EQ.NO) THEN C C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C WRITE(IPR,9782) LNNO,CDPN,NINT(2.*AXDST),JFLV,IT0MX,IFLV C9782 FORMAT(' *Y* LNNO,ICDPN,XDST =',3I5,' JFLV,IT0MX =',2I5, C + ' IFLV =',I5) C CALL PTST1R('Q ',Q ,250,IPR) C CALL PTST2R('P ',P,NT,NF1,250,1,IPR) C CALL PTST1R('A ',A(1,IX,IY),250,IPR) C ENDIF C ENDIF C DO 240 IT0=JFLV,IT0MX AIT0=IT0-1 AK=AMIN1(TEMP/AIT0+1.,ANF) KK=AK AK=AK-KK AITN=AIT0*STRCH(IKK)+1. ITN=AITN AITN=AITN-ITN PA=P(ITN ,KK)+AK*(P(ITN ,KK+1)-P(ITN ,KK)) PB=P(ITN+1,KK)+AK*(P(ITN+1,KK+1)-P(ITN+1,KK)) A(IT0,IX,IY)=A(IT0,IX,IY)+WGHT(IKK)*(PA+AITN*(PB-PA)) 240 CONTINUE C C IF(IX.EQ.1 .AND. IY.EQ.1) THEN C IF(LNNO.GE.75) THEN C CALL PTST1R('A ',A(1,IX,IY),250,IPR) C ENDIF C ENDIF C ELSE C C DO 260 IT0=JFLV,IT0MX AIT0=IT0-1 AK=AMIN1(TEMP/AIT0+1.,ANF) KK=AK AK=AK-KK ITN=WK1(IT0,IKK) AITN=WK1(IT0,IKK)-ITN PA=P(ITN ,KK)+AK*(P(ITN ,KK+1)-P(ITN ,KK)) PB=P(ITN+1,KK)+AK*(P(ITN+1,KK+1)-P(ITN+1,KK)) A(IT0,IX,IY)=A(IT0,IX,IY)+WGHT(IKK)*(PA+AITN*(PB-PA)) 260 CONTINUE C C ENDIF C CSSD CALL PUTWA(KPWRKS,OTR,JSEQDA,NT,IERR) C 280 CONTINUE C C=================================================================== C ENDIF C C=================================================================== C RETURN END