code:
module options
  !
  ! Source code copied from:
  ! http://fortranwiki.org/fortran/show/Command-line+arguments
  !
  implicit none

  type option
     !> Long name.
     character(len=100) :: name
     !> Does the option require an argument?
     logical :: has_arg
     !> Corresponding short name.
     character :: chr
     !> Description.
     character(len=500) :: descr
     !> Argument name, if required.
     character(len=20) :: argname
   contains
     procedure :: print => print_opt
  end type option

contains

  !> Parse command line options. Options and their arguments must come before
  !> all non-option arguments. Short options have the form "-X", long options
  !> have the form "--XXXX..." where "X" is any character. Parsing can be
  !> stopped with the option '--'.
  !> The following code snippet illustrates the intended use:
  !> \code
  !> do
  !>   call getopt (..., optchar=c, ...)
  !>   if (stat /= 0) then
  !>     ! optional error handling
  !>     exit
  !>   end if
  !>   select case (c)
  !>     ! process options
  !>   end select
  !> end do
  !> \endcode
  subroutine getopt (options, longopts, optchar, optarg, arglen, stat, &
       offset, remain, err)
    use iso_fortran_env, only: error_unit

    !> String containing the characters that are valid short options. If
    !> present, command line arguments are scanned for those options.
    !> If a character is followed by a colon (:) its corresponding option
    !> requires an argument. E.g. "vn:" defines two options -v and -n with -n
    !> requiring an argument.
    character(len=*), intent(in), optional :: options

    !> Array of long options. If present, options of the form '--XXXX...' are
    !> recognised. Each option has an associated option character. This can be
    !> any character of default kind, it is just an identifier. It can, but
    !> doesn't have to, match any character in the options argument. In fact it
    !> is possible to only pass long options and no short options at all.
    !> Only name, has_arg and chr need to be set.
    type(option), intent(in), optional :: longopts(:)

    !> If stat is not 1, optchar contains the option character that was parsed.
    !> Otherwise its value is undefined.
    character, intent(out), optional :: optchar

    !> If stat is 0 and the parsed option requires an argument, optarg contains
    !> the first len(optarg) (but at most 500) characters of that argument.
    !> Otherwise its value is undefined. If the arguments length exceeds 500
    !> characters and err is .true., a warning is issued.
    character(len=*), intent(out), optional :: optarg

    !> If stat is 0 and the parsed option requires an argument, arglen contains
    !> the actual length of that argument. Otherwise its value is undefined.
    !> This can be used to make sure the argument was not truncated by the
    !> limited length of optarg.
    integer, intent(out), optional :: arglen

    !> Status indicator. Can have the following values:
    !>   -  0: An option was successfully parsed.
    !>   -  1: Parsing stopped successfully because a non-option or '--' was
    !>         encountered.
    !>   - -1: An unrecognised option was encountered.
    !>   - -2: A required argument was missing.
    !>   .
    !> Its value is never undefined.
    integer, intent(out), optional :: stat

    !> If stat is 1, offset contains the number of the argument before the
    !> first non-option argument, i.e. offset+n is the nth non-option argument.
    !> If stat is not 1, offset contains the number of the argument that would
    !> be parsed in the next call to getopt. This number can be greater than
    !> the actual number of arguments.
    integer, intent(out), optional :: offset

    !> If stat is 1, remain contains the number of remaining non-option
    !> arguments, i.e. the non-option arguments are in the range
    !> (offset+1:offset+remain). If stat is not 1, remain is undefined.
    integer, intent(out), optional :: remain

    !> If err is present and .true., getopt prints messages to the standard
    !> error unit if an error is encountered (i.e. whenever stat would be set
    !> to a negative value).
    logical, intent(in), optional :: err

    integer, save :: pos = 1, cnt = 0
    character(len=500), save :: arg

    integer :: chrpos, length, st, id = 0
    character :: chr
    logical :: long

    if (cnt == 0) cnt = command_argument_count()
    long = .false.

    ! no more arguments left
    if (pos > cnt) then
       pos = pos - 1
       st = 1
       goto 10
    end if

    call get_command_argument (pos, arg, length)

    ! is argument an option?
    if (arg(1:1) == '-') then

       chr = arg(2:2)

       ! too long ('-xxxx...') for one dash?
       if (chr /= '-' .and. len_trim(arg) > 2) then
          st = -1
          goto 10
       end if

       ! forced stop ('--')
       if (chr == '-' .and. arg(3:3) == ' ') then
          st = 1
          goto 10
       end if

       ! long option ('--xxx...')
       if (chr == '-') then

          long = .true.

          ! check if valid
          id = lookup(arg(3:))

          ! option is invalid, stop
          if (id == 0) then
             st = -1
             goto 10
          end if

          chr = longopts(id)%chr

          ! check if option requires an argument
          if (.not. longopts(id)%has_arg) then
             st = 0
             goto 10
          end if

          ! check if there are still arguments left
          if (pos == cnt) then
             st = -2
             goto 10
          end if

          ! go to next position
          pos = pos + 1

          ! get argument
          call get_command_argument (pos, arg, length)

          ! make sure it is not an option
          if (arg(1:1) == '-') then
             st = -2
             pos = pos - 1
             goto 10
          end if

       end if

       ! short option
       ! check if valid
       if (present(options)) then
          chrpos = scan(options, chr)
       else
          chrpos = 0
       end if

       ! option is invalid, stop
       if (chrpos == 0) then
          st = -1
          goto 10
       end if

       ! look for argument requirement
       if (chrpos < len_trim(options)) then
          if (options(chrpos+1:chrpos+1) == ':') then

             ! check if there are still arguments left
             if (pos == cnt) then
                st = -2
                goto 10
             end if

             ! go to next position
             pos = pos + 1

             ! get argument
             call get_command_argument (pos, arg, length)

             ! make sure it is not an option
             if (arg(1:1) == '-') then
                st = -2
                pos = pos - 1
                goto 10
             end if

          end if
       end if

       ! if we get to this point, no error happened
       ! return option and the argument (if there is one)
       st = 0
       goto 10
    end if

    ! not an option, parsing stops
    st = 1
    ! we are already at the first non-option argument
    ! go one step back to the last option or option argument
    pos = pos - 1

    ! error handling and setting of return values
10  continue

    if (present(err)) then
       if (err) then

          select case (st)
          case (-1)
             write (error_unit, *) "error: unrecognised option: " // trim(arg)
          case (-2)
             if (.not. long) then
                write (error_unit, *) "error: option -" // chr &
                     // " requires an argument"
             else
                write (error_unit, *) "error: option --" &
                     // trim(longopts(id)%name) // " requires an argument"
             end if
          end select

       end if
    end if

    if (present(optchar)) optchar = chr
    if (present(optarg))  optarg  = arg
    if (present(arglen))  arglen  = length
    if (present(stat))    stat    = st
    if (present(offset))  offset  = pos
    if (present(remain))  remain  = cnt-pos

    ! setup pos for next call to getopt
    pos = pos + 1

  contains

    integer function lookup (name)
      character(len=*), intent(in) :: name
      integer :: i

      ! if there are no long options, skip the loop
      if (.not. present(longopts)) goto 10

      do i = 1, size(longopts)
         if (name == longopts(i)%name) then
            lookup = i
            return
         end if
      end do
      ! if we get to this point, the option was not found

10    lookup = 0
    end function lookup

  end subroutine getopt

  !============================================================================

  !> Print an option in the style of a man page. I.e.
  !> \code
  !> -o arg
  !> --option arg
  !>    description.................................................................
  !>    ............................................................................
  !> \endcode
  subroutine print_opt (opt, unit)
    !> the option
    class(option), intent(in) :: opt
    !> logical unit number
    integer, intent(in) :: unit

    integer :: l, c1, c2

    if (opt%has_arg) then
       write (unit, '(1x,"-",a,1x,a)') opt%chr, trim(opt%argname)
       write (unit, '(1x,"--",a,1x,a)') trim(opt%name), trim(opt%argname)
    else
       write (unit, '(1x,"-",a)') opt%chr
       write (unit, '(1x,"--",a)') trim(opt%name)
    end if
    l = len_trim(opt%descr)

    ! c1 is the first character of the line
    ! c2 is one past the last character of the line
    c1 = 1
    do
       if (c1 > l) exit
       ! print at maximum 4+76 = 80 characters
       c2 = min(c1 + 76, 500)
       ! if not at the end of the whole string
       if (c2 /= 500) then
          ! find the end of a word
          do
             if (opt%descr(c2:c2) == ' ') exit
             c2 = c2-1
          end do
       end if
       write (unit, '(4x,a)') opt%descr(c1:c2-1)
       c1 = c2+1
    end do

  end subroutine print_opt

end module options


code:
module packjt

! These variables are accessible from outside via "use packjt":
  integer jt_itype,jt_nc1,jt_nc2,jt_ng,jt_k1,jt_k2
  character*6 jt_c1,jt_c2,jt_c3

  contains

