code:
/*
 * sleep.h  1.0 02/03/10
 *
 * Defines cross-platform sleep, usleep, etc.
 *
 * By Wu Yongwei
 *
 */

#ifndef _SLEEP_H
#define _SLEEP_H

#ifdef _WIN32
# if defined(_NEED_SLEEP_ONLY) && (defined(_MSC_VER) || defined(__MINGW32__))
#  include <stdlib.h>
#  define sleep(t) _sleep((t) * 1000)
# else
#  include <windows.h>
#  define sleep(t)  Sleep((t) * 1000)
# endif
# ifndef _NEED_SLEEP_ONLY
#  define msleep(t) Sleep(t)
#  define usleep(t) Sleep((t) / 1000)
# endif
#else
# include <unistd.h>
# ifndef _NEED_SLEEP_ONLY
#  define msleep(t) usleep((t) * 1000)
# endif
#endif

#endif /* _SLEEP_H */

license: mit


code:
subroutine sleep_msec(n)
  call usleep(n*1000)
  return
end subroutine sleep_msec


code:
subroutine slope(y,npts,xpk)

! Remove best-fit slope from data in y(i).  When fitting the straight line,
! ignore the peak around xpk +/- 4 bins

  real y(npts)

  sumw=0.
  sumx=0.
  sumy=0.
  sumx2=0.
  sumxy=0.
  sumy2=0.
  do i=1,npts
     if(abs(i-xpk).gt.4.0) then
        sumw=sumw + 1.0
        x=i
        sumx=sumx + x
        sumy=sumy + y(i)
        sumx2=sumx2 + x*x
        sumxy=sumxy + x*y(i)
        sumy2=sumy2 + y(i)**2
     endif
  enddo

  delta=sumw*sumx2 - sumx**2
  a=(sumx2*sumy - sumx*sumxy) / delta
  b=(sumw*sumxy - sumx*sumy) / delta

  sq=0.
  do i=1,npts
     y(i)=y(i)-(a + b*i)
     if(abs(i-xpk).gt.4.0) sq=sq + y(i)**2
  enddo
  rms=sqrt(sq/(sumw-4.0))
  y=y/rms

  return
end subroutine slope


code:
subroutine smo(x,npts,y,nadd)

  real x(npts)
  real y(npts)

  nh=nadd/2
  do i=1+nh,npts-nh
     sum=0.
     do j=-nh,nh
        sum=sum + x(i+j)
     enddo
     y(i)=sum
  enddo
  x=y
  x(:nh)=0.
  x(npts-nh+1:)=0.

  return
end subroutine smo


code:
subroutine smo121(x,nz)

  real x(nz)

  x0=x(1)
  do i=2,nz-1
     x1=x(i)
     x(i)=0.5*x(i) + 0.25*(x0+x(i+1))
     x0=x1
  enddo

  return
end subroutine smo121


code:
!program smolorentz
subroutine smolorentz(s1,nz,w,s2)

!  parameter (nz=512)
  real s1(nz),s2(nz)
  real y(-50:50)
!  character*8 arg

!  s1=0.
!  s1(256)=1.

!  call getarg(1,arg)
!  read(arg,*) w

  do i=-50,50
     x=i
     z=x/(0.5*w)
     y(i)=0.
     if(abs(z).lt.3.0) then
        d=1.0 + z*z
        y(i)=(1.0/d - 0.1)*10.0/9.0
     endif
  enddo

  jz=nint(1.5*w)
  if(jz.gt.50) jz=50
  do i=1,nz
     s=0.
     sy=0.
     do j=-jz,jz
        k=i+j
        if(k.ge.1 .and. k.le.nz) then
           s=s + s1(k)*y(j)
           sy=sy+y(j)
        endif
     enddo
     s2(i)=s/sy
!     write(52,3002) i-256,s1(i),s2(i)
!3002 format(i5,2f10.4)
  enddo

end subroutine smolorentz


code:
subroutine softsym(id2,npts8,nsps8,newdat,fpk,syncpk,snrdb,xdt,        &
     freq,drift,a3,schk,i1SoftSymbols)

! Compute the soft symbols

  use timer_module, only: timer

  parameter (NZ2=1512,NZ3=1360)
  logical, intent(inout) :: newdat
  complex c2(0:NZ2-1)
  complex c3(0:NZ3-1)
  complex c5(0:NZ3-1)
  real a(3)
  integer*1 i1SoftSymbolsScrambled(207)
  integer*1 i1SoftSymbols(207)
  include 'jt9sync.f90'

  nspsd=16
  ndown=nsps8/nspsd

! Mix, low-pass filter, and downsample to 16 samples per symbol
  call timer('downsam9',0)
  call downsam9(id2,npts8,nsps8,newdat,nspsd,fpk,c2)
  call timer('downsam9',1)

  call peakdt9(c2,nsps8,nspsd,c3,xdt)  !Find DT

  fsample=1500.0/ndown
  a=0.
  call timer('afc9    ',0)
  call afc9(c3,nz3,fsample,a,syncpk)  !Find deltaF, fDot, extra DT
  call timer('afc9    ',1)
  freq=fpk - a(1)
  drift=-2.0*a(2)
!  write(*,3301) fpk,freq,a
!3301 format(2f9.3,3f10.4)
  a3=a(3)
  a(3)=0.

  call timer('twkfreq ',0)
  call twkfreq(c3,c5,nz3,fsample,a)   !Correct for delta f, f1, f2 ==> a(1:3)
  call timer('twkfreq ',1)

! Compute soft symbols (in scrambled order)
  call timer('symspec2',0)
  call symspec2(c5,nz3,nsps8,nspsd,fsample,freq,drift,snrdb,schk,      &
       i1SoftSymbolsScrambled)
  call timer('symspec2',1)

! Remove interleaving
  call interleave9(i1SoftSymbolsScrambled,-1,i1SoftSymbols)

  return
end subroutine softsym


code:
subroutine softsym9f(ss2,ss3,i1SoftSymbols)

! Compute soft symbols and S/N

  real ss2(0:8,85)
  real ss3(0:7,69)
  integer*1 i1SoftSymbolsScrambled(207)
  integer*1 i1SoftSymbols(207)

  ss=0.
  sig=0.
  if(ss2(0,1).eq.-999.0) return       !Silence compiler warning
  do j=1,69
     smax=0.
     do i=0,7
        smax=max(smax,ss3(i,j))
        ss=ss+ss3(i,j)
     enddo
     sig=sig+smax
     ss=ss-smax
  enddo
  ave=ss/(69*7)                       !Baseline
!  call pctile(ss2,9*85,35,xmed)                       !### better? ###
  ss3=ss3/ave
  sig=sig/69.                         !Signal

  m0=3
  k=0
  scale=10.0
  do j=1,69
     do m=m0-1,0,-1                   !Get bit-wise soft symbols
        if(m.eq.2) then
           r1=max(ss3(4,j),ss3(5,j),ss3(6,j),ss3(7,j))
           r0=max(ss3(0,j),ss3(1,j),ss3(2,j),ss3(3,j))
        else if(m.eq.1) then
           r1=max(ss3(2,j),ss3(3,j),ss3(4,j),ss3(5,j))
           r0=max(ss3(0,j),ss3(1,j),ss3(6,j),ss3(7,j))
        else
           r1=max(ss3(1,j),ss3(2,j),ss3(4,j),ss3(7,j))
           r0=max(ss3(0,j),ss3(3,j),ss3(5,j),ss3(6,j))
        endif

        k=k+1
        i4=nint(scale*(r1-r0))
        if(i4.lt.-127) i4=-127
        if(i4.gt.127) i4=127
        i1SoftSymbolsScrambled(k)=i4
     enddo
  enddo

  call interleave9(i1SoftSymbolsScrambled,-1,i1SoftSymbols)

  return
end subroutine softsym9f


code:
subroutine softsym9w(id2,npts,xdt0,f0,width,nsubmode,xdt1,snrdb,i1softsymbols)

  parameter (NFFT=6912,NH=NFFT/2,NQ=NH/2)
  real s(NQ)
  real s2(0:8,85)
  real s3(0:7,69)
  real x(NFFT)
  complex cx(0:NH)
  integer*2 id2(60*12000)
  integer*1 i1SoftSymbolsScrambled(207)
  integer*1 i1softsymbols(207)
  include 'jt9sync.f90'
  equivalence (x,cx)

  if(npts.eq.-99) stop                     !Silence compiler warning
  df=12000.0/NFFT
  i0a=max(1.0,(xdt0-1.0)*12000.0)
  i0b=(xdt0+1.0)*12000.0
  k1=max(1,nint((f0-0.5*width)/df))
  k2=min(NQ,nint((f0+0.5*width)/df))
  smax=0.
  i0pk=1
  i1softsymbols=0

  do i0=i0a,i0b,432
     s=0.
     ssum=0.
     do j=1,16
        ia=i0 + (ii(j)-1)*nfft
        ib=ia+NFFT-1
        x=1.e-6*id2(ia:ib)
        call four2a(x,nfft,1,-1,0)        !r2c FFT
        do k=1,NQ
           s(k)=s(k) + real(cx(k))**2 + aimag(cx(k))**2
        enddo
     enddo
     ssum=ssum + sum(s(k1:k2))
     if(ssum.gt.smax) then
        smax=ssum
        i0pk=i0
     else
        if(ssum.lt.0.7*smax) exit
     endif
  end do
  xdt1=(i0pk-1)/12000.0

  if(i0pk.le.0) go to 999

  m=0
  do j=1,85
     ia=i0pk + (j-1)*nfft
     ib=ia+NFFT-1
     x=1.e-6*id2(ia:ib)
     call four2a(x,nfft,1,-1,0)        !r2c FFT
     do k=1,NQ
        s(k)=real(cx(k))**2 + aimag(cx(k))**2
     enddo

     dtone=df*(2**nsubmode)
     do i=0,8
        f=f0 + i*dtone
        k1=max(1,nint((f-0.5*width)/df))
        k2=min(NQ,nint((f+0.5*width)/df))
        s2(i,j)=sum(s(k1:k2))                  !Symbol spectra, including sync
     enddo

     if(isync(j).eq.0) then
        m=m+1
        s3(0:7,m)=s2(1:8,j)                   !Symbol spectra, data only
     endif

!     write(19,3101) j,s2(0:8,j)
!3101 format(i2,9f8.2)
  enddo

  ss=0.
  sig=0.
  do j=1,69
     smax=0.
     do i=0,7
        smax=max(smax,s3(i,j))
        ss=ss+s3(i,j)
     enddo
     sig=sig+smax
     ss=ss-smax
  enddo
  ave=ss/(69*7)                           !Baseline
  call pctile(s2,9*85,35,xmed)
  s3=s3/ave
  sig=sig/69.                             !Signal
  snrdb=db(sig/xmed) - 28.0

  m0=3
  k=0
  do j=1,69
     smax=0.
     do i=0,7
        if(s3(i,j).gt.smax) smax=s3(i,j)
     enddo

     do m=m0-1,0,-1                   !Get bit-wise soft symbols
        if(m.eq.2) then
           r1=max(s3(4,j),s3(5,j),s3(6,j),s3(7,j))
           r0=max(s3(0,j),s3(1,j),s3(2,j),s3(3,j))
        else if(m.eq.1) then
           r1=max(s3(2,j),s3(3,j),s3(4,j),s3(5,j))
           r0=max(s3(0,j),s3(1,j),s3(6,j),s3(7,j))
        else
           r1=max(s3(1,j),s3(2,j),s3(4,j),s3(7,j))
           r0=max(s3(0,j),s3(3,j),s3(5,j),s3(6,j))
        endif

        k=k+1
        i4=nint(10.0*(r1-r0))
        if(i4.lt.-127) i4=-127
        if(i4.gt.127) i4=127
        i1SoftSymbolsScrambled(k)=i4
     enddo
  enddo

! Remove interleaving
  call interleave9(i1SoftSymbolsScrambled,-1,i1SoftSymbols)

999  return
end subroutine softsym9w


code:
subroutine sort(n,arr)

  integer n,m,nstack
  real arr(n)
  parameter (m=7,nstack=50)
  integer i,ir,j,jstack,k,l,istack(nstack)
  real a,temp

  jstack=0
  l=1
  ir=n
  n0=n

1 if(ir-l.lt.m) then
     do j=l+1,ir
        a=arr(j)
        do i=j-1,1,-1
           if(arr(i).le.a) goto 2
           arr(i+1)=arr(i)
        enddo
        i=0
2       arr(i+1)=a
     enddo

     if(jstack.eq.0) return

     ir=istack(jstack)
     l=istack(jstack-1)
     jstack=jstack-2

  else
     k=(l+ir)/2
     temp=arr(k)
     arr(k)=arr(l+1)
     arr(l+1)=temp

     if(arr(l+1).gt.arr(ir)) then
        temp=arr(l+1)
        arr(l+1)=arr(ir)
        arr(ir)=temp
     endif

     if(arr(l).gt.arr(ir)) then
        temp=arr(l)
        arr(l)=arr(ir)
        arr(ir)=temp
     endif

     if(arr(l+1).gt.arr(l)) then
        temp=arr(l+1)
        arr(l+1)=arr(l)
        arr(l)=temp
     endif

     i=l+1
     j=ir
     a=arr(l)
3    i=i+1
     if(i.gt.n0) then
        do jj=1,n0
           write(99,3001) jj,arr(jj),i,n,ir
3001       format(i10,e12.3,3i10)
        enddo
        close(99)
        stop 'Bounds error in sort.f90'
     endif
     if(arr(i).lt.a) goto 3

4    j=j-1
     if(arr(j).gt.a) goto 4

     if(j.lt.i) goto 5
     temp=arr(i)
     arr(i)=arr(j)
     arr(j)=temp
     goto 3

5    arr(l)=arr(j)
     arr(j)=a
     jstack=jstack+2
     if(jstack.gt.nstack) stop 'nstack too small in sort'

     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

end subroutine sort


code:
subroutine spec64(c0,npts2,mode64,jpk,s3,LL,NN)

  parameter (NSPS=3456)                      !Samples per symbol at 6000 Hz
  complex c0(0:360000)                       !Complex spectrum of dd()
  complex cs(0:NSPS-1)                       !Complex symbol spectrum
  real s3(LL,NN)                             !Synchronized symbol spectra
  real xbase0(LL),xbase(LL)

  nfft=nsps
  fac=1.0/nfft
  do j=1,NN
     jj=j+7                                  !Skip first Costas array
     if(j.ge.33) jj=j+14                     !Skip middle Costas array
     ja=jpk + (jj-1)*nfft
     jb=ja+nfft-1
     cs(0:nfft-1)=fac*c0(ja:jb)
     call four2a(cs,nfft,1,-1,1)
     do ii=1,LL
        i=ii-65
        if(i.lt.0) i=i+nfft
        s3(ii,j)=real(cs(i))**2 + aimag(cs(i))**2
     enddo
  enddo

  df=6000.0/nfft
  do i=1,LL
     call pctile(s3(i,1:NN),NN,45,xbase0(i)) !Get baseline for passband shape
  enddo

  nh=25
  xbase(1:nh-1)=sum(xbase0(1:nh-1))/(nh-1.0)
  xbase(LL-nh+1:LL)=sum(xbase0(LL-nh+1:LL))/(nh-1.0)
  do i=nh,LL-nh
     xbase(i)=sum(xbase0(i-nh+1:i+nh))/(2*nh+1)  !Smoothed passband shape
  enddo

  do i=1,LL
     s3(i,1:NN)=s3(i,1:NN)/(xbase(i)+0.001) !Apply frequency equalization
  enddo

  return
