/Proj4/nad_cvt.c

 1#ifndef lint
 2static const char SCCSID[]="@(#)nad_cvt.c	4.3   95/09/23 GIE REL";
 3#endif
 4#define PJ_LIB__
 5#include "projects.h"
 6#define MAX_TRY 9
 7#define TOL 1e-12
 8	LP
 9nad_cvt(LP in, int inverse, struct CTABLE *ct) {
10	LP t, tb;
11
12	if (in.lam == HUGE_VAL)
13		return in;
14	/* normalize input to ll origin */
15	tb = in;
16	tb.lam -= ct->ll.lam;
17	tb.phi -= ct->ll.phi;
18	tb.lam = adjlon(tb.lam - PI) + PI;
19	t = nad_intr(tb, ct);
20	if (inverse) {
21		LP del, dif;
22		int i = MAX_TRY;
23
24		if (t.lam == HUGE_VAL) return t;
25		t.lam = tb.lam + t.lam;
26		t.phi = tb.phi - t.phi;
27
28		do {
29			del = nad_intr(t, ct);
30
31                        /* This case used to return failure, but I have
32                           changed it to return the first order approximation
33                           of the inverse shift.  This avoids cases where the
34                           grid shift *into* this grid came from another grid.
35                           While we aren't returning optimally correct results
36                           I feel a close result in this case is better than
37                           no result.  NFW
38                           To demonstrate use -112.5839956 49.4914451 against
39                           the NTv2 grid shift file from Canada. */
40			if (del.lam == HUGE_VAL) 
41                        {
42                            if( getenv( "PROJ_DEBUG" ) != NULL )
43                                fprintf( stderr, 
44                                         "Inverse grid shift iteration failed, presumably at grid edge.\n"
45                                         "Using first approximation.\n" );
46                            /* return del */;
47                            break;
48                        }
49
50			t.lam -= dif.lam = t.lam - del.lam - tb.lam;
51			t.phi -= dif.phi = t.phi + del.phi - tb.phi;
52		} while (i-- && fabs(dif.lam) > TOL && fabs(dif.phi) > TOL);
53		if (i < 0) {
54                    if( getenv( "PROJ_DEBUG" ) != NULL )
55                        fprintf( stderr, 
56                                 "Inverse grid shift iterator failed to converge.\n" );
57                    t.lam = t.phi = HUGE_VAL;
58                    return t;
59		}
60		in.lam = adjlon(t.lam + ct->ll.lam);
61		in.phi = t.phi + ct->ll.phi;
62	} else {
63		if (t.lam == HUGE_VAL)
64			in = t;
65		else {
66			in.lam -= t.lam;
67			in.phi += t.phi;
68		}
69	}
70	return in;
71}