subroutine packbits(dbits,nsymd,m0,sym)

 ! Pack 0s and 1s from dbits() into sym() with m0 bits per word.
 ! NB: nsymd is the number of packed output words.

   integer sym(:)
   integer*1 dbits(:)

   k=0
   do i=1,nsymd
      n=0
      do j=1,m0
         k=k+1
         m=dbits(k)
         n=ior(ishft(n,1),m)
      enddo
      sym(i)=n
   enddo

   return
 end subroutine packbits

 subroutine unpackbits(sym,nsymd,m0,dbits)

 ! Unpack bits from sym() into dbits(), one bit per byte.
 ! NB: nsymd is the number of input words, and m0 their length.
 ! there will be m0*nsymd output bytes, each 0 or 1.

   integer sym(:)
   integer*1 dbits(:)

   k=0
   do i=1,nsymd
      mask=ishft(1,m0-1)
      do j=1,m0
         k=k+1
         dbits(k)=0
         if(iand(mask,sym(i)).ne.0) dbits(k)=1
         mask=ishft(mask,-1)
      enddo
   enddo

   return
 end subroutine unpackbits

 subroutine packcall(callsign,ncall,text)

 ! Pack a valid callsign into a 28-bit integer.

   parameter (NBASE=37*36*10*27*27*27)
   character callsign*6,c*1,tmp*6
   logical text

   text=.false.

 ! Work-around for Swaziland prefix:
   if(callsign(1:4).eq.'3DA0') callsign='3D0'//callsign(5:6)

 ! Work-around for Guinea prefixes:
   if(callsign(1:2).eq.'3X' .and. callsign(3:3).ge.'A' .and.          &
        callsign(3:3).le.'Z') callsign='Q'//callsign(3:6)

   if(callsign(1:3).eq.'CQ ') then
      ncall=NBASE + 1
      if(callsign(4:4).ge.'0' .and. callsign(4:4).le.'9' .and.        &
           callsign(5:5).ge.'0' .and. callsign(5:5).le.'9' .and.      &
           callsign(6:6).ge.'0' .and. callsign(6:6).le.'9') then
         read(callsign(4:6),*) nfreq
         ncall=NBASE + 3 + nfreq
      endif
      return
   else if(callsign(1:4).eq.'QRZ ') then
      ncall=NBASE + 2
      return
   else if(callsign(1:3).eq.'DE ') then
      ncall=267796945
      return
   endif

   tmp='      '
   if(callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') then
      tmp=callsign
   else if(callsign(2:2).ge.'0' .and. callsign(2:2).le.'9') then
      if(callsign(6:6).ne.' ') then
         text=.true.
         return
      endif
      tmp=' '//callsign(:5)
   else
      text=.true.
      return
   endif

   do i=1,6
      c=tmp(i:i)
      if(c.ge.'a' .and. c.le.'z')                                &
           tmp(i:i)=char(ichar(c)-ichar('a')+ichar('A'))
   enddo

   n1=0
   if((tmp(1:1).ge.'A'.and.tmp(1:1).le.'Z').or.tmp(1:1).eq.' ') n1=1
   if(tmp(1:1).ge.'0' .and. tmp(1:1).le.'9') n1=1
   n2=0
   if(tmp(2:2).ge.'A' .and. tmp(2:2).le.'Z') n2=1
   if(tmp(2:2).ge.'0' .and. tmp(2:2).le.'9') n2=1
   n3=0
   if(tmp(3:3).ge.'0' .and. tmp(3:3).le.'9') n3=1
   n4=0
   if((tmp(4:4).ge.'A'.and.tmp(4:4).le.'Z').or.tmp(4:4).eq.' ') n4=1
   n5=0
   if((tmp(5:5).ge.'A'.and.tmp(5:5).le.'Z').or.tmp(5:5).eq.' ') n5=1
   n6=0
   if((tmp(6:6).ge.'A'.and.tmp(6:6).le.'Z').or.tmp(6:6).eq.' ') n6=1

   if(n1+n2+n3+n4+n5+n6 .ne. 6) then
      text=.true.
      return
   endif

   ncall=nchar(tmp(1:1))
   ncall=36*ncall+nchar(tmp(2:2))
   ncall=10*ncall+nchar(tmp(3:3))
   ncall=27*ncall+nchar(tmp(4:4))-10
   ncall=27*ncall+nchar(tmp(5:5))-10
   ncall=27*ncall+nchar(tmp(6:6))-10

   return
 end subroutine packcall

 subroutine unpackcall(ncall,word,iv2,psfx)

   parameter (NBASE=37*36*10*27*27*27)
   character word*12,c*37,psfx*4

   data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ '/

   word='......'
   psfx='    '
   n=ncall
   iv2=0
   if(n.ge.262177560) go to 20
   word='......'
 !  if(n.ge.262177560) go to 999            !Plain text message ...
   i=mod(n,27)+11
   word(6:6)=c(i:i)
   n=n/27
   i=mod(n,27)+11
   word(5:5)=c(i:i)
   n=n/27
   i=mod(n,27)+11
   word(4:4)=c(i:i)
   n=n/27
   i=mod(n,10)+1
   word(3:3)=c(i:i)
   n=n/10
   i=mod(n,36)+1
   word(2:2)=c(i:i)
   n=n/36
   i=n+1
   word(1:1)=c(i:i)
   do i=1,4
      if(word(i:i).ne.' ') go to 10
   enddo
   go to 999
 10 word=word(i:)
   go to 999

 20 if(n.ge.267796946) go to 999

 ! We have a JT65v2 message
   if((n.ge.262178563) .and. (n.le.264002071)) then
 ! CQ with prefix
      iv2=1
      n=n-262178563
      i=mod(n,37)+1
      psfx(4:4)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if((n.ge.264002072) .and. (n.le.265825580)) then
 ! QRZ with prefix
      iv2=2
      n=n-264002072
      i=mod(n,37)+1
      psfx(4:4)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if((n.ge.265825581) .and. (n.le.267649089)) then
 ! DE with prefix
      iv2=3
      n=n-265825581
      i=mod(n,37)+1
      psfx(4:4)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if((n.ge.267649090) .and. (n.le.267698374)) then
 ! CQ with suffix
      iv2=4
      n=n-267649090
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if((n.ge.267698375) .and. (n.le.267747659)) then
 ! QRZ with suffix
      iv2=5
      n=n-267698375
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if((n.ge.267747660) .and. (n.le.267796944)) then
 ! DE with suffix
      iv2=6
      n=n-267747660
      i=mod(n,37)+1
      psfx(3:3)=c(i:i)
      n=n/37
      i=mod(n,37)+1
      psfx(2:2)=c(i:i)
      n=n/37
      i=n+1
      psfx(1:1)=c(i:i)

   else if(n.eq.267796945) then
 ! DE with no prefix or suffix
      iv2=7
      psfx = '    '
   endif

999 if(word(1:3).eq.'3D0') word='3DA0'//word(4:)
   if(word(1:1).eq.'Q' .and. word(2:2).ge.'A' .and.                 &
        word(2:2).le.'Z') word='3X'//word(2:)

   return
 end subroutine unpackcall

 subroutine packgrid(grid,ng,text)

   parameter (NGBASE=180*180)
   character*4 grid
   character*1 c1
   logical text

   text=.false.
   if(grid.eq.'    ') go to 90               !Blank grid is OK

 ! First, handle signal reports in the original range, -01 to -30 dB
   if(grid(1:1).eq.'-') then
      read(grid(2:3),*,err=800,end=800) n
      if(n.ge.1 .and. n.le.30) then
         ng=NGBASE+1+n
         go to 900
      endif
      go to 10
   else if(grid(1:2).eq.'R-') then
      read(grid(3:4),*,err=800,end=800) n
      if(n.ge.1 .and. n.le.30) then
         ng=NGBASE+31+n
         go to 900
      endif
      go to 10
 ! Now check for RO, RRR, or 73 in the message field normally used for grid
   else if(grid(1:4).eq.'RO  ') then
      ng=NGBASE+62
      go to 900
   else if(grid(1:4).eq.'RRR ') then
      ng=NGBASE+63
      go to 900
   else if(grid(1:4).eq.'73  ') then
      ng=NGBASE+64
      go to 900
   endif

 ! Now check for extended-range signal reports: -50 to -31, and 0 to +49.
 10 n=99
   c1=grid(1:1)
   read(grid,*,err=20,end=20) n
   go to 30
 20 read(grid(2:4),*,err=30,end=30) n
 30 if(n.ge.-50 .and. n.le.49) then
      if(c1.eq.'R') then
         write(grid,1002) n+50
 1002    format('LA',i2.2)
      else
         write(grid,1003) n+50
 1003    format('KA',i2.2)
      endif
      go to 40
   endif

 ! Maybe it's free text ?
   if(grid(1:1).lt.'A' .or. grid(1:1).gt.'R') text=.true.
   if(grid(2:2).lt.'A' .or. grid(2:2).gt.'R') text=.true.
   if(grid(3:3).lt.'0' .or. grid(3:3).gt.'9') text=.true.
   if(grid(4:4).lt.'0' .or. grid(4:4).gt.'9') text=.true.
   if(text) go to 900

 ! OK, we have a properly formatted grid locator
 40 call grid2deg(grid//'mm',dlong,dlat)
   long=int(dlong)
   lat=int(dlat+ 90.0)
   ng=((long+180)/2)*180 + lat
   go to 900

 90 ng=NGBASE + 1
   go to 900

 800 text=.true.
 900 continue

   return
 end subroutine packgrid

 subroutine unpackgrid(ng,grid)

   parameter (NGBASE=180*180)
   character grid*4,grid6*6

   grid='    '
   if(ng.ge.32400) go to 10
   dlat=mod(ng,180)-90
   dlong=(ng/180)*2 - 180 + 2
   call deg2grid(dlong,dlat,grid6)
   grid=grid6(:4)
   if(grid(1:2).eq.'KA') then
      read(grid(3:4),*) n
      n=n-50
      write(grid,1001) n
 1001 format(i3.2)
      if(grid(1:1).eq.' ') grid(1:1)='+'
   else if(grid(1:2).eq.'LA') then
      read(grid(3:4),*) n
      n=n-50
      write(grid,1002) n
 1002 format('R',i3.2)
      if(grid(2:2).eq.' ') grid(2:2)='+'
   endif
   go to 900

 10 n=ng-NGBASE-1
   if(n.ge.1 .and.n.le.30) then
      write(grid,1012) -n
 1012 format(i3.2)
   else if(n.ge.31 .and.n.le.60) then
      n=n-30
      write(grid,1022) -n
 1022 format('R',i3.2)
   else if(n.eq.61) then
      grid='RO'
   else if(n.eq.62) then
      grid='RRR'
   else if(n.eq.63) then
      grid='73'
   endif

 900 return
 end subroutine unpackgrid

 subroutine packmsg(msg0,dat,itype,bcontest)

 ! Packs a JT4/JT9/JT65 message into twelve 6-bit symbols

 ! itype Message Type
 !--------------------
 !   1   Standardd message
 !   2   Type 1 prefix
 !   3   Type 1 suffix
 !   4   Type 2 prefix
 !   5   Type 2 suffix
 !   6   Free text
 !  -1   Does not decode correctly

   parameter (NBASE=37*36*10*27*27*27)
   parameter (NBASE2=262178562)
   character*22 msg0,msg
   integer dat(:)
   character*12 c1,c2
   character*4 c3
   character*6 grid6
   logical text1,text2,text3,bcontest

   itype=1
   if(bcontest) then
      !call to_contest_msg(msg0,msg)
      ! this causes problems with freetext ala, KN4CRD DE KN4CRD -13 R
      msg=msg0
   else
      msg=msg0
   end if

   call fmtmsg(msg,iz)
   if(msg(1:3).eq.'CQ ' .and. msg(4:4).ge.'0' .and. msg(4:4).le.'9'   &
        .and. msg(5:5).eq.' ') msg='CQ 00'//msg(4:)

   if(msg(1:6).eq.'CQ DX ') msg(3:3)='9'
   if(msg(1:3).eq.'CQ ' .and.                                         &
        msg(4:4).ge.'A' .and. msg(4:4).le.'Z' .and.                   &
        msg(5:5).ge.'A' .and. msg(5:5).le.'Z' .and.                   &
        msg(6:6).eq.' ') msg='E9'//msg(4:)

 ! See if it's a CQ message
   if(msg(1:3).eq.'CQ ') then
      i=3
 ! ... and if so, does it have a reply frequency?
      if(msg(4:4).ge.'0' .and. msg(4:4).le.'9' .and.                  &
           msg(5:5).ge.'0' .and. msg(5:5).le.'9' .and.                &
           msg(6:6).ge.'0' .and. msg(6:6).le.'9') i=7
      go to 1
   endif

   do i=1,22
      if(msg(i:i).eq.' ') go to 1       !Get 1st blank
   enddo
   go to 10                             !Consider msg as plain text

 1 ia=i
   c1=msg(1:ia-1)
   do i=ia+1,22
      if(msg(i:i).eq.' ') go to 2       !Get 2nd blank
   enddo
   go to 10                             !Consider msg as plain text

 2 ib=i
   c2=msg(ia+1:ib-1)

   do i=ib+1,22
      if(msg(i:i).eq.' ') go to 3       !Get 3rd blank
   enddo
   go to 10                             !Consider msg as plain text

 3 ic=i
   c3='    '
   if(ic.ge.ib+1) c3=msg(ib+1:ic)
   if(c3.eq.'OOO ') c3='    '           !Strip out the OOO flag
   call getpfx1(c1,k1,nv2a)
   if(nv2a.ge.4) go to 10
   call packcall(c1,nc1,text1)
   if(text1) go to 10
   call getpfx1(c2,k2,nv2b)
   call packcall(c2,nc2,text2)
   if(text2) go to 10
   if(nv2a.eq.2 .or. nv2a.eq.3 .or. nv2b.eq.2 .or. nv2b.eq.3) then
      if(k1.lt.0 .or. k2.lt.0 .or. k1*k2.ne.0) go to 10
      if(k2.gt.0) k2=k2+450
      k=max(k1,k2)
      if(k.gt.0) then
         call k2grid(k,grid6)
         c3=grid6(:4)
      endif
   endif
   call packgrid(c3,ng,text3)

   if(nv2a.lt.4 .and. nv2b.lt.4 .and. (.not.text1) .and. (.not.text2) .and.  &
        (.not.text3)) go to 20

   nc1=0
   if(nv2b.eq.4) then
      if(c1(1:3).eq.'CQ ' .and. (.not.text3))  nc1=262178563 + k2
      if(c1(1:4).eq.'QRZ ' .and. (.not.text3)) nc1=264002072 + k2
      if(c1(1:3).eq.'DE ' .and. (.not.text3))  nc1=265825581 + k2
   else if(nv2b.eq.5) then
      if(c1(1:3).eq.'CQ ' .and. (.not.text3))  nc1=267649090 + k2
      if(c1(1:4).eq.'QRZ ' .and. (.not.text3)) nc1=267698375 + k2
      if(c1(1:3).eq.'DE ' .and. (.not.text3))  nc1=267747660 + k2
   endif
   if(nc1.ne.0) go to 20

 ! The message will be treated as plain text.
 10 itype=6
   call packtext(msg,nc1,nc2,ng)
   ng=ng+32768

 ! Encode data into 6-bit words
20 continue
   if(itype.ne.6) itype=max(nv2a,nv2b)
   jt_itype=itype
   jt_c1=c1(1:6)
   jt_c2=c2(1:6)
   jt_c3=c3
   jt_k1=k1
   jt_k2=k2
   jt_nc1=nc1
   jt_nc2=nc2
   jt_ng=ng
   dat(1)=iand(ishft(nc1,-22),63)                !6 bits
   dat(2)=iand(ishft(nc1,-16),63)                !6 bits
   dat(3)=iand(ishft(nc1,-10),63)                !6 bits
   dat(4)=iand(ishft(nc1, -4),63)                !6 bits
   dat(5)=4*iand(nc1,15)+iand(ishft(nc2,-26),3)  !4+2 bits
   dat(6)=iand(ishft(nc2,-20),63)                !6 bits
   dat(7)=iand(ishft(nc2,-14),63)                !6 bits
   dat(8)=iand(ishft(nc2, -8),63)                !6 bits
   dat(9)=iand(ishft(nc2, -2),63)                !6 bits
   dat(10)=16*iand(nc2,3)+iand(ishft(ng,-12),15) !2+4 bits
   dat(11)=iand(ishft(ng,-6),63)
   dat(12)=iand(ng,63)

   return
 end subroutine packmsg

 subroutine unpackmsg(dat,msg,bcontest,mygrid)

   parameter (NBASE=37*36*10*27*27*27)
   parameter (NGBASE=180*180)
   integer dat(:)
   character c1*12,c2*12,grid*4,msg*22,grid6*6,psfx*4,junk2*4,mygrid*6
   logical cqnnn,bcontest

   cqnnn=.false.
   nc1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+         &
        ishft(dat(4),4) + iand(ishft(dat(5),-2),15)

   nc2=ishft(iand(dat(5),3),26) + ishft(dat(6),20) +                   &
        ishft(dat(7),14) + ishft(dat(8),8) + ishft(dat(9),2) +         &
        iand(ishft(dat(10),-4),3)

   ng=ishft(iand(dat(10),15),12) + ishft(dat(11),6) + dat(12)

   if(ng.ge.32768) then
      call unpacktext(nc1,nc2,ng,msg)
      go to 100
   endif

   call unpackcall(nc1,c1,iv2,psfx)
   if(iv2.eq.0) then
 ! This is an "original JT65" message
      if(nc1.eq.NBASE+1) c1='CQ    '
      if(nc1.eq.NBASE+2) c1='QRZ   '
      nfreq=nc1-NBASE-3
      if(nfreq.ge.0 .and. nfreq.le.999) then
         write(c1,1002) nfreq
 1002    format('CQ ',i3.3)
         cqnnn=.true.
      endif
   endif

   call unpackcall(nc2,c2,junk1,junk2)
   call unpackgrid(ng,grid)

   if(iv2.gt.0) then
 ! This is a JT65v2 message
      do i=1,4
         if(ichar(psfx(i:i)).eq.0) psfx(i:i)=' '
      enddo

      n1=len_trim(psfx)
      n2=len_trim(c2)
      if(iv2.eq.1) msg='CQ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.2) msg='QRZ '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.3) msg='DE '//psfx(:n1)//'/'//c2(:n2)//' '//grid
      if(iv2.eq.4) msg='CQ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.5) msg='QRZ '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.6) msg='DE '//c2(:n2)//'/'//psfx(:n1)//' '//grid
      if(iv2.eq.7) then
         grid6=grid//'ma'
         call grid2k(grid6,k)
         if(k.ge.451 .and. k.le.900) then
            call getpfx2(k,c2)
            n2=len_trim(c2)
            msg='DE '//c2(:n2)
         else
            msg='DE '//c2(:n2)//' '//grid
         endif
      endif
      if(iv2.eq.8) msg=' '
      go to 100
   else

   endif

   grid6=grid//'ma'
   call grid2k(grid6,k)
   if(k.ge.1 .and. k.le.450)   call getpfx2(k,c1)
   if(k.ge.451 .and. k.le.900) call getpfx2(k,c2)

   i=index(c1,char(0))
   if(i.ge.3) c1=c1(1:i-1)//'            '
   i=index(c2,char(0))
   if(i.ge.3) c2=c2(1:i-1)//'            '

   msg='                      '
   j=0
   if(cqnnn) then
      msg=c1//'          '
      j=7                                  !### ??? ###
      go to 10
   endif

   do i=1,12
      j=j+1
      msg(j:j)=c1(i:i)
      if(c1(i:i).eq.' ') go to 10
   enddo
   j=j+1
   msg(j:j)=' '

 10 do i=1,12
      if(j.le.21) j=j+1
      msg(j:j)=c2(i:i)
      if(c2(i:i).eq.' ') go to 20
   enddo
   if(j.le.21) j=j+1
   msg(j:j)=' '

 20 if(k.eq.0) then
      do i=1,4
         if(j.le.21) j=j+1
         msg(j:j)=grid(i:i)
      enddo
      if(j.le.21) j=j+1
      msg(j:j)=' '
   endif

 100 continue
   if(msg(1:6).eq.'CQ9DX ') msg(3:3)=' '
   if(msg(1:2).eq.'E9' .and.                                          &
        msg(3:3).ge.'A' .and. msg(3:3).le.'Z' .and.                   &
        msg(4:4).ge.'A' .and. msg(4:4).le.'Z' .and.                   &
        msg(5:5).eq.' ') msg='CQ '//msg(3:)

   if(bcontest) call fix_contest_msg(mygrid,msg)

   if(msg(1:5).eq.'CQ 00' .and. msg(6:6).ge.'0' .and.                 &
        msg(6:6).le.'9') msg='CQ '//msg(6:)

   return
 end subroutine unpackmsg

 subroutine packtext(msg,nc1,nc2,nc3)

   parameter (MASK28=2**28 - 1)
   character*22 msg
   character*42 c
   data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +-./?'/

   nc1=0
   nc2=0
   nc3=0

   do i=1,5                                !First 5 characters in nc1
      do j=1,42                            !Get character code
         if(msg(i:i).eq.c(j:j)) go to 10
      enddo
      j=37
 10   j=j-1                                !Codes should start at zero
      nc1=42*nc1 + j
   enddo

   do i=6,10                               !Characters 6-10 in nc2
      do j=1,42                            !Get character code
         if(msg(i:i).eq.c(j:j)) go to 20
      enddo
      j=37
 20   j=j-1                                !Codes should start at zero
      nc2=42*nc2 + j
   enddo

   do i=11,13                              !Characters 11-13 in nc3
      do j=1,42                            !Get character code
         if(msg(i:i).eq.c(j:j)) go to 30
      enddo
      j=37
 30   j=j-1                                !Codes should start at zero
      nc3=42*nc3 + j
   enddo

 ! We now have used 17 bits in nc3.  Must move one each to nc1 and nc2.
   nc1=nc1+nc1
   if(iand(nc3,32768).ne.0) nc1=nc1+1
   nc2=nc2+nc2
   if(iand(nc3,65536).ne.0) nc2=nc2+1
   nc3=iand(nc3,32767)

   return
 end subroutine packtext

 subroutine unpacktext(nc1,nc2,nc3,msg)

   character*22 msg
   character*44 c
   data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +-./?'/

   nc3=iand(nc3,32767)                      !Remove the "plain text" bit
   if(iand(nc1,1).ne.0) nc3=nc3+32768
   nc1=nc1/2
   if(iand(nc2,1).ne.0) nc3=nc3+65536
   nc2=nc2/2

   do i=5,1,-1
      j=mod(nc1,42)+1
      msg(i:i)=c(j:j)
      nc1=nc1/42
   enddo

   do i=10,6,-1
      j=mod(nc2,42)+1
      msg(i:i)=c(j:j)
      nc2=nc2/42
   enddo

   do i=13,11,-1
      j=mod(nc3,42)+1
      msg(i:i)=c(j:j)
      nc3=nc3/42
   enddo
   msg(14:22) = '         '

   return
 end subroutine unpacktext

 subroutine getpfx1(callsign,k,nv2)

   character*12 callsign0,callsign,lof,rof
   character*8 c
   character addpfx*8,tpfx*4,tsfx*3
   logical ispfx,issfx,invalid
   common/pfxcom/addpfx
   include 'pfx.f90'

   callsign0=callsign
   nv2=1
   iz=index(callsign,' ') - 1
   if(iz.lt.0) iz=12
   islash=index(callsign(1:iz),'/')
   k=0
 !  if(k.eq.0) go to 10     !Tnx to DL9RDZ for reminder:this was for tests only!
   c='   '
   if(islash.gt.0 .and. islash.le.(iz-4)) then
 ! Add-on prefix
      c=callsign(1:islash-1)
      callsign=callsign(islash+1:iz)
      do i=1,NZ
         if(pfx(i)(1:4).eq.c) then
            k=i
            nv2=2
            go to 10
         endif
      enddo
      if(addpfx.eq.c) then
         k=449
         nv2=2
         go to 10
      endif

   else if(islash.eq.(iz-1)) then
 ! Add-on suffix
      c=callsign(islash+1:iz)
      callsign=callsign(1:islash-1)
      do i=1,NZ2
         if(sfx(i).eq.c(1:1)) then
            k=400+i
            nv2=3
            go to 10
         endif
      enddo
   endif

 10 if(islash.ne.0 .and.k.eq.0) then
 ! Original JT65 would force this compound callsign to be treated as
 ! plain text.  In JT65v2, we will encode the prefix or suffix into nc1.
 ! The task here is to compute the proper value of k.
      lof=callsign0(:islash-1)
      rof=callsign0(islash+1:)
      llof=len_trim(lof)
      lrof=len_trim(rof)
      ispfx=(llof.gt.0 .and. llof.le.4)
      issfx=(lrof.gt.0 .and. lrof.le.3)
      invalid=.not.(ispfx.or.issfx)
      if(ispfx.and.issfx) then
         if(llof.lt.3) issfx=.false.
         if(lrof.lt.3) ispfx=.false.
         if(ispfx.and.issfx) then
            i=ichar(callsign0(islash-1:islash-1))
            if(i.ge.ichar('0') .and. i.le.ichar('9')) then
               issfx=.false.
            else
               ispfx=.false.
            endif
         endif
      endif

      if(invalid) then
         k=-1
      else
         if(ispfx) then
            tpfx=lof(1:4)
            k=nchar(tpfx(1:1))
            k=37*k + nchar(tpfx(2:2))
            k=37*k + nchar(tpfx(3:3))
            k=37*k + nchar(tpfx(4:4))
            nv2=4
            i=index(callsign0,'/')
            callsign=callsign0(:i-1)
            callsign=callsign0(i+1:)
         endif
         if(issfx) then
            tsfx=rof(1:3)
            k=nchar(tsfx(1:1))
            k=37*k + nchar(tsfx(2:2))
            k=37*k + nchar(tsfx(3:3))
            nv2=5
            i=index(callsign0,'/')
            callsign=callsign0(:i-1)
         endif
      endif
   endif

   return
 end subroutine getpfx1

 subroutine getpfx2(k0,callsign)

   character callsign*12
   include 'pfx.f90'
   character addpfx*8
   common/pfxcom/addpfx

   k=k0
   if(k.gt.450) k=k-450
   if(k.ge.1 .and. k.le.NZ) then
      iz=index(pfx(k),' ') - 1
      callsign=pfx(k)(1:iz)//'/'//callsign
   else if(k.ge.401 .and. k.le.400+NZ2) then
      iz=index(callsign,' ') - 1
      callsign=callsign(1:iz)//'/'//sfx(k-400)
   else if(k.eq.449) then
      iz=index(addpfx,' ') - 1
      if(iz.lt.1) iz=8
      callsign=addpfx(1:iz)//'/'//callsign
   endif

   return
 end subroutine getpfx2

 subroutine grid2k(grid,k)

   character*6 grid

   call grid2deg(grid,xlong,xlat)
   nlong=nint(xlong)
   nlat=nint(xlat)
   k=0
   if(nlat.ge.85) k=5*(nlong+179)/2 + nlat-84

   return
 end subroutine grid2k

 subroutine k2grid(k,grid)
   character grid*6

   nlong=2*mod((k-1)/5,90)-179
   if(k.gt.450) nlong=nlong+180
   nlat=mod(k-1,5)+ 85
   dlat=nlat
   dlong=nlong
   call deg2grid(dlong,dlat,grid)

   return
 end subroutine k2grid

 subroutine grid2n(grid,n)
   character*4 grid

   i1=ichar(grid(1:1))-ichar('A')
   i2=ichar(grid(3:3))-ichar('0')
   i=10*i1 + i2
   n=-i - 31

   return
 end subroutine grid2n

 subroutine n2grid(n,grid)
   character*4 grid

   if(n.gt.-31 .or. n.lt.-70) stop 'Error in n2grid'
   i=-(n+31)                           !NB: 0 <= i <= 39
   i1=i/10
   i2=mod(i,10)
   grid(1:1)=char(ichar('A')+i1)
   grid(2:2)='A'
   grid(3:3)=char(ichar('0')+i2)
   grid(4:4)='0'

   return
 end subroutine n2grid

 function nchar(c)

 ! Convert ascii number, letter, or space to 0-36 for callsign packing.

   character c*1

   n=0                                    !Silence compiler warning
   if(c.ge.'0' .and. c.le.'9') then
      n=ichar(c)-ichar('0')
   else if(c.ge.'A' .and. c.le.'Z') then
      n=ichar(c)-ichar('A') + 10
   else if(c.ge.'a' .and. c.le.'z') then
      n=ichar(c)-ichar('a') + 10
   else if(c.ge.' ') then
      n=36
   else
      Print*,'Invalid character in callsign ',c,' ',ichar(c)
      stop
   endif
   nchar=n

   return
 end function nchar

 subroutine pack50(n1,n2,dat)

   integer*1 dat(:),i1

   i1=iand(ishft(n1,-20),255)                !8 bits
   dat(1)=i1
   i1=iand(ishft(n1,-12),255)                 !8 bits
   dat(2)=i1
   i1=iand(ishft(n1, -4),255)                 !8 bits
   dat(3)=i1
   i1=16*iand(n1,15)+iand(ishft(n2,-18),15)   !4+4 bits
   dat(4)=i1
   i1=iand(ishft(n2,-10),255)                 !8 bits
   dat(5)=i1
   i1=iand(ishft(n2, -2),255)                 !8 bits
   dat(6)=i1
   i1=64*iand(n2,3)                           !2 bits
   dat(7)=i1
   dat(8)=0
   dat(9)=0
   dat(10)=0
   dat(11)=0

   return
 end subroutine pack50

subroutine packpfx(call1,n1,ng,nadd)

  character*12 call1,call0
  character*3 pfx
  logical text

  i1=index(call1,'/')
  if(call1(i1+2:i1+2).eq.' ') then
! Single-character add-on suffix (maybe also fourth suffix letter?)
     call0=call1(:i1-1)
     call packcall(call0,n1,text)
     nadd=1
     nc=ichar(call1(i1+1:i1+1))
     if(nc.ge.48 .and. nc.le.57) then
        n=nc-48
     else if(nc.ge.65 .and. nc.le.90) then
        n=nc-65+10
     else
        n=38
     endif
     nadd=1
     ng=60000-32768+n
  else if(call1(i1+3:i1+3).eq.' ') then
! Two-character numerical suffix, /10 to /99
     call0=call1(:i1-1)
     call packcall(call0,n1,text)
     nadd=1
     n=10*(ichar(call1(i1+1:i1+1))-48) + ichar(call1(i1+2:i1+2)) - 48
     nadd=1
     ng=60000 + 26 + n
  else
! Prefix of 1 to 3 characters
     pfx=call1(:i1-1)
     if(pfx(3:3).eq.' ') pfx=' '//pfx(1:2)
     if(pfx(3:3).eq.' ') pfx=' '//pfx(1:2)
     call0=call1(i1+1:)
     call packcall(call0,n1,text)

     ng=0
     do i=1,3
        nc=ichar(pfx(i:i))
        if(nc.ge.48 .and. nc.le.57) then
           n=nc-48
        else if(nc.ge.65 .and. nc.le.90) then
           n=nc-65+10
        else
           n=36
        endif
        ng=37*ng + n
     enddo
     nadd=0
     if(ng.ge.32768) then
        ng=ng-32768
        nadd=1
     endif
  endif

  return
end subroutine packpfx

end module packjt

license: mit


code:
subroutine pctile(x,npts,npct,xpct)

  parameter (NMAX=128*1024)
  real*4 x(npts)
  real*4 tmp(NMAX)

  if(npts.le.0) then
     xpct=1.0
     go to 900
  endif
  if(npts.gt.NMAX) stop

  tmp(1:npts)=x
  call shell(npts,tmp)
  j=nint(npts*0.01*npct)
  if(j.lt.1) j=1
  if(j.gt.npts) j=npts
  xpct=tmp(j)

900 continue
  return
end subroutine pctile


code:
subroutine peakdt9(c2,nsps8,nspsd,c3,xdt)

  parameter (NZ2=1512,NZ3=1360)
  complex c2(0:NZ2-1)
  complex c3(0:NZ3-1)
  complex z
  real p(0:3300)
  include 'jt9sync.f90'

  p=0.
  i0=5*nspsd
  do i=0,NZ2-1
     z=1.e-3*sum(c2(max(i-(nspsd-1),0):i))
     p(i0+i)=real(z)**2 + aimag(z)**2      !Integrated symbol power at freq=0
  enddo

  call getlags(nsps8,lag0,lag1,lag2)
  tsymbol=nsps8/1500.0
  dtlag=tsymbol/nspsd
  smax=0.
  lagpk=0
  do lag=lag1,lag2
     sum0=0.
     sum1=0.
     j=-nspsd
     do i=1,85
        j=j+nspsd
        if(isync(i).eq.1) then
           sum1=sum1+p(j+lag)
        else
           sum0=sum0+p(j+lag)
        endif
     enddo
     ss=(sum1/16.0)/(sum0/69.0) - 1.0
     xdt=(lag-lag0)*dtlag
     if(ss.gt.smax) then
        smax=ss
        lagpk=lag
     endif
  enddo

  xdt=(lagpk-lag0)*dtlag

  do i=0,NZ3-1
     j=i+lagpk-i0-nspsd+1
     if(j.ge.0 .and. j.lt.NZ2) then
        c3(i)=c2(j)
     else
        c3(i)=0.
     endif
  enddo

  return
end subroutine peakdt9


code:
subroutine peakup(ym,y0,yp,dx)

  b=(yp-ym)/2.0
  c=(yp+ym-2.0*y0)/2.0
  dx=-b/(2.0*c)

  return
end subroutine peakup


code:
subroutine plotsave(swide,nw,nh,irow)

  real, dimension(:,:), allocatable :: sw
  real swide(0:nw-1)
  data nw0/-1/,nh0/-1/
  save nw0,nh0,sw

  if(irow.eq.-99) then
     if(allocated(sw)) deallocate(sw)
     go to 900
  endif

  if(nw.ne.nw0 .or. nh.ne.nh0 .or. (.not.allocated(sw))) then
     if(allocated(sw)) deallocate(sw)
!     if(nw0.ne.-1) deallocate(sw)
     allocate(sw(0:nw-1,0:nh-1))
     sw=0.
     nw0=nw
     nh0=nh
  endif
  df=12000.0/16384
  if(irow.lt.0) then
! Push a new row of data into sw
     do j=nh-1,1,-1
        sw(0:nw-1,j)=sw(0:nw-1,j-1)
     enddo
     sw(0:nw-1,0)=swide
  else
! Return the saved "irow" as swide(), for a waterfall replot.
     swide=sw(0:nw-1,irow)
  endif

900 return
end subroutine plotsave


code:
subroutine polyfit(x,y,sigmay,npts,nterms,mode,a,chisqr)
  implicit real*8 (a-h,o-z)
  real*8 x(npts), y(npts), sigmay(npts), a(nterms)
  real*8 sumx(10), sumy(10), array(10,10)

! Accumulate weighted sums
  nmax = 2*nterms-1
  sumx=0.
  sumy=0.
  chisq=0.
  do i=1,npts
     xi=x(i)
     yi=y(i)
     if(mode.lt.0) then
        weight=1./abs(yi)
     else if(mode.eq.0) then
        weight=1
     else
        weight=1./sigmay(i)**2
     end if
     xterm=weight
     do n=1,nmax
        sumx(n)=sumx(n)+xterm
        xterm=xterm*xi
     enddo
     yterm=weight*yi
     do n=1,nterms
        sumy(n)=sumy(n)+yterm
        yterm=yterm*xi
     enddo
     chisq=chisq+weight*yi**2
  enddo

! Construct matrices and calculate coefficients
  do j=1,nterms
     do k=1,nterms
        n=j+k-1
        array(j,k)=sumx(n)
     enddo
  enddo

  delta=determ(array,nterms)
  if(delta.eq.0) then
     chisqr=0.
     a=0.
  else
     do l=1,nterms
        do j=1,nterms
           do k=1,nterms
              n=j+k-1
              array(j,k)=sumx(n)
           enddo
           array(j,l)=sumy(j)
        enddo
        a(l)=determ(array,nterms)/delta
     enddo

! Calculate chi square

     do j=1,nterms
        chisq=chisq-2*a(j)*sumy(j)
        do k=1,nterms
           n=j+k-1
           chisq=chisq+a(j)*a(k)*sumx(n)
        enddo
     enddo
     free=npts-nterms
     chisqr=chisq/free
  end if

  return
end subroutine polyfit


code:
MODULE prog_args
  CHARACTER(len=80) :: shm_key
  CHARACTER(len=500) :: exe_dir = '.', data_dir = '.', temp_dir = '.'
END MODULE prog_args


code:
subroutine ps4(dat,nfft,s)

  parameter (NMAX=2520+2)
  parameter (NHMAX=NMAX/2-1)
  real dat(nfft)
  real dat2(NMAX)
  real s(NHMAX)
  complex c(0:NMAX)
  equivalence(dat2,c)

  nh=nfft/2
  do i=1,nh
     dat2(i)=dat(i)/128.0       !### Why 128 ??
  enddo
  do i=nh+1,nfft
     dat2(i)=0.
  enddo

  call four2a(c,nfft,1,-1,0)

  fac=1.0/nfft
  do i=1,nh
     s(i)=fac*(real(c(i))**2 + aimag(c(i))**2)
  enddo

  return
end subroutine ps4


code:
program psk_parse

  character line*256,callsign*12,callsign0*12,progname*30
  integer nc(6),ntot(6),nsingle(6)
  logical zap
  data callsign0/'            '/

  open(10,file='jt65-2',status='old')
  nc=0
  ntot=0
  ncalls=0
  nsingle=0
  zap=.false.

  do iline=1,9999999
     read(10,'(a256)',end=900) line
     n=len(trim(line))
     i1=0
     i2=0
     i3=0
     do i=1,n
        if(ichar(line(i:i)).eq.9) then
           if(i1.eq.0) then
              i1=i
              cycle
           endif
           if(i2.eq.0) then
              i2=i
              cycle
           endif
           if(i3.eq.0) then
              i3=i
              exit
           endif
        endif
     enddo
     callsign=line(1:i1-1)

     if(zap) then
        if(callsign(1:1).ge.'0' .and. callsign(1:1).le.'9' .and.               &
           callsign(2:2).ge.'0' .and. callsign(2:2).le.'9' .and.               &
           callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') cycle

        if(callsign(1:1).ge.'A' .and. callsign(1:1).le.'Z' .and.               &
           callsign(2:2).ge.'A' .and. callsign(2:2).le.'Z' .and.               &
           callsign(3:3).ge.'A' .and. callsign(3:3).le.'Z' .and.               &
           callsign(4:4).ge.'A' .and. callsign(4:4).le.'Z') cycle

        if(callsign(1:1).eq.'N' .and.                                        &
           callsign(2:2).ge.'0' .and. callsign(2:2).le.'9' .and.               &
           callsign(3:3).ge.'0' .and. callsign(3:3).le.'9') cycle

     endif

     progname=line(i1+1:i2-1)
     do i=1,len(trim(progname))
        if(progname(i:i).eq.' ') progname(i:i)='_'
     enddo
     read(line(i2+1:i3-1),*) ndecodes
     read(line(i3+1:),*) nreporters

     j=6
     if(index(progname,'WSJT-X').gt.0) j=1
     if(index(progname,'HB9HQX').gt.0) j=2
     if(index(progname,'JTDX').gt.0) j=3
     if(index(progname,'Comfort').gt.0) j=4
     if(index(progname,'COMFORT').gt.0) j=4
     if(index(progname,'JT65-HF').gt.0 .and. j.eq.6) j=5

     nctot=sum(nc)
     if(callsign.ne.callsign0 .and. nctot.gt.0) then
        write(13,1000) callsign0,nc,nctot
1000    format(a12,6i8,i10)
        if(nctot.eq.1) nsingle(j)=nsingle(j)+1
        nc=0
        callsign0=callsign
        ncalls=ncalls+1
     endif
     nc(j)=nc(j) + ndecodes
     ntot(j)=ntot(j) + ndecodes
     write(12,1010) iline,callsign,ndecodes,nreporters,progname
1010 format(i8,2x,a12,2x,2i8,2x,a30)
  enddo

900 nctot=sum(nc)
  if(nctot.gt.0) write(13,1000) callsign0,nc,nctot

  write(*,1018)
1018 format('   Total    WSJT-X  HB9HQX   JTDX  Comfort JT65-HF   Other'/  &
            '----------------------------------------------------------')
  write(*,1019)
1019 format('Spots reported; fraction of total:')
  write(*,1020) sum(ntot),ntot
1020 format(i8,2x,6i8)
  fac=1.0/sum(ntot)
  write(*,1030) 1.0000,fac*ntot
1030 format(f8.4,2x,6f8.4,f10.4)

  write(*,1038)
1038 format(/'Singletons; as fraction of spots by same program')
  write(*,1040) sum(nsingle),nsingle
1040 format(i8,2x,6i8)
  write(*,1030) float(sum(nsingle))/sum(ntot),float(nsingle)/ntot

  write(*,1050) ncalls
1050 format(/'Distinct calls spotted:',i6)

end program psk_parse


code:
#include <windows.h>
#include <stdio.h>

int ptt_(int *nport, int *ntx, int *iptt)
{
  static HANDLE hFile;
  static int open=0;
  char s[10];
  int i3=1,i4=1,i5=1,i6=1,i9=1,i00=1;

  if(*nport==0) {
    *iptt=*ntx;
    return(0);
  }

  if(*ntx && (!open)) {
    sprintf(s,"\\\\.\\COM%d",*nport);
    hFile=CreateFile(TEXT(s),GENERIC_WRITE,0,NULL,OPEN_EXISTING,
		     FILE_ATTRIBUTE_NORMAL,NULL);
    if(hFile==INVALID_HANDLE_VALUE) {
      //      printf("PTT: Cannot open COM port %d.\n",*nport);
      return 1;
    }
    open=1;
  }

  if(*ntx && open) {
    i3=EscapeCommFunction(hFile,3);
    i5=EscapeCommFunction(hFile,5);
    *iptt=1;
  }

  else {
    i4=EscapeCommFunction(hFile,4);
    i6=EscapeCommFunction(hFile,6);
    i9=EscapeCommFunction(hFile,9);
    i00=CloseHandle(hFile);
    *iptt=0;
    open=0;
  }
  /*
  if(i3==0) return 3;
  if(i4==0) return 4;
  if(i5==0) return 5;
  if(i6==0) return 6;
  if(i9==0) return 9;
  if(i00==0) return 10;
  */
  return 0;
}


code:
FC = gfortran
CC = gcc
CFLAGS = -O2 -Wall -I. -D_WIN32

# Default rules
%.o: %.c
	${CC} ${CFLAGS} -c $<
%.o: %.f
	${FC} ${FFLAGS} -c $<
%.o: %.F
	${FC} ${FFLAGS} -c $<
%.o: %.f90
	${FC} ${FFLAGS} -c $<
%.o: %.F90
	${FC} ${FFLAGS} -c $<

all:	qra64.exe

OBJS1 = main.o qra64.o
qra64.exe: $(OBJS1)
	${CC} -o qra64.exe $(OBJS1) ../qracodes/libqra64.a -lm

OBJS2 = qra64sim.o options.o wavhdr.o
qra64sim.exe: $(OBJS2)
	${FC} -o qra64sim.exe $(OBJS2) ../qracodes/libqra64.a -lm

.PHONY : clean

clean:
	$(RM) *.o qra64.exe qra64sim.exe


code:
/*
main.c
QRA64 mode encode/decode tests

(c) 2016 - Nico Palermo, IV3NWV

Thanks to Andrea Montefusco IW0HDV for his help on adapting the sources
to OSs other than MS Windows

------------------------------------------------------------------------------
This file is part of the qracodes project, a Forward Error Control
encoding/decoding package based on Q-ary RA (Repeat and Accumulate) LDPC codes.

Files in this package:
   main.c		 - this file
   qra64.c/.h     - qra64 mode encode/decoding functions

   ../qracodes/normrnd.{c,h}   - random gaussian number generator
   ../qracodes/npfwht.{c,h}    - Fast Walsh-Hadamard Transforms
   ../qracodes/pdmath.{c,h}    - Elementary math on probability distributions
   ../qracodes/qra12_63_64_irr_b.{c,h} - Tables for a QRA(12,63) irregular RA
                                         code over GF(64)
   ../qracodes/qra13_64_64_irr_e.{c,h} - Tables for a QRA(13,64) irregular RA
                                         code over GF(64)
   ../qracodes/qracodes.{c,h}  - QRA codes encoding/decoding functions

-------------------------------------------------------------------------------

   qracodes is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
   qracodes is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with qracodes source distribution.
   If not, see <http://www.gnu.org/licenses/>.

-----------------------------------------------------------------------------

The code used by the QRA64 mode is the code: QRA13_64_64_IRR_E: K=13
N=64 Q=64 irregular QRA code (defined in qra13_64_64_irr_e.{h,c}).

This code has been designed to include a CRC as the 13th information
symbol and improve the code UER (Undetected Error Rate).  The CRC
symbol is not sent along the channel (the codes are punctured) and the
resulting code is still a (12,63) code with an effective code rate of
R = 12/63.
*/

// OS dependent defines and includes ------------------------------------------

#if _WIN32 // note the underscore: without it, it's not msdn official!
// Windows (x64 and x86)
#include <windows.h>   // required only for GetTickCount(...)
#include <process.h>   // _beginthread
#endif

#if __linux__
#include <unistd.h>
#include <time.h>

unsigned GetTickCount(void) {
    struct timespec ts;
    unsigned theTick = 0U;
    clock_gettime( CLOCK_REALTIME, &ts );
    theTick  = ts.tv_nsec / 1000000;
    theTick += ts.tv_sec * 1000;
    return theTick;
}
#endif

#if __APPLE__
#endif

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include "qra64.h"
#include "../qracodes/normrnd.h"		   // gaussian numbers generator

// ----------------------------------------------------------------------------

// channel types
#define CHANNEL_AWGN     0
#define CHANNEL_RAYLEIGH 1
#define CHANNEL_FASTFADE 2

#define JT65_SNR_EBNO_OFFSET 29.1f		// with the synch used in JT65
#define QRA64_SNR_EBNO_OFFSET 31.0f		// with the costas array synch

void printwordd(char *msg, int *x, int size)
{
  int k;
  printf("\n%s ",msg);
  for (k=0;k<size;k++)
    printf("%2d ",x[k]);
  printf("\n");
}
void printwordh(char *msg, int *x, int size)
{
  int k;
  printf("\n%s ",msg);
  for (k=0;k<size;k++)
    printf("%02hx ",x[k]);
  printf("\n");
}

#define NSAMPLES (QRA64_N*QRA64_M)

static float rp[NSAMPLES];
static float rq[NSAMPLES];
static float chp[NSAMPLES];
static float chq[NSAMPLES];
static float r[NSAMPLES];

float *mfskchannel(int *x, int channel_type, float EbNodB)
{
/*
Simulate an MFSK channel, either AWGN or Rayleigh.

x is a pointer to the transmitted codeword, an array of QRA64_N
integers in the range 0..63.

Returns the received symbol energies (squared amplitudes) as an array of
(QRA64_M*QRA64_N) floats.  The first QRA64_M entries of this array are
the energies of the first symbol in the codeword.  The second QRA64_M
entries are those of the second symbol, and so on up to the last codeword
symbol.
*/
  const float No = 1.0f;		        // noise spectral density
  const float sigma   = (float)sqrt(No/2.0f);	// std dev of noise I/Q components
  const float sigmach = (float)sqrt(1/2.0f);	// std dev of channel I/Q gains
  const float R = 1.0f*QRA64_K/QRA64_N;

  float EbNo = (float)pow(10,EbNodB/10);
  float EsNo = 1.0f*QRA64_m*R*EbNo;
  float Es = EsNo*No;
  float A = (float)sqrt(Es);
  int k;

  normrnd_s(rp,NSAMPLES,0,sigma);
  normrnd_s(rq,NSAMPLES,0,sigma);

  if(EbNodB>-15)
	if (channel_type == CHANNEL_AWGN)
	    for (k=0;k<QRA64_N;k++)
		  rp[k*QRA64_M+x[k]]+=A;
	else
		if (channel_type == CHANNEL_RAYLEIGH) {
			normrnd_s(chp,QRA64_N,0,sigmach);
			normrnd_s(chq,QRA64_N,0,sigmach);
			for (k=0;k<QRA64_N;k++) {
				rp[k*QRA64_M+x[k]]+=A*chp[k];
				rq[k*QRA64_M+x[k]]+=A*chq[k];
				}
			}
		else {
		return 0;	// unknown channel type
		}

  // compute the squares of the amplitudes of the received samples
  for (k=0;k<NSAMPLES;k++)
    r[k] = rp[k]*rp[k] + rq[k]*rq[k];

  return r;
}

// These defines are some packed fields as computed by JT65
#define CALL_IV3NWV		0x7F85AE7
#define CALL_K1JT		0xF70DDD7
#define GRID_JN66		0x3AE4		// JN66
#define GRID_73 		0x7ED0		// 73

char decode_type[12][32] = {
  "[?    ?    ?] AP0",
  "[CQ   ?    ?] AP27",
  "[CQ   ?     ] AP42",
  "[CALL ?    ?] AP29",
  "[CALL ?     ] AP44",
  "[CALL CALL ?] AP57",
  "[?    CALL ?] AP29",
  "[?    CALL  ] AP44",
  "[CALL CALL G] AP72",
  "[CQ   CALL ?] AP55",
  "[CQ   CALL  ] AP70",
  "[CQ   CALL G] AP70"
};
char apmode_type[3][32] = {
  "NO AP",
  "AUTO AP",
  "USER AP"
};

int test_proc_1(int channel_type, float EbNodB, int mode)
{
/*
Here we simulate the following (dummy) QSO:

1) CQ IV3NWV
2)                 IV3NWV K1JT
3) K1JT IV3NWV 73
4)                 IV3NWV K1JT 73

No message repetition is attempted

The QSO is counted as successfull if IV3NWV received the last message
When mode=QRA_AUTOAP each decoder attempts to decode the message sent
by the other station using the a-priori information derived by what
has been already decoded in a previous phase of the QSO if decoding
with no a-priori information has not been successful.

Step 1) K1JT's decoder first attempts to decode msgs of type [? ? ?]
and if this attempt fails, it attempts to decode [CQ/QRZ ? ?]  or
[CQ/QRZ ?] msgs

Step 2) if IV3NWV's decoder is unable to decode K1JT's without AP it
attempts to decode messages of the type [IV3NWV ? ?] and [IV3NWV ?].

Step 3) K1JT's decoder attempts to decode [? ? ?] and [K1JT IV3NWV ?]
(this last decode type has been enabled by K1JT's encoder at step 2)

Step 4) IV3NWV's decoder attempts to decode [? ? ?] and [IV3NWV K1JT
?] (this last decode type has been enabled by IV3NWV's encoder at step
3)

At each step the simulation reports if a decode was successful.  In
this case it also reports the type of decode (see table decode_type
above)

When mode=QRA_NOAP, only [? ? ?] decodes are attempted and no a-priori
information is used by the decoder

The function returns 0 if all of the four messages have been decoded
by their recipients (with no retries) and -1 if any of them could not
be decoded
*/

  int x[QRA64_K], xdec[QRA64_K];
  int y[QRA64_N];
  float *rx;
  int rc;

// Each simulated station must use its own codec since it might work with
// different a-priori information.
  qra64codec *codec_iv3nwv = qra64_init(mode);  // codec for IV3NWV
  qra64codec *codec_k1jt   = qra64_init(mode);    // codec for K1JT

// Step 1a: IV3NWV makes a CQ call (with no grid)
  printf("IV3NWV tx: CQ IV3NWV\n");
  encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_BLANK);
  qra64_encode(codec_iv3nwv, y, x);
  rx = mfskchannel(y,channel_type,EbNodB);

// Step 1b: K1JT attempts to decode [? ? ?], [CQ/QRZ ? ?] or [CQ/QRZ ?]
  rc = qra64_decode(codec_k1jt, 0, xdec,rx);
  if (rc>=0) { // decoded
    printf("K1JT   rx: received with apcode=%d %s\n",rc, decode_type[rc]);

// Step 2a: K1JT replies to IV3NWV (with no grid)
    printf("K1JT   tx: IV3NWV K1JT\n");
    encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_BLANK);
    qra64_encode(codec_k1jt, y, x);
    rx = mfskchannel(y,channel_type,EbNodB);

// Step 2b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?] or [IV3NWV ?]
    rc = qra64_decode(codec_iv3nwv, 0, xdec,rx);
    if (rc>=0) { // decoded
      printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]);

// Step 3a: IV3NWV replies to K1JT with a 73
      printf("IV3NWV tx: K1JT   IV3NWV 73\n");
      encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV, GRID_73);
      qra64_encode(codec_iv3nwv, y, x);
      rx = mfskchannel(y,channel_type,EbNodB);

// Step 3b: K1JT attempts to decode [? ? ?] or [K1JT IV3NWV ?]
      rc = qra64_decode(codec_k1jt, 0, xdec,rx);
      if (rc>=0) { // decoded
	printf("K1JT   rx: received with apcode=%d %s\n",rc, decode_type[rc]);

// Step 4a: K1JT replies to IV3NWV with a 73
	printf("K1JT   tx: IV3NWV K1JT   73\n");
	encodemsg_jt65(x,CALL_IV3NWV,CALL_K1JT, GRID_73);
	qra64_encode(codec_k1jt, y, x);
	rx = mfskchannel(y,channel_type,EbNodB);

// Step 4b: IV3NWV attempts to decode [? ? ?], [IV3NWV ? ?], or [IV3NWV ?]
	rc = qra64_decode(codec_iv3nwv, 0, xdec,rx);
	if (rc>=0) { // decoded
	  printf("IV3NWV rx: received with apcode=%d %s\n",rc, decode_type[rc]);
	  return 0;
	}
      }
    }
  }
  printf("no decode\n");
  return -1;
}

