C     ************
C     *  CALF3N  *
C     ************
C
      PROGRAM CALF3N
C
      PARAMETER (ND1=100,ND2=50)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 F(ND1),S,AA,FAULTN
      DIMENSION FXX(ND2),FYY(ND2),DOFF(ND2,ND2),T0(ND2,ND2)
      CHARACTER FNAME*20
      EXTERNAL FAULTN
C
      COMMON /FPOSI/ X0,Y0,Z0
     *       /SPOSI/ XS,YS,ZS
     *       /CAL01/ DX,DY,DZ,XN,YN,ZN
     *       /CAL02/ SINT,COST,SIND,COSD
     *       /PCONS/ PMU,PLAM,VR,PAI
     *       /WAVEN/ AK,AK2,AK3,BK,TBK1,TBK2
     *       /OMEGA/ OM
C
      PAI =4.D0*DATAN(1.D0)
C
      WRITE(*,6000)
      READ(*,'(A)')FNAME
C
      OPEN(4,FILE=FNAME)
      READ (4,4010)VS,VP
      WRITE(*,6010)VS,VP
      READ (4,4020)DENS
      WRITE(*,6020)DENS
      READ (4,4030)X0,Y0,Z0
      WRITE(*,6030)X0,Y0,Z0
      READ (4,4040)FL,FW
      WRITE(*,6040)FL,FW
      READ (4,4050)DELTA,THETA,RAKE
      WRITE(*,6050)DELTA,THETA,RAKE
      DELTA=DELTA/180.D0*PAI
      THETA=THETA/180.D0*PAI
      RAKE =RAKE /180.D0*PAI
      READ (4,4060)NFX,NFY
      WRITE(*,6060)NFX,NFY
      IX=1
      IY=1
      READ (4,4070)IRS
      WRITE(*,6070)IRS
      IF(IRS.EQ.1) IY=0
      IF(IRS.EQ.2) IX=0
      IF(IRS.NE.-1)THEN
        READ (4,4080)DOFF0
        WRITE(*,6080)DOFF0
        DO 9 I=1,NFX
          DO 8 J=1,NFY
            DOFF(I,J)=DOFF0
    8     CONTINUE
    9   CONTINUE
      ELSE
        DO 15 I=1,NFX
          DO 10 J=1,NFY
            READ (4,4085)DOFF(I,J),T0(I,J)
            WRITE(*,6085)I,J,DOFF(I,J),T0(I,J)
   10     CONTINUE
   15   CONTINUE
      END IF
      READ (4,4090)VR,TAU
      WRITE(*,6090)VR,TAU
      READ (4,4100)NFX0,NFY0
      WRITE(*,6100)NFX0,NFY0
      READ (4,4110)XS,YS,ZS
      WRITE(*,6110)XS,YS,ZS
      READ (4,4120)DOM
      WRITE(*,6120)DOM
      READ (4,4130)NO
      WRITE(*,6130)NO
      CLOSE(4)
C
 4010 FORMAT(2F10.3)
 4020 FORMAT(F10.3)
 4030 FORMAT(3F10.2)
 4040 FORMAT(2F10.3)
 4050 FORMAT(3F10.3)
 4060 FORMAT(2I10)
 4070 FORMAT(I10)
 4080 FORMAT(F10.3)
 4085 FORMAT(2F10.3)
 4090 FORMAT(2F10.3)
 4100 FORMAT(2I10)
 4110 FORMAT(3F10.2)
 4120 FORMAT(F10.3)
 4130 FORMAT(I10) 