end subroutine spec64


code:
subroutine spec9f(id2,npts,nsps,s1,jz,nq)

! Compute symbol spectra at quarter-symbol steps.

  integer*2 id2(0:npts)
  real s1(nq,jz)
  real x(960)
  complex c(0:480)
  equivalence (x,c)

  nfft=2*nsps                               !FFTs at twice the symbol length
  nh=nfft/2
  do j=1,jz
     ia=(j-1)*nsps/4
     ib=ia+nsps-1
     if(ib.gt.npts) exit
     x(1:nh)=id2(ia:ib)
     x(nh+1:)=0.
     call four2a(x,nfft,1,-1,0)           !r2c
     k=mod(j-1,340)+1
     do i=1,NQ
        s1(i,j)=1.e-10*(real(c(i))**2 + aimag(c(i))**2)
     enddo
  enddo

!### Reference spectrum should be applied here (or possibly earlier?) ###
!### Normalize so that rms (or level?) is 1.0 ?  ###

  return
end subroutine spec9f


code:
program stats

  character*8 arg
  character*40 infile
  character decoded*22

  nargs=iargc()
  if(nargs.lt.1) then
     print*,'Usage: stats file1 ...'
     go to 999
  endif

  ttol=0.1
  nftol=1
  write(*,1000)
1000 format(' SNR  Files  Sync    BM    FT  Hint Total  False BadSync'/  &
          56('-'))

  do ifile=1,nargs
     call getarg(ifile,infile)
     open(10,file=infile,status='old')
     i1=index(infile,".")+1
     i2=40
     if(index(infile,"_").gt.1) i2=index(infile,"_") - 1
     snrgen=0.
     read(infile(i1:i2),*,err=1) snrgen
1     snrgen=-snrgen
     nsynced=0
     nbm=0
     nftok=0
     nhint=0
     ngood=0
     nbad=0
     nbadsync=0

     do iline=1,999999
        read(10,1010,end=100) nutc,sync,nsnr,dt,nfreq,ncandidates,nhard,  &
             ntotal,rtt,ntry,nft,nqual,decoded
1010    format(i4.4,f5.1,i4,f5.1,i5,i6,i3,i4,f6.3,i8,i2,i3,1x,a22)

        ndfreq=abs(nfreq-1500)
        if(sync.ge.1.0 .and. abs(dt).le.ttol .and. ndfreq.le.nftol) then
           nsynced=nsynced+1
        else
           nbadsync=nbadsync+1
        endif

        if(decoded.eq.'                      ') cycle
        if(nft.eq.2 .or. (ntotal.le.81 .and. rtt.le.0.87)) then       !nag=0
           if(decoded(1:11).eq.'K1ABC W9XYZ') then
              ngood=ngood+1
              if(nft.eq.1 .and. ncandidates.eq.0) nbm=nbm+1
              if(nft.eq.1) nftok=nftok+1
              if(nft.ge.2) nhint=nhint+1
           else
              nbad=nbad+1
           endif
        endif
     enddo

100  write(*,1100) snrgen,nutc,nsynced,nbm,nftok,nhint,ngood,nbad,   &
          nbadsync
1100 format(f5.1,8i6)
  enddo

999 end program stats


code:
function stdmsg(msg0,bcontest,mygrid)

  ! Is msg0 a standard "JT-style" message?

  use iso_c_binding, only: c_bool
  use packjt
  character*22 msg0,msg1,msg
  character*6 mygrid
  integer dat(12)
  logical(c_bool), value :: bcontest
  logical(c_bool) :: stdmsg

  msg1=msg0
  i0=index(msg1,' OOO ')
  if(i0.gt.10) msg1=msg0(1:i0)
  call packmsg(msg0,dat,itype,logical(bcontest))
  call unpackmsg(dat,msg,logical(bcontest),mygrid)
  stdmsg=(msg.eq.msg1) .and. (itype.ge.0) .and. itype.ne.6

  return
end function stdmsg


code:
subroutine subtract65(dd,npts,f0,dt)

! Subtract a jt65 signal
!
! Measured signal  : dd(t)    = a(t)cos(2*pi*f0*t+theta(t))
! Reference signal : cref(t)  = exp( j*(2*pi*f0*t+phi(t)) )
! Complex amp      : cfilt(t) = LPF[ dd(t)*CONJG(cref(t)) ]
! Subtract         : dd(t)    = dd(t) - 2*REAL{cref*cfilt}

  use packjt
  use timer_module, only: timer

  integer correct(63)
  parameter (NMAX=60*12000) !Samples per 60 s
  parameter (NFILT=1600)
  real*4  dd(NMAX), window(-NFILT/2:NFILT/2)
  complex cref,camp,cfilt,cw
  integer nprc(126)
  real*8 dphi,phi
  logical first
  data nprc/                                   &
    1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
    0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
    0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
    0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
    1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
    0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
    1,1,1,1,1,1/
  data first/.true./
  common/chansyms65/correct
  common/heap1/cref(NMAX),camp(NMAX),cfilt(NMAX),cw(NMAX)
  save first

  pi=4.0*atan(1.0)

! Symbol duration is 4096/11025 s.
! Sample rate is 12000/s, so 12000*(4096/11025)=4458.23 samples/symbol.
! For now, call it 4458 samples/symbol. Over the message duration, we'll be off
! by about (4458.23-4458)*126=28.98 samples; 29 samples, or 0.7% of 1 symbol.
! Could eliminate accumulated error by injecting one extra sample every
! 5 or so symbols... Maybe try this later.

  nstart=dt*12000+1;
  nsym=126
  ns=4458
  nref=nsym*ns
  nend=nstart+nref-1
  phi=0.0
  iref=1
  ind=1
  isym=1
  call timer('subtr_1 ',0)
  do k=1,nsym
    if( nprc(k) .eq. 1 ) then
        omega=2*pi*f0
    else
        omega=2*pi*(f0+2.6917*(correct(isym)+2))
        isym=isym+1
    endif
    dphi=omega/12000.0
    do i=1,ns
        cref(ind)=cexp(cmplx(0.0,phi))
        phi=modulo(phi+dphi,2*pi)
        id=nstart-1+ind
        if(id.ge.1) camp(ind)=dd(id)*conjg(cref(ind))
        ind=ind+1
     enddo
  enddo
  call timer('subtr_1 ',1)

  call timer('subtr_2 ',0)
! Smoothing filter: do the convolution by means of FFTs. Ignore end-around
! cyclic effects for now.

  nfft=564480

  if(first) then
! Create and normalize the filter
     sum=0.0
     do j=-NFILT/2,NFILT/2
        window(j)=cos(pi*j/NFILT)**2
        sum=sum+window(j)
     enddo
     cw=0.
     do i=-NFILT/2,NFILT/2
        j=i+1
        if(j.lt.1) j=j+nfft
        cw(j)=window(i)/sum
     enddo
     call four2a(cw,nfft,1,-1,1)
     first=.false.
  endif

  nz=561708
  cfilt(1:nz)=camp(1:nz)
  cfilt(nz+1:nfft)=0.
  call four2a(cfilt,nfft,1,-1,1)
  fac=1.0/float(nfft)
  cfilt(1:nfft)=fac*cfilt(1:nfft)*cw(1:nfft)
  call four2a(cfilt,nfft,1,1,1)
  call timer('subtr_2 ',1)

! Subtract the reconstructed signal
  call timer('subtr_3 ',0)
  do i=1,nref
     j=nstart+i-1
     if(j.ge.1 .and. j.le.npts) dd(j)=dd(j)-2*REAL(cfilt(i)*cref(i))
  enddo
  call timer('subtr_3 ',1)

  return
end subroutine subtract65


code:
subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)

  implicit none

  integer y                         !Year
  integer m                         !Month
  integer DD                        !Day
  integer mjd                       !Modified Julian Date
  real UT                           !UT!in hours
  real RA,Dec                       !RA and Dec of sun

! NB: Double caps here are single caps in the writeup.

! Orbital elements of the Sun (also N=0, i=0, a=1):
  real w                            !Argument of perihelion
  real e                            !Eccentricity
  real MM                           !Mean anomaly
  real Ls                           !Mean longitude

! Other standard variables:
  real v                            !True anomaly
  real EE                           !Eccentric anomaly
  real ecl                          !Obliquity of the ecliptic
  real d                            !Ephemeris time argument in days
  real r                            !Distance to sun, AU
  real xv,yv                        !x and y coords in ecliptic
  real lonsun                       !Ecliptic long and lat of sun
!Ecliptic coords of sun (geocentric)
  real xs,ys
!Equatorial coords of sun (geocentric)
  real xe,ye,ze
  real lon,lat
  real GMST0,LST,HA
  real xx,yy,zz
  real xhor,yhor,zhor
  real Az,El

  real day
  real rad
  data rad/57.2957795/

! Time in days, with Jan 0, 2000 equal to 0.0:
  d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
  mjd=d + 51543
  ecl = 23.4393 - 3.563e-7 * d

! Compute updated orbital elements for Sun:
  w = 282.9404 + 4.70935e-5 * d
  e = 0.016709 - 1.151e-9 * d
  MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
  Ls = mod(w+MM+720.0,360.0)

  EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
  EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))

  xv = cos(EE/rad) - e
  yv = sqrt(1.0-e*e) * sin(EE/rad)
  v = rad*atan2(yv,xv)
  r = sqrt(xv*xv + yv*yv)
  lonsun = mod(v + w + 720.0,360.0)
! Ecliptic coordinates of sun (rectangular):
  xs = r * cos(lonsun/rad)
  ys = r * sin(lonsun/rad)

! Equatorial coordinates of sun (rectangular):
  xe = xs
  ye = ys * cos(ecl/rad)
  ze = ys * sin(ecl/rad)

! RA and Dec in degrees:
  RA = rad*atan2(ye,xe)
  Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))

  GMST0 = (Ls + 180.0)/15.0
  LST = mod(GMST0+UT+lon/15.0+48.0,24.0)    !LST in hours
  HA = 15.0*LST - RA                        !HA in degrees
  xx = cos(HA/rad)*cos(Dec/rad)
  yy = sin(HA/rad)*cos(Dec/rad)
  zz =             sin(Dec/rad)
  xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
  yhor = yy
  zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
  Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
  El = rad*asin(zhor)
  day=d-1.5

  return
end subroutine sun


code:
subroutine symspec(shared_data,k,ntrperiod,nsps,ingain,nminw,pxdb,s,   &
     df3,ihsym,npts8,pxdbmax)

! Input:
!  k         pointer to the most recent new data
!  ntrperiod T/R sequence length, minutes
!  nsps      samples per symbol, at 12000 Hz
!  ndiskdat  0/1 to indicate if data from disk
!  nb        0/1 status of noise blanker (off/on)
!  nbslider  NB setting, 0-100

! Output:
!  pxdb      raw power (0-90 dB)
!  s()       current spectrum for waterfall display
!  ihsym     index number of this half-symbol (1-184) for 60 s modes

! jt9com
!  ss()      JT9 symbol spectra at half-symbol steps
!  savg()    average spectra for waterfall display

  use, intrinsic :: iso_c_binding, only: c_int, c_short, c_float, c_char
  include 'jt9com.f90'

  type(dec_data) :: shared_data
  real*4 w3(MAXFFT3)
  real*4 s(NSMAX)
  real*4 ssum(NSMAX)
  real*4 xc(0:MAXFFT3-1)
  real*4 tmp(NSMAX)
  complex cx(0:MAXFFT3/2)
  integer nch(7)

  common/spectra/syellow(NSMAX),ref(0:3456),filter(0:3456)
  data k0/99999999/,nfft3z/0/
  data nch/1,2,4,9,18,36,72/
  equivalence (xc,cx)
  save

  if(ntrperiod.eq.-999) stop                   !Silence compiler warning
  nfft3=16384                                  !df=12000.0/16384 = 0.732422
  jstep=nsps/2                                 !Step size = half-symbol in id2()
  if(k.gt.NMAX) go to 900
  if(k.lt.2048) then                !(2048 was nfft3)  (Any need for this ???)
     ihsym=0
     go to 900                                 !Wait for enough samples to start
  endif

  if(nfft3.ne.nfft3z) then
! Compute new window
     pi=4.0*atan(1.0)
     width=0.25*nsps
     do i=1,nfft3
        z=(i-nfft3/2)/width
        w3(i)=exp(-z*z)
     enddo
     nfft3z=nfft3
  endif

  if(k.lt.k0) then                             !Start a new data block
     ja=0
     ssum=0.
     ihsym=0
     if(.not. shared_data%params%ndiskdat) shared_data%id2(k+1:)=0   !Needed to prevent "ghosts". Not sure why.
  endif
  gain=10.0**(0.1*ingain)
  sq=0.
  pxmax=0.;
  do i=k0+1,k
     x1=shared_data%id2(i)
     if (abs(x1).gt.pxmax) pxmax = abs(x1);
     sq=sq + x1*x1
  enddo
  pxdb = 0.
  if(sq.gt.0.0) pxdb=10*log10(sq/(k-k0))
  pxdbmax=0.
  if(pxmax.gt.0) pxdbmax = 20*log10(pxmax)

  k0=k
  ja=ja+jstep                         !Index of first sample

  fac0=0.1
  do i=0,nfft3-1                      !Copy data into cx
     j=ja+i-(nfft3-1)
     xc(i)=0.
     if(j.ge.1 .and.j.le.NMAX) xc(i)=fac0*shared_data%id2(j)
  enddo
  ihsym=ihsym+1

!  xc(0:nfft3-1)=w3(1:nfft3)*xc(0:nfft3-1)    !Apply window w3
  call four2a(xc,nfft3,1,-1,0)               !Real-to-complex FFT

  df3=12000.0/nfft3                   !JT9-1: 0.732 Hz = 0.42 * tone spacing
  iz=min(NSMAX,nint(5000.0/df3))
  fac=(1.0/nfft3)**2
  do i=1,iz
     j=i-1
     if(j.lt.0) j=j+nfft3
     sx=fac*(real(cx(j))**2 + aimag(cx(j))**2)
     if(ihsym.le.184) shared_data%ss(ihsym,i)=sx
     ssum(i)=ssum(i) + sx
     s(i)=1000.0*gain*sx
  enddo

  shared_data%savg=ssum/ihsym

  if(mod(ihsym,10).eq.0) then
     mode4=nch(nminw+1)
     nsmo=min(10*mode4,150)
     nsmo=4*nsmo
     call flat1(shared_data%savg,iz,nsmo,syellow)
     if(mode4.ge.2) call smo(syellow,iz,tmp,mode4)
     if(mode4.ge.2) call smo(syellow,iz,tmp,mode4)
     syellow(1:250)=0.
     ia=500./df3
     ib=2700.0/df3
     smin=minval(syellow(ia:ib))
     smax=maxval(syellow(1:iz))
     syellow=(50.0/(smax-smin))*(syellow-smin)
     where(syellow<0) syellow=0.
  endif

900 npts8=k/8

  return
end subroutine symspec


code:
subroutine symspec2(c5,nz3,nsps8,nspsd,fsample,freq,drift,snrdb,schk,    &
     i1SoftSymbolsScrambled)

