Skip to content

Commit 91ba576

Browse files
jswhitjswhit
committed
bug fixes from proj4 svn
1 parent ef5c76b commit 91ba576

File tree

2 files changed

+90
-50
lines changed

2 files changed

+90
-50
lines changed

‎nad2bin.c

+19-7
Original file line numberDiff line numberDiff line change
@@ -77,7 +77,7 @@ int main(int argc, char **argv) {
7777
/* ==================================================================== */
7878
for( i = 1; i < argc; i++ )
7979
{
80-
if( strcmp(argv[i],"-f") && i < argc-1 )
80+
if( strcmp(argv[i],"-f") == 0 && i < argc-1 )
8181
{
8282
format = argv[++i];
8383
}
@@ -98,9 +98,15 @@ int main(int argc, char **argv) {
9898
/* Read the ASCII Table */
9999
/* ==================================================================== */
100100

101-
fgets(ct.id, MAX_TAB_ID, stdin);
102-
scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
103-
&ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi);
101+
if ( NULL == fgets(ct.id, MAX_TAB_ID, stdin) ) {
102+
perror("fgets");
103+
exit(1);
104+
}
105+
if ( EOF == scanf("%d %d %*d %lf %lf %lf %lf", &ct.lim.lam, &ct.lim.phi,
106+
&ct.ll.lam, &ct.del.lam, &ct.ll.phi, &ct.del.phi) ) {
107+
perror("scanf");
108+
exit(1);
109+
}
104110
if (!(ct.cvs = (FLP *)malloc(tsize = ct.lim.lam * ct.lim.phi *
105111
sizeof(FLP)))) {
106112
perror("mem. alloc");
@@ -112,7 +118,10 @@ int main(int argc, char **argv) {
112118
ct.del.phi *= DEG_TO_RAD;
113119
/* load table */
114120
for (p = ct.cvs, i = 0; i < ct.lim.phi; ++i) {
115-
scanf("%d:%ld %ld", &ichk, &laml, &phil);
121+
if ( EOF == scanf("%d:%ld %ld", &ichk, &laml, &phil) ) {
122+
perror("scanf on row");
123+
exit(1);
124+
}
116125
if (ichk != i) {
117126
fprintf(stderr,"format check on row\n");
118127
exit(1);
@@ -121,7 +130,10 @@ int main(int argc, char **argv) {
121130
t.phi = phil * U_SEC_TO_RAD;
122131
*p++ = t;
123132
for (j = 1; j < ct.lim.lam; ++j) {
124-
scanf("%ld %ld", &lam, &phi);
133+
if ( EOF == scanf("%ld %ld", &lam, &phi) ) {
134+
perror("scanf on column");
135+
exit(1);
136+
}
125137
t.lam = (laml += lam) * U_SEC_TO_RAD;
126138
t.phi = (phil += phi) * U_SEC_TO_RAD;
127139
*p++ = t;
@@ -183,7 +195,7 @@ int main(int argc, char **argv) {
183195
{
184196
swap_words( header + 96, 8, 4 );
185197
swap_words( header + 128, 4, 2 );
186-
swap_words( ct.cvs, 4, ct.lim.lam * ct.lim.phi );
198+
swap_words( ct.cvs, 4, ct.lim.lam * 2 * ct.lim.phi );
187199
}
188200

189201
if( fwrite( header, sizeof(header), 1, fp ) != 1 ) {

‎src/PJ_healpix.c

+71-43
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,30 @@ double standardize_lat(double x){
137137
}
138138
return x;
139139
}
140+
/**
141+
* Returns the index of the 2d array in rot.
142+
* @param index range from -3 to 3.
143+
* @return the index into the rot 3d array.
144+
*/
145+
static int get_rotate_index(int index){
146+
switch(index){
147+
case 0:
148+
return 0;
149+
case 1:
150+
return 1;
151+
case 2:
152+
return 2;
153+
case 3:
154+
return 3;
155+
case -1:
156+
return 4;
157+
case -2:
158+
return 5;
159+
case -3:
160+
return 6;
161+
}
162+
return 0;
163+
}
140164
/**
141165
* Calculates if the point lies on or within the polygon.
142166
* Very good explination of how this works: http://paulbourke.net/geometry/insidepoly/
@@ -369,7 +393,7 @@ static void dot_product(double a[2][2], double b[], double * ret){
369393
for(i = 0; i < length; i++){
370394
ret[i] = 0;
371395
for(j = 0; j < length; j++){
372-
ret[i] += a[i][j]*b[i];
396+
ret[i] += a[i][j]*b[j];
373397
}
374398
}
375399
}
@@ -426,12 +450,12 @@ static CapMap get_cap(double x, double y, double R, int npole, int spole, int in
426450
double eps;
427451
if(y > R*PI/4.0){
428452
capmap.region = north;
429-
capmap.x = -1*R*3.0*PI/4.0 + npole*R*PI/2.0;
453+
capmap.x = R*(-3.0*PI/4.0 + npole*PI/2.0);
430454
capmap.y = R*PI/2.0;
431455
x = x - npole*R*PI/2.0;
432456
}else if(y < -1*R*PI/4.0){
433457
capmap.region = south;
434-
capmap.x = -1*R*3.0*PI/4.0 + spole*R*PI/2;
458+
capmap.x = R*(-3.0*PI/4.0 + spole*PI/2);
435459
capmap.y = -1*R*PI/2.0;
436460
x = x - spole*R*PI/2.0;
437461
}else{
@@ -444,23 +468,23 @@ static CapMap get_cap(double x, double y, double R, int npole, int spole, int in
444468
eps = R*1e-15; // Kludge. Fuzz to avoid some rounding errors.
445469
if(capmap.region == north){
446470
if(y >= -1*x - R*PI/4.0 - eps && y < x + R*5.0*PI/4.0 - eps){
447-
capmap.cn = 1;
471+
capmap.cn = (npole + 1) % 4;
448472
}else if(y > -1*x -1*R*PI/4.0 + eps && y >= x + R*5.0*PI/4.0 - eps){
449-
capmap.cn = 2;
473+
capmap.cn = (npole + 2) % 4;
450474
}else if(y <= -1*x -1*R*PI/4.0 + eps && y > x + R*5.0*PI/4.0 + eps){
451-
capmap.cn = 3;
475+
capmap.cn = (npole + 3) % 4;
452476
}else{
453-
capmap.cn = 0;
477+
capmap.cn = npole;
454478
}
455479
}else if(capmap.region == south){
456480
if(y <= x + R*PI/4.0 + eps && y > -1*x - R*5.0*PI/4 + eps){
457-
capmap.cn = 1;
481+
capmap.cn = (spole + 1) % 4;
458482
}else if(y < x + R*PI/4.0 - eps && y <= -1*x - R*5.0*PI/4.0 + eps){
459-
capmap.cn = 2;
483+
capmap.cn = (spole + 2) % 4;
460484
}else if(y >= x + R*PI/4.0 - eps && y < -1*x - R*5.0*PI/4.0 - eps){
461-
capmap.cn = 3;
485+
capmap.cn = (spole + 3) % 4;
462486
}else {
463-
capmap.cn = 0;
487+
capmap.cn = spole;
464488
}
465489
}
466490
return capmap;
@@ -505,53 +529,57 @@ static XY combine_caps(double x, double y, double R, int npole, int spole, int i
505529
double c[2] = {capmap.x,capmap.y};
506530
if(capmap.region == north){
507531
pole = npole;
508-
tmpRot = rot[capmap.cn];
509-
a[0] = R*-3.0*PI/4.0;
510-
a[1] = PI/2.0;
532+
a[0] = R * (-3.0*PI/4.0 + pole * PI/2);
533+
a[1] = R * (PI/2.0 + pole * 0);
534+
535+
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
536+
vector_sub(v,c,v_min_c);
537+
dot_product(tmpRot,v_min_c,ret_dot);
538+
vector_add(ret_dot, a, vector);
539+
511540
}else {
512541
pole = spole;
513-
tmpRot = rot[capmap.cn+RFACTOR];
514-
a[0] = R*-3.0*PI/4.0;
515-
a[1] = PI/-2.0;
542+
a[0] = R * (-3.0*PI/4.0 + pole * PI/2);
543+
a[1] = R * (PI/-2.0 + pole * 0);
544+
545+
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
546+
vector_sub(v,c,v_min_c);
547+
dot_product(tmpRot,v_min_c,ret_dot);
548+
549+
vector_add(ret_dot, a, vector);
516550
}
517551

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);
525552
xy.x = vector[0];
526553
xy.y = vector[1];
527554
return xy;
528555
}else{
529556
// compute inverse function.
530557
// 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];
535558
int cn;
536-
// translate polar square to position 0
537-
vector_sub(coord,tmpVect,v);
559+
int pole = 0;
560+
double (*tmpRot)[2];
561+
double c[2] = {capmap.x,capmap.y};
538562
// disassemble
539563
if(capmap.region == north){
540-
cn = capmap.cn + RFACTOR;
541-
a[0] = R*-3*PI/4.0;
542-
a[1] = PI/2.0;
564+
pole = npole;
565+
a[0] = R * (-3.0*PI/4.0 + capmap.cn * PI/2);
566+
a[1] = R * (PI/2.0 + capmap.cn * 0);
567+
568+
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
569+
vector_sub(v,c,v_min_c);
570+
dot_product(tmpRot,v_min_c,ret_dot);
571+
vector_add(ret_dot, a, vector);
572+
543573
}else{
544-
cn = capmap.cn;
545-
a[0] = R*-3*PI/4.0;
546-
a[1] = PI/-2.0;
574+
pole = spole;
575+
a[0] = R * (-3.0*PI/4.0 + capmap.cn * PI/2);
576+
a[1] = R * (PI/-2.0 + capmap.cn * 0);
577+
578+
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
579+
vector_sub(v,c,v_min_c);
580+
dot_product(tmpRot,v_min_c,ret_dot);
581+
vector_add(ret_dot, a, vector);
547582
}
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);
555583
xy.x = vector[0];
556584
xy.y = vector[1];
557585
return xy;

0 commit comments

Comments
 (0)