code:
/* Initialize a RS codec
 *
 * Copyright 2002 Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#include <stdlib.h>

#ifdef CCSDS
#include "ccsds.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

void FREE_RS(void *p){
  struct rs *rs = (struct rs *)p;

  free(rs->alpha_to);
  free(rs->index_of);
  free(rs->genpoly);
  free(rs);
}

/* Initialize a Reed-Solomon codec
 * symsize = symbol size, bits (1-8)
 * gfpoly = Field generator polynomial coefficients
 * fcr = first root of RS code generator polynomial, index form
 * prim = primitive element to generate polynomial roots
 * nroots = RS code generator polynomial degree (number of roots)
 */
void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned fcr,unsigned prim,
		unsigned int nroots){
  struct rs *rs;
  int i, j, sr,root,iprim;

  if(symsize > 8*sizeof(DTYPE))
    return NULL; /* Need version with ints rather than chars */

  if(fcr >= (1<<symsize))
    return NULL;
  if(prim == 0 || prim >= (1<<symsize))
    return NULL;
  if(nroots >= (1<<symsize))
    return NULL; /* Can't have more roots than symbol values! */

  rs = (struct rs *)calloc(1,sizeof(struct rs));
  rs->mm = symsize;
  rs->nn = (1<<symsize)-1;

  rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
  if(rs->alpha_to == NULL){
    free(rs);
    return NULL;
  }
  rs->index_of = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
  if(rs->index_of == NULL){
    free(rs->alpha_to);
    free(rs);
    return NULL;
  }

  /* Generate Galois field lookup tables */
  rs->index_of[0] = A0; /* log(zero) = -inf */
  rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */
  sr = 1;
  for(i=0;i<rs->nn;i++){
    rs->index_of[sr] = i;
    rs->alpha_to[i] = sr;
    sr <<= 1;
    if(sr & (1<<symsize))
      sr ^= gfpoly;
    sr &= rs->nn;
  }
  if(sr != 1){
    /* field generator polynomial is not primitive! */
    free(rs->alpha_to);
    free(rs->index_of);
    free(rs);
    return NULL;
  }

  /* Form RS code generator polynomial from its roots */
  rs->genpoly = (DTYPE *)malloc(sizeof(DTYPE)*(nroots+1));
  if(rs->genpoly == NULL){
    free(rs->alpha_to);
    free(rs->index_of);
    free(rs);
    return NULL;
  }
  rs->fcr = fcr;
  rs->prim = prim;
  rs->nroots = nroots;

  /* Find prim-th root of 1, used in decoding */
  for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
    ;
  rs->iprim = iprim / prim;

  rs->genpoly[0] = 1;
  for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) {
    rs->genpoly[i+1] = 1;

    /* Multiply rs->genpoly[] by  @**(root + x) */
    for (j = i; j > 0; j--){
      if (rs->genpoly[j] != 0)
	rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)];
      else
	rs->genpoly[j] = rs->genpoly[j-1];
    }
    /* rs->genpoly[0] can never be zero */
    rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)];
  }
  /* convert rs->genpoly[] to index form for quicker encoding */
  for (i = 0; i <= nroots; i++)
    rs->genpoly[i] = rs->index_of[rs->genpoly[i]];

  return rs;
}


code:
/* Include file to configure the RS codec for integer symbols
 *
 * Copyright 2002, Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#define DTYPE int

/* Reed-Solomon codec control block */
struct rs {
  unsigned int mm;   /* Bits per symbol */
  unsigned int nn;   /* Symbols per block (= (1<<mm)-1) */
  int *alpha_to;      /* log lookup table */
  int *index_of;      /* Antilog lookup table */
  int *genpoly;       /* Generator polynomial */
  unsigned int nroots;     /* Number of generator roots = number of parity symbols */
  unsigned int fcr;        /* First consecutive root, index form */
  unsigned int prim;       /* Primitive element, index form */
  unsigned int iprim;      /* prim-th root of 1, index form */
};

static inline int modnn(struct rs *rs,int x){
  while (x >= rs->nn) {
    x -= rs->nn;
    x = (x >> rs->mm) + (x & rs->nn);
  }
  return x;
}
#define MODNN(x) modnn(rs,x)

#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)
#define FCR (rs->fcr)
#define PRIM (rs->prim)
#define IPRIM (rs->iprim)
#define A0 (NN)

#define ENCODE_RS encode_rs_int
#define DECODE_RS decode_rs_int
#define INIT_RS init_rs_int
#define FREE_RS free_rs_int

void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras, int calc_syn);
void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned int fcr,
		   unsigned int prim,unsigned int nroots);
void FREE_RS(void *p);

license: mit


code:
/* User include file for the Reed-Solomon codec
 * Copyright 2002, Phil Karn KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */

/* General purpose RS codec, integer symbols */
void encode_rs_int(void *rs,int *data,int *parity);
int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras, int calc_syn);
void *init_rs_int(int symsize,int gfpoly,int fcr,
		  int prim,int nroots,int pad);
void free_rs_int(void *rs);

/* Tables to map from conventional->dual (Taltab) and
 * dual->conventional (Tal1tab) bases
 */
extern unsigned char Taltab[],Tal1tab[];


code:
program rsdtest

  real s3(64,63)
  character msg*22,arg*12
  integer param(0:7)

  nargs=iargc()
  if(nargs.ne.2) then
     print*,'Usage: rsdtest ntrials nfiles'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) ntrials
  call getarg(2,arg)
  read(arg,*) nfiles

  open(10,file='s3_1000.bin',access='stream', status='old')
  open(22,file='kvasd.dat',access='direct',recl=1024,status='unknown')

  nadd=1
  ifile0=0
  if(nfiles.lt.0) then
     ifile0=-nfiles
     nfiles=99999
  endif

  do ifile=1,nfiles
     read(10,end=999) s3
     if(ifile.lt.ifile0) cycle
     call extract2(s3,nadd,ntrials,param,msg)
     if(ifile.eq.ifile0) exit
  enddo

999 end program rsdtest


code:
/*
 sfrsd.c

 A soft-decision decoder for the JT65 (63,12) Reed-Solomon code.

 This decoding scheme is built around Phil Karn's Berlekamp-Massey
 errors and erasures decoder. The approach is inspired by a number of
 publications, including the stochastic Chase decoder described
 in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
 IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
 "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
 and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.

 Steve Franke K9AN, Urbana IL, September 2015
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <string.h>
#include "sfrsd2.h"

//***************************************************************************
void usage(void)
{
    printf("Usage: sfrsd [options...] <path to kvasd.dat>\n");
    printf("       input file should be in kvasd format\n");
    printf("\n");
    printf("Options:\n");
    printf("       -n number of random erasure vectors to try\n");
    printf("       -v verbose\n");
}

int main(int argc, char *argv[]){

    extern char *optarg;
    extern int optind;

    int correct[63], indx[63], param[8];
    int c,i;
    char *infile;

    FILE *datfile, *logfile;
    int nsec, maxe, nads;
    float xlambda;
    int mrsym[63],mrprob[63],mr2sym[63],mr2prob[63];
    int nsec2,ncount,dat4[12];
    int ntrials, nverbose, ntry;
    int nhard;
    double tt;

    ntrials=10000;
    nverbose=1;

    while ( (c = getopt(argc, argv, "n:qv")) !=-1 ) {
        switch (c) {
                case 'n':
                ntrials=(int)strtof(optarg,NULL);
                printf("ntrials set to %d\n",ntrials);
                break;
                case 'v':
                nverbose=1;
                break;
                case 'q': //accept (and ignore) -q option for WSJT10 compatibility
                break;
                case '?':
                usage();
                exit(1);
        }
    }

    if( optind+1 > argc) {
        //        usage();
        //        exit(1);
        infile="kvasd.dat";
    } else {
        infile=argv[optind];
    }

    logfile=fopen("/tmp/sfrsd.log","a");
    if( !logfile ) {
        printf("Unable to open sfrsd.log\n");
        exit(1);
    }

    datfile=fopen(infile,"rb");
    if( !datfile ) {
        printf("Unable to open kvasd.dat\n");
        exit(1);
    } else {
        fread(&nsec,sizeof(int),1,datfile);
        fread(&xlambda,sizeof(float),1,datfile);
        fread(&maxe,sizeof(int),1,datfile);
        fread(&nads,sizeof(int),1,datfile);
        fread(&mrsym,sizeof(int),63,datfile);
        fread(&mrprob,sizeof(int),63,datfile);
        fread(&mr2sym,sizeof(int),63,datfile);
        fread(&mr2prob,sizeof(int),63,datfile);
        fread(&nsec2,sizeof(int),1,datfile);
        fread(&ncount,sizeof(int),1,datfile);
        fread(&dat4,sizeof(int),12,datfile);
        fclose(datfile);
    }

    sfrsd2_(mrsym,mrprob,mr2sym,mr2prob,&ntrials,&nverbose,correct,param,indx,&tt,&ntry);
    nhard=param[1];
    if( nhard>=0 ) {
        for (i=0; i<12; i++) {
            dat4[i]=correct[11-i];
        }
    } else {
        nhard=-1;
        memset(dat4,0,12*sizeof(int));
    }
    datfile=fopen(infile,"wb");
    if( !datfile ) {
        printf("Unable to open kvasd.dat\n");
        return 1;
    } else {
        fwrite(&nsec,sizeof(int),1,datfile);
        fwrite(&xlambda,sizeof(float),1,datfile);
        fwrite(&maxe,sizeof(int),1,datfile);
        fwrite(&nads,sizeof(int),1,datfile);
        fwrite(&mrsym,sizeof(int),63,datfile);
        fwrite(&mrprob,sizeof(int),63,datfile);
        fwrite(&mr2sym,sizeof(int),63,datfile);
        fwrite(&mr2prob,sizeof(int),63,datfile);
        fwrite(&nsec2,sizeof(int),1,datfile);
        fwrite(&nhard,sizeof(int),1,datfile);
        fwrite(&dat4,sizeof(int),12,datfile);
        fclose(datfile);
    }
    exit(0);
}


code:
void ftrsd2_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[],
	     int* ntrials0, int* verbose0, int correct[], int param[],
	     int indexes[], double tt[], int ntry[]);


code:
/*
 sfrsd2.c

 A soft-decision decoder for the JT65 (63,12) Reed-Solomon code.

 This decoding scheme is built around Phil Karn's Berlekamp-Massey
 errors and erasures decoder. The approach is inspired by a number of
 publications, including the stochastic Chase decoder described
 in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
 IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
 "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
 and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.

 Steve Franke K9AN and Joe Taylor K1JT
 */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <string.h>
#include "rs2.h"

static void *rs;

void sfrsd2_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[],
	     int* ntrials0, int* verbose0, int correct[], int param[],
	     int indexes[], double tt[], int ntry[])
{
  int rxdat[63], rxprob[63], rxdat2[63], rxprob2[63];
  int workdat[63],workdat2[63];
  int era_pos[51];
  int c, i, j, numera, nmr2, nerr, nn=63, kk=12;
  FILE *datfile, *logfile;
  int ntrials = *ntrials0;
  int verbose = *verbose0;
  int nhard=0,nhard_min=32768,nsoft=0,nsoft_min=32768, ncandidates;
  int ngmd,nera_best;
  clock_t t0=0,t1=0;
  int perr[8][8] = {
     12,     31,     44,     52,     60,     57,     50,     50,
     28,     38,     49,     58,     65,     69,     64,     80,
     40,     41,     53,     62,     66,     73,     76,     81,
     50,     53,     53,     64,     70,     76,     77,     81,
     50,     50,     52,     60,     71,     72,     77,     84,
     50,     50,     56,     62,     67,     73,     81,     85,
     50,     50,     71,     62,     70,     77,     80,     85,
     50,     50,     62,     64,     71,     75,     82,     87};

  int pmr2[8][8] = {
      4,      8,      9,      7,      6,      0,      0,      0,
     13,     18,     15,     11,      9,      7,      5,      0,
      0,     23,     21,     15,     12,     10,      7,      4,
      0,     34,     28,     20,     16,     14,     11,      7,
      0,     20,     26,     25,     19,     14,     12,      9,
      0,      0,     28,     27,     22,     19,     14,     11,
      0,      0,     40,     29,     29,     23,     18,     12,
      0,      0,     40,     35,     31,     21,     20,     13};

  if(verbose) {
    logfile=fopen("sfrsd.log","a");
    if( !logfile ) {
      printf("Unable to open sfrsd.log\n");
      exit(1);
    }
  }

// Initialize the KA9Q Reed-Solomon encoder/decoder
  unsigned int symsize=6, gfpoly=0x43, fcr=3, prim=1, nroots=51;
  rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 0);