! Compute soft symbols from the final downsampled data

  complex c5(0:4096-1)
  complex z
  integer*1 i1SoftSymbolsScrambled(207)
  real aa(3)
  real ss2(0:8,85)
  real ss3(0:7,69)
  include 'jt9sync.f90'
  data scale/10.0/

  aa(1)=-1500.0/nsps8
  aa(2)=0.
  aa(3)=0.
  do i=0,8                                         !Loop over the 9 tones
     if(i.ge.1) call twkfreq(c5,c5,nz3,fsample,aa)
     m=0
     k=-1
     do j=1,85                                     !Loop over all symbols
        z=0.
        do n=1,nspsd                               !Sum over 16 samples
           k=k+1
           z=z+c5(k)
        enddo
        ss2(i,j)=real(z)**2 + aimag(z)**2        !Symbol speactra, data and sync
        if(i.ge.1 .and. isync(j).eq.0) then
           m=m+1
           ss3(i-1,m)=ss2(i,j)                   !Symbol speactra, data only
        endif
     enddo
  enddo

  call chkss2(ss2,freq,drift,schk)

  ss=0.
  sig=0.
  do j=1,69
     smax=0.
     do i=0,7
        smax=max(smax,ss3(i,j))
        ss=ss+ss3(i,j)
     enddo
     sig=sig+smax
     ss=ss-smax
  enddo
  ave=ss/(69*7)                           !Baseline
  call pctile(ss2,9*85,35,xmed)
  ss3=ss3/ave
  sig=sig/69.                             !Signal
  t=max(1.0,sig - 1.0)
  snrdb=db(t) - 61.3

  m0=3
  k=0
  do j=1,69
        smax=0.
        do i=0,7
           if(ss3(i,j).gt.smax) smax=ss3(i,j)
        enddo

     do m=m0-1,0,-1                   !Get bit-wise soft symbols
        if(m.eq.2) then
           r1=max(ss3(4,j),ss3(5,j),ss3(6,j),ss3(7,j))
           r0=max(ss3(0,j),ss3(1,j),ss3(2,j),ss3(3,j))
        else if(m.eq.1) then
           r1=max(ss3(2,j),ss3(3,j),ss3(4,j),ss3(5,j))
           r0=max(ss3(0,j),ss3(1,j),ss3(6,j),ss3(7,j))
        else
           r1=max(ss3(1,j),ss3(2,j),ss3(4,j),ss3(7,j))
           r0=max(ss3(0,j),ss3(3,j),ss3(5,j),ss3(6,j))
        endif

        k=k+1
        i4=nint(scale*(r1-r0))
        if(i4.lt.-127) i4=-127
        if(i4.gt.127) i4=127
!        i4=i4+128
!        if(i4.le.127) i1SoftSymbolsScrambled(k)=i4
!        if(i4.ge.128) i1SoftSymbolsScrambled(k)=i4-256
        i1SoftSymbolsScrambled(k)=i4
     enddo
  enddo

  return
end subroutine symspec2


code:
!subroutine symspec65(dd,npts,ss,nqsym,savg)
subroutine symspec65(dd,npts,nqsym,savg)

! Compute JT65 symbol spectra at quarter-symbol steps

  parameter (NFFT=8192)
  parameter (NSZ=3413)                         !NFFT*5000/12000
  parameter (MAXHSYM=322)
  parameter (MAXQSYM=552)
  real*8 hstep
  real*4 dd(npts)
!  real*4 ss(MAXHSYM,NSZ)
  real*4 ss(MAXQSYM,NSZ)
  real*4 savg(NSZ)
  real*4 x(NFFT)
  real*4 w(NFFT)
  complex c(0:NFFT/2)
  logical first
  common/refspec/dfref,ref(NSZ)
  equivalence (x,c)
  data first/.true./
  save /refspec/,first,w
  common/sync/ss

  hstep=2048.d0*12000.d0/11025.d0              !half-symbol = 2229.116 samples
  qstep=hstep/2.0              !quarter-symbol = 1114.558 samples
  nsps=nint(2*hstep)
  df=12000.0/NFFT
  nhsym=(npts-NFFT)/hstep
  nqsym=(npts-NFFT)/qstep
  savg=0.
  fac1=1.e-3

  if(first) then
! Compute the FFT window
!     width=0.25*nsps
     do i=1,NFFT
!        z=(i-NFFT/2)/width
        w(i)=1
        if(i.gt.4458) w(i)=0
!        w(i)=exp(-z*z)
      enddo
     first=.false.
  endif

  do j=1,nqsym
     i0=(j-1)*qstep
     x=fac1*w*dd(i0+1:i0+NFFT)
     call four2a(c,NFFT,1,-1,0)                !r2c forward FFT
     do i=1,NSZ
        s=real(c(i))**2 + aimag(c(i))**2
        ss(j,i)=s
        savg(i)=savg(i)+s
     enddo
  enddo
  savg=savg/nhsym

!  call flat65(ss,nhsym,MAXQSYM,NSZ,ref)  !Flatten the 2d spectrum, saving
  call flat65(ss,nqsym,MAXQSYM,NSZ,ref)  !Flatten the 2d spectrum, saving
  dfref=df                               ! the reference spectrum ref()
  savg=savg/ref
!  do j=1,nhsym
  do j=1,nqsym
     ss(j,1:NSZ)=ss(j,1:NSZ)/ref
  enddo

  return
end subroutine symspec65


code:
subroutine sync4(dat,jz,ntol,NFreeze,MouseDF,mode,mode4,minwidth,    &
     dtx,dfx,snrx,snrsync,ccfblue,ccfred1,flip,width,ps0)

! Synchronizes JT4 data, finding the best-fit DT and DF.

  parameter (NFFTMAX=2520)         !Max length of FFTs
  parameter (NHMAX=NFFTMAX/2)      !Max length of power spectra
  parameter (NSMAX=525)            !Max number of half-symbol steps
  integer ntol                     !Range of DF search
  real dat(jz)
  real psavg(NHMAX)                !Average spectrum of whole record
  real ps0(450)                    !Avg spectrum for plotting
  real s2(NHMAX,NSMAX)             !2d spectrum, stepped by half-symbols
  real ccfblue(-5:540)             !CCF with pseudorandom sequence
  real ccfred(-450:450)            !Peak of ccfblue, as function of freq
  real red(-450:450)               !Peak of ccfblue, as function of freq
  real ccfred1(-224:224)           !Peak of ccfblue, as function of freq
  real tmp(1260)
  integer ipk1(1)
  integer nch(7)
  logical savered
  equivalence (ipk1,ipk1a)
  data nch/1,2,4,9,18,36,72/
  save

!  write(*,3001) 'A',ntol,nfreeze,mousedf,mode,mode4,minwidth
!3001 format(a1,6i6)

! Do FFTs of twice symbol length, stepped by half symbols.  Note that
! we have already downsampled the data by factor of 2.
  nsym=207
  nfft=2520
  nh=nfft/2
  nq=nfft/4
  nsteps=jz/nq - 1
  df=0.5*11025.0/nfft
  psavg(1:nh)=0.
  if(mode.eq.-999) width=0.                        !Silence compiler warning

  do j=1,nsteps                     !Compute spectrum for each step, get average
     k=(j-1)*nq + 1
     call ps4(dat(k),nfft,s2(1,j))
     psavg(1:nh)=psavg(1:nh) + s2(1:nh,j)
  enddo

  nsmo=min(10*mode4,150)
  call flat1b(psavg,nsmo,s2,nh,nsteps,NHMAX,NSMAX)        !Flatten spectra

  if(mode4.ge.9) call smo(psavg,nh,tmp,mode4/4)

  i0=132
  do i=1,450
     ps0(i)=5.0*(psavg(i0+2*i) + psavg(i0+2*i+1) - 2.0)
  enddo

! Set freq and lag ranges
  famin=200.0 + 3*mode4*df
  fbmax=2700.0 - 3*mode4*df
  fa=famin
  fb=fbmax
  if(NFreeze.eq.1) then
     fa=max(famin,1270.46+MouseDF-ntol)
     fb=min(fbmax,1270.46+MouseDF+ntol)
  else
     fa=max(famin,1270.46+MouseDF-600)
     fb=min(fbmax,1270.46+MouseDF+600)
  endif
  ia=fa/df - 3*mode4                   !Index of lowest tone, bottom of range
  ib=fb/df - 3*mode4                   !Index of lowest tone, top of range
  i0=nint(1270.46/df)
  irange=450
  if(ia-i0.lt.-irange) ia=i0-irange
  if(ib-i0.gt.irange)  ib=i0+irange
  lag1=-5
  lag2=59
  syncbest=-1.e30
  ccfred=0.
  jmax=-1000
  jmin=1000

  do ich=minwidth,7                       !Find best width
     kz=nch(ich)/2
     savered=.false.
     do i=ia+kz,ib-kz                     !Find best frequency channel for CCF
        call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0,   &
             lagpk0,flip)
        j=i-i0 + 3*mode4
        if(j.ge.-372 .and. j.le.372) then
           ccfred(j)=ccf0
           jmax=max(j,jmax)
           jmin=min(j,jmin)
        endif

! Find rms of the CCF, without main peak
        call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0)
        sync=abs(ccfblue(lagpk0))

! Find best sync value
        if(sync.gt.syncbest) then
           ipk=i
           lagpk=lagpk0
           ichpk=ich
           syncbest=sync
           savered=.true.
        endif
     enddo
     if(savered) red=ccfred
  enddo

  ccfred=red
!  width=df*nch(ichpk)
  dfx=(ipk-i0 + 3*mode4)*df

! Peak up in time, at best whole-channel frequency
  call xcor4(s2,ipk,nsteps,nsym,lag1,lag2,ichpk,mode4,ccfblue,ccfmax,   &
       lagpk,flip)
  xlag=lagpk
  if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then
     call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
     xlag=lagpk+dx2
  endif

! Find rms of the CCF, without the main peak
  call slope(ccfblue(lag1),lag2-lag1+1,xlag-lag1+1.0)
  sq=0.
  nsq=0
  do lag=lag1,lag2
     if(abs(lag-xlag).gt.2.0) then
        sq=sq+ccfblue(lag)**2
        nsq=nsq+1
     endif
  enddo
  rms=sqrt(sq/nsq)
  snrsync=max(0.0,db(abs(ccfblue(lagpk)/rms - 1.0)) - 4.5)
  snrx=-26.
  if(mode4.eq.2)  snrx=-25.
  if(mode4.eq.4)  snrx=-24.
  if(mode4.eq.9)  snrx=-23.
  if(mode4.eq.18) snrx=-22.
  if(mode4.eq.36) snrx=-21.
  if(mode4.eq.72) snrx=-20.
  snrx=snrx + snrsync

  dt=2.0/11025.0
  istart=xlag*nq
  dtx=istart*dt
  ccfred1=0.
  jmin=max(jmin,-224)
  jmax=min(jmax,224)
  do i=jmin,jmax
     ccfred1(i)=ccfred(i)
  enddo

  ipk1=maxloc(ccfred1) - 225
  ns=0
  s=0.
  iw=min(mode4,(ib-ia)/4)
  do i=jmin,jmax
     if(abs(i-ipk1a).gt.iw) then
        s=s+ccfred1(i)
        ns=ns+1
     endif
  enddo
  base=s/ns
  ccfred1=ccfred1-base
  ccf10=0.5*maxval(ccfred1)
  do i=ipk1a,jmin,-1
     if(ccfred1(i).le.ccf10) exit
  enddo
  i1=i
  do i=ipk1a,jmax
     if(ccfred1(i).le.ccf10) exit
  enddo
  width=(i-i1)*df

  return
end subroutine sync4

include 'flat1b.f90'


code:
subroutine sync64(c0,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk,sync,  &
     sync2,width)

  use timer_module, only: timer

  parameter (NMAX=60*12000)                  !Max size of raw data at 12000 Hz
  parameter (NSPS=3456)                      !Samples per symbol at 6000 Hz
  parameter (NSPC=7*NSPS)                    !Samples per Costas array
  real s1(0:NSPC-1)                          !Power spectrum of Costas 1
  real s2(0:NSPC-1)                          !Power spectrum of Costas 2
  real s3(0:NSPC-1)                          !Power spectrum of Costas 3
  real s0(0:NSPC-1)                          !Sum of s1+s2+s3
  real s0a(0:NSPC-1)                         !Best synchromized spectrum (saved)
  real s0b(0:NSPC-1)                         !tmp
  real a(5)
  integer icos7(0:6)                         !Costas 7x7 tones
  integer ipk0(1)
  complex cc(0:NSPC-1)                       !Costas waveform
  complex c0(0:720000)                       !Complex spectrum of dd()
  complex c1(0:NSPC-1)                       !Complex spectrum of Costas 1
  complex c2(0:NSPC-1)                       !Complex spectrum of Costas 2
  complex c3(0:NSPC-1)                       !Complex spectrum of Costas 3
  data icos7/4,2,5,6,1,3,0/                  !Costas 7x7 tone pattern
  data mode64z/-1/
  save

  if(mode64.ne.mode64z) then
     twopi=8.0*atan(1.0)
     dfgen=mode64*12000.0/6912.0
     k=-1
     phi=0.
     do j=0,6                               !Compute complex Costas waveform
        dphi=twopi*10.0*icos7(j)*dfgen/6000.0
        do i=1,NSPS
           phi=phi + dphi
           if(phi.gt.twopi) phi=phi-twopi
           k=k+1
           cc(k)=cmplx(cos(phi),sin(phi))
        enddo
     enddo
     mode64z=mode64
  endif

  nfft3=NSPC
  nh3=nfft3/2
  df3=6000.0/nfft3

  fa=max(nf1,nfqso-ntol)
  fb=min(nf2,nfqso+ntol)
  iaa=max(0,nint(fa/df3))
  ibb=min(NSPC-1,nint(fb/df3))

  maxtol=max(ntol,500)
  fa=max(nf1,nfqso-maxtol)
  fb=min(nf2,nfqso+maxtol)
  ia=max(0,nint(fa/df3))
  ib=min(NSPC-1,nint(fb/df3))
  id=0.1*(ib-ia)
  iz=ib-ia+1
  sync=-1.e30
  smaxall=0.
  jpk=0
  ja=0
  jb=(5.0+emedelay)*6000
  jstep=100
  ipk=0
  kpk=0
  nadd=10*mode64
  if(mod(nadd,2).eq.0) nadd=nadd+1       !Make nadd odd
  nskip=max(49,nadd)

  do j1=ja,jb,jstep
     call timer('sync64_1',0)
     j2=j1 + 39*NSPS
     j3=j1 + 77*NSPS
     c1=1.e-4*c0(j1:j1+NSPC-1) * conjg(cc)
     c2=1.e-4*c0(j2:j2+NSPC-1) * conjg(cc)
     c3=1.e-4*c0(j3:j3+NSPC-1) * conjg(cc)
     call four2a(c1,nfft3,1,-1,1)
     call four2a(c2,nfft3,1,-1,1)
     call four2a(c3,nfft3,1,-1,1)
     s1=0.
     s2=0.
     s3=0.
     s0b=0.
     do i=ia,ib
        freq=i*df3
        s1(i)=real(c1(i))**2 + aimag(c1(i))**2
        s2(i)=real(c2(i))**2 + aimag(c2(i))**2
        s3(i)=real(c3(i))**2 + aimag(c3(i))**2
     enddo
     call timer('sync64_1',1)

     call timer('sync64_2',0)
     s0(ia:ib)=s1(ia:ib) + s2(ia:ib) + s3(ia:ib)
     s0(:ia-1)=0.
     s0(ib+1:)=0.
     if(nadd.ge.3) then
        do ii=1,3
           s0b(ia:ib)=s0(ia:ib)
           call smo(s0b(ia:ib),iz,s0(ia:ib),nadd)
        enddo
     endif
     call averms(s0(ia+id:ib-id),iz-2*id,nskip,ave,rms)
     s=(maxval(s0(iaa:ibb))-ave)/rms
     ipk0=maxloc(s0(iaa:ibb))
     ip=ipk0(1) + iaa - 1
     if(s.gt.sync) then
        jpk=j1
        s0a=(s0-ave)/rms
        sync=s
        dtx=jpk/6000.0 - 1.0
        ipk=ip
        f0=ip*df3
     endif
     call timer('sync64_2',1)
  enddo

  s0a=s0a+2.0
