program ENZFILES
      use nas_system
      parameter(np=16500)
      dimension f(np,6),t(np)
      character zeile*72,zeilanf*49,zeilend*36,fnam*18,ext(3)*1
      data ext /'e','n','z'/      
      
      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            

      open(17,status='old',file='refseis.sig')
      open(7,file='winplote.par')
      open(8,file='winplotn.par')
      open(9,file='winplotz.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)
        torad=atan(1.)/45.
        cosaz=cos(az*torad)
        sinaz=sin(az*torad)
        
c read six synthetic seismograms (far field r,t,v and near field r,t,v)
        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    Transform r,t to e,n
c => Stations-Azimut vom Epizentrum zur Station, von Nord ueber Ost (cw)
c => Radialkomponente weg vom Epizentrum
c => Transversalkomponente rechts (90 Grad cw) von der radialen
c => Vertikalkomponente nach oben (also alles wie ueblich)
c => let's hope that the moment-tensor components are understood the same way
c => order of components is e, n, z

       do k=1,3
         do i=1,nn
           f(i,k)=f(i,k)+f(i,k+3)
         enddo
       enddo  
       
       do i=1,nn
         fx    =f(i,1)*sinaz+f(i,2)*cosaz
         f(i,2)=f(i,1)*cosaz-f(i,2)*sinaz
         f(i,1)=fx
       enddo
 
c  loop over three components
        do k=1,3
c  remove trend. The array t (time) is now used for the signal.        
	  do i=1,nn
	    t(i)=f(i,k)
          enddo
          if(ntrend.gt.0) call trend(t,ntrend,nn)
          write(fnam,4) 20000+nrad,2000+naz,ext(k)
    4     format(7x,i5,i4,".",a1)
          fnam(1:8)='rdemux\r'
          fnam(13:13)='a'
          open(18,file=fnam)
          write(6,*) 'writing 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(k+6,'(a)') fnam
        enddo
      enddo
      write(7,'(a)') 'synthetic E seismograms from REFSEIS'
      write(8,'(a)') 'synthetic N seismograms from REFSEIS'
      write(9,'(a)') 'synthetic Z seismograms from REFSEIS'
      write(6,*)
      write(6,*) 
     & 'plot-parameter files for the E, N, Z components were written'
      write(6,*) 'copy one of the files winplote.par, winplotn.par,'
      write(6,*) 'winplotz.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