« back to all changes in this revision
Viewing changes to src/PJ_healpix.c
- browse files at revision 11
- compare with another revision
- download diff
- download tarball
- view history from revision 11
- Committer: Package Import Robot
- Author(s): Francesco Paolo Lovergine, Jerome Villeneuve Larouche, Francesco Paolo Lovergine
- Date: 2013-11-25 15:11:25 UTC
- mfrom: (1.2.7)
- Revision ID: package-import@ubuntu.com-20131125151125-mvcw144wvgep1hev
Tags: 4.8.0-1
[ Jerome Villeneuve Larouche ]
* New upstream release
* Modified libproj-dev.install to remove nad_list.h and projects.h
* Modified proj-bin.install to remove nad2nad
* Modified proj-bin.manpages to remove nad2nad man
* Added symbols for libproj0
[ Francesco Paolo Lovergine ]
* Properly merged with current git master and sid version 4.7.0-2.
* Removed proj transitional package and obsolete conflicts/replaces against
series 4.6 and older.
* New upstream release
* Modified libproj-dev.install to remove nad_list.h and projects.h
* Modified proj-bin.install to remove nad2nad
* Modified proj-bin.manpages to remove nad2nad man
* Added symbols for libproj0
[ Francesco Paolo Lovergine ]
* Properly merged with current git master and sid version 4.7.0-2.
* Removed proj transitional package and obsolete conflicts/replaces against
series 4.6 and older.
- files added:
- .pc
- files removed:
- man/man1/nad2nad.1
- files modified:
- ChangeLog
added
removed
src/PJ_healpix.c
1
/******************************************************************************
2
* $Id: PJ_healpix.c 1504 2011-10-18 14:58:57Z landcare $
3
*
4
* Project: PROJ.4
5
* Purpose: Implementation of the healpix projection.
6
* Definition: http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/scenzgrid.pdf
7
* Author: Alex Raichev & Michael Speth , spethm@landcareresearch.co.nz
8
*
9
******************************************************************************
10
* Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
11
*
12
* Permission is hereby granted, free of charge, to any person obtaining a
13
* copy of this software and associated documentation files (the "Software"),
14
* to deal in the Software without restriction, including without limitation
15
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
16
* and/or sell copies of the Software, and to permit persons to whom the
17
* Software is furnished to do so, subject to the following conditions:
18
*
19
* The above copyright notice and this permission notice shall be included
20
* in all copies or substantial portions of the Software.
21
*
22
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
23
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
24
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
25
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
26
* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
27
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
28
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
29
* SOFTWARE.
30
*****************************************************************************/
31
32
#define PROJ_PARMS__ \
33
int npole;\
34
int spole;
35
36
#define PJ_LIB__
37
# include <projects.h>
38
PROJ_HEAD(healpix, "HEALPix") "\n\tSph., Ellps.";
39
PROJ_HEAD(rhealpix, "rHEALPix") "\n\tSph., Ellps.\n\tnpole= spole=";
40
# include <stdio.h>
41
# define R1 {{ 0,-1},{ 1, 0}} /** Matrix for anticlockwise rotation by pi/2 **/
42
# define R2 {{-1, 0},{ 0,-1}} /** Matrix for anticlockwise rotation by pi (R1 X R1) X = dot product **/
43
# define R3 {{ 0, 1},{-1, 0}} /** Matrix for anticlockwise rotation by 3*pi/2 (R2 X R1) **/
44
# define IDENT {{1,0},{0,1}}
45
/**
46
* 0 - Identity matrix<br>
47
* 1 - Counter-clockwise rotation by PI/2<br>
48
* 2 - Counter-clockwise rotation by PI<br>
49
* 3 - Counter-clockwise rotation by 3*PI/2<br>
50
* 4 - Counter-clockwise rotation by 3*PI/2<br>
51
* 5 - Counter-clockwise rotation by PI<br>
52
* 6 - Counter-clockwise rotation by PI/2<br>
53
**/
54
# define ROT { IDENT, R1, R2, R3, R3, R2, R1}
55
# define RFACTOR 3 /** Used for returning the rotation matrix **/
56
/** Used for calculating if a point is within the HEALPix projection for sphere. **/
57
# define EPS 1e-12
58
typedef struct {
59
int cn; // the number 0 -> 4 indicating the position of the polar cap.
60
double x,y; // the coordinates of the pole points (point of most extreme latitude on the polar caps).
61
enum Region { north, south, equatorial } region;
62
} CapMap;
63
typedef struct {
64
double x,y;
65
} Point;
66
double rot[7][2][2] = ROT;
67
68
/**
69
NOTES: Alex Raichev implemented the math in python and this is a port of his work.
70
The healpix projection is a Lambert cylindrical equal-area projection for
71
equaltorial latitudes and an interrupted Colignon projection for polar
72
latitudes.
73
**/
74
75
/**
76
* Returns the sign of the double.
77
* @param v the parameter whose sign is returned.
78
* @return 1 for positive number, -1 for negative, and 0 for zero.
79
**/
80
double sign (double v) {
81
return v > 0 ? 1 : (v < 0 ? -1 : 0);
82
}
83
/**
84
* Scales the number by a factor.
85
* @param num the number to be scaled.
86
* @param factor the factor to scale the number by.
87
* @param isInverse 1 for scaling the number by 1 / factor and 0 for scaling by the factor.
88
* @return the scaled number.
89
**/
90
double scale_number(double num, double factor, int isInverse){
91
if(isInverse == 1){
92
return num * 1.0/factor;
93
}
94
return num * factor;
95
}
96
/**
97
* Scales all the items of the array by a factor.
98
* @param xy
99
**/
100
void scale_array(XY *array, double k, int inverse){
101
double c = 0;
102
if (inverse == 1) {
103
c = 1.0/k;
104
}else{
105
c = k;
106
}
107
array->x *= c;
108
array->y *= c;
109
}
110
/**
111
* Given an angle return its equivalent angle.
112
* @param x the angle to convert
113
* @return the equivalent angle such that -PI <= the angle returend <= PI
114
**/
115
double standardize_lon(double x){
116
if(x < -1*PI || x >= PI){
117
x = x - 2*PI*floor(x/(2*PI));
118
if(x >= PI){
119
x = x - 2*PI;
120
}
121
}
122
return x;
123
}
124
/**
125
* Given an angle, return its unit-circle equivalent angle.
126
* @param x the angel to convert.
127
* @return the equivalent angle such that -PI/2 <= the angle returned <= PI/2.
128
**/
129
double standardize_lat(double x){
130
if( x < -PI/2.0 || x > PI/2){
131
x = x-2.0*PI*floor(x/(2.0*PI));
132
if(x > PI/2.0 && x <= 3.0*PI/2){
133
x = PI - x;
134
}else{
135
x = x - 2*PI;
136
}
137
}
138
return x;
139
}
140
/**
141
* Calculates if the point lies on or within the polygon.
142
* Very good explination of how this works: http://paulbourke.net/geometry/insidepoly/
143
* @param nvert the number of vertices in the polygon.
144
* @param vert the x,y-coordinates of the polygon's vertices
145
* @param testx the x-coordinate of the test point.
146
* @param testy the y-coordinate of the test point.
147
* @return 1 if on or within the bounds of the polygon, and 0 otherwise.
148
**/
149
static
150
int pnpoly(int nvert, double vert[][2], double testx, double testy){
151
152
int i,j,c = 0;
153
int counter = 0;
154
double xinters;
155
Point p1,p2;
156
157
// check for boundrary cases
158
for(i = 0; i < nvert; i++){
159
if(testx == vert[i][0] && testy == vert[i][1]){
160
return 1;
161
}
162
}
163
164
// initialize p1
165
p1.x = vert[0][0];
166
p1.y = vert[0][1];
167
168
for(i = 1; i < nvert; i++){
169
p2.x = vert[i % nvert][0];
170
p2.y = vert[i % nvert][1];
171
172
if(testy > MIN(p1.y,p2.y)){
173
if (testy <= MAX(p1.y,p2.y)) {
174
if (testx <= MAX(p1.x,p2.x)) {
175
if (p1.y != p2.y) {
176
xinters = (testy-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
177
if (p1.x == p2.x || testx <= xinters){
178
counter++;
179
}
180
}
181
}
182
}
183
}
184
p1 = p2;
185
}
186
if(counter % 2 == 0){
187
return 0;
188
}else{
189
return 1;
190
}
191
return c;
192
}
193
/**
194
* Calculates if the coordinates are within the image of projection.
195
* @param x the x-coordinate to check.
196
* @param y the y-coordinate to check.
197
* @param proj 0 for healpix and 1 for rhealpix.
198
* @param npole the positions of the polar squares, only used for rhealpix.
199
* @param spole the positions of the polar squares, only used for rhealpix.
200
* @return 1 if the coordinate is within the projection and 0 otherwise.
201
**/
202
int in_image(double x, double y, int proj, int npole, int spole){
203
if(proj == 0){
204
double healpixVertsJit[][2] = {
205
{-1.0*PI-EPS ,PI/4.0},
206
{-3.0*PI/4.0 ,PI/2.0+EPS},
207
{-1.0*PI/2.0 ,PI/4.0+EPS},
208
{-1.0*PI/4.0 ,PI/2.0+EPS},
209
{0.0 ,PI/4.0+EPS},
210
{PI/4.0 ,PI/2.0+EPS},
211
{PI/2.0 ,PI/4.0+EPS},
212
{3.0*PI/4.0 ,PI/2.0+EPS},
213
{PI+EPS ,PI/4.0},
214
{PI+EPS ,-1.0*PI/4.0},
215
{3.0*PI/4.0 ,-1.0*PI/2.0-EPS},
216
{PI/2.0 ,-1.0*PI/4.0-EPS},
217
{PI/4.0 ,-1.0*PI/2.0-EPS},
218
{0.0 ,-1.0*PI/4.0-EPS},
219
{-1.0*PI/4.0 ,-1.0*PI/2.0-EPS},
220
{-1.0*PI/2.0 ,-1.0*PI/4.0-EPS},
221
{-3.0*PI/4.0 ,-1.0*PI/2.0-EPS},
222
{-1.0*PI-EPS ,-1.0*PI/4.0}};
223
return pnpoly((int)sizeof(healpixVertsJit)/sizeof(healpixVertsJit[0]),
224
healpixVertsJit,x,y);
225
}else{
226
// Used for calculating if a point is within the rHEALPix projection for sphere.
227
double rhealpixVertsJit[][2] = {
228
{-1.0*PI-EPS ,PI/4.0+EPS},
229
{-1.0*PI + npole*PI/2.0-EPS ,PI/4.0+EPS},
230
{-1.0*PI + npole*PI/2.0-EPS ,3*PI/4.0+EPS},
231
{-1.0*PI + (npole + 1.0)*PI/2.0+EPS ,3*PI/4.0+EPS},
232
{-1.0*PI + (npole + 1.0)*PI/2.0+EPS ,PI/4.0+EPS},
233
{PI+EPS ,PI/4.0+EPS},
234
{PI+EPS ,-1.0*PI/4.0-EPS},
235
{-1.0*PI + (spole + 1.0)*PI/2.0+EPS ,-1.0*PI/4.0-EPS},
236
{-1.0*PI + (spole + 1.0)*PI/2.0+EPS ,-3.0*PI/4.0-EPS},
237
{-1.0*PI + spole*PI/2.0-EPS ,-3.0*PI/4.0-EPS},
238
{-1.0*PI + spole*PI/2.0-EPS ,-1.0*PI/4.0-EPS},
239
{-1.0*PI-EPS ,-1.0*PI/4.0-EPS}};
240
return pnpoly((int)sizeof(rhealpixVertsJit)/sizeof(rhealpixVertsJit[0]),
241
rhealpixVertsJit,x,y);
242
}
243
}
244
/**
245
* Returns an authalic latitude of the point given a point of geographic
246
* latitude phi on an ellipse of eccentricity e.
247
* pj_authlat is the inverse of the alex's auth_lat.
248
* @param phi
249
* @param e
250
* @param inverse 1 for inverse or 0 otherwise.
251
* @return the authalic latitude of the point.
252
**/
253
double auth_lat(double phi, double e, int inverse){
254
if(inverse == 0){
255
double q_numerator = ((1.0 - pow(e,2.0)) * sin(phi));
256
double q_demonitor = (1.0 - (pow(e*sin(phi),2.0)));
257
double q_subtractor = - (1.0 - pow(e,2.0)) / (2.0*e) * log((1.0 - e*sin(phi)) / (1.0+e*sin(phi)));
258
double q = ((1.0 - pow(e,2.0)) * sin(phi)) / (1.0 - (pow(e*sin(phi),2.0))) -
259
(1.0 - pow(e,2.0)) / (2.0*e) * log((1.0 - e*sin(phi)) / (1.0+e*sin(phi)));
260
261
double qp = 1.0 - (1.0-pow(e,2.0)) / (2.0*e)*log((1.0 - e) / (1.0 + e));
262
double ratio = q/qp;
263
// Rounding errors
264
if( fabsl(ratio) > 1){
265
ratio = sign(ratio);
266
}
267
return asin(ratio);
268
}
269
return phi + (pow(e,2) / 3.0 + 31*pow(e,4) / 180.0 + 517.0*pow(e,6)/5040.0) * sin(2.0*phi)
270
+ (23.0*pow(e,4)/360.0 + 251.0*pow(e,6)/3780.0)*sin(4.0*phi)
271
+ 761.0*pow(e,6)/45360.0 * sin(6.0*phi);
272
}
273
/**
274
* Compute the forward signature functions of the HEALPix
275
* projection of a sphere with radius `R` and central meridian `lon0`.
276
**/
277
XY healpix_sphere(LP lp, PJ *P){
278
double lam = standardize_lon(lp.lam);
279
double phi = standardize_lat(lp.phi);
280
double phi0 = aasin(P->ctx, 2.0/3.0);
281
XY xy;
282
// equatorial region
283
if( fabsl(phi) <= phi0) {
284
xy.x = lam;
285
xy.y = 3.0*PI/8.0*sin(phi);
286
} else {
287
double lamc;
288
double sigma = sqrt(3.0 * (1 - fabsl(sin(phi))));
289
double cn = floor(2 * lam / PI + 2);
290
if (cn >= 4) {
291
cn = 3;
292
}
293
lamc = -3*PI/4 + (PI/2)*cn;
294
xy.x = lamc + (lam - lamc) * sigma;
295
xy.y = sign(phi)*PI/4 * (2 - sigma);
296
}
297
xy.x = scale_number(xy.x,P->a,0);
298
xy.y = scale_number(xy.y,P->a,0);
299
return xy;
300
}
301
/**
302
* Compute the inverse signature functions of the HEALPix
303
* projection of a sphere with radius `R` and central meridian `lon0`.
304
**/
305
LP healpix_sphere_inv(XY xy, PJ *P){
306
double x,y,y0;
307
double cn;
308
double xc;
309
double tau;
310
LP lp;
311
// Scale down to radius 1 sphere
312
x = scale_number(xy.x,P->a,1);
313
y = scale_number(xy.y,P->a,1);
314
y0 = PI/4.0;
315
// Equatorial region.
316
if(fabsl(y) <= y0){
317
lp.lam = x;
318
lp.phi = asin(8.0*y/(3.0*PI));
319
} else if(fabsl(y) < PI/2.0){
320
cn = floor(2.0 * x/PI + 2.0);
321
if(cn >= 4){
322
cn = 3;
323
}
324
xc = -3.0 * PI/4.0 + (PI/2.0)*cn;
325
tau = 2.0 - 4.0*fabsl(y)/PI;
326
lp.lam = xc + (x - xc)/tau;
327
lp.phi = sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
328
} else {
329
lp.lam = -1.0*PI - P->lam0;
330
lp.phi = sign(y)*PI/2.0;
331
}
332
return (lp);
333
}
334
/**
335
* Adds one vector to another of length 2.
336
* @param a the first term.
337
* @param b the second term.
338
* @param ret holds the summation of the vectors.
339
**/
340
static void vector_add(double a[], double b[],double * ret){
341
int i;
342
for(i = 0; i < 2; i++){
343
ret[i] = a[i] + b[i];
344
}
345
}
346
/**
347
* Subs tracts one vector from another of length 2.
348
* @param a the minuend.
349
* @param b the subtrahend.
350
* @param ret the difference of the vectors where the difference is the result of a minus b.
351
**/
352
static void vector_sub(double a[], double b[], double * ret){
353
int i;
354
for(i = 0; i < 2; i++){
355
ret[i] = a[i] - b[i];
356
}
357
}
358
/**
359
* Calculates the dot product of the arrays.
360
* @param a the array that will be used to calculate the dot product.
361
* Must contain the same number of columns as b's rows. Must be a matrix with equal lengthed rows and columns.
362
* @param b the array that will be used to calculate the dot product; must contain the same number of rows as a's columns.
363
* @param length the size of the b array. Note, a's column size must equal b's length.
364
* @param ret the dot product of a and b.
365
**/
366
static void dot_product(double a[2][2], double b[], double * ret){
367
int i,j;
368
int length = 2;
369
for(i = 0; i < length; i++){
370
ret[i] = 0;
371
for(j = 0; j < length; j++){
372
ret[i] += a[i][j]*b[i];
373
}
374
}
375
}
376
/**
377
* Returns the polar cap number, pole point coordinates, and region
378
* for x,y in the HEALPix projection of the sphere of radius R.
379
* @param x coordinate in the HEALPix or rHEALPix.
380
* @param y coordinate in the HEALPix or rHEALPix.
381
* @param npole integer between 0 and 3 indicating the position of the north pole.
382
* @param spole integer between 0 and 3 indicating teh position of the south pole.
383
* @param inverse 1 computes the rHEALPix projection and 0 computes forward.
384
* @return a structure containing the cap poles.
385
**/
386
static CapMap get_cap(double x, double y, double R, int npole, int spole, int inverse){
387
CapMap capmap;
388
double c;
389
390
capmap.x = x;
391
capmap.y = y;
392
393
if(inverse == 0){
394
if(y > R*PI/4.0){
395
capmap.region = north;
396
c = R*PI/2.0;
397
}else if(y < -1*R*PI/4.0){
398
capmap.region = south;
399
c = -1*R*PI/2.0;
400
}else{
401
capmap.region = equatorial;
402
capmap.cn = 0;
403
return capmap;
404
}
405
// polar region
406
if(x < -1*R*PI/2.0){
407
capmap.cn = 0;
408
capmap.x = (-1*R*3.0*PI/4.0);
409
capmap.y = c;
410
}else if(x >= -1*R*PI/2.0 && x < 0){
411
capmap.cn = 1;
412
capmap.x = -1*R*PI/4.0;
413
capmap.y = c;
414
}else if(x >= 0 && x < R*PI/2.0){
415
capmap.cn = 2;
416
capmap.x = R*PI/4.0;
417
capmap.y = c;
418
}else{
419
capmap.cn = 3;
420
capmap.x = R*3.0*PI/4.0;
421
capmap.y = c;
422
}
423
return capmap;
424
}else{
425
double c;
426
double eps;
427
if(y > R*PI/4.0){
428
capmap.region = north;
429
capmap.x = -1*R*3.0*PI/4.0 + npole*R*PI/2.0;
430
capmap.y = R*PI/2.0;
431
x = x - npole*R*PI/2.0;
432
}else if(y < -1*R*PI/4.0){
433
capmap.region = south;
434
capmap.x = -1*R*3.0*PI/4.0 + spole*R*PI/2;
435
capmap.y = -1*R*PI/2.0;
436
x = x - spole*R*PI/2.0;
437
}else{
438
capmap.region = equatorial;
439
capmap.cn = 0;
440
return capmap;
441
}
442
// Polar Region, find # of HEALPix polar cap number that
443
// x,y moves to when rHEALPix polar square is disassembled.
444
eps = R*1e-15; // Kludge. Fuzz to avoid some rounding errors.
445
if(capmap.region == north){
446
if(y >= -1*x - R*PI/4.0 - eps && y < x + R*5.0*PI/4.0 - eps){
447
capmap.cn = 1;
448
}else if(y > -1*x -1*R*PI/4.0 + eps && y >= x + R*5.0*PI/4.0 - eps){
449
capmap.cn = 2;
450
}else if(y <= -1*x -1*R*PI/4.0 + eps && y > x + R*5.0*PI/4.0 + eps){
451
capmap.cn = 3;
452
}else{
453
capmap.cn = 0;
454
}
455
}else if(capmap.region == south){
456
if(y <= x + R*PI/4.0 + eps && y > -1*x - R*5.0*PI/4 + eps){
457
capmap.cn = 1;
458
}else if(y < x + R*PI/4.0 - eps && y <= -1*x - R*5.0*PI/4.0 + eps){
459
capmap.cn = 2;
460
}else if(y >= x + R*PI/4.0 - eps && y < -1*x - R*5.0*PI/4.0 - eps){
461
capmap.cn = 3;
462
}else {
463
capmap.cn = 0;
464
}
465
}
466
return capmap;
467
}
468
}
469
/**
470
* Rearrange point x,y in the HEALPix projection by
471
* combining the polar caps into two polar squares.
472
* Put the north polar square in position npole and
473
* the south polar square in position spole.
474
* @param x coordinate in the HEALPix projection of the sphere.
475
* @param y coordinate in the HEALPix projection of the sphere.
476
* @param R - the Sphere's radius.
477
* @param npole integer between 0 and 3 indicating the position
478
* of the north polar square.
479
* @param spole integer between 0 and 3 indicating the position
480
* of the south polar square.
481
* @param inverse 1 to uncombine the polar caps and 0 to combine.
482
**/
483
static XY combine_caps(double x, double y, double R, int npole, int spole, int inverse){
484
XY xy;
485
double v[2];
486
double a[2];
487
double vector[2];
488
double tmpVect[2];
489
double v_min_c[2];
490
double ret_dot[2];
491
double ret_add[2];
492
CapMap capmap = get_cap(x,y,R,npole,spole,inverse);
493
494
if(capmap.region == equatorial){
495
xy.x = capmap.x;
496
xy.y = capmap.y;
497
return xy;
498
}
499
v[0] = x;
500
v[1] = y;
501
if(inverse == 0){
502
// compute forward function by rotating, translating, and shifting xy.
503
int pole = 0;
504
double (*tmpRot)[2];
505
double c[2] = {capmap.x,capmap.y};
506
if(capmap.region == north){
507
pole = npole;
508
tmpRot = rot[capmap.cn];
509
a[0] = R*-3.0*PI/4.0;
510
a[1] = PI/2.0;
511
}else {
512
pole = spole;
513
tmpRot = rot[capmap.cn+RFACTOR];
514
a[0] = R*-3.0*PI/4.0;
515
a[1] = PI/-2.0;
516
}
517
518
tmpVect[0] = R*pole*PI/2.0;
519
tmpVect[1] = 0;
520
// translate, rotate, then shift
521
vector_sub(v,c,v_min_c);
522
dot_product(tmpRot,v_min_c,ret_dot);
523
vector_add(a,tmpVect,ret_add);
524
vector_add(ret_dot, ret_add, vector);
525
xy.x = vector[0];
526
xy.y = vector[1];
527
return xy;
528
}else{
529
// compute inverse function.
530
// get the current position of rHEALPix polar squares
531
int pole = floor( (capmap.x + R*3.0*PI/4.0) / (R*PI/2.0));
532
double tmpVect[2] = {R*pole*PI/2.0,0};
533
double coord[2] = {x,y};
534
double (*tmpRot)[2];
535
int cn;
536
// translate polar square to position 0
537
vector_sub(coord,tmpVect,v);
538
// disassemble
539
if(capmap.region == north){
540
cn = capmap.cn + RFACTOR;
541
a[0] = R*-3*PI/4.0;
542
a[1] = PI/2.0;
543
}else{
544
cn = capmap.cn;
545
a[0] = R*-3*PI/4.0;
546
a[1] = PI/-2.0;
547
}
548
tmpVect[0] = R*capmap.cn*PI/2.0;
549
tmpVect[1] = 0;
550
// Math: Rotate Matrix * v-a + a + R*CN*{PI/2,0}
551
vector_sub(v,a,v_min_c);
552
dot_product(rot[cn],v_min_c,ret_dot);
553
vector_add(ret_dot,a,ret_add);
554
vector_add(ret_add,tmpVect,vector);
555
xy.x = vector[0];
556
xy.y = vector[1];
557
return xy;
558
}
559
}
560
FORWARD(e_healpix_forward); /* ellipsoidal */
561
//int r1[][2] = R1;
562
double bet = auth_lat(lp.phi, P->e, 0);
563
lp.phi = bet;
564
P->a = P->ra;
565
return healpix_sphere(lp,P);
566
}
567
FORWARD(s_healpix_forward); /* spheroid */
568
return healpix_sphere(lp, P);
569
}
570
INVERSE(e_healpix_inverse); /* ellipsoidal */
571
double bet, x, y;
572
P->a = P->ra;
573
574
// Scale down to radius 1 sphere before checking x,y
575
x = scale_number(xy.x,P->a,1);
576
y = scale_number(xy.y,P->a,1);
577
// check if the point is in the image
578
if(in_image(x,y,0,0,0) == 0){
579
lp.lam = HUGE_VAL;
580
lp.phi = HUGE_VAL;
581
pj_ctx_set_errno( P->ctx, -15);
582
return lp;
583
}
584
585
lp = healpix_sphere_inv(xy, P);
586
587
lp.phi = auth_lat(lp.phi,P->e,1);
588
589
return (lp);
590
}
591
INVERSE(s_healpix_inverse); /* spheroid */
592
double x = xy.x;
593
double y = xy.y;
594
// Scale down to radius 1 sphere before checking x,y
595
x = scale_number(x,P->a,1);
596
y = scale_number(y,P->a,1);
597
// check if the point is in the image
598
if(in_image(x,y,0,0,0) == 0){
599
lp.lam = HUGE_VAL;
600
lp.phi = HUGE_VAL;
601
pj_ctx_set_errno( P->ctx, -15);
602
return lp;
603
}
604
return healpix_sphere_inv(xy, P);
605
}
606
FORWARD(e_rhealpix_forward); /* ellipsoidal */
607
double bet = auth_lat(lp.phi,P->e,0);
608
lp.phi = bet;
609
xy = healpix_sphere(lp,P);
610
return combine_caps(xy.x, xy.y, P->a, P->npole, P->spole, 0);
611
}
612
FORWARD(s_rhealpix_forward); /* spheroid */
613
// Compute forward function.
614
xy = healpix_sphere(lp,P);
615
return combine_caps(xy.x, xy.y, P->a, P->npole, P->spole, 0);
616
}
617
INVERSE(e_rhealpix_inverse); /* ellipsoidal */
618
double x = scale_number(xy.x,P->a,1);
619
double y = scale_number(xy.y,P->a,1);
620
// check for out of bounds coordinates
621
if(in_image(x,y,1,P->npole,P->spole) == 0){
622
lp.lam = HUGE_VAL;
623
lp.phi = HUGE_VAL;
624
pj_ctx_set_errno( P->ctx, -15);
625
return lp;
626
}
627
628
xy = combine_caps(xy.x,xy.y,P->a,P->npole,P->spole,1);
629
lp = healpix_sphere_inv(xy, P);
630
lp.phi = auth_lat(lp.phi,P->e,1);
631
return lp;
632
}
633
INVERSE(s_rhealpix_inverse); /* spheroid */
634
double x = scale_number(xy.x,P->a,1);
635
double y = scale_number(xy.y,P->a,1);
636
// check for out of bounds coordinates
637
if(in_image(x,y,1,P->npole,P->spole) == 0){
638
lp.lam = HUGE_VAL;
639
lp.phi = HUGE_VAL;
640
pj_ctx_set_errno( P->ctx, -15);
641
return lp;
642
}
643
xy = combine_caps(xy.x,xy.y,P->a,P->npole,P->spole,1);
644
return healpix_sphere_inv(xy, P);
645
}
646
FREEUP;
647
if (P) {
648
pj_dalloc(P);
649
}
650
}
651
ENTRY0(healpix)
652
if(P->es){
653
P->inv = e_healpix_inverse; P->fwd = e_healpix_forward;
654
}else{
655
P->inv = s_healpix_inverse; P->fwd = s_healpix_forward;
656
}
657
ENDENTRY(P)
658
ENTRY0(rhealpix)
659
P->npole = pj_param(P->ctx, P->params,"inpole").i;
660
P->spole = pj_param(P->ctx,P->params,"ispole").i;
661
662
// check for valid npole and spole inputs
663
if(P->npole < 0 || P->npole > 3){
664
E_ERROR(-47);
665
}
666
if(P->spole < 0 || P->spole > 3){
667
E_ERROR(-47);
668
}
669
670
if(P->es){
671
P->inv = e_rhealpix_inverse; P->fwd = e_rhealpix_forward;
672
}else{
673
P->inv = s_rhealpix_inverse; P->fwd = s_rhealpix_forward;
674
}
675
ENDENTRY(P)
Loggerhead is a web-based interface for Breezy
Version: 2.0.1