code:
// qra13_64_64_irr_e.c
// Encoding/Decoding tables for Q-ary RA code (13,64) over GF(64)
// Code Name: qra13_64_64_irr_e
// (13,64) RA Code over GF(64) RF=[3x4 4x4 6x1 3x2 5x1 7x1]/18

// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy

// 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.
//
//    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/>.

#include "qra13_64_64_irr_e.h"

#define qra_K       13 // number of information symbols
#define qra_N       64 // codeword length in symbols
#define qra_m        6 // bits/symbol
#define qra_M       64 // Symbol alphabet cardinality
#define qra_a		 1 // grouping factor
#define qra_NC      51 // number of check symbols (N-K)

// Defines used by the message passing decoder --------

#define qra_V       64 // number of variables in the code graph (N)
#define qra_C      116 // number of factors in the code graph (N +(N-K)+1)
#define qra_NMSG   218 // number of msgs in the code graph
#define qra_MAXVDEG    8   // maximum variable degree
#define qra_MAXCDEG    3   // maximum factor degree
#define qra_R     0.20313f // code rate (K/N)
#define CODE_NAME "qra_13_64_64_irr_e"

// table of the systematic symbols indexes in the accumulator chain
static const int qra_acc_input_idx[qra_NC+1] = {
 12,   4,   3,   9,   0,  11,   6,   8,  12,   1,
  2,   7,   4,  11,   3,   5,   9,   8,  12,   7,
  2,   4,  10,   3,   5,  11,  12,   8,   9,   6,
  7,   2,   5,   4,  12,   8,  11,   1,   6,   7,
  0,  10,  12,   8,  11,   5,   6,   1,   0,  10,
 12,   8
};

// table of the systematic symbols weight logarithms over GF(M)
static const  int qra_acc_input_wlog[qra_NC+1] = {
  0,  27,   0,   0,   0,  31,  28,  61,  31,   0,
  0,  52,  22,   7,  19,  47,  44,  62,  32,  50,
 52,  42,  48,  56,  40,  50,  51,  37,  37,   0,
  5,  14,   0,   0,  18,   2,   0,  45,  21,   0,
 62,   8,  11,  60,  36,  32,  17,   9,   5,   0,
 53,   0
};

// table of the logarithms of the elements of GF(M) (log(0) never used)
static const int qra_log[qra_M] = {
 -1,   0,   1,   6,   2,  12,   7,  26,   3,  32,
 13,  35,   8,  48,  27,  18,   4,  24,  33,  16,
 14,  52,  36,  54,   9,  45,  49,  38,  28,  41,
 19,  56,   5,  62,  25,  11,  34,  31,  17,  47,
 15,  23,  53,  51,  37,  44,  55,  40,  10,  61,
 46,  30,  50,  22,  39,  43,  29,  60,  42,  21,
 20,  59,  57,  58
};

// table of GF(M) elements given their logarithm
static const  int qra_exp[qra_M-1] = {
  1,   2,   4,   8,  16,  32,   3,   6,  12,  24,
 48,  35,   5,  10,  20,  40,  19,  38,  15,  30,
 60,  59,  53,  41,  17,  34,   7,  14,  28,  56,
 51,  37,   9,  18,  36,  11,  22,  44,  27,  54,
 47,  29,  58,  55,  45,  25,  50,  39,  13,  26,
 52,  43,  21,  42,  23,  46,  31,  62,  63,  61,
 57,  49,  33
};

// table of the messages weight logarithms over GF(M)
static const  int qra_msgw[qra_NMSG] = {
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,  27,   0,   0,   0,  31,
 28,  61,  31,   0,   0,  52,  22,   7,  19,  47,
 44,  62,  32,  50,  52,  42,  48,  56,  40,  50,
 51,  37,  37,   0,   5,  14,   0,   0,  18,   2,
  0,  45,  21,   0,  62,   8,  11,  60,  36,  32,
 17,   9,   5,   0,  53,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
  0,   0,   0,   0,   0,   0,   0,   0
};

// table of the degrees of the variable nodes
static const  int qra_vdeg[qra_V] = {
  4,   4,   4,   4,   5,   5,   5,   5,   7,   4,
  4,   6,   8,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3
};

// table of the degrees of the factor nodes
static const  int qra_cdeg[qra_C] = {
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   1,   1,   1,   1,   1,   1,
  1,   1,   1,   1,   2,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   3,   3,   3,   3,   3,
  3,   3,   3,   3,   3,   2
};

// table (uncompressed) of the v->c message indexes (-1=unused entry)
static const int qra_v2cmidx[qra_V*qra_MAXVDEG] = {
  0,  68, 104, 112,  -1,  -1,  -1,  -1,
  1,  73, 101, 111,  -1,  -1,  -1,  -1,
  2,  74,  84,  95,  -1,  -1,  -1,  -1,
  3,  66,  78,  87,  -1,  -1,  -1,  -1,
  4,  65,  76,  85,  97,  -1,  -1,  -1,
  5,  79,  88,  96, 109,  -1,  -1,  -1,
  6,  70,  93, 102, 110,  -1,  -1,  -1,
  7,  75,  83,  94, 103,  -1,  -1,  -1,
  8,  71,  81,  91,  99, 107, 115,  -1,
  9,  67,  80,  92,  -1,  -1,  -1,  -1,
 10,  86, 105, 113,  -1,  -1,  -1,  -1,
 11,  69,  77,  89, 100, 108,  -1,  -1,
 12,  64,  72,  82,  90,  98, 106, 114,
 13, 116, 117,  -1,  -1,  -1,  -1,  -1,
 14, 118, 119,  -1,  -1,  -1,  -1,  -1,
 15, 120, 121,  -1,  -1,  -1,  -1,  -1,
 16, 122, 123,  -1,  -1,  -1,  -1,  -1,
 17, 124, 125,  -1,  -1,  -1,  -1,  -1,
 18, 126, 127,  -1,  -1,  -1,  -1,  -1,
 19, 128, 129,  -1,  -1,  -1,  -1,  -1,
 20, 130, 131,  -1,  -1,  -1,  -1,  -1,
 21, 132, 133,  -1,  -1,  -1,  -1,  -1,
 22, 134, 135,  -1,  -1,  -1,  -1,  -1,
 23, 136, 137,  -1,  -1,  -1,  -1,  -1,
 24, 138, 139,  -1,  -1,  -1,  -1,  -1,
 25, 140, 141,  -1,  -1,  -1,  -1,  -1,
 26, 142, 143,  -1,  -1,  -1,  -1,  -1,
 27, 144, 145,  -1,  -1,  -1,  -1,  -1,
 28, 146, 147,  -1,  -1,  -1,  -1,  -1,
 29, 148, 149,  -1,  -1,  -1,  -1,  -1,
 30, 150, 151,  -1,  -1,  -1,  -1,  -1,
 31, 152, 153,  -1,  -1,  -1,  -1,  -1,
 32, 154, 155,  -1,  -1,  -1,  -1,  -1,
 33, 156, 157,  -1,  -1,  -1,  -1,  -1,
 34, 158, 159,  -1,  -1,  -1,  -1,  -1,
 35, 160, 161,  -1,  -1,  -1,  -1,  -1,
 36, 162, 163,  -1,  -1,  -1,  -1,  -1,
 37, 164, 165,  -1,  -1,  -1,  -1,  -1,
 38, 166, 167,  -1,  -1,  -1,  -1,  -1,
 39, 168, 169,  -1,  -1,  -1,  -1,  -1,
 40, 170, 171,  -1,  -1,  -1,  -1,  -1,
 41, 172, 173,  -1,  -1,  -1,  -1,  -1,
 42, 174, 175,  -1,  -1,  -1,  -1,  -1,
 43, 176, 177,  -1,  -1,  -1,  -1,  -1,
 44, 178, 179,  -1,  -1,  -1,  -1,  -1,
 45, 180, 181,  -1,  -1,  -1,  -1,  -1,
 46, 182, 183,  -1,  -1,  -1,  -1,  -1,
 47, 184, 185,  -1,  -1,  -1,  -1,  -1,
 48, 186, 187,  -1,  -1,  -1,  -1,  -1,
 49, 188, 189,  -1,  -1,  -1,  -1,  -1,
 50, 190, 191,  -1,  -1,  -1,  -1,  -1,
 51, 192, 193,  -1,  -1,  -1,  -1,  -1,
 52, 194, 195,  -1,  -1,  -1,  -1,  -1,
 53, 196, 197,  -1,  -1,  -1,  -1,  -1,
 54, 198, 199,  -1,  -1,  -1,  -1,  -1,
 55, 200, 201,  -1,  -1,  -1,  -1,  -1,
 56, 202, 203,  -1,  -1,  -1,  -1,  -1,
 57, 204, 205,  -1,  -1,  -1,  -1,  -1,
 58, 206, 207,  -1,  -1,  -1,  -1,  -1,
 59, 208, 209,  -1,  -1,  -1,  -1,  -1,
 60, 210, 211,  -1,  -1,  -1,  -1,  -1,
 61, 212, 213,  -1,  -1,  -1,  -1,  -1,
 62, 214, 215,  -1,  -1,  -1,  -1,  -1,
 63, 216, 217,  -1,  -1,  -1,  -1,  -1
};

// table (uncompressed) of the c->v message indexes (-1=unused entry)
static const int qra_c2vmidx[qra_C*qra_MAXCDEG] = {
  0,  -1,  -1,   1,  -1,  -1,   2,  -1,  -1,   3,  -1,  -1,
  4,  -1,  -1,   5,  -1,  -1,   6,  -1,  -1,   7,  -1,  -1,
  8,  -1,  -1,   9,  -1,  -1,  10,  -1,  -1,  11,  -1,  -1,
 12,  -1,  -1,  13,  -1,  -1,  14,  -1,  -1,  15,  -1,  -1,
 16,  -1,  -1,  17,  -1,  -1,  18,  -1,  -1,  19,  -1,  -1,
 20,  -1,  -1,  21,  -1,  -1,  22,  -1,  -1,  23,  -1,  -1,
 24,  -1,  -1,  25,  -1,  -1,  26,  -1,  -1,  27,  -1,  -1,
 28,  -1,  -1,  29,  -1,  -1,  30,  -1,  -1,  31,  -1,  -1,
 32,  -1,  -1,  33,  -1,  -1,  34,  -1,  -1,  35,  -1,  -1,
 36,  -1,  -1,  37,  -1,  -1,  38,  -1,  -1,  39,  -1,  -1,
 40,  -1,  -1,  41,  -1,  -1,  42,  -1,  -1,  43,  -1,  -1,
 44,  -1,  -1,  45,  -1,  -1,  46,  -1,  -1,  47,  -1,  -1,
 48,  -1,  -1,  49,  -1,  -1,  50,  -1,  -1,  51,  -1,  -1,
 52,  -1,  -1,  53,  -1,  -1,  54,  -1,  -1,  55,  -1,  -1,
 56,  -1,  -1,  57,  -1,  -1,  58,  -1,  -1,  59,  -1,  -1,
 60,  -1,  -1,  61,  -1,  -1,  62,  -1,  -1,  63,  -1,  -1,
 64, 116,  -1,  65, 117, 118,  66, 119, 120,  67, 121, 122,
 68, 123, 124,  69, 125, 126,  70, 127, 128,  71, 129, 130,
 72, 131, 132,  73, 133, 134,  74, 135, 136,  75, 137, 138,
 76, 139, 140,  77, 141, 142,  78, 143, 144,  79, 145, 146,
 80, 147, 148,  81, 149, 150,  82, 151, 152,  83, 153, 154,
 84, 155, 156,  85, 157, 158,  86, 159, 160,  87, 161, 162,
 88, 163, 164,  89, 165, 166,  90, 167, 168,  91, 169, 170,
 92, 171, 172,  93, 173, 174,  94, 175, 176,  95, 177, 178,
 96, 179, 180,  97, 181, 182,  98, 183, 184,  99, 185, 186,
100, 187, 188, 101, 189, 190, 102, 191, 192, 103, 193, 194,
104, 195, 196, 105, 197, 198, 106, 199, 200, 107, 201, 202,
108, 203, 204, 109, 205, 206, 110, 207, 208, 111, 209, 210,
112, 211, 212, 113, 213, 214, 114, 215, 216, 115, 217,  -1
};

