w4kek/WSJT-X
1
0
Fork 0
mirror of https://github.com/saitohirga/WSJT-X.git synced 2025-03-22 20:18:35 -04:00
Code Issues Packages Projects Releases Wiki Activity

Prune further...

Browse Source
This commit is contained in:
Joe Taylor 2022-12-12 11:35:23 -05:00
parent ecaa0b8861
commit 39318e9d21
20 changed files with 0 additions and 2082 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
q65w/libm65/CMakeLists.txt
View File

@ -41,7 +41,6 @@ set (libm65_FSRCS
recvpkt.f90
rfile3a.f90
s3avg.f90
sec_midn.f90
set.f90
shell.f90
sleep_msec.f90
@ -59,36 +58,10 @@ set (libm65_FSRCS
f77_wisdom.f
)
set (libm65_ka9q_CSRCS
decode_rs.c
encode_rs.c
init_rs.c
)
set_source_files_properties (${libm65_ka9q_CSRCS} PROPERTIES COMPILE_FLAGS -Wno-sign-compare)
set (libm65_CSRCS
${libm65_ka9q_CSRCS}
ftrsd2.c
# gran.c
igray.c
tmoonsub.c
usleep.c
wrapkarn.c
)
if (WIN32)
set (libm65_CSRCS ${libm65_CSRCS} ptt.c)
else ()
set (libm65_CSRCS ${libm65_CSRCS} ptt_unix.c)
endif ()
set (libm65_CXXSRCS
ipcomm.cpp
)
add_definitions (-DBIGSYM=1)
set_source_files_properties (sec_midn.f90 PROPERTIES COMPILE_FLAGS -fno-second-underscore)
set (libm65_C_and_CXXSRCS
${libm65_CSRCS}
${libm65_CXXSRCS}

93
q65w/libm65/cutil.c
View File

@ -1,93 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
// #include <sys/times.h>
// #include <time.h>
// #include <sys/time.h>
#include "sleep.h"
#include "timeval.h"
/* FORTRAN: fd = close(filedes) */
int close_(int *filedes)
{
return(close(*filedes));
}
/* FORTRAN: fd = open(filnam,mode) */
int open_(char filnam[], int *mode)
{
return(open(filnam,*mode));
}
/* FORTRAN: fd = creat(filnam,mode) */
int creat_(char filnam[],int *mode)
{
return(creat(filnam,*mode));
}
/* FORTRAN: nread = read(fd,buf,n) */
int read_(int *fd, char buf[], int *n)
{
return(read(*fd,buf,*n));
}
/* FORTRAN: nwrt = write(fd,buf,n) */
int write_(int *fd, char buf[], int *n)
{
return(write(*fd,buf,*n));
}
/* FORTRAN: ns = lseek(fd,offset,origin) */
int lseek_(int *fd,int *offset, int *origin)
{
return(lseek(*fd,*offset,*origin));
}
/* times(2) */
//int times_(struct tms *buf)
//{
// return (times(buf));
//}
/* ioperm(2) */
//ioperm_(from,num,turn_on)
//unsigned long *from,*num,*turn_on;
//{
// return (ioperm(*from,*num,*turn_on));
// return (i386_get_ioperm(*from,*num,*turn_on));
//}
/* usleep(3) */
void usleep_(unsigned long *microsec)
{
usleep(*microsec);
}
/* returns random numbers between 0 and 32767 to FORTRAN program */
int iran_(int *arg)
{
return (rand());
}
int exit_(int *n)
{
printf("\n\n");
exit(*n);
}
/*
struct tm *
gmtime_r_(const time_t *clock, struct tm *result)
{
gmtime_r(clock, result);
}
*/
time_t time_(void)
{
return time(0);
}
/* hrtime() */
double hrtime_(void)
{
struct timeval tv;
gettimeofday(&tv,NULL);
return(tv.tv_sec+1.e-6*tv.tv_usec);
}

268
q65w/libm65/decode_rs.c
View File

@ -1,268 +0,0 @@
/* Reed-Solomon decoder
* Copyright 2002 Phil Karn, KA9Q
* May be used under the terms of the GNU General Public License (GPL)
* Modified by Steve Franke, K9AN, for use in a soft-symbol RS decoder
*/
#ifdef DEBUG
#include <stdio.h>
#endif
#include <string.h>
#define NULL ((void *)0)
#define min(a,b) ((a) < (b) ? (a) : (b))
#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif
int DECODE_RS(
#ifndef FIXED
void *p,
#endif
DTYPE *data, int *eras_pos, int no_eras, int calc_syn){
#ifndef FIXED
struct rs *rs = (struct rs *)p;
#endif
int deg_lambda, el, deg_omega;
int i, j, r,k;
DTYPE u,q,tmp,num1,num2,den,discr_r;
DTYPE lambda[NROOTS+1]; // Err+Eras Locator poly
static DTYPE s[51]; // and syndrome poly
DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
int syn_error, count;
if( calc_syn ) {
/* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
for(i=0;i<NROOTS;i++)
s[i] = data[0];
for(j=1;j<NN;j++){
for(i=0;i<NROOTS;i++){
if(s[i] == 0){
s[i] = data[j];
} else {
s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
}
}
}
/* Convert syndromes to index form, checking for nonzero condition */
syn_error = 0;
for(i=0;i<NROOTS;i++){
syn_error |= s[i];
s[i] = INDEX_OF[s[i]];
}
if (!syn_error) {
/* if syndrome is zero, data[] is a codeword and there are no
* errors to correct. So return data[] unmodified
*/
count = 0;
goto finish;
}
}
memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
lambda[0] = 1;
if (no_eras > 0) {
/* Init lambda to be the erasure locator polynomial */
lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
for (i = 1; i < no_eras; i++) {
u = MODNN(PRIM*(NN-1-eras_pos[i]));
for (j = i+1; j > 0; j--) {
tmp = INDEX_OF[lambda[j - 1]];
if(tmp != A0)
lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
}
}
#if DEBUG >= 1
/* Test code that verifies the erasure locator polynomial just constructed
Needed only for decoder debugging. */
/* find roots of the erasure location polynomial */
for(i=1;i<=no_eras;i++)
reg[i] = INDEX_OF[lambda[i]];
count = 0;
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
q = 1;
for (j = 1; j <= no_eras; j++)
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
if (q != 0)
continue;
/* store root and error location number indices */
root[count] = i;
loc[count] = k;
count++;
}
if (count != no_eras) {
printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
count = -1;
goto finish;
}
#if DEBUG >= 2
printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
for (i = 0; i < count; i++)
printf("%d ", loc[i]);
printf("\n");
#endif
#endif
}
for(i=0;i<NROOTS+1;i++)
b[i] = INDEX_OF[lambda[i]];
/*
* Begin Berlekamp-Massey algorithm to determine error+erasure
* locator polynomial
*/
r = no_eras;
el = no_eras;
while (++r <= NROOTS) { /* r is the step number */
/* Compute discrepancy at the r-th step in poly-form */
discr_r = 0;
for (i = 0; i < r; i++){
if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
}
}
discr_r = INDEX_OF[discr_r]; /* Index form */
if (discr_r == A0) {
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1],b,NROOTS*sizeof(b[0]));
b[0] = A0;
} else {
/* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
t[0] = lambda[0];
for (i = 0 ; i < NROOTS; i++) {
if(b[i] != A0)
t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
else
t[i+1] = lambda[i+1];
}
if (2 * el <= r + no_eras - 1) {
el = r + no_eras - el;
/*
* 2 lines below: B(x) <-- inv(discr_r) *
* lambda(x)
*/
for (i = 0; i <= NROOTS; i++)
b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
} else {
/* 2 lines below: B(x) <-- x*B(x) */
memmove(&b[1],b,NROOTS*sizeof(b[0]));
b[0] = A0;
}
memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
}
}
/* Convert lambda to index form and compute deg(lambda(x)) */
deg_lambda = 0;
for(i=0;i<NROOTS+1;i++){
lambda[i] = INDEX_OF[lambda[i]];
if(lambda[i] != A0)
deg_lambda = i;
}
/* Find roots of the error+erasure locator polynomial by Chien search */
memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
count = 0; /* Number of roots of lambda(x) */
for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
q = 1; /* lambda[0] is always 0 */
for (j = deg_lambda; j > 0; j--){
if (reg[j] != A0) {
reg[j] = MODNN(reg[j] + j);
q ^= ALPHA_TO[reg[j]];
}
}
if (q != 0)
continue; /* Not a root */
/* store root (index-form) and error location number */
#if DEBUG>=2
printf("count %d root %d loc %d\n",count,i,k);
#endif
root[count] = i;
loc[count] = k;
/* If we've already found max possible roots,
* abort the search to save time
*/
if(++count == deg_lambda)
break;
}
if (deg_lambda != count) {
/*
* deg(lambda) unequal to number of roots => uncorrectable
* error detected
*/
count = -1;
goto finish;
}
/*
* Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
* x**NROOTS). in index form. Also find deg(omega).
*/
deg_omega = 0;
for (i = 0; i < NROOTS;i++){
tmp = 0;
j = (deg_lambda < i) ? deg_lambda : i;
for(;j >= 0; j--){
if ((s[i - j] != A0) && (lambda[j] != A0))
tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
}
if(tmp != 0)
deg_omega = i;
omega[i] = INDEX_OF[tmp];
}
omega[NROOTS] = A0;
/*
* Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
* inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
*/
for (j = count-1; j >=0; j--) {
num1 = 0;
for (i = deg_omega; i >= 0; i--) {
if (omega[i] != A0)
num1 ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
}
num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
den = 0;
/* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
if(lambda[i+1] != A0)
den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
}
if (den == 0) {
#if DEBUG >= 1
printf("\n ERROR: denominator = 0\n");
#endif
count = -1;
goto finish;
}
/* Apply error to data */
if (num1 != 0) {
data[loc[j]] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
}
}
finish:
if(eras_pos != NULL){
for(i=0;i<count;i++)
eras_pos[i] = loc[i];
}
return count;
}

