code:
# Makefile for MinGW on Windows
#
# Needed libraries are located using the following variables:
#
#   QT_DIR - point this at the root of your Qt installation
#   FFTW3_DIR - point this at the fftw v3 installation directory
#
# e.g.
#   make -f Makefile.MinGW \
#    QT_DIR=c:/Qt/5.2.1/mingw48_32 \
#    FFTW3_DIR = c:/fftw-3.3.3-dll32-2
#
# Windows re-direct:
#   C> make > junk1 2>&1

# Set paths
EXE_DIR = ..\\..\\wsjtx_install
QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
FFTW3_DIR = ..

INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include

# Compilers
CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
RANLIB = ranlib
MKDIR = mkdir -p
CP = cp
RM = rm -f

FFLAGS = -O2 -fbounds-check -Wall -Wno-precision-loss -fno-second-underscore
CFLAGS = -I. -fbounds-check -mno-stack-arg-probe

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

all:    libjt9.a libastro.a jt9.exe jt9code.exe jt65code.exe

OBJS1 = prog_args.o options.o pctile.o graycode.o sort.o ssort.o chkmsg.o \
	unpackmsg.o igray.o unpackcall.o unpackgrid.o \
	fmtmsg.o grid2k.o unpacktext.o getpfx2.o packmsg.o deg2grid.o \
	packtext.o getpfx1.o packcall.o k2grid.o packgrid.o \
	nchar.o four2a.o grid2deg.o pfxdump.o f77_wisdom.o \
	symspec.o analytic.o db.o genjt9.o flat1.o smo.o \
	packbits.o unpackbits.o encode232.o interleave9.o \
	entail.o fano232.o gran.o sync9.o decode9.o \
	fil3.o decoder.o grid2n.o n2grid.o timer.o \
	softsym.o getlags.o afc9.o fchisq.o twkfreq.o downsam9.o \
	peakdt9.o symspec2.o stdmsg.o morse.o azdist.o geodist.o \
	fillcom.o chkss2.o zplot9.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat3.o polfit.o determ.o baddata.o

libjt9.a: $(OBJS1)
	$(AR) libjt9.a $(OBJS1)
	$(RANLIB) libjt9.a

OBJS2 = jt9.o jt9a.o jt9b.o jt9c.o ipcomm.o sec_midn.o usleep.o
LIBS2 = -L${QT_DIR}/lib -lQt5Core
jt9.exe: $(OBJS2) libjt9.a
	$(CXX) -o jt9.exe $(OBJS2) $(LIBS2) libjt9.a \
	-L$(FFTW3_DIR) -lfftw3f-3 $(shell $(FC) -print-file-name=lib$(FC).a)
	-$(MKDIR) $(EXE_DIR)
	$(CP) jt9.exe $(EXE_DIR)

OBJS3 = jt9sim.o
jt9sim.exe: $(OBJS3) libjt9.a
	$(FC) -o jt9sim.exe $(OBJS3) libjt9.a

OBJS4 = jt9code.o
jt9code.exe: $(OBJS4) libjt9.a
	$(FC) -o jt9code.exe $(OBJS4) libjt9.a
	$(CP) jt9code.exe $(EXE_DIR)

OBJS5 = jt65.o
jt65.exe: $(OBJS5) libjt9.a
	$(FC) -o jt65.exe $(OBJS5) libjt9.a -L$(FFTW3_DIR) -lfftw3f-3

OBJS7 = astrosub.o astro0.o astro.o tm2.o grid2deg.o sun.o moondop.o \
	coord.o dot.o moon2.o tmoonsub.o toxyz.o geocentric.o \
	dcoord.o

libastro.a: $(OBJS7)
	$(AR) libastro.a $(OBJS7)
	$(RANLIB) libastro.a

OBJS6 = jt65code.o
jt65code.exe: $(OBJS6) libjt9.a
	$(FC) -o jt65code.exe $(OBJS6) libjt9.a
	$(CP) jt65code.exe $(EXE_DIR)

sync9.o: sync9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c sync9.f90

spec9.o: spec9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c spec9.f90

peakdt9.o: peakdt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c peakdt9.f90

jt9sim.o: jt9sim.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c jt9sim.f90

genjt9.o: genjt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c genjt9.f90

redsync.o: redsync.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c redsync.f90

unpackmsg.o: unpackmsg.f90
	$(FC) -c -O0 -fbounds-check -Wall -Wno-precision-loss unpackmsg.f90

ipcomm.o: ipcomm.cpp
	$(CXX) -c $(INCPATH) ipcomm.cpp

sec_midn.o: sec_midn.f90
	$(FC) -c -fno-second-underscore sec_midn.f90

#rig_control.o: rig_control.c
#	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include rig_control.c

tstrig.o: tstrig.c
	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include tstrig.c

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	$(RM) *.o libjt9.a wsjtx.exe jt9sim.exe jt9.exe jt65.exe


code:
# Makefile for MinGW on Windows
# Windows re-direct:
#   C> make > junk1 2>&1

# Set paths
EXE_DIR = ../../wsjtx_install
INCPATH = -I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/include/QtCore' \
	 -I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/include' \
	 -I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/include/ActiveQt' \
	 -I'release' -I'.' -I'c:/QtSDK/Desktop/Qt/4.7.4/mingw/mkspecs/win32-g++'

# Compilers
CC = gcc
CXX = g++
FC = g95

FFLAGS = -O2 -fbounds-check -Wall -Wno-precision-loss -fno-second-underscore
CFLAGS = -I. -fbounds-check -mno-stack-arg-probe

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

all:    libjt9.a jt9sim.exe jt9.exe jt9code.exe jt65.exe

OBJS1 = pctile.o graycode.o sort.o ssort.o \
	unpackmsg.o igray.o unpackcall.o unpackgrid.o \
	grid2k.o unpacktext.o getpfx2.o packmsg.o deg2grid.o \
	packtext.o getpfx1.o packcall.o k2grid.o packgrid.o \
	nchar.o four2a.o grid2deg.o pfxdump.o f77_wisdom.o \
	symspec.o analytic.o db.o genjt9.o \
	packbits.o unpackbits.o encode232.o interleave9.o \
	entail.o fano232.o gran.o sync9.o decode9.o \
	fil3.o decoder.o grid2n.o n2grid.o timer.o \
	softsym.o getlags.o afc9.o fchisq.o twkfreq.o downsam9.o \
	peakdt9.o symspec2.o stdmsg.o morse.o azdist.o geodist.o \
	fillcom.o chkss2.o zplot9.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o

libjt9.a: $(OBJS1)
	ar cr libjt9.a $(OBJS1)
	ranlib libjt9.a

OBJS2 = jt9.o jt9a.o jt9b.o jt9c.o ipcomm.o sec_midn.o usleep.o
LIBS2 = -L'c:/QtSDK/Desktop/Qt/4.7.4/mingw/lib' -lQtCore4
jt9.exe: $(OBJS2) libjt9.a
	$(CXX) -o jt9.exe $(OBJS2) $(LIBS2) libjt9.a ../libfftw3f_win.a \
	c:/MinGW/lib/libf95.a
	mkdir -p $(EXE_DIR)
	cp jt9.exe $(EXE_DIR)

OBJS3 = jt9sim.o
jt9sim.exe: $(OBJS3) libjt9.a
	$(FC) -o jt9sim.exe $(OBJS3) libjt9.a

OBJS4 = jt9code.o
jt9code.exe: $(OBJS4) libjt9.a
	$(FC) -o jt9code.exe $(OBJS4) libjt9.a

OBJS5 = jt65.o
jt65.exe: $(OBJS5) libjt9.a
	$(FC) -o jt65.exe $(OBJS5) libjt9.a ../libfftw3f_win.a

sync9.o: sync9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c sync9.f90

spec9.o: spec9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c spec9.f90

peakdt9.o: peakdt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c peakdt9.f90

jt9sim.o: jt9sim.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c jt9sim.f90

genjt9.o: genjt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c genjt9.f90

redsync.o: redsync.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c redsync.f90

unpackmsg.o: unpackmsg.f90
	$(FC) -c -O0 -fbounds-check -Wall -Wno-precision-loss unpackmsg.f90

ipcomm.o: ipcomm.cpp
	$(CXX) -c $(INCPATH) ipcomm.cpp

sec_midn.o: sec_midn.f90
	$(FC) -c -fno-second-underscore sec_midn.f90

#rig_control.o: rig_control.c
#	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include rig_control.c

tstrig.o: tstrig.c
	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include tstrig.c

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	rm -f *.o libjt9.a wsjtx.exe jt9sim.exe jt9.exe jt65.exe

license: mit


code:
# Set paths
EXE_DIR = ../../wsjtx_exp_install_latest
QT_DIR = /usr/include/qt5
INCPATH = -I${QT_DIR} -I${QT_DIR}/QtCore

CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
MKDIR = mkdir -p
CP = cp
RANLIB = ranlib
RM = rm -f

FFLAGS = -I/opt/local/include -O3 -funroll-loops -Wall -Wno-conversion -fno-second-underscore -DUNIX
CFLAGS = -I. -fbounds-check -fPIE

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

all:    libjt9.a jt65 jt65sim t3

OBJS1 = astrosub.o astro0.o astro.o sun.o coord.o tmoonsub.o \
        fmtmsg.o deg2grid.o\
        prog_args.o options.o pctile.o graycode.o sort.o chkmsg.o \
	igray.o fftw3mod.o packjt.o\
	four2a.o grid2deg.o wisdom.o \
	symspec.o analytic.o db.o \
	encode232.o interleave9.o\
	entail.o fano232.o gran.o sync9.o \
	fil3.o decoder.o fqso_first.o \
	twkfreq.o symspec2.o shell.o sync65.o peakup.o slope.o xcor.o\
	chkss2.o zplot9.o flat1.o flat2.o \
	symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat4.o determ.o baddata.o subtract65.o

libjt9.a: $(OBJS1)
	$(AR) libjt9.a $(OBJS1)
	$(RANLIB) libjt9.a

OBJS7 = jt65.o
jt65: $(OBJS7) libjt9.a libsfrsd.a
	$(FC) -o jt65 $(OBJS7) -L. -L/opt/local/lib -L./sfrsd2 -ljt9 -lsfrsd -lfftw3f_threads -lfftw3f
	$(CP) jt65 $(EXE_DIR)

OBJS2 = jt65sim.o wavhdr.o
jt65sim: $(OBJS2) libjt9.a
	$(FC) -o jt65sim $(OBJS2) -L. -L/opt/local/lib -ljt9 -lfftw3f
	$(CP) jt65sim $(EXE_DIR)

OBJS3 = t3.o
t3:   $(OBJS3) libjt9.a
	$(FC) -o t3 $(OBJS3) -L. -L/opt/local/lib -ljt9

OBJS6 = jt65code.o
jt65code: $(OBJS6) libjt9.a
	$(FC) -o jt65code $(OBJS6) libjt9.a

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	$(RM) *.o libjt9.a jt65


code:
# Set paths
EXE_DIR = C:/JTSDK/wsjtx_exp/install/Debug/bin
QT_DIR = /usr/include/qt5
INCPATH = -I${QT_DIR} -I${QT_DIR}/QtCore

CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
MKDIR = mkdir -p
CP = cp
RANLIB = ranlib
RM = rm -f

FFLAGS = -O2 -Wall -Wno-conversion -fno-second-underscore
CFLAGS = -I. -fbounds-check -fPIE -DBIGSYM

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

all:    libjt9.a jt65 jt65sim

OBJS1 = astrosub.o astro0.o astro.o sun.o coord.o tmoonsub.o \
        fmtmsg.o deg2grid.o\
        prog_args.o options.o pctile.o graycode.o sort.o chkmsg.o \
	igray.o fftw3mod.o packjt.o\
	four2a.o grid2deg.o wisdom.o \
	symspec.o analytic.o db.o \
	encode232.o interleave9.o \
	entail.o fano232.o gran.o sync9.o decjt9.o \
	fil3.o decoder.o timer.o exp_decode65.o fqso_first.o \
	twkfreq.o symspec2.o shell.o sync65.o peakup.o slope.o xcor.o\
	fillcom.o chkss2.o zplot9.o flat1.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat4.o determ.o baddata.o subtract65.o

libjt9.a: $(OBJS1)
	$(AR) libjt9.a $(OBJS1)
	$(RANLIB) libjt9.a

OBJS7 = jt65.o
jt65: $(OBJS7) libjt9.a libsfrsd.a
	$(FC) -o jt65 $(OBJS7) -L. -L./sfrsd2 -ljt9 -lsfrsd \
	 C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS2 = jt65sim.o wavhdr.o
jt65sim: $(OBJS2) libjt9.a
	$(FC) -o jt65sim $(OBJS2) -L. -ljt9 C:\JTSDK\fftw3f\libfftw3f-3.dll
	$(CP) jt65sim.exe $(EXE_DIR)

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	$(RM) *.o libjt9.a jt65


code:
# Set paths
EXE_DIR = ../../wsjtx_install

CC = clang
CXX = g++
FC = gfortran
AR = ar cr
MKDIR = mkdir -p
CP = cp
RANLIB = ranlib
RM = rm -f

FFLAGS = -I/opt/local/include -Wall -Wno-conversion -fno-second-underscore -DUNIX
CFLAGS = -I. -fbounds-check -fPIE

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

all:    libjt9.a jt65 jt65sim

OBJS1 = astrosub.o astro0.o astro.o sun.o coord.o tmoonsub.o \
        fmtmsg.o deg2grid.o\
        prog_args.o options.o pctile.o graycode.o sort.o chkmsg.o \
	igray.o fftw3mod.o packjt.o\
	four2a.o grid2deg.o wisdom.o \
	symspec.o analytic.o db.o \
	encode232.o interleave9.o\
	entail.o fano232.o sgran.o gran.o sync9.o decjt9.o \
	fil3.o decoder.o timer.o exp_decode65.o fqso_first.o\
	twkfreq.o symspec2.o shell.o sync65.o peakup.o slope.o xcor.o\
	fillcom.o chkss2.o zplot9.o flat1.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat4.o determ.o baddata.o subtract65.o

libjt9.a: $(OBJS1)
	$(AR) libjt9.a $(OBJS1)
	$(RANLIB) libjt9.a

OBJS7 = jt65.o
jt65: $(OBJS7) libjt9.a libftrsd.a
	$(FC) -o jt65 $(OBJS7) -L. -L/opt/local/lib -L./ftrsd -ljt9 -lftrsd -lfftw3f_threads -lfftw3f
	$(CP) jt65 $(EXE_DIR)

