program rfiles
      use nas_system
      parameter(np=16500)
      dimension f(np,6),t(np)
      character zeile*72,zeilanf*49,zeilend*36,fnam*19,ext(9)*2
      
      write(*,*)
      write(*,*) '>>> Reformatting synthetic seismograms from REFSEIS'
      if(iargc().eq.0) then
      write(*,*) 'To remove a baseline trend before the first arrival,'
      write(*,*) 'enter number of samples before arrival as a runstring'
      write(*,*) 'parameter. A second runstring parameter, if present,'
      write(*,*) 'limits the number of samples.'
      write(*,*)
        ntrend=0
        nmax=0
      else if(iargc().eq.1) then
        call getarg(1,zeile)
        read(zeile,*) ntrend    
      else if(iargc().eq.2) then
        call getarg(1,zeile)
        read(zeile,*) ntrend    
        call getarg(2,zeile)
        read(zeile,*) nmax    
      endif            

c  code: f=farfield, q=nearfield, r=radial, t=transverse, v=vertical
      data ext /'fr','ft','fv','qr','qt','qv','tr','tt','tv'/
      open(17,status='old',file='refseis.sig')
      open(7,file='winplotr.par')
      open(8,file='winplott.par')
      open(9,file='winplotv.par')
      read(17,'(a)') zeile
      write(6,'(a)') zeile
      read(17,'(a,a)') zeilanf,zeilend
      write(6,'(a)') zeilanf
      write(6,'(a)') zeilend
      write(6,*)
      read(zeilanf,1) n,nstat,t0,tnm1
    1 format(4x,i5,4x,i3,6x,f10.4,7x,f10.4)
      dt=(tnm1-t0)/(n-1)
      write(7,8) nstat,24.,15.,n
      write(8,8) nstat,24.,15.,n
      write(9,8) nstat,24.,15.,n
    8 format("0,",i3,",",f5.1,",",f5.1,",","  1,",i8,", 0.8")
    
c  loop over the stations
      do m=1,nstat
        read(17,'(a)') zeile
        write(6,'(a)') zeile
        read(zeile,2) rad,az
    2   format(4x,f10.4,5x,f6.2)
        nrad=nint(10.*rad)
        naz=nint(az)
        
c loop over time
        do i=1,n
          read(17,3) t(i),(f(i,k),k=1,6)
    3     format(1x,f8.4,6(1x,e10.4))
        enddo
        
        nn=n
        if(nmax.gt.0) nn=min(n,nmax)
 
c  loop over three components
        do k=7,9
c  remove trend. The array t (time) is no longer used an serves as 
c  a buffer for the total field.        
	  do i=1,nn
	    t(i)=f(i,k-6)+f(i,k-3)
          enddo
          if(ntrend.gt.0) call trend(t,ntrend,nn)
          write(fnam,4) 20000+nrad,2000+naz,ext(k)
    4     format(7x,i5,i4,".",a2)
          fnam(1:8)='rdemux\r'
          fnam(13:13)='a'
          open(18,file=fnam)
          write(18,*) 'synthetic seismogram from refseis ',fnam
          write(18,5) nn,dt
    5     format(i10,"(5e10.4)",12x,f10.6)
          write(18,'(5e10.4)') (t(i),i=1,nn)
          close(18)
          write(6,*) 
     &'total-field data were demultiplexed and written to file ',fnam
          write(k,'(a)') fnam
        enddo
      enddo
      write(7,'(a)') 'synthetic radial seismograms from REFSEIS'
      write(8,'(a)') 'synthetic transverse seismograms from REFSEIS'
      write(9,'(a)') 'synthetic vertical seismograms from REFSEIS'
      write(6,*)
      write(6,*) 
     & 'plot-parameter files for the r, t, v components were written'
      write(6,*) 'copy one of the files winplotr.par, winplott.par,'
      write(6,*) 'winplotv.par to winplot.par, then type winplot'
      end
	
      subroutine trend(x,ntrend,n)
c  determine trend from samples 1..ntrend, remove from all samples
      parameter(np=16500)
      dimension x(np)
      gn = ntrend
      alpha = 0.5*gn*(gn+1.)
      beta = (2.*gn+1.)*(gn+1.)*gn/6.
      det = gn*beta-alpha*alpha
      sx = 0.
      sjx = 0.
      do 1001 j=1,ntrend
        sx = sx+x(j)
 1001   sjx = sjx+j*x(j)
      a = (sx*beta-sjx*alpha)/det
      b = (sjx*gn-sx*alpha)/det
      do 1002 j=1,n
 1002   x(j) = x(j)-a-b*j
      return
      end