// permutation matrix to compute Prob(x*alfa^logw)
static const int qra_pmat[qra_M*qra_M] = {
  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,  14,  15,
 16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,  30,  31,
 32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,  44,  45,  46,  47,
 48,  49,  50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,
  0,  33,   1,  32,   2,  35,   3,  34,   4,  37,   5,  36,   6,  39,   7,  38,
  8,  41,   9,  40,  10,  43,  11,  42,  12,  45,  13,  44,  14,  47,  15,  46,
 16,  49,  17,  48,  18,  51,  19,  50,  20,  53,  21,  52,  22,  55,  23,  54,
 24,  57,  25,  56,  26,  59,  27,  58,  28,  61,  29,  60,  30,  63,  31,  62,
  0,  49,  33,  16,   1,  48,  32,  17,   2,  51,  35,  18,   3,  50,  34,  19,
  4,  53,  37,  20,   5,  52,  36,  21,   6,  55,  39,  22,   7,  54,  38,  23,
  8,  57,  41,  24,   9,  56,  40,  25,  10,  59,  43,  26,  11,  58,  42,  27,
 12,  61,  45,  28,  13,  60,  44,  29,  14,  63,  47,  30,  15,  62,  46,  31,
  0,  57,  49,   8,  33,  24,  16,  41,   1,  56,  48,   9,  32,  25,  17,  40,
  2,  59,  51,  10,  35,  26,  18,  43,   3,  58,  50,  11,  34,  27,  19,  42,
  4,  61,  53,  12,  37,  28,  20,  45,   5,  60,  52,  13,  36,  29,  21,  44,
  6,  63,  55,  14,  39,  30,  22,  47,   7,  62,  54,  15,  38,  31,  23,  46,
  0,  61,  57,   4,  49,  12,   8,  53,  33,  28,  24,  37,  16,  45,  41,  20,
  1,  60,  56,   5,  48,  13,   9,  52,  32,  29,  25,  36,  17,  44,  40,  21,
  2,  63,  59,   6,  51,  14,  10,  55,  35,  30,  26,  39,  18,  47,  43,  22,
  3,  62,  58,   7,  50,  15,  11,  54,  34,  31,  27,  38,  19,  46,  42,  23,
  0,  63,  61,   2,  57,   6,   4,  59,  49,  14,  12,  51,   8,  55,  53,  10,
 33,  30,  28,  35,  24,  39,  37,  26,  16,  47,  45,  18,  41,  22,  20,  43,
  1,  62,  60,   3,  56,   7,   5,  58,  48,  15,  13,  50,   9,  54,  52,  11,
 32,  31,  29,  34,  25,  38,  36,  27,  17,  46,  44,  19,  40,  23,  21,  42,
  0,  62,  63,   1,  61,   3,   2,  60,  57,   7,   6,  56,   4,  58,  59,   5,
 49,  15,  14,  48,  12,  50,  51,  13,   8,  54,  55,   9,  53,  11,  10,  52,
 33,  31,  30,  32,  28,  34,  35,  29,  24,  38,  39,  25,  37,  27,  26,  36,
 16,  46,  47,  17,  45,  19,  18,  44,  41,  23,  22,  40,  20,  42,  43,  21,
  0,  31,  62,  33,  63,  32,   1,  30,  61,  34,   3,  28,   2,  29,  60,  35,
 57,  38,   7,  24,   6,  25,  56,  39,   4,  27,  58,  37,  59,  36,   5,  26,
 49,  46,  15,  16,  14,  17,  48,  47,  12,  19,  50,  45,  51,  44,  13,  18,
  8,  23,  54,  41,  55,  40,   9,  22,  53,  42,  11,  20,  10,  21,  52,  43,
  0,  46,  31,  49,  62,  16,  33,  15,  63,  17,  32,  14,   1,  47,  30,  48,
 61,  19,  34,  12,   3,  45,  28,  50,   2,  44,  29,  51,  60,  18,  35,  13,
 57,  23,  38,   8,   7,  41,  24,  54,   6,  40,  25,  55,  56,  22,  39,   9,
  4,  42,  27,  53,  58,  20,  37,  11,  59,  21,  36,  10,   5,  43,  26,  52,
  0,  23,  46,  57,  31,   8,  49,  38,  62,  41,  16,   7,  33,  54,  15,  24,
 63,  40,  17,   6,  32,  55,  14,  25,   1,  22,  47,  56,  30,   9,  48,  39,
 61,  42,  19,   4,  34,  53,  12,  27,   3,  20,  45,  58,  28,  11,  50,  37,
  2,  21,  44,  59,  29,  10,  51,  36,  60,  43,  18,   5,  35,  52,  13,  26,
  0,  42,  23,  61,  46,   4,  57,  19,  31,  53,   8,  34,  49,  27,  38,  12,
 62,  20,  41,   3,  16,  58,   7,  45,  33,  11,  54,  28,  15,  37,  24,  50,
 63,  21,  40,   2,  17,  59,   6,  44,  32,  10,  55,  29,  14,  36,  25,  51,
  1,  43,  22,  60,  47,   5,  56,  18,  30,  52,   9,  35,  48,  26,  39,  13,
  0,  21,  42,  63,  23,   2,  61,  40,  46,  59,   4,  17,  57,  44,  19,   6,
 31,  10,  53,  32,   8,  29,  34,  55,  49,  36,  27,  14,  38,  51,  12,  25,
 62,  43,  20,   1,  41,  60,   3,  22,  16,   5,  58,  47,   7,  18,  45,  56,
 33,  52,  11,  30,  54,  35,  28,   9,  15,  26,  37,  48,  24,  13,  50,  39,
  0,  43,  21,  62,  42,   1,  63,  20,  23,  60,   2,  41,  61,  22,  40,   3,
 46,   5,  59,  16,   4,  47,  17,  58,  57,  18,  44,   7,  19,  56,   6,  45,
 31,  52,  10,  33,  53,  30,  32,  11,   8,  35,  29,  54,  34,   9,  55,  28,
 49,  26,  36,  15,  27,  48,  14,  37,  38,  13,  51,  24,  12,  39,  25,  50,
  0,  52,  43,  31,  21,  33,  62,  10,  42,  30,   1,  53,  63,  11,  20,  32,
 23,  35,  60,   8,   2,  54,  41,  29,  61,   9,  22,  34,  40,  28,   3,  55,
 46,  26,   5,  49,  59,  15,  16,  36,   4,  48,  47,  27,  17,  37,  58,  14,
 57,  13,  18,  38,  44,  24,   7,  51,  19,  39,  56,  12,   6,  50,  45,  25,
  0,  26,  52,  46,  43,  49,  31,   5,  21,  15,  33,  59,  62,  36,  10,  16,
 42,  48,  30,   4,   1,  27,  53,  47,  63,  37,  11,  17,  20,  14,  32,  58,
 23,  13,  35,  57,  60,  38,   8,  18,   2,  24,  54,  44,  41,  51,  29,   7,
 61,  39,   9,  19,  22,  12,  34,  56,  40,  50,  28,   6,   3,  25,  55,  45,
  0,  13,  26,  23,  52,  57,  46,  35,  43,  38,  49,  60,  31,  18,   5,   8,
 21,  24,  15,   2,  33,  44,  59,  54,  62,  51,  36,  41,  10,   7,  16,  29,
 42,  39,  48,  61,  30,  19,   4,   9,   1,  12,  27,  22,  53,  56,  47,  34,
 63,  50,  37,  40,  11,   6,  17,  28,  20,  25,  14,   3,  32,  45,  58,  55,
  0,  39,  13,  42,  26,  61,  23,  48,  52,  19,  57,  30,  46,   9,  35,   4,
 43,  12,  38,   1,  49,  22,  60,  27,  31,  56,  18,  53,   5,  34,   8,  47,
 21,  50,  24,  63,  15,  40,   2,  37,  33,   6,  44,  11,  59,  28,  54,  17,
 62,  25,  51,  20,  36,   3,  41,  14,  10,  45,   7,  32,  16,  55,  29,  58,
  0,  50,  39,  21,  13,  63,  42,  24,  26,  40,  61,  15,  23,  37,  48,   2,
 52,   6,  19,  33,  57,  11,  30,  44,  46,  28,   9,  59,  35,  17,   4,  54,
 43,  25,  12,  62,  38,  20,   1,  51,  49,   3,  22,  36,  60,  14,  27,  41,
 31,  45,  56,  10,  18,  32,  53,   7,   5,  55,  34,  16,   8,  58,  47,  29,
  0,  25,  50,  43,  39,  62,  21,  12,  13,  20,  63,  38,  42,  51,  24,   1,
 26,   3,  40,  49,  61,  36,  15,  22,  23,  14,  37,  60,  48,  41,   2,  27,
 52,  45,   6,  31,  19,  10,  33,  56,  57,  32,  11,  18,  30,   7,  44,  53,
 46,  55,  28,   5,   9,  16,  59,  34,  35,  58,  17,   8,   4,  29,  54,  47,
  0,  45,  25,  52,  50,  31,  43,   6,  39,  10,  62,  19,  21,  56,  12,  33,
 13,  32,  20,  57,  63,  18,  38,  11,  42,   7,  51,  30,  24,  53,   1,  44,
 26,  55,   3,  46,  40,   5,  49,  28,  61,  16,  36,   9,  15,  34,  22,  59,
 23,  58,  14,  35,  37,   8,  60,  17,  48,  29,  41,   4,   2,  47,  27,  54,
  0,  55,  45,  26,  25,  46,  52,   3,  50,   5,  31,  40,  43,  28,   6,  49,
 39,  16,  10,  61,  62,   9,  19,  36,  21,  34,  56,  15,  12,  59,  33,  22,
 13,  58,  32,  23,  20,  35,  57,  14,  63,   8,  18,  37,  38,  17,  11,  60,
 42,  29,   7,  48,  51,   4,  30,  41,  24,  47,  53,   2,   1,  54,  44,  27,
  0,  58,  55,  13,  45,  23,  26,  32,  25,  35,  46,  20,  52,  14,   3,  57,
 50,   8,   5,  63,  31,  37,  40,  18,  43,  17,  28,  38,   6,  60,  49,  11,
 39,  29,  16,  42,  10,  48,  61,   7,  62,   4,   9,  51,  19,  41,  36,  30,
 21,  47,  34,  24,  56,   2,  15,  53,  12,  54,  59,   1,  33,  27,  22,  44,
  0,  29,  58,  39,  55,  42,  13,  16,  45,  48,  23,  10,  26,   7,  32,  61,
 25,   4,  35,  62,  46,  51,  20,   9,  52,  41,  14,  19,   3,  30,  57,  36,
 50,  47,   8,  21,   5,  24,  63,  34,  31,   2,  37,  56,  40,  53,  18,  15,
 43,  54,  17,  12,  28,   1,  38,  59,   6,  27,  60,  33,  49,  44,  11,  22,
  0,  47,  29,  50,  58,  21,  39,   8,  55,  24,  42,   5,  13,  34,  16,  63,
 45,   2,  48,  31,  23,  56,  10,  37,  26,  53,   7,  40,  32,  15,  61,  18,
 25,  54,   4,  43,  35,  12,  62,  17,  46,   1,  51,  28,  20,  59,   9,  38,
 52,  27,  41,   6,  14,  33,  19,  60,   3,  44,  30,  49,  57,  22,  36,  11,
  0,  54,  47,  25,  29,  43,  50,   4,  58,  12,  21,  35,  39,  17,   8,  62,
 55,   1,  24,  46,  42,  28,   5,  51,  13,  59,  34,  20,  16,  38,  63,   9,
 45,  27,   2,  52,  48,   6,  31,  41,  23,  33,  56,  14,  10,  60,  37,  19,
 26,  44,  53,   3,   7,  49,  40,  30,  32,  22,  15,  57,  61,  11,  18,  36,
  0,  27,  54,  45,  47,  52,  25,   2,  29,   6,  43,  48,  50,  41,   4,  31,
 58,  33,  12,  23,  21,  14,  35,  56,  39,  60,  17,  10,   8,  19,  62,  37,
 55,  44,   1,  26,  24,   3,  46,  53,  42,  49,  28,   7,   5,  30,  51,  40,
 13,  22,  59,  32,  34,  57,  20,  15,  16,  11,  38,  61,  63,  36,   9,  18,
  0,  44,  27,  55,  54,  26,  45,   1,  47,   3,  52,  24,  25,  53,   2,  46,
 29,  49,   6,  42,  43,   7,  48,  28,  50,  30,  41,   5,   4,  40,  31,  51,
 58,  22,  33,  13,  12,  32,  23,  59,  21,  57,  14,  34,  35,  15,  56,  20,
 39,  11,  60,  16,  17,  61,  10,  38,   8,  36,  19,  63,  62,  18,  37,   9,
  0,  22,  44,  58,  27,  13,  55,  33,  54,  32,  26,  12,  45,  59,   1,  23,
 47,  57,   3,  21,  52,  34,  24,  14,  25,  15,  53,  35,   2,  20,  46,  56,
 29,  11,  49,  39,   6,  16,  42,  60,  43,  61,   7,  17,  48,  38,  28,  10,
 50,  36,  30,   8,  41,  63,   5,  19,   4,  18,  40,  62,  31,   9,  51,  37,
  0,  11,  22,  29,  44,  39,  58,  49,  27,  16,  13,   6,  55,  60,  33,  42,
 54,  61,  32,  43,  26,  17,  12,   7,  45,  38,  59,  48,   1,  10,  23,  28,
 47,  36,  57,  50,   3,   8,  21,  30,  52,  63,  34,  41,  24,  19,  14,   5,
 25,  18,  15,   4,  53,  62,  35,  40,   2,   9,  20,  31,  46,  37,  56,  51,
  0,  36,  11,  47,  22,  50,  29,  57,  44,   8,  39,   3,  58,  30,  49,  21,
 27,  63,  16,  52,  13,  41,   6,  34,  55,  19,  60,  24,  33,   5,  42,  14,
 54,  18,  61,  25,  32,   4,  43,  15,  26,  62,  17,  53,  12,  40,   7,  35,
 45,   9,  38,   2,  59,  31,  48,  20,   1,  37,  10,  46,  23,  51,  28,  56,
  0,  18,  36,  54,  11,  25,  47,  61,  22,   4,  50,  32,  29,  15,  57,  43,
 44,  62,   8,  26,  39,  53,   3,  17,  58,  40,  30,  12,  49,  35,  21,   7,
 27,   9,  63,  45,  16,   2,  52,  38,  13,  31,  41,  59,   6,  20,  34,  48,
 55,  37,  19,   1,  60,  46,  24,  10,  33,  51,   5,  23,  42,  56,  14,  28,
  0,   9,  18,  27,  36,  45,  54,  63,  11,   2,  25,  16,  47,  38,  61,  52,
 22,  31,   4,  13,  50,  59,  32,  41,  29,  20,  15,   6,  57,  48,  43,  34,
 44,  37,  62,  55,   8,   1,  26,  19,  39,  46,  53,  60,   3,  10,  17,  24,
 58,  51,  40,  33,  30,  23,  12,   5,  49,  56,  35,  42,  21,  28,   7,  14,
  0,  37,   9,  44,  18,  55,  27,  62,  36,   1,  45,   8,  54,  19,  63,  26,
 11,  46,   2,  39,  25,  60,  16,  53,  47,  10,  38,   3,  61,  24,  52,  17,
 22,  51,  31,  58,   4,  33,  13,  40,  50,  23,  59,  30,  32,   5,  41,  12,
 29,  56,  20,  49,  15,  42,   6,  35,  57,  28,  48,  21,  43,  14,  34,   7,
  0,  51,  37,  22,   9,  58,  44,  31,  18,  33,  55,   4,  27,  40,  62,  13,
 36,  23,   1,  50,  45,  30,   8,  59,  54,   5,  19,  32,  63,  12,  26,  41,
 11,  56,  46,  29,   2,  49,  39,  20,  25,  42,  60,  15,  16,  35,  53,   6,
 47,  28,  10,  57,  38,  21,   3,  48,  61,  14,  24,  43,  52,   7,  17,  34,
  0,  56,  51,  11,  37,  29,  22,  46,   9,  49,  58,   2,  44,  20,  31,  39,
 18,  42,  33,  25,  55,  15,   4,  60,  27,  35,  40,  16,  62,   6,  13,  53,
 36,  28,  23,  47,   1,  57,  50,  10,  45,  21,  30,  38,   8,  48,  59,   3,
 54,  14,   5,  61,  19,  43,  32,  24,  63,   7,  12,  52,  26,  34,  41,  17,
  0,  28,  56,  36,  51,  47,  11,  23,  37,  57,  29,   1,  22,  10,  46,  50,
  9,  21,  49,  45,  58,  38,   2,  30,  44,  48,  20,   8,  31,   3,  39,  59,
 18,  14,  42,  54,  33,  61,  25,   5,  55,  43,  15,  19,   4,  24,  60,  32,
 27,   7,  35,  63,  40,  52,  16,  12,  62,  34,   6,  26,  13,  17,  53,  41,
  0,  14,  28,  18,  56,  54,  36,  42,  51,  61,  47,  33,  11,   5,  23,  25,
 37,  43,  57,  55,  29,  19,   1,  15,  22,  24,  10,   4,  46,  32,  50,  60,
  9,   7,  21,  27,  49,  63,  45,  35,  58,  52,  38,  40,   2,  12,  30,  16,
 44,  34,  48,  62,  20,  26,   8,   6,  31,  17,   3,  13,  39,  41,  59,  53,
  0,   7,  14,   9,  28,  27,  18,  21,  56,  63,  54,  49,  36,  35,  42,  45,
 51,  52,  61,  58,  47,  40,  33,  38,  11,  12,   5,   2,  23,  16,  25,  30,
 37,  34,  43,  44,  57,  62,  55,  48,  29,  26,  19,  20,   1,   6,  15,   8,
 22,  17,  24,  31,  10,  13,   4,   3,  46,  41,  32,  39,  50,  53,  60,  59,
  0,  34,   7,  37,  14,  44,   9,  43,  28,  62,  27,  57,  18,  48,  21,  55,
 56,  26,  63,  29,  54,  20,  49,  19,  36,   6,  35,   1,  42,   8,  45,  15,
 51,  17,  52,  22,  61,  31,  58,  24,  47,  13,  40,  10,  33,   3,  38,   4,
 11,  41,  12,  46,   5,  39,   2,  32,  23,  53,  16,  50,  25,  59,  30,  60,
  0,  17,  34,  51,   7,  22,  37,  52,  14,  31,  44,  61,   9,  24,  43,  58,
 28,  13,  62,  47,  27,  10,  57,  40,  18,   3,  48,  33,  21,   4,  55,  38,
 56,  41,  26,  11,  63,  46,  29,  12,  54,  39,  20,   5,  49,  32,  19,   2,
 36,  53,   6,  23,  35,  50,   1,  16,  42,  59,   8,  25,  45,  60,  15,  30,
  0,  41,  17,  56,  34,  11,  51,  26,   7,  46,  22,  63,  37,  12,  52,  29,
 14,  39,  31,  54,  44,   5,  61,  20,   9,  32,  24,  49,  43,   2,  58,  19,
 28,  53,  13,  36,  62,  23,  47,   6,  27,  50,  10,  35,  57,  16,  40,   1,
 18,  59,   3,  42,  48,  25,  33,   8,  21,  60,   4,  45,  55,  30,  38,  15,
  0,  53,  41,  28,  17,  36,  56,  13,  34,  23,  11,  62,  51,   6,  26,  47,
  7,  50,  46,  27,  22,  35,  63,  10,  37,  16,  12,  57,  52,   1,  29,  40,
 14,  59,  39,  18,  31,  42,  54,   3,  44,  25,   5,  48,  61,   8,  20,  33,
  9,  60,  32,  21,  24,  45,  49,   4,  43,  30,   2,  55,  58,  15,  19,  38,
  0,  59,  53,  14,  41,  18,  28,  39,  17,  42,  36,  31,  56,   3,  13,  54,
 34,  25,  23,  44,  11,  48,  62,   5,  51,   8,   6,  61,  26,  33,  47,  20,
  7,  60,  50,   9,  46,  21,  27,  32,  22,  45,  35,  24,  63,   4,  10,  49,
 37,  30,  16,  43,  12,  55,  57,   2,  52,  15,   1,  58,  29,  38,  40,  19,
  0,  60,  59,   7,  53,   9,  14,  50,  41,  21,  18,  46,  28,  32,  39,  27,
 17,  45,  42,  22,  36,  24,  31,  35,  56,   4,   3,  63,  13,  49,  54,  10,
 34,  30,  25,  37,  23,  43,  44,  16,  11,  55,  48,  12,  62,   2,   5,  57,
 51,  15,   8,  52,   6,  58,  61,   1,  26,  38,  33,  29,  47,  19,  20,  40,
  0,  30,  60,  34,  59,  37,   7,  25,  53,  43,   9,  23,  14,  16,  50,  44,
 41,  55,  21,  11,  18,  12,  46,  48,  28,   2,  32,  62,  39,  57,  27,   5,
 17,  15,  45,  51,  42,  52,  22,   8,  36,  58,  24,   6,  31,   1,  35,  61,
 56,  38,   4,  26,   3,  29,  63,  33,  13,  19,  49,  47,  54,  40,  10,  20,
  0,  15,  30,  17,  60,  51,  34,  45,  59,  52,  37,  42,   7,   8,  25,  22,
 53,  58,  43,  36,   9,   6,  23,  24,  14,   1,  16,  31,  50,  61,  44,  35,
 41,  38,  55,  56,  21,  26,  11,   4,  18,  29,  12,   3,  46,  33,  48,  63,
 28,  19,   2,  13,  32,  47,  62,  49,  39,  40,  57,  54,  27,  20,   5,  10,
  0,  38,  15,  41,  30,  56,  17,  55,  60,  26,  51,  21,  34,   4,  45,  11,
 59,  29,  52,  18,  37,   3,  42,  12,   7,  33,   8,  46,  25,  63,  22,  48,
 53,  19,  58,  28,  43,  13,  36,   2,   9,  47,   6,  32,  23,  49,  24,  62,
 14,  40,   1,  39,  16,  54,  31,  57,  50,  20,  61,  27,  44,  10,  35,   5,
  0,  19,  38,  53,  15,  28,  41,  58,  30,  13,  56,  43,  17,   2,  55,  36,
 60,  47,  26,   9,  51,  32,  21,   6,  34,  49,   4,  23,  45,  62,  11,  24,
 59,  40,  29,  14,  52,  39,  18,   1,  37,  54,   3,  16,  42,  57,  12,  31,
  7,  20,  33,  50,   8,  27,  46,  61,  25,  10,  63,  44,  22,   5,  48,  35,
  0,  40,  19,  59,  38,  14,  53,  29,  15,  39,  28,  52,  41,   1,  58,  18,
 30,  54,  13,  37,  56,  16,  43,   3,  17,  57,   2,  42,  55,  31,  36,  12,
 60,  20,  47,   7,  26,  50,   9,  33,  51,  27,  32,   8,  21,  61,   6,  46,
 34,  10,  49,  25,   4,  44,  23,  63,  45,   5,  62,  22,  11,  35,  24,  48,
  0,  20,  40,  60,  19,   7,  59,  47,  38,  50,  14,  26,  53,  33,  29,   9,
 15,  27,  39,  51,  28,   8,  52,  32,  41,  61,   1,  21,  58,  46,  18,   6,
 30,  10,  54,  34,  13,  25,  37,  49,  56,  44,  16,   4,  43,  63,   3,  23,
 17,   5,  57,  45,   2,  22,  42,  62,  55,  35,  31,  11,  36,  48,  12,  24,
  0,  10,  20,  30,  40,  34,  60,  54,  19,  25,   7,  13,  59,  49,  47,  37,
 38,  44,  50,  56,  14,   4,  26,  16,  53,  63,  33,  43,  29,  23,   9,   3,
 15,   5,  27,  17,  39,  45,  51,  57,  28,  22,   8,   2,  52,  62,  32,  42,
 41,  35,  61,  55,   1,  11,  21,  31,  58,  48,  46,  36,  18,  24,   6,  12,
  0,   5,  10,  15,  20,  17,  30,  27,  40,  45,  34,  39,  60,  57,  54,  51,
 19,  22,  25,  28,   7,   2,  13,   8,  59,  62,  49,  52,  47,  42,  37,  32,
 38,  35,  44,  41,  50,  55,  56,  61,  14,  11,   4,   1,  26,  31,  16,  21,
 53,  48,  63,  58,  33,  36,  43,  46,  29,  24,  23,  18,   9,  12,   3,   6,
  0,  35,   5,  38,  10,  41,  15,  44,  20,  55,  17,  50,  30,  61,  27,  56,
 40,  11,  45,  14,  34,   1,  39,   4,  60,  31,  57,  26,  54,  21,  51,  16,
 19,  48,  22,  53,  25,  58,  28,  63,   7,  36,   2,  33,  13,  46,   8,  43,
 59,  24,  62,  29,  49,  18,  52,  23,  47,  12,  42,   9,  37,   6,  32,   3,
  0,  48,  35,  19,   5,  53,  38,  22,  10,  58,  41,  25,  15,  63,  44,  28,
 20,  36,  55,   7,  17,  33,  50,   2,  30,  46,  61,  13,  27,  43,  56,   8,
 40,  24,  11,  59,  45,  29,  14,  62,  34,  18,   1,  49,  39,  23,   4,  52,
 60,  12,  31,  47,  57,   9,  26,  42,  54,   6,  21,  37,  51,   3,  16,  32,
  0,  24,  48,  40,  35,  59,  19,  11,   5,  29,  53,  45,  38,  62,  22,  14,
 10,  18,  58,  34,  41,  49,  25,   1,  15,  23,  63,  39,  44,  52,  28,   4,
 20,  12,  36,  60,  55,  47,   7,  31,  17,   9,  33,  57,  50,  42,   2,  26,
 30,   6,  46,  54,  61,  37,  13,  21,  27,   3,  43,  51,  56,  32,   8,  16,
  0,  12,  24,  20,  48,  60,  40,  36,  35,  47,  59,  55,  19,  31,  11,   7,
  5,   9,  29,  17,  53,  57,  45,  33,  38,  42,  62,  50,  22,  26,  14,   2,
 10,   6,  18,  30,  58,  54,  34,  46,  41,  37,  49,  61,  25,  21,   1,  13,
 15,   3,  23,  27,  63,  51,  39,  43,  44,  32,  52,  56,  28,  16,   4,   8,
  0,   6,  12,  10,  24,  30,  20,  18,  48,  54,  60,  58,  40,  46,  36,  34,
 35,  37,  47,  41,  59,  61,  55,  49,  19,  21,  31,  25,  11,  13,   7,   1,
  5,   3,   9,  15,  29,  27,  17,  23,  53,  51,  57,  63,  45,  43,  33,  39,
 38,  32,  42,  44,  62,  56,  50,  52,  22,  16,  26,  28,  14,   8,   2,   4,
  0,   3,   6,   5,  12,  15,  10,   9,  24,  27,  30,  29,  20,  23,  18,  17,
 48,  51,  54,  53,  60,  63,  58,  57,  40,  43,  46,  45,  36,  39,  34,  33,
 35,  32,  37,  38,  47,  44,  41,  42,  59,  56,  61,  62,  55,  52,  49,  50,
 19,  16,  21,  22,  31,  28,  25,  26,  11,   8,  13,  14,   7,   4,   1,   2,
  0,  32,   3,  35,   6,  38,   5,  37,  12,  44,  15,  47,  10,  42,   9,  41,
 24,  56,  27,  59,  30,  62,  29,  61,  20,  52,  23,  55,  18,  50,  17,  49,
 48,  16,  51,  19,  54,  22,  53,  21,  60,  28,  63,  31,  58,  26,  57,  25,
 40,   8,  43,  11,  46,  14,  45,  13,  36,   4,  39,   7,  34,   2,  33,   1,
  0,  16,  32,  48,   3,  19,  35,  51,   6,  22,  38,  54,   5,  21,  37,  53,
 12,  28,  44,  60,  15,  31,  47,  63,  10,  26,  42,  58,   9,  25,  41,  57,
 24,   8,  56,  40,  27,  11,  59,  43,  30,  14,  62,  46,  29,  13,  61,  45,
 20,   4,  52,  36,  23,   7,  55,  39,  18,   2,  50,  34,  17,   1,  49,  33,
  0,   8,  16,  24,  32,  40,  48,  56,   3,  11,  19,  27,  35,  43,  51,  59,
  6,  14,  22,  30,  38,  46,  54,  62,   5,  13,  21,  29,  37,  45,  53,  61,
 12,   4,  28,  20,  44,  36,  60,  52,  15,   7,  31,  23,  47,  39,  63,  55,
 10,   2,  26,  18,  42,  34,  58,  50,   9,   1,  25,  17,  41,  33,  57,  49,
  0,   4,   8,  12,  16,  20,  24,  28,  32,  36,  40,  44,  48,  52,  56,  60,
  3,   7,  11,  15,  19,  23,  27,  31,  35,  39,  43,  47,  51,  55,  59,  63,
  6,   2,  14,  10,  22,  18,  30,  26,  38,  34,  46,  42,  54,  50,  62,  58,
  5,   1,  13,   9,  21,  17,  29,  25,  37,  33,  45,  41,  53,  49,  61,  57,
  0,   2,   4,   6,   8,  10,  12,  14,  16,  18,  20,  22,  24,  26,  28,  30,
 32,  34,  36,  38,  40,  42,  44,  46,  48,  50,  52,  54,  56,  58,  60,  62,
  3,   1,   7,   5,  11,   9,  15,  13,  19,  17,  23,  21,  27,  25,  31,  29,
 35,  33,  39,  37,  43,  41,  47,  45,  51,  49,  55,  53,  59,  57,  63,  61
};