OBJS2 = jt65sim.o wavhdr.o
jt65sim:    $(OBJS2) libjt9.a
	$(FC) -o jt65sim $(OBJS2) -L. -L/opt/local/lib -ljt9 -lfftw3f
	$(CP) jt65sim $(EXE_DIR)

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	$(RM) *.o libjt9.a jt65


code:
# Set paths
EXE_DIR = ..\\..\\wsjtx_install
QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
FFTW3_DIR = ..

INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include

# Compilers
CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
RANLIB = ranlib
MKDIR = mkdir -p
CP = cp
RM = rm -f

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

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

all:	 jt9w

OBJS1 = jt9w.o smo.o sync9w.o pctile.o shell.o lorentzian.o fchisq0.o \
      softsym9w.o four2a.o interleave9.o jt9fano.o fano232.o packjt.o \
      deg2grid.o grid2deg.o fmtmsg.o db.o decode9w.o

jt9w: $(OBJS1)
	$(FC) -o jt9w $(OBJS1) -lfftw3f

OBJS2 = t1.o four2a.o db.o
t1:   $(OBJS2)
	$(FC) -o t1 $(OBJS2) -lfftw3f

OBJS3 = t2.o four2a.o db.o
t2:   $(OBJS3)
	$(FC) -o t2 $(OBJS3) -lfftw3f

OBJS4 = t3.o four2a.o db.o
t3:   $(OBJS4)
	$(FC) -o t3 $(OBJS4) -lfftw3f

.PHONY : clean

clean:
	$(RM) *.o JTMSKcode JTMSKcode.exe


code:
# Set paths
EXE_DIR = ..\\..\\wsjtx_install
QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
FFTW3_DIR = ..

INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include

# Compilers
CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
RANLIB = ranlib
MKDIR = mkdir -p
CP = cp
RM = rm -f

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

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

all:	 jt9w

OBJS1 = jt9w.o smo.o sync9w.o pctile.o shell.o lorentzian.o fchisq0.o \
      softsym9w.o four2a.o interleave9.o jt9fano.o fano232.o packjt.o \
      deg2grid.o grid2deg.o fmtmsg.o db.o decode9w.o

jt9w: $(OBJS1)
	$(FC) -o jt9w $(OBJS1) -lfftw3f

OBJS2 = t1.o four2a.o db.o
t1:   $(OBJS2)
	$(FC) -o t1 $(OBJS2) -lfftw3f

OBJS3 = t2.o four2a.o db.o
t2:   $(OBJS3)
	$(FC) -o t2 $(OBJS3) -lfftw3f

OBJS4 = t3.o
t3:   $(OBJS4)
	$(FC) -o t3 $(OBJS4) -L. -ljt9 C:\JTSDK\fftw3f\libfftw3f-3.dll

.PHONY : clean

clean:
	$(RM) *.o JTMSKcode JTMSKcode.exe


code:
# Makefile for MinGW on Windows
# Windows re-direct:
#   C> make > junk1 2>&1

# Set paths
EXE_DIR = ../../wsjtx_install

INCPATH = -I'C:/JTSDK/Qt5/5.2.1/mingw48_32/include/QtCore' \
	-I'C:/JTSDK/Qt5/5.2.1/mingw48_32/include/'

# Compilers
CC  = c:/JTSDK/Qt5/Tools/mingw48_32/bin/gcc
FC = c:/JTSDK/Qt5/Tools/mingw48_32/bin/gfortran
CXX = c:/JTSDK/Qt5/Tools/mingw48_32/bin/g++

FFLAGS = -O3 -Wall -Wno-conversion -fno-second-underscore -DWIN32
CFLAGS = -I. -DWIN32

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

all:    libjt9.a libastro.a jt9.exe jt9code.exe jt65code.exe jt9sim.exe

OBJS1 = pctile.o graycode.o sort.o chkmsg.o \
	unpackmsg.o igray.o unpackcall.o unpackgrid.o \
	grid2k.o unpacktext.o getpfx2.o packmsg.o deg2grid.o \
	packtext.o getpfx1.o packcall.o k2grid.o packgrid.o \
	nchar.o four2a.o grid2deg.o pfxdump.o wisdom.o \
	symspec.o analytic.o db.o genjt9.o flat1.o smo.o \
	packbits.o unpackbits.o encode232.o interleave9.o \
	entail.o fano232.o gran.o sync9.o jt9fano.o \
	fil3.o decoder.o grid2n.o n2grid.o timer.o \
	softsym.o getlags.o afc9.o fchisq.o twkfreq.o downsam9.o \
	peakdt9.o symspec2.o stdmsg.o morse.o azdist.o geodist.o \
	fillcom.o chkss2.o zplot9.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat4.o polfit.o determ.o baddata.o prog_args.o \
	options.o fmtmsg.o decjt9.o

libjt9.a: $(OBJS1)
	ar cr libjt9.a $(OBJS1)
	ranlib libjt9.a

OBJS2 = jt9.o jt9a.o jt9b.o jt9c.o ipcomm.o sec_midn.o usleep.o
LIBS2 = -L'C:/JTSDK/Qt5/5.2.1/mingw48_32/lib' -lQt5Core
jt9.exe: $(OBJS2) libjt9.a
	$(CXX) -o jt9.exe -static $(OBJS2) $(LIBS2) libjt9.a \
	C:\JTSDK\fftw3f\libfftw3f-3.dll -lgfortran
#	mkdir -p $(EXE_DIR)
	cp jt9.exe $(EXE_DIR)

OBJS3 = jt9sim.o
jt9sim.exe: $(OBJS3) libjt9.a
	$(FC) -o jt9sim.exe $(OBJS3) libjt9.a

OBJS4 = jt9code.o
jt9code.exe: $(OBJS4) libjt9.a
	$(FC) -o jt9code.exe $(OBJS4) libjt9.a
	cp jt9code.exe $(EXE_DIR)

OBJS5 = jt65.o
jt65.exe: $(OBJS5) libjt9.a
	$(FC) -o jt65.exe $(OBJS5) libjt9.a ../libfftw3f_win.a

OBJS7 = astrosub.o astro0.o astro.o tm2.o grid2deg.o sun.o moondop.o \
	coord.o dot.o moon2.o tmoonsub.o toxyz.o geocentric.o \
	dcoord.o

libastro.a: $(OBJS7)
	ar cr libastro.a $(OBJS7)
	ranlib libastro.a

OBJS6 = jt65code.o
jt65code.exe: $(OBJS6) libjt9.a
	$(FC) -o jt65code.exe $(OBJS6) libjt9.a
	cp jt65code.exe $(EXE_DIR)

sync9.o: sync9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c sync9.f90

spec9.o: spec9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c spec9.f90

peakdt9.o: peakdt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c peakdt9.f90

jt9sim.o: jt9sim.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c jt9sim.f90

genjt9.o: genjt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c genjt9.f90

redsync.o: redsync.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c redsync.f90

unpackmsg.o: unpackmsg.f90
	$(FC) -c -O0 -fbounds-check -Wall -Wno-precision-loss unpackmsg.f90

ipcomm.o: ipcomm.cpp
	$(CXX) -c $(INCPATH) ipcomm.cpp

sec_midn.o: sec_midn.f90
	$(FC) -c -fno-second-underscore sec_midn.f90

#rig_control.o: rig_control.c
#	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include rig_control.c

tstrig.o: tstrig.c
	$(CC) -c -Wall -I..\..\..\hamlib-1.2.15.3\include tstrig.c

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	rm -f *.o libjt9.a wsjtx.exe jt9sim.exe jt9.exe jt65.exe


code:
# Set paths
EXE_DIR = ../../wsjtx_install
QT_DIR = /usr/include/qt5
INCPATH = -I${QT_DIR} -I${QT_DIR}/QtCore

CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
MKDIR = mkdir -p
CP = cp
RANLIB = ranlib
RM = rm -f

FFLAGS = -O3 -funroll-loops -Wall -Wno-conversion -fno-second-underscore -DUNIX
CFLAGS = -I. -fbounds-check -fPIE

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

all:    libjt9.a jt9sim jt9 jt9code jt65code

OBJS1 = astrosub.o astro0.o astro.o tm2.o sun.o moondop.o coord.o tmoonsub.o \
        fmtmsg.o geocentric.o moon2.o toxyz.o dot.o dcoord.o \
        prog_args.o options.o pctile.o graycode.o sort.o chkmsg.o \
	unpackmsg.o igray.o unpackcall.o unpackgrid.o \
	grid2k.o unpacktext.o getpfx2.o packmsg.o deg2grid.o \
	packtext.o getpfx1.o packcall.o k2grid.o packgrid.o \
	nchar.o four2a.o grid2deg.o pfxdump.o wisdom.o \
	symspec.o analytic.o db.o genjt9.o jt9fano.o \
	packbits.o unpackbits.o encode232.o interleave9.o \
	entail.o fano232.o gran.o sync9.o decjt9.o \
	fil3.o decoder.o grid2n.o n2grid.o timer.o \
	softsym.o peakdt9.o getlags.o afc9.o fchisq.o \
	twkfreq.o downsam9.o symspec2.o ipcomm.o sleep_msec.o \
	stdmsg.o sec_midn.o usleep.o azdist.o geodist.o morse.o \
	fillcom.o chkss2.o zplot9.o flat1.o flat2.o \
	jt65a.o symspec65.o flat65.o ccf65.o decode65a.o \
	filbig.o fil6521.o afc65b.o decode65b.o setup65.o \
	extract.o fchisq65.o demod64a.o chkhist.o interleave63.o ccf2.o \
	move.o indexx.o graycode65.o twkfreq65.o smo.o smo121.o \
	wrapkarn.o init_rs.o encode_rs.o decode_rs.o gen65.o fil4.o \
	flat4.o polfit.o determ.o baddata.o

libjt9.a: $(OBJS1)
	$(AR) libjt9.a $(OBJS1)
	$(RANLIB) libjt9.a

OBJS2 = jt9.o jt9a.o jt9b.o jt9c.o

jt9: 	$(OBJS2) libjt9.a
	$(CXX) -o jt9 $(OBJS2) -L. -ljt9 -lQt5Core -lfftw3f_threads -lfftw3f `$(FC) -print-file-name=libgfortran.so`
	$(MKDIR) -p $(EXE_DIR)
	$(CP) jt9 $(EXE_DIR)

OBJS3 = jt9sim.o
jt9sim: $(OBJS3) libjt9.a
	$(FC) -o jt9sim $(OBJS3) -L. -ljt9

OBJS4 = jt9code.o
jt9code: $(OBJS4) libjt9.a
	$(FC) -o jt9code $(OBJS4) -L. -ljt9
	$(CP) jt9code $(EXE_DIR)

OBJS6 = jt65code.o
jt65code: $(OBJS6) libjt9.a
	$(FC) -o jt65code $(OBJS6) libjt9.a
	$(CP) jt65code $(EXE_DIR)

sync9.o: sync9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c sync9.f90

peakdf9.o: peakdf9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c peakdf9.f90

peakdt9.o: peakdt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c peakdt9.f90

jt9sim.o: jt9sim.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c jt9sim.f90

genjt9.o: genjt9.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c genjt9.f90

redsync.o: redsync.f90 jt9sync.f90
	$(FC) $(FFLAGS) -c redsync.f90

ipcomm.o: ipcomm.cpp
	$(CXX) -c $(INCPATH) -fPIE ipcomm.cpp

sec_midn.o: sec_midn.f90
	$(FC) -c -fno-second-underscore sec_midn.f90

init_rs.o: init_rs.c
	$(CC) -c -DBIGSYM=1 -o init_rs.o init_rs.c

encode_rs.o: encode_rs.c
	$(CC) -c -DBIGSYM=1 -o encode_rs.o encode_rs.c

decode_rs.o: decode_rs.c
	$(CC) -c -DBIGSYM=1 -o decode_rs.o decode_rs.c

.PHONY : clean

clean:
	$(RM) *.o libjt9.a wsjtx jt9sim jt9 jt9code


code:
# Set paths
EXE_DIR = ..\\..\\wsjtx_install
QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
FFTW3_DIR = ..

INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include

# Compilers
CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
RANLIB = ranlib
MKDIR = mkdir -p
CP = cp
RM = rm -f

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

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

#all:    jt9code JTMSKcode
all:	 jtmsk JTMSKcode

OBJS4 = jt9code.o packjt.o fmtmsg.o gen9.o deg2grid.o grid2deg.o \
      entail.o encode232.o interleave9.o graycode.o igray.o
jt9code: $(OBJS4)
	$(FC) -o jt9code $(OBJS4)

OBJS5 = JTMSKcode.o packjt.o fmtmsg.o genmsk.o deg2grid.o grid2deg.o \
      entail.o nhash.o tab.o vit213.o
JTMSKcode: $(OBJS5)
	$(FC) -o JTMSKcode $(OBJS5)

OBJS6 = jtmsk.o jtmsk.o analytic.o four2a.o db.o mskdf.o pctile.o \
      sort.o tweak1.o syncmsk.o genmsk.o packjt.o fmtmsg.o indexx.o \
      deg2grid.o grid2deg.o entail.o nhash.o tab.o vit213.o
jtmsk: $(OBJS6)
	$(FC) -o jtmsk $(OBJS6) -lfftw3f

OBJS1 = t1.o four2a.o db.o
t1:   $(OBJS1)
	$(FC) -o t1 $(OBJS1) -lfftw3f

OBJS2 = t2.o four2a.o db.o
t2:   $(OBJS2)
	$(FC) -o t2 $(OBJS2) -lfftw3f

OBJS6 = t6.o four2a.o db.o
t6:   $(OBJS6)
	$(FC) -o t6 $(OBJS6) -lfftw3f

nhash.o: wsprd/nhash.h wsprd/nhash.c
	 $(CC) -c -O2 wsprd/nhash.c

.PHONY : clean

clean:
	$(RM) *.o JTMSKcode JTMSKcode.exe


code:
# Set paths
EXE_DIR = ..\\..\\wsjtx_install
QT_DIR = C:/wsjt-env/Qt5/5.2.1/mingw48_32
FFTW3_DIR = ..

INCPATH = -I${QT_DIR}/include/QtCore -I${QT_DIR}/include

# Compilers
CC = gcc
CXX = g++
FC = gfortran
AR = ar cr
RANLIB = ranlib
MKDIR = mkdir -p
CP = cp
RM = rm -f

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

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

#all:    jt9code JTMSKcode.exe
all:	 jtmsk.exe JTMSKsim.exe JTMSKcode.exe fixwav.exe