// Reverse the received symbol vector for BM decoder
  for (i=0; i<63; i++) {
    rxdat[i]=mrsym[62-i];
    rxprob[i]=mrprob[62-i];
    rxdat2[i]=mr2sym[62-i];
    rxprob2[i]=mr2prob[62-i];
  }

// Sort the mrsym probabilities to find the least reliable symbols
  int k, pass, tmp, nsym=63;
  int probs[63];
  for (i=0; i<63; i++) {
    indexes[i]=i;
    probs[i]=rxprob[i];
  }
  for (pass = 1; pass <= nsym-1; pass++) {
    for (k = 0; k < nsym - pass; k++) {
      if( probs[k] < probs[k+1] ) {
	tmp = probs[k];
	probs[k] = probs[k+1];
	probs[k+1] = tmp;
	tmp = indexes[k];
	indexes[k] = indexes[k+1];
	indexes[k+1] = tmp;
      }
    }
  }

// See if we can decode using BM HDD, and calculate the syndrome vector.
  memset(era_pos,0,51*sizeof(int));
  numera=0;
  memcpy(workdat,rxdat,sizeof(rxdat));
  nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
  if( nerr >= 0 ) {
    if(verbose) fprintf(logfile,"   BM decode nerrors= %3d : ",nerr);
    memcpy(correct,workdat,63*sizeof(int));
    ngmd=-1;
    param[0]=0;
    param[1]=0;
    param[2]=0;
    param[3]=0;
    param[4]=0;
    return;
  }

/*
Generate random erasure-locator vectors and see if any of them
decode. This will generate a list of potential codewords. The
"soft" distance between each codeword and the received word is
used to decide which codeword is "best".
*/

#ifdef WIN32
  srand(0xdeadbeef);
#else
  srandom(0xdeadbeef);
#endif

  float ratio, ratio0[63];
  int threshe, thresh2, nsum;
  int thresh0[63],thresh1[63], mr2flag;
  ncandidates=0;
  nsum=0;
  int ii,jj;
  for (i=0; i<nn; i++) {
    nsum=nsum+rxprob[i];
    j = indexes[62-i];
    ratio = (float)rxprob2[j]/(float)rxprob[j];
    ratio0[i]=ratio;
    ii = 7.999*ratio;
    jj = (62-i)/8;
    thresh0[i] = 1.3*perr[ii][jj];
    thresh1[i] = 0.4*pmr2[ii][jj];
  }
  if(nsum==0) return;

  for( k=0; k<ntrials; k++) {
    memset(era_pos,0,51*sizeof(int));
    memcpy(workdat,rxdat,sizeof(rxdat));

/*
Mark a subset of the symbols as erasures.
Run through the ranked symbols, starting with the worst, i=0.
NB: j is the symbol-vector index of the symbol with rank i.
*/
    numera=0;
    nmr2=0;
    for (i=0; i<nn; i++) {
      j = indexes[62-i];
      threshe=thresh0[i];
      thresh2=thresh1[i];
      long int ir, ir2;
#ifdef WIN32
      ir=rand();
      ir2=rand();
#else
      ir=random();
      ir2=random();
#endif
      if( ((ir % 100) < threshe ) && (numera+2*nmr2) < 51 ) {
          era_pos[numera]=j;
	  numera=numera+1;
      }
      if( ((ir2 % 100) < thresh2) && (numera+2*nmr2)<51) {
          workdat[j]=rxdat2[j];
          nmr2=nmr2+1;
      }
    }
    t0=clock();
//    rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 1);
    nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
    t1=clock();
    tt[0]+=(double)(t1-t0)/CLOCKS_PER_SEC;

    if( nerr >= 0 ) {
      ncandidates=ncandidates+1;
      nhard=0;
      nsoft=0;
      for (i=0; i<63; i++) {
	if(workdat[i] != rxdat[i]) {
	  nhard=nhard+1;
	  if(workdat[i] != rxdat2[i]) {
	    nsoft=nsoft+rxprob[i];
	  }
	}
      }
      nsoft=63*nsoft/nsum;
      if((nsoft < 33) && (nhard < 43) && (nhard+nsoft) < 74) {  //???
	if( (nsoft < nsoft_min) ) {
	  nsoft_min=nsoft;
	  nhard_min=nhard;
	  memcpy(correct,workdat,63*sizeof(int));
	  ngmd=0;
	  nera_best=numera;
	  ntry[0]=k;
	}
      }
      if(nsoft_min < 27) break;
      if((nsoft_min < 32) && (nhard_min < 43) &&
	 (nhard_min+nsoft_min) < 74) break;
    }
    if(k == ntrials-1) ntry[0]=k+1;
  }

  if(verbose) fprintf(logfile,
     "%d trials and %d candidates after stochastic loop\n",k,ncandidates);

  if( (ncandidates >= 0) && (nsoft_min < 36) && (nhard_min < 44) ) {
    if(verbose) {
      for (i=0; i<63; i++) {
	fprintf(logfile,"%3d %3d %3d %3d %3d %3d\n",i,correct[i],
		rxdat[i],rxprob[i],rxdat2[i],rxprob2[i]);
      }
      fprintf(logfile,"**** ncandidates %d nhard %d nsoft %d nsum %d\n",
	      ncandidates,nhard_min,nsoft_min,nsum);
    }
  } else {
    nhard_min=-1;
  }

  if(verbose) {
    fprintf(logfile,"exiting sfrsd\n");
    fclose(logfile);
  }
  param[0]=ncandidates;
  param[1]=nhard_min;
  param[2]=nsoft_min;
  param[3]=nera_best;
  param[4]=ngmd;
  if(param[0]==0) param[2]=-1;
  return;
}


code:
gcc -c gran.c
gfortran -c four2a.f90
gfortran -c f77_wisdom.f90
gfortran -o chkfft chkfft.f90 four2a.o f77_wisdom.o gran.o -lfftw3f


code:
gcc -c gran.c
gfortran -c four2a.f90
gfortran -c f77_wisdom.f90
gfortran -o chkfft chkfft.f90 four2a.o f77_wisdom.o gran.o /JTSDK-QT/appsupport/runtime/libfftw3f-3.dll


code:
gcc -c gran.c
gfortran -c fftw3mod.f90
gfortran -c f77_wisdom.f90
gfortran -o chkfft2 chkfft2.f90 f77_wisdom.o gran.o /JTSDK-QT/appsupport/runtime/libfftw3f-3.dll


code:
gcc -c gran.c
gfortran -c fftw3mod.f90
gfortran -c f77_wisdom.f90
gfortran -o chkfft3 chkfft3.f90 f77_wisdom.o gran.o /JTSDK-QT/appsupport/runtime/libfftw3f-3.dll


code:
gcc -c gran.c
gfortran -c -Wall fftw3mod.f90
gfortran -o timefft -Wall timefft.f90 timefft_opts.f90 gran.o libfftw3f-3.dll


code:
subroutine gen4(msg0,ichk,msgsent,itone,itype)

! Encode a JT4 message.  Returns msgsent, the message as it will be
! decoded, an integer array itone(206) of 4-FSK tons values in the
! range 0-3; and itype, the JT message type.

  use jt4
  use packjt
  character*22 msg0
  character*22 message          !Message to be generated
  character*22 msgsent          !Message as it will be received
  character*1 c
  integer itone(206)
  integer*4 i4Msg6BitWords(13)            !72-bit message as 6-bit words
  integer mettab(-128:127,0:1)
  save

  if(msg0(1:1).eq.'@') then
     read(msg0(2:5),*,end=1,err=1) nfreq
     go to 2
1    nfreq=1000
2    itone(1)=nfreq
     msgsent=msg0
  else
     call getmet4(mettab,ndelta)

     message=msg0
     call fmtmsg(message,iz)
     call packmsg(message,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
     call unpackmsg(i4Msg6BitWords,msgsent,.false.,'      ') !Unpack to get msgsent
     if(ichk.ne.0) go to 999
     call encode4(message,itone)                 !Encode the information bits
     i1=index(message,'-')
     c=message(i1+1:i1+1)
     if(i1.ge.9 .and. c.ge.'0' .and. c.le.'3') then
        itone=2*itone + (1-npr(2:))             !Inverted '#' sync
     else
        itone=2*itone + npr(2:)                 !Data = MSB, sync = LSB
     endif
  endif

999 return
end subroutine gen4


code:
subroutine gen65(msg0,ichk,msgsent,itone,itype)

! Encodes a JT65 message to yieild itone(1:126)
! Temporarily, does not implement EME shorthands

  use packjt
  character*22 msg0
  character*22 message          !Message to be generated
  character*22 msgsent          !Message as it will be received
  integer itone(126)
  character*3 cok               !'   ' or 'OOO'
  integer dgen(13)
  integer sent(63)
  integer nprc(126)
  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/
  save

  if(msg0(1:1).eq.'@') then
     read(msg0(2:5),*,end=1,err=1) nfreq
     go to 2
1    nfreq=1000
2    itone(1)=nfreq
  else
     message=msg0
     do i=1,22
        if(ichar(message(i:i)).eq.0) then
           message(i:)='                      '
           exit
        endif
     enddo

     do i=1,22                               !Strip leading blanks
        if(message(1:1).ne.' ') exit
        message=message(i+1:)
     enddo

     call chkmsg(message,cok,nspecial,flip)
     ntest=0
     if(flip.lt.0.0) ntest=1
     if(nspecial.eq.0) then
        call packmsg(message,dgen,itype,.false.) !Pack message into 72 bits
        call unpackmsg(dgen,msgsent,.false.,'      ') !Unpack to get message sent
        msgsent(20:22)=cok
        call fmtmsg(msgsent,iz)
        if(ichk.ne.0) go to 999             !Return if checking only

        call rs_encode(dgen,sent)           !Apply Reed-Solomon code
        call interleave63(sent,1)           !Apply interleaving
        call graycode65(sent,63,1)          !Apply Gray code
        nsym=126                            !Symbols per transmission
        k=0
        do j=1,nsym
           if(nprc(j).eq.ntest) then
              k=k+1
              itone(j)=sent(k)+2
           else
              itone(j)=0
           endif
        enddo
     else
        nsym=32
        k=0
        do j=1,nsym
           do n=1,4
              k=k+1
              if(iand(j,1).eq.1) itone(k)=0
              if(iand(j,1).eq.0) itone(k)=10*nspecial
              if(k.eq.126) go to 10
           enddo
        enddo
10      msgsent=message
        itype=7
     endif
  endif

999 return
end subroutine gen65


code:
subroutine gen9(msg0,ichk,msgsent,i4tone,itype)

! Encodes a JT9 message and returns msgsent, the message as it will
! be decoded, and an integer array i4tone(85) of 9-FSK tone values
! in the range 0-8.

  use packjt
  character*22 msg0
  character*22 message                    !Message to be generated
  character*22 msgsent                    !Message as it will be received
  integer*4 i4Msg6BitWords(13)            !72-bit message as 6-bit words
  integer*1 i1Msg8BitBytes(13)            !72 bits and zero tail as 8-bit bytes
  integer*1 i1EncodedBits(207)            !Encoded information-carrying bits
  integer*1 i1ScrambledBits(207)          !Encoded bits after interleaving
  integer*4 i4DataSymbols(69)             !Data symbols (values 0-7)
  integer*4 i4GrayCodedSymbols(69)        !Gray-coded symbols (values 0-7)
  integer*4 i4tone(85)                    !Tone #s, data and sync (values 0-8)
  include 'jt9sync.f90'
  save

  if(msg0(1:1).eq.'@') then
     read(msg0(2:5),*,end=1,err=1) nfreq
     go to 2
1    nfreq=1000
2    i4tone(1)=nfreq
  else
     message=msg0
     do i=1,22
        if(ichar(message(i:i)).eq.0) then
           message(i:)='                      '
           exit
        endif
     enddo

     do i=1,22                               !Strip leading blanks
        if(message(1:1).ne.' ') exit
        message=message(i+1:)
     enddo

     call packmsg(message,i4Msg6BitWords,itype,.false.) !Pack into 12 6-bit bytes
     call unpackmsg(i4Msg6BitWords,msgsent,.false.,'      ') !Unpack to get msgsent
     if(ichk.ne.0) go to 999
     call entail(i4Msg6BitWords,i1Msg8BitBytes)  !Add tail, make 8-bit bytes
     nsym2=206
     call encode232(i1Msg8BitBytes,nsym2,i1EncodedBits)   !Encode K=32, r=1/2
     i1EncodedBits(207)=0
     call interleave9(i1EncodedBits,1,i1ScrambledBits)    !Interleave bits
     i1ScrambledBits(207)=0
     call packbits(i1ScrambledBits,69,3,i4DataSymbols)    !Pk 3-bits into words
     call graycode(i4DataSymbols,69,1,i4GrayCodedSymbols) !Apply Gray code

! Insert sync symbols at ntone=0 and add 1 to the data-tone numbers.
     j=0
     do i=1,85
        if(isync(i).eq.1) then
           i4tone(i)=0
        else
           j=j+1
           i4tone(i)=i4GrayCodedSymbols(j)+1
        endif
     enddo
  endif

999 return
end subroutine gen9


code:
subroutine geniscat(msg,msgsent,itone)

! Generate an ISCAT waveform.

  parameter (NSZ=1291)
  character msg*28,msgsent*28                    !Message to be transmitted
  integer imsg(30)
  integer itone(NSZ)
  real*8 sps
  character c*42
  integer icos(4)                                !Costas array
  data icos/0,1,3,2/
  data nsync/4/,nlen/2/,ndat/18/
  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ /.?@-'/

  sps=256.d0*12000.d0/11025.d0
  nsym=int(30*12000.d0/sps)
  nblk=nsync+nlen+ndat

  do i=22,1,-1
     if(msg(i:i).ne.' ' .and. msg(i:i).ne.char(0)) exit
  enddo
  nmsg=i
  msglen=nmsg+1
  k=0
  kk=1
  imsg(1)=40                                  !Always start with BOM char: '@'
  do i=1,nmsg                                 !Define the tone sequence
     imsg(i+1)=36                             !Illegal char set to blank
     do j=1,42
        if(msg(i:i).eq.c(j:j)) imsg(i+1)=j-1
     enddo
  enddo

  do i=1,nsym                                 !Total symbols in 30 s
     j=mod(i-1,nblk)+1
     if(j.le.nsync) then
        itone(i)=icos(j)                      !Insert 4x4 Costas array
     else if(j.gt.nsync .and. j.le.nsync+nlen) then
        itone(i)=msglen                       !Insert message-length indicator
        if(j.ge.nsync+2) then
           n=msglen + 5*(j-nsync-1)
           if(n.gt.41) n=n-42
           itone(i)=n
        endif
     else
        k=k+1
        kk=mod(k-1,msglen)+1
        itone(i)=imsg(kk)
     endif
  enddo
  msgsent=msg

  return
end subroutine geniscat


code:
program genmet

! Generate metric table for M-FSK modulation.

  character*12 arg
  real*4 r(0:255)
  integer hist(2,-128:128)
  data hist/514*0/,idum/-1/

  lim(x)=min(127,max(-128,nint(scale*x)))

  nargs=iargc()
  if(nargs.ne.5) then
     print*,'Usage: genmet ncoh nadd m0 snr iters'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) ncoh
  call getarg(2,arg)
  read(arg,*) nadd
  call getarg(3,arg)
  read(arg,*) m0
  call getarg(4,arg)
  read(arg,*) snr
  call getarg(5,arg)
  read(arg,*) iters

  ntones=2**m0
  xm0=m0
  scale=5.0
  fac=sqrt(1.0/nadd)
  s=sqrt(10.0**(0.1*snr))
  hist=0
  nerr=0

  call sgran()

  do iter=1,iters
     do i=0,ntones-1
        r(i)=0.
        do n=1,nadd
           x1=0.707*gran()
           y1=0.707*gran()
           if(i.eq.0) x1=x1+s
           if(ncoh.eq.0) r(i)=r(i) + x1*x1 + y1*y1
           if(ncoh.ne.0) r(i)=r(i) + x1
        enddo
        r(i)=fac*r(i)
     enddo
     do m=0,m0-1
        n=2**m
        r1=0.
        r2=0.
        do i=0,ntones-1
           if(iand(i,n).ne.0) r1=max(r1,r(i))
           if(iand(i,n).eq.0) r2=max(r2,r(i))
        enddo
        don=r2-r1
        doff=r1-r2
        if(don.lt.0.0) nerr=nerr+1
        j1=lim(doff)
        hist(1,j1)=hist(1,j1)+1
        j2=lim(don)
        hist(2,j2)=hist(2,j2)+1
     enddo
  enddo

  do i=-128,127
     write(13,1010) i/scale,hist(1,i)/(xm0*iters),hist(2,i)/(xm0*iters)
