/***********************************************************************
 *                                                                     *
 * Program to control LORAN-C radio                                    *
 *                                                                     *
 * Utility subroutines and input/output routines                       *
 *                                                                     *
 ***********************************************************************
 */

#include "loran.h"
#include <stdlib.h>

/*
 * Local function declarations
 */
double gauss(double);
double envelope(double);
struct station *receive(void);
void loadreg(struct station *);
void loadgri(struct station *);
void timerq(struct station *, long);
void adj_freq(struct station *, double);
char station_type(struct station *);
double local_time(struct station *);

/*
 * Imported from loran.c
 */
extern int vco;					/* vco dac signal */
extern int vcodac;				/* vco dac bias (dac a) */
extern int agc;					/* agc dac signal */
extern int agcdac;				/* agc dac bias (dac b) */
extern int peak_detect;			/* s-signal (adc chan 2) */
extern double rcvr_delay;		/* antenna-receiver cable delay (us) */

/*
 * Imported from gri.c
 */
extern int nstation;			/* number of stations */
extern struct station *chain[];	/* vector of station pointers */
extern double delay[NSTA];		/* vector of station offsets */
extern int par;					/* mode register */
extern long gri;				/* group repetition interval (cycle) */
extern long fri;				/* frame interval (2 * gri) (cycle) */
extern long offset;				/* current frame offset (cycle) */
extern double master_delay;		/* master pulse-group delay */

/*
 * Local data declarations.
 */
double epoch = 0;				/* elapsed time since startup (us) */
struct station *lorsta[NSTA];	/* vector of station pointers */
long deltat[NSTA];				/* vector of increments to next alarm */
int wstation = 0;				/* number of alarms queued */

/*
 * The folowing stuff is used only for simulation. all signal levels are
 * defined relative to the peak of the master pulse-group envelope. The
 * present values are chosen at the lower limit of the operating
 * envelope. For simulating the US North East LORAN-C chain (gri 9960),
 * both the slave pulse-group peak and noise level are set -6 dB (0.5)
 * below the master pulse-group peak. In addition, the initial master-
 * clock frequency offset is set at 100 ppm. The master offset is set
 * arbitrarily within the frame interval. These choices result in a
 * marvellous bit of clanking and banging upon initial synchronization,
 * but all five stations eventually synchronize and settle down to less
 * than one master-clock period (200 ns) from the predicted delays. With
 * noise levels higher than this, the master fails to reliably
 * synchronize. Once synchronized, the noise level can be increased to
 * 6 dB (2) with only occasional slave-clock slips of 200 ns. In this
 * regime the snr of the reference cycle, about 6 dB below the pulse
 * peak, for the slaves is -18 dB. Allowing 12 dB for synchronous
 * integration of the 16-pulse LORAN-C signal over the frame interval,
 * the slaves are operating some -30 dB down in the muck, which is most
 * respectable. When the noise is increased to 12 dB (4) some strobe
 * slips of 10 us occur with frequent slave-clock slips up to 5 us. When
 * the noise is increased to 18 dB (8), severe strobe slips occur up to
 * +-50 us. The actual receiver in Newark, DE, operates with the 9960
 * chain at an equivalent noise level of about .25 (-12 dB) relative to
 * the master at Seneca, NY.
 */
#if defined(DEBUG)
double signal_master = 1;		/* master pulse-group peak */
double signal_slave = .5;		/* slave pulse-group peak (-6 dB) */
double noise = .01;				/* noise level (-40 dB) */
double test_slew = -.1;			/* frequency offset (0.1ppm) */
double master_offset = 12000;	/* master (m) offset (us) */

double freq = 0;				/* vco frequency offset (ppm) */
double utc = 0;					/* universal time (us) */
double gate_offset;				/* gate offset in frame (us) */
double strobe_offset;			/* strobe offset in pulse gate (us) */

