program testsig

c  common white noise and individual white noise
c
      parameter(ndim= 120000 )
c
      real*4 x(ndim),gauss(ndim)
      character name*10,par*20,text*35,noisfi*10,ja
      character*50 prolog(64),yes
      data par/'                    '/
      zpi=8.*atan(1.)
      text='     Test signal'
      write(6,*) 'Number of samples?'
      read(5,*) n
      if(n.gt.ndim) stop
      write(6,*) 'Sampling interval?'
      read(5,*) dt
      do 20 j=1,n
   20 gauss(j)=gran(j)
    6 write(6,*) 'Filename for the test signal?'
      read(5,'(a)') name
      npro=0
    1 continue
      write(6,*) 'Amplitude of common gaussian noise?'
      read(5,*) amp
      npro=npro+1
      write(prolog(npro),5) amp
    5 format('% common gaussian noise with amplitude ',f10.3)
      do 3 j=1,n
    3 x(j)=amp*gauss(j)
    2 write(6,*) 'Add Gaussian noise with rms Amplitude:'
      read(5,*) amp
      do 18 j=1,n
   18 x(j)=x(j)+amp*gran(j)
      npro=npro+1
      write(prolog(npro),9) amp
    9 format('%  white noise, amp=',f10.3)
      write(6,*) 'Apply sine taper? (y/n):'
      read(5,'(a)') yes
      if(yes.eq.'y') then
        f=zpi/(2.*n)
        do 19 j=1,n
   19     x(j)=x(j)*sin(j*f)
        npro=npro+1
        write(prolog(npro),8) 
    8   format('%  sine taper applied')
      endif

    7 write(6,'(a)') 'Fortran format, default <cr>= (8f9.3) :'
      read(5,'(a20)') par                           
      call output(name,prolog,npro,par,text,n,dt,x)
      write(6,'(a)') 'Done. Another signal? (y/n)'
      read(5,'(a)') ja
      if(ja.eq.'y') goto 6
      stop
      end

      subroutine output(name,prolog,npro,par,text,n,dt,x)
      real*4 x(n)
      character name*10,text*35,par*20,iform*20
      character*50 prolog(64)
      iform='(5f12.3)            '
      if(par(1:1).eq.'(') read(par,'(a20)') iform
      write(6,5) name,iform
    5 format('Signal is written into file ',a10/' format: ',a20)
      open(8,file=name)
      write(8,'(a)') text
      do j=1,npro
      write(8,'(a)') prolog(j)
      enddo
      write(8,1) n,iform,dt
    1 format(i10,a20,3f10.3,2i5,a10)
      write(8,iform) (x(j),j=1,n)
      close(8)
      return
      end

      function gran(j)
c  return a gaussian random number
      dimension x(12)
      call random_number(x)
      gr=x(1)
      do k=2,12
        gr=gr+x(k)
      enddo
      gran=gr-6.
      return
      end