/Proj4/nad_intr.c
C | 65 lines | 63 code | 1 blank | 1 comment | 18 complexity | 612a4dae9ba79d99761175d3a9ec1337 MD5 | raw file
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}