const qracode qra_13_64_64_irr_e = {
		qra_K,
		qra_N,
		qra_m,
		qra_M,
		qra_a,
		qra_NC,
		qra_V,
		qra_C,
		qra_NMSG,
		qra_MAXVDEG,
		qra_MAXCDEG,
		QRATYPE_CRCPUNCTURED,
		qra_R,
		CODE_NAME,
		qra_acc_input_idx,
		qra_acc_input_wlog,
		qra_log,
		qra_exp,
		qra_msgw,
		qra_vdeg,
		qra_cdeg,
		qra_v2cmidx,
		qra_c2vmidx,
		qra_pmat
};

#undef qra_K
#undef qra_N
#undef qra_m
#undef qra_M
#undef qra_a
#undef qra_NC
#undef qra_V
#undef qra_C
#undef qra_NMSG
#undef qra_MAXVDEG
#undef qra_MAXCDEG
#undef qra_R
#undef CODE_NAME


code:
// qra13_64_64_irr_e.h
// Code tables and defines for Q-ary RA code (13,64) over GF(64)
// Code Name: qra13_64_64_irr_e
// (13,64) RA Code over GF(64) RF=[3x4 4x4 6x1 3x2 5x1 7x1]/18

// (c) 2016 - Nico Palermo - IV3NWV - Microtelecom Srl, Italy

// 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.
//
//    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/>.

#ifndef _qra13_64_64_irr_e_h
#define _qra13_64_64_irr_e_h

#include "qracodes.h"

#ifdef __cplusplus
extern "C" {
#endif

extern const qracode qra_13_64_64_irr_e;

#ifdef __cplusplus
}
#endif

#endif // _qra13_64_64_irr_e_h

license: mit


code:
// qracodes.c
// Q-ary RA codes encoding/decoding functions
//
// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy
// ------------------------------------------------------------------------------
// 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.
//
//    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/>.

#include <stdio.h>
#include <math.h>

#include "npfwht.h"
#include "pdmath.h"

#include "qracodes.h"