int test_proc_2(int channel_type, float EbNodB, int mode)
{
/*
Here we simulate the decoder of K1JT after K1JT has sent a msg [IV3NWV K1JT]
and IV3NWV sends him the msg [K1JT IV3NWV JN66].

If mode=QRA_NOAP, K1JT decoder attempts to decode only msgs of type [? ? ?].

If mode=QRA_AUTOP, K1JT decoder will attempt to decode also the msgs
[K1JT IV3NWV] and [K1JT IV3NWV ?].

In the case a decode is successful the return code of the qra64_decode function
indicates the amount of a-priori information required to decode the received
message according to this table:

 rc=0    [?    ?    ?]		AP0
 rc=1    [CQ   ?    ?]		AP27
 rc=2    [CQ   ?     ]		AP42
 rc=3    [CALL ?    ?]		AP29
 rc=4    [CALL ?     ]		AP44
 rc=5    [CALL CALL ?]		AP57
 rc=6    [?    CALL ?]		AP29
 rc=7    [?    CALL  ]		AP44
 rc=8    [CALL CALL GRID]	AP72
 rc=9    [CQ   CALL ?]		AP55
 rc=10   [CQ   CALL  ]		AP70
 rc=11   [CQ   CALL GRID]	AP70

The return code is <0 when decoding is unsuccessful

This test simulates the situation ntx times and reports how many times
a particular type decode among the above 6 cases succeded.
*/

  int x[QRA64_K], xdec[QRA64_K];
  int y[QRA64_N];
  float *rx;
  float ebnodbest, ebnodbavg=0;
  int rc,k;

  int ndecok[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  int nundet = 0;
  int ntx = 200,ndec=0;

  qra64codec *codec_iv3nwv = qra64_init(mode);   // codec for IV3NWV
  qra64codec *codec_k1jt   = qra64_init(mode);     // codec for K1JT

  printf("\nQRA64 Test #2 - Decoding with AP knowledge (SNR-Eb/No offset = %.1f dB)\n\n",
		   QRA64_SNR_EBNO_OFFSET);

// This will enable K1JT's decoder to look for calls directed to him [K1JT ? ?/b]
//  printf("K1JT decoder enabled for [K1JT  ?     ?/blank]\n");
//  qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_MYCALL);

// This will enable K1JT's decoder to look for IV3NWV calls directed to him [K1JT IV3NWV ?/b]
//  printf("K1JT decoder enabled for [K1JT IV3NWV ?]\n");
//  qra64_apset(codec_k1jt, CALL_CQ,CALL_IV3NWV,0,APTYPE_BOTHCALLS);

// This will enable K1JT's decoder to look for msges sent by IV3NWV [? IV3NWV ?]
//  printf("K1JT decoder enabled for [?    IV3NWV ?/blank]\n");
//  qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_BLANK,APTYPE_HISCALL);

// This will enable K1JT's decoder to look for full-knowledge [K1JT IV3NWV JN66] msgs
  printf("K1JT decoder enabled for [K1JT IV3NWV JN66]\n");
  qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,GRID_JN66,APTYPE_FULL);

// This will enable K1JT's decoder to look for calls from IV3NWV [CQ IV3NWV ?/b] msgs
  printf("K1JT decoder enabled for [CQ   IV3NWV ?/b/JN66]\n");
  qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_JN66,APTYPE_CQHISCALL);

  // Dx station IV3NWV calls
  printf("\nIV3NWV encoder sends msg: [K1JT IV3NWV JN66]\n\n");
  encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_JN66);