47
q65w/libm65/encode_rs.c
View File

@ -1,47 +0,0 @@
/* Reed-Solomon encoder
* Copyright 2002, Phil Karn, KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
#include <string.h>
#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif
void ENCODE_RS(
#ifndef FIXED
void *p,
#endif
DTYPE *data, DTYPE *bb){
#ifndef FIXED
struct rs *rs = (struct rs *)p;
#endif
int i, j;
DTYPE feedback;
memset(bb,0,NROOTS*sizeof(DTYPE));
for(i=0;i<NN-NROOTS;i++){
feedback = INDEX_OF[data[i] ^ bb[0]];
if(feedback != A0){ /* feedback term is non-zero */
#ifdef UNNORMALIZED
/* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must
* always be for the polynomials constructed by init_rs()
*/
feedback = MODNN(NN - GENPOLY[NROOTS] + feedback);
#endif
for(j=1;j<NROOTS;j++)
bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
}
/* Shift */
memmove(&bb[0],&bb[1],sizeof(DTYPE)*(NROOTS-1));
if(feedback != A0)
bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
else
bb[NROOTS-1] = 0;
}
}

213
q65w/libm65/ftrsd2.c
View File

@ -1,213 +0,0 @@
/*
ftrsd2.c
A soft-decision decoder for the JT65 (63,12) Reed-Solomon code.
This decoding scheme is built around Phil Karn's Berlekamp-Massey
errors and erasures decoder. The approach is inspired by a number of
publications, including the stochastic Chase decoder described
in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
"Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.
Steve Franke K9AN and Joe Taylor K1JT
*/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <string.h>
#include "rs2.h"
static void *rs;
void getpp_(int workdat[], float *pp);
void ftrsd2_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[],
int* ntrials0, int correct[], int param[], int ntry[])
{
int rxdat[63], rxprob[63], rxdat2[63], rxprob2[63];
int workdat[63];
int indexes[63];
int era_pos[51];
int i, j, numera, nerr, nn=63;
int ntrials = *ntrials0;
int nhard=0,nhard_min=32768,nsoft=0,nsoft_min=32768;
int ntotal=0,ntotal_min=32768,ncandidates;
int nera_best=0;
float pp,pp1,pp2;
static unsigned int nseed;
// Power-percentage symbol metrics - composite gnnf/hf
int perr[8][8] = {
{ 4, 9, 11, 13, 14, 14, 15, 15},
{ 2, 20, 20, 30, 40, 50, 50, 50},
{ 7, 24, 27, 40, 50, 50, 50, 50},
{13, 25, 35, 46, 52, 70, 50, 50},
{17, 30, 42, 54, 55, 64, 71, 70},
{25, 39, 48, 57, 64, 66, 77, 77},
{32, 45, 54, 63, 66, 75, 78, 83},
{51, 58, 57, 66, 72, 77, 82, 86}};
// Initialize the KA9Q Reed-Solomon encoder/decoder
unsigned int symsize=6, gfpoly=0x43, fcr=3, prim=1, nroots=51;
rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 0);
// Reverse the received symbol vectors for BM decoder
for (i=0; i<63; i++) {
rxdat[i]=mrsym[62-i];
rxprob[i]=mrprob[62-i];
rxdat2[i]=mr2sym[62-i];
rxprob2[i]=mr2prob[62-i];
}
// Sort rxprob to find indexes of the least reliable symbols
int k, pass, tmp, nsym=63;
int probs[63];
for (i=0; i<63; i++) {
indexes[i]=i;
probs[i]=rxprob[i];
}
for (pass = 1; pass <= nsym-1; pass++) {
for (k = 0; k < nsym - pass; k++) {
if( probs[k] < probs[k+1] ) {
tmp = probs[k];
probs[k] = probs[k+1];
probs[k+1] = tmp;
tmp = indexes[k];
indexes[k] = indexes[k+1];
indexes[k+1] = tmp;
}
}
}
// See if we can decode using BM HDD, and calculate the syndrome vector.
memset(era_pos,0,51*sizeof(int));
numera=0;
memcpy(workdat,rxdat,sizeof(rxdat));
nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
if( nerr >= 0 ) {
// Hard-decision decoding succeeded. Save codeword and some parameters.
nhard=0;
for (i=0; i<63; i++) {
if( workdat[i] != rxdat[i] ) nhard=nhard+1;
}
memcpy(correct,workdat,63*sizeof(int));
param[0]=0;
param[1]=nhard;
param[2]=0;
param[3]=0;
param[4]=0;
param[5]=0;
param[7]=1000*1000;
ntry[0]=0;
return;
}
/*
Hard-decision decoding failed. Try the FT soft-decision method.
Generate random erasure-locator vectors and see if any of them
decode. This will generate a list of "candidate" codewords. The
soft distance between each candidate codeword and the received
word is estimated by finding the largest (pp1) and second-largest
(pp2) outputs from a synchronized filter-bank operating on the
symbol spectra, and using these to decide which candidate
codeword is "best".
*/
nseed=1; //Seed for random numbers
float ratio;
int thresh, nsum;
int thresh0[63];
ncandidates=0;
nsum=0;
int ii,jj;
for (i=0; i<nn; i++) {
nsum=nsum+rxprob[i];
j = indexes[62-i];
ratio = (float)rxprob2[j]/((float)rxprob[j]+0.01);
ii = 7.999*ratio;
jj = (62-i)/8;
thresh0[i] = 1.3*perr[ii][jj];
}
if(nsum<=0) return;
pp1=0.0;
pp2=0.0;
for (k=1; k<=ntrials; k++) {
memset(era_pos,0,51*sizeof(int));
memcpy(workdat,rxdat,sizeof(rxdat));
/*
Mark a subset of the symbols as erasures.
Run through the ranked symbols, starting with the worst, i=0.
NB: j is the symbol-vector index of the symbol with rank i.
*/
numera=0;
for (i=0; i<nn; i++) {
j = indexes[62-i];
thresh=thresh0[i];
long int ir;
// Generate a random number ir, 0 <= ir < 100 (see POSIX.1-2001 example).
nseed = nseed * 1103515245 + 12345;
ir = (unsigned)(nseed/65536) % 32768;
ir = (100*ir)/32768;
if((ir < thresh ) && numera < 51) {
era_pos[numera]=j;
numera=numera+1;
}
}
nerr=decode_rs_int(rs,workdat,era_pos,numera,0);
if( nerr >= 0 ) {
// We have a candidate codeword. Find its hard and soft distance from
// the received word. Also find pp1 and pp2 from the full array
// s3(64,63) of synchronized symbol spectra.
ncandidates=ncandidates+1;
nhard=0;
nsoft=0;
for (i=0; i<63; i++) {
if(workdat[i] != rxdat[i]) {
nhard=nhard+1;
if(workdat[i] != rxdat2[i]) {
nsoft=nsoft+rxprob[i];
}
}
}
nsoft=63*nsoft/nsum;
ntotal=nsoft+nhard;
getpp_(workdat,&pp);
if(pp>pp1) {
pp2=pp1;
pp1=pp;
nsoft_min=nsoft;
nhard_min=nhard;
ntotal_min=ntotal;
memcpy(correct,workdat,63*sizeof(int));
nera_best=numera;
ntry[0]=k;
} else {
if(pp>pp2 && pp!=pp1) pp2=pp;
}
if(nhard_min <= 41 && ntotal_min <= 71) break;
}
if(k == ntrials) ntry[0]=k;
}
param[0]=ncandidates;
param[1]=nhard_min;
param[2]=nsoft_min;
param[3]=nera_best;
param[4]= pp1 > 0 ? 1000.0*pp2/pp1 : 1000.0;
param[5]=ntotal_min;
param[6]=ntry[0];
param[7]=1000.0*pp2;
param[8]=1000.0*pp1;
if(param[0]==0) param[2]=-1;
return;
}