!  write(17) ia,ib,s0a(ia:ib)                !Save data for red curve
!  close(17)

  nskip=50
  call lorentzian(s0a(ia+nskip:ib-nskip),iz-2*nskip,a)
  f0a=(a(3)+ia+49)*df3
  w1=df3*a(4)
  w2=2*nadd*df3
  width=w1
  if(w1.gt.1.2*w2) width=sqrt(w1**2 - w2**2)

  sq=0.
  do i=1,20
     j=ia+nskip+1
     k=ib-nskip-21+i
     sq=sq + (s0a(j)-a(1))**2 + (s0a(k)-a(1))**2
  enddo
  rms2=sqrt(sq/40.0)
  sync2=10.0*log10(a(2)/rms2)

  slimit=6.0
  rewind 17
  write(17,1110) 0.0,0.0
  rewind 17
!  rewind 76
  do i=2,iz-2*nskip-1,3
     x=i
     z=(x-a(3))/(0.5*a(4))
     yfit=a(1)
     if(abs(z).lt.3.0) then
        d=1.0 + z*z
        yfit=a(1) + a(2)*(1.0/d - 0.1)
     endif
     j=i+ia+49
     freq=j*df3
     ss=(s0a(j-1)+s0a(j)+s0a(j+1))/3.0
     if(ss.gt.slimit) write(17,1110) freq,ss
1110 format(3f10.3)
!     write(76,1110) freq,ss,yfit
  enddo
  flush(17)
  close(17)
!  flush(76)

  return
end subroutine sync64


code:
subroutine sync65(nfa,nfb,naggressive,ntol,nqsym,ca,ncand,nrobust,   &
     bVHF)

  parameter (NSZ=3413,NFFT=8192,MAXCAND=300)
  real ss(552,NSZ)
  real ccfblue(-32:82)             !CCF with pseudorandom sequence
  real ccfred(NSZ)                  !Peak of ccfblue, as function of freq
  logical bVHF

  type candidate
     real freq
     real dt
     real sync
     real flip
  end type candidate
  type(candidate) ca(MAXCAND)

  common/steve/thresh0
  common/sync/ss

  if(ntol.eq.-99) stop                       !Silence compiler warning
  call setup65

  df=12000.0/NFFT                            !df = 12000.0/8192 = 1.465 Hz
  ia=max(2,nint(nfa/df))
  ib=min(NSZ-1,nint(nfb/df))
!  lag1=-11
!  lag2=59
!  lag1=-22
!  lag2=118
  lag1=-32
  lag2=82        !may need to be extended for EME
  nsym=126
  ncand=0
  fdot=0.
  ccfred=0.
  ccfblue=0.
  ccfmax=0.
  ipk=0
  do i=ia,ib
     call xcor(i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk0,flip,fdot,nrobust)
! Remove best-fit slope from ccfblue and normalize so baseline rms=1.0
     if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1,      &
          lagpk0-lag1+1.0)
     ccfred(i)=ccfblue(lagpk0)
     if(ccfred(i).gt.ccfmax) then
        ccfmax=ccfred(i)
        ipk=i
     endif
  enddo
  call pctile(ccfred(ia:ib),ib-ia+1,35,xmed)
  ccfred(ia:ib)=ccfred(ia:ib)-xmed
  ccfred(ia-1)=ccfred(ia)
  ccfred(ib+1)=ccfred(ib)

  do i=ia,ib
     freq=i*df
     itry=0
!     if(naggressive.gt.0 .and. ntol.lt.1000 .and. ccfmax.ge.thresh0) then
     if(naggressive.gt.0 .and. ccfmax.ge.thresh0) then
        if(i.ne.ipk) cycle
        itry=1
        ncand=ncand+1
     else
        if(ccfred(i).ge.thresh0 .and. ccfred(i).gt.ccfred(i-1) .and.       &
             ccfred(i).gt.ccfred(i+1)) then
           itry=1
           ncand=ncand+1
        endif
     endif
     if(itry.ne.0) then
        call xcor(i,nqsym,nsym,lag1,lag2,ccfblue,ccf0,lagpk,flip,fdot,nrobust)
        if(.not.bVHF) call slope(ccfblue(lag1),lag2-lag1+1,       &
             lagpk-lag1+1.0)
        xlag=lagpk
        if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then
           call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
           xlag=lagpk+dx2
        endif
        dtx=xlag*1024.0/11025.0
        ccfblue(lag1)=0.
        ccfblue(lag2)=0.
        ca(ncand)%freq=freq
        ca(ncand)%dt=dtx
        ca(ncand)%flip=flip
        if(bVHF) then
           ca(ncand)%sync=db(ccfred(i)) - 16.0
        else
           ca(ncand)%sync=ccfred(i)
        endif
     endif
     if(ncand.eq.MAXCAND) exit
  enddo

  return
end subroutine sync65


code:
subroutine sync9(ss,nzhsym,lag1,lag2,ia,ib,ccfred,red2,ipkbest)

  include 'constants.f90'
  real ss(184,NSMAX)
  real ss1(184)
  real ccfred(NSMAX)
  real savg(NSMAX)
  real savg2(NSMAX)
  real smo(-5:25)
  real sq(NSMAX)
  real red2(NSMAX)
  character*27 cr
  data cr/'(C) 2016, Joe Taylor - K1JT'/
  include 'jt9sync.f90'

  ipk=0
  ipkbest=0
  sbest=0.
  ccfred=0.

  do i=ia,ib                         !Loop over freq range
     ss1=ss(1:184,i)
     call pctile(ss1,nzhsym,40,xmed)

     ss1=ss1/xmed - 1.0
     do j=1,nzhsym
        if(ss1(j).gt.3.0) ss1(j)=3.0
     enddo

     call pctile(ss1,nzhsym,45,sbase)
     ss1=ss1-sbase
     sq0=dot_product(ss1(1:nzhsym),ss1(1:nzhsym))
     rms=sqrt(sq0/(nzhsym-1))

     smax=0.
     do lag=lag1,lag2                !DT = 2.5 to 5.0 s
        sum1=0.
        sq2=sq0
        nsum=nzhsym
        do j=1,16                    !Sum over 16 sync symbols
           k=ii2(j) + lag
           if(k.ge.1 .and. k.le.nzhsym) then
              sum1=sum1 + ss1(k)
              sq2=sq2 - ss1(k)*ss1(k)
              nsum=nsum-1
           endif
        enddo
        if(sum1.gt.smax) then
           smax=sum1
           ipk=i
        endif
        rms=sqrt(sq2/(nsum-1))
     enddo
     ccfred(i)=smax                        !Best at this freq, over all lags
     if(smax.gt.sbest) then
        sbest=smax
        ipkbest=ipk
     endif
  enddo

  call pctile(ccfred(ia),ib-ia+1,50,xmed)
  if(xmed.le.0.0) xmed=1.0
  ccfred=2.0*ccfred/xmed

  savg=0.
  do j=1,nzhsym
     savg(ia:ib)=savg(ia:ib) + ss(j,ia:ib)
  enddo
  savg(ia:ib)=savg(ia:ib)/nzhsym
  smo(0:20)=1.0/21.0
  smo(-5:-1)=-(1.0/21.0)*(21.0/10.0)
  smo(21:25)=smo(-5)

  do i=ia,ib
     sm=0.
     do j=-5,25
        if(i+j.ge.1 .and. i+j.lt.NSMAX) sm=sm + smo(j)*savg(i+j)
     enddo
     savg2(i)=sm
     sq(i)=sm*sm
  enddo

  call pctile(sq(ia:ib),ib-ia+1,20,sq0)
  rms=sqrt(sq0)
  savg2(ia:ib)=savg2(ia:ib)/(5.0*rms)

  red2=0.
  do i=ia+11,ib-10
     ref=max(savg2(i-10),savg2(i+10))
     red2(i)=savg2(i)-ref
     if(red2(i).lt.-99.0) red2(i)=-99.0
     if(red2(i).gt.99.0) red2(i)=99.0
  enddo

  return
end subroutine sync9


code:
subroutine sync9f(s2,nq,nfa,nfb,ss2,ss3,lagpk,ipk,ccfbest)

! Look for JT9 sync pattern in the folded symbol spectra, s2.
! Frequency search extends from nfa to nfb.  Synchronized symbol
! spectra are put into ss2() and ss3().

  integer ii4(16)
  real s2(240,340)
  real ss2(0:8,85)
  real ss3(0:7,69)
  include 'jt9sync.f90'

  ii4=4*ii-3
  ccf=0.
  ccfbest=0.
  nfft=4*nq
  df=12000.0/nfft
  ia=nfa/df
  ib=nfb/df + 0.9999

  do i=ia,ib
     do lag=0,339
        t=0.
        do n=1,16
           j=ii4(n)+lag
           if(j.gt.340) j=j-340
           t=t + s2(i,j)
        enddo
        if(t.gt.ccfbest) then
           lagpk=lag
           ipk=i
           ccfbest=t
        endif
     enddo
  enddo

  do i=0,8
     j4=lagpk-4
     i2=2*i + ipk
     if(i2.lt.1) i2=1
     m=0
     do j=1,85
        j4=j4+4
        if(j4.gt.340) j4=j4-340
        if(j4.lt.1) j4=j4+340
        ss2(i,j)=s2(i2,j4)
        if(i.ge.1 .and. isync(j).eq.0) then
           m=m+1
           ss3(i-1,m)=ss2(i,j)
        endif
     enddo
  enddo

  return
end subroutine sync9f


code:
subroutine sync9w(ss,nzhsym,lag1,lag2,ia,ib,ccfred,ccfblue,ipkbest,lagpk,nadd)

  include 'constants.f90'
  real ss(184,NSMAX)
  real ss1(184),ss1save(184)
  real ccfred(NSMAX)
  real ccfblue(-9:18)
  real sa(NSMAX),sb(NSMAX)
  include 'jt9sync.f90'

! Smooth the symbol spectra (by an amount consistent with measured width??)
  do j=1,nzhsym
     sa=ss(j,1:NSMAX)
     call smo(sa,NSMAX,sb,nadd)
     call smo(sb,NSMAX,sa,nadd)
     ss(j,1:NSMAX)=sa
  enddo

  ipk=0
  ipkbest=0
  sbest=0.
  ccfred=0.
  df=12000.0/16384.0

  do i=ia,ib                                 !Loop over specified freq range
     ss1=ss(1:184,i)                         !Symbol amplitudes at this freq
     call pctile(ss1,nzhsym,50,xmed)         !Median level at this freq
     ss1=ss1/xmed - 1.0

     smax=0.                                 !Find DT in specified range
     do lag=lag1,lag2
        sum1=0.
        nsum=nzhsym
        do j=1,16                            !Sum over 16 sync symbols
           k=ii2(j) + lag
           if(k.ge.1 .and. k.le.nzhsym) then
              sum1=sum1 + ss1(k)
              nsum=nsum-1
           endif
        enddo
        if(sum1.gt.smax) then
           smax=sum1
           ipk=i
        endif
     enddo

     ccfred(i)=smax                        !Best at this freq, over all lags
     if(smax.gt.sbest) then
        sbest=smax
        ipkbest=ipk
        ss1save=ss1
     endif
  enddo

  call pctile(ccfred(ia),ib-ia+1,50,xmed)
  if(xmed.le.0.0) xmed=1.0
  ccfred=ccfred/xmed

  ss1=ss1save
  smax=0.                                 !Find DT in specified range
  do lag=lag1,lag2
     sum1=0.
     nsum=nzhsym
     do j=1,16                            !Sum over 16 sync symbols
        k=ii2(j) + lag
        if(k.ge.1 .and. k.le.nzhsym) then
           sum1=sum1 + ss1(k)
           nsum=nsum-1
        endif
     enddo
     ccfblue(lag)=sum1
     if(sum1.gt.smax) then
        smax=sum1
        lagpk=lag
     endif
  enddo
  if(lagpk.eq.-9) lagpk=-8                !Protect the ends of ccfblue()
  if(lagpk.eq.18) lagpk=17

  return
end subroutine sync9w


code:
subroutine synciscat(cdat,npts,nh,npct,s0,jsym,df,ntol,NFreeze,    &
     MouseDF,mousebutton,mode4,nafc,psavg,xsync,sig,ndf0,msglen,         &
     ipk,jpk,idf,df1)

! Synchronize an ISCAT signal
! cdat() is the downsampled analytic signal.
! Sample rate = fsample = BW = 11025 * (9/32) = 3100.78125 Hz
! npts, nsps, etc., are all reduced by 9/32

  parameter (NMAX=30*3101)
  parameter (NSZ=4*1400)
  complex cdat(NMAX)
  complex c(288)
  real s0(288,NSZ)
  real fs0(288,96)                        !108 = 96 + 3*4
  real savg(288)
  real sref(288)
  real psavg(72)                          !Average spectrum of whole file
  integer icos(4)
  data icos/0,1,3,2/
  data nsync/4/,nlen/2/,ndat/18/

