program echosim

! Generate simulated echo-mode files -- self-echo or "measure" 

  use wavhdr
  parameter (NWAVE=27648,NMAX=32768,NZ=36000)
  type(hdr) h                            !Header for .wav file
  character arg*12,fname*17
  complex c0(0:NMAX-1)
  complex c(0:NMAX-1)
!  complex cwave(0:NWAVE-1)
  real*8 f0,dt,twopi,phi,dphi
  real wave(NZ)
  integer*2 iwave(NZ)                  !Generated full-length waveform
  equivalence (nDop0,iwave(1))
  equivalence (nDopAudio0,iwave(3))
  equivalence (nfrit0,iwave(5))
  equivalence (f10,iwave(7))
  equivalence (fspread0,iwave(9))

! Get command-line argument(s)
  nargs=iargc()
  if(nargs.ne.5) then
     print*,'Usage:    echosim   f0   fdop fspread nfiles snr'
     print*,'Examples: echosim  1500   0.0   4.0     10   -22'
     go to 999
  endif

  call getarg(1,arg)
  read(arg,*) f0                         !Tone frequency
  call getarg(2,arg)
  read(arg,*) fdop                       !Doppler shift (Hz)
  call getarg(3,arg)
  read(arg,*) fspread             !Frequency spread (Hz) (JHT Lorentzian model)
  call getarg(4,arg)
  read(arg,*) nfiles                     !Number of files
  call getarg(5,arg)
  read(arg,*) snrdb                      !SNR_2500

  twopi=8.d0*atan(1.d0)
  fs=12000.0                             !Sample rate (Hz)
  dt=1.d0/fs                              !Sample interval (s)
  bandwidth_ratio=2500.0/(fs/2.0)
  sig=sqrt(2*bandwidth_ratio) * 10.0**(0.05*snrdb)
  if(snrdb.gt.90.0) sig=1.0
  dphi=twopi*(f0+fdop)*dt

  write(*,1000)
1000 format('   N   f0     fDop fSpread   SNR  File name'/51('-'))

  do ifile=1,nfiles
     phi=0.d0
     do i=0,NWAVE-1
        phi=phi + dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        c0(i)=cmplx(cos(xphi),sin(xphi))
     enddo
     c0(NWAVE:)=0.
     if(fspread.gt.0.0) call fspread_lorentz(c0,fspread)
     c=sig*c0
     wave(1:NWAVE)=imag(c(1:NWAVE))
     peak=maxval(abs(wave))
     if(snrdb.lt.90) then
        do i=1,NWAVE                   !Add gaussian noise at specified SNR
           xnoise=gran()
           wave(i)=wave(i) + xnoise
        enddo
        do i=NWAVE+1,NZ
           xnoise=gran()
           wave(i)=xnoise
        enddo
     endif

     gain=100.0
     if(snrdb.lt.90.0) then
       wave=gain*wave
     else
       datpk=maxval(abs(wave))
       fac=32766.9/datpk
       wave=fac*wave
     endif
     if(any(abs(wave).gt.32767.0)) print*,"Warning - data will be clipped."
     iwave=nint(wave)

     nDop0=nint(fdop)
     nDopAudio0=0
     nfrit0=0
     f10=f0 + fdop
     fspread0=fspread
     
     h=default_header(12000,NMAX)
     write(fname,1102) 3*ifile
1102 format('000000_',i6.6,'.wav')
     open(10,file=fname,status='unknown',access='stream')
     write(10) h,iwave                !Save to *.wav file
     close(10)
     write(*,1110) ifile,f0,fdop,fspread,snrdb,fname
1110 format(i4,4f7.1,2x,a17)
  enddo

999 end program echosim