/Proj4/nad_intr.c

http://github.com/route-me/route-me · C · 65 lines · 63 code · 1 blank · 1 comment · 18 complexity · 612a4dae9ba79d99761175d3a9ec1337 MD5 · raw file 

    

  1. /* Determine nad table correction value */
    2. 

  2. #ifndef lint
    3. 

  3. static const char SCCSID[]="@(#)nad_intr.c	4.2   95/09/23 GIE REL";
    4. 

  4. #endif
    5. 

  5. #define PJ_LIB__
    6. 

  6. #include "projects.h"
    7. 

  7. 	LP
    8. 

  8. nad_intr(LP t, struct CTABLE *ct) {
    9. 

  9. 	LP val, frct;
    10. 

  10. 	ILP indx;
    11. 

  11. 	double m00, m10, m01, m11;
    12. 

  12. 	FLP *f00, *f10, *f01, *f11;
    13. 

  13. 	long index;
    14. 

  14. 	int in;
    15. 

  15. 

  16. 	indx.lam = floor(t.lam /= ct->del.lam);
    17. 

  17. 	indx.phi = floor(t.phi /= ct->del.phi);
    18. 

  18. 	frct.lam = t.lam - indx.lam;
    19. 

  19. 	frct.phi = t.phi - indx.phi;
    20. 

  20. 	val.lam = val.phi = HUGE_VAL;
    21. 

  21. 	if (indx.lam < 0) {
    22. 

  22. 		if (indx.lam == -1 && frct.lam > 0.99999999999) {
    23. 

  23. 			++indx.lam;
    24. 

  24. 			frct.lam = 0.;
    25. 

  25. 		} else
    26. 

  26. 			return val;
    27. 

  27. 	} else if ((in = indx.lam + 1) >= ct->lim.lam) {
    28. 

  28. 		if (in == ct->lim.lam && frct.lam < 1e-11) {
    29. 

  29. 			--indx.lam;
    30. 

  30. 			frct.lam = 1.;
    31. 

  31. 		} else
    32. 

  32. 			return val;
    33. 

  33. 	}
    34. 

  34. 	if (indx.phi < 0) {
    35. 

  35. 		if (indx.phi == -1 && frct.phi > 0.99999999999) {
    36. 

  36. 			++indx.phi;
    37. 

  37. 			frct.phi = 0.;
    38. 

  38. 		} else
    39. 

  39. 			return val;
    40. 

  40. 	} else if ((in = indx.phi + 1) >= ct->lim.phi) {
    41. 

  41. 		if (in == ct->lim.phi && frct.phi < 1e-11) {
    42. 

  42. 			--indx.phi;
    43. 

  43. 			frct.phi = 1.;
    44. 

  44. 		} else
    45. 

  45. 			return val;
    46. 

  46. 	}
    47. 

  47. 	index = indx.phi * ct->lim.lam + indx.lam;
    48. 

  48. 	f00 = ct->cvs + index++;
    49. 

  49. 	f10 = ct->cvs + index;
    50. 

  50. 	index += ct->lim.lam;
    51. 

  51. 	f11 = ct->cvs + index--;
    52. 

  52. 	f01 = ct->cvs + index;
    53. 

  53. 	m11 = m10 = frct.lam;
    54. 

  54. 	m00 = m01 = 1. - frct.lam;
    55. 

  55. 	m11 *= frct.phi;
    56. 

  56. 	m01 *= frct.phi;
    57. 

  57. 	frct.phi = 1. - frct.phi;
    58. 

  58. 	m00 *= frct.phi;
    59. 

  59. 	m10 *= frct.phi;
    60. 

  60. 	val.lam = m00 * f00->lam + m10 * f10->lam +
    61. 

  61. 			  m01 * f01->lam + m11 * f11->lam;
    62. 

  62. 	val.phi = m00 * f00->phi + m10 * f10->phi +
    63. 

  63. 			  m01 * f01->phi + m11 * f11->phi;
    64. 

  64. 	return val;
    65. 

  65. }
    66.