! Silence compiler warnings:
  sigbest=-20.0
  ndf0best=0
  msglenbest=0
  ipkbest=0
  jpkbest=0
  ipk2=0
  idfbest=mousebutton

  fsample=3100.78125                   !New sample rate
  nsps=144/mode4
  nsym=npts/nsps - 1
  nblk=nsync+nlen+ndat
  nfft=2*nsps                          !FFTs at twice the symbol length,

  kstep=nsps/4                         !  stepped by 1/4 symbol
  df=fsample/nfft
  fac=1.0/1000.0                       !Somewhat arbitrary
  savg=0.
  s0=0.

  ia=1-kstep
  do j=1,4*nsym                                   !Compute symbol spectra
     ia=ia+kstep
     ib=ia+nsps-1
     if(ib.gt.npts) exit
     c(1:nsps)=fac*cdat(ia:ib)
     c(nsps+1:nfft)=0.
     call four2a(c,nfft,1,-1,1)
     do i=1,nfft
        s0(i,j)=real(c(i))**2 + aimag(c(i))**2
        savg(i)=savg(i) + s0(i,j)                 !Accumulate avg spectrum
     enddo
     i0=40
  enddo

  jsym=4*nsym
  savg=savg/jsym

  do i=1,71                                   !Compute spectrum in dB, for plot
     if(mode4.eq.1) then
        psavg(i)=2*db(savg(4*i)+savg(4*i-1)+savg(4*i-2)+savg(4*i-3)) + 1.0
     else
        psavg(i)=2*db(savg(2*i)+savg(2*i-1)) + 7.0
     endif
  enddo

  do i=nh+1,nfft-nh
     call pctile(savg(i-nh),2*nh+1,npct,sref(i))
  enddo
  sref(1:nh)=sref(nh+11)
  sref(nfft-nh+1:nfft)=sref(nfft-nh)

  do i=1,nfft                                 !Normalize the symbol spectra
     fac=1.0/sref(i)
     if(i.lt.11) fac=1.0/savg(11)
     do j=1,jsym
        s0(i,j)=fac*s0(i,j)
     enddo
  enddo

  nfold=jsym/96
  jb=96*nfold

  ttot=npts/fsample                         !Length of record (s)
  df1=df/ttot                               !Step size for f1=fdot
  idf1=-25.0/df1
  idf2=5.0/df1
  if(nafc.eq.0) then
     idf1=0
     idf2=0
  else if(mod(-idf1,2).eq.1) then
     idf1=idf1-1
  endif

  xsyncbest=0.
  do idf=idf1,idf2                         !Loop over fdot
     fs0=0.
     do j=1,jb                             !Fold s0 into fs0, modulo 4*nblk
        k=mod(j-1,4*nblk)+1
        ii=nint(idf*float(j-jb/2)/float(jb))
        ia=max(1-ii,1)
        ib=min(nfft-ii,nfft)
        do i=ia,ib
           fs0(i,k)=fs0(i,k) + s0(i+ii,j)
        enddo
     enddo
     ref=nfold*4

     i0=27
     if(mode4.eq.1) i0=95
     ia=i0-nint(ntol/df)
     ib=i0+nint(ntol/df)
     if(ia.lt.1) ia=1
     if(ib.gt.nfft-3) ib=nfft-3

     smax=0.
     ipk=1
     jpk=1
     do j=0,4*nblk-1                            !Find sync pattern: lags 0-95
        do i=ia,ib                              !Search specified freq range
           ss=0.
           do n=1,4                             !Sum over 4 sync tones
              k=j+4*n-3
              if(k.gt.96) k=k-96
              ss=ss + fs0(i+2*icos(n),k)
           enddo
           if(ss.gt.smax) then
              smax=ss
              ipk=i                             !Frequency offset, DF
              jpk=j+1                           !Time offset, DT
           endif
        enddo
     enddo

     xsync=smax/ref - 1.0
     if(nfold.lt.26) xsync=xsync * sqrt(nfold/26.0)
     xsync=xsync-0.5                           !Empirical

     sig=db(smax/ref - 1.0) - 15.0
     if(mode4.eq.1) sig=sig-5.0
!     if(sig.lt.-20 .or. xsync.lt.1.0) sig=-20.0
!     if(sig.lt.-20) sig=-20.0
     ndf0=nint(df*(ipk-i0))

     smax=0.
     ja=jpk+16
     if(ja.gt.4*nblk) ja=ja-4*nblk
     jj=jpk+20
     if(jj.gt.4*nblk) jj=jj-4*nblk
     do i=ipk,ipk+60,2                         !Find User's message length
        ss=fs0(i,ja) + fs0(i+10,jj)
        if(ss.gt.smax) then
           smax=ss
           ipk2=i
        endif
     enddo

     msglen=(ipk2-ipk)/2
     if(msglen.lt.2 .or. msglen.gt.29) cycle

     if(xsync.ge.xsyncbest) then
        xsyncbest=xsync
        sigbest=sig
        ndf0best=ndf0
        msglenbest=msglen
        ipkbest=ipk
        jpkbest=jpk
        idfbest=idf
     endif
  enddo

  xsync=xsyncbest
  sig=sigbest
  ndf0=ndf0best
  msglen=msglenbest
  ipk=ipkbest
  jpk=jpkbest
  idf=idfbest
  if(nafc.eq.0) idf=0

  return
end subroutine synciscat


code:
subroutine syncmsk(cdat,npts,jpk,ipk,idf,rmax,snr,metric,decoded)

! Attempt synchronization, and if successful decode using Viterbi algorithm.

  use iso_c_binding, only: c_loc,c_size_t
  use packjt
  use hashing
  use timer_module, only: timer

  parameter (NSPM=1404,NSAVE=2000)
  complex cdat(npts)                    !Analytic signal
  complex cb(66)                        !Complex waveform for Barker-11 code
  complex cd(0:11,0:3)
  complex c(0:NSPM-1)                   !Complex data for one message length
  complex c2(0:NSPM-1)
  complex cb3(1:NSPM,3)
  real r(12000)
  real rdat(12000)
  real ss1(12000)
  real symbol(234)
  real rdata(198)
  real rd2(198)
  real rsave(NSAVE)
  real xp(29)
  complex z,z0,z1,z2,z3,cfac
  integer*1 e1(198)
  integer*1, target :: d8(13)
  integer*1 i1hash(4)
  integer*1 i1
  integer*4 i4Msg6BitWords(12)            !72-bit message as 6-bit words
  integer mettab(0:255,0:1)               !Metric table for BPSK modulation
  integer ipksave(NSAVE)
  integer jpksave(NSAVE)
  integer indx(NSAVE)
  integer b11(11)                      !Barker-11 code
  character*22 decoded
  character*72 c72
  logical first
  equivalence (i1,i4)
  equivalence (ihash,i1hash)
  data xp/0.500000, 0.401241, 0.309897, 0.231832, 0.168095,    &
          0.119704, 0.083523, 0.057387, 0.039215, 0.026890,    &
          0.018084, 0.012184, 0.008196, 0.005475, 0.003808,    &
          0.002481, 0.001710, 0.001052, 0.000789, 0.000469,    &
          0.000329, 0.000225, 0.000187, 0.000086, 0.000063,    &
          0.000017, 0.000091, 0.000032, 0.000045/
  data first/.true./
  data b11/1,1,1,0,0,0,1,0,0,1,0/
  save first,cb,cd,twopi,dt,f0,f1,mettab

  phi=0.
  if(first) then
! Get the metric table
     bias=0.0
     scale=20.0
     xln2=log(2.0)
     mettab=0
     do i=128,156
        x0=log(max(0.001,2.0*xp(i-127)))/xln2
        x1=log(max(0.001,2.0*(1.0-xp(i-127))))/xln2
        mettab(i,0)=nint(scale*(x0-bias))
        mettab(i,1)=nint(scale*(x1-bias))
        mettab(256-i,0)=mettab(i,1)
        mettab(256-i,1)=mettab(i,0)
     enddo
     do i=157,255
        mettab(i,0)=mettab(156,0)
        mettab(i,1)=mettab(156,1)
        mettab(256-i,0)=mettab(i,1)
        mettab(256-i,1)=mettab(i,0)
     enddo
     j=0
     twopi=8.0*atan(1.0)
     dt=1.0/12000.0
     f0=1000.0
     f1=2000.0
     dphi=0
     do i=1,11
        if(b11(i).eq.0) dphi=twopi*f0*dt
        if(b11(i).eq.1) dphi=twopi*f1*dt
        do n=1,6
           j=j+1
           phi=phi+dphi
           cb(j)=cmplx(cos(phi),sin(phi))
        enddo
     enddo
     cb3=0.
     cb3(1:66,1)=cb
     cb3(283:348,1)=cb
     cb3(769:834,1)=cb

     cb3(1:66,2)=cb
     cb3(487:552,2)=cb
     cb3(1123:1188,2)=cb

     cb3(1:66,3)=cb
     cb3(637:702,3)=cb
     cb3(919:984,3)=cb

     phi=0.
     do n=0,3
        k=-1
        dphi=twopi*f0*dt
        if(n.ge.2) dphi=twopi*f1*dt
        do i=0,5
           k=k+1
           phi=phi+dphi
           if(phi.gt.twopi) phi=phi-twopi
           cd(k,n)=cmplx(cos(phi),sin(phi))
        enddo

        dphi=twopi*f0*dt
        if(mod(n,2).eq.1) dphi=twopi*f1*dt
        do i=6,11
           k=k+1
           phi=phi+dphi
           if(phi.gt.twopi) phi=phi-twopi
           cd(k,n)=cmplx(cos(phi),sin(phi))
        enddo
     enddo

     first=.false.
  endif

  nfft=NSPM
  jz=npts-nfft
  decoded="                      "
  ipk=0
  jpk=0
  metric=-9999
  r=0.

  call timer('sync1   ',0)
  do j=1,jz                               !Find the Barker-11 sync vectors
     z=0.
     ss=0.
     do i=1,66
        ss=ss + real(cdat(j+i-1))**2 + aimag(cdat(j+i-1))**2
        z=z + cdat(j+i-1)*conjg(cb(i))    !Signal matching Barker 11
     enddo
     ss=sqrt(ss/66.0)*66.0
     r(j)=abs(z)/(0.908*ss)               !Goodness-of-fit to Barker 11
     ss1(j)=ss
  enddo
  call timer('sync1   ',1)

  call timer('sync2   ',0)
  jz=npts-nfft
  rmax=0.
!  n1=35, n2=69, n3=94
  k=0
  do j=1,jz                               !Find best full-message sync
     if(ss1(j).lt.85.0) cycle
     r1=r(j) + r(j+282) + r(j+768)        ! 6*(12+n1) 6*(24+n1+n2)
     r2=r(j) + r(j+486) + r(j+1122)       ! 6*(12+n2) 6*(24+n2+n3)
     r3=r(j) + r(j+636) + r(j+918)        ! 6*(12+n3) 6*(24+n3+n1)
     if(r1.gt.rmax) then
        rmax=r1
        jpk=j
        ipk=1
     endif
     if(r2.gt.rmax) then
        rmax=r2
        jpk=j
        ipk=2
     endif
     if(r3.gt.rmax) then
        rmax=r3
        jpk=j
        ipk=3
     endif
     rrmax=max(r1,r2,r3)
     if(rrmax.gt.1.9) then
        k=min(k+1,NSAVE)
        if(r1.eq.rrmax) ipksave(k)=1
        if(r2.eq.rrmax) ipksave(k)=2
        if(r3.eq.rrmax) ipksave(k)=3
        jpksave(k)=j
        rsave(k)=rrmax
     endif
  enddo
  call timer('sync2   ',1)
  kmax=k

  call indexx(rsave,kmax,indx)

  call timer('sync3   ',0)
  do kk=1,kmax
     k=indx(kmax+1-kk)
     ipk=ipksave(k)
     jpk=jpksave(k)
     rmax=rsave(k)

     c=conjg(cb3(1:NSPM,ipk))*cdat(jpk:jpk+nfft-1)
     smax=0.
     dfx=0.
     idfbest=0
     do itry=1,25
        idf=itry/2
        if(mod(itry,2).eq.0) idf=-idf
        idf=4*idf
        twk=idf
        call tweak1(c,NSPM,-twk,c2)
        z=sum(c2)
        if(abs(z).gt.smax) then
           dfx=twk
           smax=abs(z)
           phi=atan2(aimag(z),real(z))            !Carrier phase offset
           idfbest=idf
        endif
     enddo
     idf=idfbest
     call tweak1(cdat,npts,-dfx,cdat)
     cfac=cmplx(cos(phi),-sin(phi))
     cdat=cfac*cdat

     sig=0.
     ref=0.
     rdat(1:npts)=cdat
     iz=11
     do k=1,234                                !Compute soft symbols
        j=jpk+6*(k-1)

        z0=2.0*dot_product(cdat(j:j+iz),cd(0:iz,0))
        z1=2.0*dot_product(cdat(j:j+iz),cd(0:iz,1))
        z2=2.0*dot_product(cdat(j:j+iz),cd(0:iz,2))
        z3=2.0*dot_product(cdat(j:j+iz),cd(0:iz,3))

!### Maybe these should be weighted by yellow() ?
        if(j+1404+iz.lt.npts) then
           z0=z0 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,0))
           z1=z1 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,1))
           z2=z2 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,2))
           z3=z3 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,3))
        endif

        if(j-1404.ge.1) then
           z0=z0 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,0))
           z1=z1 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,1))
           z2=z2 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,2))
           z3=z3 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,3))
        endif

        sym=max(abs(real(z2)),abs(real(z3))) - max(abs(real(z0)),abs(real(z1)))

        if(sym.lt.0.0) then
           phi=atan2(aimag(z0),real(z0))
           sig=sig + real(z0)**2
           ref=ref + aimag(z0)**2
        else
           phi=atan2(aimag(z1),real(z1))
           sig=sig + real(z1)**2
           ref=ref + aimag(z1)**2
        endif
        n=k
        if(ipk.eq.2) n=k+47
        if(ipk.eq.3) n=k+128
        if(n.gt.234) n=n-234
        ibit=0
        if(sym.ge.0) ibit=1
        symbol(n)=sym
     enddo
     snr=db(sig/ref-1.0)

     rdata(1:35)=symbol(12:46)
     rdata(36:104)=symbol(59:127)
     rdata(105:198)=symbol(140:233)

! Re-order the symbols and make them i*1
     j=0
     do i=1,99
        i4=128+rdata(i)                       !### Should be nint() ??? ###
        if(i4.gt.255)  i4=255
        if(i4.lt.0) i4=0
        j=j+1
        e1(j)=i1
        rd2(j)=rdata(i)
        i4=128+rdata(i+99)
        if(i4.gt.255)  i4=255
        if(i4.lt.0) i4=0
        j=j+1
        e1(j)=i1
        rd2(j)=rdata(i+99)
     enddo

! Decode the message
     nb1=87
     call vit213(e1,nb1,mettab,d8,metric)
     ihash=nhash(c_loc(d8),int(9,c_size_t),146)
     ihash=2*iand(ihash,32767)
     decoded='                      '
     if(d8(10).eq.i1hash(2) .and. d8(11).eq.i1hash(1)) then
        write(c72,1012) d8(1:9)
1012    format(9b8.8)
        read(c72,1014) i4Msg6BitWords
1014    format(12b6.6)
        call unpackmsg(i4Msg6BitWords,decoded,.false.,'      ') !Unpack to get msgsent
     endif
     if(decoded.ne.'                      ') exit
  enddo
  call timer('sync3   ',1)

  return
end subroutine syncmsk


code:
program t6

  parameter (MAXFFT=1404)
  complex c(0:MAXFFT-1)
  real s(0:MAXFFT-1)

  m1=45
  m2=67
  m3=89
  nsym=3*11 + m1 + m2 + m3
  nfft=6*nsym
  nh=nfft/2

  best=9999.
!  do m1=22,67
!     do m2=37,97
  do m1=30,67
     do m2=26,100
        m3=201-m2-m1
        if(m3.lt.13) cycle
        c=0.
        n1=6*(11+m1)
        n2=n1+6*(11+m2)
        c(1:66)=1.
        c(1+n1:66+n1)=1.
        c(1+n2:66+n2)=1.

        call four2a(c,nfft,1,-1,1)            !c2c FFT

        df=12000.0/nfft
        smax=0.
        do i=0,nfft-1
           s(i)=real(c(i))**2 + aimag(c(i))**2
           if(i.ne.0) smax=max(s(i),smax)
        enddo
        sidelobe=db(smax/s(0))

        if(sidelobe.lt.best) then
           write(*,1000) m1,m2,m3,sidelobe
1000       format(3i5,f8.2)
           best=sidelobe
           s=s/s(0)
           rewind 13
           do j=0,nfft-1
              i=mod(j+nh,nfft)
              f=i*df
              if(i.gt.nh) f=f-12000.0
              write(13,1020) f,s(i)
1020          format(2f12.4)
           enddo
        endif
     enddo
  enddo

end program t6


code:
program test_init_random_seed
  real :: r(10,4)
  call init_random_seed()
  call random_number(r)
  do i =1,10
     write (*, *) (r(i,j),j=1,4)
  end do