1010 format(f8.3,2f12.9)
  enddo

  ber=nerr/(xm0*iters)
  write(*,1020) nadd,m0,snr,ber
1020 format('nadd:',i3,'   m0:',i2,'   snr: 'f5.1,'   BER:',f8.3)

  xln2=log(2.0)
  do i=-128,127
     p1=hist(2,i)/(xm0*iters)
     p0=hist(1,i)/(xm0*iters)
     if(p0+p1.eq.0.0 .and. i.lt.0) p0=1.e-6
     if(p0+p1.eq.0.0 .and. i.gt.0) p1=1.e-6
     xlhd0=log(max(0.001,2.0*p0/(p0+p1)))/xln2
     xlhd1=log(max(0.001,2.0*p1/(p0+p1)))/xln2
     write(14,1012) i/scale,xlhd0,xlhd1,p0/(p0+p1),p1/(p0+p1)
1012 format(f7.1,2f8.3,2f9.6)
  enddo

999 end program genmet


code:
subroutine genmsk144(msg0,mygrid,ichk,bcontest,msgsent,i4tone,itype)
! s8 + 48bits + s8 + 80 bits = 144 bits (72ms message duration)
!
! Encode an MSK144 message
! Input:
!   - msg0     requested message to be transmitted
!   - ichk     if ichk=1, return only msgsent
!              if ichk.ge.10000, set imsg=ichk-10000 for short msg
!   - msgsent  message as it will be decoded
!   - i4tone   array of audio tone values, 0 or 1
!   - itype    message type
!                 1 = standard message  "Call_1 Call_2 Grid/Rpt"
!                 2 = type 1 prefix
!                 3 = type 1 suffix
!                 4 = type 2 prefix
!                 5 = type 2 suffix
!                 6 = free text (up to 13 characters)
!                 7 = short message     "<Call_1 Call2> Rpt"

  use iso_c_binding, only: c_loc,c_size_t
  use packjt
  use hashing
  character*22 msg0
  character*22 message                    !Message to be generated
  character*22 msgsent                    !Message as it will be received
  character*6 mygrid
  integer*4 i4Msg6BitWords(13)            !72-bit message as 6-bit words
  integer*4 i4tone(144)                   !
  integer*1, target:: i1Msg8BitBytes(10)  !80 bits represented in 10 bytes
  integer*1 codeword(128)                 !Encoded bits before re-ordering
  integer*1 msgbits(80)                   !72-bit message + 8-bit hash
  integer*1 bitseq(144)                   !Tone #s, data and sync (values 0-1)
  integer*1 i1hash(4)
  integer*1 s8(8)
  logical bcontest
  real*8 pp(12)
  real*8 xi(864),xq(864),pi,twopi
  data s8/0,1,1,1,0,0,1,0/
  equivalence (ihash,i1hash)
  logical first
  data first/.true./
  save

  if(first) then
    first=.false.
    nsym=128
    pi=4.0*atan(1.0)
    twopi=8.*atan(1.0)
    do i=1,12
      pp(i)=sin((i-1)*pi/12)
    enddo
  endif

  if(msg0(1:1).eq.'@') then                    !Generate a fixed tone
     read(msg0(2:5),*,end=1,err=1) nfreq       !at specified frequency
     go to 2
1    nfreq=1000
2    i4tone(1)=nfreq
  else
     message=msg0
     do i=1,22
        if(ichar(message(i:i)).eq.0) then
           message(i:)='                      '
           exit
        endif
     enddo

     do i=1,22                               !Strip leading blanks
        if(message(1:1).ne.' ') exit
        message=message(i+1:)
     enddo

     if(message(1:1).eq.'<') then
        call genmsk40(message,msgsent,ichk,i4tone,itype)
        if(itype.lt.0) go to 999
        i4tone(41)=-40
        go to 999
     endif

     call packmsg(message,i4Msg6BitWords,itype,bcontest) !Pack into 12 6-bit bytes
     call unpackmsg(i4Msg6BitWords,msgsent,bcontest,mygrid) !Unpack to get msgsent

     if(ichk.eq.1) go to 999
     i4=0
     ik=0
     im=0
     do i=1,12
       nn=i4Msg6BitWords(i)
       do j=1, 6
         ik=ik+1
         i4=i4+i4+iand(1,ishft(nn,j-6))
         i4=iand(i4,255)
         if(ik.eq.8) then
           im=im+1
           i1Msg8BitBytes(im)=i4
           ik=0
         endif
       enddo
     enddo

     ihash=nhash(c_loc(i1Msg8BitBytes),int(9,c_size_t),146)
     ihash=2*iand(ihash,32767)                   !Generate the 8-bit hash
     i1Msg8BitBytes(10)=i1hash(1)                !Hash code to byte 10

     mbit=0
     do i=1, 10
       i1=i1Msg8BitBytes(i)
         do ibit=1,8
           mbit=mbit+1
           msgbits(mbit)=iand(1,ishft(i1,ibit-8))
         enddo
     enddo

     call encode_msk144(msgbits,codeword)

!Create 144-bit channel vector:
!8-bit sync word + 48 bits + 8-bit sync word + 80 bits
     bitseq=0
     bitseq(1:8)=s8
     bitseq(9:56)=codeword(1:48)
     bitseq(57:64)=s8
     bitseq(65:144)=codeword(49:128)
     bitseq=2*bitseq-1

     xq(1:6)=bitseq(1)*pp(7:12)   !first bit is mapped to 1st half-symbol on q
     do i=1,71
       is=(i-1)*12+7
       xq(is:is+11)=bitseq(2*i+1)*pp
     enddo
     xq(864-5:864)=bitseq(1)*pp(1:6)   !last half symbol
     do i=1,72
       is=(i-1)*12+1
       xi(is:is+11)=bitseq(2*i)*pp
     enddo
! Map I and Q  to tones.
    i4tone=0
    do i=1,72
      i4tone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2;
      i4tone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,144)+1)-1)/2;
    enddo
  endif

! Flip polarity
  i4tone=-i4tone+1

999 return
end subroutine genmsk144


code:
subroutine genmsk40(msg,msgsent,ichk,itone,itype)

  use hashing
  character*22 msg,msgsent,hashmsg
  character*4 crpt,rpt(0:15)
  logical first
  integer*4 itone(144)
  integer*1 message(16),codeword(32),bitseq(40)
  integer*1 s8r(8)
  data s8r/1,0,1,1,0,0,0,1/   ! Sync word is reversed wrt msk144 sync word.
  data first/.true./
  data rpt/"-03 ","+00 ","+03 ","+06 ","+10 ","+13 ","+16 ", &
           "R-03","R+00","R+03","R+06","R+10","R+13","R+16", &
           "RRR ","73  "/
  save first,rpt

  itype=-1
  msgsent='*** bad message ***'
  itone=0
  i1=index(msg,'>')
  if(i1.lt.9) go to 900
  call fmtmsg(msg,iz)
  crpt=msg(i1+2:i1+5)
  do i=0,15
     if(crpt.eq.rpt(i)) go to 10
  enddo
  go to 900

10 irpt=i                                   !Report index, 0-15
  if(ichk.lt.10000) then
     hashmsg=msg(2:i1-1)
     call hash(hashmsg,22,ihash)
     ihash=iand(ihash,4095)                 !12-bit hash
     ig=16*ihash + irpt                     !4-bit report
  else
     ig=ichk-10000
  endif

  do i=1,16
    message(i)=iand(1,ishft(ig,1-i))
  enddo

  call encode_msk40(message,codeword)
!  write(*,'(a6,i6,2x,a6,i6,2x,a6,i6)') ' msg: ',ig,'rprt: ',irpt,'hash: ',ihash
!  write(*,'(a6,32i1)') '  cw: ',codeword

  bitseq(1:8)=s8r
  bitseq(9:40)=codeword
  bitseq=2*bitseq-1