28
q65w/libm65/gran.c
View File

@ -1,28 +0,0 @@
#include <stdlib.h>
#include <math.h>
/* Generate gaussian random float with mean=0 and std_dev=1 */
float gran_()
{
float fac,rsq,v1,v2;
static float gset;
static int iset;
if(iset){
/* Already got one */
iset = 0;
return gset;
}
/* Generate two evenly distributed numbers between -1 and +1
* that are inside the unit circle
*/
do {
v1 = 2.0 * (float)rand() / RAND_MAX - 1;
v2 = 2.0 * (float)rand() / RAND_MAX - 1;
rsq = v1*v1 + v2*v2;
} while(rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0*log(rsq)/rsq);
gset = v1*fac;
iset++;
return v2*fac;
}

22
q65w/libm65/igray.c
View File

@ -1,22 +0,0 @@
#ifdef CVF
extern int __stdcall IGRAY(int *n0, int *idir)
#else
int igray_(int *n0, int *idir)
#endif
{
int n;
unsigned long sh;
unsigned long nn;
n=*n0;
if(*idir>0) return (n ^ (n >> 1));
sh = 1;
nn = (n >> sh);
while (nn > 0) {
n ^= nn;
sh <<= 1;
nn = (n >> sh);
}
return (n);
}

120
q65w/libm65/init_rs.c
View File

@ -1,120 +0,0 @@
/* Initialize a RS codec
*
* Copyright 2002 Phil Karn, KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
#include <stdlib.h>
#ifdef CCSDS
#include "ccsds.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif
void FREE_RS(void *p){
struct rs *rs = (struct rs *)p;
free(rs->alpha_to);
free(rs->index_of);
free(rs->genpoly);
free(rs);
}
/* Initialize a Reed-Solomon codec
* symsize = symbol size, bits (1-8)
* gfpoly = Field generator polynomial coefficients
* fcr = first root of RS code generator polynomial, index form
* prim = primitive element to generate polynomial roots
* nroots = RS code generator polynomial degree (number of roots)
*/
void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned fcr,unsigned prim,
unsigned int nroots){
struct rs *rs;
int i, j, sr,root,iprim;
/* Check parameter ranges */
if(symsize < 0 || symsize > (int)(8*sizeof(DTYPE)))
return NULL; /* Need version with ints rather than chars */
if(fcr >= (1<<symsize))
return NULL;
if(prim == 0 || prim >= (1<<symsize))
return NULL;
if(nroots >= (1<<symsize))
return NULL; /* Can't have more roots than symbol values! */
rs = (struct rs *)calloc(1,sizeof(struct rs));
rs->mm = symsize;
rs->nn = (1<<symsize)-1;
rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
if(rs->alpha_to == NULL){
free(rs);
return NULL;
}
rs->index_of = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
if(rs->index_of == NULL){
free(rs->alpha_to);
free(rs);
return NULL;
}
/* Generate Galois field lookup tables */
rs->index_of[0] = A0; /* log(zero) = -inf */
rs->alpha_to[A0] = 0; /* alpha**-inf = 0 */
sr = 1;
for(i=0;i<rs->nn;i++){
rs->index_of[sr] = i;
rs->alpha_to[i] = sr;
sr <<= 1;
if(sr & (1<<symsize))
sr ^= gfpoly;
sr &= rs->nn;
}
if(sr != 1){
/* field generator polynomial is not primitive! */
free(rs->alpha_to);
free(rs->index_of);
free(rs);
return NULL;
}
/* Form RS code generator polynomial from its roots */
rs->genpoly = (DTYPE *)malloc(sizeof(DTYPE)*(nroots+1));
if(rs->genpoly == NULL){
free(rs->alpha_to);
free(rs->index_of);
free(rs);
return NULL;
}
rs->fcr = fcr;
rs->prim = prim;
rs->nroots = nroots;
/* Find prim-th root of 1, used in decoding */
for(iprim=1;(iprim % prim) != 0;iprim += rs->nn)
;
rs->iprim = iprim / prim;
rs->genpoly[0] = 1;
for (i = 0,root=fcr*prim; i < nroots; i++,root += prim) {
rs->genpoly[i+1] = 1;
/* Multiply rs->genpoly[] by @**(root + x) */
for (j = i; j > 0; j--){
if (rs->genpoly[j] != 0)
rs->genpoly[j] = rs->genpoly[j-1] ^ rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[j]] + root)];
else
rs->genpoly[j] = rs->genpoly[j-1];
}
/* rs->genpoly[0] can never be zero */
rs->genpoly[0] = rs->alpha_to[modnn(rs,rs->index_of[rs->genpoly[0]] + root)];
}
/* convert rs->genpoly[] to index form for quicker encoding */
for (i = 0; i <= nroots; i++)
rs->genpoly[i] = rs->index_of[rs->genpoly[i]];
return rs;
}

