Bring the sort routines up-to-date with WSJT-X.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/map65@7497 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

libm65/CMakeLists.txt

@@ -111,7 +111,6 @@ set (FSRCS
   noisegen.f90
   packjt.f90
   pctile.f90
-  pctile2.f90
   pfxdump.f90
   qra64b.f90
   qra64c.f90
@@ -124,7 +123,6 @@ set (FSRCS
   shell.f90
   sleep_msec.f90
   smo.f90
-  sort.f90
   spec64.f90
   sun.f90
   symspec.f90
@@ -139,7 +137,6 @@ set (FSRCS
   zplot.f90
 
   f77_wisdom.f
-  ssort.f
   )
 
 set (CSRCS

libm65/ccf65.f90

@@ -12,7 +12,6 @@ subroutine ccf65(ss,nhsym,nfast,ssmax,sync1,ipol1,jpz,dt1,flipk,      &
   real pr2(NFFT)                   !JT65 shorthand pattern
   complex cpr2(0:NH)               !Complex FT of pr2
   real tmp1(322)
-  real tmp2(322)
   real ccf(-11:54,4)
   logical first
   integer npr(126)
@@ -62,7 +61,7 @@ subroutine ccf65(ss,nhsym,nfast,ssmax,sync1,ipol1,jpz,dt1,flipk,      &
 !        s(i)=ss(ip,i)+ss(ip,i+1)
         s(i)=min(ssmax,ss(ip,i)+ss(ip,i+1))
      enddo
-     call pctile(s,tmp1,nhsym-1,50,base)
+     call pctile(s,nhsym-1,50,base)
      s(1:nhsym-1)=s(1:nhsym-1)-base
      s(nhsym:NFFT)=0.
      call four2a(s,NFFT,1,-1,0)                !Real-to-complex FFT
@@ -115,7 +114,7 @@ subroutine ccf65(ss,nhsym,nfast,ssmax,sync1,ipol1,jpz,dt1,flipk,      &
   do i=1,nhsym
      tmp1(i)=ss(ipol2,i)
   enddo
-  call pctile(tmp1,tmp2,nhsym,40,base)
+  call pctile(tmp1,nhsym,40,base)
   snr2=0.398107*ccfbest2/base                !### empirical
   syncshort=0.5*ccfbest2/rms - 4.0           !### better normalizer than rms?
   dt2=(2.5 + lagpk2*(2048.0/11025.0))/nfast

libm65/display.f90

@@ -58,7 +58,7 @@ subroutine display(nkeep,ftol)
   endif
 
   call flush(26)
-  call indexx(nz,freqkHz,indx)
+  call indexx(freqkHz,nz,indx)
 
   nstart=1
   k3=0
@@ -81,7 +81,7 @@ subroutine display(nkeep,ftol)
         endif
         kz=k
         if(nstart.eq.1) then
-           call indexx(kz,utc2,indx2)
+           call indexx(utc2,kz,indx2)
            k3=0
            do k=1,kz
               k3=min(k3+1,400)
@@ -89,7 +89,7 @@ subroutine display(nkeep,ftol)
            enddo
            nstart=0
         else
-           call indexx(kz,utc2,indx2)
+           call indexx(utc2,kz,indx2)
            do k=1,kz
               k3=min(k3+1,400)
               line3(k3)=line2(indx2(k))
@@ -108,7 +108,7 @@ subroutine display(nkeep,ftol)
      j0=j
   enddo
   kz=k
-  call indexx(kz,utc2,indx2)
+  call indexx(utc2,kz,indx2)
   do k=1,kz
      k3=min(k3+1,400)
      line3(k3)=line2(indx2(k))

libm65/extract.f90

@@ -2,7 +2,6 @@ subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
 
   use packjt
   real s3(64,63)
-  real tmp(4032)
   character decoded*22
   integer era(51),dat4(12),indx(64)
   integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
@@ -23,7 +22,7 @@ subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
   save
 
   nfail=0
-  call pctile(s3,tmp,4032,50,base)     ! ### or, use ave from demod64a
+  call pctile(s3,4032,50,base)     ! ### or, use ave from demod64a
   s3=s3/base
   s3a=s3
 1 call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)

libm65/indexx.f90

@@ -1,19 +1,91 @@
-subroutine indexx(n,arr,indx)
+subroutine indexx(arr,n,indx)
 
-  parameter (NMAX=3000)
+  parameter (M=7,NSTACK=50)
-  integer indx(n)
+  integer n,indx(n)
   real arr(n)
-  real brr(NMAX)
+  integer i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
-  if(n.gt.NMAX) then
+  real a
-     print*,'n=',n,' too big in indexx.'
+
-     stop
+  do j=1,n
-  endif
+     indx(j)=j
-  do i=1,n
+  enddo
-     brr(i)=arr(i)
+
-     indx(i)=i
+  jstack=0
-  enddo
+  l=1
-  call ssort(brr,indx,n,2)
+  ir=n
+1 if(ir-l.lt.M) then
+     do j=l+1,ir
+        indxt=indx(j)
+        a=arr(indxt)
+        do i=j-1,1,-1
+           if(arr(indx(i)).le.a) goto 2
+           indx(i+1)=indx(i)
+        enddo
+        i=0
+2       indx(i+1)=indxt
+     enddo
+     if(jstack.eq.0) return
+
+     ir=istack(jstack)
+     l=istack(jstack-1)
+     jstack=jstack-2
+
+  else
+     k=(l+ir)/2
+     itemp=indx(k)
+     indx(k)=indx(l+1)
+     indx(l+1)=itemp
+
+     if(arr(indx(l+1)).gt.arr(indx(ir))) then
+        itemp=indx(l+1)
+        indx(l+1)=indx(ir)
+        indx(ir)=itemp
+     endif
+
+     if(arr(indx(l)).gt.arr(indx(ir))) then
+        itemp=indx(l)
+        indx(l)=indx(ir)
+        indx(ir)=itemp
+     endif
+
+     if(arr(indx(l+1)).gt.arr(indx(l))) then
+        itemp=indx(l+1)
+        indx(l+1)=indx(l)
+        indx(l)=itemp
+     endif
+
+     i=l+1
+     j=ir
+     indxt=indx(l)
+     a=arr(indxt)
+3    continue
+     i=i+1
+     if(arr(indx(i)).lt.a) goto 3
+
+4    continue
+     j=j-1
+     if(arr(indx(j)).gt.a) goto 4
+     if(j.lt.i) goto 5
+     itemp=indx(i)
+     indx(i)=indx(j)
+     indx(j)=itemp
+     goto 3
+
+5    indx(l)=indx(j)
+     indx(j)=indxt
+     jstack=jstack+2
+     if(jstack.gt.NSTACK) stop 'NSTACK too small in indexx'
+     if(ir-i+1.ge.j-l)then
+        istack(jstack)=ir
+        istack(jstack-1)=i
+        ir=j-1
+     else
+        istack(jstack)=j-1
+        istack(jstack-1)=l
+        l=i
+     endif
+  endif
+  goto 1
 
-  return
 end subroutine indexx
 

libm65/map65a.f90

@@ -12,7 +12,6 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb,        &
   real*4 ss(4,322,NFFT),savg(4,NFFT)
   real tavg(-50:50)                  !Temp for finding local base level
   real base(4)                       !Local basel level at 4 pol'ns
-  real tmp (200)                     !Temp storage for pctile sorting
   real sig(MAXMSG,30)                !Parameters of detected signals
   real a(5)
   real*8 fcenter
@@ -102,7 +101,7 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb,        &
                     go to 999
                  endif
               enddo
-              call pctile(tavg,tmp,101,50,base(jp))
+              call pctile(tavg,101,50,base(jp))
            enddo
         endif
 

libm65/pctile.f90

@@ -1,13 +1,22 @@
-subroutine pctile(x,tmp,nmax,npct,xpct)
+subroutine pctile(x,npts,npct,xpct)
-  real x(nmax),tmp(nmax)
+
-
+  parameter (NMAX=100000)
-  do i=1,nmax
+  real*4 x(npts)
-     tmp(i)=x(i)
+  real*4 tmp(NMAX)
-  enddo
+
-  call sort(nmax,tmp)
+  if(npts.le.0) then
-  j=nint(nmax*0.01*npct)
+     xpct=1.0
-  if(j.lt.1) j=1
+     go to 900
-  xpct=tmp(j)
+  endif
-
+  if(npts.gt.NMAX) stop
-  return
+
-end subroutine pctile
+  tmp(1:npts)=x
+  call shell(npts,tmp)
+  j=nint(npts*0.01*npct)
+  if(j.lt.1) j=1
+  if(j.gt.npts) j=npts
+  xpct=tmp(j)
+
+900 continue
+  return
+end subroutine pctile

libm65/pctile2.f90

@@ -1,22 +0,0 @@
-subroutine pctile2(x,npts,npct,xpct)
-
-  parameter (NMAX=100000)
-  real*4 x(npts)
-  real*4 tmp(NMAX)
-
-  if(npts.le.0) then
-     xpct=1.0
-     go to 900
-  endif
-  if(npts.gt.NMAX) stop
-
-  tmp(1:npts)=x
-  call shell(npts,tmp)
-  j=nint(npts*0.01*npct)
-  if(j.lt.1) j=1
-  if(j.gt.npts) j=npts
-  xpct=tmp(j)
-
-900 continue
-  return
-end subroutine pctile2

libm65/qra64c.f90

@@ -97,7 +97,7 @@ subroutine qra64c(cx,cy,nutc,nqd,ikhz,nfqso,ntol,xpol,mycall_12,     &
      if(mod(itry0,2).eq.0) idt=-idt
      jpk=jpk0 + 750*idt
      call spec64(c0,npts2,mode64,jpk,s3a,LL,NN)
-     call pctile2(s3a,LL*NN,40,base)
+     call pctile(s3a,LL*NN,40,base)
      s3a=s3a/base
      where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim
      do iter=itz,0,-2

libm65/spec64.f90

@@ -24,7 +24,7 @@ subroutine spec64(c0,npts2,mode64,jpk,s3,LL,NN)
 
   df=6000.0/nfft
   do i=1,LL
-     call pctile2(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
+     call pctile(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
   enddo
   
   nh=25

libm65/trimlist.f90

@@ -7,7 +7,7 @@ subroutine trimlist(sig,km,ftol,indx,nsiz,nz)
 !    1      2     3    4    5    6     7     8
 !  nfile  nutc  freq  snr  dt  ipol  flip  sync
 
-  call indexx(km,sig(1,3),indx)            !Sort list by frequency
+  call indexx(sig(1,3),km,indx)            !Sort list by frequency
 
   n=1
   i0=1

mainwindow.cpp

@@ -1,4 +1,4 @@
-//-------------------------------------------------------------- MainWindow
+//------------------------------------------------------------- MainWindow
 #include "mainwindow.h"
 #include "ui_mainwindow.h"
 #include "devsetup.h"