int qra_encode(const qracode *pcode, int *y, const int *x)
{
	int k,j,kk,jj;
	int t, chk = 0;

	const int K = pcode->K;
	const int M = pcode->M;
	const int NC= pcode->NC;
	const int a = pcode->a;
	const int  *acc_input_idx  = pcode->acc_input_idx;
	const int *acc_input_wlog = pcode->acc_input_wlog;
	const int  *gflog		   = pcode->gflog;
	const int *gfexp          = pcode->gfexp;

	// copy the systematic symbols to destination
	memcpy(y,x,K*sizeof(int));

	y = y+K;	// point to check symbols

	// compute the code check symbols as a weighted accumulation of a permutated
	// sequence of the (repeated) systematic input symbols:
	// chk(k+1) = x(idx(k))*alfa^(logw(k)) + chk(k)
	// (all operations performed over GF(M))

	if (a==1) { // grouping factor = 1
		for (k=0;k<NC;k++) {
			t = x[acc_input_idx[k]];
			if (t) {
				// multiply input by weight[k] and xor it with previous check
				t = (gflog[t] + acc_input_wlog[k])%(M-1);
				t = gfexp[t];
				chk ^=t;
				}
			y[k] = chk;
			}

		#ifdef QRA_DEBUG
			// verify that the encoder accumulator is terminated to 0
			// (we designed the code this way so that Iex = 1 when Ia = 1)
			t = x[acc_input_idx[k]];
			if (t) {
				t = (gflog[t] + acc_input_wlog[k])%(M-1);
				t = gfexp[t];
				// accumulation
				chk ^=t;
				}
			return (chk==0);
		#else
			return 1;
		#endif // QRA_DEBUG
		}
	else { // grouping factor > 1
		for (k=0;k<NC;k++) {
			kk = a*k;
			for (j=0;j<a;j++) {
				jj = kk+j;
				// irregular grouping support
				if (acc_input_idx[jj]<0)
					continue;
				t = x[acc_input_idx[jj]];
				if (t) {
					// multiply input by weight[k] and xor it with previous check
					t = (gflog[t] + acc_input_wlog[jj])%(M-1);
					t = gfexp[t];
					chk ^=t;
					}
				}
			y[k] = chk;
			}
		#ifdef QRA_DEBUG
			// verify that the encoder accumulator is terminated to 0
			// (we designed the code this way so that Iex = 1 when Ia = 1)
			kk = a*k;
			for (j=0;j<a;j++) {
				jj = kk+j;
				if (acc_input_idx[jj]<0)
					continue;
				t = x[acc_input_idx[jj]];
				if (t) {
					// multiply input by weight[k] and xor it with previous check
					t = (gflog[t] + acc_input_wlog[jj])%(M-1);
					t = gfexp[t];
					chk ^=t;
					}
				}
			return (chk==0);
		#else
			return 1;
		#endif // QRA_DEBUG
		}
}

static void qra_ioapprox(float *src, float C, int nitems)
{
	// In place approximation of the modified bessel function I0(x*C)
	// Computes src[k] = Io(src[k]*C) where Io() is the modified Bessel function of first kind and order 0

	float v;
	float vsq;

	while (nitems--) {
		v = src[nitems]*C;

		// rational approximation of log(Io(v))
		vsq = v*v;
		v = vsq*(v+0.039f)/(vsq*.9931f+v*2.6936f+0.5185f);

		if (v>80.f) // avoid floating point exp() overflows
			v=80.f;

		src[nitems] = (float)exp(v);
		}
}

float qra_mfskbesselmetric(float *pix, const float *rsq, const int m, const int N, float EsNoMetric)
{
	// Computes the codeword symbols intrinsic probabilities
	// given the square of the received input amplitudes.

	// The input vector rqs must be a linear array of size M*N, where M=2^m,
	// containing the squared amplitudes (rp*rp+rq*rq) of the input samples

	// First symbol amplitudes should be stored in the first M positions,
	// second symbol amplitudes stored at positions [M ... 2*M-1], and so on.

	// Output vector is the intrinsic symbol metric (the probability distribution)
	// assuming that symbols are transmitted using a M-FSK modulation
	// and incoherent demodulation.

	// As the input Es/No is generally unknown (as it cannot be exstimated accurately
	// when the codeword length is few tens symbols) but an exact metric requires it
	// we simply fix it to a predefined EsNoMetric value so that the metric is what
	// expected at that specific value.
	// The metric computed in this way is optimal only at this predefined Es/No value,
	// nevertheless it is usually better than a generic parameter-free metric which
	// makes no assumptions on the input Es/No.

	// returns the estimated noise standard deviation

	int k;
	float rsum = 0.f;
	float sigmaest, cmetric;

	const int M = 1<<m;
	const int nsamples = M*N;

	// compute total power and modulus of input signal
	for (k=0;k<nsamples;k++) {
		rsum = rsum+rsq[k];
		pix[k] = (float)sqrt(rsq[k]);
		}

	rsum = rsum/nsamples;	// average S+N

	// IMPORTANT NOTE: in computing the noise stdev it is assumed that
	// in the input amplitudes there's no strong interference!
	// A more robust estimation can be done evaluating the histogram of the input amplitudes

	sigmaest = (float)sqrt(rsum/(1.0f+EsNoMetric/M)/2); // estimated noise stdev
	cmetric = (float)sqrt(2*EsNoMetric)/sigmaest;

	for (k=0;k<N;k++) {
		// compute bessel metric for each symbol in the codeword
		qra_ioapprox(PD_ROWADDR(pix,M,k),cmetric,M);
		// normalize to a probability distribution
		pd_norm(PD_ROWADDR(pix,M,k),m);
		}

	return sigmaest;
}

#ifdef QRA_DEBUG
void pd_print(int imsg,float *ppd,int size)
{
	int k;
	printf("imsg=%d\n",imsg);
	for (k=0;k<size;k++)
		printf("%7.1e ",ppd[k]);
	printf("\n");
}
#endif

#define ADDRMSG(fp, msgidx)    PD_ROWADDR(fp,qra_M,msgidx)
#define C2VMSG(msgidx)         PD_ROWADDR(qra_c2vmsg,qra_M,msgidx)
#define V2CMSG(msgidx)         PD_ROWADDR(qra_v2cmsg,qra_M,msgidx)
#define MSGPERM(logw)          PD_ROWADDR(qra_pmat,qra_M,logw)

#define QRACODE_MAX_M	256	// Maximum alphabet size handled by qra_extrinsic

int qra_extrinsic(const qracode *pcode,
				  float *pex,
				  const float *pix,
				  int maxiter,
				  float *qra_v2cmsg,
				  float *qra_c2vmsg)
{
	const int qra_M		= pcode->M;
	const int qra_m		= pcode->m;
	const int qra_V		= pcode->V;
	const int qra_MAXVDEG  = pcode->MAXVDEG;
	const int  *qra_vdeg    = pcode->vdeg;
	const int qra_C		= pcode->C;
	const int qra_MAXCDEG  = pcode->MAXCDEG;
	const int *qra_cdeg    = pcode->cdeg;
	const int  *qra_v2cmidx = pcode->v2cmidx;
	const int  *qra_c2vmidx = pcode->c2vmidx;
	const int  *qra_pmat    = pcode->gfpmat;
	const int *qra_msgw    = pcode->msgw;

//	float msgout[qra_M];		 // buffer to store temporary results
	float msgout[QRACODE_MAX_M]; // we use a fixed size in order to avoid mallocs

	float totex;	// total extrinsic information
	int nit;		// current iteration
	int nv;		// current variable
	int nc;		// current check
	int k,kk;		// loop indexes

	int ndeg;		// current node degree
	int msgbase;	// current offset in the table of msg indexes
	int imsg;		// current message index
	int wmsg;		// current message weight

	int rc     = -1; // rc>=0  extrinsic converged to 1 at iteration rc (rc=0..maxiter-1)
					 // rc=-1  no convergence in the given number of iterations
					 // rc=-2  error in the code tables (code checks degrees must be >1)
					 // rc=-3  M is larger than QRACODE_MAX_M

	if (qra_M>QRACODE_MAX_M)
		return -3;

	// message initialization -------------------------------------------------------

	// init c->v variable intrinsic msgs
	pd_init(C2VMSG(0),pix,qra_M*qra_V);

	// init the v->c messages directed to code factors (k=1..ndeg) with the intrinsic info
	for (nv=0;nv<qra_V;nv++) {

		ndeg = qra_vdeg[nv];		// degree of current node
		msgbase = nv*qra_MAXVDEG;	// base to msg index row for the current node

		// copy intrinsics on v->c
		for (k=1;k<ndeg;k++) {
			imsg = qra_v2cmidx[msgbase+k];
			pd_init(V2CMSG(imsg),ADDRMSG(pix,nv),qra_M);
			}
		}

	// message passing algorithm iterations ------------------------------

	for (nit=0;nit<maxiter;nit++) {

		// c->v step -----------------------------------------------------
		// Computes messages from code checks to code variables.
		// As the first qra_V checks are associated with intrinsic information
		// (the code tables have been constructed in this way)
		// we need to do this step only for code checks in the range [qra_V..qra_C)

		// The convolutions of probability distributions over the alphabet of a finite field GF(qra_M)
		// are performed with a fast convolution algorithm over the given field.
		//
		// I.e. given the code check x1+x2+x3 = 0 (with x1,x2,x3 in GF(2^m))
		// and given Prob(x2) and Prob(x3), we have that:
		// Prob(x1=X1) = Prob((x2+x3)=X1) = sum((Prob(x2=X2)*Prob(x3=(X1+X2))) for all the X2s in the field
		// This translates to Prob(x1) = IWHT(WHT(Prob(x2))*WHT(Prob(x3)))
		// where WHT and IWHT are the direct and inverse Walsh-Hadamard transforms of the argument.
		// Note that the WHT and the IWHF differs only by a multiplicative coefficent and since in this step
		// we don't need that the output distribution is normalized we use the relationship
		// Prob(x1) =(proportional to) WH(WH(Prob(x2))*WH(Prob(x3)))

		// In general given the check code x1+x2+x3+..+xm = 0
		// the output distribution of a variable given the distributions of the other m-1 variables
		// is the inverse WHT of the product of the WHTs of the distribution of the other m-1 variables
		// The complexity of this algorithm scales with M*log2(M) instead of the M^2 complexity of
		// the brute force approach (M=size of the alphabet)

		for (nc=qra_V;nc<qra_C;nc++) {

			ndeg = qra_cdeg[nc];		// degree of current node

			if (ndeg==1) 				// this should never happen (code factors must have deg>1)
				return -2;				// bad code tables

			msgbase = nc*qra_MAXCDEG;	// base to msg index row for the current node

			// transforms inputs in the Walsh-Hadamard "frequency" domain
			// v->c  -> fwht(v->c)
			for (k=0;k<ndeg;k++) {
				imsg = qra_c2vmidx[msgbase+k];		// msg index
				np_fwht(qra_m,V2CMSG(imsg),V2CMSG(imsg)); // compute fwht
				}

			// compute products and transform them back in the WH "time" domain
			for (k=0;k<ndeg;k++) {

				// init output message to uniform distribution
				pd_init(msgout,pd_uniform(qra_m),qra_M);

				// c->v = prod(fwht(v->c))
				// TODO: we assume that checks degrees are not larger than three but
				// if they are larger the products can be computed more efficiently
				for (kk=0;kk<ndeg;kk++)
					if (kk!=k) {
						imsg = qra_c2vmidx[msgbase+kk];
						pd_imul(msgout,V2CMSG(imsg),qra_m);
						}

				// transform product back in the WH "time" domain

				// Very important trick:
				// we bias WHT[0] so that the sum of output pd components is always strictly positive
				// this helps avoiding the effects of underflows in the v->c steps when multipling
				// small fp numbers
				msgout[0]+=1E-7f;	// TODO: define the bias accordingly to the field size

				np_fwht(qra_m,msgout,msgout);

				// inverse weight and output
				imsg = qra_c2vmidx[msgbase+k]; // current output msg index
				wmsg = qra_msgw[imsg];		   // current msg weight

				if (wmsg==0)
					pd_init(C2VMSG(imsg),msgout,qra_M);
				else
					// output p(alfa^(-w)*x)
					pd_bwdperm(C2VMSG(imsg),msgout, MSGPERM(wmsg), qra_M);

				} // for (k=0;k<ndeg;k++)

			} // for (nc=qra_V;nc<qra_C;nc++)

		// v->c step -----------------------------------------------------
		for (nv=0;nv<qra_V;nv++) {

			ndeg = qra_vdeg[nv];		// degree of current node
			msgbase = nv*qra_MAXVDEG;	// base to msg index row for the current node

			for (k=0;k<ndeg;k++) {
				// init output message to uniform distribution
				pd_init(msgout,pd_uniform(qra_m),qra_M);

				// v->c msg = prod(c->v)
				// TODO: factor factors to reduce the number of computations for high degree nodes
				for (kk=0;kk<ndeg;kk++)
					if (kk!=k) {
						imsg = qra_v2cmidx[msgbase+kk];
						pd_imul(msgout,C2VMSG(imsg),qra_m);
						}

#ifdef QRA_DEBUG
// normalize and check if product of messages v->c are null
				// normalize output to a probability distribution
				if (pd_norm(msgout,qra_m)<=0) {
					// dump msgin;
					printf("warning: v->c pd with invalid norm. nit=%d nv=%d k=%d\n",nit,nv,k);
					for (kk=0;kk<ndeg;kk++) {
						imsg = qra_v2cmidx[msgbase+kk];
						pd_print(imsg,C2VMSG(imsg),qra_M);
						}
					printf("-----------------\n");
					}
#else
				// normalize the result to a probability distribution
				pd_norm(msgout,qra_m);
#endif
				// weight and output
				imsg = qra_v2cmidx[msgbase+k]; // current output msg index
				wmsg = qra_msgw[imsg];		   // current msg weight

				if (wmsg==0) {
					pd_init(V2CMSG(imsg),msgout,qra_M);
					}
				else {
					// output p(alfa^w*x)
					pd_fwdperm(V2CMSG(imsg),msgout, MSGPERM(wmsg), qra_M);
					}

				} // for (k=0;k<ndeg;k++)
			} // for (nv=0;nv<qra_V;nv++)

		// check extrinsic information ------------------------------
		// We assume that decoding is successful if each of the extrinsic
		// symbol probability is close to ej, where ej = [0 0 0 1(j-th position) 0 0 0 ]
		// Therefore, for each symbol k in the codeword we compute max(prob(Xk))
		// and we stop the iterations if sum(max(prob(xk)) is close to the codeword length
		// Note: this is a more restrictive criterium than that of computing the a
		// posteriori probability of each symbol, making a hard decision and then check
		// if the codeword syndrome is null.
		// WARNING: this is tricky and probably works only for the particular class of RA codes
		// we are coping with (we designed the code weights so that for any input symbol the
		// sum of its weigths is always 0, thus terminating the accumulator trellis to zero
		// for every combination of the systematic symbols).
		// More generally we should instead compute the max a posteriori probabilities
		// (as a product of the intrinsic and extrinsic information), make a symbol by symbol hard
		// decision and then check that the syndrome of the result is indeed null.

		totex = 0;
		for (nv=0;nv<qra_V;nv++)
			totex += pd_max(V2CMSG(nv),qra_M);

		if (totex>(1.*(qra_V)-0.01)) {
			// the total maximum extrinsic information of each symbol in the codeword
			// is very close to one. This means that we have reached the (1,1) point in the
			// code EXIT chart(s) and we have successfully decoded the input.
			rc = nit;
			break;	// remove the break to evaluate the decoder speed performance as a function of the max iterations number)
			}

		} // for (nit=0;nit<maxiter;nit++)

	// copy extrinsic information to output to do the actual max a posteriori prob decoding
	pd_init(pex,V2CMSG(0),(qra_M*qra_V));
	return rc;
}

void qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix)
{
// Maximum a posteriori probability decoding.
// Given the intrinsic information (pix) and extrinsic information (pex) (computed with qra_extrinsic(...))
// compute pmap = pex*pix and decode each (information) symbol of the received codeword
// as the symbol which maximizes pmap

// Returns:
//	xdec[k] = decoded (information) symbols k=[0..qra_K-1]

//  Note: pex is destroyed and overwritten with mapp

	const int qra_M		= pcode->M;
	const int qra_m		= pcode->m;
	const int qra_K		= pcode->K;

	int k;

	for (k=0;k<qra_K;k++) {
		// compute a posteriori prob
		pd_imul(PD_ROWADDR(pex,qra_M,k),PD_ROWADDR(pix,qra_M,k),qra_m);
		xdec[k]=pd_argmax(NULL, PD_ROWADDR(pex,qra_M,k), qra_M);
		}
}


code:
// qracodes.h
// Q-ary RA codes encoding/decoding functions
//
// (c) 2016 - Nico Palermo, IV3NWV - Microtelecom Srl, Italy
// ------------------------------------------------------------------------------
// 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.
//
//    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/>.

#ifndef _qracodes_h_
#define _qracodes_h_

// type of codes
#define QRATYPE_NORMAL			0x00 // normal code
#define QRATYPE_CRC 			0x01 // code with crc - last information symbol is a CRC
#define QRATYPE_CRCPUNCTURED 	0x02 // the CRC symbol is punctured (not sent along the channel)

typedef struct {
	// code parameters
	const int K;			// number of information symbols
	const int N;			// codeword length in symbols
	const int m;			// bits/symbol
	const int M;			// Symbol alphabet cardinality (2^m)
	const int a;			// code grouping factor
	const int NC;			// number of check symbols (N-K)
	const int V;			// number of variables in the code graph (N)
	const int C;			// number of factors in the code graph (N +(N-K)+1)
	const int NMSG;		// number of msgs in the code graph
	const int MAXVDEG;		// maximum variable degree
	const int MAXCDEG;		// maximum factor degree
	const int type;		// see QRATYPE_xx defines
	const float R;			// code rate (K/N)
	const char  name[64];	// code name
	// tables used by the encoder
	const int	 *acc_input_idx;
	const int   *acc_input_wlog;
	const int	 *gflog;
	const int   *gfexp;
	// tables used by the decoder -------------------------
	const int *msgw;
	const int *vdeg;
	const int *cdeg;
	const int  *v2cmidx;
	const int  *c2vmidx;
	const int  *gfpmat;
} qracode;
// Uncomment the header file of the code which needs to be tested

//#include "qra12_63_64_irr_b.h"  // irregular code (12,63) over GF(64)
//#include "qra13_64_64_irr_e.h"  // irregular code with good performance and best UER protection at AP56
//#include "qra13_64_64_reg_a.h"  // regular code with good UER but perf. inferior to that of code qra12_63_64_irr_b

#ifdef __cplusplus
extern "C" {
#endif

int   qra_encode(const qracode *pcode, int *y, const int *x);
float qra_mfskbesselmetric(float *pix, const float *rsq, const int m, const int N, float EsNoMetric);
int   qra_extrinsic(const qracode *pcode, float *pex, const float *pix, int maxiter,float *qra_v2cmsg,float *qra_c2vmsg);
void  qra_mapdecode(const qracode *pcode, int *xdec, float *pex, const float *pix);

#ifdef __cplusplus
}
#endif

#endif // _qracodes_h_


code:
subroutine qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
     emedelay,mycall_12,hiscall_12,hisgrid_6,sync,nsnr,dtx,nfreq,decoded,nft)

  use packjt
  use timer_module, only: timer

  parameter (NMAX=60*12000,LN=1152*63)
  character decoded*22
  character*12 mycall_12,hiscall_12
  character*6 mycall,hiscall,hisgrid_6
  character*4 hisgrid
  logical ltext
  complex c00(0:720000)                      !Complex spectrum of dd()
  complex c0(0:720000)                       !Complex data for dd()
  real a(3)
  real dd(NMAX)                              !Raw data sampled at 12000 Hz
  real s3(LN)                                !Symbol spectra
  real s3a(LN)                               !Symbol spectra
  integer dat4(12)                           !Decoded message (as 12 integers)
  integer dat4x(12)
  integer nap(0:11)
  data nap/0,2,3,2,3,4,2,3,6,4,6,6/
  data nc1z/-1/,nc2z/-1/,ng2z/-1/,maxaptypez/-1/
  save

  call timer('qra64a  ',0)
  irc=-1
  decoded='                      '
  nft=99
  if(nfqso.lt.nf1 .or. nfqso.gt.nf2) go to 900

  mycall=mycall_12(1:6)                     !### May need fixing ###
  hiscall=hiscall_12(1:6)
  hisgrid=hisgrid_6(1:4)
  call packcall(mycall,nc1,ltext)
  call packcall(hiscall,nc2,ltext)
  call packgrid(hisgrid,ng2,ltext)
  nSubmode=0
  if(mode64.eq.2) nSubmode=1
  if(mode64.eq.4) nSubmode=2
  if(mode64.eq.8) nSubmode=3
  if(mode64.eq.16) nSubmode=4
  b90=1.0
  nFadingModel=1
  maxaptype=4
  if(iand(ndepth,64).ne.0) maxaptype=5
  if(nc1.ne.nc1z .or. nc2.ne.nc2z .or. ng2.ne.ng2z .or.            &
     maxaptype.ne.maxaptypez) then
     do naptype=0,maxaptype
        if(naptype.eq.2 .and. maxaptype.eq.4) cycle
        call qra64_dec(s3,nc1,nc2,ng2,naptype,1,nSubmode,b90,      &
             nFadingModel,dat4,snr2,irc)
     enddo
     nc1z=nc1
     nc2z=nc2
     ng2z=ng2
     maxaptypez=maxaptype
  endif
  naptype=maxaptype

  call ana64(dd,npts,c00)
  npts2=npts/2

  call timer('sync64  ',0)
  call sync64(c00,nf1,nf2,nfqso,ntol,mode64,emedelay,dtx,f0,jpk0,sync,  &
       sync2,width)
  call timer('sync64  ',1)
  nfreq=nint(f0)
  if(mode64.eq.1 .and. minsync.ge.0 .and. (sync-7.0).lt.minsync) go to 900
!  if((sync-3.4).lt.float(minsync) .or.width.gt.340.0) go to 900
  a=0.
  a(1)=-f0
  call twkfreq(c00,c0,npts2,6000.0,a)

  irc=-99
  s3lim=20.
  itz=11
  if(mode64.eq.4) itz=9
  if(mode64.eq.2) itz=7
  if(mode64.eq.1) itz=5

  LL=64*(mode64+2)
  NN=63
  napmin=99
  do itry0=1,5
     idt=itry0/2
     if(mod(itry0,2).eq.0) idt=-idt
     jpk=jpk0 + 750*idt
     call spec64(c0,npts2,mode64,jpk,s3a,LL,NN)
     call pctile(s3a,LL*NN,40,base)
     s3a=s3a/base
     where(s3a(1:LL*NN)>s3lim) s3a(1:LL*NN)=s3lim
     do iter=itz,0,-2
        b90=1.728**iter
        if(b90.gt.230.0) cycle
        if(b90.lt.0.15*width) exit
        s3(1:LL*NN)=s3a(1:LL*NN)
        call timer('qra64_de',0)
        call qra64_dec(s3,nc1,nc2,ng2,naptype,0,nSubmode,b90,      &
             nFadingModel,dat4,snr2,irc)
        call timer('qra64_de',1)
        if(irc.eq.0) go to 10
        if(irc.gt.0) call badmsg(irc,dat4,nc1,nc2,ng2)
        iirc=max(0,min(irc,11))
        if(irc.gt.0 .and. nap(iirc).lt.napmin) then
           dat4x=dat4
           b90x=b90
           snr2x=snr2
           napmin=nap(iirc)
           irckeep=irc
           dtxkeep=jpk/6000.0 - 1.0
           itry0keep=itry0
           iterkeep=iter
        endif
     enddo
     if(irc.eq.0) exit
  enddo

  if(napmin.ne.99) then
     dat4=dat4x
     b90=b90x
     snr2=snr2x
     irc=irckeep
     dtx=dtxkeep
     itry0=itry0keep
     iter=iterkeep
  endif
10 decoded='                      '

  if(irc.ge.0) then
     call unpackmsg(dat4,decoded,.false.,'      ') !Unpack the user message
     call fmtmsg(decoded,iz)
     if(index(decoded,"000AAA ").ge.1) then
        ! Suppress a certain type of garbage decode.
        decoded='                      '
        irc=-1
     endif
     nft=100 + irc
     nsnr=nint(snr2)
  else
     snr2=0.
  endif

900 if(irc.lt.0) then
     sy=max(1.0,sync)
     if(nSubmode.eq.0) nsnr=nint(10.0*log10(sy)-35.0)   !A
     if(nSubmode.eq.1) nsnr=nint(10.0*log10(sy)-34.0)   !B
     if(nSubmode.eq.2) nsnr=nint(10.0*log10(sy)-29.0)   !C
     if(nSubmode.eq.3) nsnr=nint(10.0*log10(sy)-29.0)   !D
     if(nSubmode.eq.4) nsnr=nint(10.0*log10(sy)-24.0)   !E
  endif
  call timer('qra64a  ',1)

!  write(71,3001) nutc,dtx,f0,sync,sync2,width,minsync,decoded
!3001 format(i4.4,f7.2,4f8.1,i3,2x,a22)

  return
end subroutine qra64a


code:
program QRA64code

! Provides examples of message packing, bit and symbol ordering,
! QRA (63,12) encoding, and other necessary details of the QRA64
! protocol.

  use packjt
  character*22 msg,msg0,msg1,decoded,cok*3,msgtype*10,arg*12
  character*6 mycall
  logical ltext
  integer dgen(12),sent(63),dec(12)
  integer icos7(0:6)
  data icos7/4,2,5,6,1,3,0/                   !Costas 7x7 tone pattern

  include 'testmsg.f90'

  nargs=iargc()
  if(nargs.lt.1) then
     print*,'Usage: qra64code "message"'
     print*,'       qra64code -t'
     go to 999
  endif

  call getarg(1,msg)                     !Get message from command line
  nmsg=1
  if(msg(1:2).eq."-t") nmsg=NTEST

  write(*,1010)
1010 format("     Message                 Decoded                Err? Type"/74("-"))

  do imsg=1,nmsg
     if(nmsg.gt.1) msg=testmsg(imsg)
    call fmtmsg(msg,iz)                    !To upper, collapse mult blanks
     msg0=msg                               !Input message
     call chkmsg(msg,cok,nspecial,flip)     !See if it includes "OOO" report
     msg1=msg                               !Message without "OOO"
     call packmsg(msg1,dgen,itype,.false.) !Pack message into 12 six-bit bytes
     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 pfx"
     if(itype.eq.3) msgtype="Type 1 sfx"
     if(itype.eq.4) msgtype="Type 2 pfx"
     if(itype.eq.5) msgtype="Type 2 sfx"
     if(itype.eq.6) msgtype="Free text"

     call qra64_enc(dgen,sent)              !Encode using QRA64

     call unpackmsg(dgen,decoded,.false.,'      ') !Unpack the user message
     call fmtmsg(decoded,iz)
     ii=imsg
     write(*,1020) ii,msg0,decoded,itype,msgtype
1020 format(i4,1x,a22,2x,a22,4x,i3,": ",a13)
  enddo

  if(nmsg.eq.1) then
     write(*,1030) dgen
1030 format(/'Packed message, 6-bit symbols ',12i3) !Display packed symbols

     write(*,1040) sent
1040 format(/'Information-carrying channel symbols'/(i5,29i3))

     write(*,1050) 10*icos7,sent(1:32),10*icos7,sent(33:63),10*icos7
1050 format(/'Channel symbols including sync'/(i5,29i3))
  endif

999 end program QRA64code


code:
program qratest

  parameter (NMAX=60*12000)
  real dd(NMAX)
  character arg*8,mycall*12,hiscall*12,hisgrid*6,decoded*22
  character c*1
  logical loop

  nargs=iargc()
  if(nargs.lt.1 .or. nargs.gt.4) then
     print*,'Usage: qratest nfile [sync f0 fTol]'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) nfile
  loop=arg(1:1).eq.'+'
  minsync0=-1
  nfqso0=-1
  ntol0=-1
  if(nargs.gt.1) then
     call getarg(2,arg)
     read(arg,*) minsync0
     call getarg(3,arg)
     read(arg,*) nfqso0
     call getarg(4,arg)
     read(arg,*) ntol0
  endif
  ndepth=3
  nft=99

  open(60,file='qra64_data.bin',access='stream')

  do ifile=1,999
     read(60,end=999) dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth, &
          mycall,hiscall,hisgrid
     if(ifile.lt.nfile) cycle

     if(minsync0.ne.-1) minsync=minsync0
     if(nfqso0.ne.-1) nfqso=nfqso0
     if(ntol0.ne.-1) ntol=ntol0

     call qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,      &
          mycall,hiscall,hisgrid,sync,nsnr,dtx,nfreq,decoded,nft)
     c='a'
     if(mode64.eq.2) c='b'
     if(mode64.eq.4) c='c'
     if(mode64.eq.8) c='d'
     if(mode64.eq.16) c='e'
     write(*,1000) ifile,c,nutc,nsnr,dtx,nfreq,decoded,nft-100,sync-3.4
1000 format(i4,1x,a1,1x,i4.4,i4,f6.2,i5,1x,a22,i3,f6.2)
     if(ifile.eq.nfile .and. (.not.loop)) exit
  enddo

999 end program qratest


code:
gfortran -o twq -Wall -Wno-conversion -fbounds-check twq.f90 \
	 ../packjt.f90 wqenc.f90 wqdec.f90 packprop.f90 \
	 packname.f90 packtext2.f90 unpackprop.f90 unpackname.f90 \
	 unpacktext2.f90 unpackpfx.f90 pack50.f90 unpack50.f90 \
	 ../hash.f90 ../deg2grid.f90 ../grid2deg.f90 \
	 ../fix_contest_msg.f90 ../to_contest_msg.f90 \
	 ../fmtmsg.f90 ../azdist.f90 ../geodist.f90 ../wsprd/nhash.c


code:
subroutine pack50(n1,n2,dat)

  integer*1 dat(11),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


code:
subroutine packname(name,len,n1,n2)

  character*9 name
  real*8 dn

  dn=0
  do i=1,len
     n=ichar(name(i:i))
     if(n.ge.97 .and. n.le.122) n=n-32
     dn=27*dn + n-64
  enddo
  if(len.lt.9) then
     do i=len+1,9
        dn=27*dn
     enddo
  endif

  n2=mod(dn,32768.d0)
  dn=dn/32768.d0
  n1=dn

  return
end subroutine packname


code:
subroutine packprop(k,muf,ccur,cxp,n1)

! Pack propagation indicators into a 21-bit number.

! k      k-index, 0-9; 10="N/A"
! muf    muf, 2-60 MHz; 0=N/A, 1="none", 61=">60 MHz"
! ccur   up to two current events, each indicated by single
!        or double letter.
! cxp    zero or one expected event, indicated by single or
!        double letter

  character ccur*4,cxp*2

  j=ichar(ccur(1:1))-64
  if(j.lt.0) j=0
  n1=j
  do i=2,4
     if(ccur(i:i).eq.' ') go to 10
     if(ccur(i:i).eq.ccur(i-1:i-1)) then
        n1=n1+26
     else
        j=ichar(ccur(i:i))-64
        if(j.lt.0) j=0
        n1=53*n1 + j
     endif
  enddo

