Lensing Example using Healpy

#include
#include

int main(int argc, char **argv)
{
	if (argc !=2)
	{
		printf("  Please provide the path to a lens as argument");
		return -1;
	}
	parsec::MagneticLens lens;
	lens.loadLens(argv[1]);
  lens.normalizeLens();

	size_t counter = 0;

  double energy = 30; // EeV

	std::printf("# Energy [EeV], lon_in [rad] lat_in [rad] lon_out [rad] lat_out [rad]\n");
	while (counter < 1000)
	{
		double lon_in = -M_PI + 2*M_PI * float(rand()) / RAND_MAX; 
		double lat_in = asin(2*float(rand())/RAND_MAX-1);
    
		double lon_out = lon_in; 
		double lat_out = lat_in ;
		if (lens.transformCosmicRay(energy, lon_out, lat_out));
		{
      std::printf("%8.4f  %8.4f  %8.4f  %8.4f\n",lon_in, lat_in,
			lon_out, lat_out);
      counter+=1;
		}
	
	} 
}