end program test_init_random_seed


code:
program testfast9

  parameter (NMAX=359424)
  integer*2 id2(NMAX)
  integer narg(0:11)
  character*80 line(100)
  character submode*1,infile*80

  nargs=iargc()
  if(nargs.ne.2) then
     print*,'Usage:    testfast9 submode infile'
     print*,'Example: testfast9 E /data/VE1SKY/K1JT/JT9E/150806_123300.wav'
     go to 999
  endif
  call getarg(1,submode)
  call getarg(2,infile)

  open(10,file=infile,access='stream',status='old')
  read(10) id2(1:22)                     !Skip 44 header bytes
  npts=NMAX
  read(10,end=1) id2(1:npts)                   !Read the raw data

1 i1=index(infile,'.wav')
  read(infile(i1-6:i1-1),*) narg(0)
  narg(1)=NMAX
  n=ichar(submode)
  narg(2)=n-ichar('A')
  if(n.ge.97 .and. n.le.104) narg(2)=n-ichar('a')
  narg(3)=1
  narg(4)=0
  narg(5)=0
  narg(6)=0
  narg(7)=29951
  narg(8)=1
  narg(9)=102
  narg(10)=700
  narg(11)=500

  call fast9(id2,narg,line)
  print*,line(1)

999 end program testfast9


code:
  parameter (MAXTEST=75,NTEST=68)
  character*22 testmsg(MAXTEST)
  character*22 testmsgchk(MAXTEST)
  ! Test msgs should include the extremes for the different types
  ! See pfx.f90
  ! Type 1 P & A
  ! Type 1 1A & E5
  data testmsg(1:NTEST)/         &
       "CQ WB9XYZ EN34",         &
       "CQ DX WB9XYZ EN34",      &
       "QRZ WB9XYZ EN34",        &
       "KA1ABC WB9XYZ EN34",     &
       "KA1ABC WB9XYZ RO",       &
       "KA1ABC WB9XYZ -21",      &
       "KA1ABC WB9XYZ R-19",     &
       "KA1ABC WB9XYZ RRR",      &
       "KA1ABC WB9XYZ 73",       &
       "KA1ABC WB9XYZ",          &
       "CQ 000 WB9XYZ EN34",     &
       "CQ 999 WB9XYZ EN34",     &
       "CQ EU WB9XYZ EN34",      &
       "CQ WY WB9XYZ EN34",      &
       "1A/KA1ABC WB9XYZ",       &
       "E5/KA1ABC WB9XYZ",       &
       "KA1ABC 1A/WB9XYZ",       &
       "KA1ABC E5/WB9XYZ",       &
       "KA1ABC/P WB9XYZ",        &
       "KA1ABC/A WB9XYZ",        &
       "KA1ABC WB9XYZ/P",        &
       "KA1ABC WB9XYZ/A",        &
       "CQ KA1ABC/P",            &
       "CQ WB9XYZ/A",            &
       "QRZ KA1ABC/P",           &
       "QRZ WB9XYZ/A",           &
       "DE KA1ABC/P",            &
       "DE WB9XYZ/A",            &
       "CQ 1A/KA1ABC",           &
       "CQ E5/KA1ABC",           &
       "DE 1A/KA1ABC",           &
       "DE E5/KA1ABC",           &
       "QRZ 1A/KA1ABC",          &
       "QRZ E5/KA1ABC",          &
       "CQ WB9XYZ/1A",           &
       "CQ WB9XYZ/E5",           &
       "QRZ WB9XYZ/1A",          &
       "QRZ WB9XYZ/E5",          &
       "DE WB9XYZ/1A",           &
       "DE WB9XYZ/E5",           &
       "CQ A000/KA1ABC FM07",    &
       "CQ ZZZZ/KA1ABC FM07",    &
       "QRZ W4/KA1ABC FM07",     &
       "DE W4/KA1ABC FM07",      &
       "CQ W4/KA1ABC -22",       &
       "DE W4/KA1ABC -22",       &
       "QRZ W4/KA1ABC -22",      &
       "CQ W4/KA1ABC R-22",      &
       "DE W4/KA1ABC R-22",      &
       "QRZ W4/KA1ABC R-22",     &
       "DE W4/KA1ABC 73",        &
       "CQ KA1ABC FM07",         &
       "QRZ KA1ABC FM07",        &
       "DE KA1ABC/VE6 FM07",     &
       "CQ KA1ABC/VE6 -22",      &
       "DE KA1ABC/VE6 -22",      &
       "QRZ KA1ABC/VE6 -22",     &
       "CQ KA1ABC/VE6 R-22",     &
       "DE KA1ABC/VE6 R-22",     &
       "QRZ KA1ABC/VE6 R-22",    &
       "DE KA1ABC 73",           &
       "HELLO WORLD",            &
       "ZL4/KA1ABC 73",          &
       "KA1ABC XL/WB9XYZ",       &
       "KA1ABC WB9XYZ/W4",       &
       "DE KA1ABC/QRP 2W",       &
       "KA1ABC/1 WB9XYZ/1",      &
       "123456789ABCDEFGH"/
  data testmsgchk(1:NTEST)/      &
       "CQ WB9XYZ EN34",         &
       "CQ DX WB9XYZ EN34",      &
       "QRZ WB9XYZ EN34",        &
       "KA1ABC WB9XYZ EN34",     &
       "KA1ABC WB9XYZ RO",       &
       "KA1ABC WB9XYZ -21",      &
       "KA1ABC WB9XYZ R-19",     &
       "KA1ABC WB9XYZ RRR",      &
       "KA1ABC WB9XYZ 73",       &
       "KA1ABC WB9XYZ",          &
       "CQ 000 WB9XYZ EN34",     &
       "CQ 999 WB9XYZ EN34",     &
       "CQ EU WB9XYZ EN34",      &
       "CQ WY WB9XYZ EN34",      &
       "1A/KA1ABC WB9XYZ",       &
       "E5/KA1ABC WB9XYZ",       &
       "KA1ABC 1A/WB9XYZ",       &
       "KA1ABC E5/WB9XYZ",       &
       "KA1ABC/P WB9XYZ",        &
       "KA1ABC/A WB9XYZ",        &
       "KA1ABC WB9XYZ/P",        &
       "KA1ABC WB9XYZ/A",        &
       "CQ KA1ABC/P",            &
       "CQ WB9XYZ/A",            &
       "QRZ KA1ABC/P",           &
       "QRZ WB9XYZ/A",           &
       "DE KA1ABC/P",            &
       "DE WB9XYZ/A",            &
       "CQ 1A/KA1ABC",           &
       "CQ E5/KA1ABC",           &
       "DE 1A/KA1ABC",           &
       "DE E5/KA1ABC",           &
       "QRZ 1A/KA1ABC",          &
       "QRZ E5/KA1ABC",          &
       "CQ WB9XYZ/1A",           &
       "CQ WB9XYZ/E5",           &
       "QRZ WB9XYZ/1A",          &
       "QRZ WB9XYZ/E5",          &
       "DE WB9XYZ/1A",           &
       "DE WB9XYZ/E5",           &
       "CQ A000/KA1ABC FM07",    &
       "CQ ZZZZ/KA1ABC FM07",    &
       "QRZ W4/KA1ABC FM07",     &
       "DE W4/KA1ABC FM07",      &
       "CQ W4/KA1ABC -22",       &
       "DE W4/KA1ABC -22",       &
       "QRZ W4/KA1ABC -22",      &
       "CQ W4/KA1ABC R-22",      &
       "DE W4/KA1ABC R-22",      &
       "QRZ W4/KA1ABC R-22",     &
       "DE W4/KA1ABC 73",        &
       "CQ KA1ABC FM07",         &
       "QRZ KA1ABC FM07",        &
       "DE KA1ABC/VE6 FM07",     &
       "CQ KA1ABC/VE6 -22",      &
       "DE KA1ABC/VE6 -22",      &
       "QRZ KA1ABC/VE6 -22",     &
       "CQ KA1ABC/VE6 R-22",     &
       "DE KA1ABC/VE6 R-22",     &
       "QRZ KA1ABC/VE6 R-22",    &
       "DE KA1ABC 73",           &
       "HELLO WORLD",            &
       "ZL4/KA1ABC 73",          &
       "KA1ABC XL/WB9",          &
       "KA1ABC WB9XYZ",          &
       "DE KA1ABC/QRP",          &
       "KA1ABC/1 WB9X",          &
       "123456789ABCD"/


code:
program timefft

! Tests and times one-dimensional FFTs computed by FFTW3

  use, intrinsic :: iso_c_binding
  use FFTW3

  complex(C_FLOAT_COMPLEX),pointer :: a(:),b(:),c(:)
  real(C_FLOAT),pointer :: ar(:),br(:)
  type(C_PTR) :: plan1,plan2              !Pointers to FFTW plans
  type(C_PTR) :: pa,pb,pc
  integer(C_INT) iret
  integer*8 count0,count1,clkfreq
  character problem*9
  logical linplace,lcomplex,lthreading

! Get command-line parameters
  call timefft_opts(npatience,maxthreads,linplace,lcomplex,nfft,problem,nflags)
  lthreading=maxthreads.ge.1
  maxthreads=max(1,maxthreads)

  call sgran()                  ! see C rand generator (used in gran)

! Allocate data arrays
  pa=fftwf_alloc_complex(int(nfft,C_SIZE_T))
  call c_f_pointer(pa,a,[nfft])
  call c_f_pointer(pa,ar,[nfft])

  pb=fftwf_alloc_complex(int(nfft,C_SIZE_T))
  call c_f_pointer(pb,b,[nfft])
  call c_f_pointer(pb,br,[nfft])

  pc=fftwf_alloc_complex(int(nfft,C_SIZE_T))
  call c_f_pointer(pc,c,[nfft])

! Initialize FFTW threading
  if(lthreading) iret=fftwf_init_threads()