/*
 * Subroutine gauss()
 *
 * This subroutine generates a random number uniformly distributed
 * over the interval [-1, 1], then transforms the distribution to
 * a Gaussian distribution with zero mean and specified standard
 * deviation. See the header file for additional data.
 *
 * Calling sequence: x = gauss(sigma)
 * 	sigma		standard deviation of noise
 *	x			noise sample
 *
 * Variables and functions used (math library)
 *	rand()		generate uniform random sample over [0, 32767]
 *	sqrt()		square rootsie
 *	log()		log base e
 */
double gauss(sigma)
	double sigma;				/* standard deviation of noise */
	{
 	double x, y;				/* double temps */

	x = ((double)rand() / 16384 - 1);
	if (x > 0)
		y = sigma / sqrt(2) * log(1 / x);
	else if (x < 0)
		y = -sigma / sqrt(2) * log(1 / -x);
	else
		y = 0;
	return (y);
	}

/*
 * Subroutine envelope()
 *
 * This subroutine generates the standard LORAN-C envelope pulse
 * normalized for unit peak amplitude. By numeric methods the integral
 * of the normalized envelope is 120.072 and the integral of the square
 * of the normalized envelope is 83.1769. See the header file for
 * additional details. Note that 65 * 65 * exp(-2) = 571.792.
 *
 * Calling sequence: x = envelope(t)
 * 	t			simulation time (us)
 *	x			pulse amplitude
 *
 * Variables and functions used (math library)
 *	sin()		sine function
 *	cos()		cosine function
 *	exp()		exponent base e
 */
double envelope(t)
	double t;					/* simulation time */
	{

	if (t > 0 && t < 300)
		return (t * t * exp(-2 * t / 65) / 571.792);
	else
		return (0);
	}

#endif /* DEBUG */

/*
 * Subroutine receive()
 *
 * This subroutine assumes the adc has been enabled at the end of the
 * next gri. It sets up for the next gri, waits for completion and reads
 * the i integrator (adc channel 0), q integrator (adc channel 1) and
 * s detector (adc channel 2).
 *
 * For the purposes of simulation this subroutine can generate a
 * faithful representation of the LORAN-C envelope, carrier and noise
 * characteristics, depending on receiver mode. It calculates the
 * receiver gain, computes the time position of the pulse to generate
 * the i and q signals and draws a sample from a Gaussian noise process.
 * The resulting signals are treated as increments in the integration
 * process and are normally accumulated over the a and b intervals of
 * the frame.
 *
 * Calling sequence: sptr = receive()
 * 	sptr		station structure pointer
 *
 * Variables and functions used
 *	gri			group repetition interval (cycle)
 *	fri			frame interval (cycle)
 *	wstation	count of alarms queued
 *	lorsta[]	vector of station pointers
 *	deltat[]	vector of increments to next alarm
 * 	sptr->isig	i channel integrator accumulator
 * 	sptr->qsig	q channel integrator accumulator
 *	sptr->mode	program mode
 *	sptr->codesw a/b interval indicator (> 0: a, <= 0: b)
 *	offset		gri phase offset (cycke)
 *	peak_detect	peak-reading agc accumulator
 *	loadreg()	load registers for next gri interval
 *	loadgri()	load gri interval
 *	station_type() returns 'm' for master, 's' for slave
 */
struct station *receive() {		/* station structure pointer */

	struct station *sptr;		/* station pointer */
	int i, j;					/* int temps */
	long ltemp;					/* long temp */
	struct station *rptr;		/* temp station pointer */
#if defined( DEBUG)
	double t;					/* simulation time (us) */
	double gainfac;				/* receiver gain */
	double sigfac;				/* received signal level*/
	double u, dtemp;			/* double temporaries */
#endif

	/*
	 * Fetch the next station on the alarm queue and load registers
	 */
	sptr = lorsta[0];
	ltemp = deltat[0];
	for (i = 1; i < wstation; i++) {
		lorsta[i - 1] = lorsta[i];
		deltat[i - 1] = deltat[i];
		}
	wstation--;
	loadgri(sptr);