//  printf("\nIV3NWV encoder sends msg: [CQ IV3NWV JN66]\n\n");
//  encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_JN66);

//  printf("\nIV3NWV encoder sends msg: [CQ IV3NWV]\n\n");
//  encodemsg_jt65(x,CALL_CQ,CALL_IV3NWV,GRID_BLANK);
  qra64_encode(codec_iv3nwv, y, x);

  printf("Simulating K1JT decoder up to AP72\n");

  for (k=0;k<ntx;k++) {
    printf(".");
    rx = mfskchannel(y,channel_type,EbNodB);
    rc = qra64_decode(codec_k1jt, &ebnodbest, xdec,rx);
	if (rc>=0) {
	  ebnodbavg +=ebnodbest;
	  if (memcmp(xdec,x,12*sizeof(int))==0)
		ndecok[rc]++;
	  else
	    nundet++;
	}
  }
  printf("\n\n");

  printf("Transimtted msgs:%d\nDecoded msgs:\n\n",ntx);
  for (k=0;k<12;k++) {
    printf("%3d with %s\n",ndecok[k],decode_type[k]);
    ndec += ndecok[k];
  }
  printf("\nTotal: %d/%d (%d undetected errors)\n\n",ndec,ntx,nundet);
  printf("");

  ebnodbavg/=(ndec+nundet);
  printf("Estimated SNR (average in dB) = %.2f dB\n\n",ebnodbavg-QRA64_SNR_EBNO_OFFSET);

  return 0;
}