OBJS3 = JTMSKsim.o wavhdr.o gran.o four2a.o db.o
JTMSKsim.exe: $(OBJS3)
	$(FC) -o JTMSKsim.exe $(OBJS3) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS4 = jt9code.o packjt.o fmtmsg.o gen9.o deg2grid.o grid2deg.o \
      entail.o encode232.o interleave9.o graycode.o igray.o
jt9code: $(OBJS4)
	$(FC) -o jt9code $(OBJS4)

OBJS5 = JTMSKcode.o packjt.o fmtmsg.o genmsk.o deg2grid.o grid2deg.o \
      entail.o tab.o vit213.o hashing.o nhash.o
JTMSKcode.exe: $(OBJS5)
	$(FC) -o JTMSKcode.exe $(OBJS5)

OBJS6 = jtmsk.o analytic.o four2a.o db.o pctile.o \
      shell.o tweak1.o syncmsk.o genmsk.o packjt.o fmtmsg.o indexx.o \
      deg2grid.o grid2deg.o entail.o hashing.o nhash.o tab.o vit213.o \
      mskdt.o rectify_msk.o timer.o jtmsk_decode.o genmsk_short.o \
      jtmsk_short.o golay24_table.o hash.o

jtmsk.exe: $(OBJS6)
	$(FC) -o jtmsk.exe $(OBJS6) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS1 = fixwav.o wavhdr.o
fixwav.exe: $(OBJS1)
	$(FC) -o fixwav.exe $(OBJS1)

OBJS2 = t2.o four2a.o db.o
t2:   $(OBJS2)
	$(FC) -o t2 $(OBJS2) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS6 = t6.o four2a.o db.o
t6:   $(OBJS6)
	$(FC) -o t6 $(OBJS6) C:\JTSDK\fftw3f\libfftw3f-3.dll

.PHONY : clean

clean:
	$(RM) *.o JTMSKcode JTMSKcode.exe


code:
subroutine addit(itone,nfsample,nsym,nsps,ifreq,sig,dat)

  integer itone(nsym)
  real dat(60*12000)
  real*8 f,dt,twopi,phi,dphi,baud,fsample,freq,tsym,t

  tsym=nsps*1.d0/nfsample            !Symbol duration
  baud=1.d0/tsym
  fsample=12000.d0                   !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  twopi=8.d0*atan(1.d0)
  dphi=0.

  iters=1
  if(nsym.eq.79) iters=2
  do iter=1,iters
     f=ifreq
     phi=0.
     ntot=nsym*tsym/dt
     k=12000                             !Start audio at t = 1.0 s
     t=0.
     if(nsym.eq.79) k=12000 + (iter-1)*12000*30  !Special case for FT8
     isym0=-1
     do i=1,ntot
        t=t+dt
        isym=nint(t/tsym) + 1
        if(isym.ne.isym0) then
           freq=f + itone(isym)*baud
           dphi=twopi*freq*dt
           isym0=isym
        endif
        phi=phi + dphi
        if(phi.gt.twopi) phi=phi-twopi
        xphi=phi
        k=k+1
        dat(k)=dat(k) + sig*sin(xphi)
     enddo
  enddo

  return
end subroutine addit

subroutine addcw(icw,ncw,ifreq,sig,dat)

  integer icw(ncw)
  real dat(60*12000)
  real s(60*12000)
  real x(60*12000)
  real y(60*12000)
  real*8 dt,twopi,phi,dphi,fsample,tdit,t

  wpm=25.0
  nspd=nint(1.2*12000.0/wpm)
  fsample=12000.d0                   !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  tdit=nspd*dt
  twopi=8.d0*atan(1.d0)
  dphi=twopi*ifreq*dt
  phi=0.
  k=12000                             !Start audio at t = 1.0 s
  t=0.
  npts=60*12000
  x=0.
  do i=1,npts
     t=t+dt
     j=nint(t/tdit) + 1
     j=mod(j-1,ncw) + 1
     phi=phi + dphi
     if(phi.gt.twopi) phi=phi-twopi
     xphi=phi
     k=k+1
     x(k)=icw(j)
     s(k)=sin(xphi)
     if(t.ge.59.5) exit
  enddo

  nadd=0.004/dt
  call smo(x,npts,y,nadd)
  y=y/nadd
  dat=dat + sig*y*s

  return
end subroutine addcw


code:
subroutine afc65b(cx,npts,fsample,nflip,mode65,a,ccfbest,dtbest)

! Find delta f, f1, f2 ==> a(1:3)

  complex cx(npts)
  real a(5),deltaa(5)

  a=0.
  i2=8*mode65
  i1=-i2
  j2=8*mode65
  j1=-j2
  ccfmax=0.
  istep=2*mode65
  do iter=1,2
     do i=i1,i2,istep
        a(1)=i
        do j=j1,j2,istep
           a(2)=j
           chisq=fchisq65(cx,npts,fsample,nflip,a,ccf,dtmax)
           if(ccf.gt.ccfmax) then
              a1=a(1)
              a2=a(2)
              ccfmax=ccf
           endif
!           write(81,3081) istep,i1,i2,j1,j2,i,j,ccf,ccfmax,dtmax,a1,a2
!3081       format(7i4,5f8.2)
        enddo
     enddo
     i1=a1-istep
     i2=a1+istep
     j1=a2-istep
     j2=a2+istep
     istep=1
  enddo

!  a(1)=0.
!  a(2)=0.
  a(1)=a1
  a(2)=a2
  a(3)=0.
  a(4)=0.
  deltaa(1)=2.0*mode65
  deltaa(2)=2.0*mode65
  deltaa(3)=1.0
  nterms=2                                  !Maybe 2 is enough?

!  Start the iteration
  chisqr=0.
  chisqr0=1.e6
  do iter=1,100                              !How many iters is enough?
     do j=1,nterms
        chisq1=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
        fn=0.
        delta=deltaa(j)
10      a(j)=a(j)+delta
        chisq2=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
        if(chisq2.eq.chisq1) go to 10
        if(chisq2.gt.chisq1) then
           delta=-delta                      !Reverse direction
           a(j)=a(j)+delta
           tmp=chisq1
           chisq1=chisq2
           chisq2=tmp
        endif
20      fn=fn+1.0
        a(j)=a(j)+delta
        chisq3=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
        if(chisq3.lt.chisq2) then
           chisq1=chisq2
           chisq2=chisq3
           go to 20
        endif

! Find minimum of parabola defined by last three points
        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
        a(j)=a(j)-delta
        deltaa(j)=deltaa(j)*fn/3.
!        write(*,4000) iter,j,a(1:2),-chisq2
!4000    format(2i2,4f9.4)
     enddo
     chisqr=fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)
     fdiff=chisqr/chisqr0-1.0
!     write(*,4000) 0,0,a(1:2),-chisqr,fdiff
     if(abs(fdiff).lt.0.0001) exit
     chisqr0=chisqr
  enddo
  ccfbest=ccfmax * (1378.125/fsample)**2
  dtbest=dtmax

  return
end subroutine afc65b


code:
subroutine afc9(c3a,npts,fsample,a,syncpk)

  parameter (NZ2=1512)
  complex c3a(0:NZ2-1)
  complex c3(0:NZ2-1)
  real a(3),deltaa(3)

  a(1)=0.                                   !f0
  a(2)=0.                                   !f1
  a(3)=0.                                   !f2
  deltaa(1)=1.736
  deltaa(2)=1.736
  deltaa(3)=1.0
  nterms=3

! Start the iteration
  chisqr=0.
  chisqr0=1.e6
  c3=c3a
  a3=a(3)
  do iter=1,4
     do j=1,nterms
        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
        chisq1=fchisq(c3,npts,fsample,a)
        fn=0.
        delta=deltaa(j)
10      a(j)=a(j)+delta
        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
        chisq2=fchisq(c3,npts,fsample,a)
        if(chisq2.eq.chisq1) go to 10
        if(chisq2.gt.chisq1) then
           delta=-delta                      !Reverse direction
           a(j)=a(j)+delta
           tmp=chisq1
           chisq1=chisq2
           chisq2=tmp
        endif
20      fn=fn+1.0
        a(j)=a(j)+delta
        if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
        chisq3=fchisq(c3,npts,fsample,a)
        if(chisq3.lt.chisq2) then
           chisq1=chisq2
           chisq2=chisq3
           go to 20
        endif

! Find minimum of parabola defined by last three points
        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
        a(j)=a(j)-delta
        if(j.lt.3) deltaa(j)=deltaa(j)*fn/3.
!        write(*,4000) iter,j,a,-chisq2
!4000    format(i1,i2,3f10.4,f11.3)
     enddo
     if(a(3).ne.a3) call shft(c3a,a(3),a3,c3)
     chisqr=fchisq(c3,npts,fsample,a)
!     write(*,4000) 0,0,a,-chisqr
     if(chisqr/chisqr0.gt.0.99) exit
     chisqr0=chisqr
  enddo

  syncpk=-chisqr
  c3a=c3
!  write(*,4001) a,syncpk
!4001 format(3x,3f10.4,f11.3)

  return
end subroutine afc9

subroutine shft(c3a,a3a,a3,c3)
  complex c3a(0:1359)
  complex  c3(0:1359)

  a3=a3a
  n=nint(a3)
  c3=cshift(c3a,n)
  if(n.gt.0) c3(1360-n:)=0.0
  if(n.lt.0) c3(:n-1)=0.0

  return
end subroutine shft


code:
program allsim

! Generate simulated data for WSJT-X slow modes: JT4, JT9, JT65, QRA64,
! and WSPR.  Also unmodulated carrier and 20 WPM CW.

  use wavhdr
  use packjt
  parameter (NMAX=60*12000)
  type(hdr) h
  integer*2 iwave(NMAX)                  !Generated waveform (no noise)
  integer itone(206)                     !Channel symbols (values 0-8)
  integer icw(250)
  integer*1 msgbits(87)
  logical*1 bcontest
  real*4 dat(NMAX)
  character message*22,msgsent*22,arg*8,mygrid*6

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage:   allsim <snr>'
     go to 999
  endif

  call getarg(1,arg)
  read(arg,*) snrdb                  !S/N in dB (2500 hz reference BW)

  message='CQ KA2ABC FN20'
  mygrid='FN20  '
  bcontest=.false.
  rmsdb=25.
  rms=10.0**(0.05*rmsdb)
  sig=10.0**(0.05*snrdb)
  npts=NMAX

  call init_random_seed()       !Seed Fortran RANDOM_NUMBER generator
  call sgran()                  !Seed C rand generator (used in gran)

  h=default_header(12000,npts)
  open(10,file='000000_0000.wav',access='stream',status='unknown')
  do i=1,npts                   !Generate gaussian noise
     dat(i)=gran()
  enddo

  itone=0
  call addit(itone,12000,85,6912,400,sig,dat)    !Unmodulated carrier

  call morse('CQ CQ DE KA2ABC KA2ABC',icw,ncw)
!  print*,ncw
!  write(*,3001) icw(1:ncw)
!3001 format(50i1)
  call addcw(icw,ncw,600,sig,dat)                !CW

  call genwspr(message,msgsent,itone)
  call addit(itone,12000,86,8192,800,sig,dat)    !WSPR (only 59 s of data)

  call gen9(message,0,msgsent,itone,itype)
  call addit(itone,12000,85,6912,1000,sig,dat)   !JT9

  call gen4(message,0,msgsent,itone,itype)
  call addit(itone,11025,206,2520,1200,sig,dat)  !JT4

  i3bit=0                                        ! ### TEMPORARY ??? ###
  call genft8(message,mygrid,bcontest,i3bit,msgsent,msgbits,itone)
  call addit(itone,12000,79,1920,1400,sig,dat)   !FT8

  call genqra64(message,0,msgsent,itone,itype)
  call addit(itone,12000,84,6912,1600,sig,dat)   !QRA64

  call gen65(message,0,msgsent,itone,itype)
  call addit(itone,11025,126,4096,1800,sig,dat)  !JT65

  iwave(1:npts)=nint(rms*dat(1:npts))

  write(10) h,iwave(1:npts)
  close(10)

999 end program allsim


code:
subroutine ana64(dd,npts,c0)

  use timer_module, only: timer

  parameter (NMAX=60*12000)                  !Max size of raw data at 12000 Hz
  parameter (NSPS=3456)                      !Samples per symbol at 6000 Hz
  parameter (NSPC=7*NSPS)                    !Samples per Costas array
  real dd(NMAX)                              !Raw data
  complex c0(0:720000)                       !Complex spectrum of dd()
  save

  nfft1=672000
  nfft2=nfft1/2
  df1=12000.0/nfft1
  fac=2.0/nfft1
  c0(0:npts-1)=fac*dd(1:npts)
  c0(npts:nfft1)=0.
  call four2a(c0,nfft1,1,-1,1)             !Forward c2c FFT
  c0(nfft2/2+1:nfft2)=0.
  c0(0)=0.5*c0(0)
  call four2a(c0,nfft2,1,1,1)              !Inverse c2c FFT; c0 is analytic sig

  return
end subroutine ana64


code:
subroutine ana932(dat,npts0,cdat,npts)

  real dat(npts0)
  complex cdat(262145)

  n=log(float(npts0))/log(2.0)
  nfft1=2**(n+1)
  nfft2=9*nfft1/32
  df932=11025.0/nfft1
  fac=2.0/nfft1
  do i=1,npts0/2
     cdat(i)=fac*cmplx(dat(2*i-1),dat(2*i))
  enddo
  cdat(npts0/2+1:nfft1/2)=0.
  call four2a(cdat,nfft1,1,-1,0)               !Forward r2c FFT
  call four2a(cdat,nfft2,1,1,1)                !Inverse c2c FFT
  npts=npts0*9.0/32.0                          !Downsampled data length
  npts2=npts

  return
end subroutine ana932


code:
subroutine analytic(d,npts,nfft,c,pc,beq)

! Convert real data to analytic signal

  parameter (NFFTMAX=1024*1024)

  real d(npts)              ! passband signal
  real h(NFFTMAX/2)         ! real BPF magnitude
  real*8 pc(5),pclast(5)    ! static phase coeffs
  real*8 ac(5),aclast(5)    ! amp coeffs
  real*8 fp

  complex corr(NFFTMAX/2)  ! complex frequency-dependent correction
  complex c(NFFTMAX)        ! analytic signal

  logical*1 beq            ! boolean static equalizer flag

  data nfft0/0/
  data aclast/0.0,0.0,0.0,0.0,0.0/
  data pclast/0.0,0.0,0.0,0.0,0.0/