	/*
	 * Wait for the next gri, update the master clock and read the adc
	 * i, q and s channels
	 */
#if !defined(DEBUG)
	par = outp(PAR, par | ENG); outp(ADC, I); /* i */
	while ((inp(ADCGO) & DONE) == 0);
#endif /* DEBUG */
	epoch += ltemp * CYCLE;
	offset += ltemp;
	while (offset >= fri)
		offset -= fri;
	while (offset < 0)
		offset += fri;
	if (sptr->codesw <= 0) {
		sptr->isig = 0; sptr->qsig = 0; peak_detect = 0;
		}
#if defined( DEBUG)
	/*
	 * Calculate the signal level, gain factor and time offset of the
	 * pulse-group gate (in microseconds). Establish the adc offset for
	 * the i, q and s channels. Note that in all except the mode 5/6
	 * cases the receiver needs only the rms amplitude, so only the i
	 * channel is calculated. The s channel is presently not used, so
	 * its calculation is a fake.
	 */
	utc += ltemp * CYCLE * (1 + freq / 1e6);
	dtemp = CYCLE * fri;
	while (utc > dtemp)
		utc -= dtemp;
	while (utc < 0)
		utc += dtemp;
	t = utc + strobe_offset - rcvr_delay - master_offset - RCVDELAY;
	if (sptr->codesw <= 0)
		t -= CYCLE * gri;
	while (t > dtemp - NENV * CYCLE)
		t -= dtemp;
	while (t < -NENV * CYCLE)
		t += dtemp;
	for (i = 0; i < nstation; i++) {
		rptr = chain[i];
		if (t < rptr->ems_delay + rptr->path_delay + NENV * CYCLE)
			break;
		}
	t -= rptr->ems_delay + rptr->path_delay;
	gainfac = (1 - (agc - AGC) / 400) * 50;
	if (station_type(sptr) == 'm' && i == 0)
		sigfac = signal_master;
	else
		sigfac = signal_slave;
	sptr->isig += ADCOFS; sptr->qsig += ADCOFS; peak_detect += ADCOFS;
	switch(sptr->mode) {

		/*
		 * In mode 1 the LORAN-C signal is suppressed and the noise
		 * reduced by 40 dB to simulate the receiver at minimum gain.
		 */
		case MODE_MIN:
		sptr->isig += (int)(gainfac * gauss(noise / 100));
		break;

		/*
		 * In mode 2 the LORAN-C signal is suppressed and the noise
		 * increased by 20 dB to simulate the receiver under no-saignal
		 * conditions at maximum gain.
		 */
		case MODE_MAX:
		sptr->isig += (int)(gainfac * gauss(noise * 10));
		break;

		/*
		 * In mode 3 the LORAN-C signal is suppressed and the noise
		 * maintained at normal level to simulate the receiver under no-
		 * signal conditions at normal gain.
		 */
		case MODE_AGC:
		sptr->isig += (int)(gainfac * gauss(noise));
		break;

		/*
		 * In mode 4 the LORAN-C envelope signal plus noise is produced
		 * to simulate the receiver during the pulse-group search
		 * phase. The operation of the envelop integrators is
		 * simulated by numerical integration. The noise is of course
		 * statistically independent of the signal, so we cheat and
		 * calculate it only once and reduce it by the square root of
		 * the number of samples.
		 */
		case MODE_PSRCH:
		dtemp = gauss(noise) / sqrt(NENV);
		if (t > -300 && t < NENV * CYCLE) {
			for (u = 0; u < NENV * CYCLE; u += CYCLE)
				if (t + u > 0 && t + u < NENV * CYCLE)
					dtemp += sigfac * envelope(u + t) / NENV;
			}
		sptr->isig += (int)(gainfac * dtemp);
 		break;

		/*
		 * In mode 5 the LORAN-C envelope, carrier and noise are
		 * produced for both the i and q channels.
		 */
		default:
		sptr->isig += (int)(gainfac * (gauss(noise) + sigfac *
			envelope(t) * sin(2 * PI / CYCLE * t)));
		sptr->qsig += (int)(gainfac * (gauss(noise) + sigfac *
			envelope(t) * cos(2 * PI / CYCLE * t)));
		}
#else