! Map I and Q  to tones.
  itone=0
  do i=1, 20
    itone(2*i-1)=(bitseq(2*i)*bitseq(2*i-1)+1)/2;
    itone(2*i)=-(bitseq(2*i)*bitseq(mod(2*i,40)+1)-1)/2;
  enddo

! Flip polarity
  itone=-itone+1

  msgsent=msg
  itype=7

900 return
end subroutine genmsk40


code:
subroutine genmsk_short(msg,msgsent,ichk,itone,itype)

  use hashing
  character*22 msg,msgsent
  character*3 crpt,rpt(0:7)
  logical first
  integer itone(35)
  integer ig24(0:4096-1)                  !Codewords for Golay (24,12) code
  integer b11(11)
  data b11/1,1,1,0,0,0,1,0,0,1,0/         !Barker 11 code
  data rpt /'26 ','27 ','28 ','R26','R27','R28','RRR','73 '/
  data first/.true./
  save first,ig24

  if(first) then
     call golay24_table(ig24)             !Define the Golay(24,12) codewords
     first=.false.
  endif

  itype=-1
  msgsent='*** bad message ***'
  itone=0
  i1=index(msg,'>')
  if(i1.lt.9) go to 900
  call fmtmsg(msg,iz)
  crpt=msg(i1+2:i1+5)
  do i=0,7
     if(crpt.eq.rpt(i)) go to 10
  enddo
  go to 900

10 irpt=i                               !Report index, 0-7
  if(ichk.lt.10000) then
     call hash(msg(2:i1-1),i1-2,ihash)
     ihash=iand(ihash,511)                 !9-bit hash for the two callsigns
     ig=8*ihash + irpt                     !12-bit message information
  else
     ig=ichk-10000
  endif
  ncodeword=ig24(ig)
  itone(1:11)=b11                       !Insert the Barker-11 code
  n=2**24
  do i=12,35                            !Insert codeword into itone array
     n=n/2
     itone(i)=0
     if(iand(ncodeword,n).ne.0) itone(i)=1
  enddo
  msgsent=msg
  itype=7

  n=count(itone(1:35).eq.0)
  if(mod(n,2).ne.0) stop 'Parity error in genmsk_short.'

900 return
end subroutine genmsk_short

subroutine hash_calls(calls,ih9)

  use hashing
  character*(*) calls
  i1=index(calls,'>')
  call hash(calls(2:i1-1),i1-2,ih9)
  ih9=iand(ih9,511)                      !9-bit hash for the two callsigns

  return
end subroutine hash_calls


code:
subroutine genqra64(msg0,ichk,msgsent,itone,itype)

! Encodes a QRA64 message to yield itone(1:84)

  use packjt
  character*22 msg0
  character*22 message          !Message to be generated
  character*22 msgsent          !Message as it will be received
  integer itone(84)
  character*3 cok               !'   ' or 'OOO'
  logical old_qra_sync
  integer dgen(13)
  integer sent(63)
  integer icos7(0:6)
  data icos7/4,2,5,6,1,3,0/                   !Costas 7x7 tone pattern
  save

  if(msg0(1:1).eq.'@') then
     read(msg0(2:5),*,end=1,err=1) nfreq
     go to 2
1    nfreq=1000
2    itone(1)=nfreq
     write(msgsent,1000) nfreq
1000 format(i5,' Hz')
  else
     message=msg0
     do i=1,22
        if(ichar(message(i:i)).eq.0) then
           message(i:)='                      '
           exit
        endif
     enddo

     do i=1,22                               !Strip leading blanks
        if(message(1:1).ne.' ') exit
        message=message(i+1:)
     enddo

     call chkmsg(message,cok,nspecial,flip)
     call packmsg(message,dgen,itype,.false.) !Pack message into 72 bits
     call unpackmsg(dgen,msgsent,.false.,'      ') !Unpack to get message sent
     if(ichk.ne.0) go to 999             !Return if checking only
     call qra64_enc(dgen,sent)           !Encode using QRA64

     nsync=10
     inquire(file='old_qra_sync',exist=old_qra_sync)
     if(old_qra_sync) nsync=1

     itone(1:7)=nsync*icos7              !Insert 7x7 Costas array in 3 places
     itone(8:39)=sent(1:32)
     itone(40:46)=nsync*icos7
     itone(47:77)=sent(33:63)
     itone(78:84)=nsync*icos7
  endif

999 return
end subroutine genqra64


code:
subroutine geodist(Eplat,Eplon,Stlat,Stlon,Az,Baz,Dist)
  implicit none
  real eplat, eplon, stlat, stlon, az, baz, dist

! JHT: In actual fact, I use the first two arguments for "My Location",
!     the second two for "His location"; West longitude is positive.

!      Taken directly from:
!      Thomas, P.D., 1970, Spheroidal geodesics, reference systems,
!      & local geometry, U.S. Naval Oceanographi!Office SP-138,
!      165 pp.
!      assumes North Latitude and East Longitude are positive

!      EpLat, EpLon = End point Lat/Long
!      Stlat, Stlon = Start point lat/long
!      Az, BAz = direct & reverse azimuith
!      Dist = Dist (km); Deg = central angle, discarded

  real BOA, F, P1R, P2R, L1R, L2R, DLR, T1R, T2R, TM,          &
       DTM, STM, CTM, SDTM,CDTM, KL, KK, SDLMR, L,                &
       CD, DL, SD, T, U, V, D, X, E, Y, A, FF64, TDLPM,           &
       HAPBR, HAMBR, A1M2, A2M1

  real AL,BL,D2R,Pi2

  data AL/6378206.4/              ! Clarke 1866 ellipsoid
  data BL/6356583.8/
!      real pi /3.14159265359/
  data D2R/0.01745329251994/      ! degrees to radians conversion factor
  data Pi2/6.28318530718/

  if(abs(Eplat-Stlat).lt.0.02 .and. abs(Eplon-Stlon).lt.0.02) then
     Az=0.
     Baz=180.0
     Dist=0
     go to 999
  endif

  BOA = BL/AL
  F = 1.0 - BOA
! Convert st/end pts to radians
  P1R = Eplat * D2R
  P2R = Stlat * D2R
  L1R = Eplon * D2R
  L2R = StLon * D2R
  DLR = L2R - L1R                 ! DLR = Delta Long in Rads
  T1R = ATan(BOA * Tan(P1R))
  T2R = ATan(BOA * Tan(P2R))
  TM = (T1R + T2R) / 2.0
  DTM = (T2R - T1R) / 2.0
  STM = Sin(TM)
  CTM = Cos(TM)
  SDTM = Sin(DTM)
  CDTM = Cos(DTM)
  KL = STM * CDTM
  KK = SDTM * CTM
  SDLMR = Sin(DLR/2.0)
  L = SDTM * SDTM + SDLMR * SDLMR * (CDTM * CDTM - STM * STM)
  CD = 1.0 - 2.0 * L
  DL = ACos(CD)
  SD = Sin(DL)
  T = DL/SD
  U = 2.0 * KL * KL / (1.0 - L)
  V = 2.0 * KK * KK / L
  D = 4.0 * T * T
  X = U + V
  E = -2.0 * CD
  Y = U - V
  A = -D * E
  FF64 = F * F / 64.0
  Dist = AL*SD*(T -(F/4.0)*(T*X-Y)+FF64*(X*(A+(T-(A+E)                 &
       /2.0)*X)+Y*(-2.0*D+E*Y)+D*X*Y))/1000.0
  TDLPM = Tan((DLR+(-((E*(4.0-X)+2.0*Y)*((F/2.0)*T+FF64*               &
       (32.0*T+(A-20.0*T)*X-2.0*(D+2.0)*Y))/4.0)*Tan(DLR)))/2.0)
  HAPBR = ATan2(SDTM,(CTM*TDLPM))
  HAMBR = Atan2(CDTM,(STM*TDLPM))
  A1M2 = Pi2 + HAMBR - HAPBR
  A2M1 = Pi2 - HAMBR - HAPBR

1 If ((A1M2 .ge. 0.0) .AND. (A1M2 .lt. Pi2)) GOTO 5
  If (A1M2 .lt. Pi2) GOTO 4
  A1M2 = A1M2 - Pi2
  GOTO 1
4 A1M2 = A1M2 + Pi2
  GOTO 1

! All of this gens the proper az, baz (forward and back azimuth)

5 If ((A2M1 .ge. 0.0) .AND. (A2M1 .lt. Pi2)) GOTO 9
  If (A2M1 .lt. Pi2) GOTO 8
  A2M1 = A2M1 - Pi2
  GOTO 5
8 A2M1 = A2M1 + Pi2
  GOTO 5

9 Az = A1M2 / D2R
  BAZ = A2M1 / D2R

!Fix the mirrored coords here.

  az = 360.0 - az
  baz = 360.0 - baz

999 return
end subroutine geodist


code:
subroutine getlags(nsps8,lag0,lag1,lag2)
  if(nsps8.eq.864) then
     lag1=39
     lag2=291
     lag0=123
  else if(nsps8.eq.1920) then
     lag1=70
     lag2=184
     lag0=108
  else if(nsps8.eq.5120) then
     lag1=84
     lag2=129
     lag0=99
  else if(nsps8.eq.10368) then
     lag1=91
     lag2=112
     lag0=98
  else if(nsps8.eq.31500) then
     lag1=93
     lag2=102
     lag0=96
  else
     stop 'Error in getlags'
  endif

  return
end subroutine getlags


code: gfortran -o chkfft chkfft3.f90 four2a.f90 gran.c /jtsdk/fftw3f/libfftw3f-3.dll


code:
void golay24_table_(int cw[])
{
  // Compute arnd return the table of 4096 codewords for the Golay (24,12) code.

  // Array y contains the twelve rows (columns) of the parity-check matrix
  int y[12] = { 0x7ff, 0xee2, 0xdc5, 0xb8b, 0xf16, 0xe2d,
		0xc5b, 0x8b7, 0x96e, 0xadc, 0xdb8, 0xb71 };

  int c[2];          /* Codeword composed of 12-bit info and 12-bit parity */

  int i,j,k;
  int aux;
  int weight(int vector);

  for(k=0; k<4096; k++) {
    c[0] = k;
    c[1] = 0;
    for (i=0; i<12; i++) {
      aux = 0;
      for (j=0; j<12; j++) {
	aux = aux ^ ((c[0] & y[i]) >> j & 1);
      }
      c[1] = (c[1] << 1) ^ aux;
    }
    cw[k]=4096*c[0] + c[1];
  }
}


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

/* Generate gaussian random float with mean=0 and std_dev=1 */
float gran_()
{
  float fac,rsq,v1,v2;
  static float gset;
  static int iset;

  if(iset){
    /* Already got one */
    iset = 0;
    return gset;
  }
  /* Generate two evenly distributed numbers between -1 and +1
   * that are inside the unit circle
   */
  do {
    v1 = 2.0 * (float)rand() / RAND_MAX - 1;
    v2 = 2.0 * (float)rand() / RAND_MAX - 1;
    rsq = v1*v1 + v2*v2;
  } while(rsq >= 1.0 || rsq == 0.0);
  fac = sqrt(-2.0*log(rsq)/rsq);
  gset = v1*fac;
  iset++;
  return v2*fac;
}


code:
subroutine graycode(ia,n,idir,ib)

  integer ia(n),ib(n)
  do i=1,n
     ib(i)=igray(ia(i),idir)
  enddo

  return
end subroutine graycode


code:
subroutine graycode65(dat,n,idir)

  integer dat(n)
  do i=1,n
     dat(i)=igray(dat(i),idir)
  enddo

  return
end subroutine graycode65


code:
subroutine grayline(nyear,month,nday,uth,mygrid,nduration,isun)

  character*6 mygrid
  real LST
  real lat,lon

  call grid2deg(MyGrid,elon,lat)
  lon=-elon

  uth0=uth-0.5*nduration/60.0
  uth1=uth+0.5*nduration/60.0

  call sun(nyear,month,nday,uth0,lon,lat,RASun,DecSun,LST,        &
       AzSun,ElSun0,mjd,day)
  call sun(nyear,month,nday,uth1,lon,lat,RASun,DecSun,LST,        &
       AzSun,ElSun1,mjd,day)

  elchk=-0.8333
  isun=-1
  if(elsun0.lt.elchk .and. elsun1.ge.elchk) then
     isun=0
  else if(elsun0.gt.elchk .and. elsun1.le.elchk) then
     isun=2
  else if(elsun1.gt.elchk) then
     isun=1
  else
     isun=3
  endif

  return
end subroutine grayline


code:
subroutine grid2deg(grid0,dlong,dlat)

