From 87dcde7564bc34899344509bab4272ee435227ea Mon Sep 17 00:00:00 2001
From: Joe Taylor <joe@princeton.edu>
Date: Fri, 22 Apr 2022 11:43:28 -0400
Subject: [PATCH] Move calculation of echo snr into a separate subroutine.

---
 lib/test_snr.f90 | 89 ++++++++++++++++++++++++++++++------------------
 1 file changed, 56 insertions(+), 33 deletions(-)

diff --git a/lib/test_snr.f90 b/lib/test_snr.f90
index 04870a1a6..d1ea1685f 100644
--- a/lib/test_snr.f90
+++ b/lib/test_snr.f90
@@ -12,9 +12,12 @@ program test_snr
   complex c(0:NH)
   real x(NFFT)
   real s(8192)
+  real sa(4096)
+  real sb(4096)
+  real red(4096)
+  real blue(4096)
   character*80 infile
-  integer ipkv(1)
-  equivalence (x,c),(ipk,ipkv)
+  equivalence (x,c)
 
   nargs=iargc()
   if(nargs.lt.1) then
@@ -45,36 +48,11 @@ program test_snr
      enddo
   enddo
 
-  wh=0.5*fspread0+10.0
-  i1=nint((1500.0 - 2.0*wh)/df)
-  i2=nint((1500.0 - wh)/df)
-  i3=nint((1500.0 + wh)/df)
-  i4=nint((1500.0 + 2.0*wh)/df)
-
-  baseline=(sum(s(i1:i2-1)) + sum(s(i3+1:i4)))/(i2+i4-i1-i3)
-  s=s/baseline
-  psig=sum(s(i2:i3)-1.0)
-  pnoise_2500 = 2500.0/df
-  snrdb=db(psig/pnoise_2500)
-
-  smax=0.
-  mh=max(1,nint(0.2*fspread0/df))
-  do i=i2,i3
-     ssum=sum(s(i-mh:i+mh))
-     if(ssum.gt.smax) then
-        smax=ssum
-        ipk=i
-     endif
-  enddo
-  fpeak=ipk*df - 1500.0
-
-  call averms(s(i1:i2-1),i2-i1,-1,ave1,rms1)
-  call averms(s(i3+1:i4),i4-i3,-1,ave2,rms2)
-  perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3))
-  snr_detect=psig/perr
-  db_err=99.0
-  if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500)
+  sa=s(2049:6144)
+  sb=s(2049:6144)
 
+  call echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
+  
   nsum=npings
   nqual=0
   if(nsum.ge.2 .and. nsum.lt.4)  nqual=(snr_detect-4)/5
@@ -84,10 +62,10 @@ program test_snr
   if(nqual.lt.0)  nqual=0
   if(nqual.gt.10) nqual=10
 
-  write(*,1010) fspread0,wh,snrdb0,snrdb,snrdb-snrdb0,db_err,fpeak,   &
+  write(*,1010) fspread0,snrdb0,snrdb,snrdb-snrdb0,db_err,fpeak,   &
        snr_detect,nqual
 
-1010 format(6f6.1,2f7.1,i3)
+1010 format(5f6.1,2f7.1,i4)
 
   do i=1,8192
      write(12,1100) i*df,s(i)
@@ -95,3 +73,48 @@ program test_snr
   enddo
 
 999 end program test_snr
+
+subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
+
+  parameter (NZ=4096)
+  real sa(NZ)
+  real sb(NZ)
+  real blue(NZ)
+  real red(NZ)
+  integer ipkv(1)
+  equivalence (ipk,ipkv)
+
+  df=12000.0/32768.0
+  wh=0.5*fspread0+10.0
+  i1=nint((1500.0 - 2.0*wh)/df) - 2048
+  i2=nint((1500.0 - wh)/df) - 2048
+  i3=nint((1500.0 + wh)/df) - 2048
+  i4=nint((1500.0 + 2.0*wh)/df) - 2048
+
+  baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3)
+  blue=sa/baseline
+  red=sb/baseline
+  psig=sum(red(i2:i3)-1.0)
+  pnoise_2500 = 2500.0/df
+  snrdb=db(psig/pnoise_2500)
+
+  smax=0.
+  mh=max(1,nint(0.2*fspread0/df))
+  do i=i2,i3
+     ssum=sum(red(i-mh:i+mh))
+     if(ssum.gt.smax) then
+        smax=ssum
+        ipk=i
+     endif
+  enddo
+  fpeak=ipk*df - 750.0
+
+  call averms(red(i1:i2-1),i2-i1,-1,ave1,rms1)
+  call averms(red(i3+1:i4),i4-i3,-1,ave2,rms2)
+  perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3))
+  snr_detect=psig/perr
+  db_err=99.0
+  if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500)
+
+  return
+end subroutine echo_snr