!  data ac/1.0,0.05532,0.11438,0.12918,0.09274/ ! amp coeffs for TS2000
  data ac/1.0,0.0,0.0,0.0,0.0/

  save corr,nfft0,h,ac,aclast,pclast,pi,t,beta

  df=12000.0/nfft
  nh=nfft/2
  if( nfft.ne.nfft0 ) then
     pi=4.0*atan(1.0)
     t=1.0/2000.0
     beta=0.1
     do i=1,nh+1
        ff=(i-1)*df
        f=ff-1500.0
        h(i)=1.0
        if(abs(f).gt.(1-beta)/(2*t) .and. abs(f).le.(1+beta)/(2*t)) then
           h(i)=h(i)*0.5*(1+cos((pi*t/beta )*(abs(f)-(1-beta)/(2*t))))
        elseif( abs(f) .gt. (1+beta)/(2*t) ) then
           h(i)=0.0
        endif
     enddo
     nfft0=nfft
  endif

  if( any(aclast .ne. ac) .or. any(pclast .ne. pc) ) then
     aclast=ac
     pclast=pc
!     write(*,3001) pc
!3001 format('Phase coeffs:',5f12.6)
     do i=1,nh+1
        ff=(i-1)*df
        f=ff-1500.0
        fp=f/1000.0
        corr(i)=ac(1)+fp*(ac(2)+fp*(ac(3)+fp*(ac(4)+fp*ac(5))))
        pd=fp*fp*(pc(3)+fp*(pc(4)+fp*pc(5))) ! ignore 1st two terms
        corr(i)=corr(i)*cmplx(cos(pd),sin(pd))
     enddo
  endif

  fac=2.0/nfft
  c(1:npts)=fac*d(1:npts)
  c(npts+1:nfft)=0.
  call four2a(c,nfft,1,-1,1)               !Forward c2c FFT

  if( beq ) then
    c(1:nh+1)=h(1:nh+1)*corr(1:nh+1)*c(1:nh+1)
  else
    c(1:nh+1)=h(1:nh+1)*c(1:nh+1)
  endif

  c(1)=0.5*c(1)                            !Half of DC term
  c(nh+2:nfft)=0.                          !Zero the negative frequencies
  call four2a(c,nfft,1,1,1)                !Inverse c2c FFT
  return
end subroutine analytic


code:
subroutine astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid,              &
     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,    &
     dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0,    &
     width1,width2,xlst8,techo8)

  parameter (DEGS=57.2957795130823d0)
  character*6 mygrid,hisgrid
  real*8 AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8
  real*8 dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,xnr8,dfdt,dfdt0,dt
  real*8 sd8,poloffset8,width1,width2,xlst8
  real*8 uth8,techo8,freq8
  real*8 xl,b
  common/librcom/xl(2),b(2)
  data uth8z/0.d0/
  save

  uth=uth8
  call astro(nyear,month,nday,uth,freq8,hisgrid,2,1,                 &
       AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler,            &
       dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr,               &
       day,xlon2,xlat2,xlst,techo)
  AzMoonB8=AzMoon
  ElMoonB8=ElMoon
  xl2=xl(1)
  xl2a=xl(2)
  b2=b(1)
  b2a=b(2)
  call astro(nyear,month,nday,uth,freq8,mygrid,1,1,                  &
       AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler,            &
       dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr,               &
       day,xlon1,xlat1,xlst,techo)
  xl1=xl(1)
  xl1a=xl(2)
  b1=b(1)
  b1a=b(2)
  techo8=techo

  fghz=1.d-9*freq8
  dldt1=DEGS*(xl1a-xl1)
  dbdt1=DEGS*(b1a-b1)
  dldt2=DEGS*(xl2a-xl2)
  dbdt2=DEGS*(b2a-b2)
  rate1=2.0*sqrt(dldt1**2 + dbdt1**2)
  width1=0.5*6741*fghz*rate1
  rate2=sqrt((dldt1+dldt2)**2 + (dbdt1+dbdt2)**2)
  width2=0.5*6741*fghz*rate2

  AzSun8=AzSun
  ElSun8=ElSun
  AzMoon8=AzMoon
  ElMoon8=ElMoon
  dbMoon8=dbMoon
  RAMoon8=RAMoon/15.0
  DecMoon8=DecMoon
  HA8=HA
  xlst8=xlst
  Dgrd8=Dgrd
  sd8=sd
  poloffset8=poloffset
  xnr8=xnr
  ndop=nint(doppler)
  ndop00=nint(doppler00)

  if(uth8z.eq.0.d0) then
     uth8z=uth8-1.d0/3600.d0
     dopplerz=doppler
     doppler00z=doppler00
  endif

  dt=60.0*(uth8-uth8z)
  if(dt.le.0) dt=1.d0/60.d0
  dfdt=(doppler-dopplerz)/dt
  dfdt0=(doppler00-doppler00z)/dt
  uth8z=uth8
  dopplerz=doppler
  doppler00z=doppler00

  return
end subroutine astro0


code:
subroutine astrosub(nyear,month,nday,uth8,freq8,mygrid,hisgrid,          &
     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,  &
     RAMoon8,DecMoon8,Dgrd8,poloffset8,xnr8,techo8,width1,width2,bTx,    &
     AzElFileName,jpleph)

  implicit real*8 (a-h,o-z)
  character*6 mygrid,hisgrid,c1*1
  character*6 AzElFileName*(*),jpleph*(*)
  character*256 jpleph_file_name
  logical*1 bTx
  common/jplcom/jpleph_file_name

  jpleph_file_name=jpleph

  call astro0(nyear,month,nday,uth8,freq8,mygrid,hisgrid,                &
     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,  &
     dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0,  &
     width1,width2,xlst8,techo8)

  if (len_trim(AzElFileName) .eq. 0) go to 999
  imin=60*uth8
  isec=3600*uth8
  ih=uth8
  im=mod(imin,60)
  is=mod(isec,60)
  open(15,file=AzElFileName,status='unknown',err=900)
  c1='R'
  nRx=1
  if(bTx) then
     c1='T'
     nRx=0
  endif
  AzAux=0.
  ElAux=0.
  nfreq=freq8/1000000
  doppler=ndop
  doppler00=ndop00
  write(15,1010,err=10) ih,im,is,AzMoon8,ElMoon8,                     &
       ih,im,is,AzSun8,ElSun8,                                        &
       ih,im,is,AzAux,ElAux,                                          &
       nfreq,doppler,dfdt,doppler00,dfdt0,c1
!       TXFirst,TRPeriod,poloffset,Dgrd,xnr,ave,rms,nRx
1010 format(                                                          &
      i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Moon'/               &
      i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Sun'/                &
      i2.2,':',i2.2,':',i2.2,',',f5.1,',',f5.1,',Source'/             &
      i5,',',f8.1,',',f8.2,',',f8.1,',',f8.2,',Doppler, ',a1)
!      i1,',',i3,',',f8.1,','f8.1,',',f8.1,',',f12.3,',',f12.3,',',i1,',RPol')
10 close(15)
  go to 999

900 print*,'Error opening azel.dat'

999 return
end subroutine astrosub


code:
subroutine avecho(id2,ndop,nfrit,nqual,f1,xlevel,sigdb,snr,dfreq,width)

  integer TXLENGTH
  parameter (TXLENGTH=27648)           !27*1024
  parameter (NFFT=32768,NH=NFFT/2)
  parameter (NZ=4096)
  integer*2 id2(34560)                 !Buffer for Rx data
  real sa(NZ)      !Avg spectrum relative to initial Doppler echo freq
  real sb(NZ)      !Avg spectrum with Dither and changing Doppler removed
  integer nsum       !Number of integrations
  real dop0          !Doppler shift for initial integration (Hz)
  real dop           !Doppler shift for current integration (Hz)
  real s(8192)
  real x(NFFT)
  integer ipkv(1)
  complex c(0:NH)
  equivalence (x,c),(ipk,ipkv)
  common/echocom/nclearave,nsum,blue(NZ),red(NZ)
  save dop0,sa,sb

  dop=ndop
  sq=0.
  do i=1,TXLENGTH
     x(i)=id2(i)
     sq=sq + x(i)*x(i)
  enddo
  xlevel=10.0*log10(sq/TXLENGTH)

  if(nclearave.ne.0) nsum=0
  if(nsum.eq.0) then
     dop0=dop                             !Remember the initial Doppler
     sa=0.                                !Clear the average arrays
     sb=0.
  endif

  x(TXLENGTH+1:)=0.
  x=x/TXLENGTH
  call four2a(x,NFFT,1,-1,0)
  df=12000.0/NFFT
  do i=1,8192                             !Get spectrum 0 - 3 kHz
     s(i)=real(c(i))**2 + aimag(c(i))**2
  enddo

  fnominal=1500.0           !Nominal audio frequency w/o doppler or dither
  ia=nint((fnominal+dop0-nfrit)/df)
  ib=nint((f1+dop-nfrit)/df)
  if(ia.lt.600 .or. ib.lt.600) go to 900
  if(ia.gt.7590 .or. ib.gt.7590) go to 900

  nsum=nsum+1

  do i=1,NZ
     sa(i)=sa(i) + s(ia+i-2048)    !Center at initial doppler freq
     sb(i)=sb(i) + s(ib+i-2048)    !Center at expected echo freq
  enddo

  call pctile(sb,200,50,r0)
  call pctile(sb(1800),200,50,r1)

  sum=0.
  sq=0.
  do i=1,NZ
     y=r0 + (r1-r0)*(i-100.0)/1800.0
     blue(i)=sa(i)/y
     red(i)=sb(i)/y
     if(i.le.500 .or. i.ge.3597) then
        sum=sum+red(i)
        sq=sq + (red(i)-1.0)**2
     endif
  enddo
  ave=sum/1000.0
  rms=sqrt(sq/1000.0)

  redmax=maxval(red)
  ipkv=maxloc(red)
  fac=10.0/max(redmax,10.0)
  dfreq=(ipk-2048)*df
  snr=(redmax-ave)/rms

  sigdb=-99.0
  if(ave.gt.0.0) sigdb=10.0*log10(redmax/ave - 1.0) - 35.7

  nqual=0
  if(nsum.ge.2 .and. nsum.lt.4)  nqual=(snr-4)/5
  if(nsum.ge.4 .and. nsum.lt.8)  nqual=(snr-3)/4
  if(nsum.ge.8 .and. nsum.lt.12) nqual=(snr-3)/3
  if(nsum.ge.12) nqual=(snr-2.5)/2.5
  if(nqual.lt.0)  nqual=0
  if(nqual.gt.10) nqual=10

! Scale for plotting
  blue=fac*blue
  red=fac*red

  sum=0.
  do i=ipk,ipk+300
     if(i.gt.NZ) exit
     if(red(i).lt.1.0) exit
     sum=sum+(red(i)-1.0)
  enddo
  do i=ipk-1,ipk-300,-1
     if(i.lt.1) exit
     if(red(i).lt.1.0) exit
     sum=sum+(red(i)-1.0)
  enddo
  bins=sum/(red(ipk)-1.0)
  width=df*bins
  nsmo=max(0.0,0.25*bins)

  do i=1,nsmo
     call smo121(red,NZ)
     call smo121(blue,NZ)
  enddo

900  return
end subroutine avecho


code:
subroutine averms(x,n,nskip,ave,rms)
  real x(n)
  integer ipk(1)

  ns=0
  s=0.
  sq=0.
  ipk=maxloc(x)
  do i=1,n
     if(abs(i-ipk(1)).gt.nskip) then
        s=s + x(i)
        sq=sq + x(i)**2
        ns=ns+1
     endif
  enddo
  ave=s/ns
  rms=sqrt(sq/ns - ave*ave)

  return
end subroutine averms


code: ! The contents of this file have been migrated to lib/jt4_decode.f90


code:
subroutine azdist(grid1,grid2,utch,nAz,nEl,nDmiles,nDkm,nHotAz,nHotABetter)

  character*(*) grid1,grid2
  character*6 MyGrid,HisGrid,mygrid0,hisgrid0
  real*8 utch,utch0
  logical HotABetter,IamEast
  real eltab(22),daztab(22)
  data eltab/18.,15.,13.,11.,9.,8.,7.,6.,5.3,4.7,4.,3.3,2.7,    &
       2.,1.5,1.,0.8,0.6,0.4,0.2,0.0,0.0/
  data daztab/21.,18.,16.,15.,14.,13.,12.,11.,10.7,10.3,10.,    &
       10.,10.,10.,10.,10.,10.,9.,9.,9.,8.,8./
  data mygrid0/"      "/,hisgrid0/"      "/,utch0/-999.d0/
  save

  MyGrid=grid1//'  '
  HisGrid=grid2//'  '
  if(ichar(MyGrid(5:5)).eq.0)  MyGrid(5:6)='  '
  if(ichar(HisGrid(5:5)).eq.0) HisGrid(5:6)='  '

  if(MyGrid.eq.HisGrid) then
     naz=0
     nel=0
     ndmiles=0
     ndkm=0
     nhotaz=0
     nhotabetter=1
     go to 999
  endif

  if(mygrid.eq.mygrid0 .and. hisgrid.eq.hisgrid0 .and.              &
       abs(utch-utch0).lt.0.1666667d0) go to 900
  utch0=utch
  mygrid0=mygrid
  hisgrid0=hisgrid
  utchours=utch

  if(MyGrid(5:5).eq.' ') MyGrid(5:5)='m'
  if(MyGrid(6:6).eq.' ') MyGrid(6:6)='m'
  if(HisGrid(5:5).eq.' ') HisGrid(5:5)='m'
  if(HisGrid(6:6).eq.' ') HisGrid(6:6)='m'

  if(MyGrid.eq.HisGrid) then
     Az=0.
     Dmiles=0.
     Dkm=0.0
     El=0.
     HotA=0.
     HotB=0.
     HotABetter=.true.
     go to 900
  endif
  call grid2deg(MyGrid,dlong1,dlat1)
  call grid2deg(HisGrid,dlong2,dlat2)
  eps=1.e-6
  Az=0.
  Dmiles=0.
  Dkm=0.0
  El=0.
  HotA=0.
  HotB=0.
  HotABetter=.true.
  if(abs(dlat1-dlat2).lt.eps .and. abs(dlong1-dlong2).lt.eps) go to 900

  difflong=mod(dlong1-dlong2+720.0,360.0)
  if(abs(dlat1+dlat2).lt.eps .and. abs(difflong-180.0).lt.eps) then