! Import FFTW wisdom, if available
  iret=fftwf_import_wisdom_from_filename(C_CHAR_'wis.dat' // C_NULL_CHAR)

  do i=1,nfft                           !Generate random data
     x=gran()
     y=gran()
     b(i)=cmplx(x,y)
  enddo
  iters=100

  write(*,1000)
1000 format(/'Problem  Threads Plan    Time    Gflops     RMS   iters'/    &
             '--------------------------------------------------------')

! Try nthreads = 1,maxthreads
  do nthreads=1,maxthreads
     a(1:nfft)=b(1:nfft)                             !Copy test data into a()
     call system_clock(count0,clkfreq)
! Make the plans
     if(lthreading) call fftwf_plan_with_nthreads(nthreads)
     if(lcomplex) then
        if(linplace) then
           plan1=fftwf_plan_dft_1d(nfft,a,a,-1,nflags)
           plan2=fftwf_plan_dft_1d(nfft,a,a,+1,nflags)
        else
           plan1=fftwf_plan_dft_1d(nfft,a,c,-1,nflags)
           plan2=fftwf_plan_dft_1d(nfft,c,a,+1,nflags)
        endif
     else
        if(linplace) then
           plan1=fftwf_plan_dft_r2c_1d(nfft,ar,a,nflags)
           plan2=fftwf_plan_dft_c2r_1d(nfft,a,ar,nflags)
        else
           plan1=fftwf_plan_dft_r2c_1d(nfft,ar,c,nflags)
           plan2=fftwf_plan_dft_c2r_1d(nfft,c,ar,nflags)
        endif
     endif
     call system_clock(count1,clkfreq)
     tplan=0.5*float(count1-count0)/float(clkfreq)    !Plan time for one transform

     total=0.
     do iter=1,iters                             !Do many iterations
        a=b                                      !Copy test data into a()
        call system_clock(count0,clkfreq)
! Compute the transforms
        if(lcomplex) then
           if(linplace) then
              call fftwf_execute_dft(plan1,a,a)
              call fftwf_execute_dft(plan2,a,a)
           else
              call fftwf_execute_dft(plan1,a,c)
              call fftwf_execute_dft(plan2,c,a)
           endif
        else
           if(linplace) then
              call fftwf_execute_dft_r2c(plan1,ar,a)
              call fftwf_execute_dft_c2r(plan2,a,ar)
           else
              call fftwf_execute_dft_r2c(plan1,ar,c)
              call fftwf_execute_dft_c2r(plan2,c,ar)
           endif
        endif
        call system_clock(count1,clkfreq)
        total=total + float(count1-count0)/float(clkfreq)
        if(total>=1.0 .and. iter>=10) go to 40     !Cut iterations short ?
     enddo
     iter=iters

40   time=0.5*total/iter                         !Time for one FFT
     gflops=5.0/(1.e9*time/(nfft*log(float(nfft))/log(2.0)))
     a(1:nfft)=a(1:nfft)/nfft              !Normalize the back-transformed data

! Compute RMS difference between original data and back-transformed data.
     sq=0.
     if(lcomplex) then
        do i=1,nfft
           sq=sq + real(a(i)-b(i))**2 + aimag(a(i)-b(i))**2
        enddo
     else
        do i=1,nfft
           sq=sq + (ar(i)-br(i))**2
        enddo
     endif
     rms=sqrt(sq/nfft)

! Display results
     write(*,1050) problem,nthreads,tplan,time,gflops,rms,iter
1050 format(a9,i4,f8.3,f10.6,f7.2,f11.7,i5)
  enddo

! Export accumulated FFTW wisdom
  iret=fftwf_export_wisdom_to_filename(C_CHAR_'wis.dat' // C_NULL_CHAR)

! Clean up
  call fftwf_destroy_plan(plan1)
  call fftwf_destroy_plan(plan2)
  call fftwf_free(pa)
  call fftwf_free(pb)
  call fftwf_free(pc)
  call fftwf_cleanup_threads()
  call fftwf_cleanup()

end program timefft


code:
subroutine timefft_opts(npatience,nthreads,linplace,lcomplex,nfft,    &
     problem,nflags)

  use FFTW3

  logical linplace,lcomplex
  character problem*9,arg*12

  nargs=iargc()
  if(nargs.lt.3) then
     print*,'Usage: timefft npatience maxthreads [[o|i][r|c]]nfft'
     print*,'       npatience  - 0 to 4'
     print*,'       maxthreads - suggest #CPUs or #CPUs-1'
     print*,'       o,i        - out-of-place or in-place (default=in-place)'
     print*,'       r,c        - real or complex (default=complex)'
     print*,' '
     print*,'Examples:'
     print*,'       timefft 1 1 32768     (1 thread, in-place, complex)'
     print*,'       timefft 2 3 or32768   (more patient, 3 threads,'
     print*,'                              out-of-place, real)'
     stop
  endif

  call getarg(1,arg)
  read(arg,*) npatience
  call getarg(2,arg)
  read(arg,*) nthreads
  call getarg(3,arg)
  linplace=arg(1:1).ne.'o' .and. arg(2:2).ne.'o'
  lcomplex=arg(1:1).ne.'r' .and. arg(2:2).ne.'r'
  k=3
  if(ichar(arg(2:2)).ge.48 .and. ichar(arg(2:2)).le.57) k=2
  if(ichar(arg(1:1)).ge.48 .and. ichar(arg(1:1)).le.57) k=1
  read(arg(k:),*) nfft

  write(problem,'(i9)') nfft
  problem='ic'//adjustl(problem)
  if(.not.linplace) problem(1:1)='o'
  if(.not.lcomplex) problem(2:2)='r'

  nflags=FFTW_ESTIMATE
  if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
  if(npatience.eq.2) nflags=FFTW_MEASURE
  if(npatience.eq.3) nflags=FFTW_PATIENT
  if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE

  return
end subroutine timefft_opts


code:
module timer_c_wrapper
  use :: iso_c_binding, only: c_ptr
  use timer_module, only: timer, null_timer
  implicit none

  !
  ! C interoperable callback setup
  !
  abstract interface
     subroutine c_timer_callback (context, dname, k)
       use, intrinsic :: iso_c_binding, only: c_ptr, c_char
       implicit none
       type(c_ptr), value, intent(in) :: context
       character(c_char), intent(in) :: dname(*)
       integer, intent(in), value :: k
     end subroutine c_timer_callback
  end interface

  public :: init, fini

  private

  !
  ! the following are singleton items which assumes that any timer
  ! implementation should only assume one global instance, probably a
  ! struct or class object whose address is stored the context below
  !
  type(c_ptr), private :: the_context
  procedure(C_timer_callback), pointer, private :: the_callback

contains
  subroutine timer_callback_wrapper (dname, k)
    use, intrinsic :: iso_c_binding, only: c_null_char
    implicit none
    character(len=8), intent(in) :: dname
    integer, intent(in) :: k
    call the_callback (the_context, trim (dname) // c_null_char, k)
  end subroutine timer_callback_wrapper

  subroutine init (context, callback)
    use, intrinsic :: iso_c_binding, only: c_ptr, c_funptr, c_f_procpointer
    use iso_c_utilities, only: c_to_f_string
    use timer_module, only: timer
    implicit none
    type(c_ptr), value, intent(in) :: context
    type(c_funptr), value, intent(in) :: callback
    the_context=context
    call c_f_procpointer (callback, the_callback)
    timer => timer_callback_wrapper
  end subroutine init

  subroutine fini ()
    implicit none
    timer => null_timer
  end subroutine fini

end module timer_c_wrapper


code:
  integer :: level, onlevel(0:10)
  common/timer_private/ level, onlevel
  !$omp threadprivate(/timer_private/)


code:
module timer_impl
  !$ use omp_lib
  use :: iso_c_binding, only: c_ptr
  use timer_module, only: timer_callback
  implicit none

  public :: init_timer, fini_timer
  integer, public :: limtrace=0

  private

  integer, parameter :: MAXCALL=100
  integer :: lu=6
  real :: dut
  integer :: i,nmax=0,ncall(MAXCALL),nlevel(MAXCALL),nparent(MAXCALL)
  character(len=8) :: name(MAXCALL),space='        '
  logical :: on(MAXCALL)
  real :: total,sum,sumf,ut(MAXCALL),ut0(MAXCALL)
  !$ integer :: j,l,m,ntid(MAXCALL)

  !
  ! C interoperable callback setup
  !
  public :: C_init_timer
  abstract interface
     subroutine C_timer_callback (context, dname, k)
       use, intrinsic :: iso_c_binding, only: c_ptr
       implicit none
       type(c_ptr), intent(in) :: context
       character(len=8), intent(in) :: dname
       integer, intent(in) :: k
     end subroutine C_timer_callback
  end interface
  type(c_ptr), private :: the_context
  procedure(C_timer_callback), pointer, private :: the_C_callback

contains
  subroutine timer_callback_wrapper (dname, k)
    implicit none
    character(len=8), intent(in) :: dname
    integer, intent(in) :: k
    call the_C_callback (the_context, dname, k)
  end subroutine timer_callback_wrapper

  subroutine C_init_timer (context, callback) bind(C)
    use, intrinsic :: iso_c_binding, only: c_ptr, c_funptr, c_f_procpointer
    use iso_c_utilities, only: c_to_f_string
    use timer_module, only: timer
    implicit none
    type(c_ptr), intent(in) :: context
    type(c_funptr), intent(in) :: callback
    the_context=context
    call c_f_procpointer (callback, the_C_callback)
    timer => timer_callback_wrapper
  end subroutine C_init_timer

  !
  ! default Fortran implementation which is thread safe using OpenMP
  !
  subroutine default_timer (dname, k)

    ! Times procedure number n between a call with k=0 (tstart) and with
    ! k=1 (tstop). Accumulates sums of these times in array ut (user time).
    ! Also traces all calls (for debugging purposes) if limtrace.gt.0
    !
    ! If this is used with OpenMP than the /timer_private/ common
    ! block must be copyed into each thread of a thread team by using
    ! the copyin() clause on the !$omp parallel directive that creates
    ! the team.

    implicit none

    character(len=8), intent(in) :: dname
    integer, intent(in) :: k

    real :: ut1,eps=0.000001
    integer :: n,ndiv,ntrace=0
    !$ integer :: tid
    character(len=8) :: tname
    include 'timer_common.inc'

    !$omp critical(timer)
    if(limtrace.lt.0) go to 999
    if(k.gt.1) go to 40                        !Check for "all done" (k>1)
    onlevel(0)=0

    !$ tid=omp_get_thread_num()
    do n=1,nmax                                !Check for existing name/parent[/thread]
       if(name(n).eq.dname &
                                !$ .and.ntid(n).eq.tid &
            ) then
          if (on(n)) then
             if (nparent(n).eq.onlevel(level-1)) goto 20
          else
             if (nparent(n).eq.onlevel(level)) goto 20
          end if
       end if
    enddo

    nmax=nmax+1                                !This is a new one
    n=nmax
    !$ ntid(n)=tid
    ncall(n)=0
    on(n)=.false.
    ut(n)=eps
    name(n)=dname

20  if(k.eq.0) then                                !Get start times (k=0)
       if(on(n)) then
          print*,'Error in timer: ',dname,' already on.'
       end if
       level=level+1 !Increment the level
       on(n)=.true.
       !     call system_clock(icount,irate)
       !     ut0(n)=float(icount)/irate
       !     call cpu_time(ut0(n))
       ut0(n)=secnds(0.0)

       ncall(n)=ncall(n)+1
       if(ncall(n).gt.1.and.nlevel(n).ne.level) then
          !recursion is happening
          !
          !TODO: somehow need to account for this deeper call at the
          !shallowest instance in the call chain and this needs to be
          !done without incrementing anything here other than counters
          !and timers
          !
          nlevel(n)=-1
       else
          nlevel(n)=level
       endif
       nparent(n)=onlevel(level-1)
       onlevel(level)=n

    else if(k.eq.1) then        !Get stop times and accumulate sums. (k=1)
       if(on(n)) then
          on(n)=.false.
          !        call system_clock(icount,irate)
          !        ut1=float(icount)/irate
          !        call cpu_time(ut1)
          ut1=secnds(0.0)

          ut(n)=ut(n)+ut1-ut0(n)
       endif
       level=level-1
    endif

    ntrace=ntrace+1
    if(ntrace.lt.limtrace) write(lu,1020) ntrace,tname,k,level,nparent(n)
1020 format(i8,': ',a8,3i5)
    go to 998

    ! Write out the timer statistics

40  write(lu,1040)
1040 format(/' Name                 Time  Frac     dTime',       &
         ' dFrac    Calls'/58('-'))

    !$ !walk backwards through the database rolling up thread data by call chain
    !$ do i=nmax,1,-1
    !$    do j=1,i-1
    !$       l=j
    !$       m=i
    !$       do while (name(l).eq.name(m))
    !$          l=nparent(l)
    !$          m=nparent(m)
    !$          if (l.eq.0.or.m.eq.0) exit
    !$       end do
    !$       if (l.eq.0.and.m.eq.0) then
    !$          !same call chain so roll up data
    !$          ncall(j)=ncall(j)+ncall(i)
    !$          ut(j)=ut(j)+ut(i)
    !$          do n=1,nmax
    !$            if (nparent(n).eq.i) nparent(n)=j
    !$          end do
    !$          name(i)=space
    !$          exit
    !$       end if
    !$    end do
    !$ end do

    if(k.gt.100) then
       ndiv=k-100
       do i=1,nmax
          ncall(i)=ncall(i)/ndiv
          ut(i)=ut(i)/ndiv
       enddo
    endif

    total=ut(1)
    sum=0.
    sumf=0.
    call print_root(1)
    write(lu,1070) sum,sumf
1070 format(58('-')/32x,f10.3,f6.2)
    nmax=0
    eps=0.000001
    ntrace=0
    level=0
    onlevel(0)=0

998 flush(lu)

999 continue

    !$omp end critical(timer)
    return
  end subroutine default_timer

  recursive subroutine print_root(i)
    implicit none
    integer, intent(in) :: i
    character(len=16) :: sname
    real :: dutf, utf
    integer :: j, kk

    if (i.le.nmax) then
       if (name(i).ne.space) then
          dut=ut(i)
          do j=i,nmax
             if (name(j).ne.space.and.nparent(j).eq.i) dut=dut-ut(j)
          enddo
          if(dut.lt.0.0) dut=0.0
          utf=ut(i)/total
          dutf=dut/total
          sum=sum+dut
          sumf=sumf+dutf
          kk=nlevel(i)
          sname=space(1:kk)//name(i)//space(1:8-kk)
          write(lu,2000) sname,ut(i),utf,dut,dutf,ncall(i)
2000      format(a16,2(f10.3,f6.2),i9)
          do j=i,nmax
             if(nparent(j).eq.i) call print_root(j)
          enddo
       end if
    end if
    return
  end subroutine print_root

  subroutine init_timer (filename)
    use, intrinsic :: iso_c_binding, only: c_char
    use timer_module, only: timer
    implicit none
    character(len=*), optional, intent(in) :: filename
    include 'timer_common.inc'
    data level/0/, onlevel/11 * 0/
    if (present (filename)) then
       open (newunit=lu, file=filename, status='unknown')
    else
       open (newunit=lu, file='timer.out', status='unknown')
    end if
    timer => default_timer
  end subroutine init_timer

  subroutine fini_timer ()
    use timer_module, only: timer, null_timer
    implicit none
    timer => null_timer
    close (lu)
  end subroutine fini_timer

end module timer_impl


code:
module timer_module
  implicit none

  abstract interface
     subroutine timer_callback (dname, k)
       character(len=8), intent(in) :: dname
       integer, intent(in) :: k
     end subroutine timer_callback
  end interface

  public :: null_timer
  procedure(timer_callback), pointer :: timer => null_timer

contains
  !
  ! default Fortran implementation which does nothing
  !
  subroutine null_timer (dname, k)
    implicit none
    character(len=8), intent(in) :: dname
    integer, intent(in) :: k
    if(dname.eq.'99999999' .and. k.eq.9999) stop  !Silence compiler warnings
  end subroutine null_timer
end module timer_module


code:
subroutine timf2(x0,k,nfft,nwindow,nb,peaklimit,x1,     &
     slimit,lstrong,px,nzap)

! Sequential processing of time-domain I/Q data, using Linrad-like
! "first FFT" and "first backward FFT", treating frequencies with
! strong signals differently.  Noise blanking is applied to weak
! signals only.

!  x0       - real input data
!  nfft     - length of FFTs
!  nwindow  - 0 for no window, 2 for sin^2 window
!  x1       - real output data

! Non-windowed processing means no overlap, so kstep=nfft.
! Sin^2 window has 50% overlap, kstep=nfft/2.

! Frequencies with strong signals are identified and separated.  Back
! transforms are done separately for weak and strong signals, so that
! noise blanking can be applied to the weak-signal portion.  Strong and
! weak are finally re-combined, in the time domain.

  parameter (MAXFFT=1024,MAXNH=MAXFFT/2)
  parameter (MAXSIGS=100)
  real x0(0:nfft-1),x1(0:nfft-1)
  real x(0:MAXFFT-1),xw(0:MAXFFT-1),xs(0:MAXFFT-1)
  real xwov(0:MAXNH-1),xsov(0:MAXNH-1)
  complex cx(0:MAXFFT-1),cxt(0:MAXFFT-1)
  complex cxs(0:MAXFFT-1)                     !Strong signals
  complex cxw(0:MAXFFT-1)                     !Weak signals
  real*4 w(0:MAXFFT-1)
  real*4 s(0:MAXNH)
  logical*1 lstrong(0:MAXNH),lprev
  integer ia(MAXSIGS),ib(MAXSIGS)
  logical first
  equivalence (x,cx),(xw,cxw),(xs,cxs)
  data first/.true./
  data k0/99999999/
  save

  if(first) then
     pi=4.0*atan(1.0)
     do i=0,nfft-1
        w(i)=(sin(i*pi/nfft))**2
     enddo
     s=0.
     nh=nfft/2
     kstep=nfft
     if(nwindow.eq.2) kstep=nh
     fac=1.0/nfft
     slimit=1.e30
     first=.false.
  endif

  if(k.lt.k0) then
     xsov=0.
     xwov=0.
  endif
  k0=k

  x(0:nfft-1)=x0
  if(nwindow.eq.2) x(0:nfft-1)=w(0:nfft-1)*x(0:nfft-1)
  call four2a(x,nfft,1,-1,0)                       !First forward FFT, r2c
  cxt(0:nh)=cx(0:nh)

! Identify frequencies with strong signals.
  do i=0,nh
     p=real(cxt(i))**2 + aimag(cxt(i))**2
     s(i)=p
  enddo
  ave=sum(s(0:nh))/nh
  lstrong(0:nh)=s(0:nh).gt.10.0*ave

  nsigs=0
  lprev=.false.
  iwid=1
  ib=-99
  do i=0,nh
     if(lstrong(i) .and. (.not.lprev)) then
        if(nsigs.lt.MAXSIGS) nsigs=nsigs+1
        ia(nsigs)=i-iwid
        if(ia(nsigs).lt.0) ia(nsigs)=0
     endif
     if(.not.lstrong(i) .and. lprev) then
        ib(nsigs)=i-1+iwid
        if(ib(nsigs).gt.nh) ib(nsigs)=nh
     endif
     lprev=lstrong(i)
  enddo

  if(nsigs.gt.0) then
     do i=1,nsigs
        ja=ia(i)
        jb=ib(i)
        if(ja.lt.0 .or. ja.gt.nh .or. jb.lt.0 .or. jb.gt.nh) then
           cycle
        endif
        if(jb.eq.-99) jb=ja + min(2*iwid,nh)
        lstrong(ja:jb)=.true.
     enddo
  endif

! Copy frequency-domain data into array cs (strong) or cw (weak).
  do i=0,nh
     if(lstrong(i)) then
        cxs(i)=fac*cxt(i)
        cxw(i)=0.
     else
        cxw(i)=fac*cxt(i)
        cxs(i)=0.
     endif
  enddo

  call four2a(cxw,nfft,1,1,-1)           !Transform weak and strong back
  call four2a(cxs,nfft,1,1,-1)           !to time domain, separately (c2r)

  if(nwindow.eq.2) then
     xw(0:nh-1)=xw(0:nh-1)+xwov(0:nh-1)     !Add previous segment's 2nd half
     xwov(0:nh-1)=xw(nh:nfft-1)             !Save 2nd half
     xs(0:nh-1)=xs(0:nh-1)+xsov(0:nh-1)     !Ditto for strong signals
     xsov(0:nh-1)=xs(nh:nfft-1)
  endif

! Apply noise blanking to weak data
  if(nb.ne.0) then
     do i=0,kstep-1
        peak=abs(xw(i))
        if(peak.gt.peaklimit) then
           xw(i)=0.
           nzap=nzap+1
        endif
     enddo
  endif

! Compute power levels from weak data only
  do i=0,kstep-1
     px=px + xw(i)**2
  enddo

  x1(0:kstep-1)=xw(0:kstep-1) + xs(0:kstep-1)     !Recombine weak + strong

  return
end subroutine timf2


code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define RADS 0.0174532925199433
#define DEGS 57.2957795130823
#define TPI 6.28318530717959
#define PI 3.1415927

/* ratio of     earth radius to astronomical unit */
#define ER_OVER_AU 0.0000426352325194252

/* all prototypes here */

double getcoord(int coord);
void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
double range(double y);
double rangerad(double y);
double days(int y, int m, int dn, double hour);
double days_(int *y, int *m, int *dn, double *hour);
void moonpos(double, double *, double *, double *);
void sunpos(double , double *, double *, double *);
double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
double atan22(double y, double x);
double epsilon(double d);
void equatorial(double d, double *lon, double *lat, double *r);
void ecliptic(double d, double *lon, double *lat, double *r);
double gst(double d);
void topo(double lst, double glat, double *alp, double *dec, double *r);
double alt(double glat, double ha, double dec);
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
int daysinmonth(int y, int m);
int isleap(int y);
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
   double *mrv, double *l, double *b, double *paxis);

static const char
usage[] = "  Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
            "example: tmoon 200009 0 -00155 5230\n";

/*
  getargs() gets the arguments from the command line, does some basic error
  checking, and converts arguments into numerical form. Arguments are passed
  back in pointers. Error messages print to stderr so re-direction of output
  to file won't leave users blind. Error checking prints list of all errors
  in a command line before quitting.
*/
void getargs(int argc, char *argv[], int *y, int *m, double *tz,
             double *glong, double *glat) {

  int date, latitude, longitude;
  int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
  int longminflag = 0, latminflag = 0, dflag = 0;

  /* if not right number of arguments, then print example command line */

  if (argc !=5) {
    fprintf(stderr, usage);
    exit(EXIT_FAILURE);
  }

  date = atoi(argv[1]);
  *y = date / 100;
  *m = date - *y * 100;
  *tz = (double) atof(argv[2]);
  longitude = atoi(argv[3]);
  latitude = atoi(argv[4]);
  *glong = RADS * getcoord(longitude);
  *glat = RADS * getcoord(latitude);

  /* set a flag for each error found */

  if (*m > 12 || *m < 1) mflag = 1;
  if (*y > 2500) yflag = 1;
  if (date < 150001) dflag = 1;
  if (fabs((float) *glong) > 180 * RADS) longflag = 1;
  if (abs(longitude) % 100 > 59) longminflag = 1;
  if (fabs((float) *glat) > 90 * RADS) latflag = 1;
  if (abs(latitude) % 100 > 59) latminflag = 1;
  if (fabs((float) *tz) > 12) tzflag = 1;

  /* print all the errors found */

  if (dflag == 1) {
    fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
  }
  if (yflag == 1) {
    fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
  }
  if (mflag == 1) {
    fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
  }
  if (tzflag == 1) {
    fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
  }
  if (longflag == 1) {
    fprintf(stderr, "long: must be in range +/- 180 degrees\n");
  }
  if (longminflag == 1) {
    fprintf(stderr, "long: last two digits are arcmin - max 59\n");
  }
  if (latflag == 1) {
    fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
  }
  if (latminflag == 1) {
    fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
  }

  /* quits if one or more flags set */

  if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
    exit(EXIT_FAILURE);
  }

}

/*
   returns coordinates in decimal degrees given the
   coord as a ddmm value stored in an integer.
*/
double getcoord(int coord) {
  int west = 1;
  double glg, deg;
  if (coord < 0) west = -1;
  glg = fabs((double) coord/100);
  deg = floor(glg);
  glg = west* (deg + (glg - deg)*100 / 60);
  return(glg);
}

/*
  days() takes the year, month, day in the month and decimal hours
  in the day and returns the number of days since J2000.0.
  Assumes Gregorian calendar.
*/
double days(int y, int m, int d, double h) {
  int a, b;
  double day;

  /*
    The lines below work from 1900 march to feb 2100
    a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
    day = (double)a - 730531.5 + hour / 24;
  */

  /*  These lines work for any Gregorian date since 0 AD */
  if (m ==1 || m==2) {
    m +=12;
    y -= 1;
  }
  a = y / 100;
  b = 2 - a + a/4;
  day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
    + d + b - 1524.5 - 2451545 + h/24;
  return(day);
}
double days_(int *y0, int *m0, int *d0, double *h0)
{
  return days(*y0,*m0,*d0,*h0);
}

/*
Returns 1 if y a leap year, and 0 otherwise, according
to the Gregorian calendar
*/
int isleap(int y) {
  int a = 0;
  if(y % 4 == 0) a = 1;
  if(y % 100 == 0) a = 0;
  if(y % 400 == 0) a = 1;
  return(a);
}

/*
Given the year and the month, function returns the
number of days in the month. Valid for Gregorian
calendar.
*/
int daysinmonth(int y, int m) {
  int b = 31;
  if(m == 2) {
    if(isleap(y) == 1) b= 29;
    else b = 28;
  }
  if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
  return(b);
}

/*
moonpos() takes days from J2000.0 and returns ecliptic coordinates
of moon in the pointers. Note call by reference.
This function is within a couple of arcminutes most of the time,
and is truncated from the Meeus Ch45 series, themselves truncations of
ELP-2000. Returns moon distance in earth radii.
Terms have been written out explicitly rather than using the
table based method as only a small number of terms is
retained.
*/
void moonpos(double d, double *lambda, double *beta, double *rvec) {
  double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;

  t = d / 36525;

  L = range(218.3164591  + 481267.88134236  * t) * RADS;
  D = range(297.8502042  + 445267.1115168  * t) * RADS;
  M = range(357.5291092  + 35999.0502909  * t) * RADS;
  M1 = range(134.9634114  + 477198.8676313  * t - .008997 * t * t) * RADS;
  F = range(93.27209929999999  + 483202.0175273  * t - .0034029*t*t)*RADS;
  e = 1 - .002516 * t;

  dl =      6288774 * sin(M1);
  dl +=     1274027 * sin(2 * D - M1);
  dl +=      658314 * sin(2 * D);
  dl +=      213618 * sin(2 * M1);
  dl -=  e * 185116 * sin(M);
  dl -=      114332 * sin(2 * F) ;
  dl +=       58793 * sin(2 * D - 2 * M1);
  dl +=   e * 57066 * sin(2 * D - M - M1) ;
  dl +=       53322 * sin(2 * D + M1);
  dl +=   e * 45758 * sin(2 * D - M);
  dl -=   e * 40923 * sin(M - M1);
  dl -=       34720 * sin(D) ;
  dl -=   e * 30383 * sin(M + M1) ;
  dl +=       15327 * sin(2 * D - 2 * F) ;
  dl -=       12528 * sin(M1 + 2 * F);
  dl +=       10980 * sin(M1 - 2 * F);
  lm = rangerad(L + dl / 1000000 * RADS);

  dB =   5128122 * sin(F);
  dB +=   280602 * sin(M1 + F);
  dB +=   277693 * sin(M1 - F);
  dB +=   173237 * sin(2 * D - F);
  dB +=    55413 * sin(2 * D - M1 + F);
  dB +=    46271 * sin(2 * D - M1 - F);
  dB +=    32573 * sin(2 * D + F);
  dB +=    17198 * sin(2 * M1 + F);
  dB +=     9266 * sin(2 * D + M1 - F);
  dB +=     8822 * sin(2 * M1 - F);
  dB += e * 8216 * sin(2 * D - M - F);
  dB +=     4324 * sin(2 * D - 2 * M1 - F);
  bm = dB / 1000000 * RADS;

  dR =    -20905355 * cos(M1);
  dR -=     3699111 * cos(2 * D - M1);
  dR -=     2955968 * cos(2 * D);
  dR -=      569925 * cos(2 * M1);
  dR +=   e * 48888 * cos(M);
  dR -=        3149 * cos(2 * F);
  dR +=      246158 * cos(2 * D - 2 * M1);
  dR -=  e * 152138 * cos(2 * D - M - M1) ;
  dR -=      170733 * cos(2 * D + M1);
  dR -=  e * 204586 * cos(2 * D - M);
  dR -=  e * 129620 * cos(M - M1);
  dR +=      108743 * cos(D);
  dR +=  e * 104755 * cos(M + M1);
  dR +=       79661 * cos(M1 - 2 * F);
  rm = 385000.56  + dR / 1000;

  *lambda = lm;
  *beta = bm;
  /* distance to Moon must be in Earth radii */
  *rvec = rm / 6378.14;
}

/*
topomoon() takes the local siderial time, the geographical
latitude of the observer, and pointers to the geocentric
equatorial coordinates. The function overwrites the geocentric
coordinates with topocentric coordinates on a simple spherical
earth model (no polar flattening). Expects Moon-Earth distance in
Earth radii.    Formulas scavenged from Astronomical Almanac 'low
precision formulae for Moon position' page D46.
*/

void topo(double lst, double glat, double *alp, double *dec, double *r) {
  double x, y, z, r1;
  x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
  y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
  z = *r * sin(*dec)  - sin(glat);
  r1 = sqrt(x*x + y*y + z*z);
  *alp = atan22(y, x);
  *dec = asin(z / r1);
  *r = r1;
}

/*
moontransit() takes date, the time zone and geographic longitude
of observer and returns the time (decimal hours) of lunar transit
on that day if there is one, and sets the notransit flag if there
isn't. See Explanatory Supplement to Astronomical Almanac
section 9.32 and 9.31 for the method.
*/

double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
  double hm, ht, ht1, lon, lat, rv, dnew, lst;
  int itcount;

  ht1 = 180 * RADS;
  ht = 0;
  itcount = 0;
  *notransit = 0;
  do {
    ht = ht1;
    itcount++;
    dnew = days(y, m, d, ht * DEGS/15) - tz/24;
    lst = gst(dnew) + glong;
    /* find the topocentric Moon ra (hence hour angle) and dec */
    moonpos(dnew, &lon, &lat, &rv);
    equatorial(dnew, &lon, &lat, &rv);
    topo(lst, glat, &lon, &lat, &rv);
    hm = rangerad(lst -  lon);
    ht1 = rangerad(ht - hm);
    /* if no convergence, then no transit on that day */
    if (itcount > 30) {
      *notransit = 1;
      break;
    }
  }
  while (fabs(ht - ht1) > 0.04 * RADS);
  return(ht1);
}