54
q65w/libm65/int.h
View File

@ -1,54 +0,0 @@
/* Include file to configure the RS codec for integer symbols
*
* Copyright 2002, Phil Karn, KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
#define DTYPE int
/* Reed-Solomon codec control block */
struct rs {
unsigned int mm; /* Bits per symbol */
unsigned int nn; /* Symbols per block (= (1<<mm)-1) */
int *alpha_to; /* log lookup table */
int *index_of; /* Antilog lookup table */
int *genpoly; /* Generator polynomial */
unsigned int nroots; /* Number of generator roots = number of parity symbols */
unsigned int fcr; /* First consecutive root, index form */
unsigned int prim; /* Primitive element, index form */
unsigned int iprim; /* prim-th root of 1, index form */
};
static inline int modnn(struct rs *rs,int x){
while (x >= rs->nn) {
x -= rs->nn;
x = (x >> rs->mm) + (x & rs->nn);
}
return x;
}
#define MODNN(x) modnn(rs,x)
#define MM (rs->mm)
#define NN (rs->nn)
#define ALPHA_TO (rs->alpha_to)
#define INDEX_OF (rs->index_of)
#define GENPOLY (rs->genpoly)
#define NROOTS (rs->nroots)
#define FCR (rs->fcr)
#define PRIM (rs->prim)
#define IPRIM (rs->iprim)
#define A0 (NN)
#define ENCODE_RS encode_rs_int
#define DECODE_RS decode_rs_int
#define INIT_RS init_rs_int
#define FREE_RS free_rs_int
void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras, int calc_syn);
void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned int fcr,
unsigned int prim,unsigned int nroots);
void FREE_RS(void *p);

1
q65w/libm65/map65a.f90
View File

@ -121,7 +121,6 @@ subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb, &
short=0. !Zero the whole short array
jpz=1
print*,'AAA',mode65
call timer('filbig ',0)
call filbig(dd,NSMAX,f0,newdat,nfsample,xpol,cx,cy,n5)
call timer('filbig ',1)

43
q65w/libm65/ptt.c
View File

@ -1,43 +0,0 @@
#include <windows.h>
#include <stdio.h>
int ptt_(int *nport, int *ntx, int *iptt)
{
static HANDLE hFile;
static int open=0;
char s[10];
int i3=0,i4=0,i5=0,i6=0,i9=0,i00=0;
if(*nport==0) {
*iptt=*ntx;
return(0);
}
if(*ntx && (!open)) {
sprintf(s,"COM%d",*nport);
hFile=CreateFile(TEXT(s),GENERIC_WRITE,0,NULL,OPEN_EXISTING,
FILE_ATTRIBUTE_NORMAL,NULL);
if(hFile==INVALID_HANDLE_VALUE) {
// printf("PTT: Cannot open COM port %d.\n",*nport);
return 1;
}
open=1;
}
if(*ntx && open) {
EscapeCommFunction(hFile,3);
EscapeCommFunction(hFile,5);
*iptt=1;
}
else {
EscapeCommFunction(hFile,4);
EscapeCommFunction(hFile,6);
EscapeCommFunction(hFile,9);
i00=CloseHandle(hFile);
*iptt=0;
open=0;
}
if((i00+i3+i4+i5+i6+i9)==-99) return -1; //Silence compiler warning
return 0;
}

405
q65w/libm65/ptt_unix.c
View File