	/*
	 * This is the input/output code sequence which retrieves the
	 * hardware i-, q- and s-channel values for the real radio. Not
	 * nearly as dramatic as the simulation code. This stuff is
	 * presently done by spinning on busy bits; it should be using
	 * interrupts. Later.
	 */
	sptr->isig += inp(ADC);
	par = outp(PAR, par & ~ENG);		/* q */
	outp(ADC, Q); outp(ADCGO, START);
	while ((inp(ADCGO) & DONE) == 0);
	sptr->qsig += inp(ADC);
	if (sptr->mode < MODE_PSRCH) {
		outp(ADC, S); outp(ADCGO, START); /* agc (calibration only) */
		while ((inp(ADCGO) & DONE) == 0);
		peak_detect += inp(ADC);
		}
#endif /* DEBUG */
	if (wstation > 0)
		rptr = lorsta[0];
	else
		rptr = sptr;
	loadreg(rptr);
	return(sptr);
	}

/*
 * Subroutine loadreg()
 *
 * This subroutine loads the registers and dacs for the next station
 * gri interval. This must be done before the beginning of the next
 * pulse-group window. Before exit the adc is set to sample the i
 * channel automatically at the end of the next gri interval.
 *
 * Calling sequence: loadreg(sptr)
 * 	sptr		station structure pointer
 *
 * Variables and functions used
 *
 *	sptr->level	agc level
 *	sptr->env	cycle strobe (cycle)
 *	sptr->slew	cycle phase (clock)
 *	sptr->mode	station mode
 *	sptr->codesw a/b interval indicator (> 0: a, <= 0: b)
 *	vco			vco dac
 *	agc			agc dac
 *	par			parameter register
 *	station_type() returns 'm' for master, 's' for slave
 */
void loadreg(sptr)
	struct station *sptr;		/* station structure pointer */
	{
	long ltemp;					/* long temp */
	double dtemp;				/* double temp */

	dtemp = sptr->level;		/* agc */
	if (dtemp > 4095)
		dtemp = 4095;
	if (dtemp < 0)
		dtemp = 0;
	agc = (int)dtemp;
	if (sptr->mode > MODE_PSRCH) { /* parameter register and strobe */
		par = par & 0xf0 | 0x0c;
		}
	else {
		par = par & 0xf0 | 0x0a;
		}
	ltemp = sptr->env * CLOCK + sptr->slew + RCVDELAY * (CLOCK / CYCLE);
#if defined(DEBUG)
	strobe_offset = ltemp * 1. / (CLOCK / CYCLE);
	freq = test_slew + Ko * CYCLE * (vco - vcodac);
#else
	outp(TGC, LOADDP + 0x0c);
	outp(TGD, (int)ltemp); outp(TGD, (int)ltemp >> 8);
	if (station_type(sptr) == 'm') /* code reg */
		if (sptr->codesw <= 0)
			outp(CODE, MPCB);
		else
			outp(CODE, MPCA);
	else
		if (sptr->codesw <= 0)
			outp(CODE, SPCB);
		else
			outp(CODE, SPCA);
	par = outp(PAR, par & ~MSB); outp(DACB, agc); /* agc dac */
	par = outp(PAR, par | MSB); outp(DACB, agc >> 8);
	par = outp(PAR, par & ~MSB); outp(DACA, vco); /* vco dac */
	par = outp(PAR, par | MSB); outp(DACA, vco >> 8);
/*	par = outp(PAR, par | ENG); outp(ADC, I); /* i */
#endif
	}

/*
 * Subroutine loadgri()
 *
 * This subroutine loads the gri counter buffer for the next station in
 * sequence. This must be done before the end of the next pulse-group
 * gate, at which time the buffer is loaded into the gri counter. At the
 * time of this call the registers for the first station on the alarm
 * queue have already been loaded by the software and the gri counter
 * loaded automatically from the gri counter buffer. Therefore, this
 * subroutine loads the gri counter buffer for the next (first) station
 * on the alarm queue. However, the station now to run has to be checked
 * for a pending display statement. If so, the gri counter buffer is
 * increased by one frame time (fri) to allow time for the display
 * operation to complete. The same trick could be used in case a disk
 * access could happen as the result of future monitoring code, for
 * example.
 *
 * Calling sequence: loadgri(sptr)
 * 	sptr		station structure pointer
 *
 * Variables and functions used
 *	gri			group repetition interval (cycle)
 *	fri			frame interval (cycle)
 *	wstation	count of alarms queued
 *	lorsta[]	vector of station pointers
 *	deltat[]	vector of increments to next alarm
 *	sptr->codesw a/b interval indicator (> 0: a, <= 0: b)
 *	sptr->report[] text string to be displayed (\0 if empty)
 */
