program testsig
c
c   Generats a random-telegraph signal

      parameter(ndim= 120000 )
c
      real*4 x(ndim)
      character name*10,par*20,text*35,noisfi*10,ja
      character*50 prolog(2),yes
      data par/'                    '/
      text='     random-telegraph signal'
      zpi=8.*atan(1.)
      write(6,*) 'Number of samples?'
      read(5,*) n
      if(n.gt.ndim) stop
      write(6,*) 'Sampling interval?'
      read(5,*) dt
      write(6,*)'clock interval as a multiple of the sampling interval:'
      read(5,*) nper
      if(nper.lt.1) then
        write(*,*) 'impossible clock interval'
        stop
      endif
      write(6,*) 'Amplitude?'
      read(5,*) amp
    6 write(6,*) 'Filename for the test signal?'
      read(5,'(a)') name
      npro=1
      write(prolog(npro),5) per*dt,amp
    5 format('%   clock interval=',f10.3,'sec, amp=',f10.3)
      
      do i=1,n
        x(i)=amp
        if(mod(i,nper).eq.1) then
          call random_number(z)
          if(z.gt.0.5) then
            amp=-amp
          endif
        endif
      enddo
 
      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=2
        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)

      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(2)
      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 i=1,npro
      write(8,'(a)') prolog(i)
      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