From 5693be1ee149240b17b43d9f8a5b3d7f7b49058b Mon Sep 17 00:00:00 2001
From: Joe Taylor <joe@princeton.edu>
Date: Thu, 7 May 2020 13:54:02 -0400
Subject: [PATCH] Fix bugs for generating and decoding JT65 shorthand messages.

---
 lib/jt65sim.f90 |  4 +++-
 lib/sh65.f90    | 51 +++++++++++++++++++++++++++----------------------
 2 files changed, 31 insertions(+), 24 deletions(-)

diff --git a/lib/jt65sim.f90 b/lib/jt65sim.f90
index 26e558f44..623a0fca9 100644
--- a/lib/jt65sim.f90
+++ b/lib/jt65sim.f90
@@ -196,9 +196,11 @@ program jt65sim
           if(msg(1:2).eq.'73') nshorthand=4
           if(nshorthand.gt.0) then
             ntoggle=0
-            do i=1,nsym,2
+            do i=1,nsym,4
               itone(i)=ntoggle*10*nshorthand
               if(i+1.le.126) itone(i+1)=ntoggle*10*nshorthand
+              if(i+2.le.126) itone(i+2)=ntoggle*10*nshorthand
+              if(i+3.le.126) itone(i+3)=ntoggle*10*nshorthand
               ntoggle=mod(ntoggle+1,2)
             enddo
           endif
diff --git a/lib/sh65.f90 b/lib/sh65.f90
index 8f13fdc3a..09bef9c24 100644
--- a/lib/sh65.f90
+++ b/lib/sh65.f90
@@ -1,20 +1,19 @@
 subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
-  parameter(NFFT=2048,NH=NFFT/2,MAXSTEPS=150)
-  complex cx(90000)
+  parameter(NFFT=2048,NH=NFFT/2)
+  complex cx(n5)               !Centered on nfqso, sample rate 1378.125
   complex c(0:NFFT-1)
   real s(-NH+1:NH)
-  real s2(-NH+1:NH,MAXSTEPS)
-  real ss(-NH+1:NH,8)
-  real sigmax(8)
-  integer ipk(8)
+  real ss(-NH+1:NH,16)
+  real sigmax(16)
+  integer ipk(16)
 
-  s=0.
   ss=0.
 
-  jstep=NFFT/4
-  nblks=n5/jstep - 3
+  jstep=NFFT/8
+  nblks=272
   ia=-jstep+1
   do iblk=1,nblks
+     n=mod(iblk-1,16) + 1
      ia=ia+jstep
      ib=ia+NFFT-1
      c=cx(ia:ib)
@@ -22,34 +21,37 @@ subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
      do i=0,NFFT-1
         j=i
         if(j.gt.NH) j=j-NFFT
-        p=real(c(i))**2 + aimag(c(i))**2
-        s(j)=s(j) + p
-        s2(j,iblk)=p
+        ss(j,n)=ss(j,n) + real(c(i))**2 + aimag(c(i))**2
      enddo
-     n=mod(iblk-1,8) +1
-     ss(-NH+1:NH,n)=ss(-NH+1:NH,n) + s2(-NH+1:NH,iblk)
   enddo
 
-  s=1.e-6*s
-  ss=1.e-6*ss
+  s=1.e-5*s
+  ss=1.e-5*ss
   df=1378.1285/NFFT
   nfac=40*mode65
   dtstep=0.25/df
 
-! Define range of frequencies to be searched
+  do i=1,2*mode65
+     call smo121(ss,16*NFFT)
+  enddo
+
+!  do i=-NH+1,NH
+!     write(72,3072) i*df,(ss(i,j),j=1,16)
+!3072 format(17f7.1)
+!  enddo
+
+! Define freq range to be searched. Upper tone is at sync freq  + 4*nfac*df Hz
   fa=-ntol
   fb=ntol
   ia2=max(-NH+1,nint(fa/df))
-! Upper tone is above sync tone by 4*nfac*df Hz
   ib2=min(NH,nint(fb/df + 4.1*nfac))
 
-! Find strongest line in each of the 4 phases, repeating for each drift rate.
+! Find strongest line in each of the 16 phases
   sbest=0.
   snrbest=0.
   nbest=1
   ipk=0
-
-  do n=1,8
+  do n=1,16
      sigmax(n)=0.
      do i=ia2,ib2
         sig=ss(i,n)
@@ -63,8 +65,8 @@ subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
         endif
      enddo
   enddo
-  n2best=nbest+4
-  if(n2best.gt.8) n2best=nbest-4
+  n2best=nbest+8
+  if(n2best.gt.16) n2best=nbest-8
   xdf=min(ipk(nbest),ipk(n2best))*df
   nspecial=0
   if(abs(xdf).le.ntol) then
@@ -74,6 +76,9 @@ subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
      iderr=nint((xk-k)*nfac)
 !     maxerr=nint(0.008*abs(idiff) + 0.51)
      maxerr=nint(0.02*abs(idiff) + 0.51)     !### Better test ??? ###
+!     write(71,3001) nbest,n2best,idiff,iderr,maxerr,k,      &
+!          ipk(nbest)*df,ipk(n2best)*df,sbest
+!3001 format(6i4,2f7.1,f7.2)
      if(abs(iderr).le.maxerr .and. k.ge.2 .and. k.le.4) nspecial=k
      snrdb=-30.0
      if(nspecial.gt.0) then