@ -1,405 +0,0 @@
/*
* WSJT is Copyright (c) 2001-2006 by Joseph H. Taylor, Jr., K1JT,
* and is licensed under the GNU General Public License (GPL).
*
* Code used from cwdaemon for parallel port ptt only.
*
* cwdaemon - morse sounding daemon for the parallel or serial port
* Copyright (C) 2002 -2005 Joop Stakenborg <pg4i@amsat.org>
* and many authors, see the AUTHORS file.
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Library General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
# if HAVE_STDIO_H
# include <stdio.h>
#endif
#if STDC_HEADERS
# include <stdlib.h>
# include <stddef.h>
#else
# if HAVE_STDLIB_H
# include <stdlib.h>
# endif
#endif
#if HAVE_UNISTD_H
# include <unistd.h>
#endif
#if HAVE_SYS_STAT_H
# include <sys/stat.h>
#endif
#if HAVE_SYS_IOCTL_H
# include <sys/ioctl.h>
#endif
#if HAVE_FCNTL_H
# include <fcntl.h>
#endif
#ifdef HAVE_LINUX_PPDEV_H
# include <linux/ppdev.h>
# include <linux/parport.h>
#endif
#ifdef HAVE_DEV_PPBUS_PPI_H
# include <dev/ppbus/ppi.h>
# include <dev/ppbus/ppbconf.h>
#endif
int lp_reset (int fd);
int lp_ptt (int fd, int onoff);
#ifdef HAVE_SYS_STAT_H
# include <sys/stat.h>
#endif
#if (defined(__unix__) || defined(unix)) && !defined(USG)
# include <sys/param.h>
#endif
#include <string.h>
/* parport functions */
int dev_is_parport(int fd);
int ptt_parallel(int fd, int *ntx, int *iptt);
int ptt_serial(int fd, int *ntx, int *iptt);
int fd=-1; /* Used for both serial and parallel */
/*
* ptt_
*
* generic unix PTT routine called from Fortran
*
* Inputs
* unused Unused, to satisfy old windows calling convention
* ptt_port device name serial or parallel
* ntx pointer to fortran command on or off
* iptt pointer to fortran command status on or off
* Returns - non 0 if error
*/
/* Tiny state machine */
#define STATE_PORT_CLOSED 0
#define STATE_PORT_OPEN_PARALLEL 1
#define STATE_PORT_OPEN_SERIAL 2
//int ptt_(int *unused, char *ptt_port, int *ntx, int *iptt)
int ptt_(int *unused, int *ntx, int *iptt)
{
static int state=0;
char *p;
// ### Temporary:
char* ptt_port;
if(*unused != -99) {
*iptt=*ntx;
return 0;
}
// ###
/* In the very unlikely event of a NULL pointer, just return.
* Yes, I realise this should not be possible in WSJT.
*/
if (ptt_port == NULL) {
*iptt = *ntx;
return (0);
}
switch (state) {
case STATE_PORT_CLOSED:
/* Remove trailing ' ' */
if ((p = strchr(ptt_port, ' ')) != NULL)
*p = '\0';
/* If all that is left is a '\0' then also just return */
if (*ptt_port == '\0') {
*iptt = *ntx;
return(0);
}
if ((fd = open(ptt_port, O_RDWR|O_NONBLOCK)) < 0) {
fprintf(stderr, "Can't open %s.\n", ptt_port);
return (1);
}
if (dev_is_parport(fd)) {
state = STATE_PORT_OPEN_PARALLEL;
lp_reset(fd);
ptt_parallel(fd, ntx, iptt);
} else {
state = STATE_PORT_OPEN_SERIAL;
ptt_serial(fd, ntx, iptt);
}
break;
case STATE_PORT_OPEN_PARALLEL:
ptt_parallel(fd, ntx, iptt);
break;
case STATE_PORT_OPEN_SERIAL:
ptt_serial(fd, ntx, iptt);
break;
default:
close(fd);
fd = -1;
state = STATE_PORT_CLOSED;
break;
}
return(0);
}
/*
* ptt_serial
*
* generic serial unix PTT routine called indirectly from Fortran
*
* fd - already opened file descriptor
* ntx - pointer to fortran command on or off
* iptt - pointer to fortran command status on or off
*/
int
ptt_serial(int fd, int *ntx, int *iptt)
{
int control = TIOCM_RTS | TIOCM_DTR;
#if defined (TIOCMBIS) && defined (TIOCMBIS)
if(*ntx) {
ioctl(fd, TIOCMBIS, &control); /* Set DTR and RTS */
*iptt = 1;
} else {
ioctl(fd, TIOCMBIC, &control);
*iptt = 0;
}
#else
unsigned y;
ioctl(fd, TIOCMGET, &y);
if (*ntx) {
y |= control;
} else {
y &= ~control;
}
ioctl(fd, TIOCMSET, &y);
#endif
return(0);
}
/* parport functions */
/*
* dev_is_parport(fd):
*
* inputs - Already open fd
* output - 1 if parallel port, 0 if not
* side effects - Unfortunately, this is platform specific.
*/
#if defined(HAVE_LINUX_PPDEV_H) /* Linux (ppdev) */
int
dev_is_parport(int fd)
{
struct stat st;
int m;
if ((fstat(fd, &st) == -1) ||
((st.st_mode & S_IFMT) != S_IFCHR) ||
(ioctl(fd, PPGETMODE, &m) == -1))
return(0);
return(1);
}
#elif defined(HAVE_DEV_PPBUS_PPI_H) /* FreeBSD (ppbus/ppi) */
int
dev_is_parport(int fd)
{
struct stat st;
unsigned char c;
if ((fstat(fd, &st) == -1) ||
((st.st_mode & S_IFMT) != S_IFCHR) ||
(ioctl(fd, PPISSTATUS, &c) == -1))
return(0);
return(1);
}
#else /* Fallback (nothing) */
int
dev_is_parport(int fd)
{
return(0);
}
#endif
/* Linux wrapper around PPFCONTROL */
#ifdef HAVE_LINUX_PPDEV_H
static void
parport_control (int fd, unsigned char controlbits, int values)
{
struct ppdev_frob_struct frob;
frob.mask = controlbits;
frob.val = values;
if (ioctl (fd, PPFCONTROL, &frob) == -1)
{
fprintf(stderr, "Parallel port PPFCONTROL");
exit (1);
}
}
#endif
/* FreeBSD wrapper around PPISCTRL */
#ifdef HAVE_DEV_PPBUS_PPI_H
static void
parport_control (int fd, unsigned char controlbits, int values)
{
unsigned char val;
if (ioctl (fd, PPIGCTRL, &val) == -1)
{
fprintf(stderr, "Parallel port PPIGCTRL");
exit (1);
}
val &= ~controlbits;
val |= values;
if (ioctl (fd, PPISCTRL, &val) == -1)
{
fprintf(stderr, "Parallel port PPISCTRL");
exit (1);
}
}
#endif
/* Initialise a parallel port, given open fd */
int
lp_init (int fd)
{
#ifdef HAVE_LINUX_PPDEV_H
int mode;
#endif
#ifdef HAVE_LINUX_PPDEV_H
mode = PARPORT_MODE_PCSPP;
if (ioctl (fd, PPSETMODE, &mode) == -1)
{
fprintf(stderr, "Setting parallel port mode");
close (fd);
return(-1);
}
if (ioctl (fd, PPEXCL, NULL) == -1)
{
fprintf(stderr, "Parallel port is already in use.\n");
close (fd);
return(-1);
}
if (ioctl (fd, PPCLAIM, NULL) == -1)
{
fprintf(stderr, "Claiming parallel port.\n");
fprintf(stderr, "HINT: did you unload the lp kernel module?");
close (fd);
return(-1);
}
/* Enable CW & PTT - /STROBE bit (pin 1) */
parport_control (fd, PARPORT_CONTROL_STROBE, PARPORT_CONTROL_STROBE);
#endif
#ifdef HAVE_DEV_PPBUS_PPI_H
parport_control (fd, STROBE, STROBE);
#endif
lp_reset (fd);
return(0);
}
/* release ppdev and close port */
int
lp_free (int fd)
{
#ifdef HAVE_LINUX_PPDEV_H
lp_reset (fd);
/* Disable CW & PTT - /STROBE bit (pin 1) */
parport_control (fd, PARPORT_CONTROL_STROBE, 0);
ioctl (fd, PPRELEASE);
#endif
#ifdef HAVE_DEV_PPBUS_PPI_H
/* Disable CW & PTT - /STROBE bit (pin 1) */
parport_control (fd, STROBE, 0);
#endif
close (fd);
return(0);
}
/* set to a known state */
int
lp_reset (int fd)
{
#if defined (HAVE_LINUX_PPDEV_H) || defined (HAVE_DEV_PPBUS_PPI_H)
lp_ptt (fd, 0);
#endif
return(0);
}
/* SSB PTT keying - /INIT bit (pin 16) (inverted) */
int
lp_ptt (int fd, int onoff)
{
#ifdef HAVE_LINUX_PPDEV_H
if (onoff == 1)
parport_control (fd, PARPORT_CONTROL_INIT,
PARPORT_CONTROL_INIT);
else
parport_control (fd, PARPORT_CONTROL_INIT, 0);
#endif
#ifdef HAVE_DEV_PPBUS_PPI_H
if (onoff == 1)
parport_control (fd, nINIT,
nINIT);
else
parport_control (fd, nINIT, 0);
#endif
return(0);
}
/*
* ptt_parallel
*
* generic parallel unix PTT routine called indirectly from Fortran
*
* fd - already opened file descriptor
* ntx - pointer to fortran command on or off
* iptt - pointer to fortran command status on or off
*/
int
ptt_parallel(int fd, int *ntx, int *iptt)
{
if(*ntx) {
lp_ptt(fd, 1);
*iptt=1;
} else {
lp_ptt(fd, 0);
*iptt=0;
}
return(0);
}

35
q65w/libm65/rs.h
View File

