From 701d517e6e0fb1b8c84c7a0f0598bd03f54aa7cb Mon Sep 17 00:00:00 2001
From: Steven Franke <s.j.franke@icloud.com>
Date: Tue, 22 Dec 2020 12:33:24 -0600
Subject: [PATCH] For fst4sim, use Lorentzian fading spectrum when fspread is
 negative.

---
 lib/fst4/fst4sim.f90           |  3 ++-
 lib/fst4/lorentzian_fading.f90 | 43 ++++++++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+), 1 deletion(-)
 create mode 100644 lib/fst4/lorentzian_fading.f90

diff --git a/lib/fst4/fst4sim.f90 b/lib/fst4/fst4sim.f90
index 5b8f33bc3..e42bbbb02 100644
--- a/lib/fst4/fst4sim.f90
+++ b/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
diff --git a/lib/fst4/lorentzian_fading.f90 b/lib/fst4/lorentzian_fading.f90
new file mode 100644
index 000000000..c1bd5c325
--- /dev/null
+++ b/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