int test_fastfading(float EbNodB, float B90, int fadingModel, int submode, int apmode, int olddec, int channel_type, int ntx)
{
  int x[QRA64_K], xdec[QRA64_K];
  int y[QRA64_N];
  float *rx;
  float ebnodbest, ebnodbavg=0;
  int rc,k;
  float rxolddec[QRA64_N*QRA64_M];	// holds the energies at nominal tone freqs

  int ndecok[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  int nundet = 0;
  int ndec=0;

  qra64codec *codec_iv3nwv;
  qra64codec *codec_k1jt;

  codec_iv3nwv=qra64_init(QRA_NOAP);
  codec_k1jt  =qra64_init(apmode);

  if (channel_type==2) {	// fast-fading case
	printf("Simulating the fast-fading channel\n");
	printf("B90=%.2f Hz - Fading Model=%s - Submode=QRA64%c\n",B90,fadingModel?"Lorentz":"Gauss",submode+'A');
	printf("Decoder metric = %s\n",olddec?"AWGN":"Matched to fast-fading signal");
	}
  else {
    printf("Simulating the %s channel\n",channel_type?"Rayleigh block fading":"AWGN");
	printf("Decoder metric = AWGN\n");
	}

  printf("\nEncoding msg [K1JT IV3NWV JN66]\n");
  encodemsg_jt65(x,CALL_K1JT,CALL_IV3NWV,GRID_JN66);
// printf("[");
//  for (k=0;k<11;k++) printf("%02hX ",x[k]); printf("%02hX]\n",x[11]);

  qra64_encode(codec_iv3nwv, y, x);
  printf("%d transmissions will be simulated\n\n",ntx);

  if (apmode==QRA_USERAP) {
	// This will enable K1JT's decoder to look for cq/qrz calls [CQ/QRZ ? ?/b]
	printf("K1JT decoder enabled for [CQ    ?     ?/blank]\n");
	qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_CQQRZ);

	// This will enable K1JT's decoder to look for calls directed to him [K1JT ? ?/b]
	printf("K1JT decoder enabled for [K1JT  ?     ?/blank]\n");
	qra64_apset(codec_k1jt, CALL_K1JT,0,0,APTYPE_MYCALL);

	// This will enable K1JT's decoder to look for msges sent by IV3NWV [? IV3NWV ?]
	printf("K1JT decoder enabled for [?    IV3NWV ?/blank]\n");
	qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_BLANK,APTYPE_HISCALL);

	// This will enable K1JT's decoder to look for IV3NWV calls directed to him [K1JT IV3NWV ?/b]
	printf("K1JT decoder enabled for [K1JT IV3NWV ?]\n");
	qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,0,APTYPE_BOTHCALLS);

	// This will enable K1JT's decoder to look for full-knowledge [K1JT IV3NWV JN66] msgs
	printf("K1JT decoder enabled for [K1JT IV3NWV JN66]\n");
	qra64_apset(codec_k1jt, CALL_K1JT,CALL_IV3NWV,GRID_JN66,APTYPE_FULL);

	// This will enable K1JT's decoder to look for calls from IV3NWV [CQ IV3NWV ?/b] msgs
	printf("K1JT decoder enabled for [CQ   IV3NWV ?/b/JN66]\n");
	qra64_apset(codec_k1jt, 0,CALL_IV3NWV,GRID_JN66,APTYPE_CQHISCALL);

  }

  printf("\nNow decoding with K1JT's decoder...\n");
