program CATFILES
c this version concatenates the three output files for each station
      use nas_system
      parameter(np=16500)
      parameter(ncat=3*np)
      dimension f(np,6),t(np),cat(ncat)
      character zeile*72,zeilanf*49,zeilend*36,fnam*20    
      
      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='winplotenz.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) 
      nn=n
      if(nmax.gt.0) nn=min(n,nmax) 
      write(7,8) nstat,24.,15.,3*nn
    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        
        
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  concatenate 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)
          do i=1,nn
            cat((k-1)*nn+i)=t(i)
          enddo
        enddo

c store concatenated components        
        write(fnam,4) 20000+nrad,2000+naz,'.enz'
  4     format(7x,i5,i4,a4)
        fnam(1:8)='rdemux\r'
        fnam(13:13)='a'
        open(18,file=fnam)
        write(6,*) 'writing file ',fnam
        write(18,*) 'concatenated e,n,z seismograms from REFSEIS'
        write(18,5) 3*nn,dt
  5     format(i10,"(5e12.4)",12x,f10.6)
        write(18,'(5e12.4)') (cat(i),i=1,3*nn)
        close(18)
        write(7,'(a)') fnam                
      enddo
      write(7,'(a)') 'concatenated e,n,z seismograms from REFSEIS'
      write(6,*)
      write(6,*) 'plot-parameter file for the concatenated E, N, Z compo
     &nents was written'
      write(6,*) 
     & 'copy file winplotenz.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