@ -1,35 +0,0 @@
/* User include file for the Reed-Solomon codec
* Copyright 2002, Phil Karn KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
/* General purpose RS codec, 8-bit symbols */
void encode_rs_char(void *rs,unsigned char *data,unsigned char *parity);
int decode_rs_char(void *rs,unsigned char *data,int *eras_pos,
int no_eras);
void *init_rs_char(int symsize,int gfpoly,
int fcr,int prim,int nroots,
int pad);
void free_rs_char(void *rs);
/* General purpose RS codec, integer symbols */
void encode_rs_int(void *rs,int *data,int *parity);
int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras);
void *init_rs_int(int symsize,int gfpoly,int fcr,
int prim,int nroots,int pad);
void free_rs_int(void *rs);
/* CCSDS standard (255,223) RS codec with conventional (*not* dual-basis)
* symbol representation
*/
void encode_rs_8(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_8(unsigned char *data,int *eras_pos,int no_eras,int pad);
/* CCSDS standard (255,223) RS codec with dual-basis symbol representation */
void encode_rs_ccsds(unsigned char *data,unsigned char *parity,int pad);
int decode_rs_ccsds(unsigned char *data,int *eras_pos,int no_eras,int pad);
/* Tables to map from conventional->dual (Taltab) and
* dual->conventional (Tal1tab) bases
*/
extern unsigned char Taltab[],Tal1tab[];

16
q65w/libm65/rs2.h
View File

@ -1,16 +0,0 @@
/* User include file for the Reed-Solomon codec
* Copyright 2002, Phil Karn KA9Q
* May be used under the terms of the GNU General Public License (GPL)
*/
/* General purpose RS codec, integer symbols */
void encode_rs_int(void *rs,int *data,int *parity);
int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras, int calc_syn);
void *init_rs_int(int symsize,int gfpoly,int fcr,
int prim,int nroots,int pad);
void free_rs_int(void *rs);
/* Tables to map from conventional->dual (Taltab) and
* dual->conventional (Tal1tab) bases
*/
extern unsigned char Taltab[],Tal1tab[];

11
q65w/libm65/sec_midn.f90
View File

@ -1,11 +0,0 @@
real function sec_midn()
sec_midn=secnds(0.0)
return
end function sec_midn
subroutine sleep_msec(n)
call usleep(1000*n)
return
end subroutine sleep_msec

32
q65w/libm65/sleep.h
View File

@ -1,32 +0,0 @@
/*
* sleep.h 1.0 02/03/10
*
* Defines cross-platform sleep, usleep, etc.
*
* By Wu Yongwei
*
*/
#ifndef _SLEEP_H
#define _SLEEP_H
#ifdef _WIN32
# if defined(_NEED_SLEEP_ONLY) && (defined(_MSC_VER) || defined(__MINGW32__))
# include <stdlib.h>
# define sleep(t) _sleep((t) * 1000)
# else
# include <windows.h>
# define sleep(t) Sleep((t) * 1000)
# endif
# ifndef _NEED_SLEEP_ONLY
# define msleep(t) Sleep(t)
# define usleep(t) Sleep((t) / 1000)
# endif
#else
# include <unistd.h>
# ifndef _NEED_SLEEP_ONLY
# define msleep(t) usleep((t) * 1000)
# endif
#endif
#endif /* _SLEEP_H */

76
q65w/libm65/timeval.h
View File

@ -1,76 +0,0 @@
/*
* timeval.h 1.0 01/12/19
*
* Defines gettimeofday, timeval, etc. for Win32
*
* By Wu Yongwei
*
*/
#ifndef _TIMEVAL_H
#define _TIMEVAL_H
#ifdef _WIN32
#define WIN32_LEAN_AND_MEAN
#include <windows.h>
#include <time.h>
#ifndef __GNUC__
#define EPOCHFILETIME (116444736000000000i64)
#else
#define EPOCHFILETIME (116444736000000000LL)
#endif
//struct timeval {
// long tv_sec; /* seconds */
// long tv_usec; /* microseconds */
//};
/*
struct timezone {
int tz_minuteswest; // minutes W of Greenwich
int tz_dsttime; // type of dst correction
};
*/
__inline int gettimeofday(struct timeval *tv, struct timezone *tz)
{
FILETIME ft;
LARGE_INTEGER li;
__int64 t;
static int tzflag;
if (tv)
{
GetSystemTimeAsFileTime(&ft);
li.LowPart = ft.dwLowDateTime;
li.HighPart = ft.dwHighDateTime;
t = li.QuadPart; /* In 100-nanosecond intervals */
t -= EPOCHFILETIME; /* Offset to the Epoch time */
t /= 10; /* In microseconds */
tv->tv_sec = (long)(t / 1000000);
tv->tv_usec = (long)(t % 1000000);
}
if (tz)
{
if (!tzflag)
{
_tzset();
tzflag++;
}
tz->tz_minuteswest = _timezone / 60;
tz->tz_dsttime = _daylight;
}
return 0;
}
#else /* _WIN32 */
#include <sys/time.h>
#endif /* _WIN32 */
#endif /* _TIMEVAL_H */

514
q65w/libm65/tmoonsub.c
View File