! Antipodes
     Dkm=20400
     go to 900
  endif

  call geodist(dlat1,dlong1,dlat2,dlong2,Az,Baz,Dkm)

  ndkm=Dkm/100
  j=ndkm-4
  if(j.lt.1) j=1
  if(j.gt.21)j=21
  if(Dkm.lt.500.0) then
     El=18.0
  else
     u=(Dkm-100.0*ndkm)/100.0
     El=(1.0-u)*eltab(j) + u*eltab(j+1)
  endif

  daz=daztab(j) + u * (daztab(j+1)-daztab(j))
  Dmiles=Dkm/1.609344

  tmid=mod(UTChours-0.5*(dlong1+dlong2)/15.0+48.0,24.0)
  IamEast=.false.
  if(dlong1.lt.dlong2) IamEast=.true.
  if(dlong1.eq.dlong2 .and. dlat1.gt.dlat2) IamEast=.false.
  azEast=baz
  if(IamEast) azEast=az
  if((azEast.ge.45.0 .and. azEast.lt.135.0) .or.                    &
       (azEast.ge.225.0 .and. azEast.lt.315.0)) then
! The path will be taken as "east-west".
     HotABetter=.true.
     if(abs(tmid-6.0).lt.6.0) HotABetter=.false.
     if((dlat1+dlat2)/2.0 .lt. 0.0) HotABetter=.not.HotABetter
  else
! The path will be taken as "north-south".
     HotABetter=.false.
     if(abs(tmid-12.0).lt.6.0) HotABetter=.true.
  endif
  if(IamEast) then
     HotA = Az - daz
     HotB = Az + daz
  else
     HotA = Az + daz
     HotB = Az - daz
  endif
  if(HotA.lt.0.0)   HotA=HotA+360.0
  if(HotA.gt.360.0) HotA=HotA-360.0
  if(HotB.lt.0.0)   HotB=HotB+360.0
  if(HotB.gt.360.0) HotB=HotB-360.0

900 continue
  naz=nint(Az)
  nel=nint(el)
  nDmiles=nint(Dmiles)
  nDkm=nint(Dkm)
  nHotAz=nint(HotB)
  nHotABetter=0
  if(HotABetter) then
     nHotAz=nint(HotA)
     nHotABetter=1
  endif

999 return
end subroutine azdist


code:
logical function baddata(id2,nz)

  integer*2 id2(nz)

  nadd=1200
  j=0
  smin=1.e30
  smax=-smin
  iz=49*12000/nadd
  do i=1,iz
     sq=0.
     do n=1,nadd
        j=j+1
        x=id2(j)
        sq=sq + x*x
     enddo
     rms=sqrt(sq/nadd)
     if(i.gt.6) then
        smin=min(smin,rms)
        smax=max(smax,rms)
     endif
  enddo

  sratio=smax/(smin+1.e-30)
  baddata=.false.
  if(sratio.gt.1.e30) baddata=.true.

  return
end function baddata


code:
subroutine badmsg(irc,dat,nc1,nc2,ng2)

! Get rid of a few QRA64 false decodes that cannot be correct messages.

  integer dat(12)                           !Decoded message (as 12 integers)

  ic1=ishft(dat(1),22) + ishft(dat(2),16) + ishft(dat(3),10)+         &
       ishft(dat(4),4) + iand(ishft(dat(5),-2),15)

! Test for "......" or "CQ 000"
  if(ic1.eq.262177560 .or. ic1.eq.262177563) then
     irc=-1
     return
  endif

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

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

! Test for blank, -01 to -30, R-01 to R-30, RO, RRR, 73
  if(ig.ge.32401 .and. ig.le.32464) return

  if(ig.ge.14220 .and. ig.le.14229) return  !-41 to -50
  if(ig.ge.14040 .and. ig.le.14049) return  !-31 to -40

  if(ig.ge.13320 .and. ig.le.13329) return  !+00 to +09
  if(ig.ge.13140 .and. ig.le.13149) return  !+10 to +19
  if(ig.ge.12960 .and. ig.le.12969) return  !+20 to +29
  if(ig.ge.12780 .and. ig.le.12789) return  !+30 to +39
  if(ig.ge.12600 .and. ig.le.12609) return  !+40 to +49

  if(ig.ge.12420 .and. ig.le.12429) return  !R-41 to R-50
  if(ig.ge.12240 .and. ig.le.12249) return  !R-31 to R-40

  if(ig.ge.11520 .and. ig.le.11529) return  !R+00 to R+09
  if(ig.ge.11340 .and. ig.le.11349) return  !R+10 to R+19
  if(ig.ge.11160 .and. ig.le.11169) return  !R+20 to R+29
  if(ig.ge.10980 .and. ig.le.10989) return  !R+30 to R+39
  if(ig.ge.10800 .and. ig.le.10809) return  !R+40 to R+49

  if(ic1.eq.nc1 .and. ic2.eq.nc2 .and. ng2.ne.32401 .and. ig.ne.ng2) irc=-1

  return
end subroutine badmsg


code:
subroutine pltanh(x,y)
  isign=+1
  z=x
  if( x.lt.0 ) then
    isign=-1
    z=abs(x)
  endif
  if( z.le. 0.8 ) then
    y=0.83*x
    return
  elseif( z.le. 1.6 ) then
    y=isign*(0.322*z+0.4064)
    return
  elseif( z.le. 3.0 ) then
    y=isign*(0.0524*z+0.8378)
    return
  elseif( z.lt. 7.0 ) then
    y=isign*(0.0012*z+0.9914)
    return
  else
    y=isign*0.9998
    return
  endif
end subroutine pltanh

subroutine platanh(x,y)
  isign=+1
  z=x
  if( x.lt.0 ) then
    isign=-1
    z=abs(x)
  endif
  if( z.le. 0.664 ) then
    y=x/0.83
    return
  elseif( z.le. 0.9217 ) then
    y=isign*(z-0.4064)/0.322
    return
  elseif( z.le. 0.9951 ) then
    y=isign*(z-0.8378)/0.0524
    return
  elseif( z.le. 0.9998 ) then
    y=isign*(z-0.9914)/0.0012
    return
  else
    y=isign*7.0
    return
  endif
end subroutine platanh

subroutine bpdecode144(llr,maxiterations,decoded,niterations)
!
! A log-domain belief propagation decoder for the msk144 code.
! The code is a regular (128,80) code with column weight 3 and row weight 8.
! k9an August, 2016
!
integer, parameter:: N=128, K=80, M=N-K
integer*1 codeword(N),cw(N)
integer*1 colorder(N)
integer*1 decoded(K)
integer Nm(8,M)  ! 8 bits per check
integer Mn(3,N)  ! 3 checks per bit
integer synd(M)
real tov(3,N)    ! single precision seems to be adequate in log-domain
real toc(8,M)
real tanhtoc(8,M)
real zn(N)
real llr(N)
real Tmn

data colorder/0,1,2,3,4,5,6,7,8,9, &
              10,11,12,13,14,15,24,26,29,30, &
              32,43,44,47,60,77,79,97,101,111, &
              96,38,64,53,93,34,59,94,74,90, &
              108,123,85,57,70,25,69,62,48,49, &
              50,51,52,33,54,55,56,21,58,36, &
              16,61,23,63,20,65,66,67,68,46, &
              22,71,72,73,31,75,76,45,78,17, &
              80,81,82,83,84,42,86,87,88,89, &
              39,91,92,35,37,95,19,27,98,99, &
              100,28,102,103,104,105,106,107,40,109, &
              110,18,112,113,114,115,116,117,118,119, &
              120,121,122,41,124,125,126,127/

data Mn/               &
   1,  14,  38, &
   2,   4,  41, &
   3,  19,  39, &
   5,  29,  34, &
   6,  35,  40, &
   7,  20,  45, &
   8,  28,  48, &
   9,  22,  25, &
  10,  24,  36, &
  11,  12,  37, &
  13,  43,  44, &
  15,  18,  46, &
  16,  17,  47, &
  21,  32,  33, &
  23,  30,  31, &
  26,  27,  42, &
   1,  12,  46, &
   2,  36,  38, &
   3,   5,  10, &
   4,   9,  23, &
   6,  13,  39, &
   7,  15,  17, &
   8,  18,  27, &
  11,  33,  40, &
  14,  28,  44, &
  16,  29,  31, &
  19,  20,  22, &
  21,  30,  42, &
  24,  26,  47, &
  25,  37,  48, &
  32,  34,  45, &
   8,  35,  41, &
  12,  31,  43, &
   1,  19,  21, &
   2,  43,  45, &
   3,   4,  11, &
   5,  18,  33, &
   6,  25,  47, &
   7,  28,  30, &
   9,  14,  34, &
  10,  35,  42, &
  13,  15,  22, &
  16,  37,  38, &
  17,  41,  44, &
  20,  24,  29, &
  18,  23,  39, &
  12,  26,  32, &
  27,  38,  40, &
  15,  36,  48, &
   2,  30,  46, &
   1,   4,  13, &
   3,  28,  32, &
   5,  43,  47, &
   6,  34,  46, &
   7,   9,  40, &
   8,  11,  45, &
  10,  17,  23, &
  14,  31,  35, &
  16,  22,  42, &
  19,  37,  44, &
  20,  33,  48, &
  21,  24,  41, &
  25,  27,  29, &
  26,  39,  48, &
  19,  31,  36, &
   1,   5,   7, &
   2,  29,  39, &
   3,  16,  46, &
   4,  26,  37, &
   6,  28,  45, &
   8,  22,  33, &
   9,  21,  43, &
  10,  25,  38, &
  11,  14,  24, &
  12,  17,  40, &
  13,  27,  30, &
  15,  32,  35, &
  18,  44,  47, &
  20,  23,  36, &
  34,  41,  42, &
   1,  32,  48, &
   2,   3,  33, &
   4,  29,  42, &
   5,  14,  37, &
   6,   7,  36, &
   8,   9,  39, &
  10,  13,  19, &
  11,  18,  30, &
  12,  16,  20, &
  15,  29,  44, &
  17,  34,  38, &
   6,  21,  22, &
  23,  32,  40, &
  24,  27,  46, &
  25,  41,  45, &
   7,  26,  43, &
  28,  31,  47, &
  20,  35,  38, &
   1,  33,  41, &
   2,  42,  44, &
   3,  23,  48, &
   4,  31,  45, &
   5,   8,  30, &
   9,  16,  36, &
  10,  40,  47, &
  11,  17,  46, &
  12,  21,  34, &
  13,  24,  28, &
  14,  18,  43, &
  15,  25,  26, &
  19,  27,  35, &
  22,  37,  39, &
   1,  16,  18, &
   2,   6,  20, &
   3,  30,  43, &
   4,  28,  33, &
   5,  22,  23, &
   7,  39,  42, &
   8,  12,  38, &
   9,  35,  46, &
  10,  27,  32, &
  11,  15,  34, &
  13,  36,  37, &
  14,  41,  47, &
  17,  21,  25, &
  19,  29,  45, &
  24,  31,  48, &
  26,  40,  44/

data Nm/               &
   1,  17,  34,  51,  66,  81,  99, 113, &
   2,  18,  35,  50,  67,  82, 100, 114, &
   3,  19,  36,  52,  68,  82, 101, 115, &
   2,  20,  36,  51,  69,  83, 102, 116, &
   4,  19,  37,  53,  66,  84, 103, 117, &
   5,  21,  38,  54,  70,  85,  92, 114, &
   6,  22,  39,  55,  66,  85,  96, 118, &
   7,  23,  32,  56,  71,  86, 103, 119, &
   8,  20,  40,  55,  72,  86, 104, 120, &
   9,  19,  41,  57,  73,  87, 105, 121, &
  10,  24,  36,  56,  74,  88, 106, 122, &
  10,  17,  33,  47,  75,  89, 107, 119, &
  11,  21,  42,  51,  76,  87, 108, 123, &
   1,  25,  40,  58,  74,  84, 109, 124, &
  12,  22,  42,  49,  77,  90, 110, 122, &
  13,  26,  43,  59,  68,  89, 104, 113, &
  13,  22,  44,  57,  75,  91, 106, 125, &
  12,  23,  37,  46,  78,  88, 109, 113, &
   3,  27,  34,  60,  65,  87, 111, 126, &
   6,  27,  45,  61,  79,  89,  98, 114, &
  14,  28,  34,  62,  72,  92, 107, 125, &
   8,  27,  42,  59,  71,  92, 112, 117, &
  15,  20,  46,  57,  79,  93, 101, 117, &
   9,  29,  45,  62,  74,  94, 108, 127, &
   8,  30,  38,  63,  73,  95, 110, 125, &
  16,  29,  47,  64,  69,  96, 110, 128, &
  16,  23,  48,  63,  76,  94, 111, 121, &
   7,  25,  39,  52,  70,  97, 108, 116, &
   4,  26,  45,  63,  67,  83,  90, 126, &
  15,  28,  39,  50,  76,  88, 103, 115, &
  15,  26,  33,  58,  65,  97, 102, 127, &
  14,  31,  47,  52,  77,  81,  93, 121, &
  14,  24,  37,  61,  71,  82,  99, 116, &
   4,  31,  40,  54,  80,  91, 107, 122, &
   5,  32,  41,  58,  77,  98, 111, 120, &
   9,  18,  49,  65,  79,  85, 104, 123, &
  10,  30,  43,  60,  69,  84, 112, 123, &
   1,  18,  43,  48,  73,  91,  98, 119, &
   3,  21,  46,  64,  67,  86, 112, 118, &
   5,  24,  48,  55,  75,  93, 105, 128, &
   2,  32,  44,  62,  80,  95,  99, 124, &
  16,  28,  41,  59,  80,  83, 100, 118, &
  11,  33,  35,  53,  72,  96, 109, 115, &
  11,  25,  44,  60,  78,  90, 100, 128, &
   6,  31,  35,  56,  70,  95, 102, 126, &
  12,  17,  50,  54,  68,  94, 106, 120, &
  13,  29,  38,  53,  78,  97, 105, 124, &
   7,  30,  49,  61,  64,  81, 101, 127/

nrw=8
ncw=3

toc=0
tov=0
tanhtoc=0

! initial messages to checks
do j=1,M
  do i=1,nrw
    toc(i,j)=llr((Nm(i,j)))
  enddo
enddo

ncnt=0

do iter=0,maxiterations

! Update bit log likelihood ratios
  do i=1,N
    zn(i)=llr(i)+sum(tov(1:ncw,i))
  enddo

! Check to see if we have a codeword
  cw=0
  where( zn .gt. 0. ) cw=1
  ncheck=0
  do i=1,M
    synd(i)=sum(cw(Nm(:,i)))
    if( mod(synd(i),2) .ne. 0 ) ncheck=ncheck+1
  enddo

  if( ncheck .eq. 0 ) then ! we have a codeword
    niterations=iter
    codeword=cw(colorder+1)
    decoded=codeword(M+1:N)
    return
  endif

  if( iter.gt.0 ) then  ! this code block implements an early stopping criterion
    nd=ncheck-nclast
    if( nd .lt. 0 ) then ! # of unsatisfied parity checks decreased
      ncnt=0  ! reset counter
    else
      ncnt=ncnt+1
    endif
