      parameter (MOB=150,MEV=21)

      implicit real*8(a-h,o-z)

      integer flag(MOB)
      dimension dlat(MOB),dlon(MOB),dist(MOB),rat(MOB)
      dimension dlat1(MOB),dlon1(MOB)
      dimension d1(MOB),As(MOB),Bs(MOB),Cs(MOB),rad(MOB)

      dimension u(3,MOB),anu(3,MOB),p(3,MOB),sv(3,MOB),sh(3,MOB)
      dimension fsp(MOB),fsv(MOB),fsh(MOB)
      character INFILE1(MEV)*100,INFILE2(MEV)*100,INFILE3(MEV)*100,
     &          OTFILE1(MEV)*100,OTFILE2(MEV)*100,OTFILE3(MEV)*100,a*10
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      NEV=MEV

      PI=dacos(-1.d0)

      a='0123456789'

      open(16,file='rad16.inp')

      do 100 iEV=1,NEV

cccc ループから外すもの（対象外） cccc
         if (iEV.eq.12) then
            goto 100
         else if (iEV.eq.13) then
            goto 100
         else if (iEV.eq.14) then
            goto 100
         else if (iEV.eq.15) then
            goto 100
         else if (iEV.eq.16) then
            goto 100
         endif
cccccccccccccccccccccccccccccccccccccc

         if (iEV.lt.10) then
            INFILE1(iEV)='sou'//a(iEV+1:iEV+1)//'.inp'
            INFILE2(iEV)='tmp.inp'
            INFILE3(iEV)='f'//a(iEV+1:iEV+1)//'.inp'
            OTFILE1(iEV)='d'//a(iEV+1:iEV+1)//'.inp'
            OTFILE2(iEV)='rad'//a(iEV+1:iEV+1)//'.inp'
            OTFILE3(iEV)='rat'//a(iEV+1:iEV+1)//'.inp'
         else if (iEV.ge.10.and.iEV.lt.20) then
            INFILE1(iEV)='sou1'//a(iEV-9:iEV-9)//'.inp'
            INFILE2(iEV)='tmp.inp'
            INFILE3(iEV)='f1'//a(iEV-9:iEV-9)//'.inp'
            OTFILE1(iEV)='d1'//a(iEV-9:iEV-9)//'.inp'
            OTFILE2(iEV)='rad1'//a(iEV-9:iEV-9)//'.inp'
            OTFILE3(iEV)='rat1'//a(iEV-9:iEV-9)//'.inp'
         else if (iEV.ge.20.and.iEV.lt.30) then
            INFILE1(iEV)='sou2'//a(iEV-19:iEV-19)//'.inp'
            INFILE2(iEV)='tmp.inp'
            INFILE3(iEV)='f2'//a(iEV-19:iEV-19)//'.inp'
            OTFILE1(iEV)='d2'//a(iEV-19:iEV-19)//'.inp'
            OTFILE2(iEV)='rad2'//a(iEV-19:iEV-19)//'.inp'
            OTFILE3(iEV)='rat2'//a(iEV-19:iEV-19)//'.inp'
         endif
         call CountNumberOfChracters(INFILE1,iEV,1,iend1,MEV)
         call CountNumberOfChracters(INFILE2,iEV,1,iend6,MEV)
         call CountNumberOfChracters(INFILE3,iEV,1,iend2,MEV)
         call CountNumberOfChracters(OTFILE1,iEV,1,iend3,MEV)
         call CountNumberOfChracters(OTFILE2,iEV,1,iend4,MEV)
         call CountNumberOfChracters(OTFILE3,iEV,1,iend5,MEV)

         open(10,file=INFILE1(iEV)(1:iend1),status='old')
         open(11,file=INFILE2(iEV)(1:iend6),status='old')
         open(12,file=INFILE3(iEV)(1:iend2),status='old')
         open(13,file=OTFILE1(iEV)(1:iend3))
         open(14,file=OTFILE2(iEV)(1:iend4))
         open(15,file=OTFILE3(iEV)(1:iend5))


c    hypocentral distance
c
c
c	open(13, file="outdata")
c	open(14, file="outdata2")
c	open(15, file="outdata3")
c	open(16, file="outdata4")
c        OPEN(10, file="indata1",status='old')
c        OPEN(11, file="indata2",status='old')
c	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
      close(10)

      call UTM(alat,alon,alat1,alon1)

c   横浜市内150点の緯度、経度
c
      do i=1,MOB
         read(11,*)dlat(i),dlon(i)
         dlat0=dlat(i)
         dlon0=dlon(i)
         call UTM(dlat0,dlon0,UTMobservLat,UTMobservLon)
         dlat1(i)=UTMobservLat
         dlon1(i)=UTMobservLon
      enddo
      close(11)
C	write(6,1001)(dlat(i),dlon(i),i=1,MOB)

      read(12,*)(flag(j),j=1,MOB)
      close(12)



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

c 	do 11 i=1,MOB
c	As(i)=(cos(dlat(i)*PI/180))*(cos(dlon(i)*PI/180))
c	Bs(i)=(cos(dlat(i)*PI/180))*(sin(dlon(i)*PI/180))
c	Cs(i)=sin(dlat(i)*PI/180)
c11	continue

