/Proj4/nad_intr.c

 1/* Determine nad table correction value */
 2#ifndef lint
 3static const char SCCSID[]="@(#)nad_intr.c	4.2   95/09/23 GIE REL";
 4#endif
 5#define PJ_LIB__
 6#include "projects.h"
 7	LP
 8nad_intr(LP t, struct CTABLE *ct) {
 9	LP val, frct;
10	ILP indx;
11	double m00, m10, m01, m11;
12	FLP *f00, *f10, *f01, *f11;
13	long index;
14	int in;
15
16	indx.lam = floor(t.lam /= ct->del.lam);
17	indx.phi = floor(t.phi /= ct->del.phi);
18	frct.lam = t.lam - indx.lam;
19	frct.phi = t.phi - indx.phi;
20	val.lam = val.phi = HUGE_VAL;
21	if (indx.lam < 0) {
22		if (indx.lam == -1 && frct.lam > 0.99999999999) {
23			++indx.lam;
24			frct.lam = 0.;
25		} else
26			return val;
27	} else if ((in = indx.lam + 1) >= ct->lim.lam) {
28		if (in == ct->lim.lam && frct.lam < 1e-11) {
29			--indx.lam;
30			frct.lam = 1.;
31		} else
32			return val;
33	}
34	if (indx.phi < 0) {
35		if (indx.phi == -1 && frct.phi > 0.99999999999) {
36			++indx.phi;
37			frct.phi = 0.;
38		} else
39			return val;
40	} else if ((in = indx.phi + 1) >= ct->lim.phi) {
41		if (in == ct->lim.phi && frct.phi < 1e-11) {
42			--indx.phi;
43			frct.phi = 1.;
44		} else
45			return val;
46	}
47	index = indx.phi * ct->lim.lam + indx.lam;
48	f00 = ct->cvs + index++;
49	f10 = ct->cvs + index;
50	index += ct->lim.lam;
51	f11 = ct->cvs + index--;
52	f01 = ct->cvs + index;
53	m11 = m10 = frct.lam;
54	m00 = m01 = 1. - frct.lam;
55	m11 *= frct.phi;
56	m01 *= frct.phi;
57	frct.phi = 1. - frct.phi;
58	m00 *= frct.phi;
59	m10 *= frct.phi;
60	val.lam = m00 * f00->lam + m10 * f10->lam +
61			  m01 * f01->lam + m11 * f11->lam;
62	val.phi = m00 * f00->phi + m10 * f10->phi +
63			  m01 * f01->phi + m11 * f11->phi;
64	return val;
65}