! Converts Maidenhead grid locator to degrees of West longitude
! and North latitude.

  character*6 grid0,grid
  character*1 g1,g2,g3,g4,g5,g6

  grid=grid0
  i=ichar(grid(5:5))
  if(grid(5:5).eq.' ' .or. i.le.64 .or. i.ge.128) grid(5:6)='mm'

  if(grid(1:1).ge.'a' .and. grid(1:1).le.'z') grid(1:1)=          &
        char(ichar(grid(1:1))+ichar('A')-ichar('a'))
  if(grid(2:2).ge.'a' .and. grid(2:2).le.'z') grid(2:2)=          &
        char(ichar(grid(2:2))+ichar('A')-ichar('a'))
  if(grid(5:5).ge.'A' .and. grid(5:5).le.'Z') grid(5:5)=          &
        char(ichar(grid(5:5))-ichar('A')+ichar('a'))
  if(grid(6:6).ge.'A' .and. grid(6:6).le.'Z') grid(6:6)=          &
        char(ichar(grid(6:6))-ichar('A')+ichar('a'))

  g1=grid(1:1)
  g2=grid(2:2)
  g3=grid(3:3)
  g4=grid(4:4)
  g5=grid(5:5)
  g6=grid(6:6)

  nlong = 180 - 20*(ichar(g1)-ichar('A'))
  n20d = 2*(ichar(g3)-ichar('0'))
  xminlong = 5*(ichar(g5)-ichar('a')+0.5)
  dlong = nlong - n20d - xminlong/60.0
  nlat = -90+10*(ichar(g2)-ichar('A')) + ichar(g4)-ichar('0')
  xminlat = 2.5*(ichar(g6)-ichar('a')+0.5)
  dlat = nlat + xminlat/60.0

  return
end subroutine grid2deg


code:
subroutine hash(string,len,ihash)
  use iso_c_binding, only: c_loc,c_size_t
  use hashing
  parameter (MASK15=32767)
!  character*(*), target :: string
  character*1, target :: string
     i=nhash(c_loc(string),int(len,c_size_t),146)
     ihash=iand(i,MASK15)
  return
end subroutine hash


code:
module hashing
  interface
     integer(c_int32_t) function nhash (key, length, initval) bind(C, name="nhash")
       use iso_c_binding, only: c_ptr, c_size_t, c_int32_t
       type(c_ptr), intent(in), value :: key
       integer(c_size_t), intent(in), value :: length
       integer(c_int32_t), intent(in), value :: initval
     end function nhash
  end interface
end module hashing