!    write(*,*) iter,ncheck,nd,ncnt
    if( ncnt .ge. 3 .and. iter .ge. 5 .and. ncheck .gt. 10) then
      niterations=-1
      return
    endif
  endif
  nclast=ncheck

! Send messages from bits to check nodes
  do j=1,M
    do i=1,nrw
      ibj=Nm(i,j)
      toc(i,j)=zn(ibj)
      do kk=1,ncw ! subtract off what the bit had received from the check
        if( Mn(kk,ibj) .eq. j ) then  ! Mn(3,128)
          toc(i,j)=toc(i,j)-tov(kk,ibj)
        endif
      enddo
    enddo
  enddo

! send messages from check nodes to variable nodes
  do i=1,M
    tanhtoc(1:nrw,i)=tanh(-toc(1:nrw,i)/2)
  enddo

  do j=1,N
    do i=1,ncw
      ichk=Mn(i,j)  ! Mn(:,j) are the checks that include bit j
      Tmn=product(tanhtoc(:,ichk),mask=Nm(:,ichk).ne.j)
      call platanh(-Tmn,y)
      tov(i,j)=2*y
    enddo
  enddo

enddo
niterations=-1
end subroutine bpdecode144


code:
subroutine calibrate(data_dir,iz,a,b,rms,sigmaa,sigmab,irc)

! Average groups of frequency-calibration measurements, then fit a
! straight line for slope and intercept.

  parameter (NZ=1000)
  implicit real*8 (a-h,o-z)
  character*(*) data_dir
  character*256 infile,outfile
  character*8 cutc,cutc1
  character*1 c1
  real*8 fd(NZ),deltaf(NZ),r(NZ),rmsd(NZ)
  integer nn(NZ)

  infile=trim(data_dir)//'fmt.all'
  outfile=trim(data_dir)//'fcal2.out'

  open(10,file=trim(infile),status='old',err=996)
  open(12,file=trim(outfile),status='unknown',err=997)

  nkhz0=0
  sum=0.d0
  sumsq=0.d0
  n=0
  j=0
  do i=1,99999
     read(10,*,end=10,err=995) cutc,nkHz,ncal,noffset,faudio,df,dblevel,snr
     if((nkHz.ne.nkHz0) .and. i.ne.1) then
        ave=sum/n
        rms=0.d0
        if(n.gt.1) then
           rms=sqrt(abs(sumsq - sum*sum/n)/(n-1.d0))
        endif
        fMHz=0.001d0*nkHz0
        j=j+1
        fd(j)=fMHz
        deltaf(j)=ave
        r(j)=0.d0
        rmsd(j)=rms
        nn(j)=n
        sum=0.d0
        sumsq=0.d0
        n=0
     endif
     dial_error=faudio-noffset
     sum=sum + dial_error
     sumsq=sumsq + dial_error**2
     n=n+1
     if(n.eq.1) then
        cutc1=cutc
        ncal0=ncal
     endif
     nkHz0=nkHz
  enddo

10 ave=sum/n
  rms=0.d0
  if(n.gt.0) then
     rms=sqrt((sumsq - sum*sum/n)/(n-1.d0))
  endif
  fMHz=0.001d0*nkHz
  j=j+1
  fd(j)=fMHz
  deltaf(j)=ave
  r(j)=0.d0
  rmsd(j)=rms
  nn(j)=n
  iz=j
  if(iz.lt.2) go to 998

  call fitcal(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)

  write(12,1002)
1002 format('    Freq      DF     Meas Freq    N    rms    Resid'/        &
            '   (MHz)     (Hz)      (MHz)          (Hz)     (Hz)'/        &
            '----------------------------------------------------')
  irc=0
  do i=1,iz
     fm=fd(i) + 1.d-6*deltaf(i)
     c1=' '
     if(rmsd(i).gt.1.0d0) c1='*'
     write(12,1012)  fd(i),deltaf(i),fm,nn(i),rmsd(i),r(i),c1
1012 format(f8.3,f9.3,f14.9,i4,f7.2,f9.3,1x,a1)
  enddo
  go to 999

995 irc=-4; iz=i; go to 999
996 irc=-1; go to 999
997 irc=-2; go to 999
998 irc=-3
999 continue
  close(10)
  close(12)

  return
end subroutine calibrate


code:
subroutine ccf2(ss,nz,nflip,ccfbest,xlagpk)

!  parameter (LAGMIN=-86,LAGMAX=258)
  parameter (LAGMIN=-112,LAGMAX=258)  ! Look for DT from -3.6s to +5.0s
  real ss(nz)
  real ccf(-LAGMAX:LAGMAX)
  integer npr(126)

! The JT65 pseudo-random sync pattern:
  data npr/                                    &
      1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0, &
      0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1, &
      0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1, &
      0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1, &
      1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1, &
      0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1, &
      1,1,1,1,1,1/
  save

  ccfbest=0.
  lag1=LAGMIN
  lag2=LAGMAX
  do lag=lag1,lag2
     s0=0.
     s1=0.
     do i=1,126
        j=16*(i-1)+1 + lag
        if(j.ge.1 .and. j.le.nz-8) then
           x=ss(j)
           if(npr(i).eq.0) then
              s0=s0 + x
           else
              s1=s1 + x
           endif
        endif
     enddo
     ccf(lag)=nflip*(s1-s0)
     if(ccf(lag).gt.ccfbest) then
        ccfbest=ccf(lag)
        lagpk=lag
        xlagpk=lagpk
     endif
  enddo
  if( lagpk.gt.-LAGMAX .and. lagpk.lt.LAGMAX) then
     call peakup(ccf(lagpk-1),ccf(lagpk),ccf(lagpk+1),dx)
     xlagpk=lagpk+dx
  endif
  return
end subroutine ccf2


code:
subroutine ccf65(ss,nhsym,ssmax,sync1,dt1,flipk,syncshort,snr2,dt2)

  parameter (NFFT=512,NH=NFFT/2)
  real ss(322)                     !Input: half-symbol normalized powers
  real s(NFFT)                     !CCF = ss*pr
  complex cs(0:NH)                 !Complex FT of s
  real s2(NFFT)                    !CCF = ss*pr2
  complex cs2(0:NH)                !Complex FT of s2
  real pr(NFFT)                    !JT65 pseudo-random sync pattern
  complex cpr(0:NH)                !Complex FT of pr
  real pr2(NFFT)                   !JT65 shorthand pattern
  complex cpr2(0:NH)               !Complex FT of pr2
  real tmp1(322)
  real ccf(-11:54)
  logical first
  integer npr(126)
  data first/.true./
  equivalence (s,cs),(pr,cpr),(s2,cs2),(pr2,cpr2)
  save

! The JT65 pseudo-random sync pattern:
  data npr/                                        &
      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/

  if(first) then
! Initialize pr, pr2; compute cpr, cpr2.
     fac=1.0/NFFT
     do i=1,NFFT
        pr(i)=0.
        pr2(i)=0.
        k=2*mod((i-1)/8,2)-1
        if(i.le.NH) pr2(i)=fac*k
     enddo
     do i=1,126
        j=2*i
        pr(j)=fac*(2*npr(i)-1)
! Not sure why, but it works significantly better without the following line:
!        pr(j-1)=pr(j)
     enddo
     call four2a(cpr,NFFT,1,-1,0)
     call four2a(cpr2,NFFT,1,-1,0)
     first=.false.
  endif

! Look for JT65 sync pattern and shorthand square-wave pattern.
  ccfbest=0.
  ccfbest2=0.
  do i=1,nhsym-1
     s(i)=min(ssmax,ss(i)+ss(i+1))
!     s(i)=ss(i)+ss(i+1)
  enddo

  call pctile(s,nhsym-1,50,base)
  s(1:nhsym-1)=s(1:nhsym-1)-base
  s(nhsym:NFFT)=0.
  call four2a(cs,NFFT,1,-1,0)               !Real-to-complex FFT
  do i=0,NH
!     cs2(i)=cs(i)*conjg(cpr2(i))            !Mult by complex FFT of pr2
     cs(i)=cs(i)*conjg(cpr(i))              !Mult by complex FFT of pr
  enddo
  call four2a(cs,NFFT,1,1,-1)               !Complex-to-real inv-FFT
!  call four2a(cs2,NFFT,1,1,-1)              !Complex-to-real inv-FFT

  do lag=-11,54                             !Check for best JT65 sync
     j=lag
     if(j.lt.1) j=j+NFFT
     ccf(lag)=s(j)
!     if(abs(ccf(lag)).gt.ccfbest) then
     if(ccf(lag).gt.ccfbest) then           !No inverted sync for use at HF
!        ccfbest=abs(ccf(lag))
        ccfbest=ccf(lag)
        lagpk=lag
        flipk=1.0
!        if(ccf(lag).lt.0.0) flipk=-1.0
     endif
  enddo

!  do lag=-11,54                             !Check for best shorthand
!     ccf2=s2(lag+28)
!     if(ccf2.gt.ccfbest2) then
!        ccfbest2=ccf2
!        lagpk2=lag
!     endif
!  enddo

! Find rms level on baseline of "ccfblue", for normalization.
  sum=0.
  do lag=-11,54
     if(abs(lag-lagpk).gt.1) sum=sum + ccf(lag)
  enddo
  base=sum/50.0
  sq=0.
  do lag=-11,54
     if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag)-base)**2
  enddo
  rms=sqrt(sq/49.0)
  sync1=ccfbest/rms - 4.0
  dt1=lagpk*(2048.0/11025.0) - 2.5

! Find base level for normalizing snr2.
  do i=1,nhsym
     tmp1(i)=ss(i)
  enddo
  call pctile(tmp1,nhsym,40,base)
  snr2=0.398107*ccfbest2/base                !### empirical
  syncshort=0.5*ccfbest2/rms - 4.0           !### better normalizer than rms?
!  dt2=(2.5 + lagpk2*(2048.0/11025.0))
  dt2=0.

  return
end subroutine ccf65


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

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

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

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

#define ENCODE_RS encode_rs_char
#define DECODE_RS decode_rs_char
#define INIT_RS init_rs_char
#define FREE_RS free_rs_char

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


code:
subroutine chkcall(w,bc,cok)

! Check "w" to see if it could be a valid standard callsign or a valid
! compound callsign.
! Return base call "bc" and a logical "cok" indicator.

  character w*13                            !A putative callsign
  character bc*6                            !Base call (tentative)
  character c*1
  logical cok,isdigit,isletter

  isdigit(c)=(ichar(c).ge.ichar('0')) .and. (ichar(c).le.ichar('9'))
  isletter(c)=(ichar(c).ge.ichar('A')) .and. (ichar(c).le.ichar('Z'))

  cok=.true.
  bc=w(1:6)
  n1=len_trim(w)
  if(n1.gt.11) go to 100
  if(index(w,'.').ge.1) go to 100
  if(index(w,'+').ge.1) go to 100
  if(index(w,'-').ge.1) go to 100
  if(index(w,'?').ge.1) go to 100
  if(n1.gt.6 .and. index(w,'/').le.0) go to 100

  i0=index(w,'/')
  if(max(i0-1,n1-i0).gt.6) go to 100      !Base call must be < 7 characters
  if(i0.ge.2 .and. i0.le.n1-1) then       !Extract base call from compound call
     if(i0-1.le.n1-i0) bc=w(i0+1:n1)//'   '
     if(i0-1.gt.n1-i0) bc=w(1:i0-1)//'   '
  endif

  nbc=len_trim(bc)
  if(nbc.gt.6) go to 100  !Base call should have no more than 6 characters

! One of first two characters (c1 or c2) must be a letter
  if((.not.isletter(bc(1:1))) .and. (.not.isletter(bc(2:2)))) go to 100
  if(bc(1:1).eq.'Q') go to 100              !Calls don't start with Q

! Must have a digit in 2nd or 3rd position
  i1=0
  if(isdigit(bc(2:2))) i1=2
  if(isdigit(bc(3:3))) i1=3
  if(i1.eq.0) go to 100

! Callsign must have a suffix of 1-3 letters
  if(i1.eq.nbc) go to 100
  n=0
  do i=i1+1,nbc
     j=ichar(bc(i:i))
     if(j.lt.ichar('A') .or. j.gt.ichar('Z')) go to 100
     n=n+1
  enddo
  if(n.ge.1 .and. n.le.3) go to 200

100 cok=.false.

200 return
end subroutine chkcall


code:
program chkfft

! Tests and times one-dimensional FFTs computed by four2a().
! An all-Fortran version of four2a() is available, but the preferred
! version uses calls to the FFTW library.

  parameter (NMAX=8*1024*1024)            !Maximum FFT length
  complex a(NMAX),b(NMAX)
  real ar(NMAX),br(NMAX)
  real mflops
  character infile*12,arg*8
  logical list
  common/patience/npatience
  equivalence (a,ar),(b,br)

  nargs=iargc()
  if(nargs.ne.5) then
     print*,'Usage: chkfft <nfft | infile> nr nw nc np'
     print*,'       nfft:   length of FFT'
     print*,'       nfft=0: do lengths 2^n, n=2^4 to 2^23'
     print*,'       infile: name of file with nfft values, one per line'
     print*,'       nr:     0/1 to not read (or read) wisdom'
     print*,'       nw:     0/1 to not write (or write) wisdom'
     print*,'       nc:     0/1 for real or complex data'
     print*,'       np:     0-4 patience for finding best algorithm'
     go to 999
  endif

  list=.false.
  nfft=-1
  call getarg(1,infile)
  open(10,file=infile,status='old',err=1)
  list=.true.                          !A valid file name was provided
  go to 2
1 read(infile,*) nfft                  !Takje first argument to be nfft
2  call getarg(2,arg)
  read(arg,*) nr
  call getarg(3,arg)
  read(arg,*) nw
  call getarg(4,arg)
  read(arg,*) ncomplex
  call getarg(5,arg)
  read(arg,*) npatience

  call sgran()

  if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
1000 format(/'infile: ',a12,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)
  if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
1002 format(/'nfft: ',i10,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)

  open(12,file='chkfft.out',status='unknown')
  open(13,file='fftwf_wisdom.dat',status='unknown')

  if(nr.ne.0) then
     call import_wisdom_from_file(isuccess,13)
     if(isuccess.eq.0) then
        write(*,1010)
