view release on metacpan or  search on metacpan

Index.xs  view on Meta::CPAN

		lat_0_idx = tmp;
	}
	
	/* Populate return values (all unsigned integers) */
	
	Inline_Stack_Reset;
	Inline_Stack_Push(sv_2mortal(newSVuv( grid_level )));
	Inline_Stack_Push(sv_2mortal(newSVuv( grid_size )));
	Inline_Stack_Push(sv_2mortal(newSVuv( max_grid_idx )));
	Inline_Stack_Push(sv_2mortal(newSVuv( lat_0_idx )));
	Inline_Stack_Push(sv_2mortal(newSVuv( lat_1_idx )));
	Inline_Stack_Push(sv_2mortal(newSVuv( lon_0_idx )));
	Inline_Stack_Push(sv_2mortal(newSVuv( lon_1_idx )));
	Inline_Stack_Done;
}




/*
   Compute distance between two points on a sphere
   
   diameter is in meters
   lat_0, lon_0, lat_1, and lon_1 are in radians
   Distance returned is in meters
   
   To compute a distance first call SetUpDistance(...) with the first point 
   then call HaversineDistance(...) to get the distance to a second point.
   
   The C version can use either floats or doubles whereas Perl uses doubles.  
   When using floats instead of doubles the loss of precision is typically less 
   than a meter (about 2 meters in the worst case).   On modern hardware there 
   should be little noticable difference between using floats and doubles.  On 
   older hardware or in embedded systems floats may give better performance.
*/


/* Functions using floats */

float f_diameter, f_lat_1, f_lon_1;
float f_cos_lat_1;

void SetUpDistance_float(float new_diameter, float new_lat_1, float new_lon_1) {
	f_diameter = new_diameter;
	f_lat_1 = new_lat_1;
	f_lon_1 = new_lon_1;
	f_cos_lat_1 = cosf( new_lat_1 );
}

float HaversineDistance_float(float lat_0, float lon_0) {
	float sin_lat_diff_over_2 = sinf( ( lat_0 - f_lat_1 ) / 2.0 );
	float sin_lon_diff_over_2 = sinf( ( lon_0 - f_lon_1 ) / 2.0 );
	
	float n = ( sin_lat_diff_over_2 * sin_lat_diff_over_2 ) 
	          + (
	              ( sin_lon_diff_over_2 * sin_lon_diff_over_2 )
	              * f_cos_lat_1
	              * cosf( lat_0 )
	            );
	
	/* The haversine formula may get messy around antipodal points so clip to the largest sane value. */
	if ( n < 0.0 ) { n = 0.0; }
	
	return f_diameter  * asinf( sqrtf(n) );
}



/* Functions using doubles */

double d_diameter, d_lat_1, d_lon_1;
double d_cos_lat_1;

void SetUpDistance_double(double new_diameter, double new_lat_1, double new_lon_1) {
	d_diameter = new_diameter;
	d_lat_1 = new_lat_1;
	d_lon_1 = new_lon_1;
	d_cos_lat_1 = cos( new_lat_1 );
}

double HaversineDistance_double(double lat_0, double lon_0) {
	double sin_lat_diff_over_2 = sin( ( lat_0 - d_lat_1 ) / 2.0 );
	double sin_lon_diff_over_2 = sin( ( lon_0 - d_lon_1 ) / 2.0 );
	
	double n = ( sin_lat_diff_over_2 * sin_lat_diff_over_2 ) 
	          + (
	              ( sin_lon_diff_over_2 * sin_lon_diff_over_2 )
	              * d_cos_lat_1
	              * cos( lat_0 )
	            );
	
	/* The haversine formula may get messy around antipodal points so clip to the largest sane value. */
	if ( n < 0.0 ) { n = 0.0; }
	
	return d_diameter  * asin( sqrt(n) );
}

/* End of C code */

MODULE = Geo::Index  PACKAGE = Geo::Index  

PROTOTYPES: DISABLE


unsigned int
GetCCodeVersion ()

unsigned int
fast_log2_double (n)
	double	n

unsigned int
fast_log2_float (n)
	float	n

void
ComputeAreaExtrema_float (tile_adjust, max_size, max_level, p_lat, p_lat_rad, p_lon, polar_circumference, search_radius)
	int	tile_adjust
	unsigned long	max_size
	unsigned int	max_level
	float	p_lat
	float	p_lat_rad
	float	p_lon
	float	polar_circumference
	float	search_radius
        PREINIT:
        I32* temp;
        PPCODE:
        temp = PL_markstack_ptr++;
        ComputeAreaExtrema_float(tile_adjust, max_size, max_level, p_lat, p_lat_rad, p_lon, polar_circumference, search_radius);
        if (PL_markstack_ptr != temp) {
          /* truly void, because dXSARGS not invoked */
          PL_markstack_ptr = temp;
          XSRETURN_EMPTY; /* return empty stack */
        }
        /* must have used dXSARGS; list context implied */
        return; /* assume stack size is correct */

void
ComputeAreaExtrema_double (tile_adjust, max_size, max_level, p_lat, p_lat_rad, p_lon, polar_circumference, search_radius)
	int	tile_adjust
	unsigned long	max_size
	unsigned int	max_level
	double	p_lat
	double	p_lat_rad
	double	p_lon
	double	polar_circumference
	double	search_radius
        PREINIT:
        I32* temp;
        PPCODE:
        temp = PL_markstack_ptr++;