#define _GNU_SOURCE
/* Tests for simultaneous transformation of pm and pos */

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

/* positions in rad */
extern void tpmpm_(double *alpha, double *delta, double *pma, double *pmd, 
	double *foreshort, double *srcEp, double *destEp, double *alphaOut, 
	double *deltaOut, double *pmaOut, double *pmdOut);


int testCase(double alpha, double delta, double pma, double pmd, double fore,
	double srcEp, double destEp)
{
	usPos pos={alpha,delta}, pm={pma, pmd};
	usPos posFromC, pmFromC, posFromFortran, pmFromFortran;

	xpa_movePmAndPm(pos, pm, fore, srcEp, destEp, &posFromC, &pmFromC);
	/* the fortran code does the *cos(delta) in the function */
	pm.alpha /= cos(pos.delta); 
	tpmpm_(&(pos.alpha), &(pos.delta), &(pm.alpha), &(pm.delta), &fore,
		&srcEp, &destEp, &(posFromFortran.alpha), &(posFromFortran.delta),
		&(pmFromFortran.alpha), &(pmFromFortran.delta));
	if (!xpa_usPosEqual(posFromFortran, posFromC)) {
		fprintf(stderr, "Pos failure for %g %g (%g %g), %g %g %g\n",
			alpha, delta, pma, pmd, fore, srcEp, destEp);
		return 1;
	}
	pmFromFortran.alpha *= cos(posFromFortran.delta);
	if (!xpa_usPosEqual(pmFromFortran, pmFromC)) {
		fprintf(stderr, "PM failure for %g %g (%g %g), %g %g %g\n",
			alpha, delta, pma, pmd, fore, srcEp, destEp);
		return 1;
	}
	return 0;
}

int main(void)
{
	double alpha, delta;
	int pmFactor;
	int res=0;

	for (alpha=0; alpha<=2*M_PI; alpha+=M_PI/10) {
		for (delta=-M_PI+1e-4; delta<M_PI-1e-10; delta+=M_PI/10) {
			for (pmFactor=1; pmFactor<100; pmFactor*=2) {
				res |= testCase(alpha, delta, .001*pmFactor, -0.001*pmFactor, 
					0, 2000.0, 2110.0);
				res |= testCase(alpha, delta, -0.001*pmFactor, -0.001*pmFactor, 
					0, 2000.0, 2110.0);
				res |= testCase(alpha, delta, -0.001*pmFactor, 0.001*pmFactor, 
					0, 2000.0, 2110.0);
				res |= testCase(alpha, delta, 0.001*pmFactor, 0.001*pmFactor, 
					0, 2000.0, 2110.0);
			}
		}
	}
	res |= testCase(0.23, -1.2, 0.0, 0.0, 0, 2000.0, 2110.0);
	res |= testCase(0.23, -1.2, 0.0, 0.0, 0, 2110.0, 200.0);
	res |= testCase(0.23, -1.2, 0.0, 0.0, 0.01, 2110.0, 200.0);
	res |= testCase(0.23, -1.2, 0.02, -0.02, 0.01, 2110.0, 200.0);
	return res;
}