@ -1,514 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define RADS 0.0174532925199433
#define DEGS 57.2957795130823
#define TPI 6.28318530717959
#define PI 3.1415927
/* ratio of earth radius to astronomical unit */
#define ER_OVER_AU 0.0000426352325194252
/* all prototypes here */
double getcoord(int coord);
void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
double range(double y);
double rangerad(double y);
double days(int y, int m, int dn, double hour);
double days_(int *y, int *m, int *dn, double *hour);
void moonpos(double, double *, double *, double *);
void sunpos(double , double *, double *, double *);
double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
double atan22(double y, double x);
double epsilon(double d);
void equatorial(double d, double *lon, double *lat, double *r);
void ecliptic(double d, double *lon, double *lat, double *r);
double gst(double d);
void topo(double lst, double glat, double *alp, double *dec, double *r);
double alt(double glat, double ha, double dec);
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
int daysinmonth(int y, int m);
int isleap(int y);
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
double *mrv, double *l, double *b, double *paxis);
static const char
*usage = " Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
"example: tmoon 200009 0 -00155 5230\n";
/*
getargs() gets the arguments from the command line, does some basic error
checking, and converts arguments into numerical form. Arguments are passed
back in pointers. Error messages print to stderr so re-direction of output
to file won't leave users blind. Error checking prints list of all errors
in a command line before quitting.
*/
void getargs(int argc, char *argv[], int *y, int *m, double *tz,
double *glong, double *glat) {
int date, latitude, longitude;
int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
int longminflag = 0, latminflag = 0, dflag = 0;
/* if not right number of arguments, then print example command line */
if (argc !=5) {
fprintf(stderr, usage);
exit(EXIT_FAILURE);
}
date = atoi(argv[1]);
*y = date / 100;
*m = date - *y * 100;
*tz = (double) atof(argv[2]);
longitude = atoi(argv[3]);
latitude = atoi(argv[4]);
*glong = RADS * getcoord(longitude);
*glat = RADS * getcoord(latitude);
/* set a flag for each error found */
if (*m > 12 || *m < 1) mflag = 1;
if (*y > 2500) yflag = 1;
if (date < 150001) dflag = 1;
if (fabs((float) *glong) > 180 * RADS) longflag = 1;
if (abs(longitude) % 100 > 59) longminflag = 1;
if (fabs((float) *glat) > 90 * RADS) latflag = 1;
if (abs(latitude) % 100 > 59) latminflag = 1;
if (fabs((float) *tz) > 12) tzflag = 1;
/* print all the errors found */
if (dflag == 1) {
fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
}
if (yflag == 1) {
fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
}
if (mflag == 1) {
fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
}
if (tzflag == 1) {
fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
}
if (longflag == 1) {
fprintf(stderr, "long: must be in range +/- 180 degrees\n");
}
if (longminflag == 1) {
fprintf(stderr, "long: last two digits are arcmin - max 59\n");
}
if (latflag == 1) {
fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
}
if (latminflag == 1) {
fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
}
/* quits if one or more flags set */
if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
exit(EXIT_FAILURE);
}
}
/*
returns coordinates in decimal degrees given the
coord as a ddmm value stored in an integer.
*/
double getcoord(int coord) {
int west = 1;
double glg, deg;
if (coord < 0) west = -1;
glg = fabs((double) coord/100);
deg = floor(glg);
glg = west* (deg + (glg - deg)*100 / 60);
return(glg);
}
/*
days() takes the year, month, day in the month and decimal hours
in the day and returns the number of days since J2000.0.
Assumes Gregorian calendar.
*/
double days(int y, int m, int d, double h) {
int a, b;
double day;
/*
The lines below work from 1900 march to feb 2100
a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
day = (double)a - 730531.5 + hour / 24;
*/
/* These lines work for any Gregorian date since 0 AD */
if (m ==1 || m==2) {
m +=12;
y -= 1;
}
a = y / 100;
b = 2 - a + a/4;
day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
+ d + b - 1524.5 - 2451545 + h/24;
return(day);
}
double days_(int *y0, int *m0, int *d0, double *h0)
{
return days(*y0,*m0,*d0,*h0);
}
/*
Returns 1 if y a leap year, and 0 otherwise, according
to the Gregorian calendar
*/
int isleap(int y) {
int a = 0;
if(y % 4 == 0) a = 1;
if(y % 100 == 0) a = 0;
if(y % 400 == 0) a = 1;
return(a);
}
/*
Given the year and the month, function returns the
number of days in the month. Valid for Gregorian
calendar.
*/
int daysinmonth(int y, int m) {
int b = 31;
if(m == 2) {
if(isleap(y) == 1) b= 29;
else b = 28;
}
if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
return(b);
}
/*
moonpos() takes days from J2000.0 and returns ecliptic coordinates
of moon in the pointers. Note call by reference.
This function is within a couple of arcminutes most of the time,
and is truncated from the Meeus Ch45 series, themselves truncations of
ELP-2000. Returns moon distance in earth radii.
Terms have been written out explicitly rather than using the
table based method as only a small number of terms is
retained.
*/
void moonpos(double d, double *lambda, double *beta, double *rvec) {
double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;
t = d / 36525;
L = range(218.3164591 + 481267.88134236 * t) * RADS;
D = range(297.8502042 + 445267.1115168 * t) * RADS;
M = range(357.5291092 + 35999.0502909 * t) * RADS;
M1 = range(134.9634114 + 477198.8676313 * t - .008997 * t * t) * RADS;
F = range(93.27209929999999 + 483202.0175273 * t - .0034029*t*t)*RADS;
e = 1 - .002516 * t;
dl = 6288774 * sin(M1);
dl += 1274027 * sin(2 * D - M1);
dl += 658314 * sin(2 * D);
dl += 213618 * sin(2 * M1);
dl -= e * 185116 * sin(M);
dl -= 114332 * sin(2 * F) ;
dl += 58793 * sin(2 * D - 2 * M1);
dl += e * 57066 * sin(2 * D - M - M1) ;
dl += 53322 * sin(2 * D + M1);
dl += e * 45758 * sin(2 * D - M);
dl -= e * 40923 * sin(M - M1);
dl -= 34720 * sin(D) ;
dl -= e * 30383 * sin(M + M1) ;
dl += 15327 * sin(2 * D - 2 * F) ;
dl -= 12528 * sin(M1 + 2 * F);
dl += 10980 * sin(M1 - 2 * F);
lm = rangerad(L + dl / 1000000 * RADS);
dB = 5128122 * sin(F);
dB += 280602 * sin(M1 + F);
dB += 277693 * sin(M1 - F);
dB += 173237 * sin(2 * D - F);
dB += 55413 * sin(2 * D - M1 + F);
dB += 46271 * sin(2 * D - M1 - F);
dB += 32573 * sin(2 * D + F);
dB += 17198 * sin(2 * M1 + F);
dB += 9266 * sin(2 * D + M1 - F);
dB += 8822 * sin(2 * M1 - F);
dB += e * 8216 * sin(2 * D - M - F);
dB += 4324 * sin(2 * D - 2 * M1 - F);
bm = dB / 1000000 * RADS;
dR = -20905355 * cos(M1);
dR -= 3699111 * cos(2 * D - M1);
dR -= 2955968 * cos(2 * D);
dR -= 569925 * cos(2 * M1);
dR += e * 48888 * cos(M);
dR -= 3149 * cos(2 * F);
dR += 246158 * cos(2 * D - 2 * M1);
dR -= e * 152138 * cos(2 * D - M - M1) ;
dR -= 170733 * cos(2 * D + M1);
dR -= e * 204586 * cos(2 * D - M);
dR -= e * 129620 * cos(M - M1);
dR += 108743 * cos(D);
dR += e * 104755 * cos(M + M1);
dR += 79661 * cos(M1 - 2 * F);
rm = 385000.56 + dR / 1000;
*lambda = lm;
*beta = bm;
/* distance to Moon must be in Earth radii */
*rvec = rm / 6378.14;
}
/*
topomoon() takes the local siderial time, the geographical
latitude of the observer, and pointers to the geocentric
equatorial coordinates. The function overwrites the geocentric
coordinates with topocentric coordinates on a simple spherical
earth model (no polar flattening). Expects Moon-Earth distance in
Earth radii. Formulas scavenged from Astronomical Almanac 'low
precision formulae for Moon position' page D46.
*/
void topo(double lst, double glat, double *alp, double *dec, double *r) {
double x, y, z, r1;
x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
z = *r * sin(*dec) - sin(glat);
r1 = sqrt(x*x + y*y + z*z);
*alp = atan22(y, x);
*dec = asin(z / r1);
*r = r1;
}
/*
moontransit() takes date, the time zone and geographic longitude
of observer and returns the time (decimal hours) of lunar transit
on that day if there is one, and sets the notransit flag if there
isn't. See Explanatory Supplement to Astronomical Almanac
section 9.32 and 9.31 for the method.
*/
double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
double hm, ht, ht1, lon, lat, rv, dnew, lst;
int itcount;
ht1 = 180 * RADS;
ht = 0;
itcount = 0;
*notransit = 0;
do {
ht = ht1;
itcount++;
dnew = days(y, m, d, ht * DEGS/15) - tz/24;
lst = gst(dnew) + glong;
/* find the topocentric Moon ra (hence hour angle) and dec */
moonpos(dnew, &lon, &lat, &rv);
equatorial(dnew, &lon, &lat, &rv);
topo(lst, glat, &lon, &lat, &rv);
hm = rangerad(lst - lon);
ht1 = rangerad(ht - hm);
/* if no convergence, then no transit on that day */
if (itcount > 30) {
*notransit = 1;
break;
}
}
while (fabs(ht - ht1) > 0.04 * RADS);
return(ht1);
}
/*
Calculates the selenographic coordinates of either the sub Earth point
(optical libration) or the sub-solar point (selen. coords of centre of
bright hemisphere). Based on Meeus chapter 51 but neglects physical
libration and nutation, with some simplification of the formulas.
*/
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
double i, f, omega, w, y, x, a, t, eps;
t = day / 36525;
i = 1.54242 * RADS;
eps = epsilon(day);
f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
w = lambda - omega;
y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
x = cos(w) * cos(beta);
a = atan22(y, x);
*l = a - f;
/* kludge to catch cases of 'round the back' angles */
if (*l < -90 * RADS) *l += TPI;
if (*l > 90 * RADS) *l -= TPI;
*b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));
/* pa pole axis - not used for Sun stuff */
x = sin(i) * sin(omega);
y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
w = atan22(x, y);
*p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
}
/*
Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
eq coords Sun
Returns: position angle of bright limb wrt NCP, percentage illumination
of Sun
*/
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
double x, y, phi, i;
y = cos(sdec) * sin(sra - lra);
x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
*pabl = atan22(y, x);
phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
i = atan22(sin(phi) , (dr - cos(phi)));
*ill = 0.5*(1 + cos(i));
}
/*
sunpos() takes days from J2000.0 and returns ecliptic longitude
of Sun in the pointers. Latitude is zero at this level of precision,
but pointer left in for consistency in number of arguments.
This function is within 0.01 degree (1 arcmin) almost all the time
for a century either side of J2000.0. This is from the 'low precision
fomulas for the Sun' from C24 of Astronomical Alamanac
*/
void sunpos(double d, double *lambda, double *beta, double *rvec) {
double L, g, ls, bs, rs;
L = range(280.461 + .9856474 * d) * RADS;
g = range(357.528 + .9856003 * d) * RADS;
ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
bs = 0;
rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
*lambda = ls;
*beta = bs;
*rvec = rs;
}
/*
this routine returns the altitude given the days since J2000.0
the hour angle and declination of the object and the latitude
of the observer. Used to find the Sun's altitude to put a letter
code on the transit time, and to find the Moon's altitude at
transit just to make sure that the Moon is visible.
*/
double alt(double glat, double ha, double dec) {
return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
}
/* returns an angle in degrees in the range 0 to 360 */
double range(double x) {
double a, b;
b = x / 360;
a = 360 * (b - floor(b));
if (a < 0)
a = 360 + a;
return(a);
}
/* returns an angle in rads in the range 0 to two pi */
double rangerad(double x) {
double a, b;
b = x / TPI;
a = TPI * (b - floor(b));
if (a < 0)
a = TPI + a;
return(a);
}
/*
gets the atan2 function returning angles in the right
order and range
*/
double atan22(double y, double x) {
double a;
a = atan2(y, x);
if (a < 0) a += TPI;
return(a);
}
/*
returns mean obliquity of ecliptic in radians given days since
J2000.0.
*/
double epsilon(double d) {
double t = d/ 36525;
return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
}
/*
replaces ecliptic coordinates with equatorial coordinates
note: call by reference destroys original values
R is unchanged.
*/
void equatorial(double d, double *lon, double *lat, double *r) {
double eps, ceps, seps, l, b;
l = *lon;
b = * lat;
eps = epsilon(d);
ceps = cos(eps);
seps = sin(eps);
*lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
*lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
}
/*
replaces equatorial coordinates with ecliptic ones. Inverse
of above, but used to find topocentric ecliptic coords.
*/
void ecliptic(double d, double *lon, double *lat, double *r) {
double eps, ceps, seps, alp, dec;
alp = *lon;
dec = *lat;
eps = epsilon(d);
ceps = cos(eps);
seps = sin(eps);
*lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
*lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
}
/*
returns the siderial time at greenwich meridian as
an angle in radians given the days since J2000.0
*/
double gst( double d) {
double t = d / 36525;
double theta;
theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
return(theta * RADS);
}
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
double *mrv, double *l, double *b, double *paxis)
{
double mlambda, mbeta;
double malpha, mdelta;
double lst, mhr;
double tlambda, tbeta, trv;
lst = gst(*day) + *glong;
/* find Moon topocentric coordinates for libration calculations */
moonpos(*day, &mlambda, &mbeta, mrv);
malpha = mlambda;
mdelta = mbeta;
equatorial(*day, &malpha, &mdelta, mrv);
topo(lst, *glat, &malpha, &mdelta, mrv);
mhr = rangerad(lst - malpha);
*moonalt = alt(*glat, mhr, mdelta);
/* Optical libration and Position angle of the Pole */
tlambda = malpha;
tbeta = mdelta;
trv = *mrv;
ecliptic(*day, &tlambda, &tbeta, &trv);
libration(*day, tlambda, tbeta, malpha, l, b, paxis);
}