C
 6000 FORMAT(//'**********  INPUT FILE   **********'
     *       //' INPUT FILE NAME...........',\)
 6010 FORMAT( /'*******  PHYSICAL CONSTANT  *******'
     *       //' VS,VP.....................',2F10.3)
 6020 FORMAT(  ' DENSITY...................',F10.3)
 6030 FORMAT( /'**********  FAULT DATA  ***********'
     *       //' POSITION (X,Y,Z)..........',3F10.2)
 6040 FORMAT(  ' LENGTH,WIDTH..............',2F10.2)
 6050 FORMAT(  ' ANGLE (DELTA,THETA,RAKE)..',3F10.2)
 6060 FORMAT(  ' NFX,NFY...................',2I10)
 6070 FORMAT(  ' RUPTURE STYLE (-1,0,1,2)..',I10)
 6080 FORMAT(  ' DISP......................',F10.3)
 6085 FORMAT(  ' DISP,TIMELAG (',I2,',',I2,')......',2F10.3)
 6090 FORMAT(  ' RUPTURE VEL,RISETIME......',2F10.3)
 6100 FORMAT(  ' POSITION OF R.S.P (NX,NY).',2I10)
 6110 FORMAT( /'********  OBSERVATION DATA  *******'
     *       //' POSITION (X,Y,Z)..........',3F10.2)
 6120 FORMAT( /'**********  FOURIER DATA  *********'
     *       //' DELTA OMEGA...............',F10.5)
 6130 FORMAT(  ' NUMBER OF OMEGA...........',I10/)
 6140 FORMAT(//'**********  OUTPUT FILE  **********'
     *       //' OUTPUT FILE NAME..........',\)
C
      PMU =VS*VS*DENS
      PLAM=VP*VP*DENS-2.D0*PMU
C
      SINT=DSIN(THETA)
      COST=DCOS(THETA)
      SIND=DSIN(DELTA)
      COSD=DCOS(DELTA)
      SINR=DSIN( RAKE)
      COSR=DCOS( RAKE)
C
      XN  =-SINT*SIND
      YN  = COST*SIND
      ZN  =-     COSD
C
      DX  =COSR*COST-SINR*COSD*SINT
      DY  =COSR*SINT+SINR*COSD*COST
      DZ  =          SINR*SIND
C
      EFX =FL/DFLOAT(NFX)
      EFY =FW/DFLOAT(NFY)
C
      DO 30 J=1,NFX
        FXX(J)=EFX*(DFLOAT(J)-0.5D0)
   30 CONTINUE
      DO 40 K=1,NFY
        FYY(K)=EFY*(DFLOAT(K)-0.5D0)
   40 CONTINUE
C
      FX0 =FXX(NFX0)
      FY0 =FYY(NFY0)
C
      IF(IRS.EQ.-1) GOTO 60
      DO 50 J=1,NFX
        DO 50 K=1,NFY
          T0(J,K)=DSQRT((FXX(J)-FX0)*(FXX(J)-FX0)*IX
     *                 +(FYY(K)-FY0)*(FYY(K)-FY0)*IY)/VR
   50 CONTINUE
C
   60 DO 150 I=1,NO
C
        S=(0.D0,0.D0)
C
        OM  =DOM*DFLOAT(I)
        AK  =OM/VP
        AK2 =AK*AK
        AK3 =AK2*AK
        BK  =OM/VS
        TBK1=1.D0/BK
        TBK2=TBK1/BK
C
        OMT=OM*TAU*0.5D0
        IF(IRS.EQ.2) THEN 
          OML=OM*EFY*0.5D0/VR
        ELSE
          OML=OM*EFX*0.5D0/VR
        END IF
C
        DO 100 J=1,NFX
          DO 100 K=1,NFY
C
            S=S+FAULTN(FXX(J),FYY(K))*CDEXP(DCMPLX(0.D0,-OM*T0(J,K)))
     *        *DOFF(J,K)
C
  100   CONTINUE
C
        AA=EFX*EFY*DSIN(OMT)/OMT*DSIN(OML)/OML
     *    *DCMPLX(0.D0,-1.D0)/OM*CDEXP(DCMPLX(0.D0,-OMT-OML))
C
        F(I)=2.D0*AA*S
C
  150 CONTINUE
C
      WRITE(*,6140)
      READ (*,'(A)')FNAME
      OPEN(7,FILE=FNAME)
      WRITE(7,7000)DOM,NO
      WRITE(7,7010)1
      DO 200 I=1,NO,5
        WRITE(7,7020)(F(I+K),K=0,4)
  200 CONTINUE
      CLOSE(7)
C
 7000 FORMAT(F15.7,I15)
 7010 FORMAT(I15///)
 7020 FORMAT(10F8.5)
C
      STOP
      END
 