Modify echosim to use the same Lorentzian model for fspread as used in q65sim.

CMakeLists.txt

@@ -420,6 +420,7 @@ set (wsjt_FSRCS
   lib/fmtmsg.f90
   lib/foldspec9f.f90
   lib/four2a.f90
+  lib/fspread_lorentz.f90
   lib/ft8/foxfilt.f90
   lib/ft8/foxgen.f90
   lib/ft8/foxgen_wrap.f90

lib/avecho.f90

@@ -87,8 +87,12 @@ subroutine avecho(id2,ndop,nfrit,nauto,nqual,f1,xlevel,snrdb,db_err,  &
      call smo121(blue,NZ)
   enddo
 
-!  write(*,3001) snrdb,db_err,dfreq,snr_detect,redmax,nqual,nsmo,nclearave,nsum
-!3001 format('A',5f10.1,4i4)
+  ia=200.0/df
+  ib=400.0/df
+  call pctile(red(ia:ib),ib-ia+1,50,bred)
+  red=red-bred
+  call pctile(blue(ia:ib),ib-ia+1,50,bblue)
+  blue=blue-bblue
 
 900  return
 end subroutine avecho

lib/echosim.f90

@@ -20,9 +20,9 @@ program echosim
 
 ! Get command-line argument(s)
   nargs=iargc()
-  if(nargs.ne.6) then
-     print*,'Usage:    echosim   f0   fdop fspread delay nfiles snr'
-     print*,'Examples: echosim  1500   0.0   4.0    5.0    10   -22'
+  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
 
@@ -31,12 +31,10 @@ program echosim
   call getarg(2,arg)
   read(arg,*) fdop                       !Doppler shift (Hz)
   call getarg(3,arg)
-  read(arg,*) fspread                    !Watterson frequency spread (Hz)
+  read(arg,*) fspread                    !Frequency spread (Hz) (JHT Lorentzian model)
   call getarg(4,arg)
-  read(arg,*) delay                      !Watterson delay (ms)
-  call getarg(5,arg)
   read(arg,*) nfiles                     !Number of files
-  call getarg(6,arg)
+  call getarg(5,arg)
   read(arg,*) snrdb                      !SNR_2500
 
   twopi=8.d0*atan(1.d0)
@@ -56,7 +54,7 @@ program echosim
         c0(i)=cmplx(cos(xphi),sin(xphi))
      enddo
      c0(NWAVE:)=0.
-     if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c0,NMAX,NWAVE,fs,delay,fspread)
+     if(fspread.gt.0.0) call fspread_lorentz(c0,fspread)
      c=sig*c0
      wave(1:NWAVE)=imag(c(1:NWAVE))
      peak=maxval(abs(wave))

lib/fspread_lorentz.f90

@@ -0,0 +1,47 @@
+subroutine fspread_lorentz(cdat,fspread)
+
+  parameter (NZ=3*12000)
+  complex cdat(0:NZ-1)
+  complex cspread(0:NZ-1)
+  complex z
+
+  twopi=8.0*atan(1.0)
+  nfft=NZ
+  nh=nfft/2
+  df=12000.0/nfft
+  cspread(0)=1.0
+  cspread(nh)=0.
+  b=6.0                       !Use truncated Lorenzian shape for fspread
+  do i=1,nh
+     f=i*df
+     x=b*f/fspread
+     z=0.
+     a=0.
+     if(x.lt.3.0) then                          !Cutoff beyond x=3
+        a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian amplitude
+        phi1=twopi*rran()                       !Random phase
+        z=a*cmplx(cos(phi1),sin(phi1))
+     endif
+     cspread(i)=z
+     z=0.
+     if(x.lt.3.0) then                !Same thing for negative freqs
+        phi2=twopi*rran()
+        z=a*cmplx(cos(phi2),sin(phi2))
+     endif
+     cspread(nfft-i)=z
+  enddo
+
+  call four2a(cspread,nfft,1,1,1)             !Transform to time domain
+  
+  sum=0.
+  do i=0,nfft-1
+     p=real(cspread(i))**2 + aimag(cspread(i))**2
+     sum=sum+p
+  enddo
+  avep=sum/nfft
+  fac=sqrt(1.0/avep)
+  cspread=fac*cspread                   !Normalize to constant avg power
+  cdat=cspread*cdat                     !Apply Rayleigh fading
+
+  return
+end subroutine fspread_lorentz