1010    format('Failed to import FFTW wisdom.')
        go to 999
     endif
  endif

  idum=-1                               !Set random seed
  ndim=1                                !One-dimensional transforms
  do i=1,NMAX                           !Set random data
     x=gran()
     y=gran()
     b(i)=cmplx(x,y)                    !Generate random data
  enddo

  iters=1000000
  if(list .or. (nfft.gt.0)) then
     n1=1
     n2=1
     if(nfft.eq.-1) n2=999999
     write(*,1020)
1020 format('    NFFT     Time        rms      MHz   MFlops  iters',    &
          '  tplan'/61('-'))
  else
     n1=4
     n2=23
     write(*,1030)
1030 format(' n   N=2^n     Time        rms      MHz   MFlops  iters',  &
          '  tplan'/63('-'))
  endif

  do ii=n1,n2                           !Test one or more FFT lengths
     if(list) then
        read(10,*,end=900) nfft         !Read nfft from file
     else if(n2.gt.n1) then
        nfft=2**ii                      !Do powers of 2
     endif

     iformf=1
     iformb=1
     if(ncomplex.eq.0) then
        iformf=0                        !Real-to-complex transform
        iformb=-1                       !Complex-to-real (inverse) transform
     endif

     if(nfft.gt.NMAX) go to 900
     a(1:nfft)=b(1:nfft)                !Copy test data into a()
     t0=second()
     call four2a(a,nfft,ndim,-1,iformf) !Get planning time for forward FFT
     call four2a(a,nfft,ndim,+1,iformb) !Get planning time for backward FFT
     t2=second()
     tplan=t2-t0                        !Total planning time for this length

     total=0.
     do iter=1,iters                    !Now do many iterations
        a(1:nfft)=b(1:nfft)             !Copy test data into a()

        t0=second()
        call four2a(a,nfft,ndim,-1,iformf) !Forward FFT
        call four2a(a,nfft,ndim,+1,iformb) !Backward FFT on same data
        t1=second()
        total=total+t1-t0
        if(total.ge.1.0) go to 40       !Cut iterations short if t>1 s
     enddo
     iter=iters

40   time=0.5*total/iter                !Time for one FFT of current length
     tplan=0.5*tplan-time               !Planning time for one FFT
     if(tplan.lt.0) tplan=0.
     a(1:nfft)=a(1:nfft)/nfft

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

     freq=1.e-6*nfft/time
     mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
     if(n2.eq.1 .or. n2.eq.999999) then
        write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
        write(12,1050) nfft,time,rms,freq,mflops,iter,tplan