code:
subroutine hint65(s3,mrs,mrs2,nadd,nflip,mycall,hiscall,hisgrid,qual,decoded)

  use packjt
  use prog_args
  parameter (MAXCALLS=10000,MAXRPT=63)
  parameter (MAXMSG=2*MAXCALLS + 2 + MAXRPT)
  real s3(64,63)
  integer*1 sym1(0:62,MAXMSG)
  integer*1 sym2(0:62,MAXMSG)
  integer mrs(63),mrs2(63)
  integer dgen(12),sym(0:62),sym_rev(0:62)
  character*6 mycall,hiscall,hisgrid,call2(MAXCALLS)
  character*4 grid2(MAXCALLS),rpt(MAXRPT)
  character callsign*12,grid*4
  character*180 line
  character ceme*3,msg*22,msg00*22
  character*22 msg0(MAXMSG),decoded
  logical*1 eme(MAXCALLS)
  logical first
  data first/.true./
  data rpt/'-01','-02','-03','-04','-05',          &
           '-06','-07','-08','-09','-10',          &
           '-11','-12','-13','-14','-15',          &
           '-16','-17','-18','-19','-20',          &
           '-21','-22','-23','-24','-25',          &
           '-26','-27','-28','-29','-30',          &
           'R-01','R-02','R-03','R-04','R-05',     &
           'R-06','R-07','R-08','R-09','R-10',     &
           'R-11','R-12','R-13','R-14','R-15',     &
           'R-16','R-17','R-18','R-19','R-20',     &
           'R-21','R-22','R-23','R-24','R-25',     &
           'R-26','R-27','R-28','R-29','R-30',     &
           'RO','RRR','73'/
  save first,sym1,nused,msg0,sym2

  first=.true.   !### For now, at least: always recompute hypothetical messages
  if(first) then
     neme=0
     open(23,file=trim(data_dir)//'/CALL3.TXT',status='unknown')
     icall=0
     j=0
     do i=1,MAXCALLS
        read(23,1002,end=10) line
1002    format(a80)
        if(line(1:4).eq.'ZZZZ') cycle
        if(line(1:2).eq.'//') cycle
        i1=index(line,',')
        if(i1.lt.4) cycle
        i2=index(line(i1+1:),',')
        if(i2.lt.5) cycle
        i2=i2+i1
        i3=index(line(i2+1:),',')
        if(i3.lt.1) i3=index(line(i2+1:),' ')
        i3=i2+i3
        callsign=line(1:i1-1)
        grid=line(i1+1:i1+4)
        ceme=line(i2+1:i3-1)
        eme(i)=ceme.eq.'EME'
        if(neme.eq.1 .and. (.not.eme(i))) cycle
        j=j+1
        call2(j)=callsign(1:6)               !### Fix for compound callsigns!
        grid2(j)=grid
     enddo
10   ncalls=j
     if(ncalls.lt.10) then
        write(*,1010) ncalls
1010    format('CALL3.TXT very short (N =',i2,') or missing?')
     endif
     close(23)

! NB: generation of test messages is not yet complete!
     j=0
     do i=-1,ncalls
        if(i.eq.0 .and. hiscall.eq.'      ' .and. hisgrid(1:4).eq.'    ') cycle
        mz=2
        if(i.eq.-1) mz=1
        if(i.eq.0) mz=65
        do m=1,mz
           j=j+1
           if(i.eq.-1) then
              msg='0123456789ABC'
           else if(i.eq.0) then
              if(m.eq.1) msg=mycall//' '//hiscall//' '//hisgrid(1:4)
              if(m.eq.2) msg='CQ '//hiscall//' '//hisgrid(1:4)
              if(m.ge.3) msg=mycall//' '//hiscall//' '//rpt(m-2)
           else
              if(m.eq.1)  msg=mycall//' '//call2(i)//' '//grid2(i)
              if(m.eq.2)  msg='CQ '//call2(i)//' '//grid2(i)
           endif
           call fmtmsg(msg,iz)
           call packmsg(msg,dgen,itype,.false.) !Pack message into 72 bits
           call rs_encode(dgen,sym_rev)            !RS encode
           sym(0:62)=sym_rev(62:0:-1)
           sym1(0:62,j)=sym

           call interleave63(sym_rev,1)            !Interleave channel symbols
           call graycode(sym_rev,63,1,sym_rev)     !Apply Gray code
           sym2(0:62,j)=sym_rev(0:62)
           msg0(j)=msg
        enddo
     enddo
     nused=j
     first=.false.
  endif

  ref0=0.
  do j=1,63
     ref0=ref0 + s3(mrs(j)+1,j)
  enddo

  u1=0.
  u1=-99.0
  u2=u1

! Find u1 and u2 (best and second-best) codeword from a list, using
! a bank of matched filters on the symbol spectra s3(i,j).
  ipk=1
  ipk2=0
  msg00='                      '
  do k=1,nused
     if(k.ge.2 .and. k.le.64 .and. nflip.lt.0) cycle
! Test all messages if nflip=+1; skip the CQ messages if nflip=-1.
     if(nflip.gt.0 .or. msg0(k)(1:3).ne.'CQ ') then
        psum=0.
        ref=ref0
        do j=1,63
           i=sym2(j-1,k)+1
           psum=psum + s3(i,j)
           if(i.eq.mrs(j)+1) ref=ref - s3(i,j) + s3(mrs2(j)+1,j)
        enddo
        p=psum/ref

        if(p.gt.u1) then
           if(msg0(k).ne.msg00) then
              ipk2=ipk
              u2=u1
           endif
           u1=p
           ipk=k
           msg00=msg0(k)
        endif
        if(msg0(k).ne.msg00 .and. p.gt.u2) then
           u2=p
           ipk2=k
        endif
     endif
  enddo

!### Just in case ???
!  rewind 77
!  write(77,*) u1,u2,ipk,ipk2
!  call flush(77)
!###

  decoded='                      '
  bias=max(1.12*u2,0.35)
  if(nadd.ge.4) bias=max(1.08*u2,0.45)
  if(nadd.ge.8) bias=max(1.04*u2,0.60)
  qual=100.0*(u1-bias)
!  write(*,3301) u1,u2,u1/u2,bias,qual,nadd,ipk,ipk2
!3301 format(5f6.2,i3,2i6)
  qmin=1.0
  if(qual.ge.qmin) decoded=msg0(ipk)

  return
end subroutine hint65


code:
subroutine hspec(id2,k,nutc0,ntrpdepth,nrxfreq,ntol,bmsk144,bcontest,  &
     btrain,pcoeffs,ingain,mycall,hiscall,bshmsg,bswl,datadir,green,s, &
     jh,pxmax,dbNoGain,line1,mygrid)

! Input:
!  k         pointer to the most recent new data
!  nutc0     UTC for display of decode(s)
!  ntrpdepth TR period and 1000*ndepth
!  nrxfreq   Rx audio center frequency
!  ntol      Decoding range is +/- ntol
!  bmsk144   Boolean, true if in MSK144 mode
!  btrain    Boolean, turns on training in MSK144 mode
!  ingain    Relative gain for spectra

! Output:
!  green()   power
!  s()       spectrum for horizontal spectrogram
!  jh        index of most recent data in green(), s()

  parameter (JZ=703)
  character*80 line1
  character*512 datadir
  character*12 mycall,hiscall
  character*6 mygrid
  integer*2 id2(0:120*12000-1)
  logical*1 bmsk144,bcontest,bshmsg,btrain,bswl
  real green(0:JZ-1)
  real s(0:63,0:JZ-1)
  real x(512)
  real*8 pcoeffs(5)
  complex cx(0:256)
  data rms/999.0/,k0/99999999/
  equivalence (x,cx)
  save ja,rms0

  ndepth=ntrpdepth/1000
  ntrperiod=ntrpdepth - 1000*ndepth
  gain=10.0**(0.1*ingain)
  nfft=512
  nstep=nfft
  nblks=7
  if(ntrperiod.lt.30) then
     nstep=256
     nblks=14
  endif

  if(k.gt.30*12000) go to 900
  if(k.lt.nfft) then
     jh=0
     go to 900                                 !Wait for enough samples to start
  endif

  if(k.lt.k0) then                             !Start a new data block
     ja=-nstep
     jh=-1
     rms0=0.0
  endif

  pxmax = 0;
  do iblk=1,nblks
     if(jh.lt.JZ-1) jh=jh+1
     ja=ja+nstep
     jb=ja+nfft-1
     x=id2(ja:jb)
     sq=dot_product(x,x)
     xmax = maxval(x);
     xmin = abs(minval(x));
     if (xmin > xmax) xmax = xmin;
     if (xmax.gt.0.0) pxmax=20.0*log10(xmax);
     rms=sqrt(gain*sq/nfft)
     rms2=sqrt(sq/nfft);
     green(jh)=0.
     if(rms.gt.0.0) then
        green(jh)=20.0*log10(rms)
        dbNoGain=20.0*log10(rms2);
     endif
     call four2a(x,nfft,1,-1,0)                   !Real-to-complex FFT
     df=12000.0/nfft
     fac=(1.0/nfft)**2
     do i=1,64
        j=2*i
        sx=real(cx(j))**2 + aimag(cx(j))**2 + real(cx(j-1))**2 +        &
             aimag(cx(j-1))**2
        s(i-1,jh)=fac*gain*sx
     enddo
     if(ja+2*nfft.gt.k) exit
  enddo
  k0=k

  if(bmsk144) then
     if(k.ge.7168) then
        tsec=(k-7168)/12000.0
        k0=k-7168
        tt1=sum(float(abs(id2(k0:k0+3583))))
        k0=k-3584
        tt2=sum(float(abs(id2(k0:k0+3583))))
        if(tt1.ne.0.0 .and. tt2.ne.0) then
           call mskrtd(id2(k-7168+1:k),nutc0,tsec,ntol,nrxfreq,ndepth,   &
                mycall,mygrid,hiscall,bshmsg,bcontest,btrain,pcoeffs,bswl,&
                datadir,line1)
        endif
     endif
  endif

900 return
end subroutine hspec


code:
int igray_(int *n0, int *idir)
{
  int n;
  unsigned long sh;
  unsigned long nn;
  n=*n0;

  if(*idir>0) return (n ^ (n >> 1));

  sh = 1;
  nn = (n >> sh);
  while (nn > 0) {
    n ^= nn;
    sh <<= 1;
    nn = (n >> sh);
  }
  return (n);
}


code:
subroutine imopen(plotfile)
  character*(*) plotfile
  common/imcom/ lu,npage

  lu=80
  open(lu,file=plotfile,status='unknown')
  write(lu,1000)
1000 format('%!PS-Adobe-2.0'/                                            &
          '/rightshow { dup stringwidth pop neg 0 rmoveto show } def'/   &
          '/centershow { dup stringwidth pop neg 2 div ',                &
          '0 rmoveto show } def'/                                        &
          '/lt { lineto } def'/'%%Page: 1 1')
  npage=1

  return
end subroutine imopen

subroutine impalette(palette)
  character*(*) palette
  integer r(0:8),g(0:8),b(0:8)
  integer rr,gg,bb
  common/imcom/ lu,npage
  common/imcom2/rr(0:255),gg(0:255),bb(0:255)

  if(palette.eq.'afmhot') then
     do i=0,255
        j=255-i
        rr(i)=min(255,2*j)
        gg(i)=max(0,min(255,2*j-128))
        bb(i)=max(0,min(255,2*j-256))
     enddo
  else if(palette.eq.'hot') then
     do i=0,255
        j=255-i
        rr(i)=min(255,3*j)
        gg(i)=max(0,min(255,3*j-256))
        bb(i)=max(0,min(255,3*j-512))
     enddo
  else
     open(11,file="Palettes/"//palette,status="old")
     do j=0,8
        read(11,*) r(j),g(j),b(j)
     enddo
     close(11)
     do i=0,255
        j0=i/32
        j1=j0+1
        k=i-32*j0
        rr(i)=r(j0) + int((k*(r(j1)-r(j0)))/31 + 0.5)
        gg(i)=g(j0) + int((k*(g(j1)-g(j0)))/31 + 0.5)
        bb(i)=b(j0) + int((k*(b(j1)-b(j0)))/31 + 0.5)
     enddo

  endif

  return
end subroutine impalette

subroutine imclose
  common/imcom/ lu,npage
  write(lu,1000)
1000 format('showpage'/'%%Trailer')
  close(lu)
  return
end subroutine imclose

subroutine imnewpage
  common/imcom/ lu,npage
  npage=npage+1
  write(lu,1000) npage,npage
1000 format('showpage'/'%%Page:',2i4)
  return
end subroutine imnewpage

subroutine imxline(x,y,dx)
! Draw a line from (x,y) to (x+dx,y)  integer r,g,b
  common/imcom/ lu,npage
  write(lu,1000) 72.0*x,72.0*y,72.0*dx
1000 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto stroke')
  return
end subroutine imxline

subroutine imyline(x,y,dy)
! Draw a line from (x,y) to (x,y+dy)
  common/imcom/ lu,npage
  write(lu,1000) 72.0*x,72.0*y,72.0*dy
1000 format('newpath',2f7.1,' moveto 0',f7.1,' rlineto stroke')
  return
end subroutine imyline

subroutine imwidth(width)
  common/imcom/ lu,npage
  write(lu,1000) width
1000 format(f7.1,' setlinewidth')
  return
end subroutine imwidth

subroutine imfont(fontname,npoints)
  character*(*) fontname
  common/imcom/ lu,npage
  write(lu,1000) fontname,npoints
1000 format('/',a,' findfont',i4,' scalefont setfont')
  return
end subroutine imfont

subroutine imstring(string,x,y,just,ndeg)
  character*(*) string
  common/imcom/ lu,npage
  write(lu,1000) 72.0*x,72.0*y,ndeg,string
1000 format(2f7.1,' moveto',i4,' rotate'/'(',a,')')
  if(just.eq.1) write(lu,*) 'rightshow'
  if(just.eq.2) write(lu,*) 'centershow'
  if(just.eq.3) write(lu,*) 'show'
  write(lu,1010) -ndeg
1010 format(i4,' rotate'/)
  return
end subroutine imstring

subroutine imr4mat(z,IP,JP,imax,jmax,zz1,zz2,x,y,dx,dy,nbox)
  real z(IP,JP)
  integer idat(2048)
  common/imcom/ lu,npage

  z1=zz1
  z2=zz2
  if(z1.eq.0.0 .and. z2.eq.0.0) then
     z1=z(1,1)
     z2=z1
     do i=1,imax
        do j=1,jmax
           z1=min(z(i,j),z1)
           z2=max(z(i,j),z2)
        enddo
     enddo
  endif
  scale=255.99/(z2-z1)

  write(lu,1002) 72.0*x,72.0*y,72.0*dx,72.0*dy
1002 format(2f7.1,' translate',2f7.1,' scale')
  write(lu,*) imax,jmax,8,' [',imax,0,0,jmax,0,0,']'
  write(lu,*) '{<'

  do j=1,jmax
     do i=1,imax
        idat(i)=scale*(z(i,j)-z1)
        idat(i)=max(idat(i),0)
        idat(i)=min(idat(i),255)
        idat(i)=255-idat(i)
     enddo
     write(lu,1004) (idat(i),i=1,imax)
1004 format(30z2.2)
  enddo
  write(lu,*) '>} image'
  write(lu,1006) 1.0/(72.0*dx),1.0/(72.0*dy),-72.0*x,-72.0*y
1006 format(2f9.6,' scale',2f7.1,' translate')

  if(nbox.ne.0) then
     write(lu,1010) 72.0*x,72.0*y,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0',             &
          f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
  endif

  return
end subroutine imr4mat

subroutine imr4mat_color(z,IP,JP,imax,jmax,zz1,zz2,x,y,dx,dy,nbox)
  real z(IP,JP)
  integer idat(2048,3)
  integer rr,gg,bb
  common/imcom/ lu,npage
  common/imcom2/rr(0:255),gg(0:255),bb(0:255)

  z1=zz1
  z2=zz2
  if(z1.eq.0.0 .and. z2.eq.0.0) then
     z1=z(1,1)
     z2=z1
     do i=1,imax
        do j=1,jmax
           z1=min(z(i,j),z1)
           z2=max(z(i,j),z2)
        enddo
     enddo
  endif
  scale=255.99/(z2-z1)

  write(lu,1002) 72.0*x,72.0*y,72.0*dx,72.0*dy
1002 format(2f7.1,' translate',2f7.1,' scale')
  write(lu,1003) imax,jmax,8,imax,0,0,jmax,0,0
1003 format(3i5,' [',6i4,']')
  write(lu,1004) imax
1004 format('{currentfile 3',i4,' mul string readhexstring pop} bind'/   &
          'false 3 colorimage')

  do j=1,jmax
     do i=1,imax
        n=scale*(z(i,j)-z1)
        n=max(n,0)
        n=min(n,255)
        idat(i,1)=rr(n)
        idat(i,2)=gg(n)
        idat(i,3)=bb(n)
     enddo
     write(lu,1005) (idat(i,1),idat(i,2),idat(i,3),i=1,imax)
1005 format(30z2.2)
  enddo

  write(lu,1006) 1.0/(72.0*dx),1.0/(72.0*dy),-72.0*x,-72.0*y
1006 format(2f9.6,' scale',2f7.1,' translate')

  if(nbox.ne.0) then
     write(lu,1010) 72.0*x,72.0*y,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0',             &
          f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
  endif

  return
end subroutine imr4mat_color

subroutine imr4pro(p,imax,yy1,yy2,x,y,dx,dy,nbox)
  real p(imax)
  common/imcom/ lu,npage

  y1=yy1
  y2=yy2
  if(y1.eq.0.0 .and. y2.eq.0.0) then
     y1=p(1)
     y2=y1
     do i=1,imax
        y1=min(p(i),y1)
        y2=max(p(i),y2)
     enddo
  endif

  xscale=72.0*dx/imax
  xoff=72.0*x
  yscale=72.0*dy
  if(y1.ne.y2) yscale=yscale/(y2-y1)
  yoff=72.0*y

  write(lu,*) '1.416 setmiterlimit'
  write(lu,1002) xoff+0.5*xscale,yoff+yscale*(p(1)-y1)
1002 format('newpath',2f7.1,' moveto')

  do i=2,imax
     write(lu,1004) xoff+(i-0.5)*xscale,yoff+yscale*(p(i)-y1)
1004 format(2f6.1,' lt')
  enddo
  write(lu,*) 'stroke'

  if(nbox.ne.0) then
     write(lu,1010) xoff,yoff,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0',               &
          f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
  endif

  return
end subroutine imr4pro

subroutine imline(x1,y1,x2,y2)
  common/imcom/ lu,npage
  write(lu,1000) 72*x1,72*y1,72*x2,72*y2
1000 format('newpath',2f7.1,' moveto',2f7.1,' lineto stroke')
  return
end subroutine imline

subroutine imcircle(x,y,radius,shade)
  common/imcom/ lu,npage
  write(lu,1000) shade
1000 format(f7.1,' setgray')
  write(lu,1002) 72*x,72*y,72*radius
1002 format('newpath',3f7.1,' 0 360 arc fill')
  write(lu,1000) 0.0
  write(lu,1004) 72*x,72*y,72*radius
1004 format('newpath',3f7.1,' 0 360 arc stroke')
  return
end subroutine imcircle

subroutine imtriangle(x,y,rr,shade)
  common/imcom/ lu,npage
  write(lu,1000) shade
1000 format(f7.1,' setgray')
  write(lu,1002) 72*x,72*(y+rr)
1002 format('newpath',2f7.1,' moveto ')
  write(lu,1004) 72*(x-rr),72*(y-rr)
1004 format(2f7.1,' lineto ')
  write(lu,1004) 72*(x+rr),72*(y-rr)
  write(lu,*) 'closepath fill 0 setgray'
  write(lu,1002) 72*x,72*(y+rr)
  write(lu,1004) 72*(x-rr),72*(y-rr)
  write(lu,1004) 72*(x+rr),72*(y-rr)
  write(lu,*) 'closepath stroke'

  return
end subroutine imtriangle

subroutine imr4prov(p,jmax,xx1,xx2,x,y,dx,dy,nbox)
  real p(jmax)
  common/imcom/ lu,npage

  x1=xx1
  x2=xx2
  if(x1.eq.0.0 .and. x2.eq.0.0) then
     x1=p(1)
     x2=x1
     do j=1,jmax
        x1=min(p(j),x1)
        x2=max(p(j),x2)
     enddo
  endif

  xscale=72.0*dx
  xoff=72.0*x
  if(x1.ne.x2) xscale=xscale/(x2-x1)

  yscale=72.0*dy/jmax
  yoff=72.0*y

  write(lu,*) '1.416 setmiterlimit'
  write(lu,1002) xoff+xscale*(x2-p(1)),yoff+0.5*yscale
1002 format('newpath',2f7.1,' moveto')

  do j=2,jmax
     write(lu,1004) xoff+xscale*(x2-p(j)),yoff+(j-0.5)*yscale
1004 format(2f6.1,' lt')
  enddo
  write(lu,*) 'stroke'

  if(nbox.ne.0) then
     write(lu,1010) xoff,yoff,72.0*dx,72.0*dy,-72*dx
1010 format('newpath',2f7.1,' moveto',f7.1,' 0 rlineto 0',            &
          f7.1,' rlineto',f7.1,' 0 rlineto closepath stroke')
  endif

  return
end subroutine imr4prov


code:
subroutine indexx(arr,n,indx)

  parameter (M=7,NSTACK=50)
  integer n,indx(n)
  real arr(n)
  integer i,indxt,ir,itemp,j,jstack,k,l,istack(NSTACK)
  real a

  do j=1,n
     indx(j)=j
  enddo

  jstack=0
  l=1
  ir=n
1 if(ir-l.lt.M) then
     do j=l+1,ir
        indxt=indx(j)
        a=arr(indxt)
        do i=j-1,1,-1
           if(arr(indx(i)).le.a) goto 2
           indx(i+1)=indx(i)
        enddo
        i=0
2       indx(i+1)=indxt
     enddo
     if(jstack.eq.0) return

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

  else
     k=(l+ir)/2
     itemp=indx(k)
     indx(k)=indx(l+1)
     indx(l+1)=itemp

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

     if(arr(indx(l)).gt.arr(indx(ir))) then
        itemp=indx(l)
        indx(l)=indx(ir)
        indx(ir)=itemp
     endif

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

     i=l+1
     j=ir
     indxt=indx(l)
     a=arr(indxt)
3    continue
     i=i+1
     if(arr(indx(i)).lt.a) goto 3

4    continue
     j=j-1
     if(arr(indx(j)).gt.a) goto 4
     if(j.lt.i) goto 5
     itemp=indx(i)
     indx(i)=indx(j)
     indx(j)=itemp
     goto 3

5    indx(l)=indx(j)
     indx(j)=indxt
     jstack=jstack+2
     if(jstack.gt.NSTACK) stop 'NSTACK too small in indexx'
     if(ir-i+1.ge.j-l)then
        istack(jstack)=ir
        istack(jstack-1)=i
        ir=j-1
     else
        istack(jstack)=j-1
        istack(jstack-1)=l
        l=i
     endif
  endif
  goto 1

end subroutine indexx


code:
#include "init_random_seed.h"

#include <stdlib.h>
#include <stdint.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <sys/time.h>
#include <unistd.h>

/* basic PRNG to use for improving the basic seed selection */
static unsigned lcg (uint64_t seed)
{
  if (0 ==seed)
    {
      seed = UINT64_C(104729);
    }
  else
    {
      seed %= UINT64_C(4294967296);
    }
  seed = (seed * UINT64_C(279470273)) % UINT64_C(4294967291);
  return seed % UINT64_MAX;
}

/* Generate a good PRNG seed value */
void init_random_seed(void)
{
  unsigned seed = 0u;
  int have_seed = 0;

  // try /dev/urandom for an initial seed
  int random_source;
  if ((random_source = open ("/dev/urandom", O_RDONLY)) >= 0)
    {
      size_t random_data_length = 0;
      have_seed = -1;
      while (random_data_length < sizeof seed)
        {
          ssize_t result = read (random_source, &seed + random_data_length, (sizeof seed) - random_data_length);
          if (result < 0)
            {
              // error, unable to read /dev/random
              have_seed = 0;
            }
          random_data_length += result;
        }
      close (random_source);
    }
  if (!have_seed)
    {
      // fallback to combining the time and PID in a fairly random way
      pid_t pid = getpid ();
      struct timeval tv;
      gettimeofday (&tv, NULL);
      seed = (unsigned)(((unsigned)pid << 16)
        ^ (unsigned)pid
        ^ (unsigned)tv.tv_sec
        ^ (unsigned)tv.tv_usec);
      seed = lcg (seed);
    }
  srand (seed);
}

#ifdef TEST
#include <stdio.h>

int main (int argc, char * argv[])
{
  init_random_seed ();
  int i, j;
  int r[10][4];
  for (i = 0; i < 10; ++i)
    {
      for (j = 0; j < 4; ++j)
        {
          printf ("%10d ", rand ());
        }
      printf ("\n");
    }
  return 0;
}
#endif


code:
!
! Generate a seed for the RANDOM_NUMBER PRNG that is guaranteed to be
! unique even if many processes are started simultaneously
!
subroutine init_random_seed()
  use iso_fortran_env, only: int64
  implicit none
  integer, allocatable :: seed(:)
  integer :: i, n, un, istat, dt(8), pid
  integer(int64) :: t

  call random_seed(size = n)
  allocate(seed(n))
  ! First try if the OS provides a random number generator
  open(newunit=un, file="/dev/urandom", access="stream", &
       form="unformatted", action="read", status="old", iostat=istat)
  if (istat == 0) then
     read(un) seed
     close(un)
  else
     ! Fallback to XOR:ing the current time and pid. The PID is
     ! useful in case one launches multiple instances of the same
     ! program in parallel.
     call system_clock(t)
     if (t == 0) then
        call date_and_time(values=dt)
        t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
             + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
             + dt(3) * 24_int64 * 60 * 60 * 1000 &
             + dt(5) * 60 * 60 * 1000 &
             + dt(6) * 60 * 1000 + dt(7) * 1000 &
             + dt(8)
     end if
     pid = getpid()
     t = ieor(t, int(pid, kind(t)))
     do i = 1, n
        seed(i) = lcg(t)
     end do
  end if
  call random_seed(put=seed)
contains
  ! This simple PRNG might not be good enough for real work, but is
  ! sufficient for seeding a better PRNG.
  function lcg(s)
    integer :: lcg
    integer(int64) :: s
    if (s == 0) then
       s = 104729
    else
       s = mod(s, 4294967296_int64)
    end if
    s = mod(s * 279470273_int64, 4294967291_int64)
    lcg = int(mod(s, int(huge(0), int64)), kind(0))
  end function lcg
end subroutine init_random_seed


code:
#ifndef INIT_RANDOM_SEED_H__
#define INIT_RANDOM_SEED_H__

#ifdef __cplusplus
extern "C" {
#endif

  /*
   * Generate a seed for the RANDOM_NUMBER PRNG that is guaranteed to
   * be unique even if many processes are started simultaneously
   *
   * Not suitable for multi-threaded requirements
   */
  void init_random_seed (void);

#ifdef __cplusplus
}
#endif

#endif


code:
/* Initialize a RS codec
 *
 * Copyright 2002 Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#include <stdlib.h>

#ifdef CCSDS
#include "ccsds.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

void FREE_RS(void *p){
  struct rs *rs = (struct rs *)p;

  free(rs->alpha_to);
  free(rs->index_of);
  free(rs->genpoly);
  free(rs);
}

/* Initialize a Reed-Solomon codec
 * symsize = symbol size, bits (1-8)
 * gfpoly = Field generator polynomial coefficients
 * fcr = first root of RS code generator polynomial, index form
 * prim = primitive element to generate polynomial roots
 * nroots = RS code generator polynomial degree (number of roots)
 * pad = padding bytes at front of shortened block
 */
void *INIT_RS(int symsize,int gfpoly,int fcr,int prim,
	int nroots,int pad){
  struct rs *rs;
  int i, j, sr,root,iprim;

  /* Check parameter ranges */
  if(symsize < 0 || symsize > (int)(8*sizeof(DTYPE)))
    return NULL; /* Need version with ints rather than chars */

  if(fcr < 0 || fcr >= (1<<symsize))
    return NULL;
  if(prim <= 0 || prim >= (1<<symsize))
    return NULL;
  if(nroots < 0 || nroots >= (1<<symsize))
    return NULL; /* Can't have more roots than symbol values! */
  if(pad < 0 || pad >= ((1<<symsize) -1 - nroots))
    return NULL; /* Too much padding */

  rs = (struct rs *)calloc(1,sizeof(struct rs));
  rs->mm = symsize;
  rs->nn = (1<<symsize)-1;
  rs->pad = pad;

  rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
  if(rs->alpha_to == NULL){
    free(rs);
    return NULL;
  }
  rs->index_of = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
  if(rs->index_of == NULL){
    free(rs->alpha_to);
    free(rs);
    return NULL;
  }

  /* Generate Galois field lookup tables */
  rs->index_of[0] = A0; /* log(zero) = -inf */
  rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */
  sr = 1;
  for(i=0;i<rs->nn;i++){
    rs->index_of[sr] = i;
    rs->alpha_to[i] = sr;
    sr <<= 1;
    if(sr & (1<<symsize))
      sr ^= gfpoly;
    sr &= rs->nn;
  }
  if(sr != 1){
    /* field generator polynomial is not primitive! */
    free(rs->alpha_to);
    free(rs->index_of);
    free(rs);
    return NULL;
  }

  /* Form RS code generator polynomial from its roots */
  rs->genpoly = (DTYPE *)malloc(sizeof(DTYPE)*(nroots+1));
  if(rs->genpoly == NULL){
    free(rs->alpha_to);
    free(rs->index_of);
    free(rs);
    return NULL;
  }
  rs->fcr = fcr;
  rs->prim = prim;
  rs->nroots = nroots;

  /* Find prim-th root of 1, used in decoding */
  for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
    ;
  rs->iprim = iprim / prim;

  rs->genpoly[0] = 1;
  for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) {
    rs->genpoly[i+1] = 1;

    /* Multiply rs->genpoly[] by  @**(root + x) */
    for (j = i; j > 0; j--){
      if (rs->genpoly[j] != 0)
	rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)];
      else
	rs->genpoly[j] = rs->genpoly[j-1];
    }
    /* rs->genpoly[0] can never be zero */
    rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)];
  }
  /* convert rs->genpoly[] to index form for quicker encoding */
  for (i = 0; i <= nroots; i++)
    rs->genpoly[i] = rs->index_of[rs->genpoly[i]];

  return rs;
}


