1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * gmt_proj.c contains the specific projection functions that convert
19 * longitude and latitude to x, y. These are used via function pointers
20 * set in gmt_map.c
21 *
22 * Map_projections include functions that will set up the transformation
23 * between xy and latlon for several map projections.
24 *
25 * A few of the core coordinate transformation functions are based on similar
26 * FORTRAN routines written by Pat Manley, Doug Shearer, and Bill Haxby, and
27 * have been rewritten in C and subsequently streamlined. The Lambert conformal
28 * was originally coded up by Bernie Coakley. The rest is coded by Wessel/Smith
29 * based on P. Snyder, "Map Projections - a working manual", USGS Prof paper 1395.
30 *
31 * Transformations supported (both forward and inverse):
32 *
33 * Non-geographic projections:
34 *
35 * Linear x/y[/z] scaling
36 * Polar (theta, r) scaling
37 *
38 * Map projections:
39 *
40 * Cylindrical:
41 * Mercator
42 * Cassini Cylindrical
43 * Cylindrical Stereographic (e.g., Gall, B.S.A.M)
44 * Miller Cylindrical
45 * Oblique Mercator
46 * TM Transverse Mercator (Ellipsoidal and Spherical)
47 * UTM Universal Transverse Mercator
48 * Cylindrical Equal-area (e.g., Gall, Behrmann)
49 * Cylindrical Equidistant (e.g., Plate Carree)
50 * Conic:
51 * Albers Equal-Area Conic
52 * Lambert Conformal Conic
53 * Equidistant Conic
54 * Azimuthal:
55 * Stereographic Conformal
56 * Lambert Azimuthal Equal-Area
57 * Orthographic
58 * Azimuthal Equidistant
59 * Gnomonic
60 * Thematic:
61 * Mollweide Equal-Area
62 * Hammer-Aitoff Equal-Area
63 * Sinusoidal Equal-Area
64 * Winkel Tripel
65 * Robinson
66 * Eckert IV
67 * Eckert IV
68 * Van der Grinten
69 *
70 * Author: Paul Wessel
71 * Date: 1-JAN-2010
72 * Version: 5.x
73 */
74
75 #include "gmt_dev.h"
76 #include "gmt_internals.h"
77
78 #define GMT_PROJ_MAX_ITERATIONS 200
79 #define GMT_PROJ_CONV_LIMIT 1e-9
80 #define gmt_M_proj_is_zero(x) (fabs (x) < GMT_PROJ_CONV_LIMIT)
81
gmtproj_robinson_spline(struct GMT_CTRL * GMT,double xp,double * x,double * y,double * c)82 GMT_LOCAL double gmtproj_robinson_spline (struct GMT_CTRL *GMT, double xp, double *x, double *y, double *c) {
83 /* Returns the interpolated value y(xp) from the Robinson coefficients */
84
85 int j = 0, j1;
86 double yp, a, b, h, ih, dx;
87
88 if (xp < x[0] || xp > x[GMT_N_ROBINSON-1]) /* Desired point outside data range */
89 return (GMT->session.d_NaN);
90
91 while (j < GMT_N_ROBINSON && x[j] <= xp) j++;
92 if (j == GMT_N_ROBINSON) j--;
93 if (j > 0) j--;
94
95 dx = xp - x[j];
96 switch (GMT->current.setting.interpolant) { /* gmtproj_vrobinson would not allow case 0 so only GMT_SPLINE_AKIMA | GMT_SPLINE_CUBIC is possible */
97 case GMT_SPLINE_AKIMA:
98 yp = ((c[3*j+2]*dx + c[3*j+1])*dx + c[3*j])*dx + y[j];
99 break;
100 case GMT_SPLINE_CUBIC:
101 j1 = j + 1;
102 h = x[j1] - x[j];
103 ih = 1.0 / h;
104 a = (x[j1] - xp) * ih;
105 b = dx * ih;
106 yp = a * y[j] + b * y[j1] + ((a*a*a - a) * c[j] + (b*b*b - b) * c[j1]) * (h*h) / 6.0;
107 break;
108 default:
109 yp = 0;
110 }
111 return (yp);
112 }
113
gmtproj_check_R_J(struct GMT_CTRL * GMT,double * clon)114 GMT_LOCAL void gmtproj_check_R_J (struct GMT_CTRL *GMT, double *clon) /* Make sure -R and -J agree for global plots; J given priority */ {
115 double lon0 = 0.5 * (GMT->common.R.wesn[XLO] + GMT->common.R.wesn[XHI]);
116
117 if (GMT->current.map.is_world && lon0 != *clon) {
118 GMT->common.R.wesn[XLO] = *clon - 180.0;
119 GMT->common.R.wesn[XHI] = *clon + 180.0;
120 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian set with -J (%g) implies -R%g/%g/%g/%g\n",
121 *clon, GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
122 }
123 else if (!GMT->current.map.is_world) {
124 lon0 = *clon - 360.0;
125 while (lon0 < GMT->common.R.wesn[XLO]) lon0 += 360.0;
126 if (lon0 > GMT->common.R.wesn[XHI]) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Central meridian outside region\n");
127 }
128 }
129
gmtproj_iwinkel_sub(struct GMT_CTRL * GMT,double y,double * phi)130 GMT_LOCAL void gmtproj_iwinkel_sub (struct GMT_CTRL *GMT, double y, double *phi) {
131 /* Valid only along meridian 180 degree from central meridian. Used in left/right_winkel only */
132 int n_iter = 0;
133 double c, phi0, delta, sp, cp;
134
135 c = 2.0 * y * GMT->current.proj.i_EQ_RAD;
136 *phi = y * GMT->current.proj.i_EQ_RAD;
137 do {
138 phi0 = *phi;
139 sincos (phi0, &sp, &cp);
140 *phi = phi0 - (phi0 + M_PI_2 * sp - c) / (1.0 + M_PI_2 * cp);
141 delta = fabs (*phi - phi0);
142 n_iter++;
143 }
144 while (delta > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
145 *phi *= R2D;
146 }
147
gmtproj_ipolyconic_sub(struct GMT_CTRL * GMT,double y,double lon,double * x)148 GMT_LOCAL void gmtproj_ipolyconic_sub (struct GMT_CTRL *GMT, double y, double lon, double *x) {
149 /* Used in left/right_polyconic only */
150 double E, sp, cp, phi0, phi, delta;
151 int n_iter = 0;
152
153 *x = lon;
154 gmt_M_wind_lon (GMT, *x);
155 y *= GMT->current.proj.i_EQ_RAD;
156 y += GMT->current.proj.pole * D2R;
157 if (gmt_M_proj_is_zero (y))
158 *x *= GMT->current.proj.EQ_RAD * D2R;
159 else {
160 phi = y;
161 do {
162 phi0 = phi;
163 sincos (phi, &sp, &cp);
164 E = (*x) * sp;
165 cp /= sp; /* = cot(phi) */
166 phi = y - cp * (1.0 - cosd(E));
167 delta = fabs (phi - phi0);
168 n_iter++;
169 }
170 while (delta > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
171 *x = GMT->current.proj.EQ_RAD * cp * sin(E);
172 }
173 }
174
175 /* conversion from geodetic latitude to geocentric latitude */
176
gmtproj_genper_getgeocentric(struct GMT_CTRL * GMT,double phi,double h)177 GMT_LOCAL double gmtproj_genper_getgeocentric (struct GMT_CTRL *GMT, double phi, double h) {
178 double phig, sphi, cphi, N1;
179
180 sincosd (phi, &sphi, &cphi);
181
182 N1 = GMT->current.proj.EQ_RAD/sqrt(1.0 - (GMT->current.proj.ECC2*sphi*sphi));
183
184 phig = phi - asind(N1*GMT->current.proj.ECC2*sphi*cphi/((h/GMT->current.proj.EQ_RAD+1.0)*GMT->current.proj.EQ_RAD));
185
186 return (phig);
187 }
188
gmtproj_genper_toxy(struct GMT_CTRL * P,double lat,double lon,double h,double * x,double * y)189 GMT_LOCAL void gmtproj_genper_toxy (struct GMT_CTRL *P, double lat, double lon, double h, double *x, double *y) {
190 double angle;
191 double xp, yp, rp;
192 double N, C, S, K;
193 double sphi, cphi;
194 double sdphi, cdphi;
195 double sphi1, cphi1;
196 double sdlon, cdlon;
197
198 cdphi = P->current.proj.g_cdphi;
199 sdphi = P->current.proj.g_sdphi;
200
201 cphi1 = P->current.proj.g_cphi1;
202 sphi1 = P->current.proj.g_sphi1;
203 h *= 1e3;
204
205 sincosd (lat, &sphi, &cphi);
206
207 N = P->current.proj.g_R/sqrt(1.0 - (P->current.proj.ECC2*sphi*sphi));
208 C = ((N+h)/P->current.proj.g_R)*cphi;
209 S = ((N*P->current.proj.one_m_ECC2 + h)/P->current.proj.g_R)*sphi;
210
211 sincosd (lon - P->current.proj.g_lon0, &sdlon, &cdlon);
212
213 K = P->current.proj.g_H / (P->current.proj.g_P*cdphi - S*sphi1 - C*cphi1*cdlon);
214
215 xp = K*C*sdlon;
216 yp = K*(P->current.proj.g_P*sdphi + S*cphi1 - C*sphi1*cdlon);
217
218 rp = sqrt(xp*xp + yp*yp);
219
220 if (rp > P->current.proj.g_rmax) {
221 angle = atan2(xp, yp);
222 sincos (angle, &xp, &yp);
223 xp *= P->current.proj.g_rmax;
224 yp *= P->current.proj.g_rmax;
225 }
226 *x = xp;
227 *y = yp;
228
229 if (P->current.proj.g_debug > 1) {
230 GMT_Report (P->parent, GMT_MSG_DEBUG, "\n");
231 GMT_Report (P->parent, GMT_MSG_DEBUG, "lat %12.3f\n", lat);
232 GMT_Report (P->parent, GMT_MSG_DEBUG, "lon %12.3f\n", lon);
233 GMT_Report (P->parent, GMT_MSG_DEBUG, "h %12.3f\n", h);
234 GMT_Report (P->parent, GMT_MSG_DEBUG, "N %12.1f\n", N);
235 GMT_Report (P->parent, GMT_MSG_DEBUG, "C %12.7f\n", C);
236 GMT_Report (P->parent, GMT_MSG_DEBUG, "S %12.7f\n", S);
237 GMT_Report (P->parent, GMT_MSG_DEBUG, "K %12.1f\n", K);
238 GMT_Report (P->parent, GMT_MSG_DEBUG, "x %12.1f\n", *x);
239 GMT_Report (P->parent, GMT_MSG_DEBUG, "y %12.1f\n", *y);
240 }
241 }
242
gmtproj_genper_to_xtyt(struct GMT_CTRL * GMT,double angle,double x,double y,double offset,double * xt,double * yt)243 GMT_LOCAL void gmtproj_genper_to_xtyt (struct GMT_CTRL *GMT, double angle, double x, double y, double offset, double *xt, double *yt) {
244 double A, theta, xp, yp;
245
246 theta = GMT->current.proj.g_azimuth - angle;
247
248 A = (y * GMT->current.proj.g_cos_azimuth + x * GMT->current.proj.g_sin_azimuth) * GMT->current.proj.g_sin_tilt / GMT->current.proj.g_H + GMT->current.proj.g_cos_tilt;
249
250 if (A > 0.0) {
251 xp = (x * GMT->current.proj.g_cos_azimuth - y * GMT->current.proj.g_sin_azimuth) * GMT->current.proj.g_cos_tilt / A;
252 yp = (y * GMT->current.proj.g_cos_azimuth + x * GMT->current.proj.g_sin_azimuth) / A;
253 if (fabs(yp) > fabs(GMT->current.proj.g_max_yt)) {
254 yp = -GMT->current.proj.g_max_yt;
255 xp = -yp * tand(theta);
256 }
257 }
258 else {
259 yp = -GMT->current.proj.g_max_yt;
260 xp = -yp * tand(theta);
261 }
262
263 yp -= offset;
264
265 *xt = xp * GMT->current.proj.g_cos_twist - yp * GMT->current.proj.g_sin_twist;
266 *yt = yp * GMT->current.proj.g_cos_twist + xp * GMT->current.proj.g_sin_twist;
267
268 return;
269 }
270
gmtproj_genper_tolatlong(struct GMT_CTRL * GMT,double x,double y,double h,double * lat,double * lon)271 GMT_LOCAL int gmtproj_genper_tolatlong (struct GMT_CTRL *GMT, double x, double y, double h, double *lat, double *lon) {
272 double P, H, B, D;
273 double u, v, t, Kp, X, Y;
274 double E, S;
275 double phi_last;
276 double Kp2;
277 double phi, sphi;
278 double cphig;
279 double e2, R, one_m_e2;
280 double cphi1, sphi1;
281
282 int niter;
283 int set_exit = 0;
284
285 h *= 1e3;
286
287 H = GMT->current.proj.g_H;
288 P = GMT->current.proj.g_P;
289 R = GMT->current.proj.g_R;
290
291 one_m_e2 = GMT->current.proj.one_m_ECC2;
292 e2 = GMT->current.proj.ECC2;
293
294 cphig = GMT->current.proj.g_cphig;
295 cphi1 = GMT->current.proj.g_cphi1;
296 sphi1 = GMT->current.proj.g_sphi1;
297
298 B = GMT->current.proj.g_B;
299 D = GMT->current.proj.g_D;
300
301 u = GMT->current.proj.g_BLH - GMT->current.proj.g_DG*y + GMT->current.proj.g_BJ*y + GMT->current.proj.g_DHJ;
302 v = GMT->current.proj.g_LH2 + GMT->current.proj.g_G*y*y - GMT->current.proj.g_HJ*y + one_m_e2*x*x;
303
304 if (GMT->current.proj.g_debug > 1) {
305 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong - 1\n");
306 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: x %12.1f\n", x);
307 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: y %12.1f\n", y);
308 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: \n");
309 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: P %12.7f\n", P);
310 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phig %12.7f\n", GMT->current.proj.g_phig);
311 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: \n");
312 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: B %12.7f\n", B);
313 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: D %12.7f\n", D);
314 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: u %12.1f\n", u);
315 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: v %12.6e\n", v);
316 }
317 E = 1;
318
319 t = P*P*(1.0 - e2*cphig*cphig) - E*one_m_e2;
320 Kp2 = (1.0 - 4.0*(t/u)*(v/u));
321 if (Kp2 < 0.0)
322 Kp = -u/(2.0*t);
323 else
324 Kp = (-u + sqrt(u*u-4.0*t*v))/(2.0*t);
325 X = R*((B-H/Kp)*cphi1 - (y/Kp - D)*sphi1);
326 Y = R*x/Kp;
327 S = (y/Kp-D)*cphi1 + (B-H/Kp)*sphi1;
328
329 if (gmt_M_is_dnan(Kp) || gmt_M_is_dnan(X) || gmt_M_is_dnan(Y) || gmt_M_is_dnan(S)) set_exit++;
330
331 if (set_exit == 1) {
332 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong - 2\n");
333 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: x %12.1f\n", x);
334 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: y %12.1f\n", y);
335 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: \n");
336 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: P %12.7f\n", P);
337 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phig %12.7f\n", GMT->current.proj.g_phig);
338 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: \n");
339 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: B %12.7f\n", B);
340 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: D %12.7f\n", D);
341 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: u %12.1f\n", u);
342 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: v %12.6e\n", v);
343 }
344 if (set_exit || GMT->current.proj.g_debug > 1) {
345 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: t %12.7f\n", t);
346 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: Kp %12.1f\n", Kp);
347 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: Kp2 %12.1f\n", Kp2);
348 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: X %12.1f\n", X);
349 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: Y %12.1f\n", Y);
350 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: S %12.7f\n", S);
351 }
352
353 if (h == 0) {
354 phi = atan(S/sqrt(one_m_e2*(1.0 - e2 - S*S)));
355 /* if (gmt_M_is_dnan(phi)) set_exit++; */
356 }
357 else {
358 double t1, t2;
359 niter = 0;
360 t2 = h*h/(R*R*one_m_e2);
361
362 sphi = S/(one_m_e2/sqrt(1.0 - e2*S*S) + h/R);
363 phi = asin(sphi);
364
365 t1 = (1.0/sqrt(1.0 - e2*sphi*sphi) + h/R);
366 E = t1 * t1 - e2*sphi*sphi*(1.0/(1.0 - e2*sphi*sphi) - t2);
367
368 if (gmt_M_is_dnan(E)) set_exit++;
369
370 if (set_exit == 1) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong - 3\n");
371 if (GMT->current.proj.g_debug > 1 || set_exit) {
372 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: asinS %12.7f\n", asind(S));
373 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phi %12.7f\n", R2D*phi);
374 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: E %12.7f\n", E);
375 }
376
377 do {
378 niter++;
379 phi_last = phi;
380 t = P*P*(1.0 - e2*cphig*cphig) - E*one_m_e2;
381 Kp2 = (1.0 - 4.0*(t/u)*(v/u));
382 if (Kp2 < 0.0)
383 Kp = -u/(2.0*t);
384 else
385 Kp = (-u + sqrt(u*u-4.0*t*v))/(2.0*t);
386 X = R*((B-H/Kp)*cphi1 - (y/Kp - D)*sphi1);
387 Y = R*x/Kp;
388 S = (y/Kp-D)*cphi1 + (B-H/Kp)*sphi1;
389 sphi = S/(one_m_e2/sqrt(1.0 - e2*sphi*sphi) + h/R);
390 phi = asin(sphi);
391 t1 = (1.0/sqrt(1.0 - e2*sphi*sphi) + h/R);
392 E = t1 * t1 - e2*sphi*sphi*(1.0/(1.0 - e2*sphi*sphi) - t2);
393
394 if (gmt_M_is_dnan(Kp) || gmt_M_is_dnan(X) || gmt_M_is_dnan(Y) || gmt_M_is_dnan(S) || gmt_M_is_dnan(phi) || gmt_M_is_dnan(E)) set_exit++;
395 if (set_exit == 1) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong - 4 \n");
396 if (set_exit || GMT->current.proj.g_debug > 1) {
397 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: iter %d\n", niter);
398 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: t %12.7f\n", t);
399 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: Kp %12.1f\n", Kp);
400 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: X %12.1f\n", X);
401 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: Y %12.1f\n", Y);
402 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: S %12.7f\n", S);
403 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phi %12.7f\n", phi*R2D);
404 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: E %12.7f\n", E);
405 }
406 }
407 while (fabs(phi - phi_last) > 1e-7);
408 }
409 if (set_exit == 1) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong - 5\n");
410 if (set_exit || GMT->current.proj.g_debug > 1) {
411 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: gmtproj_genper_tolatlong phi %12.7f\n", phi*R2D);
412 return GMT_PROJECTION_ERROR;
413 }
414 *lat = phi * R2D;
415 *lon = atan2d (Y, X) + GMT->current.proj.g_lon0;
416 return (GMT_OK);
417 }
418
gmtproj_genper_setup(struct GMT_CTRL * GMT,double h0,double altitude,double lat,double lon0)419 GMT_LOCAL void gmtproj_genper_setup (struct GMT_CTRL *GMT, double h0, double altitude, double lat, double lon0) {
420 /* if ellipsoid lat0 is geodetic latitude and must convert to geocentric latitude */
421
422 double N1, phig_last;
423 double R, H, P;
424 double sphi1, cphi1, sphig, cphig;
425 double a, e2, phig;
426
427 int niter;
428
429 a = GMT->current.proj.EQ_RAD;
430 e2 = GMT->current.proj.ECC2;
431
432 h0 *= 1e3;
433
434 sincosd (lat, &sphi1, &cphi1);
435 sphig = sphi1; cphig = cphi1;
436
437 N1 = a / sqrt (1.0 - (e2*sphi1*sphi1));
438 niter = 0;
439
440 if (GMT->current.proj.g_radius) { /* Altitude given as the radial distance from the center of the Earth (in km) */
441 H = fabs (altitude * METERS_IN_A_KM) - a;
442 P = H/a + 1.0;
443 phig = lat;
444 }
445 else if (GMT->current.proj.g_geosync) {/* Select implicit altitude of geosynchronous viewpoint */
446 double temp = 86164.1/TWO_PI; /* Siderial day rotation rate */
447 H = pow (3.98603e14*temp*temp, 0.3333) - a; /* Standard gravitational parameter GM for Earth */
448 P = H/a + 1.0;
449 phig = lat - asind (N1*e2*sphi1*cphi1/(P*a));
450 sincosd (phig, &sphig, &cphig);
451 if (cphi1 != 0.0)
452 H = P*a*cphig/cphi1 - N1 - h0;
453 else
454 H = P*a - N1 - h0;
455 }
456 else if (GMT->current.proj.g_earth_radius) { /* Altitude given (in Earth radii) */
457 P = altitude;
458 /* need to setup H from P equation */
459 phig = lat - asind (N1*e2*sphi1*cphi1/(P*a));
460 sincosd (phig, &sphig, &cphig);
461 if (cphi1 != 0.0)
462 H = P*a*cphig/cphi1 - N1 - h0;
463 else
464 H = P*a - N1 - h0;
465 }
466 else { /* Altitude given as normal (in km) */
467 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: altitude %f\n", altitude);
468 H = altitude*1e3;
469 /* need to setup P from iterating phig */
470 phig = lat;
471 do {
472 niter++;
473 sincosd (phig, &sphig, &cphig);
474 P = (cphi1/cphig) * (H + N1 + h0)/a;
475 phig_last = phig;
476 phig = lat - asind(N1*e2*sphi1*cphi1/(P*a));
477 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: %2d P %12.7f phig %12.7f\n", niter, P, phig);
478 }
479 while (fabs (phig - phig_last) > 1e-9);
480 sincosd (phig, &sphig, &cphig);
481 P = (cphi1/cphig)*(H + N1 + h0)/a;
482 }
483
484 /* R = H/(P-1.0); */
485 R = a;
486 /* XXX Which one is it ? */
487
488 GMT->current.proj.g_H = H;
489 GMT->current.proj.g_P = P;
490 GMT->current.proj.g_R = R;
491 GMT->current.proj.g_lon0 = lon0;
492 GMT->current.proj.g_sphi1 = sphi1;
493 GMT->current.proj.g_cphi1 = cphi1;
494
495 GMT->current.proj.g_phig = phig;
496 GMT->current.proj.g_sphig = sphig;
497 GMT->current.proj.g_cphig = cphig;
498
499 sincosd (lat-phig, &(GMT->current.proj.g_sdphi), &(GMT->current.proj.g_cdphi));
500
501 GMT->current.proj.g_L = 1.0 - e2*cphi1*cphi1;
502 GMT->current.proj.g_G = 1.0 - e2*sphi1*sphi1;
503 GMT->current.proj.g_J = 2.0*e2*sphi1*cphi1;
504 GMT->current.proj.g_B = P*GMT->current.proj.g_cdphi;
505 GMT->current.proj.g_D = P*GMT->current.proj.g_sdphi;
506 GMT->current.proj.g_BLH = -2.0*GMT->current.proj.g_B*GMT->current.proj.g_L*H;
507 GMT->current.proj.g_DG = 2.0*GMT->current.proj.g_D*GMT->current.proj.g_G;
508 GMT->current.proj.g_BJ = GMT->current.proj.g_B*GMT->current.proj.g_J;
509 GMT->current.proj.g_HJ = H*GMT->current.proj.g_J;
510 GMT->current.proj.g_DHJ = GMT->current.proj.g_D*GMT->current.proj.g_HJ;
511 GMT->current.proj.g_LH2 = GMT->current.proj.g_L*H*H;
512
513 if (GMT->current.proj.g_debug > 0) {
514 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: a %12.4f\n", a);
515 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: R %12.4f\n", R);
516 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: e^2 %12.7f\n", e2);
517 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: H %12.4f\n", H);
518 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phi1 %12.4f\n", lat);
519 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: lon0 %12.4f\n", lon0);
520 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: h0 %12.4f\n", h0);
521 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: N1 %12.1f\n", N1);
522 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: P %12.7f\n", P);
523 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: phig %12.7f\n", phig);
524 }
525
526 return;
527 }
528
529 /* LINEAR TRANSFORMATIONS */
530
gmtlib_translin(struct GMT_CTRL * GMT,double forw,double * inv)531 void gmtlib_translin (struct GMT_CTRL *GMT, double forw, double *inv) /* Linear forward */ {
532 gmt_M_unused(GMT);
533 *inv = forw;
534 }
535
gmtproj_translind(struct GMT_CTRL * GMT,double forw,double * inv)536 GMT_LOCAL void gmtproj_translind (struct GMT_CTRL *GMT, double forw, double *inv) /* Linear forward, but with degrees*/ {
537 gmt_M_wind_lon (GMT, forw) /* Make sure forw is in -180/+180 range after removing central meridian */
538 *inv = forw ;
539 }
540
gmtproj_itranslind(struct GMT_CTRL * GMT,double * forw,double inv)541 GMT_LOCAL void gmtproj_itranslind (struct GMT_CTRL *GMT, double *forw, double inv) /* Linear inverse, but with degrees*/ {
542 while (inv < -GMT_180) inv += 360.0;
543 while (inv > +GMT_180) inv -= 360.0;
544 *forw = inv + GMT->current.proj.central_meridian;
545 }
546
gmtlib_itranslin(struct GMT_CTRL * GMT,double * forw,double inv)547 void gmtlib_itranslin (struct GMT_CTRL *GMT, double *forw, double inv) /* Linear inverse */ {
548 gmt_M_unused(GMT);
549 *forw = inv;
550 }
551
gmtproj_translog10(struct GMT_CTRL * GMT,double forw,double * inv)552 GMT_LOCAL void gmtproj_translog10 (struct GMT_CTRL *GMT, double forw, double *inv) /* Log10 forward */ {
553 *inv = d_log10 (GMT, forw);
554 }
555
gmtproj_itranslog10(struct GMT_CTRL * GMT,double * forw,double inv)556 GMT_LOCAL void gmtproj_itranslog10 (struct GMT_CTRL *GMT, double *forw, double inv) /* Log10 inverse */ {
557 gmt_M_unused(GMT);
558 *forw = pow (10.0, inv);
559 }
560
gmtproj_transpowx(struct GMT_CTRL * GMT,double x,double * x_in)561 GMT_LOCAL void gmtproj_transpowx (struct GMT_CTRL *GMT, double x, double *x_in) /* pow x forward */ {
562 *x_in = pow (x, GMT->current.proj.xyz_pow[GMT_X]);
563 }
564
gmtproj_itranspowx(struct GMT_CTRL * GMT,double * x,double x_in)565 GMT_LOCAL void gmtproj_itranspowx (struct GMT_CTRL *GMT, double *x, double x_in) /* pow x inverse */ {
566 *x = pow (x_in, GMT->current.proj.xyz_ipow[GMT_X]);
567 }
568
gmtproj_transpowy(struct GMT_CTRL * GMT,double y,double * y_in)569 GMT_LOCAL void gmtproj_transpowy (struct GMT_CTRL *GMT, double y, double *y_in) /* pow y forward */ {
570 *y_in = pow (y, GMT->current.proj.xyz_pow[GMT_Y]);
571 }
572
gmtproj_itranspowy(struct GMT_CTRL * GMT,double * y,double y_in)573 GMT_LOCAL void gmtproj_itranspowy (struct GMT_CTRL *GMT, double *y, double y_in) /* pow y inverse */ {
574 *y = pow (y_in, GMT->current.proj.xyz_ipow[GMT_Y]);
575 }
576
gmtproj_transpowz(struct GMT_CTRL * GMT,double z,double * z_in)577 GMT_LOCAL void gmtproj_transpowz (struct GMT_CTRL *GMT, double z, double *z_in) /* pow z forward */ {
578 *z_in = pow (z, GMT->current.proj.xyz_pow[GMT_Z]);
579 }
580
gmtproj_itranspowz(struct GMT_CTRL * GMT,double * z,double z_in)581 GMT_LOCAL void gmtproj_itranspowz (struct GMT_CTRL *GMT, double *z, double z_in) /* pow z inverse */ {
582 *z = pow (z_in, GMT->current.proj.xyz_ipow[GMT_Z]);
583 }
584
585 /* -JP POLAR (r-theta) PROJECTION */
586
gmtproj_planet_radius(struct GMT_CTRL * GMT,char * modifier)587 GMT_LOCAL double gmtproj_planet_radius (struct GMT_CTRL *GMT, char *modifier) {
588 static char *U[2] = {"m", "km"};
589 unsigned int k = 0;
590 double r;
591 /* Set planetary radius in correct units (m or km) depending on y-range */
592 r = GMT->current.setting.ref_ellipsoid[GMT->current.setting.proj_ellipsoid].eq_radius; /* In meters */
593 if ((r/ (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO])) >= METERS_IN_A_KM) { /* -R seems given in km */
594 r /= METERS_IN_A_KM;
595 k = 1;
596 }
597 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Planetary radius (%s) automatically set to %g %s\n", modifier, r, U[k]);
598 return (r);
599 }
600
gmtproj_vpolar(struct GMT_CTRL * GMT,double lon0)601 GMT_LOCAL void gmtproj_vpolar (struct GMT_CTRL *GMT, double lon0) {
602 /* Set up a Polar (theta,r) transformation */
603
604 GMT->current.proj.p_base_angle = lon0;
605 GMT->current.proj.central_meridian = 0.5 * (GMT->common.R.wesn[XHI] + GMT->common.R.wesn[XLO]);
606
607 if (GMT->current.proj.flip) { /* Want radial direction inwards */
608 if (GMT->current.proj.flip_radius < 0.0) /* Flag to just flip z = north - r */
609 GMT->current.proj.flip_radius = GMT->common.R.wesn[YHI];
610 else if (GMT->current.proj.flip_radius == 0.0) /* Flag to just flip z = planet_radius - r */
611 GMT->current.proj.flip_radius = gmtproj_planet_radius (GMT, "+fp");
612 /* else the radius was set specifically */
613 }
614 if (GMT->current.proj.z_down) { /* Annotate a flavor of z = radius - r */
615 if (GMT->current.proj.z_down == GMT_ZDOWN_ZP) { /* Given z; annotate r = planet_radius - z */
616 if (GMT->current.proj.flip_radius > 0.0) /* Already obtained above */
617 GMT->current.proj.z_radius = GMT->current.proj.flip_radius;
618 else
619 GMT->current.proj.z_radius = GMT->current.proj.flip_radius = gmtproj_planet_radius (GMT, "+zp");
620 GMT->current.proj.flip = true;
621 }
622 else if (GMT->current.proj.z_down == GMT_ZDOWN_Z) /* z = north - r */
623 GMT->current.proj.z_radius = GMT->common.R.wesn[YHI];
624 }
625
626 /* Plus pretend that it is kind of a geographic polar projection */
627 GMT->current.proj.north_pole = GMT->current.proj.got_elevations;
628 GMT->current.proj.pole = (GMT->current.proj.got_elevations) ? 90.0 : 0.0;
629 GMT->current.proj.radial_offset /= GMT->current.proj.pars[0]; /* Convert any radial offset in inches to user units so we can use it in gmtproj_polar/ipolar */
630 }
631
gmtproj_polar(struct GMT_CTRL * GMT,double x,double y,double * x_i,double * y_i)632 GMT_LOCAL void gmtproj_polar (struct GMT_CTRL *GMT, double x, double y, double *x_i, double *y_i) {
633 /* Transform x and y to polar(cylindrical) coordinates */
634 if (GMT->current.proj.got_azimuths) x = 90.0 - x; /* Azimuths, not directions given as x */
635 if (GMT->current.proj.flip) y = GMT->current.proj.flip_radius - y; /* Depth down or elevations given as y */
636 sincosd (x - GMT->current.proj.p_base_angle, y_i, x_i); /* Change base line angle */
637 (*x_i) *= (y + GMT->current.proj.radial_offset); /* Allow for inner circle radius before we start plotting */
638 (*y_i) *= (y + GMT->current.proj.radial_offset);
639 }
640
gmtproj_ipolar(struct GMT_CTRL * GMT,double * x,double * y,double x_i,double y_i)641 GMT_LOCAL void gmtproj_ipolar (struct GMT_CTRL *GMT, double *x, double *y, double x_i, double y_i) {
642 /* Inversely transform both x and y from polar(cylindrical) coordinates */
643 *x = d_atan2d (y_i, x_i) + GMT->current.proj.p_base_angle;
644 if (GMT->current.proj.got_azimuths) *x = 90.0 - (*x); /* Azimuths, not directions for x */
645 *y = hypot (x_i, y_i) - GMT->current.proj.radial_offset; /* Allow for inner circle radius */
646 if (GMT->current.proj.flip) *y = GMT->current.proj.flip_radius - (*y); /* Depth down or elevations for y */
647 }
648
649 /* -JM MERCATOR PROJECTION */
650
gmtproj_vmerc(struct GMT_CTRL * GMT,double lon0,double lat0)651 GMT_LOCAL void gmtproj_vmerc (struct GMT_CTRL *GMT, double lon0, double lat0) {
652 /* Set up a Mercator transformation with origin at (lon0, lat0) */
653
654 double aux_lat0 = (GMT->current.proj.GMT_convert_latitudes) ? gmt_M_latg_to_latc (GMT, lat0) : lat0;
655
656 GMT->current.proj.central_meridian = lon0;
657 /* Need geodetic latitude in this expression: */
658 GMT->current.proj.j_x = cosd (lat0) / d_sqrt (1.0 - GMT->current.proj.ECC2 * sind (lat0) * sind (lat0)) * GMT->current.proj.EQ_RAD;
659 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
660 /* Need conformal latitude in this expression (same as in gmtproj_merc_sph) */
661 GMT->current.proj.j_yc = (fabs (lat0) > 0.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * aux_lat0)) : 0.0;
662 }
663
664 /* Mercator projection for the sphere */
665
gmtproj_merc_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)666 GMT_LOCAL void gmtproj_merc_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
667 /* Convert lon/lat to Mercator x/y (GMT->current.proj.EQ_RAD in GMT->current.proj.j_x) */
668
669 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
670 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
671
672 *x = GMT->current.proj.j_x * D2R * lon;
673 *y = (fabs (lat) < 90.0) ? GMT->current.proj.j_x * d_log (GMT, tand (45.0 + 0.5 * lat)) : copysign (DBL_MAX, lat);
674 }
675
gmtproj_imerc_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)676 GMT_LOCAL void gmtproj_imerc_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
677 /* Convert Mercator x/y to lon/lat (GMT->current.proj.EQ_RAD in GMT->current.proj.j_ix) */
678
679 *lon = x * GMT->current.proj.j_ix * R2D + GMT->current.proj.central_meridian;
680 *lat = atand (sinh (y * GMT->current.proj.j_ix));
681 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
682 }
683
684 /* -JY CYLINDRICAL EQUAL-AREA PROJECTION */
685
gmtproj_vcyleq(struct GMT_CTRL * GMT,double lon0,double slat)686 GMT_LOCAL void gmtproj_vcyleq (struct GMT_CTRL *GMT, double lon0, double slat) {
687 /* Set up a Cylindrical equal-area transformation */
688
689 gmtproj_check_R_J (GMT, &lon0);
690 GMT->current.proj.central_meridian = lon0;
691 GMT->current.proj.j_x = GMT->current.proj.EQ_RAD * D2R * cosd (slat);
692 GMT->current.proj.j_y = GMT->current.proj.EQ_RAD / cosd (slat);
693 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
694 GMT->current.proj.j_iy = 1.0 / GMT->current.proj.j_y;
695 }
696
gmtproj_cyleq(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)697 GMT_LOCAL void gmtproj_cyleq (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
698 /* Convert lon/lat to Cylindrical equal-area x/y */
699
700 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
701 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
702
703 *x = lon * GMT->current.proj.j_x;
704 *y = GMT->current.proj.j_y * sind (lat);
705 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
706 (*x) *= GMT->current.proj.Dx;
707 (*y) *= GMT->current.proj.Dy;
708 }
709 }
710
gmtproj_icyleq(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)711 GMT_LOCAL void gmtproj_icyleq (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
712 /* Convert Cylindrical equal-area x/y to lon/lat */
713
714 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
715 x *= GMT->current.proj.iDx;
716 y *= GMT->current.proj.iDy;
717 }
718 *lon = x * GMT->current.proj.j_ix + GMT->current.proj.central_meridian;
719 *lat = d_asind (y * GMT->current.proj.j_iy);
720 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
721 }
722
723 /* -JQ CYLINDRICAL EQUIDISTANT PROJECTION */
724
gmtproj_vcyleqdist(struct GMT_CTRL * GMT,double lon0,double slat)725 GMT_LOCAL void gmtproj_vcyleqdist (struct GMT_CTRL *GMT, double lon0, double slat) {
726 /* Set up a Cylindrical equidistant transformation */
727
728 gmtproj_check_R_J (GMT, &lon0);
729 GMT->current.proj.central_meridian = lon0;
730 GMT->current.proj.j_x = D2R * GMT->current.proj.EQ_RAD * cosd (slat);
731 GMT->current.proj.j_y = D2R * GMT->current.proj.EQ_RAD;
732 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
733 GMT->current.proj.j_iy = 1.0 / GMT->current.proj.j_y;
734 }
735
gmtproj_cyleqdist(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)736 GMT_LOCAL void gmtproj_cyleqdist (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
737 /* Convert lon/lat to Cylindrical equidistant x/y */
738
739 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
740 *x = lon * GMT->current.proj.j_x;
741 *y = lat * GMT->current.proj.j_y;
742 }
743
gmtproj_icyleqdist(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)744 GMT_LOCAL void gmtproj_icyleqdist (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
745 /* Convert Cylindrical equal-area x/y to lon/lat */
746
747 *lon = x * GMT->current.proj.j_ix + GMT->current.proj.central_meridian;
748 *lat = y * GMT->current.proj.j_iy;
749 }
750
751 /* -JJ MILLER CYLINDRICAL PROJECTION */
752
753 //#define CHRISTMAS
754 /* Turning on Christmas makes the Miller projection a triangular projection
755 * that projects 90 degrees of longitude and latitudes 45-90 into a 45-degree
756 * triangle. Doing all for quadrants results in a square map with radial
757 * meridians and lots of distortion along the boundaries. This was used to
758 * build a 3-D cube of the world with this triangle projection being used to
759 * map the top (N polar to 34N) and bottom (S pole to 45S) sides, with the
760 * remaining 4 sides just being -JQ maps. I left it here since I may want
761 * to mess with this in the future. P. Wessel, Dec. 2016.
762 */
763 #ifdef CHRISTMAS
764 /* Bypass Miller projection entirely and introduce a triangle projection */
gmtproj_vmiller(struct GMT_CTRL * GMT,double lon0,double slat)765 GMT_LOCAL void gmtproj_vmiller (struct GMT_CTRL *GMT, double lon0, double slat) {
766 /* Set up a Cylindrical equidistant transformation */
767 GMT->current.proj.north_pole = (slat > 0.0);
768 gmtproj_check_R_J (GMT, &lon0);
769 GMT->current.proj.central_meridian = lon0;
770 GMT->current.proj.j_x = 0.25 * D2R * GMT->current.proj.EQ_RAD;
771 GMT->current.proj.j_y = 0.25 * D2R * GMT->current.proj.EQ_RAD;
772 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
773 GMT->current.proj.j_iy = 1.0 / GMT->current.proj.j_y;
774 }
775
gmtproj_miller(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)776 GMT_LOCAL void gmtproj_miller (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
777 /* Convert lon/lat to Cylindrical equidistant x/y */
778
779 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
780 if (lat > 0.0) {
781 *x = (0.5 + lon * (90.0 - lat) / 4050.0) * GMT->current.proj.j_x;
782 *y = ((lat-45.0) / 90.0) * GMT->current.proj.j_y;
783 }
784 else {
785 *x = (0.5 - lon * (90.0 + lat) / 4050.0) * GMT->current.proj.j_x;
786 *y = -((lat+45) / 90.0) * GMT->current.proj.j_y;
787 }
788 }
789
gmtproj_imiller(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)790 GMT_LOCAL void gmtproj_imiller (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
791 /* Convert Cylindrical equal-area x/y to lon/lat */
792
793 if (GMT->current.proj.north_pole) {
794 *lat = 45.0 + 90.0 * y * GMT->current.proj.j_iy;
795 *lon = (4050.0 * (x * GMT->current.proj.j_ix - 0.5)) / (90.0 - *lat) + GMT->current.proj.central_meridian;
796 }
797 else {
798 *lat = -(45.0 + 90.0 * y * GMT->current.proj.j_iy);
799 *lon = -(4050.0 * (x * GMT->current.proj.j_ix - 0.5)) / (90.0 + *lat) + GMT->current.proj.central_meridian;
800 }
801 }
802 #else
gmtproj_vmiller(struct GMT_CTRL * GMT,double lon0)803 GMT_LOCAL void gmtproj_vmiller (struct GMT_CTRL *GMT, double lon0) {
804 /* Set up a Miller Cylindrical transformation */
805
806 gmtproj_check_R_J (GMT, &lon0);
807 GMT->current.proj.central_meridian = lon0;
808 GMT->current.proj.j_x = D2R * GMT->current.proj.EQ_RAD;
809 GMT->current.proj.j_y = 1.25 * GMT->current.proj.EQ_RAD;
810 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
811 GMT->current.proj.j_iy = 1.0 / GMT->current.proj.j_y;
812 }
813
gmtproj_miller(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)814 GMT_LOCAL void gmtproj_miller (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
815 /* Convert lon/lat to Miller Cylindrical x/y */
816
817 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
818 *x = lon * GMT->current.proj.j_x;
819 *y = GMT->current.proj.j_y * d_log (GMT, tand (45.0 + 0.4 * lat));
820 }
821
gmtproj_imiller(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)822 GMT_LOCAL void gmtproj_imiller (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
823 /* Convert Miller Cylindrical x/y to lon/lat */
824
825 *lon = x * GMT->current.proj.j_ix + GMT->current.proj.central_meridian;
826 *lat = 2.5 * atand (exp (y * GMT->current.proj.j_iy)) - 112.5;
827 }
828 #endif
829
830 /* -JCyl_stere CYLINDRICAL STEREOGRAPHIC PROJECTION */
831
gmtproj_vcylstereo(struct GMT_CTRL * GMT,double lon0,double slat)832 GMT_LOCAL void gmtproj_vcylstereo (struct GMT_CTRL *GMT, double lon0, double slat) {
833 /* Set up a Cylindrical Stereographic transformation */
834
835 gmtproj_check_R_J (GMT, &lon0);
836 GMT->current.proj.central_meridian = lon0;
837 GMT->current.proj.j_x = GMT->current.proj.EQ_RAD * D2R * cosd (slat);
838 GMT->current.proj.j_y = GMT->current.proj.EQ_RAD * (1.0 + cosd (slat));
839 GMT->current.proj.j_ix = 1.0 / GMT->current.proj.j_x;
840 GMT->current.proj.j_iy = 1.0 / GMT->current.proj.j_y;
841 }
842
gmtproj_cylstereo(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)843 GMT_LOCAL void gmtproj_cylstereo (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
844 /* Convert lon/lat to Cylindrical Stereographic x/y */
845
846 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
847 *x = lon * GMT->current.proj.j_x;
848 *y = GMT->current.proj.j_y * tand (0.5 * lat);
849 }
850
gmtproj_icylstereo(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)851 GMT_LOCAL void gmtproj_icylstereo (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
852 /* Convert Cylindrical Stereographic x/y to lon/lat */
853
854 *lon = x * GMT->current.proj.j_ix + GMT->current.proj.central_meridian;
855 *lat = 2.0 * atand (y * GMT->current.proj.j_iy);
856 }
857
858 /* -JS POLAR STEREOGRAPHIC PROJECTION */
859
gmtproj_vstereo(struct GMT_CTRL * GMT,double lon0,double lat0,double horizon)860 GMT_LOCAL void gmtproj_vstereo (struct GMT_CTRL *GMT, double lon0, double lat0, double horizon) {
861 /* Set up a Stereographic transformation */
862
863 double clat;
864 if (GMT->current.proj.GMT_convert_latitudes) { /* Set Conformal radius and pole latitude */
865 gmtlib_scale_eqrad (GMT);
866 clat = gmt_M_latg_to_latc (GMT, lat0);
867 }
868 else
869 clat = lat0;
870
871 GMT->current.proj.central_meridian = lon0;
872 GMT->current.proj.pole = lat0; /* This is always geodetic */
873 sincosd (clat, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp)); /* These may be conformal */
874 GMT->current.proj.north_pole = (lat0 > 0.0);
875 GMT->current.proj.s_c = 2.0 * GMT->current.proj.EQ_RAD * GMT->current.setting.proj_scale_factor;
876 GMT->current.proj.s_ic = 1.0 / GMT->current.proj.s_c;
877 GMT->current.proj.f_horizon = horizon;
878 GMT->current.proj.rho_max = tand (0.5 * horizon) * GMT->current.proj.s_c;
879 }
880
gmtproj_plrs_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)881 GMT_LOCAL void gmtproj_plrs_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
882 /* Convert lon/lat to x/y using Spherical polar projection */
883 double rho, slon, clon;
884
885 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
886 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
887
888 sincosd (lon, &slon, &clon);
889 if (GMT->current.proj.north_pole) {
890 rho = GMT->current.proj.s_c * tand (45.0 - 0.5 * lat);
891 *y = -rho * clon;
892 *x = rho * slon;
893 }
894 else {
895 rho = GMT->current.proj.s_c * tand (45.0 + 0.5 * lat);
896 *y = rho * clon;
897 *x = rho * slon;
898 }
899 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
900 (*x) *= GMT->current.proj.Dx;
901 (*y) *= GMT->current.proj.Dy;
902 }
903 }
904
gmtproj_iplrs_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)905 GMT_LOCAL void gmtproj_iplrs_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
906 /* Convert Spherical polar x/y to lon/lat */
907 double c;
908
909 if (x == 0.0 && y == 0.0) {
910 *lon = GMT->current.proj.central_meridian;
911 *lat = GMT->current.proj.pole;
912 return;
913 }
914
915 if (GMT->current.proj.GMT_convert_latitudes) { /* Undo effect of fudge factors */
916 x *= GMT->current.proj.iDx;
917 y *= GMT->current.proj.iDy;
918 }
919 c = 2.0 * atan (hypot (x, y) * GMT->current.proj.s_ic);
920
921 if (GMT->current.proj.north_pole) {
922 *lon = GMT->current.proj.central_meridian + d_atan2d (x, -y);
923 *lat = d_asind (cos (c));
924 }
925 else {
926 *lon = GMT->current.proj.central_meridian + d_atan2d (x, y);
927 *lat = d_asind (-cos (c));
928 }
929 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
930
931 }
932
gmtproj_stereo1_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)933 GMT_LOCAL void gmtproj_stereo1_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
934 /* Convert lon/lat to x/y using Spherical stereographic projection, oblique view */
935
936 double sin_dlon, cos_dlon, s, c, cc, A;
937
938 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
939 sincosd (lon - GMT->current.proj.central_meridian, &sin_dlon, &cos_dlon);
940 sincosd (lat, &s, &c);
941 cc = c * cos_dlon;
942 A = GMT->current.proj.s_c / (1.0 + GMT->current.proj.sinp * s + GMT->current.proj.cosp * cc);
943 *x = A * c * sin_dlon;
944 *y = A * (GMT->current.proj.cosp * s - GMT->current.proj.sinp * cc);
945 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
946 (*x) *= GMT->current.proj.Dx;
947 (*y) *= GMT->current.proj.Dy;
948 }
949 }
950
gmtproj_istereo_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)951 GMT_LOCAL void gmtproj_istereo_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
952 double rho, c, sin_c, cos_c;
953
954 if (GMT->current.proj.GMT_convert_latitudes) { /* Undo effect of fudge factors */
955 x *= GMT->current.proj.iDx;
956 y *= GMT->current.proj.iDy;
957 }
958 rho = hypot (x, y);
959 if (rho > GMT->current.proj.rho_max) {
960 *lat = *lon = GMT->session.d_NaN;
961 return;
962 }
963 c = 2.0 * atan (rho * GMT->current.proj.s_ic);
964 sincos (c, &sin_c, &cos_c);
965 if (rho != 0.0) sin_c /= rho;
966 *lat = asind (cos_c * GMT->current.proj.sinp + y * sin_c * GMT->current.proj.cosp);
967 *lon = d_atan2d (x * sin_c, cos_c * GMT->current.proj.cosp - y * sin_c * GMT->current.proj.sinp) + GMT->current.proj.central_meridian;
968 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
969 }
970
971 /* Spherical equatorial view */
972
gmtproj_stereo2_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)973 GMT_LOCAL void gmtproj_stereo2_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
974 /* Convert lon/lat to x/y using stereographic projection, equatorial view */
975
976 double dlon, s, c, clon, slon, A;
977
978 dlon = lon - GMT->current.proj.central_meridian;
979 if (doubleAlmostEqual (dlon, 180.0)) {
980 *x = *y = 0.0;
981 }
982 else {
983 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
984 sincosd (lat, &s, &c);
985 sincosd (dlon, &slon, &clon);
986 A = GMT->current.proj.s_c / (1.0 + c * clon);
987 *x = A * c * slon;
988 *y = A * s;
989 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
990 (*x) *= GMT->current.proj.Dx;
991 (*y) *= GMT->current.proj.Dy;
992 }
993 }
994 }
995
996 /* -JL LAMBERT CONFORMAL CONIC PROJECTION */
997
gmtproj_vlamb(struct GMT_CTRL * GMT,double rlong0,double rlat0,double pha,double phb)998 GMT_LOCAL void gmtproj_vlamb (struct GMT_CTRL *GMT, double rlong0, double rlat0, double pha, double phb) {
999 /* Set up a Lambert Conformal Conic projection (Spherical when e = 0) */
1000
1001 double sin_pha, cos_pha, sin_phb, cos_phb, t_pha, m_pha, t_phb, m_phb, t_rlat0;
1002
1003 gmtproj_check_R_J (GMT, &rlong0);
1004 GMT->current.proj.north_pole = (GMT->common.R.wesn[YHI] > 0.0 && (GMT->common.R.wesn[YLO] >= 0.0 || (-GMT->common.R.wesn[YLO]) < GMT->common.R.wesn[YHI]));
1005 GMT->current.proj.pole = (GMT->current.proj.north_pole) ? 90.0 : -90.0;
1006 sincosd (pha, &sin_pha, &cos_pha);
1007 sincosd (phb, &sin_phb, &cos_phb);
1008
1009 t_pha = tand (45.0 - 0.5 * pha) / pow ((1.0 - GMT->current.proj.ECC *
1010 sin_pha) / (1.0 + GMT->current.proj.ECC * sin_pha), GMT->current.proj.half_ECC);
1011 m_pha = cos_pha / d_sqrt (1.0 - GMT->current.proj.ECC2 * sin_pha * sin_pha);
1012 t_phb = tand (45.0 - 0.5 * phb) / pow ((1.0 - GMT->current.proj.ECC *
1013 sin_phb) / (1.0 + GMT->current.proj.ECC * sin_phb), GMT->current.proj.half_ECC);
1014 m_phb = cos_phb / d_sqrt (1.0 - GMT->current.proj.ECC2 * sin_phb * sin_phb);
1015 t_rlat0 = tand (45.0 - 0.5 * rlat0) /
1016 pow ((1.0 - GMT->current.proj.ECC * sind (rlat0)) /
1017 (1.0 + GMT->current.proj.ECC * sind (rlat0)), GMT->current.proj.half_ECC);
1018
1019 if (doubleAlmostEqualZero (pha, phb))
1020 GMT->current.proj.l_N = sind (pha);
1021 else
1022 GMT->current.proj.l_N = (d_log (GMT, m_pha) - d_log (GMT, m_phb))/(d_log (GMT, t_pha) - d_log (GMT, t_phb));
1023
1024 GMT->current.proj.l_i_N = 1.0 / GMT->current.proj.l_N;
1025 GMT->current.proj.l_F = m_pha / (GMT->current.proj.l_N * pow (t_pha, GMT->current.proj.l_N));
1026 GMT->current.proj.central_meridian = rlong0;
1027 GMT->current.proj.l_rF = GMT->current.proj.EQ_RAD * GMT->current.proj.l_F;
1028 GMT->current.proj.l_i_rF = 1.0 / GMT->current.proj.l_rF;
1029 GMT->current.proj.l_rho0 = GMT->current.proj.l_rF * pow (t_rlat0, GMT->current.proj.l_N);
1030 GMT->current.proj.l_Nr = GMT->current.proj.l_N * D2R;
1031 GMT->current.proj.l_i_Nr = 1.0 / GMT->current.proj.l_Nr;
1032 }
1033
gmtproj_lamb(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1034 GMT_LOCAL void gmtproj_lamb (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1035 double rho, theta, hold1, hold2, hold3, es, s, c;
1036
1037 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1038
1039 es = GMT->current.proj.ECC * sind (lat);
1040 hold2 = pow ((1.0 - es) / (1.0 + es), GMT->current.proj.half_ECC);
1041 hold3 = tand (45.0 - 0.5 * lat);
1042 hold1 = pow (hold3 / hold2, GMT->current.proj.l_N);
1043 rho = GMT->current.proj.l_rF * hold1;
1044 theta = GMT->current.proj.l_Nr * lon;
1045 sincos (theta, &s, &c);
1046 *x = rho * s;
1047 *y = GMT->current.proj.l_rho0 - rho * c;
1048 }
1049
gmtproj_ilamb(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1050 GMT_LOCAL void gmtproj_ilamb (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1051 int i;
1052 double theta, rho, t, tphi, phi, dy, r;
1053
1054 dy = GMT->current.proj.l_rho0 - y;
1055 theta = (GMT->current.proj.l_N < 0.0) ? d_atan2 (-x, -dy) : d_atan2 (x, dy);
1056 *lon = theta * GMT->current.proj.l_i_Nr + GMT->current.proj.central_meridian;
1057
1058 rho = copysign (hypot (x, dy), GMT->current.proj.l_N);
1059 t = pow (rho * GMT->current.proj.l_i_rF, GMT->current.proj.l_i_N);
1060 phi = 0.0; tphi = 999.0; /* Initialize phi = 0 */
1061 for (i = 0; i < GMT_PROJ_MAX_ITERATIONS && fabs (tphi - phi) > GMT_PROJ_CONV_LIMIT; i++) {
1062 tphi = phi;
1063 r = GMT->current.proj.ECC * sin (phi);
1064 phi = M_PI_2 - 2.0 * atan (t * pow ((1.0 - r) / (1.0 + r), GMT->current.proj.half_ECC));
1065 }
1066 *lat = phi * R2D;
1067 }
1068
1069 /* Spherical cases of Lambert */
1070
gmtproj_lamb_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1071 GMT_LOCAL void gmtproj_lamb_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1072 double rho, theta, t, s, c;
1073
1074 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1075 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
1076
1077 t = MAX (0.0, tand (45.0 - 0.5 * lat)); /* Guard against negative t */
1078 rho = GMT->current.proj.l_rF * pow (t, GMT->current.proj.l_N);
1079 theta = GMT->current.proj.l_Nr * lon;
1080
1081 sincos (theta, &s, &c);
1082 *x = rho * s;
1083 *y = GMT->current.proj.l_rho0 - rho * c;
1084 }
1085
gmtproj_ilamb_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1086 GMT_LOCAL void gmtproj_ilamb_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1087 double theta, rho, t, dy;
1088
1089 dy = GMT->current.proj.l_rho0 - y;
1090 theta = (GMT->current.proj.l_N < 0.0) ? d_atan2 (-x, -dy) : d_atan2 (x, dy);
1091 *lon = theta * GMT->current.proj.l_i_Nr + GMT->current.proj.central_meridian;
1092
1093 rho = copysign (hypot (x, dy), GMT->current.proj.l_N);
1094 t = pow (rho * GMT->current.proj.l_i_rF, GMT->current.proj.l_i_N);
1095 *lat = 90.0 - 2.0 * atand (t);
1096 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
1097 }
1098
1099 /* -JO OBLIQUE MERCATOR PROJECTION */
1100
gmtproj_obl(struct GMT_CTRL * GMT,double lon,double lat,double * olon,double * olat)1101 GMT_LOCAL void gmtproj_obl (struct GMT_CTRL *GMT, double lon, double lat, double *olon, double *olat) {
1102 /* Convert a longitude/latitude point to Oblique lon/lat (all in rads) */
1103 double p_cross_x[3], X[3];
1104
1105 gmt_geo_to_cart (GMT, lat, lon, X, false);
1106
1107 *olat = d_asin (gmt_dot3v (GMT, X, GMT->current.proj.o_FP));
1108
1109 gmt_cross3v (GMT, GMT->current.proj.o_FP, X, p_cross_x);
1110 gmt_normalize3v (GMT, p_cross_x);
1111
1112 *olon = copysign (d_acos (gmt_dot3v (GMT, p_cross_x, GMT->current.proj.o_FC)), gmt_dot3v (GMT, X, GMT->current.proj.o_FC));
1113 }
1114
gmtlib_iobl(struct GMT_CTRL * GMT,double * lon,double * lat,double olon,double olat)1115 void gmtlib_iobl (struct GMT_CTRL *GMT, double *lon, double *lat, double olon, double olat) {
1116 /* Convert a longitude/latitude point from Oblique lon/lat (all in rads) */
1117 double p_cross_x[3], X[3];
1118
1119 gmt_geo_to_cart (GMT, olat, olon, X, false);
1120 *lat = d_asin (gmt_dot3v (GMT, X, GMT->current.proj.o_IP));
1121
1122 gmt_cross3v (GMT, GMT->current.proj.o_IP, X, p_cross_x);
1123 gmt_normalize3v (GMT, p_cross_x);
1124
1125 *lon = copysign (d_acos (gmt_dot3v (GMT, p_cross_x, GMT->current.proj.o_IC)), gmt_dot3v (GMT, X, GMT->current.proj.o_IC));
1126
1127 while ((*lon) < 0.0) (*lon) += TWO_PI;
1128 while ((*lon) >= TWO_PI) (*lon) -= TWO_PI;
1129 }
1130
gmtproj_oblmrc(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1131 GMT_LOCAL void gmtproj_oblmrc (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1132 /* Convert a longitude/latitude point to Oblique Mercator coordinates
1133 * by way of rotation coordinates and then using regular Mercator */
1134 double tlon, tlat;
1135 /* o_shift deals with difference between user's origin and our logical origin */
1136
1137 gmtproj_obl (GMT, lon * D2R, lat * D2R, &tlon, &tlat);
1138
1139 *x = GMT->current.proj.j_x * tlon;
1140 *y = (fabs (tlat) < M_PI_2) ? GMT->current.proj.j_x * d_log (GMT, tan (M_PI_4 + 0.5 * tlat)) - GMT->current.proj.o_shift : copysign (DBL_MAX, tlat);
1141 if (GMT->current.proj.obl_flip) {
1142 /* Let oblique Equator be y-axis, so flip x and y but must let y be negative [that change takes place in map_setxy] */
1143 gmt_M_double_swap (*x, *y);
1144 if (GMT->current.proj.pars[1] < 0.0) *x = -*x, *y = -*y; /* S hemisphere must rotate 180 */
1145 }
1146 else if (GMT->current.proj.o_spole) {
1147 *x = -(*x);
1148 *y = -(*y);
1149 }
1150 }
1151
gmtproj_ioblmrc(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1152 GMT_LOCAL void gmtproj_ioblmrc (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1153 /* Convert a longitude/latitude point from Oblique Mercator coordinates
1154 * by way of regular Mercator and then rotate coordinates */
1155 double tlon, tlat;
1156 /* o_shift deals with difference between user's origin and our logical origin */
1157
1158 if (GMT->current.proj.obl_flip) {
1159 /* Had oblique Equator be y-axis */
1160 if (GMT->current.proj.pars[1] < 0.0) x = -x, y = -y; /* S hemisphere must rotate 180 */
1161 gmt_M_double_swap (x, y);
1162 }
1163 else if (GMT->current.proj.o_spole) {
1164 x = -x;
1165 y = -y;
1166 }
1167 tlon = x * GMT->current.proj.j_ix;
1168 y += GMT->current.proj.o_shift;
1169 tlat = atan (sinh (y * GMT->current.proj.j_ix));
1170
1171 gmtlib_iobl (GMT, lon, lat, tlon, tlat);
1172
1173 (*lon) *= R2D; (*lat) *= R2D;
1174 }
1175
1176 /* -JT TRANSVERSE MERCATOR PROJECTION */
1177
gmtproj_vtm(struct GMT_CTRL * GMT,double lon0,double lat0)1178 GMT_LOCAL void gmtproj_vtm (struct GMT_CTRL *GMT, double lon0, double lat0) {
1179 /* Set up an TM projection */
1180 double e1, s2, c2;
1181
1182 /* gmtproj_check_R_J (&lon0); */
1183 e1 = (1.0 - d_sqrt (GMT->current.proj.one_m_ECC2)) / (1.0 + d_sqrt (GMT->current.proj.one_m_ECC2));
1184 GMT->current.proj.t_e2 = GMT->current.proj.ECC2 * GMT->current.proj.i_one_m_ECC2;
1185 GMT->current.proj.t_c1 = 1.0 - (1.0/4.0) * GMT->current.proj.ECC2 - (3.0/64.0) * GMT->current.proj.ECC4 - (5.0/256.0) * GMT->current.proj.ECC6;
1186 GMT->current.proj.t_c2 = -((3.0/8.0) * GMT->current.proj.ECC2 + (3.0/32.0) * GMT->current.proj.ECC4 + (25.0/768.0) * GMT->current.proj.ECC6);
1187 GMT->current.proj.t_c3 = (15.0/128.0) * GMT->current.proj.ECC4 + (45.0/512.0) * GMT->current.proj.ECC6;
1188 GMT->current.proj.t_c4 = -(35.0/768.0) * GMT->current.proj.ECC6;
1189 GMT->current.proj.t_i1 = 1.0 / (GMT->current.proj.EQ_RAD * GMT->current.proj.t_c1);
1190 GMT->current.proj.t_i2 = (3.0/2.0) * e1 - (29.0/12.0) * pow (e1, 3.0);
1191 GMT->current.proj.t_i3 = (21.0/8.0) * e1 * e1 - (1537.0/128.0) * pow (e1, 4.0);
1192 GMT->current.proj.t_i4 = (151.0/24.0) * pow (e1, 3.0);
1193 GMT->current.proj.t_i5 = (1097.0/64.0) * pow (e1, 4.0);
1194 GMT->current.proj.central_meridian = lon0;
1195 GMT->current.proj.t_lat0 = lat0 * D2R; /* In radians */
1196 sincos (2.0 * GMT->current.proj.t_lat0, &s2, &c2);
1197 GMT->current.proj.t_M0 = GMT->current.proj.EQ_RAD * (GMT->current.proj.t_c1 * GMT->current.proj.t_lat0 + s2 * (GMT->current.proj.t_c2 + c2 * (GMT->current.proj.t_c3 + c2 * GMT->current.proj.t_c4)));
1198 GMT->current.proj.t_r = GMT->current.proj.EQ_RAD * GMT->current.setting.proj_scale_factor;
1199 GMT->current.proj.t_ir = 1.0 / GMT->current.proj.t_r;
1200 }
1201
1202 /* Ellipsoidal TM functions */
1203
gmtproj_tm(struct GMT_CTRL * P,double lon,double lat,double * x,double * y)1204 GMT_LOCAL void gmtproj_tm (struct GMT_CTRL *P, double lon, double lat, double *x, double *y) {
1205 /* Convert lon/lat to TM x/y */
1206 double N, T, T2, C, A, M, dlon, tan_lat, A2, A3, A5, s, c, s2, c2;
1207
1208 if (doubleAlmostEqual (fabs (lat), 90.0)) {
1209 M = P->current.proj.EQ_RAD * P->current.proj.t_c1 * M_PI_2;
1210 *x = 0.0;
1211 *y = P->current.setting.proj_scale_factor * M;
1212 }
1213 else {
1214 lat *= D2R;
1215 sincos (lat, &s, &c);
1216 sincos (2.0 * lat, &s2, &c2);
1217 tan_lat = s / c;
1218 M = P->current.proj.EQ_RAD * (P->current.proj.t_c1 * lat + s2 * (P->current.proj.t_c2 + c2 * (P->current.proj.t_c3 + c2 * P->current.proj.t_c4)));
1219 gmt_M_set_delta_lon (P->current.proj.central_meridian, lon, dlon);
1220 N = P->current.proj.EQ_RAD / d_sqrt (1.0 - P->current.proj.ECC2 * s * s);
1221 T = tan_lat * tan_lat;
1222 T2 = T * T;
1223 C = P->current.proj.t_e2 * c * c;
1224 A = dlon * D2R * c;
1225 A2 = A * A; A3 = A2 * A; A5 = A3 * A2;
1226 *x = P->current.setting.proj_scale_factor * N * (A + (1.0 - T + C) * (A3 * 0.16666666666666666667)
1227 + (5.0 - 18.0 * T + T2 + 72.0 * C - 58.0 * P->current.proj.t_e2) * (A5 * 0.00833333333333333333));
1228 A3 *= A; A5 *= A;
1229 *y = P->current.setting.proj_scale_factor * (M - P->current.proj.t_M0 + N * tan_lat * (0.5 * A2 + (5.0 - T + 9.0 * C + 4.0 * C * C) * (A3 * 0.04166666666666666667)
1230 + (61.0 - 58.0 * T + T2 + 600.0 * C - 330.0 * P->current.proj.t_e2) * (A5 * 0.00138888888888888889)));
1231 }
1232 }
1233
gmtproj_itm(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1234 GMT_LOCAL void gmtproj_itm (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1235 /* Convert TM x/y to lon/lat */
1236 double M, mu, u2, s, c, phi1, C1, C12, T1, T12, tmp, tmp2, N1, R_1, D, D2, D3, D5, tan_phi1, cp2;
1237
1238 M = y / GMT->current.setting.proj_scale_factor + GMT->current.proj.t_M0;
1239 mu = M * GMT->current.proj.t_i1;
1240
1241 u2 = 2.0 * mu;
1242 sincos (u2, &s, &c);
1243 phi1 = mu + s * (GMT->current.proj.t_i2 + c * (GMT->current.proj.t_i3 + c * (GMT->current.proj.t_i4 + c * GMT->current.proj.t_i5)));
1244
1245 sincos (phi1, &s, &c);
1246 tan_phi1 = s / c;
1247 cp2 = c * c;
1248 C1 = GMT->current.proj.t_e2 * cp2;
1249 C12 = C1 * C1;
1250 T1 = tan_phi1 * tan_phi1;
1251 T12 = T1 * T1;
1252 tmp = 1.0 - GMT->current.proj.ECC2 * (1.0 - cp2);
1253 tmp2 = d_sqrt (tmp);
1254 N1 = GMT->current.proj.EQ_RAD / tmp2;
1255 R_1 = GMT->current.proj.EQ_RAD * GMT->current.proj.one_m_ECC2 / (tmp * tmp2);
1256 D = x / (N1 * GMT->current.setting.proj_scale_factor);
1257 D2 = D * D; D3 = D2 * D; D5 = D3 * D2;
1258
1259 *lon = GMT->current.proj.central_meridian + R2D * (D - (1.0 + 2.0 * T1 + C1) * (D3 * 0.16666666666666666667)
1260 + (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C12 + 8.0 * GMT->current.proj.t_e2 + 24.0 * T12)
1261 * (D5 * 0.00833333333333333333)) / c;
1262 D3 *= D; D5 *= D;
1263 *lat = phi1 - (N1 * tan_phi1 / R_1) * (0.5 * D2 -
1264 (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C12 - 9.0 * GMT->current.proj.t_e2) * (D3 * 0.04166666666666666667)
1265 + (61.0 + 90.0 * T1 + 298 * C1 + 45.0 * T12 - 252.0 * GMT->current.proj.t_e2 - 3.0 * C12) * (D5 * 0.00138888888888888889));
1266 (*lat) *= R2D;
1267 }
1268
1269 /*Spherical TM functions */
1270
gmtproj_tm_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1271 GMT_LOCAL void gmtproj_tm_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1272 /* Convert lon/lat to TM x/y by spherical formula */
1273 double dlon, b, clat, slat, clon, slon, xx, yy;
1274
1275 gmt_M_set_delta_lon (GMT->current.proj.central_meridian, lon, dlon);
1276 if (fabs (lat) > 90.0) {
1277 /* Invalid latitude. Treat as in gmtproj_merc_sph(), but transversely: */
1278 *x = copysign (1.0e100, dlon);
1279 *y = 0.0;
1280 return;
1281 }
1282
1283 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_latc (GMT, lat);
1284
1285 sincosd (lat, &slat, &clat);
1286 sincosd (dlon, &slon, &clon);
1287 b = clat * slon;
1288 if (fabs(b) >= 1.0) {
1289 /* This corresponds to the transverse "pole"; the point at x = +-infinity, y = -lat0.
1290 Treat as in gmtproj_merc_sph(), but transversely: */
1291 *x = copysign (1.0e100, dlon);
1292 *y = -GMT->current.proj.t_r * GMT->current.proj.t_lat0;
1293 return;
1294 }
1295
1296 xx = atanh (b);
1297
1298 /* this should get us "over the pole";
1299 see not Snyder's formula but his example Fig. 10 on p. 50: */
1300
1301 yy = atan2 (slat, (clat * clon)) - GMT->current.proj.t_lat0;
1302 if (yy < -M_PI_2) yy += TWO_PI;
1303 *x = GMT->current.proj.t_r * xx;
1304 *y = GMT->current.proj.t_r * yy;
1305 }
1306
gmtproj_itm_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1307 GMT_LOCAL void gmtproj_itm_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1308 /* Convert TM x/y to lon/lat by spherical approximation. */
1309
1310 double xx, yy, sinhxx, coshxx, sind, cosd, lambda, phi;
1311
1312 xx = x * GMT->current.proj.t_ir;
1313 yy = y * GMT->current.proj.t_ir + GMT->current.proj.t_lat0;
1314
1315 sinhxx = sinh (xx);
1316 coshxx = cosh (xx);
1317
1318 sincos (yy, &sind, &cosd);
1319 phi = asind (sind / coshxx);
1320 *lat = phi;
1321
1322 lambda = atan2d (sinhxx, cosd);
1323 *lon = lambda + GMT->current.proj.central_meridian;
1324 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_latc_to_latg (GMT, *lat);
1325 }
1326
1327 /* -JU UNIVERSAL TRANSVERSE MERCATOR PROJECTION */
1328
1329 /* Ellipsoidal UTM */
1330
gmtproj_utm(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1331 GMT_LOCAL void gmtproj_utm (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1332 /* Convert lon/lat to UTM x/y */
1333
1334 if (lon < 0.0) lon += 360.0;
1335 gmtproj_tm (GMT, lon, lat, x, y);
1336 (*x) += GMT_FALSE_EASTING;
1337 if (!GMT->current.proj.north_pole) (*y) += GMT_FALSE_NORTHING; /* For S hemisphere, add 10^7 m */
1338 }
1339
gmtproj_iutm(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1340 GMT_LOCAL void gmtproj_iutm (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1341 /* Convert UTM x/y to lon/lat */
1342
1343 x -= GMT_FALSE_EASTING;
1344 if (!GMT->current.proj.north_pole) y -= GMT_FALSE_NORTHING;
1345 gmtproj_itm (GMT, lon, lat, x, y);
1346 }
1347
1348 /* Spherical UTM */
1349
gmtproj_utm_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1350 GMT_LOCAL void gmtproj_utm_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1351 /* Convert lon/lat to UTM x/y */
1352
1353 if (lon < 0.0) lon += 360.0;
1354 gmtproj_tm_sph (GMT, lon, lat, x, y);
1355 (*x) += GMT_FALSE_EASTING;
1356 if (!GMT->current.proj.north_pole) (*y) += GMT_FALSE_NORTHING; /* For S hemisphere, add 10^7 m */
1357 }
1358
gmtproj_iutm_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1359 GMT_LOCAL void gmtproj_iutm_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1360 /* Convert UTM x/y to lon/lat */
1361
1362 x -= GMT_FALSE_EASTING;
1363 if (!GMT->current.proj.north_pole) y -= GMT_FALSE_NORTHING;
1364 gmtproj_itm_sph (GMT, lon, lat, x, y);
1365 }
1366
1367 /* -JA LAMBERT AZIMUTHAL EQUAL AREA PROJECTION */
1368
gmtproj_vlambeq(struct GMT_CTRL * GMT,double lon0,double lat0,double horizon)1369 GMT_LOCAL void gmtproj_vlambeq (struct GMT_CTRL *GMT, double lon0, double lat0, double horizon) {
1370 /* Set up Spherical Lambert Azimuthal Equal-Area projection */
1371
1372 GMT->current.proj.central_meridian = lon0;
1373 GMT->current.proj.pole = lat0;
1374 if (GMT->current.proj.GMT_convert_latitudes) lat0 = gmt_M_latg_to_lata (GMT, lat0);
1375 sincosd (lat0, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp));
1376 GMT->current.proj.f_horizon = horizon;
1377 GMT->current.proj.rho_max = 2.0 * sind (0.5 * horizon) * GMT->current.proj.EQ_RAD;
1378 }
1379
gmtproj_lambeq(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1380 GMT_LOCAL void gmtproj_lambeq (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1381 /* Convert lon/lat to Spherical Lambert Azimuthal Equal-Area x/y */
1382 double k, tmp, sin_lat, cos_lat, sin_lon, cos_lon, c;
1383
1384 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1385 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
1386 sincosd (lat, &sin_lat, &cos_lat);
1387 sincosd (lon, &sin_lon, &cos_lon);
1388 c = cos_lat * cos_lon;
1389 tmp = 1.0 + GMT->current.proj.sinp * sin_lat + GMT->current.proj.cosp * c;
1390 if (tmp > 0.0) {
1391 k = GMT->current.proj.EQ_RAD * d_sqrt (2.0 / tmp);
1392 *x = k * cos_lat * sin_lon;
1393 *y = k * (GMT->current.proj.cosp * sin_lat - GMT->current.proj.sinp * c);
1394 if (GMT->current.proj.GMT_convert_latitudes) { /* Gotta fudge abit */
1395 (*x) *= GMT->current.proj.Dx;
1396 (*y) *= GMT->current.proj.Dy;
1397 }
1398 }
1399 else
1400 *x = *y = -DBL_MAX;
1401 }
1402
gmtproj_ilambeq(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1403 GMT_LOCAL void gmtproj_ilambeq (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1404 /* Convert Lambert Azimuthal Equal-Area x/y to lon/lat */
1405 double rho, a, sin_c, cos_c;
1406
1407 if (GMT->current.proj.GMT_convert_latitudes) { /* Undo effect of fudge factors */
1408 x *= GMT->current.proj.iDx;
1409 y *= GMT->current.proj.iDy;
1410 }
1411
1412 rho = hypot (x, y);
1413 if (rho > GMT->current.proj.rho_max) { /* Horizon */
1414 *lat = *lon = GMT->session.d_NaN;
1415 return;
1416 }
1417 a = 0.5 * rho * GMT->current.proj.i_EQ_RAD; /* a = sin(c/2) */
1418 a *= a; /* a = sin(c/2)**2 */
1419 cos_c = 1.0 - 2.0 * a; /* cos_c = cos(c) */
1420 sin_c = sqrt (1.0 - a) * GMT->current.proj.i_EQ_RAD; /* sin_c = sin(c)/rho */
1421 *lat = d_asind (cos_c * GMT->current.proj.sinp + y * sin_c * GMT->current.proj.cosp);
1422 *lon = d_atan2d (x * sin_c, GMT->current.proj.cosp * cos_c - y * GMT->current.proj.sinp * sin_c) + GMT->current.proj.central_meridian;
1423 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
1424 }
1425
1426 /* -JG GENERAL PERSPECTIVE PROJECTION */
1427
1428 /* Set up General Perspective projection */
1429
1430 /* Convert lon/lat to General Perspective x/y */
1431
gmtproj_genper(struct GMT_CTRL * GMT,double lon,double lat,double * xt,double * yt)1432 GMT_LOCAL void gmtproj_genper (struct GMT_CTRL *GMT, double lon, double lat, double *xt, double *yt) {
1433 double dlon, sin_lat, cos_lat, sin_dlon, cos_dlon;
1434 double cosc, sinc;
1435 double x, y, kp;
1436 double angle;
1437
1438 dlon = lon - GMT->current.proj.central_meridian;
1439 while (dlon < -GMT_180) dlon += 360.0;
1440 while (dlon > 180.0) dlon -= 360.0;
1441 dlon *= D2R;
1442
1443 lat = gmtproj_genper_getgeocentric (GMT, lat, 0.0);
1444
1445 sincosd (lat, &sin_lat, &cos_lat);
1446 sincos (dlon, &sin_dlon, &cos_dlon);
1447
1448 cosc = GMT->current.proj.sinp * sin_lat + GMT->current.proj.cosp * cos_lat * cos_dlon;
1449 sinc = d_sqrt(1.0 - cosc * cosc);
1450
1451 GMT->current.proj.g_outside = false;
1452
1453 if (cosc < GMT->current.proj.g_P_inverse) { /* over the horizon */
1454 GMT->current.proj.g_outside = true;
1455
1456 if (GMT->current.proj.polar)
1457 angle = M_PI - dlon;
1458 else if (GMT->current.proj.cosp*sinc != 0.0) {
1459 angle = d_acos((sin_lat - GMT->current.proj.sinp*cosc)/(GMT->current.proj.cosp*sinc));
1460 if (dlon < 0.0) angle = -angle;
1461 }
1462 else
1463 angle = 0.0;
1464
1465 sincos (angle, &x, &y);
1466 x *= GMT->current.proj.g_rmax;
1467 y *= GMT->current.proj.g_rmax;
1468 angle *= R2D;
1469 }
1470 else if (GMT->current.proj.ECC2 != 0.0) { /* within field of view, ellipsoidal earth */
1471 gmtproj_genper_toxy (GMT, lat, lon, 0.0, &x, &y);
1472 /* angle = GMT->current.proj.g_azimuth; */
1473 angle = atan2d(x, y);
1474 /* XXX Which one is it? Forgotten R2D. Switched x and y. */
1475 }
1476 else { /* within field of view, spherical earth */
1477 kp = GMT->current.proj.g_R * (GMT->current.proj.g_P - 1.0) / (GMT->current.proj.g_P - cosc);
1478 x = kp * cos_lat * sin_dlon;
1479 y = kp * (GMT->current.proj.cosp * sin_lat - GMT->current.proj.sinp * cos_lat * cos_dlon);
1480 /* angle = GMT->current.proj.g_azimuth; */
1481 angle = atan2d(x, y);
1482 /* XXX Which one is it? Forgotten R2D. Switched x and y. */
1483 }
1484
1485 gmtproj_genper_to_xtyt (GMT, angle, x, y, GMT->current.proj.g_yoffset, xt, yt);
1486
1487 if (gmt_M_is_dnan(*yt) || gmt_M_is_dnan(*xt)) {
1488 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: yt or xt nan\n");
1489 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: lon %6.3f lat %6.3f\n", lon, lat);
1490 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper: xt %10.3e yt %10.3e\n", *xt, *yt);
1491 }
1492 }
1493
gmtproj_vgenper(struct GMT_CTRL * GMT,double lon0,double lat0,double altitude,double azimuth,double tilt,double twist,double width,double height)1494 GMT_LOCAL void gmtproj_vgenper (struct GMT_CTRL *GMT, double lon0, double lat0, double altitude, double azimuth, double tilt, double twist, double width, double height) {
1495 double R, Req, Rpolar, Rlat0;
1496 double H, P, PP;
1497 double rho, rho2;
1498 double eca, cos_eca, sin_eca ;
1499 double yt_min, yt_max;
1500 double xt_min, xt_max;
1501 double rmax;
1502 double rmax_at_lat0, rmax_min, rmax_max;
1503 double t1, t2, t12, t22;
1504 double kp;
1505 double xt, yt;
1506 double x, y;
1507 double gamma, Omega, sinOmega, cosOmega, omega_max;
1508 double sinlatvp, coslatvp, latvp, lonvp;
1509 double xt_vp, yt_vp;
1510 double eccen;
1511 double max_yt;
1512 double sin_lat0, cos_lat0, sin_lat1, cos_lat1, dlong;
1513 double lat0_save;
1514 double vp_lat, vp_long;
1515 int az;
1516
1517 GMT->current.proj.central_meridian = lon0;
1518
1519 lat0_save = lat0;
1520
1521 Req = R = GMT->current.proj.EQ_RAD;
1522 Rpolar = Req * sqrt(GMT->current.proj.one_m_ECC2);
1523
1524 sincosd (lat0, &t2, &t1);
1525
1526 t1 *= R;
1527 t12 = R*t1;
1528 t2 *= Rpolar;
1529 t22 = Rpolar*t2;
1530
1531 Rlat0 = sqrt ((t12*t12 + t22*t22)/(t1*t1 + t2*t2));
1532
1533 lat0 = gmtproj_genper_getgeocentric (GMT, lat0, 0.0);
1534 sincosd (lat0, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp));
1535
1536 if (!gmt_M_proj_is_zero (GMT->current.proj.ECC2)) {
1537 gmtproj_genper_setup (GMT, 0.0, altitude, lat0_save, lon0);
1538 GMT->current.proj.central_meridian = lon0;
1539 GMT->current.proj.pole = GMT->current.proj.g_phig;
1540 }
1541 else {
1542 GMT->current.proj.pole = lat0;
1543
1544 if (GMT->current.proj.g_radius || (altitude < -10.0)) {
1545 /* use altitude as the radial distance from the center of the earth*/
1546 H = fabs(altitude*1e3) - R;
1547 P = H/R + 1.0;
1548 }
1549 else if (altitude <= 0.0) {
1550 /* compute altitude of geosynchronous viewpoint n*/
1551 double temp = 86164.1/TWO_PI;
1552 H = pow(3.98603e14*temp*temp, 0.3333) - R;
1553 P = H/R + 1.0;
1554 }
1555 else if (altitude < 10.0) {
1556 P = altitude;
1557 H = R * (P - 1.0);
1558 }
1559 else {
1560 H = altitude*1e3;
1561 P = H/R + 1.0;
1562 }
1563 GMT->current.proj.g_R = R;
1564 GMT->current.proj.g_H = H;
1565 GMT->current.proj.g_P = P;
1566 }
1567
1568 H = GMT->current.proj.g_H;
1569 P = GMT->current.proj.g_P;
1570 R = GMT->current.proj.g_R;
1571
1572 GMT->current.proj.g_P_inverse = P > 0.0 ? 1.0/P : 1.0;
1573
1574 if (GMT->current.proj.g_longlat_set) {
1575 double norm_long = lon0;
1576
1577 vp_lat = tilt;
1578 vp_long = azimuth;
1579
1580 if (vp_lat == lat0_save && vp_long == lon0)
1581 tilt = azimuth = 0.0;
1582 else {
1583 if (GMT->current.proj.g_debug > 0) {
1584 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: sensor point long %7.4f lat %7.4f\n", lon0, lat0);
1585 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: input view point long %7.4f lat %7.4f\n", vp_long, vp_lat);
1586 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: input twist %7.4f\n", twist);
1587 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: altitude %f H %f R %f P %7.4f\n", altitude, H/1000.0, R/1000.0,P);
1588 }
1589
1590 sincosd (90.0 - lat0, &sin_lat0, &cos_lat0);
1591 sincosd (90.0 - vp_lat, &sin_lat1, &cos_lat1);
1592
1593 while (vp_long < 0.0) vp_long += 360.0;
1594 while (norm_long < 0.0) norm_long += 360.0;
1595
1596 dlong = vp_long - norm_long;
1597 if (dlong < -GMT_180) dlong += 360.0;
1598
1599 cos_eca = cos_lat0*cos_lat1 + sin_lat0*sin_lat1*cosd(dlong);
1600 eca = d_acos (cos_eca);
1601 sin_eca = sin (eca);
1602
1603 rho2 = P*P + 1.0 - 2.0*P*cos_eca;
1604 rho = sqrt (rho2);
1605
1606 tilt = d_acosd ((rho2 + P*P - 1.0)/(2.0*rho*P));
1607
1608 azimuth = d_acosd ((cos_lat1 - cos_lat0*cos_eca)/(sin_lat0*sin_eca));
1609
1610 if (dlong < 0) azimuth = 360.0 - azimuth;
1611
1612 }
1613 if (GMT->current.proj.g_debug > 0) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: pointing at longitude %10.4f latitude %10.4f\n with computed tilt %5.2f azimuth %6.2f\n", vp_long, vp_lat, tilt, azimuth);
1614
1615 }
1616 else if (GMT->current.proj.g_debug > 1) {
1617 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: sensor point long %6.3f lat %6.3f\n", lon0, lat0);
1618 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: input azimuth %6.3f tilt %6.3f\n", azimuth, tilt);
1619 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: input twist %6.3f\n", twist);
1620 }
1621
1622 if (tilt < 0.0) tilt = d_asind (GMT->current.proj.g_P_inverse);
1623
1624 sincosd (tilt, &(GMT->current.proj.g_sin_tilt), &(GMT->current.proj.g_cos_tilt));
1625 sincosd (twist, &(GMT->current.proj.g_sin_twist), &(GMT->current.proj.g_cos_twist));
1626
1627 GMT->current.proj.g_box = !(fabs (width) < GMT_CONV4_LIMIT);
1628
1629 if (width != 0.0 && height == 0) height = width;
1630 if (height != 0.0 && width == 0) width = height;
1631 GMT->current.proj.g_width = width/2.0;
1632
1633 GMT->current.proj.g_azimuth = azimuth;
1634 sincosd (azimuth, &(GMT->current.proj.g_sin_azimuth), &(GMT->current.proj.g_cos_azimuth));
1635
1636 PP = sqrt ((P - 1.0)/(P + 1.0));
1637 rmax = R*PP;
1638 rmax_min = Rpolar*PP;
1639 rmax_max = Req*PP;
1640 rmax_at_lat0 = Rlat0*PP;
1641
1642 if (GMT->current.proj.ECC2 != 0.0) rmax = rmax_at_lat0;
1643
1644 kp = R*(P - 1.0) / (P - GMT->current.proj.g_P_inverse);
1645
1646 omega_max = d_acosd(GMT->current.proj.g_P_inverse);
1647 GMT->current.proj.g_rmax = rmax;
1648
1649 /* GMT->current.proj.f_horizon = GMT->current.proj.g_P_inverse; */
1650 GMT->current.proj.f_horizon = omega_max;
1651 /* XXX Which one is it ? */
1652
1653 max_yt = 2.0*R*sind (0.5*omega_max);
1654
1655 eccen = sind (tilt)/sqrt (1.0 - 1.0/(P*P));
1656
1657 gamma = 180.0 - d_asind (GMT->current.proj.g_sin_tilt * P);
1658 Omega = 180.0 - tilt - gamma;
1659
1660 if (GMT->current.proj.g_debug > 0)
1661 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: tilt %6.3f sin_tilt %10.6f P %6.4f gamma %6.4f\n Omega %6.4f eccen %10.4f\n",
1662 tilt, GMT->current.proj.g_sin_tilt, P, gamma, Omega, eccen);
1663
1664 if (eccen == 1.0) {
1665 max_yt = MIN (max_yt, rmax * 2.0);
1666 if (GMT->current.proj.g_debug > 1)
1667 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Projected map is a parabola with requested tilt %6.3f\n max ECA is %6.3f degrees.\n Plot truncated for projected distances > rmax %8.2f\n", tilt, omega_max, rmax/1000.0);
1668 }
1669 else if (eccen > 1.0) {
1670 if (width != 0.0) {
1671 if (GMT->current.proj.g_debug > 1)
1672 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Projected map is a hyperbola with requested tilt %6.3f\n max ECA is %6.3f degrees.\n", tilt, omega_max);
1673 }
1674 else {
1675 max_yt = MIN (max_yt, rmax * 2.0);
1676 if (GMT->current.proj.g_debug > 1)
1677 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Projected map is a hyperbola with requested tilt %6.3f\n max ECA is %6.3f degrees.\n Plot truncated for projected distances > rmax %8.2f\n", tilt, omega_max, rmax/1000.0);
1678 }
1679 }
1680 else if (eccen > 0.5) {
1681 if (width != 0.0) {
1682 double t = sind (tilt), Pecc, maxecc = 0.5;
1683 Pecc = sqrt (1.0/(1.0 - (t*t/maxecc)));
1684 max_yt = R*sqrt ((Pecc-1.0)/(Pecc+1.0));
1685 if (GMT->current.proj.g_debug > 1)
1686 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Projected map is an enlongated ellipse (eccentricity of %6.4f) with "
1687 "requested tilt %6.3f\nwill truncate plot at rmax %8.2f\n", eccen, tilt, max_yt);
1688 }
1689 else {
1690 if (max_yt > rmax *2.0) max_yt = rmax * 2.0;
1691 if (GMT->current.proj.g_debug > 1)
1692 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Projected map is an enlongated ellipse with requested tilt %6.3f\n eccentricity %6.3f\n Plot truncated for projected distances > rmax %8.2f\n", tilt, eccen, rmax/1000.0);
1693 }
1694 }
1695
1696 GMT->current.proj.g_max_yt = max_yt;
1697
1698 sincosd (Omega, &sinOmega, &cosOmega);
1699
1700 rho2 = P*P + 1.0 - 2.0*P*cosOmega;
1701 rho = sqrt (rho2);
1702
1703 sinlatvp = GMT->current.proj.sinp*cosOmega + GMT->current.proj.cosp*sinOmega*GMT->current.proj.g_cos_azimuth;
1704 latvp = d_asind (sinlatvp);
1705 coslatvp = sqrt (1.0 - sinlatvp*sinlatvp);
1706
1707 lonvp = acosd ((cosOmega - sinlatvp*GMT->current.proj.sinp)/(GMT->current.proj.cosp*coslatvp));
1708 if (azimuth > 180.0) lonvp = -lonvp;
1709 lonvp += lon0;
1710
1711 if (GMT->current.proj.g_debug > 1)
1712 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: pointing at longitude %10.4f latitude %10.4f\n "
1713 "with tilt %5.2f azimuth %6.2f at distance %6.4f\n "
1714 "with width %6.3f height %6.3f twist %6.2f\n", lonvp, latvp, tilt, azimuth, rho, width, height, twist);
1715
1716 GMT->current.proj.g_yoffset = 0.0;
1717
1718 if (height != 0.0) {
1719 GMT->current.proj.g_yoffset = GMT->current.proj.g_sin_tilt * H ;
1720 xt_max = R * rho * sind (0.5*width);
1721 xt_min = -xt_max;
1722 yt_max = R * rho * sind (0.5*height);
1723 yt_min = -yt_max;
1724 }
1725 else {
1726 FILE *fp = NULL;
1727 xt_min = 1e20;
1728 xt_max = -xt_min;
1729
1730 yt_min = 1e20;
1731 yt_max = -yt_min;
1732
1733 if (GMT->current.proj.g_debug > 2) {
1734 if ((fp = fopen("g_border.txt", "w")) == NULL) {
1735 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: Failed to write the g_border.txt file\n");
1736 }
1737 else {
1738 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: tilt %10.4f sin_tilt %10.4f cos_tilt %10.4f\n",
1739 tilt, GMT->current.proj.g_sin_tilt, GMT->current.proj.g_cos_tilt);
1740 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: azimuth %10.4f sin_azimuth %10.4f cos_azimuth %10.4f\n",
1741 azimuth, GMT->current.proj.g_sin_azimuth, GMT->current.proj.g_cos_azimuth);
1742 }
1743 }
1744
1745 for (az = 0 ; az < 360 ; az++) {
1746 sincosd ((double)az, &x, &y);
1747 x *= rmax;
1748 y *= rmax;
1749 gmtproj_genper_to_xtyt (GMT, (double)az, x, y, GMT->current.proj.g_yoffset, &xt, &yt);
1750
1751 if (fp && GMT->current.proj.g_debug > 2)
1752 fprintf (fp,"%3d x %10.2f y %10.2f xt %10.3f yt %10.3f\n", az, x/1000, y/1000, xt/1000, yt/1000);
1753
1754 xt_min = MIN (xt, xt_min);
1755 xt_max = MAX (xt, xt_max);
1756 yt_min = MIN (yt, yt_min);
1757 yt_max = MAX (yt, yt_max);
1758 }
1759 if (GMT->current.proj.g_debug > 2 && fp) fclose(fp);
1760 if (eccen > 0.5) {
1761 GMT->current.proj.g_width = atand (2.0*rmax/H);
1762 height = width = 2.0 * GMT->current.proj.g_width;
1763 xt_max = yt_max = R * rho * sind (GMT->current.proj.g_width);
1764 xt_min = yt_min = -xt_max;
1765 }
1766 }
1767 if (GMT->current.proj.g_debug > 1) {
1768 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: xt max %7.1f km\n", xt_max/1000.0);
1769 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: xt min %7.1f km\n", xt_min/1000.0);
1770 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: yt max %7.1f km\n", yt_max/1000.0);
1771 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: yt min %7.1f km\n", yt_min/1000.0);
1772 }
1773
1774 GMT->current.proj.g_xmin = xt_min;
1775 GMT->current.proj.g_xmax = xt_max;
1776 GMT->current.proj.g_ymin = yt_min;
1777 GMT->current.proj.g_ymax = yt_max;
1778
1779 if (width != 0.0) GMT->current.proj.scale[GMT_Y] = GMT->current.proj.scale[GMT_X]/width*height;
1780
1781 if (GMT->current.proj.g_debug > 0) {
1782 gmtproj_genper (GMT, lonvp, latvp, &xt_vp, &yt_vp);
1783 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: polar %d north %d\n", GMT->current.proj.polar, GMT->current.proj.north_pole);
1784 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: altitude H %7.1f km P %7.4f\n", H/1000.0, P);
1785 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: azimuth %5.1f tilt %5.1f\n", azimuth, tilt);
1786 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: viewpoint width %5.1f height %5.1f degrees\n", width, height);
1787 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: radius max %7.1f km\n", GMT->current.proj.g_rmax/1000.0);
1788 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: eccentricity %7.4f km\n", eccen);
1789 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: eq radius max %7.1f km\n", rmax_max/1000.0);
1790 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: polar radius max %7.1f km\n", rmax_min/1000.0);
1791 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: lat0 radius max %7.1f km\n", rmax_at_lat0/1000.0);
1792 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: kp %7.1f \n", kp/1000.0);
1793 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: y offset %7.1f km\n", GMT->current.proj.g_yoffset/1000.0);
1794 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: yt max %7.1f km\n", yt_max/1000.0);
1795 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: yt min %7.1f km\n", yt_min/1000.0);
1796 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: y max %7.1f km\n", GMT->current.proj.g_ymax/1000.0);
1797 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: y min %7.1f km\n", GMT->current.proj.g_ymin/1000.0);
1798 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: x max %7.1f km\n", GMT->current.proj.g_xmax/1000.0);
1799 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: x min %7.1f km\n", GMT->current.proj.g_xmin/1000.0);
1800 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: omega max %6.2f degrees\n", omega_max);
1801 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: gamma %6.3f Omega %6.3f \n", gamma, Omega);
1802 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: viewpoint lon %6.3f lat %6.3f \n", lonvp, latvp);
1803 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: viewpoint xt %6.3f yt %6.3f \n", xt_vp/1000.0, yt_vp/1000.0);
1804 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "vgenper: user viewpoint %d\n", GMT->current.proj.g_box);
1805 }
1806
1807 }
1808
1809
1810 /* Convert General Perspective x/y to lon/lat */
1811
gmtproj_igenper(struct GMT_CTRL * GMT,double * lon,double * lat,double xt,double yt)1812 GMT_LOCAL void gmtproj_igenper (struct GMT_CTRL *GMT, double *lon, double *lat, double xt, double yt) {
1813 double H, P, R;
1814 double sin_c, cos_c;
1815 double x, y;
1816 double M, Q;
1817 double con, com, rho;
1818
1819 H = GMT->current.proj.g_H;
1820 R = GMT->current.proj.g_R;
1821 P = GMT->current.proj.g_P;
1822
1823 x = xt;
1824 y = yt;
1825
1826 xt = (x * GMT->current.proj.g_cos_twist + y * GMT->current.proj.g_sin_twist);
1827 yt = (y * GMT->current.proj.g_cos_twist - x * GMT->current.proj.g_sin_twist);
1828 yt += GMT->current.proj.g_yoffset;
1829
1830 M = H * xt / (H - yt * GMT->current.proj.g_sin_tilt);
1831 Q = H * yt * GMT->current.proj.g_cos_tilt /(H - yt * GMT->current.proj.g_sin_tilt);
1832 x = M * GMT->current.proj.g_cos_azimuth + Q * GMT->current.proj.g_sin_azimuth;
1833 y = Q * GMT->current.proj.g_cos_azimuth - M * GMT->current.proj.g_sin_azimuth;
1834
1835 rho = hypot(x, y);
1836
1837 GMT->current.proj.g_outside = false;
1838
1839 if (rho < GMT_CONV4_LIMIT) {
1840 *lat = GMT->current.proj.pole;
1841 *lon = GMT->current.proj.central_meridian;
1842 return;
1843 }
1844 if (rho > GMT->current.proj.g_rmax) {
1845 x *= GMT->current.proj.g_rmax/rho;
1846 y *= GMT->current.proj.g_rmax/rho;
1847 rho = GMT->current.proj.g_rmax;
1848 GMT->current.proj.g_outside = true;
1849 }
1850
1851 con = P - 1.0;
1852 com = P + 1.0;
1853
1854 if (GMT->current.proj.ECC2 != 0.0)
1855 gmtproj_genper_tolatlong (GMT, x, y, 0.0, lat, lon);
1856 else {
1857 sin_c = (P - d_sqrt(1.0 - (rho * rho * com)/(R*R*con))) / (R * con / rho + rho/(R*con));
1858 cos_c = d_sqrt(1.0 - sin_c * sin_c);
1859 sin_c /= rho;
1860 *lat = d_asind (cos_c * GMT->current.proj.sinp + y * sin_c * GMT->current.proj.cosp);
1861 *lon = d_atan2d (x * sin_c, cos_c * GMT->current.proj.cosp - y * sin_c * GMT->current.proj.sinp) + GMT->current.proj.central_meridian;
1862 }
1863 return;
1864 }
1865
gmtlib_genper_map_clip_path(struct GMT_CTRL * GMT,uint64_t np,double * work_x,double * work_y)1866 int gmtlib_genper_map_clip_path (struct GMT_CTRL *GMT, uint64_t np, double *work_x, double *work_y) {
1867 uint64_t i;
1868 double da, angle;
1869 double x, y, xt, yt;
1870
1871 if (GMT->current.proj.g_debug > 0) {
1872 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper_map_clip_path: np %" PRIu64 "\n", np);
1873 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "genper_map_clip_path: x_scale %e y_scale %e, x0 %e y0 %e\n", GMT->current.proj.scale[GMT_X], GMT->current.proj.scale[GMT_Y], GMT->current.proj.origin[GMT_X], GMT->current.proj.origin[GMT_Y]);
1874 }
1875 assert (np > 1);
1876 da = TWO_PI/(np-1);
1877
1878 for (i = 0; i < np; i++) {
1879 angle = i * da;
1880 sincos (angle, &x, &y);
1881 x *= GMT->current.proj.g_rmax;
1882 y *= GMT->current.proj.g_rmax;
1883
1884 /* XXX forgotten R2D */
1885 gmtproj_genper_to_xtyt (GMT, angle*R2D, x, y, GMT->current.proj.g_yoffset, &xt, &yt);
1886
1887 if (GMT->current.proj.g_width != 0.0) {
1888 xt = MAX (GMT->current.proj.g_xmin, MIN (xt, GMT->current.proj.g_xmax));
1889 yt = MAX (GMT->current.proj.g_ymin, MIN (yt, GMT->current.proj.g_ymax));
1890 }
1891 work_x[i] = xt * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X];
1892 work_y[i] = yt * GMT->current.proj.scale[GMT_Y] + GMT->current.proj.origin[GMT_Y];
1893 }
1894 return 0;
1895 }
1896
1897 /* -JG ORTHOGRAPHIC PROJECTION */
1898
gmtproj_vortho(struct GMT_CTRL * GMT,double lon0,double lat0,double horizon)1899 GMT_LOCAL void gmtproj_vortho (struct GMT_CTRL *GMT, double lon0, double lat0, double horizon) {
1900 /* Set up Orthographic projection */
1901
1902 GMT->current.proj.central_meridian = lon0;
1903 GMT->current.proj.pole = lat0;
1904 sincosd (lat0, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp));
1905 GMT->current.proj.f_horizon = horizon;
1906 GMT->current.proj.rho_max = sind (horizon);
1907 }
1908
gmtproj_ortho(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1909 GMT_LOCAL void gmtproj_ortho (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1910 /* Convert lon/lat to Orthographic x/y */
1911 double sin_lat, cos_lat, sin_lon, cos_lon;
1912
1913 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1914
1915 sincosd (lat, &sin_lat, &cos_lat);
1916 sincosd (lon, &sin_lon, &cos_lon);
1917 *x = GMT->current.proj.EQ_RAD * cos_lat * sin_lon;
1918 *y = GMT->current.proj.EQ_RAD * (GMT->current.proj.cosp * sin_lat - GMT->current.proj.sinp * cos_lat * cos_lon);
1919 }
1920
gmtproj_iortho(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1921 GMT_LOCAL void gmtproj_iortho (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1922 /* Convert Orthographic x/y to lon/lat */
1923 double rho, cos_c;
1924
1925 x *= GMT->current.proj.i_EQ_RAD;
1926 y *= GMT->current.proj.i_EQ_RAD;
1927 rho = hypot (x, y);
1928 if (rho > GMT->current.proj.rho_max) {
1929 *lat = *lon = GMT->session.d_NaN;
1930 return;
1931 }
1932 cos_c = sqrt (1.0 - rho * rho); /* Produces NaN for rho > 1: beyond horizon */
1933 *lat = d_asind (cos_c * GMT->current.proj.sinp + y * GMT->current.proj.cosp);
1934 *lon = d_atan2d (x, cos_c * GMT->current.proj.cosp - y * GMT->current.proj.sinp) + GMT->current.proj.central_meridian;
1935 }
1936
1937 /* -JF GNOMONIC PROJECTION */
1938
gmtproj_vgnomonic(struct GMT_CTRL * GMT,double lon0,double lat0,double horizon)1939 GMT_LOCAL void gmtproj_vgnomonic (struct GMT_CTRL *GMT, double lon0, double lat0, double horizon) {
1940 /* Set up Gnomonic projection */
1941
1942 GMT->current.proj.central_meridian = lon0;
1943 GMT->current.proj.f_horizon = horizon;
1944 GMT->current.proj.rho_max = tand (GMT->current.proj.f_horizon);
1945 GMT->current.proj.pole = lat0;
1946 GMT->current.proj.north_pole = (lat0 > 0.0);
1947 sincosd (lat0, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp));
1948 }
1949
gmtproj_gnomonic(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1950 GMT_LOCAL void gmtproj_gnomonic (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1951 /* Convert lon/lat to Gnomonic x/y */
1952 double k, sin_lat, cos_lat, sin_lon, cos_lon, cc;
1953
1954 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1955 sincosd (lat, &sin_lat, &cos_lat);
1956 sincosd (lon, &sin_lon, &cos_lon);
1957
1958 cc = cos_lat * cos_lon;
1959 k = GMT->current.proj.EQ_RAD / (GMT->current.proj.sinp * sin_lat + GMT->current.proj.cosp * cc);
1960 *x = k * cos_lat * sin_lon;
1961 *y = k * (GMT->current.proj.cosp * sin_lat - GMT->current.proj.sinp * cc);
1962 }
1963
gmtproj_ignomonic(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)1964 GMT_LOCAL void gmtproj_ignomonic (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
1965 /* Convert Gnomonic x/y to lon/lat */
1966 double rho, c;
1967
1968 x *= GMT->current.proj.i_EQ_RAD;
1969 y *= GMT->current.proj.i_EQ_RAD;
1970 rho = hypot (x, y);
1971 if (rho > GMT->current.proj.rho_max) {
1972 *lat = *lon = GMT->session.d_NaN;
1973 return;
1974 }
1975 c = atan (rho);
1976 *lat = d_asind (cos(c) * (GMT->current.proj.sinp + y * GMT->current.proj.cosp));
1977 *lon = d_atan2d (x, GMT->current.proj.cosp - y * GMT->current.proj.sinp) + GMT->current.proj.central_meridian;
1978 }
1979
1980 /* -JE AZIMUTHAL EQUIDISTANT PROJECTION */
1981
gmtproj_vazeqdist(struct GMT_CTRL * GMT,double lon0,double lat0,double horizon)1982 GMT_LOCAL void gmtproj_vazeqdist (struct GMT_CTRL *GMT, double lon0, double lat0, double horizon) {
1983 /* Set up azimuthal equidistant projection */
1984
1985 GMT->current.proj.central_meridian = lon0;
1986 GMT->current.proj.pole = lat0;
1987 sincosd (lat0, &(GMT->current.proj.sinp), &(GMT->current.proj.cosp));
1988 GMT->current.proj.f_horizon = horizon;
1989 GMT->current.proj.rho_max = horizon * D2R;
1990 }
1991
gmtproj_azeqdist(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)1992 GMT_LOCAL void gmtproj_azeqdist (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
1993 /* Convert lon/lat to azimuthal equidistant x/y */
1994 double k, cc, c, clat, slon, clon, slat, t;
1995
1996 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
1997
1998 if (gmt_M_proj_is_zero (lat-GMT->current.proj.pole) && gmt_M_proj_is_zero (lon)) { /* Center of projection */
1999 *x = *y = 0.0;
2000 return;
2001 }
2002 sincosd (lat, &slat, &clat);
2003 sincosd (lon, &slon, &clon);
2004
2005 t = clat * clon;
2006 cc = GMT->current.proj.sinp * slat + GMT->current.proj.cosp * t;
2007 if (cc <= -1.0) { /* Antipode is a circle, so flag x,y as NaN and increase uint64_t */
2008 *x = *y = GMT->session.d_NaN;
2009 GMT->current.proj.n_antipoles++;
2010 }
2011 else {
2012 c = d_acos (cc);
2013 k = (gmt_M_proj_is_zero (c)) ? GMT->current.proj.EQ_RAD : GMT->current.proj.EQ_RAD * c / sin (c);
2014 *x = k * clat * slon;
2015 *y = k * (GMT->current.proj.cosp * slat - GMT->current.proj.sinp * t);
2016 }
2017 }
2018
gmtproj_iazeqdist(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2019 GMT_LOCAL void gmtproj_iazeqdist (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2020 /* Convert azimuthal equidistant x/y to lon/lat */
2021 double rho, sin_c, cos_c;
2022
2023 x *= GMT->current.proj.i_EQ_RAD;
2024 y *= GMT->current.proj.i_EQ_RAD;
2025 rho = hypot (x, y);
2026 if (rho > GMT->current.proj.rho_max) { /* Horizon */
2027 *lat = *lon = GMT->session.d_NaN;
2028 return;
2029 }
2030 sincos (rho, &sin_c, &cos_c);
2031 if (rho != 0.0) sin_c /= rho;
2032 /* Since in the rest we only have sin_c in combination with x or y, it is not necessary to set
2033 sin_c = 1 when rho == 0. In that case x*sin_c and y*sin_c or 0 anyhow. */
2034 *lat = d_asind (cos_c * GMT->current.proj.sinp + y * sin_c * GMT->current.proj.cosp);
2035 *lon = d_atan2d (x * sin_c, cos_c * GMT->current.proj.cosp - y * sin_c * GMT->current.proj.sinp) + GMT->current.proj.central_meridian;
2036 }
2037
2038 /* -JW MOLLWEIDE EQUAL AREA PROJECTION */
2039
gmtproj_vmollweide(struct GMT_CTRL * GMT,double lon0,double scale)2040 GMT_LOCAL void gmtproj_vmollweide (struct GMT_CTRL *GMT, double lon0, double scale) {
2041 /* Set up Mollweide Equal-Area projection */
2042
2043 gmtproj_check_R_J (GMT, &lon0);
2044 GMT->current.proj.central_meridian = lon0;
2045 GMT->current.proj.w_x = GMT->current.proj.EQ_RAD * D2R * d_sqrt (8.0) / M_PI;
2046 GMT->current.proj.w_y = GMT->current.proj.EQ_RAD * M_SQRT2;
2047 GMT->current.proj.w_iy = 1.0 / GMT->current.proj.w_y;
2048 GMT->current.proj.w_r = 0.25 * (scale * GMT->current.proj.M_PR_DEG * 360.0); /* = Half the minor axis */
2049 }
2050
gmtproj_mollweide(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2051 GMT_LOCAL void gmtproj_mollweide (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2052 /* Convert lon/lat to Mollweide Equal-Area x/y */
2053 int i;
2054 double phi, delta, psin_lat, c, s;
2055
2056 if (doubleAlmostEqual (fabs (lat), 90.0)) { /* Special case */
2057 *x = 0.0;
2058 *y = copysign (GMT->current.proj.w_y, lat);
2059 return;
2060 }
2061
2062 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2063 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
2064 lat *= D2R;
2065
2066 phi = lat;
2067 psin_lat = M_PI * sin (lat);
2068 i = 0;
2069 do {
2070 i++;
2071 sincos (phi, &s, &c);
2072 delta = -(phi + s - psin_lat) / (1.0 + c);
2073 phi += delta;
2074 }
2075 while (fabs (delta) > GMT_PROJ_CONV_LIMIT && i < GMT_PROJ_MAX_ITERATIONS);
2076 phi *= 0.5;
2077 sincos (phi, &s, &c);
2078 *x = GMT->current.proj.w_x * lon * c;
2079 *y = GMT->current.proj.w_y * s;
2080 }
2081
gmtproj_imollweide(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2082 GMT_LOCAL void gmtproj_imollweide (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2083 /* Convert Mollweide Equal-Area x/y to lon/lat */
2084 double phi, phi2;
2085
2086 phi = asin (y * GMT->current.proj.w_iy);
2087 *lon = x / (GMT->current.proj.w_x * cos(phi));
2088 if (fabs (*lon) > 180.0) { /* Horizon */
2089 *lat = *lon = GMT->session.d_NaN;
2090 return;
2091 }
2092 *lon += GMT->current.proj.central_meridian;
2093 phi2 = 2.0 * phi;
2094 *lat = asind ((phi2 + sin (phi2)) / M_PI);
2095 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
2096 }
2097
2098 /* -JH HAMMER-AITOFF EQUAL AREA PROJECTION */
2099
gmtproj_vhammer(struct GMT_CTRL * GMT,double lon0,double scale)2100 GMT_LOCAL void gmtproj_vhammer (struct GMT_CTRL *GMT, double lon0, double scale) {
2101 /* Set up Hammer-Aitoff Equal-Area projection */
2102
2103 gmtproj_check_R_J (GMT, &lon0);
2104 GMT->current.proj.central_meridian = lon0;
2105 GMT->current.proj.w_r = 0.25 * (scale * GMT->current.proj.M_PR_DEG * 360.0); /* = Half the minor axis */
2106 }
2107
gmtproj_hammer(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2108 GMT_LOCAL void gmtproj_hammer (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2109 /* Convert lon/lat to Hammer-Aitoff Equal-Area x/y */
2110 double slat, clat, slon, clon, D;
2111
2112 if (doubleAlmostEqual (fabs (lat), 90.0)) { /* Save time */
2113 *x = 0.0;
2114 *y = M_SQRT2 * copysign (GMT->current.proj.EQ_RAD, lat);
2115 return;
2116 }
2117
2118 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2119 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
2120 sincosd (lat, &slat, &clat);
2121 sincosd (0.5 * lon, &slon, &clon);
2122
2123 D = GMT->current.proj.EQ_RAD * sqrt (2.0 / (1.0 + clat * clon));
2124 *x = 2.0 * D * clat * slon;
2125 *y = D * slat;
2126 }
2127
gmtproj_ihammer(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2128 GMT_LOCAL void gmtproj_ihammer (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2129 /* Convert Hammer-Aitoff Equal-Area x/y to lon/lat */
2130 double rho, a, sin_c, cos_c;
2131
2132 x *= 0.5;
2133 rho = hypot (x, y);
2134 a = 0.5 * rho * GMT->current.proj.i_EQ_RAD; /* a = sin(c/2) */
2135 a *= a; /* a = sin(c/2)**2 */
2136 cos_c = 1.0 - 2.0 * a; /* cos_c = cos(c) */
2137 if (cos_c < -GMT_PROJ_CONV_LIMIT) { /* Horizon */
2138 *lat = *lon = GMT->session.d_NaN;
2139 return;
2140 }
2141 sin_c = sqrt (1.0 - a) * GMT->current.proj.i_EQ_RAD; /* sin_c = sin(c)/rho */
2142 *lat = asind (y * sin_c);
2143 *lon = 2.0 * d_atan2d (x * sin_c, cos_c) + GMT->current.proj.central_meridian;
2144 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
2145 }
2146
2147 /* -JV VAN DER GRINTEN PROJECTION */
2148
gmtproj_vgrinten(struct GMT_CTRL * GMT,double lon0,double scale)2149 GMT_LOCAL void gmtproj_vgrinten (struct GMT_CTRL *GMT, double lon0, double scale) {
2150 /* Set up van der Grinten projection; scale is unused for now(?) */
2151 gmt_M_unused(scale);
2152 gmtproj_check_R_J (GMT, &lon0);
2153 GMT->current.proj.central_meridian = lon0;
2154 GMT->current.proj.v_r = M_PI * GMT->current.proj.EQ_RAD;
2155 GMT->current.proj.v_ir = 1.0 / GMT->current.proj.v_r;
2156 }
2157
gmtproj_grinten(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2158 GMT_LOCAL void gmtproj_grinten (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2159 /* Convert lon/lat to van der Grinten x/y */
2160 double flat, A, A2, G, P, P2, Q, P2A2, i_P2A2, GP2, c, s, theta;
2161
2162 flat = fabs (lat);
2163 if (flat > (90.0 - GMT_PROJ_CONV_LIMIT)) { /* Save time */
2164 *x = 0.0;
2165 *y = M_PI * copysign (GMT->current.proj.EQ_RAD, lat);
2166 return;
2167 }
2168 if (doubleAlmostEqualZero (lon, GMT->current.proj.central_meridian)) { /* Save time */
2169 theta = d_asin (2.0 * fabs (lat) / 180.0);
2170 *x = 0.0;
2171 *y = M_PI * copysign (GMT->current.proj.EQ_RAD, lat) * tan (0.5 * theta);
2172 return;
2173 }
2174
2175 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2176
2177 if (gmt_M_proj_is_zero (flat)) { /* Save time */
2178 *x = GMT->current.proj.EQ_RAD * D2R * lon;
2179 *y = 0.0;
2180 return;
2181 }
2182
2183 theta = d_asin (2.0 * fabs (lat) / 180.0);
2184
2185 A = 0.5 * fabs (180.0 / lon - lon / 180.0);
2186 A2 = A * A;
2187 sincos (theta, &s, &c);
2188 G = c / (s + c - 1.0);
2189 P = G * (2.0 / s - 1.0);
2190 Q = A2 + G;
2191 P2 = P * P;
2192 P2A2 = P2 + A2;
2193 GP2 = G - P2;
2194 i_P2A2 = 1.0 / P2A2;
2195
2196 *x = copysign (GMT->current.proj.v_r, lon) * (A * GP2 + sqrt (A2 * GP2 * GP2 - P2A2 * (G*G - P2))) * i_P2A2;
2197 *y = copysign (GMT->current.proj.v_r, lat) * (P * Q - A * sqrt ((A2 + 1.0) * P2A2 - Q * Q)) * i_P2A2;
2198 }
2199
gmtproj_igrinten(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2200 GMT_LOCAL void gmtproj_igrinten (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2201 /* Convert van der Grinten x/y to lon/lat */
2202 double x2, c1, x2y2, y2, c2, c3, d, a1, m1, theta1, p;
2203
2204 x *= GMT->current.proj.v_ir;
2205 y *= GMT->current.proj.v_ir;
2206 x2 = x * x; y2 = y * y;
2207 x2y2 = x2 + y2;
2208 if (x2y2 > 1.0) { /* Horizon */
2209 *lat = *lon = GMT->session.d_NaN;
2210 return;
2211 }
2212 c1 = -fabs(y) * (1.0 + x2y2);
2213 c2 = c1 - 2 * y2 + x2;
2214 p = x2y2 * x2y2;
2215 c3 = -2.0 * c1 + 1.0 + 2.0 * y2 + p;
2216 d = y2 / c3 + (2 * pow (c2 / c3, 3.0) - 9.0 * c1 * c2 / (c3 * c3)) / 27.0;
2217 a1 = (c1 - c2 * c2 / (3.0 * c3)) / c3;
2218 m1 = 2.0 * sqrt (-a1 / 3.0);
2219 theta1 = d_acos (3.0 * d / (a1 * m1)) / 3.0;
2220
2221 *lat = copysign (180.0, y) * (-m1 * cos (theta1 + M_PI/3.0) - c2 / (3.0 * c3));
2222 *lon = GMT->current.proj.central_meridian;
2223 if (x != 0.0) *lon += 90.0 * (x2y2 - 1.0 + sqrt (1.0 + 2 * (x2 - y2) + p)) / x;
2224 }
2225
2226 /* -JR WINKEL-TRIPEL MODIFIED gmt_M_is_azimuthal PROJECTION */
2227
gmtproj_vwinkel(struct GMT_CTRL * GMT,double lon0,double scale)2228 GMT_LOCAL void gmtproj_vwinkel (struct GMT_CTRL *GMT, double lon0, double scale) {
2229 /* Set up Winkel Tripel projection */
2230
2231 gmtproj_check_R_J (GMT, &lon0);
2232 GMT->current.proj.r_cosphi1 = cosd (50.0+(28.0/60.0));
2233 GMT->current.proj.central_meridian = lon0;
2234 GMT->current.proj.w_r = 0.25 * (scale * GMT->current.proj.M_PR_DEG * 360.0); /* = Half the minor axis */
2235 }
2236
gmtproj_winkel(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2237 GMT_LOCAL void gmtproj_winkel (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2238 /* Convert lon/lat to Winkel Tripel x/y */
2239 double X, D, x1, x2, y1, c, s;
2240
2241 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2242 lat *= D2R;
2243 lon *= (0.5 * D2R);
2244
2245 /* Fist find Aitoff x/y */
2246
2247 sincos (lat, &s, &c);
2248 D = d_acos (c * cos (lon));
2249 if (gmt_M_proj_is_zero (D))
2250 x1 = y1 = 0.0;
2251 else {
2252 X = s / sin (D);
2253 x1 = copysign (D * d_sqrt (1.0 - X * X), lon);
2254 y1 = D * X;
2255 }
2256
2257 /* Then get equirectangular projection */
2258
2259 x2 = lon * GMT->current.proj.r_cosphi1;
2260 /* y2 = lat; use directly in averaging below */
2261
2262 /* Winkler is the average value */
2263
2264 *x = GMT->current.proj.EQ_RAD * (x1 + x2);
2265 *y = 0.5 * GMT->current.proj.EQ_RAD * (y1 + lat);
2266 }
2267
gmtproj_iwinkel(struct GMT_CTRL * P,double * lon,double * lat,double x,double y)2268 GMT_LOCAL void gmtproj_iwinkel (struct GMT_CTRL *P, double *lon, double *lat, double x, double y) {
2269 /* Convert Winkel Tripel x/y to lon/lat */
2270 /* Based on iterative solution published by:
2271 * Ipbuker, 2002, Cartography & Geographical Information Science, 20, 1, 37-42.
2272 */
2273
2274 int n_iter = 0;
2275 double phi0, lambda0, sp, cp, s2p, sl, cl, sl2, cl2, C, D, sq_C, C_32;
2276 double f1, f2, df1dp, df1dl, df2dp, df2dl, denom, delta;
2277
2278 x *= P->current.proj.i_EQ_RAD;
2279 y *= P->current.proj.i_EQ_RAD;
2280 *lat = y / M_PI; /* Initial guesses for lon and lat */
2281 *lon = x / M_PI;
2282 if (fabs (y) < GMT_PROJ_CONV_LIMIT && fabs (x) < GMT_PROJ_CONV_LIMIT) { /* At ~origin, C is ~zero so no division */
2283 *lon *= R2D;
2284 *lon += P->current.proj.central_meridian;
2285 return;
2286 }
2287 do {
2288 phi0 = *lat;
2289 lambda0 = *lon;
2290 sincos (phi0, &sp, &cp);
2291 sincos (lambda0, &sl, &cl);
2292 sincos (0.5 * lambda0, &sl2, &cl2);
2293 s2p = sin (2.0 * phi0);
2294 D = acos (cp * cl2);
2295 C = 1.0 - cp * cp * cl2 * cl2;
2296 sq_C = sqrt (C);
2297 C_32 = C * sq_C;
2298 f1 = 0.5 * ((2.0 * D * cp * sl2) / sq_C + lambda0 * P->current.proj.r_cosphi1) - x;
2299 f2 = 0.5 * (D * sp / sq_C + phi0) - y;
2300 df1dp = (0.25 * (sl * s2p) / C) - ((D * sp * sl2) / C_32);
2301 df1dl = 0.5 * ((cp * cp * sl2 * sl2) / C + (D * cp * cl2 * sp * sp) / C_32 + P->current.proj.r_cosphi1);
2302 df2dp = 0.5 * ((sp * sp * cl2) / C + (D * (1.0 - cl2 * cl2) * cp) / C_32 + 1.0);
2303 df2dl = 0.125 * ((s2p * sl2) / C - (D * sp * cp * cp * sl) / C_32);
2304 denom = df1dp * df2dl - df2dp * df1dl;
2305 *lat = phi0 - (f1 * df2dl - f2 * df1dl) / denom;
2306 *lon = lambda0 - (f2 * df1dp - f1 * df2dp) / denom;
2307 delta = fabs (*lat - phi0) + fabs (*lon - lambda0);
2308 n_iter++;
2309 }
2310 while (delta > 1e-12 && n_iter < GMT_PROJ_MAX_ITERATIONS);
2311 *lon *= R2D;
2312 if (fabs (*lon) > 180.0) { /* Horizon */
2313 *lat = *lon = P->session.d_NaN;
2314 return;
2315 }
2316 *lon += P->current.proj.central_meridian;
2317 *lat *= R2D;
2318 }
2319
gmtproj_left_winkel(struct GMT_CTRL * GMT,double y)2320 double gmtproj_left_winkel (struct GMT_CTRL *GMT, double y) {
2321 double c, phi, x;
2322
2323 y -= GMT->current.proj.origin[GMT_Y];
2324 y *= GMT->current.proj.i_scale[GMT_Y];
2325 gmtproj_iwinkel_sub (GMT, y, &phi);
2326 gmt_geo_to_xy (GMT, GMT->current.proj.central_meridian - 180.0, phi, &x, &c);
2327 return (x);
2328 }
2329
gmtproj_right_winkel(struct GMT_CTRL * GMT,double y)2330 double gmtproj_right_winkel (struct GMT_CTRL *GMT, double y) {
2331 double c, phi, x;
2332
2333 y -= GMT->current.proj.origin[GMT_Y];
2334 y *= GMT->current.proj.i_scale[GMT_Y];
2335 gmtproj_iwinkel_sub (GMT, y, &phi);
2336 gmt_geo_to_xy (GMT, GMT->current.proj.central_meridian + 180.0, phi, &x, &c);
2337 return (x);
2338 }
2339
2340 /* -JKf ECKERT IV PROJECTION */
2341
gmtproj_veckert4(struct GMT_CTRL * GMT,double lon0)2342 GMT_LOCAL void gmtproj_veckert4 (struct GMT_CTRL *GMT, double lon0) {
2343 /* Set up Eckert IV projection */
2344
2345 gmtproj_check_R_J (GMT, &lon0);
2346 GMT->current.proj.k4_x = 2.0 * GMT->current.proj.EQ_RAD / d_sqrt (M_PI * (4.0 + M_PI));
2347 GMT->current.proj.k4_y = 2.0 * GMT->current.proj.EQ_RAD * d_sqrt (M_PI / (4.0 + M_PI));
2348 GMT->current.proj.k4_ix = 1.0 / GMT->current.proj.k4_x;
2349 GMT->current.proj.k4_iy = 1.0 / GMT->current.proj.k4_y;
2350 GMT->current.proj.central_meridian = lon0;
2351 }
2352
gmtproj_eckert4(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2353 GMT_LOCAL void gmtproj_eckert4 (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2354 /* Convert lon/lat to Eckert IV x/y */
2355 int n_iter = 0;
2356 double phi, delta, s_lat, s, c;
2357
2358 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2359 phi = lat * D2R;
2360 s_lat = sin (phi);
2361 phi *= 0.5;
2362 do {
2363 sincos (phi, &s, &c);
2364 delta = -(phi + s * c + 2.0 * s - (2.0 + M_PI_2) * s_lat) / (2.0 * c * (1.0 + c));
2365 phi += delta;
2366 n_iter++;
2367 }
2368 while (fabs(delta) > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
2369
2370 sincos (phi, &s, &c);
2371 *x = GMT->current.proj.k4_x * lon * D2R * (1.0 + c);
2372 *y = GMT->current.proj.k4_y * s;
2373 }
2374
gmtproj_ieckert4(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2375 GMT_LOCAL void gmtproj_ieckert4 (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2376 /* Convert Eckert IV x/y to lon/lat */
2377
2378 double phi, s, c;
2379
2380 s = y * GMT->current.proj.k4_iy;
2381 phi = d_asin (s);
2382 c = cos (phi);
2383 *lon = R2D * x * GMT->current.proj.k4_ix / (1.0 + c);
2384 if (fabs (*lon) > 180.0) { /* Horizon */
2385 *lat = *lon = GMT->session.d_NaN;
2386 return;
2387 }
2388 *lon += GMT->current.proj.central_meridian;
2389 *lat = d_asind ((phi + s * c + 2.0 * s) / (2.0 + M_PI_2));
2390 }
2391
gmtproj_left_eckert4(struct GMT_CTRL * GMT,double y)2392 double gmtproj_left_eckert4 (struct GMT_CTRL *GMT, double y) {
2393 double x, phi;
2394
2395 y -= GMT->current.proj.origin[GMT_Y];
2396 y *= GMT->current.proj.i_scale[GMT_Y];
2397 phi = d_asin (y * GMT->current.proj.k4_iy);
2398 x = GMT->current.proj.k4_x * D2R * (GMT->common.R.wesn[XLO] - GMT->current.proj.central_meridian) * (1.0 + cos (phi));
2399 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2400 }
2401
gmtproj_right_eckert4(struct GMT_CTRL * GMT,double y)2402 double gmtproj_right_eckert4 (struct GMT_CTRL *GMT, double y) {
2403 double x, phi;
2404
2405 y -= GMT->current.proj.origin[GMT_Y];
2406 y *= GMT->current.proj.i_scale[GMT_Y];
2407 phi = d_asin (y * GMT->current.proj.k4_iy);
2408 x = GMT->current.proj.k4_x * D2R * (GMT->common.R.wesn[XHI] - GMT->current.proj.central_meridian) * (1.0 + cos (phi));
2409 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2410 }
2411
2412 /* -JKs ECKERT VI PROJECTION */
2413
gmtproj_veckert6(struct GMT_CTRL * GMT,double lon0)2414 GMT_LOCAL void gmtproj_veckert6 (struct GMT_CTRL *GMT, double lon0) {
2415 /* Set up Eckert VI projection */
2416
2417 gmtproj_check_R_J (GMT, &lon0);
2418 GMT->current.proj.k6_r = GMT->current.proj.EQ_RAD / sqrt (2.0 + M_PI);
2419 GMT->current.proj.k6_ir = 1.0 / GMT->current.proj.k6_r;
2420 GMT->current.proj.central_meridian = lon0;
2421 }
2422
gmtproj_eckert6(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2423 GMT_LOCAL void gmtproj_eckert6 (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2424 /* Convert lon/lat to Eckert VI x/y */
2425 int n_iter = 0;
2426 double phi, delta, s_lat, s, c;
2427
2428 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2429 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
2430 phi = lat * D2R;
2431 s_lat = sin (phi);
2432 do {
2433 sincos (phi, &s, &c);
2434 delta = -(phi + s - (1.0 + M_PI_2) * s_lat) / (1.0 + c);
2435 phi += delta;
2436 n_iter++;
2437 }
2438 while (fabs(delta) > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
2439
2440 *x = GMT->current.proj.k6_r * lon * D2R * (1.0 + cos (phi));
2441 *y = 2.0 * GMT->current.proj.k6_r * phi;
2442 }
2443
gmtproj_ieckert6(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2444 GMT_LOCAL void gmtproj_ieckert6 (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2445 /* Convert Eckert VI x/y to lon/lat */
2446
2447 double phi, c, s;
2448
2449 phi = 0.5 * y * GMT->current.proj.k6_ir;
2450 sincos (phi, &s, &c);
2451 *lon = R2D * x * GMT->current.proj.k6_ir / (1.0 + c);
2452 if (fabs (*lon) > 180.0) { /* Horizon */
2453 *lat = *lon = GMT->session.d_NaN;
2454 return;
2455 }
2456 *lon += GMT->current.proj.central_meridian;
2457 *lat = d_asin ((phi + s) / (1.0 + M_PI_2)) * R2D;
2458 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
2459 }
2460
gmtproj_left_eckert6(struct GMT_CTRL * GMT,double y)2461 double gmtproj_left_eckert6 (struct GMT_CTRL *GMT, double y) {
2462 double x, phi;
2463
2464 y -= GMT->current.proj.origin[GMT_Y];
2465 y *= GMT->current.proj.i_scale[GMT_Y];
2466 phi = 0.5 * y * GMT->current.proj.k6_ir;
2467 x = GMT->current.proj.k6_r * D2R * (GMT->common.R.wesn[XLO] - GMT->current.proj.central_meridian) * (1.0 + cos (phi));
2468 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2469 }
2470
gmtproj_right_eckert6(struct GMT_CTRL * GMT,double y)2471 double gmtproj_right_eckert6 (struct GMT_CTRL *GMT, double y) {
2472 double x, phi;
2473
2474 y -= GMT->current.proj.origin[GMT_Y];
2475 y *= GMT->current.proj.i_scale[GMT_Y];
2476 phi = 0.5 * y * GMT->current.proj.k6_ir;
2477 x = GMT->current.proj.k6_r * D2R * (GMT->common.R.wesn[XHI] - GMT->current.proj.central_meridian) * (1.0 + cos (phi));
2478 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2479 }
2480
2481 /* -JN ROBINSON PSEUDOCYLINDRICAL PROJECTION */
2482
gmtproj_vrobinson(struct GMT_CTRL * GMT,double lon0)2483 GMT_LOCAL void gmtproj_vrobinson (struct GMT_CTRL *GMT, double lon0) {
2484 /* Set up Robinson projection */
2485 int err_flag = 0;
2486
2487 if (GMT->current.setting.interpolant == GMT_SPLINE_LINEAR) { /* Must reset and warn */
2488 GMT_Report (GMT->parent, GMT_MSG_WARNING, "-JN requires Akima or Cubic spline interpolant, set to Akima\n");
2489 GMT->current.setting.interpolant = GMT_SPLINE_AKIMA;
2490 }
2491
2492 gmtproj_check_R_J (GMT, &lon0);
2493 GMT->current.proj.n_cx = 0.8487 * GMT->current.proj.EQ_RAD * D2R;
2494 GMT->current.proj.n_cy = 1.3523 * GMT->current.proj.EQ_RAD;
2495 GMT->current.proj.n_i_cy = 1.0 / GMT->current.proj.n_cy;
2496 GMT->current.proj.central_meridian = lon0;
2497
2498 GMT->current.proj.n_phi[0] = 0; GMT->current.proj.n_X[0] = 1.0000; GMT->current.proj.n_Y[0] = 0.0000;
2499 GMT->current.proj.n_phi[1] = 5; GMT->current.proj.n_X[1] = 0.9986; GMT->current.proj.n_Y[1] = 0.0620;
2500 GMT->current.proj.n_phi[2] = 10; GMT->current.proj.n_X[2] = 0.9954; GMT->current.proj.n_Y[2] = 0.1240;
2501 GMT->current.proj.n_phi[3] = 15; GMT->current.proj.n_X[3] = 0.9900; GMT->current.proj.n_Y[3] = 0.1860;
2502 GMT->current.proj.n_phi[4] = 20; GMT->current.proj.n_X[4] = 0.9822; GMT->current.proj.n_Y[4] = 0.2480;
2503 GMT->current.proj.n_phi[5] = 25; GMT->current.proj.n_X[5] = 0.9730; GMT->current.proj.n_Y[5] = 0.3100;
2504 GMT->current.proj.n_phi[6] = 30; GMT->current.proj.n_X[6] = 0.9600; GMT->current.proj.n_Y[6] = 0.3720;
2505 GMT->current.proj.n_phi[7] = 35; GMT->current.proj.n_X[7] = 0.9427; GMT->current.proj.n_Y[7] = 0.4340;
2506 GMT->current.proj.n_phi[8] = 40; GMT->current.proj.n_X[8] = 0.9216; GMT->current.proj.n_Y[8] = 0.4958;
2507 GMT->current.proj.n_phi[9] = 45; GMT->current.proj.n_X[9] = 0.8962; GMT->current.proj.n_Y[9] = 0.5571;
2508 GMT->current.proj.n_phi[10] = 50; GMT->current.proj.n_X[10] = 0.8679; GMT->current.proj.n_Y[10] = 0.6176;
2509 GMT->current.proj.n_phi[11] = 55; GMT->current.proj.n_X[11] = 0.8350; GMT->current.proj.n_Y[11] = 0.6769;
2510 GMT->current.proj.n_phi[12] = 60; GMT->current.proj.n_X[12] = 0.7986; GMT->current.proj.n_Y[12] = 0.7346;
2511 GMT->current.proj.n_phi[13] = 65; GMT->current.proj.n_X[13] = 0.7597; GMT->current.proj.n_Y[13] = 0.7903;
2512 GMT->current.proj.n_phi[14] = 70; GMT->current.proj.n_X[14] = 0.7186; GMT->current.proj.n_Y[14] = 0.8435;
2513 GMT->current.proj.n_phi[15] = 75; GMT->current.proj.n_X[15] = 0.6732; GMT->current.proj.n_Y[15] = 0.8936;
2514 GMT->current.proj.n_phi[16] = 80; GMT->current.proj.n_X[16] = 0.6213; GMT->current.proj.n_Y[16] = 0.9394;
2515 GMT->current.proj.n_phi[17] = 85; GMT->current.proj.n_X[17] = 0.5722; GMT->current.proj.n_Y[17] = 0.9761;
2516 GMT->current.proj.n_phi[18] = 90; GMT->current.proj.n_X[18] = 0.5322; GMT->current.proj.n_Y[18] = 1.0000;
2517 if (GMT->current.setting.interpolant == GMT_SPLINE_CUBIC) { /* Natural cubic spline */
2518 err_flag = gmtlib_cspline (GMT, GMT->current.proj.n_phi, GMT->current.proj.n_X, GMT_N_ROBINSON, GMT->current.proj.n_x_coeff);
2519 err_flag += gmtlib_cspline (GMT, GMT->current.proj.n_phi, GMT->current.proj.n_Y, GMT_N_ROBINSON, GMT->current.proj.n_y_coeff);
2520 err_flag += gmtlib_cspline (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_X, GMT_N_ROBINSON, GMT->current.proj.n_yx_coeff);
2521 err_flag += gmtlib_cspline (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_phi, GMT_N_ROBINSON, GMT->current.proj.n_iy_coeff);
2522 }
2523 else { /* Akimas spline */
2524 err_flag = gmtlib_akima (GMT, GMT->current.proj.n_phi, GMT->current.proj.n_X, GMT_N_ROBINSON, GMT->current.proj.n_x_coeff);
2525 err_flag += gmtlib_akima (GMT, GMT->current.proj.n_phi, GMT->current.proj.n_Y, GMT_N_ROBINSON, GMT->current.proj.n_y_coeff);
2526 err_flag += gmtlib_akima (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_X, GMT_N_ROBINSON, GMT->current.proj.n_yx_coeff);
2527 err_flag += gmtlib_akima (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_phi, GMT_N_ROBINSON, GMT->current.proj.n_iy_coeff);
2528 }
2529 if (err_flag) GMT_Report (GMT->parent, GMT_MSG_ERROR, "Interpolation failed in gmtproj_vrobinson?\n");
2530 }
2531
gmtproj_robinson(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2532 GMT_LOCAL void gmtproj_robinson (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2533 /* Convert lon/lat to Robinson x/y */
2534 double phi, X, Y;
2535
2536 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2537 phi = fabs (lat);
2538
2539 X = gmtproj_robinson_spline (GMT, phi, GMT->current.proj.n_phi, GMT->current.proj.n_X, GMT->current.proj.n_x_coeff);
2540 Y = gmtproj_robinson_spline (GMT, phi, GMT->current.proj.n_phi, GMT->current.proj.n_Y, GMT->current.proj.n_y_coeff);
2541 *x = GMT->current.proj.n_cx * X * lon; /* D2R is in n_cx already */
2542 *y = GMT->current.proj.n_cy * copysign (Y, lat);
2543 }
2544
gmtproj_irobinson(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2545 GMT_LOCAL void gmtproj_irobinson (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2546 /* Convert Robinson x/y to lon/lat */
2547 double X, Y;
2548
2549 Y = fabs (y * GMT->current.proj.n_i_cy);
2550 *lat = gmtproj_robinson_spline (GMT, Y, GMT->current.proj.n_Y, GMT->current.proj.n_phi, GMT->current.proj.n_iy_coeff);
2551 X = gmtproj_robinson_spline (GMT, *lat, GMT->current.proj.n_phi, GMT->current.proj.n_X, GMT->current.proj.n_x_coeff);
2552 *lon = x / (GMT->current.proj.n_cx * X);
2553 if ((fabs (*lon) - GMT_PROJ_CONV_LIMIT) > 180.0) { /* Horizon */
2554 *lat = *lon = GMT->session.d_NaN;
2555 return;
2556 }
2557 *lon += GMT->current.proj.central_meridian;
2558 if (y < 0.0) *lat = -(*lat);
2559 }
2560
gmtproj_left_robinson(struct GMT_CTRL * GMT,double y)2561 double gmtproj_left_robinson (struct GMT_CTRL *GMT, double y) {
2562 double x, X, Y;
2563
2564 y -= GMT->current.proj.origin[GMT_Y];
2565 y *= GMT->current.proj.i_scale[GMT_Y];
2566 Y = fabs (y * GMT->current.proj.n_i_cy);
2567 X = gmtproj_robinson_spline (GMT, Y, GMT->current.proj.n_Y, GMT->current.proj.n_X, GMT->current.proj.n_yx_coeff);
2568 x = GMT->current.proj.n_cx * X * (GMT->common.R.wesn[XLO] - GMT->current.proj.central_meridian);
2569 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2570 }
2571
gmtproj_right_robinson(struct GMT_CTRL * GMT,double y)2572 double gmtproj_right_robinson (struct GMT_CTRL *GMT, double y) {
2573 double x, X, Y;
2574
2575 y -= GMT->current.proj.origin[GMT_Y];
2576 y *= GMT->current.proj.i_scale[GMT_Y];
2577 Y = fabs (y * GMT->current.proj.n_i_cy);
2578 X = gmtproj_robinson_spline (GMT, Y, GMT->current.proj.n_Y, GMT->current.proj.n_X, GMT->current.proj.n_yx_coeff);
2579 x = GMT->current.proj.n_cx * X * (GMT->common.R.wesn[XHI] - GMT->current.proj.central_meridian);
2580 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2581 }
2582
2583 #if 0
2584 double gmtproj_left_robinson (struct GMT_CTRL *GMT, double y) {
2585 double x, X, Y;
2586
2587 y -= GMT->current.proj.origin[GMT_Y];
2588 y *= GMT->current.proj.i_scale[GMT_Y];
2589 Y = fabs (y * GMT->current.proj.n_i_cy);
2590 if (gmt_intpol (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_X, NULL, GMT_N_ROBINSON, 1, &Y, &X, 0.0, GMT->current.setting.interpolant)) {
2591 GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT Internal error in gmtproj_left_robinson!\n");
2592 return GMT->session.d_NaN;
2593 }
2594
2595 x = GMT->current.proj.n_cx * X * (GMT->common.R.wesn[XLO] - GMT->current.proj.central_meridian);
2596 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2597 }
2598
2599 double gmtproj_right_robinson (struct GMT_CTRL *GMT, double y) {
2600 double x, X, Y;
2601
2602 y -= GMT->current.proj.origin[GMT_Y];
2603 y *= GMT->current.proj.i_scale[GMT_Y];
2604 Y = fabs (y * GMT->current.proj.n_i_cy);
2605 if (gmt_intpol (GMT, GMT->current.proj.n_Y, GMT->current.proj.n_X, NULL, GMT_N_ROBINSON, 1, &Y, &X, 0.0, GMT->current.setting.interpolant)) {
2606 GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT Internal error in gmtproj_right_robinson!\n");
2607 return GMT->session.d_NaN;
2608 }
2609
2610 x = GMT->current.proj.n_cx * X * (GMT->common.R.wesn[XHI] - GMT->current.proj.central_meridian);
2611 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2612 }
2613 #endif
2614
2615 /* -JI SINUSOIDAL EQUAL AREA PROJECTION */
2616
gmtproj_vsinusoidal(struct GMT_CTRL * GMT,double lon0)2617 GMT_LOCAL void gmtproj_vsinusoidal (struct GMT_CTRL *GMT, double lon0) {
2618 /* Set up Sinusoidal projection */
2619
2620 gmtproj_check_R_J (GMT, &lon0);
2621 GMT->current.proj.central_meridian = lon0;
2622 }
2623
gmtproj_sinusoidal(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2624 GMT_LOCAL void gmtproj_sinusoidal (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2625 /* Convert lon/lat to Sinusoidal Equal-Area x/y */
2626
2627 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2628 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
2629
2630 lat *= D2R;
2631
2632 *x = GMT->current.proj.EQ_RAD * lon * D2R * cos (lat);
2633 *y = GMT->current.proj.EQ_RAD * lat;
2634 }
2635
gmtproj_isinusoidal(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2636 GMT_LOCAL void gmtproj_isinusoidal (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2637 /* Convert Sinusoidal Equal-Area x/y to lon/lat */
2638
2639 *lat = y * GMT->current.proj.i_EQ_RAD;
2640 *lon = (doubleAlmostEqual (fabs (*lat), M_PI)) ? 0.0 : R2D * x / (GMT->current.proj.EQ_RAD * cos (*lat));
2641 if (fabs (*lon) > 180.0) { /* Horizon */
2642 *lat = *lon = GMT->session.d_NaN;
2643 return;
2644 }
2645 *lon += GMT->current.proj.central_meridian;
2646 *lat *= R2D;
2647 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
2648 }
2649
gmtproj_left_sinusoidal(struct GMT_CTRL * GMT,double y)2650 double gmtproj_left_sinusoidal (struct GMT_CTRL *GMT, double y) {
2651 double x, lat;
2652 y -= GMT->current.proj.origin[GMT_Y];
2653 y *= GMT->current.proj.i_scale[GMT_Y];
2654 lat = y * GMT->current.proj.i_EQ_RAD;
2655 x = GMT->current.proj.EQ_RAD * (GMT->common.R.wesn[XLO] - GMT->current.proj.central_meridian) * D2R * cos (lat);
2656 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2657 }
2658
gmtproj_right_sinusoidal(struct GMT_CTRL * GMT,double y)2659 double gmtproj_right_sinusoidal (struct GMT_CTRL *GMT, double y) {
2660 double x, lat;
2661 y -= GMT->current.proj.origin[GMT_Y];
2662 y *= GMT->current.proj.i_scale[GMT_Y];
2663 lat = y * GMT->current.proj.i_EQ_RAD;
2664 x = GMT->current.proj.EQ_RAD * (GMT->common.R.wesn[XHI] - GMT->current.proj.central_meridian) * D2R * cos (lat);
2665 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
2666 }
2667
2668 /* -JC CASSINI PROJECTION */
2669
gmtproj_vcassini(struct GMT_CTRL * GMT,double lon0,double lat0)2670 GMT_LOCAL void gmtproj_vcassini (struct GMT_CTRL *GMT, double lon0, double lat0) {
2671 /* Set up Cassini projection */
2672
2673 double e1, s2, c2;
2674
2675 gmtproj_check_R_J (GMT, &lon0);
2676 GMT->current.proj.central_meridian = lon0;
2677 GMT->current.proj.pole = lat0;
2678 GMT->current.proj.c_p = lat0 * D2R;
2679 sincos (2.0 * GMT->current.proj.c_p, &s2, &c2);
2680
2681 e1 = (1.0 - d_sqrt (GMT->current.proj.one_m_ECC2)) / (1.0 + d_sqrt (GMT->current.proj.one_m_ECC2));
2682
2683 GMT->current.proj.c_c1 = 1.0 - (1.0/4.0) * GMT->current.proj.ECC2 - (3.0/64.0) * GMT->current.proj.ECC4 - (5.0/256.0) * GMT->current.proj.ECC6;
2684 GMT->current.proj.c_c2 = -((3.0/8.0) * GMT->current.proj.ECC2 + (3.0/32.0) * GMT->current.proj.ECC4 + (25.0/768.0) * GMT->current.proj.ECC6);
2685 GMT->current.proj.c_c3 = (15.0/128.0) * GMT->current.proj.ECC4 + (45.0/512.0) * GMT->current.proj.ECC6;
2686 GMT->current.proj.c_c4 = -(35.0/768.0) * GMT->current.proj.ECC6;
2687 GMT->current.proj.c_M0 = GMT->current.proj.EQ_RAD * (GMT->current.proj.c_c1 * GMT->current.proj.c_p + s2 * (GMT->current.proj.c_c2 + c2 * (GMT->current.proj.c_c3 + c2 * GMT->current.proj.c_c4)));
2688 GMT->current.proj.c_i1 = 1.0 / (GMT->current.proj.EQ_RAD * GMT->current.proj.c_c1);
2689 GMT->current.proj.c_i2 = (3.0/2.0) * e1 - (29.0/12.0) * pow (e1, 3.0);
2690 GMT->current.proj.c_i3 = (21.0/8.0) * e1 * e1 - (1537.0/128.0) * pow (e1, 4.0);
2691 GMT->current.proj.c_i4 = (151.0/24.0) * pow (e1, 3.0);
2692 GMT->current.proj.c_i5 = (1097.0/64.0) * pow (e1, 4.0);
2693 }
2694
gmtproj_cassini(struct GMT_CTRL * P,double lon,double lat,double * x,double * y)2695 GMT_LOCAL void gmtproj_cassini (struct GMT_CTRL *P, double lon, double lat, double *x, double *y) {
2696 /* Convert lon/lat to Cassini x/y */
2697
2698 double tany, N, T, A, C, M, s, c, s2, c2, A2, A3;
2699
2700 gmt_M_wind_lon (P, lon) /* Remove central meridian and place lon in -180/+180 range */
2701 lon *= D2R;
2702
2703 if (gmt_M_proj_is_zero (lat)) { /* Quick when lat is zero */
2704 *x = P->current.proj.EQ_RAD * lon;
2705 *y = -P->current.proj.c_M0;
2706 return;
2707 }
2708
2709 lat *= D2R;
2710 sincos (lat, &s, &c);
2711 sincos (2.0 * lat, &s2, &c2);
2712 tany = s / c;
2713 N = P->current.proj.EQ_RAD / sqrt (1.0 - P->current.proj.ECC2 * s * s);
2714 T = tany * tany;
2715 A = lon * c;
2716 A2 = A * A;
2717 A3 = A2 * A;
2718 C = P->current.proj.ECC2 * c * c * P->current.proj.i_one_m_ECC2;
2719 M = P->current.proj.EQ_RAD * (P->current.proj.c_c1 * lat + s2 * (P->current.proj.c_c2 + c2 * (P->current.proj.c_c3 + c2 * P->current.proj.c_c4)));
2720 *x = N * (A - T * A3 / 6.0 - (8.0 - T + 8 * C) * T * A3 * A2 / 120.0);
2721 *y = M - P->current.proj.c_M0 + N * tany * (0.5 * A2 + (5.0 - T + 6.0 * C) * A2 * A2 / 24.0);
2722 }
2723
gmtproj_icassini(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2724 GMT_LOCAL void gmtproj_icassini (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2725 /* Convert Cassini x/y to lon/lat */
2726
2727 double M1, u1, u2, s, c, phi1, tany, T1, N1, R_1, D, S2, D2, D3;
2728
2729 M1 = GMT->current.proj.c_M0 + y;
2730 u1 = M1 * GMT->current.proj.c_i1;
2731 u2 = 2.0 * u1;
2732 sincos (u2, &s, &c);
2733 phi1 = u1 + s * (GMT->current.proj.c_i2 + c * (GMT->current.proj.c_i3 + c * (GMT->current.proj.c_i4 + c * GMT->current.proj.c_i5)));
2734 if (doubleAlmostEqual (fabs (phi1), M_PI_2)) {
2735 *lat = copysign (M_PI_2, phi1);
2736 *lon = GMT->current.proj.central_meridian;
2737 }
2738 else {
2739 sincos (phi1, &s, &c);
2740 tany = s / c;
2741 T1 = tany * tany;
2742 S2 = 1.0 - GMT->current.proj.ECC2 * s * s;
2743 N1 = GMT->current.proj.EQ_RAD / sqrt (S2);
2744 R_1 = GMT->current.proj.EQ_RAD * GMT->current.proj.one_m_ECC2 / pow (S2, 1.5);
2745 D = x / N1;
2746 D2 = D * D;
2747 D3 = D2 * D;
2748 *lat = R2D * (phi1 - (N1 * tany / R_1) * (0.5 * D2 - (1.0 + 3.0 * T1) * D2 * D2 / 24.0));
2749 *lon = GMT->current.proj.central_meridian + R2D * (D - T1 * D3 / 3.0 + (1.0 + 3.0 * T1) * T1 * D3 * D2 / 15.0) / c;
2750 }
2751 }
2752
2753 /* Cassini functions for the Sphere */
2754
gmtproj_cassini_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2755 GMT_LOCAL void gmtproj_cassini_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2756 /* Convert lon/lat to Cassini x/y for spherical earth */
2757
2758 double slon, clon, clat, tlat, slat;
2759
2760 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2761
2762 if (gmt_M_proj_is_zero (lat)) { /* Quick when lat is zero */
2763 *x = GMT->current.proj.EQ_RAD * lon * D2R;
2764 *y = -GMT->current.proj.EQ_RAD * GMT->current.proj.c_p;
2765 return;
2766 }
2767
2768 sincosd (lon, &slon, &clon);
2769 sincosd (lat, &slat, &clat);
2770 tlat = slat / clat;
2771 *x = GMT->current.proj.EQ_RAD * d_asin (clat * slon);
2772 *y = GMT->current.proj.EQ_RAD * (atan (tlat / clon) - GMT->current.proj.c_p);
2773 }
2774
2775
gmtproj_icassini_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2776 GMT_LOCAL void gmtproj_icassini_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2777 /* Convert Cassini x/y to lon/lat */
2778
2779 double D, sD, cD, cx, tx, sx;
2780
2781 x *= GMT->current.proj.i_EQ_RAD;
2782 D = y * GMT->current.proj.i_EQ_RAD + GMT->current.proj.c_p;
2783 sincos (D, &sD, &cD);
2784 sincos (x, &sx, &cx);
2785 tx = sx / cx;
2786 *lat = d_asind (sD * cx);
2787 *lon = GMT->current.proj.central_meridian + atand (tx / cD);
2788 }
2789
2790 /* -JB ALBERS EQUAL-AREA CONIC PROJECTION */
2791
gmtproj_valbers(struct GMT_CTRL * GMT,double lon0,double lat0,double ph1,double ph2)2792 GMT_LOCAL void gmtproj_valbers (struct GMT_CTRL *GMT, double lon0, double lat0, double ph1, double ph2) {
2793 /* Set up Albers projection */
2794
2795 double s0, s1, s2, c1, c2, q0, q1, q2, m1, m2;
2796
2797 gmtproj_check_R_J (GMT, &lon0);
2798 GMT->current.proj.central_meridian = lon0;
2799 GMT->current.proj.north_pole = (GMT->common.R.wesn[YHI] > 0.0 && (GMT->common.R.wesn[YLO] >= 0.0 || (-GMT->common.R.wesn[YLO]) < GMT->common.R.wesn[YHI]));
2800 GMT->current.proj.pole = (GMT->current.proj.north_pole) ? 90.0 : -90.0;
2801
2802 s0 = sind (lat0);
2803 sincosd (ph1, &s1, &c1);
2804 sincosd (ph2, &s2, &c2);
2805
2806 m1 = c1 * c1 / (1.0 - GMT->current.proj.ECC2 * s1 * s1); /* Actually m1 and m2 squared */
2807 m2 = c2 * c2 / (1.0 - GMT->current.proj.ECC2 * s2 * s2);
2808 q0 = (gmt_M_proj_is_zero (GMT->current.proj.ECC)) ? 2.0 * s0 : GMT->current.proj.one_m_ECC2 * (s0 / (1.0 - GMT->current.proj.ECC2 * s0 * s0) - GMT->current.proj.i_half_ECC * log ((1.0 - GMT->current.proj.ECC * s0) / (1.0 + GMT->current.proj.ECC * s0)));
2809 q1 = (gmt_M_proj_is_zero (GMT->current.proj.ECC)) ? 2.0 * s1 : GMT->current.proj.one_m_ECC2 * (s1 / (1.0 - GMT->current.proj.ECC2 * s1 * s1) - GMT->current.proj.i_half_ECC * log ((1.0 - GMT->current.proj.ECC * s1) / (1.0 + GMT->current.proj.ECC * s1)));
2810 q2 = (gmt_M_proj_is_zero (GMT->current.proj.ECC)) ? 2.0 * s2 : GMT->current.proj.one_m_ECC2 * (s2 / (1.0 - GMT->current.proj.ECC2 * s2 * s2) - GMT->current.proj.i_half_ECC * log ((1.0 - GMT->current.proj.ECC * s2) / (1.0 + GMT->current.proj.ECC * s2)));
2811
2812 GMT->current.proj.a_n = (doubleAlmostEqualZero (ph1, ph2)) ? s1 : (m1 - m2) / (q2 - q1);
2813 GMT->current.proj.a_i_n = 1.0 / GMT->current.proj.a_n;
2814 GMT->current.proj.a_C = m1 + GMT->current.proj.a_n * q1;
2815 GMT->current.proj.a_rho0 = GMT->current.proj.EQ_RAD * sqrt (GMT->current.proj.a_C - GMT->current.proj.a_n * q0) * GMT->current.proj.a_i_n;
2816 GMT->current.proj.a_n2ir2 = (GMT->current.proj.a_n * GMT->current.proj.a_n) / (GMT->current.proj.EQ_RAD * GMT->current.proj.EQ_RAD);
2817 GMT->current.proj.a_test = 1.0 - (GMT->current.proj.i_half_ECC * GMT->current.proj.one_m_ECC2) * log ((1.0 - GMT->current.proj.ECC) / (1.0 + GMT->current.proj.ECC));
2818
2819 }
2820
gmtproj_valbers_sph(struct GMT_CTRL * GMT,double lon0,double lat0,double ph1,double ph2)2821 GMT_LOCAL void gmtproj_valbers_sph (struct GMT_CTRL *GMT, double lon0, double lat0, double ph1, double ph2) {
2822 /* Set up Spherical Albers projection (Snyder, page 100) */
2823
2824 double s1, c1;
2825
2826 gmtproj_check_R_J (GMT, &lon0);
2827 GMT->current.proj.central_meridian = lon0;
2828 GMT->current.proj.north_pole = (GMT->common.R.wesn[YHI] > 0.0 && (GMT->common.R.wesn[YLO] >= 0.0 || (-GMT->common.R.wesn[YLO]) < GMT->common.R.wesn[YHI]));
2829 GMT->current.proj.pole = (GMT->current.proj.north_pole) ? 90.0 : -90.0;
2830
2831 sincosd (ph1, &s1, &c1);
2832
2833 GMT->current.proj.a_n = 0.5 * (s1 + sind (ph2));
2834 GMT->current.proj.a_i_n = 1.0 / GMT->current.proj.a_n;
2835 GMT->current.proj.a_C = c1 * c1 + 2.0 * GMT->current.proj.a_n * s1;
2836 GMT->current.proj.a_rho0 = GMT->current.proj.EQ_RAD * sqrt (GMT->current.proj.a_C - 2.0 * GMT->current.proj.a_n * sind (lat0)) * GMT->current.proj.a_i_n;
2837 GMT->current.proj.a_n2ir2 = 0.5 * GMT->current.proj.a_n / (GMT->current.proj.EQ_RAD * GMT->current.proj.EQ_RAD);
2838 GMT->current.proj.a_Cin = 0.5 * GMT->current.proj.a_C / GMT->current.proj.a_n;
2839 }
2840
gmtproj_albers(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2841 GMT_LOCAL void gmtproj_albers (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2842 /* Convert lon/lat to Albers x/y */
2843
2844 double s, c, q, theta, rho, r;
2845
2846 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2847
2848 s = sind (lat);
2849 if (gmt_M_proj_is_zero (GMT->current.proj.ECC))
2850 q = 2.0 * s;
2851 else {
2852 r = GMT->current.proj.ECC * s;
2853 q = GMT->current.proj.one_m_ECC2 * (s / (1.0 - GMT->current.proj.ECC2 * s * s) - GMT->current.proj.i_half_ECC * log ((1.0 - r) / (1.0 + r)));
2854 }
2855 theta = GMT->current.proj.a_n * lon * D2R;
2856 rho = GMT->current.proj.EQ_RAD * sqrt (GMT->current.proj.a_C - GMT->current.proj.a_n * q) * GMT->current.proj.a_i_n;
2857
2858 sincos (theta, &s, &c);
2859 *x = rho * s;
2860 *y = GMT->current.proj.a_rho0 - rho * c;
2861 }
2862
gmtproj_ialbers(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2863 GMT_LOCAL void gmtproj_ialbers (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2864 /* Convert Albers x/y to lon/lat */
2865
2866 int n_iter;
2867 double theta, rho, q, phi, phi0, s, c, s2, ex_1, delta, r;
2868
2869 theta = (GMT->current.proj.a_n < 0.0) ? d_atan2 (-x, y - GMT->current.proj.a_rho0) : d_atan2 (x, GMT->current.proj.a_rho0 - y);
2870 rho = hypot (x, GMT->current.proj.a_rho0 - y);
2871 q = (GMT->current.proj.a_C - rho * rho * GMT->current.proj.a_n2ir2) * GMT->current.proj.a_i_n;
2872
2873 if (doubleAlmostEqualZero (fabs (q), GMT->current.proj.a_test))
2874 *lat = copysign (90.0, q);
2875 else {
2876 phi = d_asin (0.5 * q);
2877 n_iter = 0;
2878 do {
2879 phi0 = phi;
2880 sincos (phi0, &s, &c);
2881 r = GMT->current.proj.ECC * s;
2882 s2 = s * s;
2883 ex_1 = 1.0 - GMT->current.proj.ECC2 * s2;
2884 phi = phi0 + 0.5 * ex_1 * ex_1 * ((q * GMT->current.proj.i_one_m_ECC2) - s / ex_1
2885 + GMT->current.proj.i_half_ECC * log ((1 - r) / (1.0 + r))) / c;
2886 delta = fabs (phi - phi0);
2887 n_iter++;
2888 }
2889 while (delta > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
2890 *lat = R2D * phi;
2891 }
2892 *lon = GMT->current.proj.central_meridian + R2D * theta * GMT->current.proj.a_i_n;
2893 }
2894
2895 /* Spherical versions of Albers */
2896
gmtproj_albers_sph(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2897 GMT_LOCAL void gmtproj_albers_sph (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2898 /* Convert lon/lat to Spherical Albers x/y */
2899
2900 double s, c, theta, rho;
2901
2902 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2903
2904 if (GMT->current.proj.GMT_convert_latitudes) lat = gmt_M_latg_to_lata (GMT, lat);
2905
2906 theta = GMT->current.proj.a_n * lon * D2R;
2907 rho = GMT->current.proj.EQ_RAD * sqrt (GMT->current.proj.a_C - 2.0 * GMT->current.proj.a_n * sind (lat)) * GMT->current.proj.a_i_n;
2908
2909 sincos (theta, &s, &c);
2910 *x = rho * s;
2911 *y = GMT->current.proj.a_rho0 - rho * c;
2912 }
2913
gmtproj_ialbers_sph(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2914 GMT_LOCAL void gmtproj_ialbers_sph (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2915 /* Convert Spherical Albers x/y to lon/lat */
2916
2917 double theta, A, dy;
2918
2919 theta = (GMT->current.proj.a_n < 0.0) ? d_atan2 (-x, y - GMT->current.proj.a_rho0) : d_atan2 (x, GMT->current.proj.a_rho0 - y);
2920 dy = GMT->current.proj.a_rho0 - y;
2921 A = (x * x + dy * dy) * GMT->current.proj.a_n2ir2;
2922
2923 *lat = d_asind (GMT->current.proj.a_Cin - A);
2924 *lon = GMT->current.proj.central_meridian + R2D * theta * GMT->current.proj.a_i_n;
2925 if (GMT->current.proj.GMT_convert_latitudes) *lat = gmt_M_lata_to_latg (GMT, *lat);
2926 }
2927
2928 /* -JD EQUIDISTANT CONIC PROJECTION */
2929
gmtproj_veconic(struct GMT_CTRL * GMT,double lon0,double lat0,double lat1,double lat2)2930 GMT_LOCAL void gmtproj_veconic (struct GMT_CTRL *GMT, double lon0, double lat0, double lat1, double lat2) {
2931 /* Set up Equidistant Conic projection */
2932 double c1;
2933
2934 gmtproj_check_R_J (GMT, &lon0);
2935 GMT->current.proj.north_pole = (GMT->common.R.wesn[YHI] > 0.0 && (GMT->common.R.wesn[YLO] >= 0.0 || (-GMT->common.R.wesn[YLO]) < GMT->common.R.wesn[YHI]));
2936 c1 = cosd (lat1);
2937 GMT->current.proj.d_n = (doubleAlmostEqualZero (lat1, lat2)) ? sind (lat1) : (c1 - cosd (lat2)) / (D2R * (lat2 - lat1));
2938 GMT->current.proj.d_i_n = R2D / GMT->current.proj.d_n; /* R2D put here instead of in lon for ieconic */
2939 GMT->current.proj.d_G = (c1 / GMT->current.proj.d_n) + lat1 * D2R;
2940 GMT->current.proj.d_rho0 = GMT->current.proj.EQ_RAD * (GMT->current.proj.d_G - lat0 * D2R);
2941 GMT->current.proj.central_meridian = lon0;
2942 }
2943
gmtproj_econic(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2944 GMT_LOCAL void gmtproj_econic (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2945 /* Convert lon/lat to Equidistant Conic x/y */
2946 double rho, theta, s, c;
2947
2948 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2949
2950 rho = GMT->current.proj.EQ_RAD * (GMT->current.proj.d_G - lat * D2R);
2951 theta = GMT->current.proj.d_n * lon * D2R;
2952
2953 sincos (theta, &s, &c);
2954 *x = rho * s;
2955 *y = GMT->current.proj.d_rho0 - rho * c;
2956 }
2957
gmtproj_ieconic(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2958 GMT_LOCAL void gmtproj_ieconic (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
2959 /* Convert Equidistant Conic x/y to lon/lat */
2960
2961 double rho, theta;
2962
2963 rho = hypot (x, GMT->current.proj.d_rho0 - y);
2964 if (GMT->current.proj.d_n < 0.0) rho = -rho;
2965 theta = (GMT->current.proj.d_n < 0.0) ? d_atan2 (-x, y - GMT->current.proj.d_rho0) : d_atan2 (x, GMT->current.proj.d_rho0 - y);
2966 *lat = (GMT->current.proj.d_G - rho * GMT->current.proj.i_EQ_RAD) * R2D;
2967 *lon = GMT->current.proj.central_meridian + theta * GMT->current.proj.d_i_n;
2968 }
2969
2970 /* -JPoly POLYCONIC PROJECTION */
2971
gmtproj_vpolyconic(struct GMT_CTRL * GMT,double lon0,double lat0)2972 GMT_LOCAL void gmtproj_vpolyconic (struct GMT_CTRL *GMT, double lon0, double lat0) {
2973 /* Set up Polyconic projection */
2974
2975 gmtproj_check_R_J (GMT, &lon0);
2976 GMT->current.proj.central_meridian = lon0;
2977 GMT->current.proj.pole = lat0;
2978 }
2979
gmtproj_polyconic(struct GMT_CTRL * GMT,double lon,double lat,double * x,double * y)2980 GMT_LOCAL void gmtproj_polyconic (struct GMT_CTRL *GMT, double lon, double lat, double *x, double *y) {
2981 /* Convert lon/lat to Polyconic x/y */
2982 double sE, cE, sp, cp;
2983
2984 gmt_M_wind_lon (GMT, lon) /* Remove central meridian and place lon in -180/+180 range */
2985
2986 if (gmt_M_proj_is_zero(lat)) {
2987 *x = GMT->current.proj.EQ_RAD * lon * D2R;
2988 *y = GMT->current.proj.EQ_RAD * (lat - GMT->current.proj.pole) * D2R;
2989 }
2990 else {
2991 sincosd(lat, &sp, &cp);
2992 sincosd(lon * sp, &sE, &cE);
2993 cp /= sp; /* = cot(phi) */
2994 *x = GMT->current.proj.EQ_RAD * cp * sE;
2995 *y = GMT->current.proj.EQ_RAD * ((lat - GMT->current.proj.pole) * D2R + cp * (1.0 - cE));
2996 }
2997 }
2998
gmtproj_ipolyconic(struct GMT_CTRL * GMT,double * lon,double * lat,double x,double y)2999 GMT_LOCAL void gmtproj_ipolyconic (struct GMT_CTRL *GMT, double *lon, double *lat, double x, double y) {
3000 /* Convert Polyconic x/y to lon/lat */
3001
3002 double B, phi, phi0, tanp, delta;
3003 int n_iter = 0;
3004
3005 x *= GMT->current.proj.i_EQ_RAD;
3006 y *= GMT->current.proj.i_EQ_RAD;
3007 y += GMT->current.proj.pole * D2R;
3008 if (gmt_M_proj_is_zero (y)) {
3009 *lat = y * R2D + GMT->current.proj.pole;
3010 *lon = x * R2D + GMT->current.proj.central_meridian;
3011 }
3012 else {
3013 B = x * x + y * y;
3014 phi = y;
3015 do {
3016 phi0 = phi;
3017 tanp = tan(phi);
3018 phi -= (y * (phi * tanp + 1.0) - phi - 0.5 * (phi * phi + B) * tanp) /
3019 ((phi - y) / tanp - 1.0);
3020 delta = fabs (phi - phi0);
3021 n_iter++;
3022 }
3023 while (delta > GMT_PROJ_CONV_LIMIT && n_iter < GMT_PROJ_MAX_ITERATIONS);
3024 *lat = R2D * phi;
3025 *lon = GMT->current.proj.central_meridian + asind(x * tanp) / sin(phi);
3026 }
3027 }
3028
gmtproj_left_polyconic(struct GMT_CTRL * GMT,double y)3029 double gmtproj_left_polyconic (struct GMT_CTRL *GMT, double y) {
3030 double x;
3031 y -= GMT->current.proj.origin[GMT_Y];
3032 y *= GMT->current.proj.i_scale[GMT_Y];
3033 gmtproj_ipolyconic_sub (GMT, y, GMT->common.R.wesn[XLO], &x);
3034 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
3035 }
3036
gmtproj_right_polyconic(struct GMT_CTRL * GMT,double y)3037 double gmtproj_right_polyconic (struct GMT_CTRL *GMT, double y) {
3038 double x;
3039 y -= GMT->current.proj.origin[GMT_Y];
3040 y *= GMT->current.proj.i_scale[GMT_Y];
3041 gmtproj_ipolyconic_sub (GMT, y, GMT->common.R.wesn[XHI], &x);
3042 return (x * GMT->current.proj.scale[GMT_X] + GMT->current.proj.origin[GMT_X]);
3043 }
3044