/*
  Calculates the selenographic coordinates of either the sub Earth point
  (optical libration) or the sub-solar point (selen. coords of centre of
  bright hemisphere).  Based on Meeus chapter 51 but neglects physical
  libration and nutation, with some simplification of the formulas.
*/
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
  double i, f, omega, w, y, x, a, t, eps;
  t = day / 36525;
  i = 1.54242 * RADS;
  eps = epsilon(day);
  f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
  omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
  w = lambda - omega;
  y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
  x = cos(w) * cos(beta);
  a = atan22(y, x);
  *l = a - f;

  /*  kludge to catch cases of 'round the back' angles  */
  if (*l < -90 * RADS) *l += TPI;
  if (*l > 90 * RADS)  *l -= TPI;
  *b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));

  /*  pa pole axis - not used for Sun stuff */
  x = sin(i) * sin(omega);
  y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
  w = atan22(x, y);
  *p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
}

/*
  Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
  eq coords Sun
  Returns: position angle of bright limb wrt NCP, percentage illumination
  of Sun
*/
void illumination(double day , double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
  double x, y, phi, i;
  (void)day;
  y = cos(sdec) * sin(sra - lra);
  x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
  *pabl = atan22(y, x);
  phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
  i = atan22(sin(phi) , (dr - cos(phi)));
  *ill = 0.5*(1 + cos(i));
}

/*
sunpos() takes days from J2000.0 and returns ecliptic longitude
of Sun in the pointers. Latitude is zero at this level of precision,
but pointer left in for consistency in number of arguments.
This function is within 0.01 degree (1 arcmin) almost all the time
for a century either side of J2000.0. This is from the 'low precision
fomulas for the Sun' from C24 of Astronomical Alamanac
*/
void sunpos(double d, double *lambda, double *beta, double *rvec) {
  double L, g, ls, bs, rs;

  L = range(280.461 + .9856474 * d) * RADS;
  g = range(357.528 + .9856003 * d) * RADS;
  ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
  bs = 0;
  rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
  *lambda = ls;
  *beta = bs;
  *rvec = rs;
}

/*
this routine returns the altitude given the days since J2000.0
the hour angle and declination of the object and the latitude
of the observer. Used to find the Sun's altitude to put a letter
code on the transit time, and to find the Moon's altitude at
transit just to make sure that the Moon is visible.
*/
double alt(double glat, double ha, double dec) {
  return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
}

/* returns an angle in degrees in the range 0 to 360 */
double range(double x) {
  double a, b;
  b = x / 360;
  a = 360 * (b - floor(b));
  if (a < 0)
    a = 360 + a;
  return(a);
}

/* returns an angle in rads in the range 0 to two pi */
double rangerad(double x) {
  double a, b;
  b = x / TPI;
  a = TPI * (b - floor(b));
  if (a < 0)
    a = TPI + a;
  return(a);
}

/*
gets the atan2 function returning angles in the right
order and  range
*/
double atan22(double y, double x) {
  double a;

  a = atan2(y, x);
  if (a < 0) a += TPI;
  return(a);
}

/*
returns mean obliquity of ecliptic in radians given days since
J2000.0.
*/
double epsilon(double d) {
  double t = d/ 36525;
  return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
}

/*
replaces ecliptic coordinates with equatorial coordinates
note: call by reference destroys original values
R is unchanged.
*/
void equatorial(double d, double *lon, double *lat, double * r) {
  double  eps, ceps, seps, l, b;
  (void)r;

  l = *lon;
  b = * lat;
  eps = epsilon(d);
  ceps = cos(eps);
  seps = sin(eps);
  *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
  *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
}

/*
replaces equatorial coordinates with ecliptic ones. Inverse
of above, but used to find topocentric ecliptic coords.
*/
void ecliptic(double d, double *lon, double *lat, double * r) {
  double  eps, ceps, seps, alp, dec;
  (void)r;

  alp = *lon;
  dec = *lat;
  eps = epsilon(d);
  ceps = cos(eps);
  seps = sin(eps);
  *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
  *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
}

/*
returns the siderial time at greenwich meridian as
an angle in radians given the days since J2000.0
*/
double gst( double d) {
  double t = d / 36525;
  double theta;
  theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
  return(theta * RADS);
}

void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
   double *mrv, double *l, double *b, double *paxis)
{
  double mlambda, mbeta;
  double malpha, mdelta;
  double lst, mhr;
  double tlambda, tbeta, trv;

  lst = gst(*day) + *glong;

  /* find Moon topocentric coordinates for libration calculations */

  moonpos(*day, &mlambda, &mbeta, mrv);
  malpha = mlambda;
  mdelta = mbeta;
  equatorial(*day, &malpha, &mdelta, mrv);
  topo(lst, *glat, &malpha, &mdelta, mrv);
  mhr = rangerad(lst - malpha);
  *moonalt = alt(*glat, mhr, mdelta);

  /* Optical libration and Position angle of the Pole */

  tlambda = malpha;
  tbeta = mdelta;
  trv = *mrv;
  ecliptic(*day, &tlambda, &tbeta, &trv);
  libration(*day, tlambda, tbeta, malpha,  l, b, paxis);
}