10 j=ichar(cxp(1:1))-64
  if(j.lt.0) j=0
  if(cxp(2:2).eq.cxp(1:1)) j=j+26
  n1=53*n1 + j
  n1=11*n1 + k
  n1=62*n1 + muf

  return
end subroutine packprop


code:
subroutine packtext2(msg,n1,ng)

  character*8 msg
  real*8 dn
  character*41 c
  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +./?'/

  dn=0.
  do i=1,8
     do j=1,41
        if(msg(i:i).eq.c(j:j)) go to 10
     enddo
     j=37
10   j=j-1                                !Codes should start at zero
     dn=41.d0*dn + j
  enddo

  ng=mod(dn,32768.d0)
  n1=(dn-ng)/32768.d0

  return
end subroutine packtext2


code:
program twq

  character*22 msg0,msg
  integer*1 data0(11)

  open(10,file='wqmsg.txt',status='old')
  write(*,1000)
1000 format(4x,'Encoded message',9x,'Decoded as',12x,'itype'/55('-'))

  do line=1,9999
     read(10,*,end=999) msg0
     call wqenc(msg0,itype,data0)
     call wqdec(data0,msg,ntype)
     write(*,1100) line,msg0,msg,ntype
1100 format(i2,'.',1x,a22,2x,a22,i3)
  enddo

999 end program twq


code:
subroutine unpack50(dat,n1,n2)

  integer*1 dat(11)

  i=dat(1)
  i4=iand(i,255)
  n1=ishft(i4,20)
  i=dat(2)
  i4=iand(i,255)
  n1=n1 + ishft(i4,12)
  i=dat(3)
  i4=iand(i,255)
  n1=n1 + ishft(i4,4)
  i=dat(4)
  i4=iand(i,255)
  n1=n1 + iand(ishft(i4,-4),15)
  n2=ishft(iand(i4,15),18)
  i=dat(5)
  i4=iand(i,255)
  n2=n2 + ishft(i4,10)
  i=dat(6)
  i4=iand(i,255)
  n2=n2 + ishft(i4,2)
  i=dat(7)
  i4=iand(i,255)
  n2=n2 + iand(ishft(i4,-6),3)

  return
end subroutine unpack50


code:
subroutine unpackname(n1,n2,name,len)

  character*9 name
  real*8 dn

  dn=32768.d0*n1 + n2
  len=0
  do i=9,1,-1
     j=mod(dn,27.d0)
     if(j.ge.1) then
        name(i:i)=char(64+j)
        len=len+1
     else
        name(i:i)=' '
     endif
     dn=dn/27.d0
  enddo

  return
end subroutine unpackname


code:
subroutine unpackpfx(ng,call1)

  character*12 call1
  character*3 pfx

  if(ng.lt.60000) then
! Add-on prefix of 1 to 3 characters
     n=ng
     do i=3,1,-1
        nc=mod(n,37)
        if(nc.ge.0 .and. nc.le.9) then
           pfx(i:i)=char(nc+48)
        else if(nc.ge.10 .and. nc.le.35) then
           pfx(i:i)=char(nc+55)
        else
           pfx(i:i)=' '
        endif
        n=n/37
     enddo
     call1=pfx//'/'//call1
     if(call1(1:1).eq.' ') call1=call1(2:)
     if(call1(1:1).eq.' ') call1=call1(2:)
  else
! Add-on suffix, one character
     i1=index(call1,' ')
     nc=ng-60000
     if(nc.ge.0 .and. nc.le.9) then
        call1=call1(:i1-1)//'/'//char(nc+48)
     else if(nc.ge.10 .and. nc.le.35) then
        call1=call1(:i1-1)//'/'//char(nc+55)
     endif
  endif

  return
end subroutine unpackpfx


code:
subroutine unpackprop(n1,k,muf,ccur,cxp)

  character ccur*4,cxp*2

  muf=mod(n1,62)
  n1=n1/62

  k=mod(n1,11)
  n1=n1/11

  j=mod(n1,53)
  n1=n1/53
  if(j.eq.0) cxp='*'
  if(j.ge.1 .and. j.le.26) cxp=char(64+j)
  if(j.gt.26) cxp=char(64+j-26)//char(64+j-26)

  j=mod(n1,53)
  n1=n1/53
  if(j.eq.0) ccur(2:2)='*'
  if(j.ge.1 .and. j.le.26) ccur(2:2)=char(64+j)
  if(j.gt.26) ccur(2:3)=char(64+j-26)//char(64+j-26)
  j=n1
  if(j.eq.0) ccur(1:1)='*'
  if(j.ge.1 .and. j.le.26) ccur(1:1)=char(64+j)
  if(j.gt.26) ccur=char(64+j-26)//char(64+j-26)//ccur(2:3)

  return
end subroutine unpackprop


code:
subroutine unpacktext2(n1,ng,msg)

  character*22 msg
  real*8 dn
  character*41 c
  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +./?'/

  msg='                      '
  dn=32768.d0*n1 + ng
  do i=8,1,-1
     j=mod(dn,41.d0)
     msg(i:i)=c(j+1:j+1)
     dn=dn/41.d0
  enddo

  return
end subroutine unpacktext2


code:
subroutine wqdec(data0,message,ntype)

  use packjt
  parameter (N15=32758)
  integer*1 data0(11)
  character*22 message
  character*12 callsign
  character*3 cdbm,cf
  character*2 crpt
  character*4 grid,psfx
  character*9 name
  character*36 fmt
  character*6 cwx(4)
  character*7 cwind(5)
  character ccur*4,cxp*2
  logical first
  character*12 dcall(0:N15-1)
  data first/.true./
  data cwx/'CLEAR','CLOUDY','RAIN','SNOW'/
  data cwind/'CALM','BREEZES','WINDY','DRY','HUMID'/
  save first,dcall

  if(first) then
     dcall='            '
     first=.false.
  endif

  message='                      '
  call unpack50(data0,n1,n2)
  call unpackcall(n1,callsign,iv2,psfx)
  i1=index(callsign,' ')
  call unpackgrid(n2/128,grid)
  ntype=iand(n2,127) -64

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  nu=mod(ntype,10)
  if(ntype.ge.0 .and. ntype.le.60 .and. (nu.eq.0 .or. nu.eq.3 .or.   &
       nu.eq.7)) then
     write(cdbm,'(i3)'),ntype
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     if(cdbm(1:1).eq.' ') cdbm=cdbm(2:)
     message=callsign(1:i1)//grid//' '//cdbm
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1)

! "Best DX" WSPR response (type 1)
  else if(ntype.eq.1) then
     message=grid//' DE '//callsign

! CQ (msg 3; types 2,4,5)
  else if(ntype.eq.2) then
     message='CQ '//callsign(:i1)//grid
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1)
  else if(ntype.eq.4 .or. ntype.eq.5) then
     ng=n2/128 + 32768*(ntype-4)
     call unpackpfx(ng,callsign)
     message='CQ '//callsign
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1)