code:
/* Include file to configure the RS codec for integer symbols
 *
 * Copyright 2002, Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#define DTYPE int

/* Reed-Solomon codec control block */
struct rs {
  int mm;              /* Bits per symbol */
  int nn;              /* Symbols per block (= (1<<mm)-1) */
  DTYPE *alpha_to;     /* log lookup table */
  DTYPE *index_of;     /* Antilog lookup table */
  DTYPE *genpoly;      /* Generator polynomial */
  int nroots;     /* Number of generator roots = number of parity symbols */
  int fcr;        /* First consecutive root, index form */
  int prim;       /* Primitive element, index form */
  int iprim;      /* prim-th root of 1, index form */
  int pad;        /* Padding bytes in shortened block */
};

static int modnn(struct rs *rs,int x){
  while (x >= rs->nn) {
    x -= rs->nn;
    x = (x >> rs->mm) + (x & rs->nn);
  }
  return x;
}
#define MODNN(x) modnn(rs,x)

#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
//#define NROOTS (rs->nroots)
#define NROOTS (51)
#define FCR (rs->fcr)
#define PRIM (rs->prim)
#define IPRIM (rs->iprim)
#define PAD (rs->pad)
#define A0 (NN)

#define ENCODE_RS encode_rs_int
#define DECODE_RS decode_rs_int
#define INIT_RS init_rs_int
#define FREE_RS free_rs_int

void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras);
void *INIT_RS(int symsize,int gfpoly,int fcr,
		   int prim,int nroots,int pad);
void FREE_RS(void *p);


code:
subroutine inter_wspr(id,ndir)

! Interleave (ndir=1) or de-interleave (ndir=-1) the array id.

  integer*1 id(0:161),itmp(0:161)
  integer j0(0:161)
  logical first
  data first/.true./
  save

  if(first) then
! Compute the interleave table using bit reversal.
     k=-1
     do i=0,255
        n=0
        ii=i
        do j=0,7
           n=n+n
           if(iand(ii,1).ne.0) n=n+1
           ii=ii/2
        enddo
        if(n.le.161) then
           k=k+1
           j0(k)=n
        endif
     enddo
     first=.false.
  endif

  if(ndir.eq.1) then
     do i=0,161
        itmp(j0(i))=id(i)
     enddo
  else
     do i=0,161
        itmp(i)=id(j0(i))
     enddo
  endif

  do i=0,161
     id(i)=itmp(i)
  enddo

  return
end subroutine inter_wspr


code:
subroutine interleave4(id,ndir)
  integer*1 id(0:205),itmp(0:205)
  integer j0(0:205)
  logical first
  data first/.true./
  save first,j0

  if(first) then
     k=-1
     do i=0,255
        m=i
        n=iand(m,1)
        n=2*n + iand(m/2,1)
        n=2*n + iand(m/4,1)
        n=2*n + iand(m/8,1)
        n=2*n + iand(m/16,1)
        n=2*n + iand(m/32,1)
        n=2*n + iand(m/64,1)
        n=2*n + iand(m/128,1)
        if(n.le.205) then
           k=k+1
           j0(k)=n
        endif
     enddo
     first=.false.
  endif

  if(ndir.eq.1) then
     do i=0,205
        itmp(j0(i))=id(i)
     enddo
  else
     do i=0,205
        itmp(i)=id(j0(i))
     enddo
  endif

  do i=0,205
     id(i)=itmp(i)
  enddo

  return
end subroutine interleave4


code:
subroutine interleave63(d1,idir)

! Interleave (idir=1) or de-interleave (idir=-1) the array d1.

  integer d1(0:6,0:8)
  integer d2(0:8,0:6)

  if(idir.ge.0) then
     do i=0,6
        do j=0,8
           d2(j,i)=d1(i,j)
        enddo
     enddo
     call move(d2,d1,63)
  else
     call move(d1,d2,63)
     do i=0,6
        do j=0,8
           d1(i,j)=d2(j,i)
        enddo
     enddo
  endif

  return
