      Program splot
C     IF You input Strike, Dip, Epicenter and center of 
C     Subfaulte location, output location for each Subfault
C     (Longitude and Latitude).
C     Parameter
C     Xd : strike Direction (origin = epicenter)  [km]
C     Yd : Dip Direction ( Up + ) (origin = epicenter) [km]
C     Ay : latitude
C     Ax : longitude
C     alat : latitude for Epicenter
C     alon : longitude for Epiceter
C     Str, dip   : Strike , Dip for mainshock
C     slip : direction of slip [Degree]
C     dm   : Displacement of slip [m]
C     Ed   : Error of slip [m]
C     NN   : total number of data 
C     vel  : unit length for slip vectors
C     Input File
C     SD.dat :  location ( x,y ), slip, ed
C     Output File
C     MD.dat : location and dm
C     MED.dat : location and Ed
C     mask.xyz : location and dm
C     [Only the range of the estimated error of 80% or less.]
C     mv.dat  : location and slip vecter 
C
      parameter (nm=1000,vel=0.2)
      real xd(nm),yd(nm),ax(nm),ay(nm)
     &     ,str,dip,dm(nm),slip(nm),ed(nm)
C
C	Open File
C
      OPEN(10,file='SD.dat')
      OPEN(17,file='MD.dat')
      OPEN(18,file='mv.dat')
      OPEN(19,file='MED.dat')
      OPEN(20,file='mask.xyz')
C
C	Read SD.dat
C
      read(10,*)alat,alon,nn
      read(10,*)str,dip
      read(10,*)(xd(i),yd(i),dm(i),slip(i),ed(i),i=1,nn)
C
C       CALCULATION
C
C       Get Degree/km Facter (deglat,deglon)
C
      pideg=3.1415926/180
      call geocen(alat,alon,olt,oln,pideg)
      hlat=alat-0.1
      hlon=alon
      call geocen(hlat,hlon,hlt,hln,pideg)
      call distan(olt,oln,hlt,hln,delta)
      deglat=0.1/delta
      hlat=alat
      hlon=alon-0.1
      call geocen(hlat,hlon,hlt,hln,pideg)
      call distan(olt,oln,hlt,hln,delta)
      deglon=0.1/delta
C
C       Get Plot Location ( Degree ; lat,lon)
C
      dip=dip*pideg
      str1=(str-90.)*pideg
      do i=1,nn
         slip(i) = 90. - str + slip(i)
         yd(i) =  yd(i)*cos(dip)
         work1 =  xd(i)*cos(str1) + yd(i)*sin(str1)
         work2 = -xd(i)*sin(str1) + yd(i)*cos(str1)
         ax(i) = alon + deglon * work1
         ay(i) = alat + deglat * work2
      end do
C
C	Out Put File
C	
      write(17,200)(ay(i),ax(i),dm(i),i=1,nn)
      write(18,210)(ay(i),ax(i),slip(i),dm(i)*vel,i=1,nn)
      do i=1,nn
         if(dm(i).gt.0.)then
            write(19,200)ay(i),ax(i),ed(i)
         end if
         if (dm(i)*0.8.gt.ed(i))then
            write(20,200)ay(i),ax(i),dm(i)
         end if
      end do
      
 200  format(2f12.6,3x,f10.6)
 210  format(2f12.6,3x,2f15.6)
cc========================By yagi(97.6.29)======
      stop
      end
c*************geocen*************************
c
      subroutine geocen(atdeg,ondeg,atrad,onrad,pideg)
c
c     GEOGRAPHICAL--->GEOCENTRICAL COORDINATE
c
      obla=atdeg*pideg
      atrad=obla-0.00337267*sin(2*obla)
      onrad=ondeg*pideg
      return
      end
c
c**************distan*********************
c
      subroutine distan(soulat,soulon,obslat,obslon,r)
c
c     CALCULATING EPICENTRAL DISTANCE
c
      snsola=sin(soulat)
      cssola=cos(soulat)
      ramda=obslon-soulon
      phy=obslat-soulat
      snphy=sin(phy/2.)
      snphy=snphy*snphy
      snram=sin(ramda/2.)
      snram=snram*snram
      csobla=cos(obslat)
      ws=snphy+snram*cssola*csobla
      wc=abs(1-ws)
      wc=sqrt(ws/wc)
      r=12739*atan(wc)
      return
      end

