/Proj4/PJ_nsper.c

http://github.com/route-me/route-me · C · 152 lines · 147 code · 5 blank · 0 comment · 15 complexity · 48a66bc194b7d17201ee4e63be916f63 MD5 · raw file 

    

  1. #ifndef lint
    2. 

  2. static const char SCCSID[]="@(#)PJ_nsper.c	4.1	94/02/15	GIE	REL";
    3. 

  3. #endif
    4. 

  4. #define PROJ_PARMS__ \
    5. 

  5. 	double	height; \
    6. 

  6. 	double	sinph0; \
    7. 

  7. 	double	cosph0; \
    8. 

  8. 	double	p; \
    9. 

  9. 	double	rp; \
    10. 

  10. 	double	pn1; \
    11. 

  11. 	double	pfact; \
    12. 

  12. 	double	h; \
    13. 

  13. 	double	cg; \
    14. 

  14. 	double	sg; \
    15. 

  15. 	double	sw; \
    16. 

  16. 	double	cw; \
    17. 

  17. 	int		mode; \
    18. 

  18. 	int		tilt;
    19. 

  19. #define PJ_LIB__
    20. 

  20. #include	"projects.h"
    21. 

  21. PROJ_HEAD(nsper, "Near-sided perspective") "\n\tAzi, Sph\n\th=";
    22. 

  22. PROJ_HEAD(tpers, "Tilted perspective") "\n\tAzi, Sph\n\ttilt= azi= h=";
    23. 

  23. # define EPS10 1.e-10
    24. 

  24. # define N_POLE	0
    25. 

  25. # define S_POLE 1
    26. 

  26. # define EQUIT	2
    27. 

  27. # define OBLIQ	3
    28. 

  28. FORWARD(s_forward); /* spheroid */
    29. 

  29. 	double  coslam, cosphi, sinphi;
    30. 

  30. 

  31. 	sinphi = sin(lp.phi);
    32. 

  32. 	cosphi = cos(lp.phi);
    33. 

  33. 	coslam = cos(lp.lam);
    34. 

  34. 	switch (P->mode) {
    35. 

  35. 	case OBLIQ:
    36. 

  36. 		xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam;
    37. 

  37. 		break;
    38. 

  38. 	case EQUIT:
    39. 

  39. 		xy.y = cosphi * coslam;
    40. 

  40. 		break;
    41. 

  41. 	case S_POLE:
    42. 

  42. 		xy.y = - sinphi;
    43. 

  43. 		break;
    44. 

  44. 	case N_POLE:
    45. 

  45. 		xy.y = sinphi;
    46. 

  46. 		break;
    47. 

  47. 	}
    48. 

  48. 	if (xy.y < P->rp) F_ERROR;
    49. 

  49. 	xy.y = P->pn1 / (P->p - xy.y);
    50. 

  50. 	xy.x = xy.y * cosphi * sin(lp.lam);
    51. 

  51. 	switch (P->mode) {
    52. 

  52. 	case OBLIQ:
    53. 

  53. 		xy.y *= (P->cosph0 * sinphi -
    54. 

  54. 		   P->sinph0 * cosphi * coslam);
    55. 

  55. 		break;
    56. 

  56. 	case EQUIT:
    57. 

  57. 		xy.y *= sinphi;
    58. 

  58. 		break;
    59. 

  59. 	case N_POLE:
    60. 

  60. 		coslam = - coslam;
    61. 

  61. 	case S_POLE:
    62. 

  62. 		xy.y *= cosphi * coslam;
    63. 

  63. 		break;
    64. 

  64. 	}
    65. 

  65. 	if (P->tilt) {
    66. 

  66. 		double yt, ba;
    67. 

  67. 

  68. 		yt = xy.y * P->cg + xy.x * P->sg;
    69. 

  69. 		ba = 1. / (yt * P->sw * P->h + P->cw);
    70. 

  70. 		xy.x = (xy.x * P->cg - xy.y * P->sg) * P->cw * ba;
    71. 

  71. 		xy.y = yt * ba;
    72. 

  72. 	}
    73. 

  73. 	return (xy);
    74. 

  74. }
    75. 

  75. INVERSE(s_inverse); /* spheroid */
    76. 

  76. 	double  rh, cosz, sinz;
    77. 

  77. 

  78. 	if (P->tilt) {
    79. 

  79. 		double bm, bq, yt;
    80. 

  80. 

  81. 		yt = 1./(P->pn1 - xy.y * P->sw);
    82. 

  82. 		bm = P->pn1 * xy.x * yt;
    83. 

  83. 		bq = P->pn1 * xy.y * P->cw * yt;
    84. 

  84. 		xy.x = bm * P->cg + bq * P->sg;
    85. 

  85. 		xy.y = bq * P->cg - bm * P->sg;
    86. 

  86. 	}
    87. 

  87. 	rh = hypot(xy.x, xy.y);
    88. 

  88. 	if ((sinz = 1. - rh * rh * P->pfact) < 0.) I_ERROR;
    89. 

  89. 	sinz = (P->p - sqrt(sinz)) / (P->pn1 / rh + rh / P->pn1);
    90. 

  90. 	cosz = sqrt(1. - sinz * sinz);
    91. 

  91. 	if (fabs(rh) <= EPS10) {
    92. 

  92. 		lp.lam = 0.;
    93. 

  93. 		lp.phi = P->phi0;
    94. 

  94. 	} else {
    95. 

  95. 		switch (P->mode) {
    96. 

  96. 		case OBLIQ:
    97. 

  97. 			lp.phi = asin(cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh);
    98. 

  98. 			xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh;
    99. 

  99. 			xy.x *= sinz * P->cosph0;
    100. 

  100. 			break;
    101. 

  101. 		case EQUIT:
    102. 

  102. 			lp.phi = asin(xy.y * sinz / rh);
    103. 

  103. 			xy.y = cosz * rh;
    104. 

  104. 			xy.x *= sinz;
    105. 

  105. 			break;
    106. 

  106. 		case N_POLE:
    107. 

  107. 			lp.phi = asin(cosz);
    108. 

  108. 			xy.y = -xy.y;
    109. 

  109. 			break;
    110. 

  110. 		case S_POLE:
    111. 

  111. 			lp.phi = - asin(cosz);
    112. 

  112. 			break;
    113. 

  113. 		}
    114. 

  114. 		lp.lam = atan2(xy.x, xy.y);
    115. 

  115. 	}
    116. 

  116. 	return (lp);
    117. 

  117. }
    118. 

  118. FREEUP; if (P) pj_dalloc(P); }
    119. 

  119. 	static PJ *
    120. 

  120. setup(PJ *P) {
    121. 

  121. 	if ((P->height = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
    122. 

  122. 	if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
    123. 

  123. 		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
    124. 

  124. 	else if (fabs(P->phi0) < EPS10)
    125. 

  125. 		P->mode = EQUIT;
    126. 

  126. 	else {
    127. 

  127. 		P->mode = OBLIQ;
    128. 

  128. 		P->sinph0 = sin(P->phi0);
    129. 

  129. 		P->cosph0 = cos(P->phi0);
    130. 

  130. 	}
    131. 

  131. 	P->pn1 = P->height / P->a; /* normalize by radius */
    132. 

  132. 	P->p = 1. + P->pn1;
    133. 

  133. 	P->rp = 1. / P->p;
    134. 

  134. 	P->h = 1. / P->pn1;
    135. 

  135. 	P->pfact = (P->p + 1.) * P->h;
    136. 

  136. 	P->inv = s_inverse;
    137. 

  137. 	P->fwd = s_forward;
    138. 

  138. 	P->es = 0.;
    139. 

  139. 	return P;
    140. 

  140. }
    141. 

  141. ENTRY0(nsper)
    142. 

  142. 	P->tilt = 0;
    143. 

  143. ENDENTRY(setup(P))
    144. 

  144. ENTRY0(tpers)
    145. 

  145. 	double omega, gamma;
    146. 

  146. 

  147. 	omega = pj_param(P->params, "dtilt").f * DEG_TO_RAD;
    148. 

  148. 	gamma = pj_param(P->params, "dazi").f * DEG_TO_RAD;
    149. 

  149. 	P->tilt = 1;
    150. 

  150. 	P->cg = cos(gamma); P->sg = sin(gamma);
    151. 

  151. 	P->cw = cos(omega); P->sw = sin(omega);
    152. 

  152. ENDENTRY(setup(P))
    153.