      Program Plain2
C     IF You input Strike, Dip, Epicenter and fault plain
C     , to output Center of Subfaulte fault Plain at 
C     Longitude and Latitude.
C     You must run program [Plain] first,
C     and second run this program [Plain2].
C
      parameter (nm=1500)
      real x(nm),y(nm),xd(nm),yd(nm),ax(nm),ay(nm)
     &     ,str,dip
      pideg=3.1415926/180.
cc================================================
      OPEN(11,file='Line.dat')
      OPEN(17,file='mline.dat')
      OPEN(18,file='line.dat')
      read(5,*)olat,olon
      read(5,*)str,dip
      read(11,*)nn
      read(11,*)(xd(i),yd(i),i=1,nn)
      write(18,'(a)')'>'
      do i=1,nn,2
         write(18,200)yd(i),-xd(i)
         write(18,200)yd(i+1),-xd(i+1)
         write(18,'(a)')'>'
      end do
c***** degree/km **********************************
c
      call geocen(olat,olon,olt,oln,pideg)
      hlat=olat-0.1
      hlon=olon
      call geocen(hlat,hlon,hlt,hln,pideg)
      call distan(olt,oln,hlt,hln,delta)
      deglat=0.1/delta
      hlat=olat
      hlon=olon-0.1
      call geocen(hlat,hlon,hlt,hln,pideg)
      call distan(olt,oln,hlt,hln,delta)
      deglon=0.1/delta
      write(6,*)'deglat=',deglat
      write(6,*)'deglon=',deglon
cc==============================================
      do i=1,nn
         xd(i)=xd(i)*cos(dip*pideg)
      end do
      str=str*pideg
      do i=1,nn
         x(i)=xd(i)*cos(str)+yd(i)*sin(str)
         y(i)=-xd(i)*sin(str)+yd(i)*cos(str)
      end do
      do i=1,nn
         ax(i)=olon+deglon*x(i)
         ay(i)=olat+deglat*y(i)
      end do
      write(17,'(a)')'>'
      do i=1,nn,2
         write(17,200)ay(i),ax(i)
         write(17,200)ay(i+1),ax(i+1)
         write(17,'(a)')'>'
      end do
200    format(2f12.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