/*
  if (channel_type==2) 	// simulate a fast-faded signal
	  printf("Simulating a fast-fading channel with given B90 and spread type\n");
  else
	  printf("Simulating a %s channel\n",channel_type?"Rayleigh block fading":"AWGN");
*/
  for (k=0;k<ntx;k++) {

	  if ((k%10)==0)
	    printf("  %5.1f %%\r",100.0*k/ntx);
//		printf(".");	// work in progress

	if (channel_type==2) {
		// generate a fast-faded signal
		rc = qra64_fastfading_channel(&rx,y,submode,EbNodB,B90,fadingModel);
		if (rc<0) {
			printf("\nqra64_fastfading_channel error. rc=%d\n",rc);
			return -1;
			}
		}
	else // generate a awgn or Rayleigh block fading signal
		rx = mfskchannel(y, channel_type, EbNodB);

	if (channel_type==2)	// fast-fading case
		if (olddec==1) {
			int k, j;
			int jj = 1<<submode;
			int bps = QRA64_M*(2+jj);
			float *rxbase;
			float *out = rxolddec;
			// calc energies at nominal freqs
			for (k=0;k<QRA64_N;k++) {
				rxbase = rx + QRA64_M + k*bps;
				for (j=0;j<QRA64_M;j++) {
					*out++=*rxbase;
					rxbase+=jj;
					}
				}
			// decode with awgn decoder
			rc = qra64_decode(codec_k1jt,&ebnodbest,xdec,rxolddec);
			}
		else // use fast-fading decoder
			rc = qra64_decode_fastfading(codec_k1jt,&ebnodbest,xdec,rx,submode,B90,fadingModel);
	else // awgn or rayleigh channel. use the old decoder whatever the olddec option is
		rc = qra64_decode(codec_k1jt,&ebnodbest,xdec,rx);

	if (rc>=0) {
	  ebnodbavg +=ebnodbest;
	  if (memcmp(xdec,x,12*sizeof(int))==0)
		ndecok[rc]++;
	  else {
		fprintf(stderr,"\nUndetected error with rc=%d\n",rc);
		nundet++;
		}
	  }

  }
  printf("  %5.1f %%\r",100.0*k/ntx);

  printf("\n\n");

  printf("Msgs transmitted:%d\nMsg decoded:\n\n",ntx);
  for (k=0;k<12;k++) {
    printf("rc=%2d   %3d with %s\n",k,ndecok[k],decode_type[k]);
    ndec += ndecok[k];
  }
  printf("\nTotal: %d/%d (%d undetected errors)\n\n",ndec,ntx,nundet);
  printf("");

  if (ndec>0) {
	ebnodbavg/=(ndec+nundet);
	printf("Estimated SNR (average in dB) = %.2f dB\n\n",ebnodbavg-QRA64_SNR_EBNO_OFFSET);
	}

  return 0;
}