end subroutine interleave63


code:
subroutine interleave9(ia,ndir,ib)
  integer*1 ia(0:205),ib(0:205)
  integer j0(0:205)
  logical first
  data first/.true./
  save first,j0                     !Save not working, or j0 overwritten ???

  if(first) then
     k=-1
     do i=0,255
        m=i
        n=iand(m,1)
        n=2*n + iand(m/2,1)
        n=2*n + iand(m/4,1)
        n=2*n + iand(m/8,1)
        n=2*n + iand(m/16,1)
        n=2*n + iand(m/32,1)
        n=2*n + iand(m/64,1)
        n=2*n + iand(m/128,1)
        if(n.le.205) then
           k=k+1
           j0(k)=n
        endif
     enddo
!     first=.false.
  endif

  if(ndir.gt.0) then
     do i=0,205
        ib(j0(i))=ia(i)
     enddo
  else
     do i=0,205
        ib(i)=ia(j0(i))
     enddo
  endif

  return
end subroutine interleave9


code:
#include <QDebug>
#include <QString>
#include <QSharedMemory>
#include <QSystemSemaphore>

#include "../commons.h"

// Multiple instances: KK1D, 17 Jul 2013
QSharedMemory mem_jt9;

// Semaphore not changed, as the acquire/release calls do not
// appear to be used anywhere.
QSystemSemaphore sem_jt9("sem_jt9", 1, QSystemSemaphore::Open);

extern "C" {
  bool attach_jt9_();
  bool create_jt9_(int nsize);
  bool detach_jt9_();
  bool lock_jt9_();
  bool unlock_jt9_();
  struct jt9com * address_jt9_();
  int size_jt9_();
// Multiple instances:  wrapper for QSharedMemory::setKey()
  bool setkey_jt9_(char* mykey, int mykey_len);

  bool acquire_jt9_();
  bool release_jt9_();
}

bool attach_jt9_() {return mem_jt9.attach();}
bool create_jt9_(int nsize) {return mem_jt9.create(nsize);}
bool detach_jt9_() {return mem_jt9.detach();}
bool lock_jt9_() {return mem_jt9.lock();}
bool unlock_jt9_() {return mem_jt9.unlock();}
struct jt9com * address_jt9_() {return reinterpret_cast<struct jt9com *>(mem_jt9.data());}
int size_jt9_() {return (int)mem_jt9.size();}

// Multiple instances:
bool setkey_jt9_(char* mykey, int mykey_len) {
   char *tempstr = (char *)calloc(mykey_len+1,1);
   memset(tempstr, 0, mykey_len+1);
   strncpy(tempstr, mykey, mykey_len);
   QString s1 = QString(QLatin1String(tempstr));
   mem_jt9.setKey(s1);
   return true;}

bool acquire_jt9_() {return sem_jt9.acquire();}
bool release_jt9_() {return sem_jt9.release();}


code:
subroutine iscat(cdat0,npts0,nh,npct,t2,pick,cfile6,minsync,ntol,   &
     NFreeze,MouseDF,mousebutton,mode4,nafc,nmore,psavg,maxlines,nlines,line)

! Decode an ISCAT signal

  parameter (NMAX=30*3101)
  parameter (NSZ=4*1400)
  character cfile6*6                      !File time
  character c42*42
  character msg*29,msg1*29,msgbig*29
  character*80 line(100)
  character csync*1
  complex cdat0(NMAX)
  complex cdat(NMAX)
  real s0(288,NSZ)
  real fs1(0:41,30)
  real psavg(72)                          !Average spectrum of whole file
  integer nsum(30)
  integer ntol
  integer icos(4)
  logical pick,last
  data icos/0,1,3,2/
  data nsync/4/,nlen/2/,ndat/18/
  data c42/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ /.?@-'/
  save cdat,s0

  nlines = 0
  fsample=3100.78125                   !Sample rate after 9/32 downsampling
  nsps=144/mode4

  bigworst=-1.e30                      !Silence compiler warnings ...
  bigxsync=0.
  bigsig=-1.e30
  msglenbig=0
  ndf0big=0
  nfdotbig=0
  bigt2=0.
  bigavg=0.
  bigtana=0.
  if(nmore.eq.-999) bigsig=-1         !... to here

  last=.false.
  do inf=1,6                           !Loop over data-segment sizes
     nframes=2**inf
     if(nframes*24*nsps.gt.npts0) then
        nframes=npts0/(24*nsps)
        last=.true.
     endif
     npts=nframes*24*nsps

     do ia=1,npts0-npts,nsps*24        !Loop over start times stepped by 1 frame
        ib=ia+npts-1
        cdat(1:npts)=cdat0(ia:ib)
        t3=(ia + 0.5*npts)/fsample + 0.9
        if(pick) t3=t2+t3

! Compute symbol spectra and establish sync:
        call synciscat(cdat,npts,nh,npct,s0,jsym,df,ntol,NFreeze,     &
             MouseDF,mousebutton,mode4,nafc,psavg,xsync,sig,ndf0,msglen,    &
             ipk,jpk,idf,df1)
        nfdot=nint(idf*df1)

        isync=xsync
        if(msglen.eq.0 .or. isync.lt.max(minsync,0)) then
           msglen=0
           worst=1.
           avg=1.
           ndf0=0
           cycle
        endif

        ipk3=0                                  !Silence compiler warning
        nblk=nsync+nlen+ndat
        fs1=0.
        nsum=0
        nfold=jsym/96
        jb=96*nfold
        k=0
        n=0
        do j=jpk,jsym,4                !Fold information symbols into fs1
           k=k+1
           km=mod(k-1,nblk)+1
           if(km.gt.6) then
              n=n+1
              m=mod(n-1,msglen)+1
              ii=nint(idf*float(j-jb/2)/float(jb))
              do i=0,41
                 iii=ii+ipk+2*i
                 if(iii.ge.1 .and. iii.le.288) fs1(i,m)=fs1(i,m) + s0(iii,j)
              enddo
              nsum(m)=nsum(m)+1
           endif
        enddo

        do m=1,msglen
           fs1(0:41,m)=fs1(0:41,m)/nsum(m)
        enddo

! Read out the message contents:
        msg= '                             '
        msg1='                             '
        mpk=0
        worst=9999.
        sum=0.
        do m=1,msglen
           smax=0.
           smax2=0.
           do i=0,41
              if(fs1(i,m).gt.smax) then
                 smax=fs1(i,m)
                 ipk3=i
              endif
           enddo
           do i=0,41
              if(fs1(i,m).gt.smax2 .and. i.ne.ipk3) smax2=fs1(i,m)
           enddo
           rr=0.
           if(smax2.gt.0.0) rr=smax/smax2
           sum=sum + rr
           if(rr.lt.worst) worst=rr
           if(ipk3.eq.40) mpk=m
           msg1(m:m)=c42(ipk3+1:ipk3+1)
        enddo

        avg=sum/msglen
        if(mpk.eq.1) then
           msg=msg1(2:)
        else if(mpk.lt.msglen) then
           msg=msg1(mpk+1:msglen)//msg1(1:mpk-1)
        else
           msg=msg1(1:msglen-1)
        endif

        ttot=npts/3100.78125

        if(worst.gt.bigworst) then
           bigworst=worst
           bigavg=avg
           bigxsync=xsync
           bigsig=sig
           ndf0big=ndf0
           nfdotbig=nfdot
           msgbig=msg
           msglenbig=msglen
           bigt2=t3
           bigtana=nframes*24*nsps/fsample
        endif

        isync = xsync
        if(avg.gt.2.5 .and. xsync.ge.max(float(minsync),1.5) .and. &
             maxlines.ge.2) then
           nsig=nint(sig)
           nworst=10.0*(worst-1.0)
           navg=10.0*(avg-1.0)
           if(nworst.gt.10) nworst=10
           if(navg.gt.10) navg=10
           tana=nframes*24*nsps/fsample
           csync=' '
           if(isync.ge.1) csync='*'
           if(nlines.le.maxlines-1) nlines = nlines + 1
           write(line(nlines),1020) cfile6,isync,nsig,t2,ndf0,nfdot,csync, &
                msg(1:28),msglen,navg,nworst,tana,char(0)
        endif
     enddo
     if(last) exit
  enddo

  worst=bigworst
  avg=bigavg
  xsync=bigxsync
  sig=bigsig
  ndf0=ndf0big
  nfdot=nfdotbig
  msg=msgbig
  msglen=msglenbig
  t2=bigt2
  tana=bigtana

  isync=xsync
  nworst=10.0*(worst-1.0)
  navg=10.0*(avg-1.0)
  if(nworst.gt.10) nworst=10
  if(navg.gt.10) navg=10

  if(navg.le.0 .or. isync.lt.max(minsync,0)) then
     msg='                             '
     nworst=0
     navg=0
     ndf0=0
     nfdot=0
     sig=-20
     msglen=0
     tana=0.
     t2=0.
  endif
  csync=' '
  if(isync.ge.1) csync='*'
  nsig=nint(sig)

  if(nlines.le.maxlines-1) nlines = nlines + 1
  write(line(nlines),1020) cfile6,isync,nsig,t2,ndf0,nfdot,csync,msg(1:28),  &
       msglen,navg,nworst,tana,char(0)
1020 format(a6,2i4,f5.1,i5,i4,1x,a1,2x,a28,i4,i3,2x,i1,f5.1,a1)

  return
end subroutine iscat


code:
module iso_c_utilities

  use, intrinsic :: iso_c_binding, only: c_ptr, c_char, c_f_pointer, c_associated
  implicit none

  public :: c_to_f_string, c_f_dyn_string

  private

  character(c_char), dimension(1), save, target :: dummy_string = "?"

  interface        ! strlen is a standard C function from <string.h>
     ! int strlen(char *string)
     function strlen (string) result (len) bind (c, name="strlen")
       use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
       implicit none
       type(c_ptr), value :: string
       integer(kind=c_size_t) :: len
     end function strlen

     ! void free(void * p)
     subroutine c_free (p) bind (c, name="free")
       use, intrinsic :: iso_c_binding, only: c_ptr
       implicit none
       type(c_ptr), value :: p
     end subroutine c_free
  end interface

contains

  !
  ! Cast C string pointer to Fortran string pointer
  !
  ! Warning! - C data must outlive result scope
  !
  function c_to_f_string (c_str) result (f_str)
    use, intrinsic :: iso_c_binding, only: c_ptr, c_f_pointer, c_char
    implicit none
    type(c_ptr), intent(in) :: c_str
    character(kind=c_char, len=:), pointer :: f_str
    character(kind=c_char), pointer :: arr(:)
    interface        ! strlen is a standard C function from <string.h>
       ! int strlen(char *string)
       function strlen (string) result (len) bind (c, name="strlen")
         use, intrinsic :: iso_c_binding, only: c_ptr, c_size_t
         implicit none
         type(c_ptr), value :: string
         integer(kind=c_size_t) :: len
       end function strlen
    end interface
    call c_f_pointer (c_str, arr, [strlen (c_str)])
    call get_scalar_pointer (size (arr), arr, f_str)
  end function c_to_f_string

  subroutine get_scalar_pointer (scalar_len, scalar, fptr)
    ! Convert a null-terminated C string into a Fortran character pointer
    use, intrinsic :: iso_c_binding, only: c_char
    integer, intent(in) :: scalar_len
    character(kind=c_char, len=scalar_len), intent(in), target :: scalar(1)
    character(kind=c_char, len=:), pointer :: fptr
    fptr => scalar(1)
  end subroutine get_scalar_pointer

  function c_f_dyn_string (cptr) result (fstr)
    ! Convert a null-terminated malloc'ed C string into a Fortran character array
    type(c_ptr), intent(in) :: cptr ! The C address
    character(kind=c_char), allocatable :: fstr(:)
    character(kind=c_char), pointer :: fptr(:)
    interface        ! strlen is a standard C function from <string.h>
       ! void free(void * p)
       subroutine c_free (p) bind (c, name="free")
         use, intrinsic :: iso_c_binding, only: c_ptr
         implicit none
         type(c_ptr), value :: p
       end subroutine c_free
    end interface
    if (c_associated (cptr)) then
       call c_f_pointer (fptr=fptr, cptr=cptr, shape=[strlen(cptr)])
    else
       ! To avoid segfaults, associate FPTR with a dummy target:
       fptr => dummy_string
    end if
    fstr = fptr
    call c_free (cptr)
  end function c_f_dyn_string

end module