	dimension dlat(150),dlon(150),dist(150),flag(150),rat(150)
	dimension d1(150),As(150),Bs(150),Cs(150),rad(150)
	implicit real*8(a-h,o-z)

	dimension u(3,150),anu(3,150),p(3,150),sv(3,150),sh(3,150)
	dimension fsp(150),fsv(150),fsh(150)

c    hypocentral distance
c
c		
	open(13, file="outdata")
	open(14, file="outdata2")
	open(15, file="outdata3")
c	open(16, file="outdata4")
        OPEN(10, file="indata1",status='old')
        OPEN(11, file="indata2",status='old')
	open(12, file="indata3",status='old')
c   震源情報
c     緯度、経度
c     震源の深さ（Ｋｍ）  
c     strike-angle phi ---- strike direction from x to y
c     dip-angle    dlt ---- dip 
c     rake-angle    alm ---- rake

        read(10,*)alat
	read(10,*)alon
	read(10,*)adepth
	read(10,*)phi
        read(10,*)dlt
	read(10,*)alm

c   横浜市内150点の緯度、経度    
c
        pi=4*atan(1.0)
	do i=1,150
	read(11,*)dlat(i),dlon(i)
	enddo
C	write(6,1001)(dlat(i),dlon(i),i=1,150)
	    read(12,*)(flag(j),j=1,150)

	Ae=(cos(alat*pi/180))*(cos(alon*pi/180))
	Be=(cos(alat*pi/180))*(sin(alon*pi/180))
	Ce=sin(alat*pi/180)

 	do 11 i=1,150
	As(i)=(cos(dlat(i)*pi/180))*(cos(dlon(i)*pi/180))
	Bs(i)=(cos(dlat(i)*pi/180))*(sin(dlon(i)*pi/180))
	Cs(i)=sin(dlat(i)*pi/180)
11	continue
 
	do 12 j=1,150
12      d1(j)=(sqrt((Ae-As(j))**2+(Be-Bs(j))**2+(Ce-Cs(j))**2))*6371
	do 20 k=1,150
20	    dist(k)=sqrt(d1(k)**2+adepth**2)
	
c   ラディエーションパターン 係数の計算 
c     radiation pattern   (see Aki & Richards, 1980 )
c       originally coded by S. Sekiguchi (NIED)
c
c     input data (in radian)
c     fault parameters (p.106)
c       phi ---- strike direction from x to y
c       dlt ---- dip
c       alm ---- rake
c     ray parameter    (p.114)
c       tht ---- colatitude
c       psi ---- longitude
c     output
c       fp,fsv,fsh
c      dimension u(3),anu(3),p(3),sv(3),sh(3)
c
c     slip direction
       do 200 j=1,150

       u(1,j)=dcos(alm)*dcos(phi)+dcos(dlt)*dsin(alm)*dsin(phi)
       u(2,j)=dcos(alm)*dsin(phi)-dcos(dlt)*dsin(alm)*dcos(phi)
       u(3,j)=-dsin(alm)*dsin(dlt)
c
c     fault normal
      anu(1,j)=-dsin(dlt)*dsin(phi)
      anu(2,j)=dsin(dlt)*dcos(phi)
      anu(3,j)=-dcos(dlt)
c
c     p -wave direction
      p(1,j)  =dsin(dlat(j))*dcos(dlon(j))
      p(2,j)  =dsin(dlat(j))*dsin(dlon(j))
      p(3,j)  =dcos(dlat(j))
c
c     sv-wave direction
      sv(1,j) =dcos(dlat(j))*dcos(dlon(j))
      sv(2,j) =dcos(dlat(j))*dsin(dlon(j))
      sv(3,j) =-dsin(dlat(j))
c
c     sh-wave direction
      sh(1,j) =-dsin(dlon(j))
      sh(2,j) =dcos(dlon(j))
      sh(3,j) =0.
c
c     radiation pattern
c
      fsp(j)= 2.*( p(1,j)*anu(1,j)+ p(2,j)*anu(2,j)+ p(3,j)*anu(3,j))
     &      *( p(1,j)*  u(1,j)+ p(2,j)*  u(2,j)+ p(3,j)*  u(3,j))
      fsv(j)=   ( p(1,j)*anu(1,j)+ p(2,j)*anu(2,j)+ p(3,j)*anu(3,j))
     &      *(sv(1,j)*  u(1,j)+sv(2,j)*  u(2,j)+sv(3,j)*  u(3,j))
     &      +( p(1,j)*  u(1,j)+ p(2,j)*  u(2,j)+ p(3,j)*  u(3,j))
     &      *(sv(1,j)*anu(1,j)+sv(2,j)*anu(2,j)+sv(3,j)*anu(3,j))
      fsh(j)=   ( p(1,j)*anu(1,j)+ p(2,j)*anu(2,j)+ p(3,j)*anu(3,j))
     &      *(sh(1,j)*  u(1,j)+sh(2,j)*  u(2,j)+sh(3,j)*  u(3,j))
     &      +( p(1,j)*  u(1,j)+ p(2,j)*  u(2,j)+ p(3,j)*  u(3,j))
     &      *(sh(1,j)*anu(1,j)+sh(2,j)*anu(2,j)+sh(3,j)*anu(3,j))
c
      rad(j)=sqrt(fsv(j)**2+fsh(j)**2)
	rat(j)=(abs(fsv(j)))/(abs(fsh(j)))
 200	continue
            a1=0
	    a2=0
       do 1 k=1,150
           a1=a1+rad(k)
	   a2=a2+rad(k)**2
 1         continue
           s1=a1/150
	   st1=sqrt((150*a2-a1**2)/(150*149))
c
c
C    flag = 1 : recorded
c    flag = 0 : not recorded
       do 31 i=1,150
                if (flag(i).eq.1) go to 32
                if (flag(i).eq.0) go to 31
 32		    write(13,1000)dist(i)
 31		continue

	   write(6,*)dist(1),rad(1),s1,st1,phi,dlt,alm
       do 41 i=1,150
                if (flag(i).eq.1) go to 42
                if (flag(i).eq.0) go to 41
 42		write(14,1001)rad(i)
 41	    continue

       do 51 i=1,150
                if (flag(i).eq.1) go to 52
                if (flag(i).eq.0) go to 51
 52		write(15,1000)rat(i)
 51	    continue
1000    format(1f10.5)

1001	format(2f8.3)
        stop
	end


	 

