1 /* set ellipsoid parameters a and es */
2
3 #include <errno.h>
4 #include <math.h>
5 #include <stddef.h>
6 #include <string.h>
7
8 #include "proj.h"
9 #include "proj_internal.h"
10
11
12 /* Prototypes of the pj_ellipsoid helper functions */
13 static int ellps_ellps (PJ *P);
14 static int ellps_size (PJ *P);
15 static int ellps_shape (PJ *P);
16 static int ellps_spherification (PJ *P);
17
18 static paralist *pj_get_param (paralist *list, const char *key);
19 static char *pj_param_value (paralist *list);
20 static const PJ_ELLPS *pj_find_ellps (const char *name);
21
22
23 /***************************************************************************************/
pj_ellipsoid(PJ * P)24 int pj_ellipsoid (PJ *P) {
25 /****************************************************************************************
26 This is a replacement for the classic PROJ pj_ell_set function. The main difference
27 is that pj_ellipsoid augments the PJ object with a copy of the exact tags used to
28 define its related ellipsoid.
29
30 This makes it possible to let a new PJ object inherit the geometrical properties
31 of an existing one.
32
33 A complete ellipsoid definition comprises a size (primary) and a shape (secondary)
34 parameter.
35
36 Size parameters supported are:
37 R, defining the radius of a spherical planet
38 a, defining the semimajor axis of an ellipsoidal planet
39
40 Shape parameters supported are:
41 rf, the reverse flattening of the ellipsoid
42 f, the flattening of the ellipsoid
43 es, the eccentricity squared
44 e, the eccentricity
45 b, the semiminor axis
46
47 The ellps=xxx parameter provides both size and shape for a number of built in
48 ellipsoid definitions.
49
50 The ellipsoid definition may be augmented with a spherification flag, turning
51 the ellipsoid into a sphere with features defined by the ellipsoid.
52
53 Spherification parameters supported are:
54 R_A, which gives a sphere with the same surface area as the ellipsoid
55 R_A, which gives a sphere with the same volume as the ellipsoid
56
57 R_a, which gives a sphere with R = (a + b)/2 (arithmetic mean)
58 R_g, which gives a sphere with R = sqrt(a*b) (geometric mean)
59 R_h, which gives a sphere with R = 2*a*b/(a+b) (harmonic mean)
60
61 R_lat_a=phi, which gives a sphere with R being the arithmetic mean of
62 of the corresponding ellipsoid at latitude phi.
63 R_lat_g=phi, which gives a sphere with R being the geometric mean of
64 of the corresponding ellipsoid at latitude phi.
65
66 If R is given as size parameter, any shape and spherification parameters
67 given are ignored.
68
69 If size and shape are given as ellps=xxx, later shape and size parameters
70 are are taken into account as modifiers for the built in ellipsoid definition.
71
72 While this may seem strange, it is in accordance with historical PROJ
73 behavior. It can e.g. be used to define coordinates on the ellipsoid
74 scaled to unit semimajor axis by specifying "+ellps=xxx +a=1"
75
76 ****************************************************************************************/
77 int err = proj_errno_reset (P);
78 const char *empty = {""};
79
80 pj_dealloc(P->def_size);
81 P->def_size = nullptr;
82 pj_dealloc(P->def_shape);
83 P->def_shape = nullptr;
84 pj_dealloc(P->def_spherification);
85 P->def_spherification = nullptr;
86 pj_dealloc(P->def_ellps);
87 P->def_ellps = nullptr;
88
89 /* Specifying R overrules everything */
90 if (pj_get_param (P->params, "R")) {
91 if (0 != ellps_size (P))
92 return 1;
93 pj_calc_ellipsoid_params (P, P->a, 0);
94 if (proj_errno (P))
95 return 1;
96 return proj_errno_restore (P, err);
97 }
98
99
100 /* If an ellps argument is specified, start by using that */
101 if (0 != ellps_ellps (P))
102 return 1;
103
104 /* We may overwrite the size */
105 if (0 != ellps_size (P))
106 return 2;
107
108 /* We may also overwrite the shape */
109 if (0 != ellps_shape (P))
110 return 3;
111
112 /* When we're done with it, we compute all related ellipsoid parameters */
113 pj_calc_ellipsoid_params (P, P->a, P->es);
114
115 /* And finally, we may turn it into a sphere */
116 if (0 != ellps_spherification (P))
117 return 4;
118
119 proj_log_trace (P, "pj_ellipsoid - final: a=%.3f f=1/%7.3f, errno=%d",
120 P->a, P->f!=0? 1/P->f: 0, proj_errno (P));
121 proj_log_trace (P, "pj_ellipsoid - final: %s %s %s %s",
122 P->def_size? P->def_size: empty,
123 P->def_shape? P->def_shape: empty,
124 P->def_spherification? P->def_spherification: empty,
125 P->def_ellps? P->def_ellps: empty );
126
127 if (proj_errno (P))
128 return 5;
129
130 /* success */
131 return proj_errno_restore (P, err);
132 }
133
134
135 /***************************************************************************************/
ellps_ellps(PJ * P)136 static int ellps_ellps (PJ *P) {
137 /***************************************************************************************/
138 const PJ_ELLPS *ellps;
139 paralist *par = nullptr;
140 char *name;
141 int err;
142
143 /* Sail home if ellps=xxx is not specified */
144 par = pj_get_param (P->params, "ellps");
145 if (nullptr==par)
146 return 0;
147
148 /* Then look up the right size and shape parameters from the builtin list */
149 if (strlen (par->param) < 7)
150 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
151 name = par->param + 6;
152 ellps = pj_find_ellps (name);
153 if (nullptr==ellps)
154 return proj_errno_set (P, PJD_ERR_UNKNOWN_ELLP_PARAM);
155
156 /* Now, get things ready for ellps_size/ellps_shape, make them do their thing */
157 err = proj_errno_reset (P);
158
159 paralist* new_params = pj_mkparam (ellps->major);
160 if (nullptr == new_params)
161 return proj_errno_set (P, ENOMEM);
162 new_params->next = pj_mkparam (ellps->ell);
163 if (nullptr == new_params->next)
164 {
165 pj_dealloc(new_params);
166 return proj_errno_set (P, ENOMEM);
167 }
168 paralist* old_params = P->params;
169 P->params = new_params;
170
171 {
172 PJ empty_PJ;
173 pj_inherit_ellipsoid_def(&empty_PJ, P);
174 }
175 ellps_size (P);
176 ellps_shape (P);
177
178 P->params = old_params;
179 pj_dealloc (new_params->next);
180 pj_dealloc (new_params);
181 if (proj_errno (P))
182 return proj_errno (P);
183
184 /* Finally update P and sail home */
185 P->def_ellps = pj_strdup(par->param);
186 par->used = 1;
187
188 return proj_errno_restore (P, err);
189 }
190
191
192 /***************************************************************************************/
ellps_size(PJ * P)193 static int ellps_size (PJ *P) {
194 /***************************************************************************************/
195 paralist *par = nullptr;
196 int a_was_set = 0;
197
198 pj_dealloc(P->def_size);
199 P->def_size = nullptr;
200
201 /* A size parameter *must* be given, but may have been given as ellps prior */
202 if (P->a != 0)
203 a_was_set = 1;
204
205 /* Check which size key is specified */
206 par = pj_get_param (P->params, "R");
207 if (nullptr==par)
208 par = pj_get_param (P->params, "a");
209 if (nullptr==par)
210 return a_was_set? 0: proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
211
212 P->def_size = pj_strdup(par->param);
213 par->used = 1;
214 P->a = pj_atof (pj_param_value (par));
215 if (P->a <= 0)
216 return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
217 if (HUGE_VAL==P->a)
218 return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
219
220 if ('R'==par->param[0]) {
221 P->es = P->f = P->e = P->rf = 0;
222 P->b = P->a;
223 }
224 return 0;
225 }
226
227
228 /***************************************************************************************/
ellps_shape(PJ * P)229 static int ellps_shape (PJ *P) {
230 /***************************************************************************************/
231 const char *keys[] = {"rf", "f", "es", "e", "b"};
232 paralist *par = nullptr;
233 size_t i, len;
234
235 par = nullptr;
236 len = sizeof (keys) / sizeof (char *);
237
238 pj_dealloc(P->def_shape);
239 P->def_shape = nullptr;
240
241 /* Check which shape key is specified */
242 for (i = 0; i < len; i++) {
243 par = pj_get_param (P->params, keys[i]);
244 if (par)
245 break;
246 }
247
248 /* Not giving a shape parameter means selecting a sphere, unless shape */
249 /* has been selected previously via ellps=xxx */
250 if (nullptr==par && P->es != 0)
251 return 0;
252 if (nullptr==par && P->es==0) {
253 P->es = P->f = 0;
254 P->b = P->a;
255 return 0;
256 }
257
258 P->def_shape = pj_strdup(par->param);
259 par->used = 1;
260 P->es = P->f = P->b = P->e = P->rf = 0;
261
262 switch (i) {
263
264 /* reverse flattening, rf */
265 case 0:
266 P->rf = pj_atof (pj_param_value (par));
267 if (HUGE_VAL==P->rf)
268 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
269 if (0==P->rf)
270 return proj_errno_set (P, PJD_ERR_REV_FLATTENING_IS_ZERO);
271 P->f = 1 / P->rf;
272 P->es = 2*P->f - P->f*P->f;
273 break;
274
275 /* flattening, f */
276 case 1:
277 P->f = pj_atof (pj_param_value (par));
278 if (HUGE_VAL==P->f)
279 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
280
281 P->rf = P->f != 0.0 ? 1.0/P->f: HUGE_VAL;
282 P->es = 2*P->f - P->f*P->f;
283 break;
284
285 /* eccentricity squared, es */
286 case 2:
287 P->es = pj_atof (pj_param_value (par));
288 if (HUGE_VAL==P->es)
289 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
290 if (P->es >= 1)
291 return proj_errno_set (P, PJD_ERR_INVALID_ECCENTRICITY);
292 break;
293
294 /* eccentricity, e */
295 case 3:
296 P->e = pj_atof (pj_param_value (par));
297 if (HUGE_VAL==P->e)
298 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
299 if (P->e < 0 || P->e >= 1)
300 return proj_errno_set (P, PJD_ERR_INVALID_ECCENTRICITY);
301 P->es = P->e * P->e;
302 break;
303
304 /* semiminor axis, b */
305 case 4:
306 P->b = pj_atof (pj_param_value (par));
307 if (HUGE_VAL==P->b)
308 return proj_errno_set (P, PJD_ERR_INVALID_ARG);
309 if (P->b <= 0)
310 return proj_errno_set (P, PJD_ERR_INVALID_ECCENTRICITY);
311 if (P->b==P->a)
312 break;
313 P->f = (P->a - P->b) / P->a;
314 P->es = 2*P->f - P->f*P->f;
315 break;
316 default:
317 return PJD_ERR_INVALID_ARG;
318
319 }
320
321 // Written that way to catch NaN
322 if (!(P->es >= 0))
323 return proj_errno_set (P, PJD_ERR_ES_LESS_THAN_ZERO);
324 return 0;
325 }
326
327
328 /* series coefficients for calculating ellipsoid-equivalent spheres */
329 static const double SIXTH = 1/6.;
330 static const double RA4 = 17/360.;
331 static const double RA6 = 67/3024.;
332 static const double RV4 = 5/72.;
333 static const double RV6 = 55/1296.;
334
335 /***************************************************************************************/
ellps_spherification(PJ * P)336 static int ellps_spherification (PJ *P) {
337 /***************************************************************************************/
338 const char *keys[] = {"R_A", "R_V", "R_a", "R_g", "R_h", "R_lat_a", "R_lat_g"};
339 size_t len, i;
340 paralist *par = nullptr;
341
342 double t;
343 char *v, *endp;
344
345 len = sizeof (keys) / sizeof (char *);
346
347 /* Check which spherification key is specified */
348 for (i = 0; i < len; i++) {
349 par = pj_get_param (P->params, keys[i]);
350 if (par)
351 break;
352 }
353
354 /* No spherification specified? Then we're done */
355 if (i==len)
356 return 0;
357
358 /* Store definition */
359 P->def_spherification = pj_strdup(par->param);
360 par->used = 1;
361
362 switch (i) {
363
364 /* R_A - a sphere with same area as ellipsoid */
365 case 0:
366 P->a *= 1. - P->es * (SIXTH + P->es * (RA4 + P->es * RA6));
367 break;
368
369 /* R_V - a sphere with same volume as ellipsoid */
370 case 1:
371 P->a *= 1. - P->es * (SIXTH + P->es * (RV4 + P->es * RV6));
372 break;
373
374 /* R_a - a sphere with R = the arithmetic mean of the ellipsoid */
375 case 2:
376 P->a = (P->a + P->b) / 2;
377 break;
378
379 /* R_g - a sphere with R = the geometric mean of the ellipsoid */
380 case 3:
381 P->a = sqrt (P->a * P->b);
382 break;
383
384 /* R_h - a sphere with R = the harmonic mean of the ellipsoid */
385 case 4:
386 if (P->a + P->b == 0)
387 return proj_errno_set (P, PJD_ERR_TOLERANCE_CONDITION);
388 P->a = (2*P->a * P->b) / (P->a + P->b);
389 break;
390
391 /* R_lat_a - a sphere with R = the arithmetic mean of the ellipsoid at given latitude */
392 case 5:
393 /* R_lat_g - a sphere with R = the geometric mean of the ellipsoid at given latitude */
394 case 6:
395 v = pj_param_value (par);
396 t = proj_dmstor (v, &endp);
397 if (fabs (t) > M_HALFPI)
398 return proj_errno_set (P, PJD_ERR_REF_RAD_LARGER_THAN_90);
399 t = sin (t);
400 t = 1 - P->es * t * t;
401 if (t == 0.) {
402 return proj_errno_set(P, PJD_ERR_INVALID_ECCENTRICITY);
403 }
404 if (i==5) /* arithmetic */
405 P->a *= (1. - P->es + t) / (2 * t * sqrt(t));
406 else /* geometric */
407 P->a *= sqrt (1 - P->es) / t;
408 break;
409 }
410
411 if (P->a <= 0.) {
412 return proj_errno_set(P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
413 }
414
415 /* Clean up the ellipsoidal parameters to reflect the sphere */
416 P->es = P->e = P->f = 0;
417 P->rf = HUGE_VAL;
418 P->b = P->a;
419 pj_calc_ellipsoid_params (P, P->a, 0);
420
421 return 0;
422 }
423
424
425 /* locate parameter in list */
pj_get_param(paralist * list,const char * key)426 static paralist *pj_get_param (paralist *list, const char *key) {
427 size_t l = strlen(key);
428 while (list && !(0==strncmp(list->param, key, l) && (0==list->param[l] || list->param[l] == '=') ) )
429 list = list->next;
430 return list;
431 }
432
433
pj_param_value(paralist * list)434 static char *pj_param_value (paralist *list) {
435 char *key, *value;
436 if (nullptr==list)
437 return nullptr;
438
439 key = list->param;
440 value = strchr (key, '=');
441
442 /* a flag (i.e. a key without value) has its own name (key) as value */
443 return value? value + 1: key;
444 }
445
446
pj_find_ellps(const char * name)447 static const PJ_ELLPS *pj_find_ellps (const char *name) {
448 int i;
449 const char *s;
450 const PJ_ELLPS *ellps;
451
452 if (nullptr==name)
453 return nullptr;
454
455 ellps = proj_list_ellps();
456
457 /* Search through internal ellipsoid list for name */
458 for (i = 0; (s = ellps[i].id) && strcmp(name, s) ; ++i);
459 if (nullptr==s)
460 return nullptr;
461 return ellps + i;
462 }
463
464
465 /**************************************************************************************/
pj_inherit_ellipsoid_def(const PJ * src,PJ * dst)466 void pj_inherit_ellipsoid_def (const PJ *src, PJ *dst) {
467 /***************************************************************************************
468 Brute force copy the ellipsoidal parameters from src to dst. This code was
469 written before the actual ellipsoid setup parameters were kept available in
470 the PJ->def_xxx elements.
471 ***************************************************************************************/
472
473 /* The linear parameters */
474 dst->a = src->a;
475 dst->b = src->b;
476 dst->ra = src->ra;
477 dst->rb = src->rb;
478
479 /* The eccentricities */
480 dst->alpha = src->alpha;
481 dst->e = src->e;
482 dst->es = src->es;
483 dst->e2 = src->e2;
484 dst->e2s = src->e2s;
485 dst->e3 = src->e3;
486 dst->e3s = src->e3s;
487 dst->one_es = src->one_es;
488 dst->rone_es = src->rone_es;
489
490 /* The flattenings */
491 dst->f = src->f;
492 dst->f2 = src->f2;
493 dst->n = src->n;
494 dst->rf = src->rf;
495 dst->rf2 = src->rf2;
496 dst->rn = src->rn;
497
498 /* This one's for GRS80 */
499 dst->J = src->J;
500
501 /* es and a before any +proj related adjustment */
502 dst->es_orig = src->es_orig;
503 dst->a_orig = src->a_orig;
504 }
505
506
507 /***************************************************************************************/
pj_calc_ellipsoid_params(PJ * P,double a,double es)508 int pj_calc_ellipsoid_params (PJ *P, double a, double es) {
509 /****************************************************************************************
510 Calculate a large number of ancillary ellipsoidal parameters, in addition to
511 the two traditional PROJ defining parameters: Semimajor axis, a, and the
512 eccentricity squared, es.
513
514 Most of these parameters are fairly cheap to compute in comparison to the overall
515 effort involved in initializing a PJ object. They may, however, take a substantial
516 part of the time taken in computing an individual point transformation.
517
518 So by providing them up front, we can amortize the (already modest) cost over all
519 transformations carried out over the entire lifetime of a PJ object, rather than
520 incur that cost for every single transformation.
521
522 Most of the parameter calculations here are based on the "angular eccentricity",
523 i.e. the angle, measured from the semiminor axis, of a line going from the north
524 pole to one of the foci of the ellipsoid - or in other words: The arc sine of the
525 eccentricity.
526
527 The formulae used are mostly taken from:
528
529 Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
530 Columbus, Ohio: Dept. of Geodetic Science
531 and Surveying, Ohio State University.
532
533 ****************************************************************************************/
534
535 P->a = a;
536 P->es = es;
537
538 /* Compute some ancillary ellipsoidal parameters */
539 if (P->e==0)
540 P->e = sqrt(P->es); /* eccentricity */
541 P->alpha = asin (P->e); /* angular eccentricity */
542
543 /* second eccentricity */
544 P->e2 = tan (P->alpha);
545 P->e2s = P->e2 * P->e2;
546
547 /* third eccentricity */
548 P->e3 = (0!=P->alpha)? sin (P->alpha) / sqrt(2 - sin (P->alpha)*sin (P->alpha)): 0;
549 P->e3s = P->e3 * P->e3;
550
551 /* flattening */
552 if (0==P->f)
553 P->f = 1 - cos (P->alpha); /* = 1 - sqrt (1 - PIN->es); */
554 if (P->f == 1.0) {
555 pj_ctx_set_errno( P->ctx, PJD_ERR_INVALID_ECCENTRICITY);
556 return PJD_ERR_INVALID_ECCENTRICITY;
557 }
558 P->rf = P->f != 0.0 ? 1.0/P->f: HUGE_VAL;
559
560 /* second flattening */
561 P->f2 = (cos(P->alpha)!=0)? 1/cos (P->alpha) - 1: 0;
562 P->rf2 = P->f2 != 0.0 ? 1/P->f2: HUGE_VAL;
563
564 /* third flattening */
565 P->n = pow (tan (P->alpha/2), 2);
566 P->rn = P->n != 0.0 ? 1/P->n: HUGE_VAL;
567
568 /* ...and a few more */
569 if (0==P->b)
570 P->b = (1 - P->f)*P->a;
571 P->rb = 1. / P->b;
572 P->ra = 1. / P->a;
573
574 P->one_es = 1. - P->es;
575 if (P->one_es == 0.) {
576 pj_ctx_set_errno( P->ctx, PJD_ERR_INVALID_ECCENTRICITY);
577 return PJD_ERR_INVALID_ECCENTRICITY;
578 }
579
580 P->rone_es = 1./P->one_es;
581
582 return 0;
583 }
584
585
586
587 #ifndef KEEP_ORIGINAL_PJ_ELL_SET
588 /**************************************************************************************/
pj_ell_set(PJ_CONTEXT * ctx,paralist * pl,double * a,double * es)589 int pj_ell_set (PJ_CONTEXT *ctx, paralist *pl, double *a, double *es) {
590 /***************************************************************************************
591 Initialize ellipsoidal parameters by emulating the original ellipsoid setup
592 function by Gerald Evenden, through a call to pj_ellipsoid
593 ***************************************************************************************/
594 PJ B;
595 int ret;
596
597 B.ctx = ctx;
598 B.params = pl;
599
600 ret = pj_ellipsoid (&B);
601 if (ret)
602 return ret;
603
604 *a = B.a;
605 *es = B.es;
606 return 0;
607 }
608 #else
609
610
611 /**************************************************************************************/
pj_ell_set(projCtx ctx,paralist * pl,double * a,double * es)612 int pj_ell_set (projCtx ctx, paralist *pl, double *a, double *es) {
613 /***************************************************************************************
614 Initialize ellipsoidal parameters: This is the original ellipsoid setup
615 function by Gerald Evenden - significantly more compact than pj_ellipsoid and
616 its many helper functions, and still quite readable.
617
618 It is, however, also so tight that it is hard to modify and add functionality,
619 and equally hard to find the right place to add further commentary for improved
620 future maintainability.
621
622 Hence, when the need to store in the PJ object, the parameters actually used to
623 define the ellipsoid came up, rather than modifying this little gem of
624 "economy of expression", a much more verbose reimplementation, pj_ellipsoid,
625 was written.
626 ***************************************************************************************/
627 int i;
628 double b=0.0, e;
629 char *name;
630 paralist *start = 0;
631
632 /* clear any previous error */
633 pj_ctx_set_errno(ctx,0);
634
635 /* check for varying forms of ellipsoid input */
636 *a = *es = 0.;
637
638 /* R takes precedence */
639 if (pj_param(ctx, pl, "tR").i)
640 *a = pj_param(ctx,pl, "dR").f;
641
642 /* probable elliptical figure */
643 else {
644 /* check if ellps present and temporarily append its values to pl */
645 if ((name = pj_param(ctx,pl, "sellps").s) != NULL) {
646 char *s;
647
648 for (start = pl; start && start->next ; start = start->next) ;
649 for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i) ;
650 if (!s) {
651 pj_ctx_set_errno( ctx, PJD_ERR_UNKNOWN_ELLP_PARAM);
652 return 1;
653 }
654 start->next = pj_mkparam(pj_ellps[i].major);
655 start->next->next = pj_mkparam(pj_ellps[i].ell);
656 }
657
658 *a = pj_param(ctx,pl, "da").f;
659
660 if (pj_param(ctx,pl, "tes").i) /* eccentricity squared */
661 *es = pj_param(ctx,pl, "des").f;
662 else if (pj_param(ctx,pl, "te").i) { /* eccentricity */
663 e = pj_param(ctx,pl, "de").f;
664 if (e < 0) {
665 pj_ctx_set_errno(ctx, PJD_ERR_INVALID_ECCENTRICITY);
666 return 1;
667 }
668 *es = e * e;
669 } else if (pj_param(ctx,pl, "trf").i) { /* recip flattening */
670 *es = pj_param(ctx,pl, "drf").f;
671 if (*es == 0.0) {
672 pj_ctx_set_errno(ctx, PJD_ERR_REV_FLATTENING_IS_ZERO);
673 goto bomb;
674 }
675 *es = 1./ *es;
676 *es = *es * (2. - *es);
677 } else if (pj_param(ctx,pl, "tf").i) { /* flattening */
678 *es = pj_param(ctx,pl, "df").f;
679 *es = *es * (2. - *es);
680 } else if (pj_param(ctx,pl, "tb").i) { /* minor axis */
681 b = pj_param(ctx,pl, "db").f;
682 *es = 1. - (b * b) / (*a * *a);
683 } /* else *es == 0. and sphere of radius *a */
684 if (b == 0.0)
685 b = *a * sqrt(1. - *es);
686
687
688 /* following options turn ellipsoid into equivalent sphere */
689 if (pj_param(ctx,pl, "bR_A").i) { /* sphere--area of ellipsoid */
690 *a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
691 *es = 0.;
692 } else if (pj_param(ctx,pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
693 *a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
694 *es = 0.;
695 } else if (pj_param(ctx,pl, "bR_a").i) { /* sphere--arithmetic mean */
696 *a = .5 * (*a + b);
697 *es = 0.;
698 } else if (pj_param(ctx,pl, "bR_g").i) { /* sphere--geometric mean */
699 *a = sqrt(*a * b);
700 *es = 0.;
701 } else if (pj_param(ctx,pl, "bR_h").i) { /* sphere--harmonic mean */
702 if ( (*a + b) == 0.0) {
703 pj_ctx_set_errno(ctx, PJD_ERR_TOLERANCE_CONDITION);
704 goto bomb;
705 }
706 *a = 2. * *a * b / (*a + b);
707 *es = 0.;
708 } else if ((i = pj_param(ctx,pl, "tR_lat_a").i) || /* sphere--arith. */
709 pj_param(ctx,pl, "tR_lat_g").i) { /* or geom. mean at latitude */
710 double tmp;
711
712 tmp = sin(pj_param(ctx,pl, i ? "rR_lat_a" : "rR_lat_g").f);
713 if (fabs(tmp) > M_HALFPI) {
714 pj_ctx_set_errno(ctx, PJD_ERR_REF_RAD_LARGER_THAN_90);
715 goto bomb;
716 }
717 tmp = 1. - *es * tmp * tmp;
718 *a *= i ? .5 * (1. - *es + tmp) / ( tmp * sqrt(tmp)) :
719 sqrt(1. - *es) / tmp;
720 *es = 0.;
721 }
722 bomb:
723 if (start) { /* clean up temporary extension of list */
724 pj_dalloc(start->next->next);
725 pj_dalloc(start->next);
726 start->next = 0;
727 }
728 if (ctx->last_errno)
729 return 1;
730 }
731 /* some remaining checks */
732 if (*es < 0.) {
733 pj_ctx_set_errno(ctx, PJD_ERR_ES_LESS_THAN_ZERO);
734 return 1;
735 }
736 if (*es >= 1.) {
737 pj_ctx_set_errno(ctx, PJD_ERR_INVALID_ECCENTRICITY);
738 return 1;
739 }
740 if (*a <= 0.) {
741 pj_ctx_set_errno(ctx, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
742 return 1;
743 }
744 return 0;
745 }
746 #endif
747