1050    format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     else
        write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
        write(12,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
1060    format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     endif
     if(mod(ii,50).eq.0) call four2a(0,-1,0,0,0)
  enddo

900  continue
  if(nw.eq.1) then
     rewind 13
     call export_wisdom_to_file(13)
!     write(*,1070)
!1070 format(/'Exported FFTW wisdom')
  endif

999 call four2a(0,-1,0,0,0)
end program chkfft


code:
		Brief Description of chkfft, by K1JT
		------------------------------------

Discrete Fourier transforms (DFTs) are found at the root of most
digital signal processing tasks. In WSJT and its sister programs the
transforms are done using the FFTW library, and subroutine four2
provides a convenient interface to the library.  Program chkfft is a
command-line utility offering a convenient way to test FFT execution
times under a variety of circumstances.

To compile chkfft in Linux:

$ gfortran -o chkfft chkfft.f90 four2a.f90 f77_wisdom.f90 gran.c -lfftw3f

To compile chkfft in Windows (you may need to customize the hard-coded
path shown here for libfftw3f-3.dll):

> gfortran -o chkfft chkfft.f90 four2a.f90 f77_wisdom.f90 gran.c \
  /JTSDK-QT/appsupport/runtime/libfftw3f-3.dll

To see a brief usage message, type chkfft at the command prompt:

$ chkfft
 Usage: chkfft <nfft | infile> nr nw nc np
        nfft:   length of FFT
        nfft=0: do lengths 2^n, n=2^4 to 2^23
        infile: name of file with nfft values, one per line
        nr:     0/1 to not read (or read) wisdom
        nw:     0/1 to not write (or write) wisdom
        nc:     0/1 for real or complex data
        np:     0-4 patience for finding best algorithm

As an example, to measure the speed of a complex DFT of length 131072:

#######################################################################
$ chkfft 131072 0 1 1 2

nfft:     131072   nr: 0   nw 1   nc: 1   np: 2

    NFFT     Time        rms      MHz   MFlops  iters  tplan
-------------------------------------------------------------
  131072  0.0021948  0.00000032  59.72  5076.1     231   2.9
#######################################################################

Program output shows that on the test machine the average time for one
forward (or inverse) transform of length N=131072 is about 2.2 ms,
corresponding to slightly over 5 GFlops computing speed.  The planning
time in FFTW was 2.9 s.

Running the command again with parameter nr=1 will use the
"wisdom" already accumulated for complex N=131072 FFTs.  The execution
speed will be essentially the same, but no planning time is required:

#######################################################################
$ chkfft 131072 1 1 1 2

nfft:     131072   nr: 1   nw 1   nc: 1   np: 2

    NFFT     Time        rms      MHz   MFlops  iters  tplan
-------------------------------------------------------------
  131072  0.0021575  0.00000032  60.75  5164.0     235   0.0
#######################################################################

Optimized algorithms can compute DFTs much faster for lengths that are
the product of small integers.  Length N=131072 = 2^17 is a good
example, and FFTs should be very efficient.  For comparison, look at
the speed for N=131071, a prime number.  The average time is now about
7 times larger:

#######################################################################
C:\JTSDK-QT\src\wsjtx\lib>chkfft 131071 1 1 1 2

nfft:     131071   nr: 1   nw 1   nc: 1   np: 2

    NFFT     Time        rms      MHz   MFlops  iters  tplan
-------------------------------------------------------------
  131071  0.0153637  0.00000065   8.53   725.2      33   5.6
#######################################################################

Here's an example that measures execution times for all integral
power-of-2 lengths from 2^4 to 2^23:

#######################################################################
$ chkfft 0 1 1 1 2

nfft:          0   nr: 1   nw 1   nc: 1   np: 2

  n   N=2^n    Time        rms      MHz   MFlops  iters  tplan
---------------------------------------------------------------
 4      16  0.0000003  0.00000014  58.61  1172.2 1000000   0.0
 5      32  0.0000004  0.00000016  89.19  2229.6 1000000   0.0
 6      64  0.0000006  0.00000016 109.44  3283.2  866975   0.0
 7     128  0.0000009  0.00000021 135.92  4757.1  538369   0.0
 8     256  0.0000016  0.00000020 158.40  6335.8  313701   0.0
 9     512  0.0000032  0.00000021 162.53  7313.8  160943   0.1
10    1024  0.0000067  0.00000023 152.53  7626.5   75521   0.1
11    2048  0.0000136  0.00000025 150.42  8273.3   37239   0.2
12    4096  0.0000316  0.00000027 129.75  7784.8   16060   0.3
13    8192  0.0000720  0.00000026 113.75  7393.8    7040   0.5
14   16384  0.0001620  0.00000028 101.11  7078.0    3129   0.9
15   32768  0.0003227  0.00000030 101.53  7615.1    1571   1.7
16   65536  0.0010020  0.00000030  65.41  5232.5     506   4.1
17  131072  0.0021575  0.00000032  60.75  5164.0     235   0.0
18  262144  0.0053937  0.00000032  48.60  4374.2      94   3.6
19  524288  0.0190668  0.00000034  27.50  2612.2      27   6.8
20 1048576  0.0468001  0.00000035  22.41  2240.5      11   2.4
21 2097152  0.0936012  0.00000036  22.41  2352.5       6  31.6
22 4194304  0.1949997  0.00000037  21.51  2366.0       3   9.8
23 8388608  0.4212036  0.00000038  19.92  2290.3       2 112.9
#######################################################################

Test data for all transforms is gaussian random noise of zero mean and
standard deviation 1.  Tabulated values of "rms" are the
root-mean-square differences between the original data and the
back-transfmred data.

File nfft.dat contains all numbers between 2^3 and 2^23 with no factor
greater than 7, followed by their factors.  These numbers are good
choices for FFT lengths.  File all_fft.out gives the result on one
machine of running the command

$ chkfft nfft.dat 0 1 1 2

Take note: this task may take as much as 24 hours, or even more!


code:
program chkfft

! Tests and times one-dimensional FFTs computed by FFTW3
  use FFTW3
  parameter (NMAX=8*1024*1024)           !Maximum FFT length
  complex a(NMAX),b(NMAX),c(NMAX)
  real ar(NMAX),br(NMAX),cr(NMAX)
  real mflops
  integer*8 plan1,plan2                  !Pointers to stored plans
  character infile*12,arg*8
  logical list
  common/patience/npatience
  equivalence (a,ar),(b,br),(c,cr)
!  include 'fftw3.f90'                    !FFTW definitions

  nargs=iargc()
  if(nargs.ne.6) then
     print*,'Usage: chkfft <nfft | infile> nr nw nc np inplace'
     print*,'       nfft:    length of FFT'
     print*,'       nfft=0:  do lengths 2^n, n=2^4 to 2^23'
     print*,'       infile:  name of file with nfft values, one per line'
     print*,'       nr:      0/1 to not read (or read) wisdom'
     print*,'       nw:      0/1 to not write (or write) wisdom'
     print*,'       nc:      0/1 for real or complex data'
     print*,'       np:      0-4 patience for finding best algorithm'
     print*,'       inplace: 1 for inplace, 0 otherwise'
     go to 999
  endif

  list=.false.
  nfft=-1
  call getarg(1,infile)
  open(10,file=infile,status='old',err=1)
  list=.true.                          !A valid file name was provided
  go to 2
1 read(infile,*) nfft                  !Take first argument to be nfft
2  call getarg(2,arg)
  read(arg,*) nr
  call getarg(3,arg)
  read(arg,*) nw
  call getarg(4,arg)
  read(arg,*) ncomplex
  call getarg(5,arg)
  read(arg,*) npatience
  call getarg(6,arg)
  read(arg,*) inplace

  if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
1000 format(/'infile: ',a12,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)
  if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
1002 format(/'nfft: ',i10,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)

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

  open(12,file='chkfft.out',status='unknown')
  open(13,file='fftwf_wisdom.dat',status='unknown')

  if(nr.ne.0) then
     call import_wisdom_from_file(isuccess,13)
     if(isuccess.eq.0) then
        write(*,1010)
1010    format('Failed to import FFTW wisdom.')
        go to 999
     endif
  endif

  idum=-1                               !Set random seed
  ndim=1                                !One-dimensional transforms
  do i=1,NMAX                           !Set random data
     x=gran()
     y=gran()
     b(i)=cmplx(x,y)                    !Generate random data
  enddo

  iters=1000000
  if(list .or. (nfft.gt.0)) then
     n1=1
     n2=1
     if(nfft.eq.-1) n2=999999
     write(*,1020)
1020 format('    NFFT     Time        rms      MHz   MFlops  iters',    &
          '  tplan'/61('-'))
  else
     n1=4
     n2=23
     write(*,1030)
1030 format(' n   N=2^n     Time        rms      MHz   MFlops  iters',  &
          '  tplan'/63('-'))
  endif

  do ii=n1,n2                           !Test one or more FFT lengths
     if(list) then
        read(10,*,end=900) nfft         !Read nfft from file
     else if(n2.gt.n1) then
        nfft=2**ii                      !Do powers of 2
     endif

     iformf=1
     iformb=1
     if(ncomplex.eq.0) then
        iformf=0                        !Real-to-complex transform
        iformb=-1                       !Complex-to-real (inverse) transform
     endif

     if(nfft.gt.NMAX) go to 900
     a(1:nfft)=b(1:nfft)                !Copy test data into a()
     t0=second()
     if(inplace.ne.0) then
        if(ncomplex.ne.0) then
           call sfftw_plan_dft_1d(plan1,nfft,a,a,FFTW_FORWARD,nflags)
           call sfftw_plan_dft_1d(plan2,nfft,a,a,FFTW_BACKWARD,nflags)
        else
           call sfftw_plan_dft_r2c_1d(plan1,nfft,a,a,nflags)
           call sfftw_plan_dft_c2r_1d(plan2,nfft,a,a,nflags)
        endif
     else
        if(ncomplex.ne.0) then
           call sfftw_plan_dft_1d(plan1,nfft,a,c,FFTW_FORWARD,nflags)
           call sfftw_plan_dft_1d(plan2,nfft,c,a,FFTW_BACKWARD,nflags)
        else
           call sfftw_plan_dft_r2c_1d(plan1,nfft,a,c,nflags)
           call sfftw_plan_dft_c2r_1d(plan2,nfft,c,a,nflags)
        endif
     endif

     t2=second()
     tplan=t2-t0                        !Total planning time for this length

     total=0.
     do iter=1,iters                    !Now do many iterations
        a(1:nfft)=b(1:nfft)             !Copy test data into a()
        t0=second()
        call sfftw_execute(plan1)
        call sfftw_execute(plan2)
        t1=second()
        total=total+t1-t0
        if(total.ge.1.0) go to 40       !Cut iterations short if t>1 s
     enddo
     iter=iters

40   time=0.5*total/iter                !Time for one FFT of current length
     tplan=0.5*tplan-time               !Planning time for one FFT
     if(tplan.lt.0) tplan=0.
     a(1:nfft)=a(1:nfft)/nfft

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

     freq=1.e-6*nfft/time
     mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
     if(n2.eq.1 .or. n2.eq.999999) then
        write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
        write(12,1050) nfft,time,rms,freq,mflops,iter,tplan
1050    format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     else
        write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
        write(12,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
1060    format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     endif
  enddo

900  continue
  if(nw.eq.1) then
     rewind 13
     call export_wisdom_to_file(13)
!     write(*,1070)
!1070 format(/'Exported FFTW wisdom')
  endif

  call sfftw_destroy_plan(plan1)
  call sfftw_destroy_plan(plan2)

999 end program chkfft


code:
program chkfft

! Tests and times one-dimensional FFTs computed by FFTW3
  use, intrinsic :: iso_c_binding
  use FFTW3
  parameter (NMAX=8*1024*1024)           !Maximum FFT length
  complex a(NMAX),b(NMAX),c(NMAX)
  real ar(NMAX),br(NMAX),cr(NMAX)
  real mflops
!  integer*8 plan1,plan2                  !Pointers to stored plans
  type(C_PTR) :: plan1,plan2              !Pointers to FFTW plans
  character infile*12,arg*8
  logical list
  common/patience/npatience
  equivalence (a,ar),(b,br),(c,cr)
!  include 'fftw3.f90'                    !FFTW definitions

  nargs=iargc()
  if(nargs.ne.5) then
     print*,'Usage: chkfft3 <nfft | infile> nr nw nc np'
     print*,'       nfft:    length of FFT'
     print*,'       nfft=0:  do lengths 2^n, n=2^4 to 2^23'
     print*,'       infile:  name of file with nfft values, one per line'
     print*,'       nr:      0/1 to not read (or read) wisdom'
     print*,'       nw:      0/1 to not write (or write) wisdom'
     print*,'       nc:      0/1 for real or complex data'
     print*,'       np:      0-4 patience for finding best algorithm'
     go to 999
  endif

  list=.false.
  nfft=-1
  call getarg(1,infile)
  open(10,file=infile,status='old',err=1)
  list=.true.                          !A valid file name was provided
  go to 2
1 read(infile,*) nfft                  !Take first argument to be nfft
2  call getarg(2,arg)
  read(arg,*) nr
  call getarg(3,arg)
  read(arg,*) nw
  call getarg(4,arg)
  read(arg,*) ncomplex
  call getarg(5,arg)
  read(arg,*) npatience

  if(list) write(*,1000) infile,nr,nw,ncomplex,npatience
1000 format(/'infile: ',a12,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)
  if(.not.list) write(*,1002) nfft,nr,nw,ncomplex,npatience
1002 format(/'nfft: ',i10,'   nr:',i2,'   nw',i2,'   nc:',i2,'   np:',i2/)

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

  open(12,file='chkfft.out',status='unknown')

  if(nr.ne.0) then
     isuccess=fftwf_import_wisdom_from_filename('fftwf_wisdom.dat'//char(0))
     if(isuccess.eq.1) then
        write(*,1010)
1010    format('Imported FFTW wisdom.')
     else
        write(*,1012)
1012    format('Failed to import FFTW wisdom.')
        go to 999
     endif
  endif

  idum=-1                               !Set random seed
  ndim=1                                !One-dimensional transforms
  do i=1,NMAX                           !Set random data
     x=gran()
     y=gran()
     b(i)=cmplx(x,y)                    !Generate random data
  enddo

  iters=1000000
  if(list .or. (nfft.gt.0)) then
     n1=1
     n2=1
     if(nfft.eq.-1) n2=999999
     write(*,1020)
1020 format('    NFFT     Time        rms      MHz   MFlops  iters',    &
          '  tplan'/61('-'))
  else
     n1=4
     n2=23
     write(*,1030)
1030 format(' n   N=2^n     Time        rms      MHz   MFlops  iters',  &
          '  tplan'/63('-'))
  endif

  do ii=n1,n2                           !Test one or more FFT lengths
     if(list) then
        read(10,*,end=900) nfft         !Read nfft from file
     else if(n2.gt.n1) then
        nfft=2**ii                      !Do powers of 2
     endif

     iformf=1
     iformb=1
     if(ncomplex.eq.0) then
        iformf=0                        !Real-to-complex transform
        iformb=-1                       !Complex-to-real (inverse) transform
     endif

     if(nfft.gt.NMAX) go to 900
     a(1:nfft)=b(1:nfft)                !Copy test data into a()
     t0=second()
     if(ncomplex.ne.0) then
        plan1=fftwf_plan_dft_1d(nfft,a,c,-1,nflags)
        plan2=fftwf_plan_dft_1d(nfft,a,c,+1,nflags)
     else
        plan1=fftwf_plan_dft_r2c_1d(nfft,ar,c,nflags)
        plan2=fftwf_plan_dft_c2r_1d(nfft,c,ar,nflags)
     endif

     t2=second()
     tplan=t2-t0                        !Total planning time for this length

     total=0.
     do iter=1,iters                    !Now do many iterations
        a(1:nfft)=b(1:nfft)             !Copy test data into a()
        t0=second()
        call fftwf_execute_dft(plan1,a,c)
        call fftwf_execute_dft(plan2,c,a)
        t1=second()
        total=total+t1-t0
        if(total.ge.1.0) go to 40       !Cut iterations short if t>1 s
     enddo
     iter=iters

40   time=0.5*total/iter                !Time for one FFT of current length
     tplan=0.5*tplan-time               !Planning time for one FFT
     if(tplan.lt.0) tplan=0.
     a(1:nfft)=a(1:nfft)/nfft

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

     freq=1.e-6*nfft/time
     mflops=5.0/(1.e6*time/(nfft*log(float(nfft))/log(2.0)))
     if(n2.eq.1 .or. n2.eq.999999) then
        write(*,1050) nfft,time,rms,freq,mflops,iter,tplan
        write(12,1050) nfft,time,rms,freq,mflops,iter,tplan
1050    format(i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     else
        write(*,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
        write(12,1060) ii,nfft,time,rms,freq,mflops,iter,tplan
1060    format(i2,i8,f11.7,f12.8,f7.2,f8.1,i8,f6.1)
     endif
  enddo

900  continue
  if(nw.eq.1) then
     ierr=fftwf_export_wisdom_to_filename('fftwf_wisdom.dat'//char(0))
     write(*,1070)
1070 format(/'Exported FFTW wisdom')
  endif

  call fftwf_destroy_plan(plan1)
  call fftwf_destroy_plan(plan2)

999 end program chkfft


code:
subroutine chkhist(mrsym,nmax,ipk)

  integer mrsym(63)
  integer hist(0:63)

  hist=0
  do j=1,63
     i=mrsym(j)
     hist(i)=hist(i)+1
  enddo

  nmax=0
  do i=0,63
     if(hist(i).gt.nmax) then
        nmax=hist(i)
        ipk=i+1
     endif
  enddo

  return
end subroutine chkhist


code:
subroutine chkmsg(message,cok,nspecial,flip)

  character message*22,cok*3

  nspecial=0
  flip=1.0
  cok="   "

  do i=22,1,-1
     if(message(i:i).ne.' ') go to 10
  enddo
  i=22

10 if(i.ge.11) then
     if((message(i-3:i).eq.' OOO') .or. (message(20:22).eq.' OO')) then
        cok='OOO'
        flip=-1.0
        if(message(20:22).eq.' OO') then
           message=message(1:19)
        else
           message=message(1:i-4)
        endif
     endif
  endif

  if(message.eq.'RO                    ') nspecial=2
  if(message.eq.'RRR                   ') nspecial=3
  if(message.eq.'73                    ') nspecial=4

  return
end subroutine chkmsg


code:
subroutine chkss2(ss2,freq,drift,schk)

  real ss2(0:8,85)
  real s(0:8,85)
  include 'jt9sync.f90'

  ave=sum(ss2)/(9*85)
  if(freq+drift.eq.-999999.0) ave=0.      !To silence compiler warning
  s=ss2/ave-1.0

  s1=0.
  do i=1,16
     j=ii(i)
     if(j.le.85) s1=s1 + s(0,j)
  enddo
  schk=s1/16.0

  return
end subroutine chkss2


code:
program code426

  parameter (MZ=26)                !Number of 4-FSK symbols
  parameter (JZMAX=64)             !Desired number of codewords
  integer ic(MZ,JZMAX),icsave(MZ)
  real c(MZ)
  character*12 arg

  nargs=iargc()
  if(nargs.ne.2) then
     print*,'Usage:   code426 <nmsgs> <iters>'
     print*,'Example: code426   64   10000000'
     go to 999
  endif
  call getarg(1,arg)
  read(arg,*) nmsgs
  call getarg(2,arg)
  read(arg,*) iters

  call init_random_seed()

  open(13,file='code426.out',status='unknown')

  write(*,1002) nmsgs,iters
  write(13,1002) nmsgs,iters
1002 format('Nmsgs:',i4,'   Iters:',i10/(66('-')))

  do i=1,MZ                     !Create 4 mutually orthogonal codewords
     ic(i,1)=mod(i-1,4)
     ic(i,2)=mod(i,4)
     ic(i,3)=mod(i+1,4)
     ic(i,4)=mod(i+2,4)
  enddo

  do j=1,4                      !Write them out
     write(*,1000) j,MZ,ic(1:MZ,j)
     write(13,1000) j,MZ,ic(1:MZ,j)
1000 format(2i5,3x,26i2)
  enddo

  do j=5,nmsgs                  !Find codewords up to j=nmsgs with maximum
     npk=0                      !distance from all the rest
     do i=1,iters
        call random_number(c)   !Generate a random codeword candidate
        ic(1:MZ,j)=int(4*c)     !Convert real to integer
!        nd=MZ
!        do k=1,j-1              !Test candidate against all others in list
!           n=count(ic(1:MZ,j).ne.ic(1:MZ,k))
!           nd=min(n,nd)
!        enddo
        call dist426(ic,j,mind)
        if(mind.gt.npk) then
           npk=mind
           icsave=ic(1:MZ,j)    !Best candidate so far, save it
!           if(npk.ge.19) exit   !It won't get any better...
        endif
     enddo
     write(*,1000) j,npk,ic(1:MZ,j)
     write(13,1000) j,npk,ic(1:MZ,j)
  enddo

999 end program code426


code:
  integer, parameter :: NTMAX=300
  integer, parameter :: NMAX=NTMAX*12000 !Total sample intervals (one minute)
  integer, parameter :: NDMAX=NTMAX*1500 !Sample intervals at 1500 Hz rate
  integer, parameter :: NSMAX=6827       !Max length of saved spectra
  integer, parameter :: MAXFFT3=16384


code:
program contest72

  use packjt
  integer dat(12)
  logical text,bcontest,ok
  character*22 msg,msg0,msg1
  character*72 ct1,ct2
  character*12 callsign1,callsign2
  character*1 c0
  character*42 c
  character*6 mygrid
  data c/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ +-./?'/
  data bcontest/.true./
  data mygrid/"EM48  "/

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

  nargs=iargc()
  if(nargs.eq.0) open(10,file='contest_msgs.txt',status='old')

  nn=0
  do imsg=1,9999
     if(nargs.eq.1) then
        if(imsg.gt.1) exit
        call getarg(1,msg0)
     else
        read(10,1001,end=999) msg0
1001    format(a22)
     endif
     msg=msg0
     call packmsg(msg,dat,itype,bcontest)
     call unpackmsg(dat,msg1,bcontest,mygrid)
     ok=msg1.eq.msg0
     if(msg0.eq.'                      ') then
        write(*,1002)
     else
        if(jt_c2(1:1).eq.'W') msg0='  '//msg0(1:20)
        nn=nn+1
        write(*,1002) nn,msg0,ok,jt_itype,jt_nc1,jt_nc2,jt_ng,jt_k1,jt_k2
1002    format(i1,'. ',a22,L2,i2,2i10,i6,2i8)
        if(index(msg1,' 73 ').gt.4) nn=0
     endif
     if(.not.ok) print*,msg0,msg1
     if(itype.lt.0 .or. itype.eq.6) cycle

     if(msg(1:3).eq.'CQ ') then
        m=2
        write(ct1,1010) dat
1010    format(12b6.6)
!        write(*,1014) ct1
!1014    format(a72)
        cycle
     endif

     i1=index(msg,'<')
     if(i1.eq.1) then
        m=0
        cycle
     endif

     if(i.ge.5) then
        m=3
        cycle
     endif

     if(msg(1:6).eq.'73 CQ ') then
        m=4
        cycle
     endif

     call packmsg(msg,dat,itype,.false.)
     write(ct1,1010) dat
     call packtext(msg,nc1,nc2,ng,.false.,'')
!     write(ct2,1012) nc1,nc2,ng+32768
!1012 format(2b28.28,b16.16)
!     write(*,1014) ct1
!     write(*,1014) ct2
!     write(*,1014)
  enddo

999 end program contest72


code:
SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2)

! Examples:
! 1. From ha,dec to az,el:
!      call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
! 2. From az,el to ha,dec:
!      call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
! 3. From ra,dec to l,b
!      call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
!        ra,dec,l,b)
! 4. From l,b to ra,dec
!      call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
!        0.478220215d0,l,b,ra,dec)
! 5. From ra,dec to ecliptic latitude (eb) and longitude (el):
!      call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
! 6. From ecliptic latitude (eb) and longitude (el) to ra,dec:
!      call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec)

  SB0=sin(B0)
  CB0=cos(B0)
  SBP=sin(BP)
  CBP=cos(BP)
  SB1=sin(B1)
  CB1=cos(B1)
  SB2=SBP*SB1 + CBP*CB1*cos(AP-A1)
  CB2=SQRT(1.e0-SB2**2)
  B2=atan(SB2/CB2)
  SAA=sin(AP-A1)*CB1/CB2
  CAA=(SB1-SB2*SBP)/(CB2*CBP)
  CBB=SB0/CBP
  SBB=sin(AP-A0)*CB0
  SA2=SAA*CBB-CAA*SBB
  CA2=CAA*CBB+SAA*SBB
  TA2O2=0.0 !Shut up compiler warnings. -db
  IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2
  IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2)
  A2=2.e0*atan(TA2O2)
  IF(A2.LT.0.e0) A2=A2+6.2831853

  RETURN
END SUBROUTINE COORD


code:
program count4

  parameter(NMAX=1000)
  character*47 line
  real snr(NMAX)
  real dt(NMAX)
  real f(NMAX)

  open(10,file='/users/joe/appdata/local/wsjt-x/all.txt',status='old')

  read(10,1000,end=10) line
1000 format(a47)

  nsync1=0
  nsync2=0
  n1=0
  n2=0
  nerr=0

  do i=1,99999
     read(10,1000,end=10) line
     if(line(47:47).ne.' ') cycle                !Skip average decodes
     if(line(20:20).eq.'*') nsync1=nsync1+1
     if(line(20:20).eq.'#') nsync2=nsync2+1
     if(line(22:34).eq.'CQ K1ABC FN42') then
        n2=n2+1                                  !Correlation decode
        read(line,1002) snr(n2),dt(n2),f(n2)
1002    format(4x,f4.0,f5.2,f5.0)
        if(line(42:42).eq.'*') n1=n1+1           !Convolutional decode
     else
        if(line(22:34).ne.'             ') nerr=nerr+1
     endif
  enddo

10 call stats(snr,n2,snrave,snrdev)
  call stats(dt,n2,dtave,dtdev)
  call stats(f,n2,fave,fdev)

write(*,1010) nsync1,nsync2,n1,n2,nerr,snrave,dtave,fave,snrdev,dtdev,fdev
1010 format(5i5,f7.1,f7.2,f7.0/25x,f7.1,f7.2,f7.0)

end program count4

subroutine stats(x,nz,ave,rms)
  real x(nz)

  ave=0.
  rms=0.
  if(nz.gt.0) ave=sum(x)/nz
  x=x-ave
  if(nz.gt.1) rms=sqrt(dot_product(x,x)/(nz-1))

  return
end subroutine stats