#define _GNU_SOURCE

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

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


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, posFromFortran;

	posFromC = xpa_movePm(pos, pm, fore, srcEp, destEp);
	/* the fortran code does the *cos(delta) in the function */
	pm.alpha /= cos(pos.delta); 
	tpmpos_(&(pos.alpha), &(pos.delta), &(pm.alpha), &(pm.delta), &fore,
		&srcEp, &destEp, &(posFromFortran.alpha), &(posFromFortran.delta));
	if (!xpa_usPosEqual(posFromFortran, posFromC)) {
		fprintf(stderr, "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;

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