void loadgri(sptr)
	struct station *sptr;		/* station structure pointer */
	{
	int i;						/* int temps */
	long ltemp, guard;			/* long temps */

	guard = PGUARD;
	if (wstation > 0) {
		if (sptr->codesw <= 0 && sptr->report[0] != '\0')
			guard = DGUARD;
		while (deltat[0] < guard)
			deltat[0] += fri;
		ltemp = deltat[0];
		}
	else
		ltemp = gri;
	ltemp -= 800;
#if !defined(DEBUG)
	outp(TGC, LOADDP + 0x0a);
	outp(TGD, (int)ltemp); outp(TGD, (int)ltemp >> 8);
#endif
	}

/*
 * Subroutine timerq()
 *
 * This subroutine is called at the end of gri processing. It sets an
 * alarm to wake up the next station following a computed delay. Note
 * that the stations are always served in round-robin fashion: first the
 * master (m), then the slaves (w, x, y, z) in order of emission delay.
 * In this way all stations can be served in the same gri.
 *
 * Calling sequence: timerq(sptr, delay)
 * 	sptr		station structure pointer
 * 	delay		offset for next alarm (cycle)
 *
 * Variables and functions used
 *	gri			group repetition interval (cycle)
 *	fri			frame interval (cycle)
 *	wstation	count of alarms queued
 *	lorsta[]	vector of station pointers
 *	deltat[]	vector of increments (cycle) to next alarm
 */
void timerq(sptr, delay)
	struct station *sptr;		/* station structure pointer */
	long delay;					/* phase offset (cycle) */
	{
	long ltemp;					/* long temps */
	int i;						/* int temps */

	ltemp = delay + gri;
	for (i = 0; i < wstation; i++)
		ltemp -= deltat[i];
	while (ltemp < 0)
		ltemp += fri;
	lorsta[wstation] = sptr;
	deltat[wstation] = ltemp;
	wstation++;
	}

/*
 * Subroutine adj_freq()
 *
 * This subroutine implements a type-II, adaptive-parameter, phase-lock
 * loop, which adjusts the voltage-controlled oscillator (vco) phase and
 * frequency to agree with the master station of the LORAN-C chain.
 * Since slave stations are phase locked to the master station, only the
 * phase has to be resolved using a simple type-I loop. In this case the
 * phase adjustment is done by slewing the strobe over the cycle in 200-
 * ns (CLOCK) steps.
 *
 * Calling sequence: adj_freq(sptr, nphase)
 *	sptr		station structure pointer
 *	phase		phase detector output
 *
 * Inputs
 *	epoch		elapsed time since startup (us)
 *	sptr->update time of last pll update (us)
 *
 * Outputs
 *	sptr->vphase pll phase adjustment
 *	sptr->vfreq	pll frequency adjustment
 *	sptr->slew	digital pll adjustment (200 ns)
 *	sptr->pfac	phase time constant
 *	sptr->psnr	mean phase error
 *	station_type() returns 'm' for master, 's' for slave
 */