7
q65w/libm65/usleep.c
View File

@ -1,7 +0,0 @@
#include <unistd.h>
/* usleep(3) */
void usleep_(unsigned long *microsec)
{
usleep(*microsec);
}

70
q65w/libm65/wrapkarn.c
View File

@ -1,70 +0,0 @@
#include <math.h>
#include <stdio.h>
#include <float.h>
#include <limits.h>
#include <stdlib.h>
#include "rs.h"
static void *rs;
static int first=1;
void rs_encode_(int *dgen, int *sent)
// Encode JT65 data dgen[12], producing sent[63].
{
int dat1[12];
int b[51];
int i;
if(first) {
// Initialize the JT65 codec
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
// Reverse data order for the Karn codec.
for(i=0; i<12; i++) {
dat1[i]=dgen[11-i];
}
// Compute the parity symbols
encode_rs_int(rs,dat1,b);
// Move parity symbols and data into sent[] array, in reverse order.
for (i = 0; i < 51; i++) sent[50-i] = b[i];
for (i = 0; i < 12; i++) sent[i+51] = dat1[11-i];
}
void rs_decode_(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
// Decode JT65 received data recd0[63], producing decoded[12].
// Erasures are indicated in era0[numera]. The number of corrected
// errors is *nerr. If the data are uncorrectable, *nerr=-1 is returned.
{
int numera;
int i;
int era_pos[50];
int recd[63];
if(first) {
rs=init_rs_int(6,0x43,3,1,51,0);
first=0;
}
numera=*numera0;
for(i=0; i<12; i++) recd[i]=recd0[62-i];
for(i=0; i<51; i++) recd[12+i]=recd0[50-i];
if(numera)
for(i=0; i<numera; i++) era_pos[i]=era0[i];
*nerr=decode_rs_int(rs,recd,era_pos,numera);
for(i=0; i<12; i++) decoded[i]=recd[11-i];
}
void rs_encode__(int *dgen, int *sent)
{
rs_encode_(dgen, sent);
}
void rs_decode__(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
{
rs_decode_(recd0, era0, numera0, decoded, nerr);
}