! Reply to CQ (msg #2; type 6)
  else if(ntype.eq.6) then
     ih=(n2-64-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message='<'//dcall(ih)(:i2-1)//'> '//callsign(:i1-1)
     else
        message='<...> '//callsign
     endif
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Reply to CQ (msg #2; type 8)
  else if(ntype.eq.8) then
     message='DE '//callsign(:i1)//grid
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Reply to CQ, DE pfx/call (msg #2; types 9, 11)
  else if(ntype.eq.9 .or. ntype.eq.11) then
     ng=n2/128 + 32768*(ntype-9)/2
     call unpackpfx(ng,callsign)
     message='DE '//callsign
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Calls and report (msg #3; types -1 to -9)
  else if(ntype.le.-1 .and. ntype.ge.-9) then
     write(crpt,1010) -ntype
1010 format('S',i1)
     ih=(n2-62-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message=callsign(:i1)//'<'//dcall(ih)(:i2-1)//'> '//crpt
     else
        message=callsign(:i1)//'<...> '//crpt
     endif
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! pfx/call and report (msg #3; types -10 to -27)
  else if(ntype.le.-10 .and. ntype.ge.-27) then
     ng=n2/128
     nrpt=-ntype-9
     if(ntype.le.-19) then
        ng=ng + 32768
        nrpt=-ntype-18
     endif
     write(crpt,1010) nrpt
     call unpackpfx(ng,callsign)
     message=callsign//' '//crpt
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Calls and R and report (msg #4; types -28 to -36)
  else if(ntype.le.-28 .and. ntype.ge.-36) then
     write(crpt,1010) -(ntype+27)
     ih=(n2-64+28-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message=callsign(:i1)//'<'//dcall(ih)(:i2-1)//'> '//'R '//crpt
     else
        message=callsign(:i1)//'<...> '//'R '//crpt
     endif
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! pfx/call R and report (msg #4; types -37 to -54)
  else if(ntype.le.-37 .and. ntype.ge.-54) then
     ng=n2/128
     nrpt=-ntype-36
     if(ntype.le.-46) then
        ng=ng + 32768
        nrpt=-ntype-45
     endif
     write(crpt,1010) nrpt
     call unpackpfx(ng,callsign)
     message=callsign//' R '//crpt
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Calls and RRR (msg#5; type 12)
  else if(ntype.eq.12) then
     ih=(n2-64+28-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message=callsign(:i1)//'<'//dcall(ih)(:i2-1)//'> RRR'
     else
        message=callsign(:i1)//'<...> RRR'
     endif
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! Calls and RRR (msg#5; type 14)
  else if(ntype.eq.14) then
     ih=(n2-64+28-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message='<'//dcall(ih)(:i2-1)//'> '//callsign(:i1)//'RRR'
     else
        message='<...> '//callsign(:i1)//' RRR'
     endif
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! DE pfx/call and RRR (msg#5; types 15, 16)
  else if(ntype.eq.15 .or. ntype.eq.16) then
     ng=n2/128 + 32768*(ntype-15)
     call unpackpfx(ng,callsign)
     message='DE '//callsign//' RRR'
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! TNX [name] 73 GL (msg #6; type 18)
  else if(ntype.eq.18) then
     ng=(n2-18-64)/128
     call unpackname(n1,ng,name,len)
     message='TNX '//name(:len)//' 73 GL'

! OP [name] 73 GL (msg #6; type 18)
  else if(ntype.eq.-56) then
     ng=(n2+56-64)/128
     call unpackname(n1,ng,name,len)
     message='OP '//name(:len)//' 73 GL'

! 73 DE [call] [grid] (msg #6; type 19)
  else if(ntype.eq.19) then
     ng=(n2-19-64)/128
     message='73 DE '//callsign(:i1)//grid
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! 73 DE pfx/call (msg #6; type 21, 22)
  else if(ntype.eq.21 .or. ntype.eq.22) then
     ng=n2/128 + (ntype-21)*32768
     call unpackpfx(ng,callsign)
     i1=index(callsign,' ')
     message='73 DE '//callsign
     call hash(callsign,i1-1,ih)
     dcall(ih)=callsign(:i1-1)

! [power] W [gain] DBD 73 GL (msg#6; type 24, 25)
  else if(ntype.eq.24 .or. ntype.eq.25) then
     ng=(n2-24-64)/128 - 32
     i1=1
     if(n1.gt.0) i1=log10(float(n1)) + 1
     i2=1
     if(ng.ge.10) i2=2
     if(ng.lt.0) i2=i2+1
     if(n1.le.3000) then
        if(ntype.eq.24) fmt="(i4,' W ',i2,' DBD 73 GL')"
        if(ntype.eq.25) fmt="(i4,' W ',i2,' DBD      ')"
        fmt(3:3)=char(48+i1)
        fmt(12:12)=char(48+i2)
        if(ng.le.100) then
           write(message,fmt) n1,ng
        else
           if(ng.eq.30000) fmt=fmt(1:8)//"DIPOLE')"
           if(ng.eq.30001) fmt=fmt(1:8)//"VERTICAL')"
           write(message,fmt) n1
        endif
     else
        mw=n1-3000
        if(ntype.eq.24) fmt="('0.',i3.3,' W ',i2,' DBD 73 GL')"
        if(ntype.eq.25) fmt="('0.',i3.3,' W ',i2,' DBD      ')"
        fmt(19:19)=char(48+i2)
        if(ng.le.100) then
           write(message,fmt) mw,ng
        else
           if(ng.eq.30000) fmt=fmt(1:15)//"DIPOLE')"
           if(ng.eq.30001) fmt=fmt(1:15)//"VERTICAL')"
           write(message,fmt) n1
        endif
        if(index(message,'***').gt.0) go to 700
     endif

! QRZ call (msg #3; type 26)
  else if(ntype.eq.26) then
     ng=(n2-24-64)/128 - 32
     message='QRZ '//callsign

! PSE QSY [nnn] KHZ (msg #6; type 28)
  else if(ntype.eq.28) then
     if(n1.gt.0) i1=log10(float(n1)) + 1
     fmt="('PSE QSY ',i2,' KHZ')"
     fmt(14:14)=char(48+i1)
     write(message,fmt) n1

! WX wx temp C/F wind (msg #6; type 29)
  else if(ntype.eq.29) then
     nwx=n1/10000
     ntemp=mod(n1,10000) - 100
     cf=' F '
     if(ntemp.gt.800) then
        ntemp=ntemp-1000
        cf=' C '
     endif
     n2a=n2/128
     if(nwx.ge.1 .and. nwx.le.4 .and. n2a.ge.1 .and. n2a.le.5) then
        write(message,1020) cwx(nwx),ntemp,cf,cwind(n2/128)
1020    format('WX ',a6,i3,a3,a7)
     else
        message='WX'//' (BadMsg)'
     endif

! Hexadecimal data (type 62)
  else if(ntype.eq.62) then
     ng=n2/128
     write(message,'(z4.4,z7.7)') ng,n1

! Solar/geomagnetic/ionospheric data (type 63)
  else if(ntype.eq.63) then
     ih=(n2-64-ntype)/128
     if(dcall(ih)(1:1).ne.' ') then
        i2=index(dcall(ih),' ')
        message='<'//dcall(ih)(:i2-1)//'> '
     else
        message='<...> '
     endif
     call unpackprop(n1,k,muf,ccur,cxp)
     i2=index(message,'>')
     write(message(i2+1:),'(i3,i3)') k,muf
     message=message(:i2+7)//ccur//' '//cxp

! [plain text] (msg #6; type -57)
  else if(ntype.eq.-57) then
     ng=n2/128
     call unpacktext2(n1,ng,message)
  else
     go to 700
  endif
  go to 750

!     message='<Unknown message type>'
700 i1=index(callsign,' ')
  if(i1.lt.1) i1=12
  message=callsign(:i1)//' (BadMsg)'

750 do i=1,22
     if(ichar(message(i:i)).eq.0) message(i:i)=' '
  enddo

  do i=22,1,-1
     if(message(i:i).ne.' ') go to 800
  enddo
800 i2=i
  do n=1,20
     i1=index(message(:i2),'  ')
     if(i1.le.0) go to 900
     message=message(1:i1)//message(i1+2:)
     i2=i2-1
  enddo

900  return
end subroutine wqdec


code:
subroutine wqenc(msg,ntype,data0)

!  Parse and encode a WSPR message.

  use packjt
  parameter (MASK15=32767)
  character*22 msg
  character*12 call1,call2
  character*4 grid
  character*9 name
  character ccur*4,cxp*2
  logical lbad1,lbad2
  integer*1 data0(11)
  integer nu(0:9)
  data nu/0,-1,1,0,-1,2,1,0,-1,1/

  read(msg,1001,end=1,err=1) ng,n1
1001 format(z4,z7)
  ntype=62
  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)             !Pack 8 bits per byte, add tail
  go to 900

1  if(msg(1:6).eq.'73 DE ') go to 80
  if(index(msg,' W ').gt.0 .and. index(msg,' DBD ').gt.0) go to 90
  if(msg(1:4).eq.'QRZ ') go to 100
  if(msg(1:8).eq.'PSE QSY ') go to 110
  if(msg(1:3).eq.'WX ') go to 120

! Standard WSPR message (types 0 3 7 10 13 17 ... 60)
  i1=index(msg,' ')
  if(i1.lt.4 .or. i1.gt.7) go to 10
  call1=msg(:i1-1)
  grid=msg(i1+1:i1+4)
  call packcall(call1,n1,lbad1)
  call packgrid(grid,ng,lbad2)
  if(lbad1 .or. lbad2) go to 10
  ndbm=0
  read(msg(i1+5:),*,err=10,end=800) ndbm
  if(ndbm.lt.0 .or. ndbm.gt.60) go to 800
  ndbm=ndbm+nu(mod(ndbm,10))
  n2=128*ng + (ndbm+64)
  call pack50(n1,n2,data0)
  ntype=ndbm
  go to 900

! "BestDX" automated WSPR reply (type 1)
10 if(i1.ne.5 .or. msg(5:8).ne.' DE ') go to 20
  grid=msg(1:4)
  call packgrid(grid,ng,lbad2)
  if(lbad2) go to 800
  call1=msg(9:)
  call packcall(call1,n1,lbad1)
  if(lbad1) go to 800
  ntype=1
  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)             !Pack 8 bits per byte, add tail
  go to 900

! CQ (msg #1; types 2, 4, 5)
20  if(msg(1:3).ne.'CQ ') go to 30
  if(index(msg,'/').le.0) then
     i2=index(msg(4:),' ')
     call1=msg(4:i2+3)
     grid=msg(i2+4:)
     call packcall(call1,n1,lbad1)
     if(lbad1) go to 30
     call packgrid(grid,ng,lbad2)
     if(lbad2) go to 30
     ntype=2
     n2=128*ng + (ntype+64)
     call pack50(n1,n2,data0)
  else
     ntype=4                                     ! or 5
     call1=msg(4:)
     call packpfx(call1,n1,ng,nadd)
     ntype=ntype+nadd
     n2=128*ng + ntype + 64
     call pack50(n1,n2,data0)
  endif
  go to 900

! Reply to CQ (msg #2; types 6,8,9,11)
30 if(msg(1:1).ne.'<' .and. msg(1:3).ne.'DE ') go to 40
  if(index(msg,' RRR ').gt.0) go to 50
  if(msg(1:1).eq.'<') then
     ntype=6
     i1=index(msg,'>')
     call1=msg(2:i1-1)
     read(msg(i1+1:),*,err=31,end=31) k,muf,ccur,cxp
     go to 130
31   call2=msg(i1+2:)
     call hash(call1,i1-2,ih)
     call packcall(call2,n1,lbad1)
     n2=128*ih + (ntype+64)
     call pack50(n1,n2,data0)
  else
     i1=index(msg(4:),' ')
     call1=msg(4:i1+2)
     if(index(msg,'/').le.0) then
        ntype=8
        ih=0
        call packcall(call1,n1,lbad1)
        grid=msg(i1+4:i1+7)
        call packgrid(grid,ng,lbad2)
        n2=128*ng + (ntype+64)
        call pack50(n1,n2,data0)
     else
        ntype=9                                   ! or 11
        call1=msg(4:)
        call packpfx(call1,n1,ng,nadd)
        ntype=ntype + 2*nadd
        n2=128*ng + ntype + 64
        call pack50(n1,n2,data0)
     endif
  endif
  go to 900

! Call(s) + report (msg #3; types -1 to -27)
! Call(s) + R + report (msg #4; types -28 to -54)
40 if(index(msg,' RRR').gt.0) go to 50
  i1=index(msg,'<')
  if(i1.gt.0 .and. (i1.lt.5 .or. i1.gt.8)) go to 50
  i2=index(msg,'/')
  if(i2.gt.0 .and.i2.le.4) then
     ntype=-10                                   ! -10 to -27
     i0=index(msg,' ')
     call1=msg(:i0-1)
     call packpfx(call1,n1,ng,nadd)
     ntype=ntype - 9*nadd
     i2=index(msg,' ')
     i3=index(msg,' R ')
     if(i3.gt.0) i2=i2+2                            !-28 to -36
     read(msg(i2+2:i2+2),*,end=800,err=800) nrpt
     ntype=ntype - (nrpt-1)
     if(i3.gt.0) ntype=ntype-27
     n2=128*ng + ntype + 64
     call pack50(n1,n2,data0)
     go to 900
  else if(i1.eq.0) then
     go to 50
  endif
  call1=msg(:i1-2)                               !-1 to -9
  i2=index(msg,'>')
  call2=msg(i1+1:i2-1)
  call hash(call2,i2-i1-1,ih)
  i3=index(msg,' R ')
  if(i3.gt.0) i2=i2+2                            !-28 to -36
  read(msg(i2+3:i2+3),*,end=42,err=42) nrpt
  go to 43
42 nrpt=1
43 ntype=-nrpt
  if(i3.gt.0) ntype=-(nrpt+27)
  call packcall(call1,n1,lbad1)
  n2=128*ih + (ntype+64)
  call pack50(n1,n2,data0)
  go to 900

50 i0=index(msg,'<')
  if(i0.le.0 .and. msg(1:3).ne.'DE ') go to 60
  i3=index(msg,' RRR')
  if(i3.le.0) go to 60
! Call or calls and RRR (msg#5; type2 12,14,15,16)
  i0=index(msg,'<')
  if(i0.eq.1) then
     if(index(msg,'/').le.0) then
        ntype=14
        i1=index(msg,'>')
        call1=msg(2:i1-1)
        call2=msg(i1+2:)
        i2=index(call2,' ')
        call2=call2(:i2-1)
        call packcall(call2,n1,lbad1)
        call hash(call1,i1-2,ih)
        n2=128*ih + (ntype+64)
        call pack50(n1,n2,data0)
     else
        stop '0002'
     endif
  else if(i0.ge.5 .and. i0.le.8) then
     if(index(msg,'/').le.0) then
        ntype=12
        i1=index(msg,'>')
        call1=msg(:i0-2)
        call2=msg(i0+1:i1-1)
        call packcall(call1,n1,lbad1)
        call hash(call2,i1-i0-1,ih)
        n2=128*ih + (ntype+64)
        call pack50(n1,n2,data0)
     else
        stop '0002'
     endif
  else
     i1=index(msg(4:),' ')
     call1=msg(4:i1+2)
     if(index(msg,'/').le.0) then
        ntype=9
        grid=msg(i1+4:i1+7)
     else
        ntype=15                                   ! or 16
        call1=msg(4:)
        i0=index(call1,' ')
        call1=call1(:i0-1)
        call packpfx(call1,n1,ng,nadd)
        ntype=ntype+nadd
        n2=128*ng + ntype + 64
        call pack50(n1,n2,data0)
     endif
  endif
  go to 900

! TNX <name> 73 GL (msg #6; type 18 ...)
60 if(msg(1:4).ne.'TNX ') go to 70
  ntype=18
  n1=0
  i2=index(msg(5:),' ')
  name=msg(5:i2+4)
  call packname(name,i2-1,n1,ng)
  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)
  go to 900

! TNX name 73 GL (msg #6; type -56 ...)
70 if(msg(1:3).ne.'OP ') go to 80
  ntype=-56
  n1=0
  i2=index(msg(4:),' ')
  name=msg(4:i2+3)
  call packname(name,i2-1,n1,ng)
  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)
  go to 900

! 73 DE call grid (msg #6; type 19)
80 if(msg(1:6).ne.'73 DE ') go to 90
  ntype=19
  i1=index(msg(7:),' ')
  call1=msg(7:)
  if(index(call1,'/').le.0) then
     i1=index(call1,' ')
     grid=call1(i1+1:)
     call1=call1(:i1-1)
     call packcall(call1,n1,lbad1)
     call packgrid(grid,ng,lbad2)
     if(lbad1 .or. lbad2) go to 800
     n2=128*ng + (ntype+64)
     call pack50(n1,n2,data0)
     go to 900
  else
     ntype=21                                   ! or 22
     call packpfx(call1,n1,ng,nadd)
     ntype=ntype + nadd
     n2=128*ng + ntype + 64
     call pack50(n1,n2,data0)
     go to 900
  endif

! [pwr] W [gain] DBD [73 GL] (msg #6; types 24, 25)
90  if(index(msg,' W ').le.0) go to 140
  ntype=25
  if(index(msg,' DBD 73 GL').gt.0) ntype=24
  i1=index(msg,' ')
  read(msg(:i1-1),*,end=800,err=800) watts
  if(watts.ge.1.0) nwatts=watts
  if(watts.lt.1.0) nwatts=3000 + nint(1000.*watts)
  if(index(msg,'DIPOLE').gt.0) then
     ndbd=30000
  else if(index(msg,'VERTICAL').gt.0) then
     ndbd=30001
  else
     i2=index(msg(i1+3:),' ')
     read(msg(i1+3:i1+i2+1),*,end=800,err=800) ndbd
  endif
  n1=nwatts
  ng=ndbd + 32
  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)
  go to 900

! QRZ call (msg #3; type 26)
100 call1=msg(5:)
  call packcall(call1,n1,lbad1)
  if(lbad1) go to 800
  ntype=26
  n2=ntype+64
  call pack50(n1,n2,data0)
  go to 900

! PSE QSY [nnn] KHZ (msg #6; type 28)
110 ntype=28
  read(msg(9:),*,end=800,err=800) n1
  n2=ntype+64
  call pack50(n1,n2,data0)
  go to 900

! WX wx temp C|F wind (msg #6; type 29)
120 ntype=29
  if(index(msg,' CLEAR ').gt.0) then
     i1=10
     n1=10000
  else if(index(msg,' CLOUDY ').gt.0) then
     i1=11
     n1=20000
  else if(index(msg,' RAIN ').gt.0) then
     i1=9
     n1=30000
  else if(index(msg,' SNOW ').gt.0) then
     i1=9
     n1=40000
  endif
  read(msg(i1:),*,err=800,end=800) ntemp
  ntemp=ntemp+100
  i1=index(msg,' C ')
  if(i1.gt.0) ntemp=ntemp+1000
  n1=n1+ntemp
  if(index(msg,' CALM').gt.0) ng=1
  if(index(msg,' BREEZES').gt.0) ng=2
  if(index(msg,' WINDY').gt.0) ng=3
  if(index(msg,' DRY').gt.0) ng=4
  if(index(msg,' HUMID').gt.0) ng=5

  n2=128*ng + (ntype+64)
  call pack50(n1,n2,data0)

  go to 900

! Solar/geomagnetic/ionospheric data
130 ntype=63
  call packprop(k,muf,ccur,cxp,n1)
  call hash(call1,i1-2,ih)
  n2=128*ih + ntype + 64
  call pack50(n1,n2,data0)
  go to 900

140 continue

! Plain text
800 ntype=-57
  call packtext2(msg(:8),n1,ng)
  n2=128*ng + ntype + 64
  call pack50(n1,n2,data0)
  go to 900

900 continue
  return
end subroutine wqenc


code:
"CQ K1JT FN20"
"CQ PJ4/K1JT"
"<K1JT> W6CQZ"
"DE W6CQZ CM87"
"DE PJ4/K1JT"
"W6CQZ <K1JT> S4"
"QRZ K1JT"
"PJ4/W6CQZ S4"
"K1JT <W6CQZ> R S3"
"PJ4/K1JT R S3"
"<W6CQZ> K1JT RRR"
"W6CQZ <K1JT> RRR"
"DE PJ4/K1JT RRR"
"73 DE W6CQZ CM87"
"73 DE PJ4/K1JT"
"TNX VICTORIA 73 GL"
"OP HARRY 73 GL"
"5 W DIPOLE"
"10 W VERTICAL"
"1 W 0 DBD"
"1500 W 21 DBD 73 GL"
"PSE QSY 1811 KHZ"
"WX SNOW -5 C CALM"
"CUL JACK"
"."
"CQ K1JT FN20"
"<K1JT> W6CQZ"
"W6CQZ <K1JT> S4"
"K1JT <W6CQZ> R S3"
"<W6CQZ> K1JT RRR"
"TNX JOE 73 GL"


code:
!
! readwav - open and read the header of a WAV format file
!
! On successful exit the file is left positioned at the start of the
! data.
!
! Example of usage:
!
!  use readwav
!  integer*2 sample
!  type(wav_header) wav
!  call wav%read ('file.wav')
!  write (*,*) 'Sample rate is: ', wav%audio_format%sample_rate
!  do i=0,wav%data_size
!    read (unit=wav%lun) sample
!    ! process sample
!  end do
!
module readwav
  implicit none

  type format_chunk
     integer*2 audio_format
     integer*2 num_channels
     integer sample_rate
     integer byte_rate
     integer*2 block_align
     integer*2 bits_per_sample
  end type format_chunk

  type, public :: wav_header
     integer :: lun
     type(format_chunk) :: audio_format
     integer :: data_size
   contains
     procedure :: read
  end type wav_header

  private
contains
  subroutine read (this, filename)
    implicit none

    type riff_descriptor
       character(len=4) :: id
       integer :: size
    end type riff_descriptor

    class(wav_header), intent(inout) :: this
    character(len=*), intent(in) :: filename

    integer :: filepos
    type(riff_descriptor) :: desc
    character(len=4) :: riff_type

    open (newunit=this%lun, file=filename, access='stream', form='unformatted', status='old')
    read (unit=this%lun) desc,riff_type
    inquire (unit=this%lun, pos=filepos)
    do
       read (unit=this%lun, pos=filepos) desc
       inquire (unit=this%lun, pos=filepos)
       if (desc%id .eq. 'fmt ') then
          read (unit=this%lun) this%audio_format
       else if (desc%id .eq. 'data') then
          this%data_size = desc%size
          exit
       end if
       filepos = filepos + (desc%size + 1) / 2 * 2 ! pad to even alignment
    end do
  end subroutine read
end module readwav


code:
subroutine rectify_msk(c,msg0,imsg,freq2)

  parameter (NSPM=1404)
  complex c(0:NSPM-1)                         !Received data
  complex cmsg(0:NSPM-1)                      !Message waveform
  complex c1(0:NSPM-1)                        !Rectified signal
  complex c2(0:NSPM-1)                        !Integral of rectified signal
  complex c3(0:2*NSPM-1)                      !FFT of rectified signal
  complex cfac
  character*22 msg0,msg,msgsent
  integer i4tone(234)

  ichk=0
  msg=msg0
  nsym=234
  if(imsg.ge.0) then
     ichk=10000+imsg
     msg="<C1ALL C2ALL> 73"
     nsym=35
  endif
  call genmsk(msg,ichk,msgsent,i4tone,itype)  !Get tone sequence for msg

  twopi=8.0*atan(1.0)
  dt=1.0/12000.0
  f0=1000.0
  f1=2000.0
  phi=0.
  dphi=0.
  k=-1
  c2=0.
  do j=1,nsym                                  !Generate Tx waveform for msg
     if(i4tone(j).eq.0) dphi=twopi*f0*dt
     if(i4tone(j).eq.1) dphi=twopi*f1*dt
     do i=1,6
        k=k+1
        phi=phi+dphi
        cmsg(k)=cmplx(cos(phi),sin(phi))
        c1(k)=conjg(cmsg(k))*c(k)
        if(k.ge.1) c2(k)=c2(k-1) + c1(k)
     enddo
  enddo
  c2(0)=c2(1)
  pha=atan2(aimag(c2(NSPM-1)),real(c2(NSPM-1)))
  cfac=cmplx(cos(pha),-sin(pha))
  c1=cfac*c1
  c2=cfac*c2
  nfft=2*NSPM
  c3(0:NSPM-1)=c2
  c3(NSPM:nfft-1)=0.
  df=12000.0/nfft
  call four2a(c3,nfft,1,-1,1)
  smax=0.
  do i=0,nfft-1
     f=i*df
     if(i.gt.nfft/2) f=f-12000.0
     s=1.e-10*(real(c3(i))**2 + aimag(c3(i))**2)
     if(s.gt.smax) then
        smax=s
        freq2=1500.0 + f
     endif
  enddo

  return
end subroutine rectify_msk


code:
subroutine refspectrum(id2,bclear,brefspec,buseref,fname)

! Input:
!  id2       i*2        Raw 16-bit integer data, 12000 Hz sample rate
!  brefspec  logical    True when accumulating a reference spectrum

  parameter (NFFT=6912,NH=NFFT/2,NPOLYLOW=400,NPOLYHIGH=2600)
  integer*2 id2(NFFT)
  logical*1 bclear,brefspec,buseref,blastuse

  real xs(0:NH-1)                         !Saved upper half of input chunk convolved with h(t)
  real x(0:NFFT-1)                        !Work array
  real*4 w(0:NFFT-1)                      !Window function
  real*4 s(0:NH)                          !Average spectrum
  real*4 fil(0:NH)
  real*8 xfit(1500),yfit(1500),sigmay(1500),a(5),chisqr !Polyfit arrays
  logical first
  complex cx(0:NH)                        !Complex frequency-domain work array
  complex cfil(0:NH)
  character*(*) fname
  common/spectra/syellow(6827),ref(0:NH),filter(0:NH)
  equivalence(x,cx)
  data first/.true./,blastuse/.false./
  save

  if(first) then
     pi=4.0*atan(1.0)
     do i=0,NFFT-1
        ww=sin(i*pi/NFFT)
        w(i)=ww*ww/NFFT
     enddo
     nsave=0
     s=0.0
     filter=1.0
     xs=0.
     first=.false.
  endif
  if(bclear) s=0.

  if(brefspec) then
     x(0:NH-1)=0.001*id2(1:NH)
     x(NH:NFFT-1)=0.0
     call four2a(x,NFFT,1,-1,0)                 !r2c FFT

     do i=1,NH
        s(i)=s(i) + real(cx(i))**2 + aimag(cx(i))**2
     enddo
     nsave=nsave+1

     fac0=0.9
     if(mod(nsave,4).eq.0) then
        df=12000.0/NFFT
        ia=nint(1000.0/df)
        ib=nint(2000.0/df)
        avemid=sum(s(ia:ib))/(ib-ia+1)
        do i=0,NH
           fil(i)=0.
           if(s(i).gt.0.0) then
              fil(i)=sqrt(avemid/s(i))
           endif
        enddo

! Default range is 240 - 4000 Hz.  For narrower filters, use frequencies
! at which gain is -20 dB relative to 1500 Hz.
        ia=nint(240.0/df)
        ib=nint(4000.0/df)
        i0=nint(1500.0/df)
        do i=i0,ia,-1
           if(s(i)/s(i0).lt.0.01) exit
        enddo
        ia=i
        do i=i0,ib,1
           if(s(i)/s(i0).lt.0.01) exit
        enddo
        ib=i

        fac=fac0
        do i=ia,1,-1
           fac=fac*fac0
           fil(i)=fac*fil(i)
        enddo

        fac=fac0
        do i=ib,NH
           fac=fac*fac0
           fil(i)=fac*fil(i)
        enddo

        do iter=1,100                        !### ??? ###
           call smo121(fil,NH)
        enddo

        do i=0,NH
           filter(i)=-60.0
           if(s(i).gt.0.0) filter(i)=20.0*log10(fil(i))
        enddo

        il=nint(NPOLYLOW/df)
        ih=nint(NPOLYHIGH/df)
        nfit=ih-il+1
        mode=0
        nterms=5
        do i=1,nfit
          xfit(i)=((i+il-1)*df-1500.0)/1000.0
          yfit(i)=fil(i+il-1)
          sigmay(i)=1.0
        enddo
        call polyfit(xfit,yfit,sigmay,nfit,nterms,mode,a,chisqr)

        open(16,file=fname,status='unknown')
        write(16,1003) NPOLYLOW,NPOLYHIGH,nterms,a
1003    format(3i5,5e25.16)
        do i=1,NH
           freq=i*df
           ref(i)=db(s(i)/avemid)
           write(16,1005) freq,s(i),ref(i),fil(i),filter(i)
1005       format(f10.3,e12.3,f12.6,e12.3,f12.6)
        enddo
        close(16)
     endif
     return
  endif

  if(buseref) then
     if(blastuse.neqv.buseref) then !just enabled so read filter
        fil=1.0
        open(16,file=fname,status='old',err=110)
        read(16,1003,err=20,end=100) ndummy,ndummy,nterms,a
        goto 30
20      rewind(16)              !allow for old style refspec.dat with no header
30      do i=1,NH
           read(16,1005,err=100,end=100) freq,s(i),ref(i),fil(i),filter(i)
        enddo
! Make the filter causal for overlap and add.
        cx(0)=0.0
        cx(1:NH)=fil(1:NH)/NFFT
        call four2a(x,NFFT,1,1,-1)
        x=cshift(x,-400)
        x(800:NH)=0.0
        call four2a(cx,NFFT,1,-1,0)
        cfil=cx
100     close(16)
110     continue
     endif
! Use overlap and add method to apply causal reference filter.
     x(0:NH-1)=id2(1:NH)
     x(NH:NFFT-1)=0.0
     x=x/NFFT
     call four2a(x,NFFT,1,-1,0)
     cx=cfil*cx
     call four2a(cx,NFFT,1,1,-1)
     x(0:NH-1)=x(0:NH-1)+xs
     xs=x(NH:NFFT-1)
     id2(1:NH)=nint(x(0:NH-1))
  endif
  blastuse=buseref

  return
end subroutine refspectrum


code:
#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "tstrig.h"

RIG *my_rig;                                       // handle to rig

int set_conf(RIG *my_rig, char *conf_parms);

//------------------------------------------------------------------------
int set_conf(RIG *my_rig, char *conf_parms)
{
  char *p, *q, *n;
  int iret;

  p = conf_parms;
  while (p && *p != '\0') {
    /* FIXME: left hand value of = cannot be null */
    q = strchr(p, '=');
    if ( !q )
      return -RIG_EINVAL;
    *q++ = '\0';
    n = strchr(q, ',');
    if (n) *n++ = '\0';
    iret = rig_set_conf(my_rig, rig_token_lookup(my_rig, p), q);
    if (iret != RIG_OK)
      return iret;
    p = n;
  }
  return RIG_OK;
}

//------------------------------------------------------------------------
int rigOpen(int verbose, rig_model_t my_model, const char* rig_file,
	    int serial_rate, const char* conf_parms2)
{
  int iret;		/* generic return code from functions */
  char *civaddr = NULL;	/* NULL means no need to set conf */
  //  const char *rig_file;
  //  const char *conf_parms2;
  //  int serial_rate;

  rig_set_debug(verbose);
  my_rig=rig_init(my_model);

  if (!my_rig) {
    //    fprintf(stderr, "Unknown rig num %d, or initialization error.\n",my_model);
    return -1;
  }

  //  rig_file="COM1";
  //  serial_rate=4800;
  //  conf_parms2="data_bits=8,stop_bits=2,serial_handshake=Hardware";

  iret=set_conf(my_rig, conf_parms2);
  if (iret!=RIG_OK) {
    //    fprintf(stderr, "Config parameter error: %s\n", rigerror(iret));
    return -2;
  }

  if (rig_file)
    strncpy(my_rig->state.rigport.pathname, rig_file, FILPATHLEN - 1);

  if (serial_rate!=0)
    my_rig->state.rigport.parm.serial.rate = serial_rate;

  if (civaddr)
    rig_set_conf(my_rig, rig_token_lookup(my_rig, "civaddr"), civaddr);

  iret = rig_open(my_rig);
  if(iret!=0) return -3;
  return 0;
}

int rigSetFreq(int fHz)
{
  return rig_set_freq(my_rig,RIG_VFO_CURR,fHz);
}

int rigFreq(int *fHz)
{
  int iret=0;
  freq_t freq;
  iret=rig_get_freq(my_rig, RIG_VFO_CURR, &freq);
  *fHz=freq;
  return iret;
}

int rigSetPTT(int iptt)
{
  return rig_set_ptt(my_rig, RIG_VFO_CURR, iptt);
}

void rigClose()
{
  rig_close(my_rig);
  rig_cleanup(my_rig);
}


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, 8-bit symbols */
void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity);
int decode_rs_char(void *rs,unsigned char *data,int *eras_pos,
		   int no_eras);
void *init_rs_char(int symsize,int gfpoly,
		   int fcr,int prim,int nroots,
		   int pad);
void free_rs_char(void *rs);

/* 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);
void *init_rs_int(int symsize,int gfpoly,int fcr,
		  int prim,int nroots,int pad);
void free_rs_int(void *rs);

/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis)
 * symbol representation
 */
void encode_rs_8(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras,int pad);

/* CCSDS standard (255,223) RS codec with dual-basis symbol representation */
void encode_rs_ccsds(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras,int pad);

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


code:
integer function savec2(c2name,ntrseconds,f0m1500)

! Array c0() has complex samples at 1500 Hz sample rate.
! WSPR-2:  downsample by 1/4 to produce c2, centered at 1500 Hz
! WSPR-15: downsample by 1/32 to produce c2, centered at 1612.5 Hz

  parameter (NDMAX=120*1500)         !Sample intervals at 1500 Hz rate
  parameter (MAXFFT=256*1024)

  character*(*) c2name
  character*14 outfile
  real*8 f0m1500
  complex c0
  complex c1(0:MAXFFT-1)
  complex c2(0:65535)
  common/c0com/c0(0:NDMAX-1)

  ntrminutes=ntrseconds/60
  npts=114*1500
  nfft1=262144
  if(ntrminutes.eq.15) then
     npts=890*1500
     nfft1=MAXFFT
  endif
  df1=1500.0/nfft1
  fac=1.0/nfft1
  c1(0:npts-1)=fac*c0(0:npts-1)
  c1(npts:nfft1-1)=0.

  call four2a(c1,nfft1,1,1,1)                 !Complex FFT to frequency domain

! Select the desired frequency range
  nfft2=65536
  nh2=nfft2/2
  if(ntrminutes.eq.2) then
     c2(0:nh2)=c1(0:nh2)
     c2(nh2+1:nfft2-1)=c1(nfft1-nh2+1:nfft1-1)
  else
     i0=nint(112.5/df1)
     c2(0:nh2)=c1(i0:i0+nh2)
     c2(nh2+1:nfft2-1)=c1(i0-nh2+1:i0-1)
  endif

  call four2a(c2,nfft2,1,-1,1)      !Shorter complex FFT, back to time domain

! Write complex time-domain data to disk.
  i1=index(c2name,'.c2')
  outfile=c2name(i1-11:i1+2)
  open(18,file=c2name,status='unknown',access='stream', iostat=ioerr)
  if (ioerr.eq.0) then
     write(18) outfile,ntrminutes,f0m1500,c2(0:45000-1)
     close(18)
  endif
  savec2 = ioerr
end function savec2


code:
real function sec_midn()
  sec_midn=secnds(0.0)
  return
end function sec_midn

subroutine sleep_msec(n)

  call usleep(1000*n)

  return
end subroutine sleep_msec


code:
subroutine setup65

! Defines arrays related to the JT65 pseudo-random synchronizing pattern.
! Executed at program start.

  use jt65_mod
  integer nprc(126)

! JT65
  data nprc/                                   &
      1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
      0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
      0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
      0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
      1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
      0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
      1,1,1,1,1,1/
  data mr2/0/                !Silence compiler warning

! Put the appropriate pseudo-random sequence into pr
  nsym=126
  do i=1,nsym
     pr(i)=2*nprc(i)-1
  enddo

! Determine locations of data and reference symbols
  k=0
  mr1=0
  do i=1,nsym
     if(pr(i).lt.0.0) then
        k=k+1
        mdat(k)=i
     else
        mr2=i
        if(mr1.eq.0) mr1=i
     endif
  enddo
  nsig=k

! Determine the reference symbols for each data symbol.
  do k=1,nsig
     m=mdat(k)
     mref(k,1)=mr1
     do n=1,10                     !Get ref symbol before data
        if((m-n).gt.0) then
           if (pr(m-n).gt.0.0) go to 10
        endif
     enddo
     go to 12
10   mref(k,1)=m-n
12   mref(k,2)=mr2
     do n=1,10                     !Get ref symbol after data
        if((m+n).le.nsym) then
           if (pr(m+n).gt.0.0) go to 20
        endif
     enddo
     cycle
20   mref(k,2)=m+n
  enddo

! Now do it all again, using opposite logic on pr(i)
  k=0
  mr1=0
  do i=1,nsym
     if(pr(i).gt.0.0) then
        k=k+1
        mdat2(k)=i
     else
        mr2=i
        if(mr1.eq.0) mr1=i
     endif
  enddo
  nsig=k

  do k=1,nsig
     m=mdat2(k)
     mref2(k,1)=mr1
     do n=1,10
        if((m-n).gt.0) then
           if (pr(m-n).lt.0.0) go to 110
        endif
     enddo
     go to 112
110  mref2(k,1)=m-n
112  mref2(k,2)=mr2
     do n=1,10
        if((m+n).le.nsym) then
           if (pr(m+n).lt.0.0) go to 120
        endif
     enddo
     cycle
120  mref2(k,2)=m+n
  enddo

  return
end subroutine setup65