void adj_freq(sptr, phase)
	struct station *sptr;		/* station pointer */
	double phase;				/* phase adjustment */
	{
	double time_constant, sample_interval;
	double dtemp;				/* double temp */

	sample_interval = epoch - sptr->update;
	sptr->update = epoch;
	sptr->vphase = phase / (VGAIN * sptr->pfac);
	sptr->vfreq += phase * sample_interval /
		(16e6 * VGAIN * sptr->pfac * sptr->pfac);
	if (station_type(sptr) == 'm') {
		dtemp = vcodac - (sptr->vphase + sptr->vfreq);
		if (dtemp > 4095)
			dtemp = 4095;
		if (dtemp < 0)
			dtemp = 0;
		vco = (int)dtemp;
		master_delay = (offset + sptr->strobe) * CYCLE - delay[0];
		}
	else {
		dtemp = DGAIN / sptr->pfac;
		if (sptr->vfreq > dtemp) {
			sptr->vfreq -= dtemp;
			sptr->slew--;
			while (sptr->slew < -CLOCK / 2) {
				sptr->slew += CLOCK;
				sptr->borrow -= CLOCK;
				}
			}
		if (sptr->vfreq < -dtemp) {
			sptr->vfreq += dtemp;
			sptr->slew++;
			while (sptr->slew >= CLOCK / 2) {
				sptr->slew -= CLOCK;
				sptr->borrow += CLOCK;
				}
			}
		}
	sptr->psnr += (phase - sptr->psnr) / (sptr->pfac * 8);
	sptr->system_delay = (offset + sptr->strobe) * CYCLE +
		sptr->slew * 1. / (CLOCK / CYCLE) - master_delay;
	while (sptr->system_delay < 0)
		sptr->system_delay += fri * CYCLE;
	}

/*
 * Subroutine station_type()
 *
 * During pulse-group search mode the first two of five stations operate
 * as masters, while the next two operate as slaves and the last idles.
 * In other modes the first station operates as the master, the rest as
 * slaves.
 *
 * Calling sequence: char = station_type(sptr)
 *	sptr		station structure pointer
 *	char		station type or 's' (slave)
 *
 * Variables and functions used
 *	sptr->type	station type 'm', 'w', 'x', 'y', 'z', ' '
 *	sptr->index	station index (1-5)
 */
char station_type(sptr)
	struct station *sptr;		/* station structure pointer */
	{

	if (sptr->mode == MODE_PSRCH)
		if (sptr->index < 2)
			return ('m');
		else
			return ('s');
	else
		return(sptr->type);
	}

/*
 * Subroutine local_time()
 *
 * This subroutine collects various intricate quantities to determine
 * the precise local time relative to the LORAN-C master station of the
 * selected chain. However, there is an ambiguity relative to the epoch
 * of the transmitted pulse group which can only be resolved by outside
 * means, such as using the time-of-coincidence (TOC) values published
 * by USNO or tracking more than one chain from the same station.
 *
 * The epoch variable represents the local time determined at the end of
 * the latest gri. The local time relative to the epoch of the master
 * station pulse group can be determined for all stations and the
 * differences used to adjust the receiver position to minimize residual
 * errors, although that trick hasn't happened yet. To find the actual
 * time, subtract the following corrections using the units specified:
 *
 * -800 (GRI)		correction to the beginning of the pulse-group gate
 *					(cycle)
 * sptr->strobe		interval to the strobe (reference) cycle (cycle)
 * sptr->slew		interval to the reference phase of the reference
 *					cycle (clock) (slave stations only). This varies
 *					normally over the range [-CLOCK, CLOCK - 1].
 * sptr->borrow		add to the above to correct for occasional wander
 *					beyond one cycle when close to a cycle boundary
 *					(cycle). This is normally either zero or +-CLOCK.
 * 30 (RCVDELAY)	cycle-strobe delay (us) to compensate for the
 *					receiver bandpass characteristic beyond the pulse-
 *					group gate.
 * rcvr_delay		antenna-receiver delay (us) to compensate for the
 *					antenna, feedline and rf-stage delays to the pulse-
 *					group gate. This is specified in the data file.
 * sptr->ems_delay	emission delay (us) (slave stations only)
 * sptr->path_delay	path delay from transmitter to receiver (us)
 *
 * Calling sequence: time = local_time(sptr)
 *	sptr		station structure pointer
 *
 * Variables and functions used (see above)
 */
double local_time(sptr)
	struct station *sptr;		/* station pointer */
	{
	return (epoch - sptr->ems_delay - sptr->path_delay - RCVDELAY -
		rcvr_delay - (GRI + sptr->strobe) * CLOCK - (sptr->slew +
		sptr->borrow) * 1. / (CLOCK / CYCLE));
	}

/* end program */