For fst4sim, use Lorentzian fading spectrum when fspread is negative.

lib/fst4/fst4sim.f90

@@ -118,7 +118,8 @@ program fst4sim
 
    do ifile=1,nfiles
       c=c0
-      if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,nwave,NZ,fs,delay,fspread)
+      if(fspread.gt.0.0 .or. delay.ne.0.0) call watterson(c,nwave,NZ,fs,delay,fspread)
+      if(fspread.lt.0.0) call lorentzian_fading(c,nwave,fs,-fspread)
       c=sig*c
       wave=real(c)
       if(snrdb.lt.90) then

lib/fst4/lorentzian_fading.f90

@@ -0,0 +1,43 @@
+subroutine lorentzian_fading(c,npts,fs,fspread)
+!
+! npts is the total length of the simulated data vector
+!
+  complex c(0:npts-1)
+  complex cspread(0:npts-1)
+  complex z
+
+  twopi=8.0*atan(1.0)
+  df=fs/npts
+  nh=npts/2
+  cspread(0)=1.0
+  cspread(nh)=0.
+  b=6.0
+  do i=1,nh
+     f=i*df
+     x=b*f/fspread
+     z=0.
+     a=0.
+     if(x.lt.3.0) then
+        a=sqrt(1.111/(1.0+x*x)-0.1)
+        phi1=twopi*rran()
+        z=a*cmplx(cos(phi1),sin(phi1))
+     endif
+     cspread(i)=z
+     z=0.
+     if(x.lt.3.0) then
+        phi2=twopi*rran()
+        z=a*cmplx(cos(phi2),sin(phi2))
+     endif
+     cspread(npts-i)=z
+  enddo
+
+  call four2a(cspread,npts,1,1,1)
+
+  s=sum(abs(cspread)**2)
+  avep=s/npts
+  fac=sqrt(1.0/avep)
+  cspread=fac*cspread
+  c=cspread*c
+   
+  return
+end subroutine lorentzian_fading