void syntax(void)
{

  printf("\nQRA64 Mode Tests\n");
  printf("2016, Nico Palermo - IV3NWV\n\n");
  printf("---------------------------\n\n");
  printf("Syntax: qra64 [-s<snrdb>] [-c<channel>] [-a<ap-type>] [-t<testtype>] [-h]\n");
  printf("Options: \n");
  printf("       -s<snrdb>   : set simulation SNR in 2500 Hz BW (default:-27.5 dB)\n");
  printf("       -c<channel> : set channel type 0=AWGN (default) 1=Rayleigh 2=Fast-fading\n");
  printf("       -a<ap-type> : set decode type 0=NOAP 1=AUTOAP (default) 2=USERAP\n");
  printf("       -t<testtype>: 0=simulate seq of msgs between IV3NWV and K1JT (default)\n");
  printf("                     1=simulate K1JT receiving K1JT IV3NWV JN66\n");
  printf("                     2=simulate fast-fading/awgn/rayliegh decoders performance\n");
  printf("       -n<ntx>     : simulate the transmission of ntx codewords (default=100)\n");

  printf("Options used only for fast-fading simulations (-c2):\n");
  printf("       -b          : 90%% fading bandwidth in Hz [1..230 Hz] (default = 2.5 Hz)\n");
  printf("       -m          : fading model. 0=Gauss, 1=Lorentz (default = Lorentz)\n");
  printf("       -q          : qra64 submode. 0=QRA64A,... 4=QRA64E (default = QRA64A)\n");
  printf("       -d          : use the old awgn decoder\n");
  printf("       -h: this help\n");
  printf("Example:\n");
  printf("        qra64 -t2 -c2 -a2 -b50 -m1 -q2 -n10000 -s-26\n");
  printf("        runs the error performance test (-t2)\n");
  printf("        with USER_AP (-a2)\n");
  printf("        simulating a fast fading channel (-c2)\n");
  printf("        with B90 = 50 Hz (-b50), Lorentz Doppler (-m1), mode QRA64C (-q2)\n");
  printf("        ntx = 10000 codewords (-n10000) and SNR = -26 dB (-s-26)\n");

}

