#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

int ncalls=0;

int miss(const gsl_vector * v, void * params, gsl_vector * f){
	ncalls++;
	const double v0=150;
	const double g=9.8;
	double theta = gsl_vector_get(v,0);
	double phi = gsl_vector_get(v,1);
	double* p = (double*)params;
	double xtarget = p[0];
	double ytarget = p[1];
	double t=2*v0*cos(theta)/g;
	double x=v0*sin(theta)*sin(phi)*t;
	double y=v0*sin(theta)*cos(phi)*t;
	gsl_vector_set(f,0,x-xtarget);
	gsl_vector_set(f,1,y-ytarget);
return GSL_SUCCESS;
}

int main(void){
	double p[2]={500,-400};
	gsl_multiroot_function F;
	F.f=&miss;
	F.n=2;
	F.params=(void*)p;

	gsl_multiroot_fsolver *s =
		gsl_multiroot_fsolver_alloc(gsl_multiroot_fsolver_hybrids,2);

	gsl_vector *v = gsl_vector_alloc(2);
	gsl_vector_set(v,0,1);
	gsl_vector_set(v,1,1);
	gsl_multiroot_fsolver_set(s,&F,v);

	int status,iter=0;
	do{
		iter++;
		status = gsl_multiroot_fsolver_iterate (s);
		if (status)   /* check if solver is stuck */ break;

	printf("%i %g\n",iter,
	sqrt( pow(gsl_vector_get(s->f,0),2)+ pow(gsl_vector_get(s->f,1),2)));

		status = gsl_multiroot_test_residual (s->f, 1e-6);
	}while (status == GSL_CONTINUE && iter < 1000);

fprintf(stderr,"ncalls=%i\n",ncalls);
	
return 0;
}