program Et2SacMH
c Generate minimum-header SAC data files from PSE data on event tapes
c ** 18-Jul-00 **
c f77 PseSacGenMH.f Ntimex.f lpedit.f spedit.f wsac1.f
c limits: 720 records (10h:52m) for LP and 90 records (1h:21m) for SPZ
c
	real data(259200)
	integer in(4864),lp(360,3),sp(2882),it(4)
	integer*2 exhdr(8)
	logical newfmt
	character cmp*3(4),ext*8,fname*16
	equivalence (exhdr(1),in(1))
	data cmp/'lpx','lpy','lpz','spz'/
c
c open input PSE ET data file in Exabyte format
10	write(*,"('Open etdata.'$)")
	read(*,'(a)') ext
	if (ext(1:1).eq.' ') stop
	fname='etdata.'//ext
	open(1,file=fname,access='direct',form='unformatted',recl=19456)
c name output file
	write(*,"('Output file psac.'$)")
	read(*,'(a)') ext
	i=2
	do while (ext(i:i).ne.' ')
	 i=i+1
	 end do
	 ii=i-1
c specify data length and set sample interval
20	write(*,"('Record/block range: '$)")
	read(*,*) n1,nb1,n2,nb2
	read(1,rec=n1) in
	nlb=(exhdr(6)+1)*3
	nch=4-exhdr(6)
	npt=((n2-n1)*nlb+nb2-nb1+1)*360
	if (npt.gt.259200) then
	 write(*,"('Data too long')")
	 go to 20
	 end if
	lbs=4860/nlb
	newfmt=exhdr(6).eq.1
	ist=exhdr(2)
	it(1)=exhdr(5)
	call time03(in(5),it(2))
	if (nb1.gt.1) then
	 ih=it(2)*24+it(3)/100
	 call time03(in((nb1-1)*lbs+5),it(2))
	 if (it(2)*24+it(3)/100.lt.ih) it(1)=it(1)+1
	 end if
	sec=it(4)/1000.
	write(*,"('Data start at',3i5,f7.3)") (it(i),i=1,3),sec
	delt=.1509375
c process data
	do kp=1,nch
	 if (kp.eq.4) then
	  npt=npt*8
	  if (npt.gt.259200) then
	   npt=259200
	   n2=n1+30
	   nb2=nb1-1
	   if (nb2.lt.1) then
	    nb2=nb2+3
	    n2=n2-1
	    end if
	   end if
	  delt=.01886719
	  end if
	 j=1
	 do nr=n1,n2
	  read(1,rec=nr) in
	  l1=1
	  l2=nlb
	  if (nr.eq.n1) l1=nb1
	  if (nr.eq.n2) l2=nb2
	  do lr=l1,l2
	   ip=(lr-1)*lbs+5
	   if (kp.lt.4) then
	    call lpedit(in(ip),lp,newfmt)
	    do i=1,360
	     data(j)=lp(i,kp)
	     j=j+1
	     end do
	   else
	    call spedit(in(ip),sp,ist)
	    if (nr.eq.n1.and.lr.eq.1) sp(3)=sp(4)
	    do i=3,2882
	     data(j)=sp(i)
	     j=j+1
	     end do
	    end if
	   end do
	  end do
c write a SAC data file
	 fname='psac.'//ext(1:ii)//'.'//cmp(kp)
	 call wsac1(fname,data,npt,0.,delt,nerr)
	 write(*,*) npt,' written to '//fname
	 end do
	close(1)
	go to 10
	end