diff --git a/q65w/libm65/CMakeLists.txt b/q65w/libm65/CMakeLists.txt index cda8cb193..bc47030fb 100644 --- a/q65w/libm65/CMakeLists.txt +++ b/q65w/libm65/CMakeLists.txt @@ -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} diff --git a/q65w/libm65/cutil.c b/q65w/libm65/cutil.c deleted file mode 100644 index 69a526b07..000000000 --- a/q65w/libm65/cutil.c +++ /dev/null @@ -1,93 +0,0 @@ -#include -#include -#include -#include -#include -#include -// #include -// #include -// #include -#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); -} diff --git a/q65w/libm65/decode_rs.c b/q65w/libm65/decode_rs.c deleted file mode 100644 index 91f582ac1..000000000 --- a/q65w/libm65/decode_rs.c +++ /dev/null @@ -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 -#endif - -#include - -#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 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 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 - -#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 -#include -#include -#include -#include -#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= 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; -} diff --git a/q65w/libm65/gran.c b/q65w/libm65/gran.c deleted file mode 100644 index 24b986503..000000000 --- a/q65w/libm65/gran.c +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include - -/* 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; -} diff --git a/q65w/libm65/igray.c b/q65w/libm65/igray.c deleted file mode 100644 index 395f79712..000000000 --- a/q65w/libm65/igray.c +++ /dev/null @@ -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); -} diff --git a/q65w/libm65/init_rs.c b/q65w/libm65/init_rs.c deleted file mode 100644 index 876819f8c..000000000 --- a/q65w/libm65/init_rs.c +++ /dev/null @@ -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 - -#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<= (1<= (1<mm = symsize; - rs->nn = (1<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;inn;i++){ - rs->index_of[sr] = i; - rs->alpha_to[i] = sr; - sr <<= 1; - if(sr & (1<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; -} diff --git a/q65w/libm65/int.h b/q65w/libm65/int.h deleted file mode 100644 index ada5bfd4c..000000000 --- a/q65w/libm65/int.h +++ /dev/null @@ -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<= 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); - - - - diff --git a/q65w/libm65/map65a.f90 b/q65w/libm65/map65a.f90 index c1b86bd16..f3bd5aa5b 100644 --- a/q65w/libm65/map65a.f90 +++ b/q65w/libm65/map65a.f90 @@ -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) diff --git a/q65w/libm65/ptt.c b/q65w/libm65/ptt.c deleted file mode 100644 index 0f99d1be8..000000000 --- a/q65w/libm65/ptt.c +++ /dev/null @@ -1,43 +0,0 @@ -#include -#include - -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; -} diff --git a/q65w/libm65/ptt_unix.c b/q65w/libm65/ptt_unix.c deleted file mode 100644 index 1c583f641..000000000 --- a/q65w/libm65/ptt_unix.c +++ /dev/null @@ -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 - * 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 -#endif -#if STDC_HEADERS -# include -# include -#else -# if HAVE_STDLIB_H -# include -# endif -#endif -#if HAVE_UNISTD_H -# include -#endif -#if HAVE_SYS_STAT_H -# include -#endif -#if HAVE_SYS_IOCTL_H -# include -#endif -#if HAVE_FCNTL_H -# include -#endif - -#ifdef HAVE_LINUX_PPDEV_H -# include -# include -#endif -#ifdef HAVE_DEV_PPBUS_PPI_H -# include -# include -#endif - -int lp_reset (int fd); -int lp_ptt (int fd, int onoff); - -#ifdef HAVE_SYS_STAT_H -# include -#endif -#if (defined(__unix__) || defined(unix)) && !defined(USG) -# include -#endif - -#include -/* 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); -} diff --git a/q65w/libm65/rs.h b/q65w/libm65/rs.h deleted file mode 100644 index 06cbe344f..000000000 --- a/q65w/libm65/rs.h +++ /dev/null @@ -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[]; diff --git a/q65w/libm65/rs2.h b/q65w/libm65/rs2.h deleted file mode 100644 index c2b807d15..000000000 --- a/q65w/libm65/rs2.h +++ /dev/null @@ -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[]; diff --git a/q65w/libm65/sec_midn.f90 b/q65w/libm65/sec_midn.f90 deleted file mode 100644 index 0bbe62c2c..000000000 --- a/q65w/libm65/sec_midn.f90 +++ /dev/null @@ -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 diff --git a/q65w/libm65/sleep.h b/q65w/libm65/sleep.h deleted file mode 100644 index df60bc92a..000000000 --- a/q65w/libm65/sleep.h +++ /dev/null @@ -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 -# define sleep(t) _sleep((t) * 1000) -# else -# include -# 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 -# ifndef _NEED_SLEEP_ONLY -# define msleep(t) usleep((t) * 1000) -# endif -#endif - -#endif /* _SLEEP_H */ diff --git a/q65w/libm65/timeval.h b/q65w/libm65/timeval.h deleted file mode 100644 index 83c77d5a1..000000000 --- a/q65w/libm65/timeval.h +++ /dev/null @@ -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 -#include - -#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 - -#endif /* _WIN32 */ - -#endif /* _TIMEVAL_H */ diff --git a/q65w/libm65/tmoonsub.c b/q65w/libm65/tmoonsub.c deleted file mode 100644 index 29ac28b49..000000000 --- a/q65w/libm65/tmoonsub.c +++ /dev/null @@ -1,514 +0,0 @@ -#include -#include -#include - -#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); -} diff --git a/q65w/libm65/usleep.c b/q65w/libm65/usleep.c deleted file mode 100644 index 21d242a68..000000000 --- a/q65w/libm65/usleep.c +++ /dev/null @@ -1,7 +0,0 @@ -#include - -/* usleep(3) */ -void usleep_(unsigned long *microsec) -{ - usleep(*microsec); -} diff --git a/q65w/libm65/wrapkarn.c b/q65w/libm65/wrapkarn.c deleted file mode 100644 index 9e0a51caf..000000000 --- a/q65w/libm65/wrapkarn.c +++ /dev/null @@ -1,70 +0,0 @@ -#include -#include -#include -#include -#include -#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