File: gnomo.c

package info (click to toggle)
pgsphere 1.1.1%2B2020-10-20-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 3,700 kB
  • sloc: ansic: 12,032; sql: 6,091; cpp: 853; makefile: 216; perl: 168; yacc: 145; xml: 66; lex: 55; sh: 1
file content (70 lines) | stat: -rw-r--r-- 1,790 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include <postgres.h>
#include <fmgr.h>

#include <utils/geo_decls.h> /* Point */

#include <point.h> /* SPoint */
#include <gnomo.h>

#include <math.h>

PG_FUNCTION_INFO_V1(gnomonic_proj);
PG_FUNCTION_INFO_V1(gnomonic_inv);

/*
 * Direct gnomonic projection.
 *
 * point gnomonic_proj(spoint spherical_point, spoint tangential_point)
 */
Datum gnomonic_proj(PG_FUNCTION_ARGS)
{
	Point*  g = (Point*)  palloc(sizeof(Point));
	SPoint* p = (SPoint*) PG_GETARG_POINTER(0);
	SPoint* t = (SPoint*) PG_GETARG_POINTER(1);

	double delta_lng = p->lng - t->lng;
	double sin_lat = sin(p->lat);
	double cos_lat = cos(p->lat);
	double sin_lat_t = sin(t->lat);
	double cos_lat_t = cos(t->lat);
	double sin_delta = sin(delta_lng);
	double cos_delta = cos(delta_lng);

	double cos_lat__cos_delta = cos_lat * cos_delta;
	double cos_dist = sin_lat_t * sin_lat + cos_lat_t * cos_lat__cos_delta;

	g->x = cos_lat * sin_delta / cos_dist;
	g->y = (cos_lat_t * sin_lat - sin_lat_t * cos_lat__cos_delta) / cos_dist;

	if (p->lng == t->lng && p->lat == t->lat)
	{
		g->x = 0;
		g->y = 0;
	}
	PG_RETURN_POINTER(g);
}

/*
 * Inverse gnomonic projection.
 *
 * spoint gnomonic_inv(point plane_point, spoint tangential_point)
 */
Datum gnomonic_inv(PG_FUNCTION_ARGS)
{
	SPoint* p = (SPoint*) palloc(sizeof(SPoint));
	Point*  g = (Point*)  PG_GETARG_POINTER(0);
	SPoint* t = (SPoint*) PG_GETARG_POINTER(1);

	double rho_sq = g->x * g->x + g->y * g->y;
	double rho = sqrt(rho_sq);
	double cos_c = 1 / sqrt(1 + rho_sq);
	double sin_c = 1 / sqrt(1 + 1 / rho_sq);
	double cos_lat_t = cos(t->lat);
	double sin_lat_t = sin(t->lat);

	p->lng = t->lng + atan2(g->x * sin_c,
							rho * cos_lat_t * cos_c - g->y * sin_lat_t * sin_c);
	p->lat = asin(cos_c * sin_lat_t + g->y *sin_c * cos_lat_t / rho);

	PG_RETURN_POINTER(p);
}