c	do 12 j=1,MOB
c12      d1(j)=(sqrt((Ae-As(j))**2+(Be-Bs(j))**2+(Ce-Cs(j))**2))*6371
c	do 20 k=1,MOB
c20	    dist(k)=sqrt(d1(k)**2+adepth**2)

      do 12 j=1,MOB
         d1(j)=dsqrt((alat1-dlat1(j))**2+(alon1-dlon1(j))**2)
12    continue
      do 20 k=1,MOB
         dist(k)=dsqrt((d1(k)**2)+adepth**2)
20    continue
	
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,MOB

      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)=dsqrt(fsv(j)**2+fsh(j)**2)
      rat(j)=(dabs(fsv(j)))/(dabs(fsh(j)))
 200	continue

      a1=0.d0
      a2=0.d0
      do 1 k=1,MOB
         a1=a1+rad(k)
         a2=a2+rad(k)**2
 1    continue
      s1=a1/MOB
      st1=dsqrt((MOB*a2-a1**2)/(MOB*149))
c
c
C    flag = 1 : recorded
c    flag = 0 : not recorded
       do 31 i=1,MOB
          if (flag(i).eq.1) go to 32
          if (flag(i).eq.0) go to 31
 32	         write(13,1000)dist(i)
 31    continue
      close(13)

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

      do 51 i=1,MOB
         if (flag(i).eq.1) go to 52
         if (flag(i).eq.0) go to 51
 52         write(15,1000)rat(i)
 51   continue
      close(15)

1000  format(1f10.5)

1001  format(1f10.5)


100   continue

      close(16)

      stop
      end





C**********************************************************************
                          subroutine UTM
     &           (TransLat,TransLon,UTMLat,UTMLon)
C**********************************************************************

      implicit real*8(a-h,o-z)
      real*8 TransLat,TransLon,UTMLat,UTMLon

C     地球楕円体の長半径
      a=6377.397155d0
C     地球楕円体の短半径
      b=6356.078963d0

      ee1=((a**2)-(b**2))/(a**2)
      ee2=((a**2)-(b**2))/(b**2)

      PI=dacos(-1.0d0)

      TransLatR=(TransLat/180.0d0)*PI
      TransLonR=(TransLon/180.0d0)*PI

C     UTMゾーンの原点
      if (TransLon.ge.120.0d0.and.TransLon.lt.126.0d0) then
         OriginLon=2.14675498d0
      else if (TransLon.ge.126.0d0.and.TransLon.lt.132.0d0) then
         OriginLon=2.25147474d0
      else if (TransLon.ge.132.0d0.and.TransLon.lt.138.0d0) then
         OriginLon=2.35619449d0
      else if (TransLon.ge.138.0d0.and.TransLon.lt.144.0d0) then
         OriginLon=2.46091425d0
      else if (TransLon.ge.144.0d0.and.TransLon.lt.150.0d0) then
         OriginLon=2.56563400d0
      else if (TransLon.ge.150.0d0.and.TransLon.lt.156.0d0) then
         OriginLon=2.67035376d0
      end if

      d=TransLonR-OriginLon

      TanLat=dtan(TransLatR)

      CosLat=dcos(TransLatR)

      SinLat=dsin(TransLatR)

      rr=ee2*(CosLat**2)

      CapitalA=a/(sqrt(1.0d0-(ee1*(SinLat**2))))

      B1=0.998329312962d0*TransLatR
      B2=0.00501415801d0*SinLat*CosLat
      B3=0.00001049328d0*((SinLat*(CosLat**3))-(CosLat*(SinLat**3)))
      CapitalB=a*(B1-B2+B3)

      g1=CapitalB
      g2=(CapitalA*(d**2)/2.0d0)*SinLat*CosLat
      g3=(d**2)*(CosLat**2)/12.0d0
      g4=5.0d0-(TanLat**2)+(9.0d0*rr)+(4.0d0*(rr**2))
      g5=(d**2)*(CosLat**2)/30.0d0
      g6=61.0d0-(58.0d0*(TanLat**2))
      g7=(TanLat**4)+(270.0d0*rr)-(330.0d0*(TanLat**2)*rr)
      g=g1+g2*(1.0d0+g3*(g4+g5*(g6+g7)))

      f1=CapitalA*d*CosLat
      f2=(d**2)*(CosLat**2)/6.0d0
      f3=1.0d0-(TanLat**2)+rr
      f4=(d**2)*(CosLat**2)/20.0d0
      f5=5.0d0-(18.0d0*(TanLat**2))+(TanLat**4)
      f6=(14.0d0*rr)-(58.0d0*(TanLat**2)*rr)
      f=f1*(1.0d0+f2*(f3+f4*(f5+f6)))

      UTMLon=(0.9996d0*f)+500.0d0
      UTMLat=0.9996d0*g

      return

      end



C***********************************************************************
      subroutine CountNumberOfChracters(a,n,istart,iend,MEV)
C***********************************************************************

      implicit real*8(a-h,o-z)

      character a(MEV)*100

      do 100 i=istart,100
         if (a(n)(i:i).eq.',') then
            go to 1
         else if (a(n)(i:i).eq.'') then
            go to 1
         end if
100   continue
1     iend=i-1

      return
      end