int main(int argc, char* argv[])
{
  int k, rc, nok=0;
  float SNRdB = -27.5f;
  unsigned int channel = CHANNEL_AWGN;
  unsigned int mode    = QRA_AUTOAP;
  unsigned int testtype=0;
  int   nqso = 100;
  float EbNodB;
  float B90 = 2.5;
  int fadingModel = 1;
  int submode = 0;
  int olddec = 0;
  int ntx = 100;

// Parse the command line
  while(--argc) {
    argv++;

    if (strncmp(*argv,"-h",2)==0) {
      syntax();
      return 0;
	  }
	else
	if (strncmp(*argv,"-n",2)==0) {
		ntx = ( int)atoi((*argv)+2);
		if (ntx<100 || ntx>1000000) {
			printf("Invalid -n option. ntx must be in the range [100..1000000]\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-a",2)==0) {
		mode = ( int)atoi((*argv)+2);
		if (mode>2) {
			printf("Invalid decoding mode\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-s",2)==0) {
		SNRdB = (float)atof((*argv)+2);
		if (SNRdB>20 || SNRdB<-50) {
			printf("SNR should be in the range [-50..20]\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-t",2)==0) {
	    testtype = ( int)atoi((*argv)+2);
	    if (testtype>2) {
			printf("Invalid test type\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-c",2)==0) {
		channel = ( int)atoi((*argv)+2);
	    if (channel>CHANNEL_FASTFADE) {
			printf("Invalid channel type\n");
			syntax();
			return -1;
			}
	    }
	else
	if (strncmp(*argv,"-b",2)==0) {
	    B90 = (float)atof((*argv)+2);
	    if (B90<1 || B90>230) {
			printf("Invalid B90\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-m",2)==0) {
	    fadingModel = (int)atoi((*argv)+2);
	    if (fadingModel<0 || fadingModel>1) {
			printf("Invalid fading model\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-q",2)==0) {
	    submode = (int)atoi((*argv)+2);
	    if (submode<0 || submode>4) {
			printf("Invalid submode\n");
			syntax();
			return -1;
			}
		}
	else
	if (strncmp(*argv,"-d",2)==0) {
		olddec = 1;
		}
	else {
	    printf("Invalid option\n");
	    syntax();
	    return -1;
	    }
	}

  if (testtype<2) // old tests
	  if (channel==CHANNEL_FASTFADE) {
		printf("Invalid Option. Test type 0 and 1 supports only AWGN or Rayleigh Channel model\n");
		return -1;
		}

  EbNodB = SNRdB+QRA64_SNR_EBNO_OFFSET;

#if defined(__linux__) || defined(__unix__)
  srand48(GetTickCount());
#endif

  if (testtype==0) {
	for (k=0;k<nqso;k++) {
		printf("\n\n------------------------\n");
		rc = test_proc_1(channel, EbNodB, mode);
		if (rc==0)
			nok++;
		}
	printf("\n\n%d/%d QSOs to end without repetitions\n",nok,nqso);
	printf("Input SNR = %.1fdB channel=%s ap-mode=%s\n\n",
		SNRdB,
		channel==CHANNEL_AWGN?"AWGN":"RAYLEIGH",
		apmode_type[mode]
		);
	}
  else if (testtype==1) {
	test_proc_2(channel, EbNodB, mode);
	printf("Input SNR = %.1fdB channel=%s ap-mode=%s\n\n",
		SNRdB,
		channel==CHANNEL_AWGN?"AWGN":"RAYLEIGH",
		apmode_type[mode]
		);
	}
  else {
	printf("Input SNR = %.1fdB ap-mode=%s\n\n",
	 SNRdB,
	 apmode_type[mode]
	 );
	test_fastfading(EbNodB,B90,fadingModel,submode,mode,olddec, channel, ntx);
  }
  return 0;
}