1 /*============================================================================
2 WCSLIB 7.7 - an implementation of the FITS WCS standard.
3 Copyright (C) 1995-2021, Mark Calabretta
4
5 This file is part of WCSLIB.
6
7 WCSLIB is free software: you can redistribute it and/or modify it under the
8 terms of the GNU Lesser General Public License as published by the Free
9 Software Foundation, either version 3 of the License, or (at your option)
10 any later version.
11
12 WCSLIB is distributed in the hope that it will be useful, but WITHOUT ANY
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with WCSLIB. If not, see http://www.gnu.org/licenses.
19
20 Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
21 http://www.atnf.csiro.au/people/Mark.Calabretta
22 $Id: dis.c,v 7.7 2021/07/12 06:36:49 mcalabre Exp $
23 *===========================================================================*/
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29
30 #include "wcserr.h"
31 #include "wcsprintf.h"
32 #include "wcsutil.h"
33 #include "dis.h"
34
35 const int DISSET = 137;
36
37 const int DIS_TPD = 1;
38 const int DIS_POLYNOMIAL = 2;
39 const int DIS_DOTPD = 1024;
40
41 // Maximum number of DPja or DQia keywords.
42 int NDPMAX = 256;
43
44 // Map status return value to message.
45 const char *dis_errmsg[] = {
46 "Success",
47 "Null disprm pointer passed",
48 "Memory allocation failed",
49 "Invalid parameter value",
50 "Distort error",
51 "De-distort error"};
52
53 // Convenience macro for invoking wcserr_set().
54 #define DIS_ERRMSG(status) WCSERR_SET(status), dis_errmsg[status]
55
56 // Internal helper functions, not for general use.
57 static int polyset(int j, struct disprm *dis);
58 static int tpdset(int j, struct disprm *dis);
59
60 static int pol2tpd(int j, struct disprm *dis);
61 static int tpvset(int j, struct disprm *dis);
62 static int sipset(int j, struct disprm *dis);
63 static int dssset(int j, struct disprm *dis);
64 static int watset(int j, struct disprm *dis);
65 static int cheleg(int type, int m, int n, double coeffm[], double coeffn[]);
66
67 static int dispoly(DISP2X_ARGS);
68 static int tpd1(DISP2X_ARGS);
69 static int tpd2(DISP2X_ARGS);
70 static int tpd3(DISP2X_ARGS);
71 static int tpd4(DISP2X_ARGS);
72 static int tpd5(DISP2X_ARGS);
73 static int tpd6(DISP2X_ARGS);
74 static int tpd7(DISP2X_ARGS);
75 static int tpd8(DISP2X_ARGS);
76 static int tpd9(DISP2X_ARGS);
77
78 // The first three iparm indices have meanings common to all distortion
79 // functions. They are used by disp2x(), disx2p(), disprt(), and dishdo().
80 #define I_DTYPE 0 // Distortion type code.
81 #define I_NIPARM 1 // Full (allocated) length of iparm[].
82 #define I_NDPARM 2 // No. of parameters in dparm[], excl. work space.
83
84 //----------------------------------------------------------------------------
85
disndp(int ndpmax)86 int disndp(int ndpmax) { if (ndpmax >= 0) NDPMAX = ndpmax; return NDPMAX; }
87
88 //----------------------------------------------------------------------------
89
dpfill(struct dpkey * dp,const char * keyword,const char * field,int j,int type,int i,double f)90 int dpfill(
91 struct dpkey *dp,
92 const char *keyword,
93 const char *field,
94 int j,
95 int type,
96 int i,
97 double f)
98
99 {
100 if (keyword) {
101 if (field) {
102 if (j && 2 <= strlen(keyword)) {
103 // Fill in the axis number from the value given.
104 if (keyword[2] == '\0') {
105 sprintf(dp->field, "%s%d.%s", keyword, j, field);
106 } else {
107 // Take care not to overwrite any alternate code.
108 char axno[8];
109 sprintf(dp->field, "%s.%s", keyword, field);
110 sprintf(axno, "%d", j);
111 dp->field[2] = axno[0];
112 }
113
114 } else {
115 sprintf(dp->field, "%s.%s", keyword, field);
116 }
117 } else {
118 strcpy(dp->field, keyword);
119 }
120 } else if (field) {
121 strcpy(dp->field, field);
122 }
123
124 if (j) {
125 dp->j = j;
126 } else {
127 // The field name must either be given or preset.
128 char *cp;
129 if ((cp = strpbrk(dp->field, "0123456789")) != 0x0) {
130 sscanf(cp, "%d.", &(dp->j));
131 }
132 }
133
134 if ((dp->type = type)) {
135 dp->value.f = f;
136 } else {
137 dp->value.i = i;
138 }
139
140 return 0;
141 }
142
143 //----------------------------------------------------------------------------
144
dpkeyi(const struct dpkey * dp)145 int dpkeyi(const struct dpkey *dp)
146
147 {
148 if (dp->type != 0) {
149 return (int)dp->value.f;
150 }
151
152 return dp->value.i;
153 }
154
155 //----------------------------------------------------------------------------
156
dpkeyd(const struct dpkey * dp)157 double dpkeyd(const struct dpkey *dp)
158
159 {
160 if (dp->type == 0) {
161 return (double)dp->value.i;
162 }
163
164 return dp->value.f;
165 }
166
167 //----------------------------------------------------------------------------
168
disini(int alloc,int naxis,struct disprm * dis)169 int disini(int alloc, int naxis, struct disprm *dis)
170
171 {
172 return disinit(alloc, naxis, dis, -1);
173 }
174
175 //----------------------------------------------------------------------------
176
disinit(int alloc,int naxis,struct disprm * dis,int ndpmax)177 int disinit(int alloc, int naxis, struct disprm *dis, int ndpmax)
178
179 {
180 static const char *function = "disinit";
181
182 // Check inputs.
183 if (dis == 0x0) return DISERR_NULL_POINTER;
184
185 if (ndpmax < 0) ndpmax = disndp(-1);
186
187
188 // Initialize error message handling.
189 if (dis->flag == -1) {
190 dis->err = 0x0;
191 }
192 struct wcserr **err = &(dis->err);
193 wcserr_clear(err);
194
195
196 // Initialize pointers.
197 if (dis->flag == -1 || dis->m_flag != DISSET) {
198 if (dis->flag == -1) {
199 dis->docorr = 0x0;
200 dis->Nhat = 0x0;
201
202 dis->axmap = 0x0;
203 dis->offset = 0x0;
204 dis->scale = 0x0;
205 dis->iparm = 0x0;
206 dis->dparm = 0x0;
207
208 dis->disp2x = 0x0;
209 dis->disx2p = 0x0;
210 dis->tmpmem = 0x0;
211
212 dis->i_naxis = 0;
213 }
214
215 // Initialize memory management.
216 dis->m_flag = 0;
217 dis->m_naxis = 0;
218 dis->m_dtype = 0x0;
219 dis->m_dp = 0x0;
220 dis->m_maxdis = 0x0;
221 }
222
223 if (naxis < 0) {
224 return wcserr_set(WCSERR_SET(DISERR_MEMORY),
225 "naxis must not be negative (got %d)", naxis);
226 }
227
228
229 // Allocate memory for arrays if required.
230 if (alloc ||
231 dis->dtype == 0x0 ||
232 (ndpmax && dis->dp == 0x0) ||
233 dis->maxdis == 0x0) {
234
235 // Was sufficient allocated previously?
236 if (dis->m_flag == DISSET &&
237 (dis->m_naxis < naxis ||
238 dis->ndpmax < ndpmax)) {
239 // No, free it.
240 disfree(dis);
241 }
242
243 if (alloc || dis->dtype == 0x0) {
244 if (dis->m_dtype) {
245 // In case the caller fiddled with it.
246 dis->dtype = dis->m_dtype;
247
248 } else {
249 if ((dis->dtype = calloc(naxis, sizeof(char [72]))) == 0x0) {
250 disfree(dis);
251 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
252 }
253
254 dis->m_flag = DISSET;
255 dis->m_naxis = naxis;
256 dis->m_dtype = dis->dtype;
257 }
258 }
259
260 if (alloc || dis->dp == 0x0) {
261 if (dis->m_dp) {
262 // In case the caller fiddled with it.
263 dis->dp = dis->m_dp;
264
265 } else {
266 if (ndpmax) {
267 if ((dis->dp = calloc(ndpmax, sizeof(struct dpkey))) == 0x0) {
268 disfree(dis);
269 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
270 }
271 } else {
272 dis->dp = 0x0;
273 }
274
275 dis->ndpmax = ndpmax;
276
277 dis->m_flag = DISSET;
278 dis->m_naxis = naxis;
279 dis->m_dp = dis->dp;
280 }
281 }
282
283 if (alloc || dis->maxdis == 0x0) {
284 if (dis->m_maxdis) {
285 // In case the caller fiddled with it.
286 dis->maxdis = dis->m_maxdis;
287
288 } else {
289 if ((dis->maxdis = calloc(naxis, sizeof(double))) == 0x0) {
290 disfree(dis);
291 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
292 }
293
294 dis->m_flag = DISSET;
295 dis->m_naxis = naxis;
296 dis->m_maxdis = dis->maxdis;
297 }
298 }
299 }
300
301
302 // Set defaults.
303 dis->flag = 0;
304 dis->naxis = naxis;
305
306 if (naxis) {
307 memset(dis->dtype, 0, naxis*sizeof(char [72]));
308 }
309
310 dis->ndp = 0;
311 if (ndpmax) {
312 memset(dis->dp, 0, ndpmax*sizeof(struct dpkey));
313 }
314
315 if (naxis) {
316 memset(dis->maxdis, 0, naxis*sizeof(double));
317 }
318 dis->totdis = 0.0;
319
320 return 0;
321 }
322
323
324 //----------------------------------------------------------------------------
325
discpy(int alloc,const struct disprm * dissrc,struct disprm * disdst)326 int discpy(int alloc, const struct disprm *dissrc, struct disprm *disdst)
327
328 {
329 static const char *function = "discpy";
330
331 if (dissrc == 0x0) return DISERR_NULL_POINTER;
332 if (disdst == 0x0) return DISERR_NULL_POINTER;
333 struct wcserr **err = &(disdst->err);
334
335 int naxis = dissrc->naxis;
336 if (naxis < 1) {
337 return wcserr_set(WCSERR_SET(DISERR_MEMORY),
338 "naxis must be positive (got %d)", naxis);
339 }
340
341 int status;
342 if ((status = disinit(alloc, naxis, disdst, dissrc->ndpmax))) {
343 return status;
344 }
345
346 memcpy(disdst->dtype, dissrc->dtype, naxis*sizeof(char [72]));
347
348 disdst->ndp = dissrc->ndp;
349 memcpy(disdst->dp, dissrc->dp, dissrc->ndpmax*sizeof(struct dpkey));
350
351 memcpy(disdst->maxdis, dissrc->maxdis, naxis*sizeof(double));
352 disdst->totdis = dissrc->totdis;
353
354 return 0;
355 }
356
357 //----------------------------------------------------------------------------
358
disfree(struct disprm * dis)359 int disfree(struct disprm *dis)
360
361 {
362 if (dis == 0x0) return DISERR_NULL_POINTER;
363
364 if (dis->flag != -1) {
365 // Optionally allocated by disinit() for given parameters.
366 if (dis->m_flag == DISSET) {
367 if (dis->dtype == dis->m_dtype) dis->dtype = 0x0;
368 if (dis->dp == dis->m_dp) dis->dp = 0x0;
369 if (dis->maxdis == dis->m_maxdis) dis->maxdis = 0x0;
370
371 if (dis->m_dtype) free(dis->m_dtype);
372 if (dis->m_dp) free(dis->m_dp);
373 if (dis->m_maxdis) free(dis->m_maxdis);
374 }
375
376 // The remainder were allocated by disset().
377 if (dis->docorr) free(dis->docorr);
378 if (dis->Nhat) free(dis->Nhat);
379
380 // Recall that axmap, offset, and scale were allocated in bulk.
381 if (dis->axmap && dis->axmap[0]) free(dis->axmap[0]);
382 if (dis->offset && dis->offset[0]) free(dis->offset[0]);
383 if (dis->scale && dis->scale[0]) free(dis->scale[0]);
384
385 if (dis->axmap) free(dis->axmap);
386 if (dis->offset) free(dis->offset);
387 if (dis->scale) free(dis->scale);
388
389 if (dis->iparm) {
390 for (int j = 0; j < dis->i_naxis; j++) {
391 if (dis->iparm[j]) free(dis->iparm[j]);
392 }
393 free(dis->iparm);
394 }
395
396 if (dis->dparm) {
397 for (int j = 0; j < dis->i_naxis; j++) {
398 if (dis->dparm[j]) free(dis->dparm[j]);
399 }
400 free(dis->dparm);
401 }
402
403 if (dis->disp2x) free(dis->disp2x);
404 if (dis->disx2p) free(dis->disx2p);
405 if (dis->tmpmem) free(dis->tmpmem);
406 }
407
408 dis->m_flag = 0;
409 dis->m_naxis = 0;
410 dis->m_dtype = 0x0;
411 dis->m_dp = 0x0;
412 dis->m_maxdis = 0x0;
413
414 dis->docorr = 0x0;
415 dis->Nhat = 0x0;
416 dis->axmap = 0x0;
417 dis->offset = 0x0;
418 dis->scale = 0x0;
419 dis->iparm = 0x0;
420 dis->dparm = 0x0;
421 dis->disp2x = 0x0;
422 dis->disx2p = 0x0;
423 dis->tmpmem = 0x0;
424
425 wcserr_clear(&(dis->err));
426
427 dis->flag = 0;
428
429 return 0;
430 }
431
432 //----------------------------------------------------------------------------
433
dissize(const struct disprm * dis,int sizes[2])434 int dissize(const struct disprm *dis, int sizes[2])
435
436 {
437 if (dis == 0x0) {
438 sizes[0] = sizes[1] = 0;
439 return DISERR_NULL_POINTER;
440 }
441
442 // Base size, in bytes.
443 sizes[0] = sizeof(struct disprm);
444
445 // Total size of allocated memory, in bytes.
446 sizes[1] = 0;
447
448 int naxis = dis->naxis;
449
450 // disprm::dtype[].
451 sizes[1] += naxis * sizeof(char [72]);
452
453 // disprm::dp[].
454 sizes[1] += dis->ndpmax * sizeof(struct dpkey);
455
456 // disprm::maxdis[].
457 sizes[1] += naxis * sizeof(double);
458
459 // dis::err[].
460 int exsizes[2];
461 wcserr_size(dis->err, exsizes);
462 sizes[1] += exsizes[0] + exsizes[1];
463
464 // The remaining arrays are allocated by disset().
465 if (dis->flag != DISSET) {
466 return DISERR_SUCCESS;
467 }
468
469 // dis::docorr[].
470 sizes[1] += naxis * sizeof(int *);
471
472 // dis::Nhat[].
473 sizes[1] += naxis * sizeof(int *);
474
475 // dis::axmap[][].
476 sizes[1] += naxis * sizeof(int *);
477 sizes[1] += naxis*naxis * sizeof(int);
478
479 // dis::offset[][].
480 sizes[1] += naxis * sizeof(double *);
481 sizes[1] += naxis*naxis * sizeof(double);
482
483 // dis::scale[][].
484 sizes[1] += naxis * sizeof(double *);
485 sizes[1] += naxis*naxis * sizeof(double);
486
487 // dis::iparm[][].
488 sizes[1] += naxis * sizeof(int *);
489 for (int j = 0; j < naxis; j++) {
490 if (dis->iparm[j]) {
491 sizes[1] += dis->iparm[j][I_NIPARM] * sizeof(int);
492 }
493 }
494
495 // dis::dparm[][].
496 sizes[1] += naxis * sizeof(double *);
497 for (int j = 0; j < naxis; j++) {
498 if (dis->dparm[j]) {
499 sizes[1] += dis->dparm[j][I_NDPARM] * sizeof(double);
500 }
501 }
502
503 // dis::disp2x[].
504 sizes[1] += naxis * sizeof(int (*)(DISP2X_ARGS));
505
506 // dis::disx2p[].
507 sizes[1] += naxis * sizeof(int (*)(DISX2P_ARGS));
508
509 // dis::tmpmem[].
510 sizes[1] += 5*naxis * sizeof(double);
511
512 return DISERR_SUCCESS;
513 }
514
515 //----------------------------------------------------------------------------
516
disprt(const struct disprm * dis)517 int disprt(const struct disprm *dis)
518
519 {
520 if (dis == 0x0) return DISERR_NULL_POINTER;
521
522 if (dis->flag != DISSET) {
523 wcsprintf("The disprm struct is UNINITIALIZED.\n");
524 return 0;
525 }
526
527 int naxis = dis->naxis;
528
529
530 wcsprintf(" flag: %d\n", dis->flag);
531
532 // Parameters supplied.
533 wcsprintf(" naxis: %d\n", naxis);
534
535 WCSPRINTF_PTR(" dtype: ", dis->dtype, "\n");
536 for (int j = 0; j < naxis; j++) {
537 wcsprintf(" \"%s\"\n", dis->dtype[j]);
538 }
539
540 wcsprintf(" ndp: %d\n", dis->ndp);
541 wcsprintf(" ndpmax: %d\n", dis->ndpmax);
542 WCSPRINTF_PTR(" dp: ", dis->dp, "\n");
543 for (int i = 0; i < dis->ndp; i++) {
544 if (dis->dp[i].type) {
545 wcsprintf(" %3d%3d %#- 11.5g %.32s\n",
546 dis->dp[i].j, dis->dp[i].type, dis->dp[i].value.f, dis->dp[i].field);
547 } else {
548 wcsprintf(" %3d%3d %11d %.32s\n",
549 dis->dp[i].j, dis->dp[i].type, dis->dp[i].value.i, dis->dp[i].field);
550 }
551 }
552
553 WCSPRINTF_PTR(" maxdis: ", dis->maxdis, "\n");
554 wcsprintf(" ");
555 for (int j = 0; j < naxis; j++) {
556 wcsprintf(" %#- 11.5g", dis->maxdis[j]);
557 }
558 wcsprintf("\n");
559
560 wcsprintf(" totdis: %#- 11.5g\n", dis->totdis);
561
562 // Derived values.
563 WCSPRINTF_PTR(" docorr: ", dis->docorr, "\n");
564 wcsprintf(" ");
565 for (int j = 0; j < naxis; j++) {
566 wcsprintf("%6d", dis->docorr[j]);
567 }
568 wcsprintf("\n");
569
570 WCSPRINTF_PTR(" Nhat: ", dis->Nhat, "\n");
571 wcsprintf(" ");
572 for (int j = 0; j < naxis; j++) {
573 wcsprintf("%6d", dis->Nhat[j]);
574 }
575 wcsprintf("\n");
576
577 WCSPRINTF_PTR(" axmap: ", dis->axmap, "\n");
578 for (int j = 0; j < naxis; j++) {
579 wcsprintf(" axmap[%d][]:", j);
580 for (int jhat = 0; jhat < naxis; jhat++) {
581 wcsprintf("%6d", dis->axmap[j][jhat]);
582 }
583 wcsprintf("\n");
584 }
585
586 WCSPRINTF_PTR(" offset: ", dis->offset, "\n");
587 for (int j = 0; j < naxis; j++) {
588 wcsprintf("offset[%d][]:", j);
589 for (int jhat = 0; jhat < naxis; jhat++) {
590 wcsprintf(" %#- 11.5g", dis->offset[j][jhat]);
591 }
592 wcsprintf("\n");
593 }
594
595 WCSPRINTF_PTR(" scale: ", dis->scale, "\n");
596 for (int j = 0; j < naxis; j++) {
597 wcsprintf(" scale[%d][]:", j);
598 for (int jhat = 0; jhat < naxis; jhat++) {
599 wcsprintf(" %#- 11.5g", dis->scale[j][jhat]);
600 }
601 wcsprintf("\n");
602 }
603
604 WCSPRINTF_PTR(" iparm: ", dis->iparm, "\n");
605 for (int j = 0; j < naxis; j++) {
606 wcsprintf(" iparm[%d] : ", j);
607 WCSPRINTF_PTR("", dis->iparm[j], "\n");
608
609 if (dis->iparm[j]) {
610 wcsprintf(" iparm[%d][]:", j);
611 for (int k = 0; k < dis->iparm[j][I_NIPARM]; k++) {
612 if (k && k%5 == 0) {
613 wcsprintf("\n ");
614 }
615 wcsprintf(" %11d", dis->iparm[j][k]);
616 }
617 wcsprintf("\n");
618 }
619 }
620
621 WCSPRINTF_PTR(" dparm: ", dis->dparm, "\n");
622 for (int j = 0; j < naxis; j++) {
623 wcsprintf(" dparm[%d] : ", j);
624 WCSPRINTF_PTR("", dis->dparm[j], "\n");
625
626 if (dis->dparm[j]) {
627 wcsprintf(" dparm[%d][]:", j);
628 for (int k = 0; k < dis->iparm[j][I_NDPARM]; k++) {
629 if (k && k%5 == 0) {
630 wcsprintf("\n ");
631 }
632 wcsprintf(" %#- 11.5g", dis->dparm[j][k]);
633 }
634 wcsprintf("\n");
635 }
636 }
637
638 wcsprintf(" i_naxis: %d\n", dis->i_naxis);
639 wcsprintf(" ndis: %d\n", dis->ndis);
640
641 // Error handling.
642 WCSPRINTF_PTR(" err: ", dis->err, "\n");
643 if (dis->err) {
644 wcserr_prt(dis->err, " ");
645 }
646
647 // Work arrays.
648 char hext[32];
649 WCSPRINTF_PTR(" disp2x: ", dis->disp2x, "\n");
650 for (int j = 0; j < naxis; j++) {
651 wcsprintf(" disp2x[%d]: %s", j,
652 wcsutil_fptr2str((void (*)(void))dis->disp2x[j], hext));
653 if (dis->disp2x[j] == dispoly) {
654 wcsprintf(" (= dispoly)\n");
655 } else if (dis->disp2x[j] == tpd1) {
656 wcsprintf(" (= tpd1)\n");
657 } else if (dis->disp2x[j] == tpd2) {
658 wcsprintf(" (= tpd2)\n");
659 } else if (dis->disp2x[j] == tpd3) {
660 wcsprintf(" (= tpd3)\n");
661 } else if (dis->disp2x[j] == tpd4) {
662 wcsprintf(" (= tpd4)\n");
663 } else if (dis->disp2x[j] == tpd5) {
664 wcsprintf(" (= tpd5)\n");
665 } else if (dis->disp2x[j] == tpd6) {
666 wcsprintf(" (= tpd6)\n");
667 } else if (dis->disp2x[j] == tpd7) {
668 wcsprintf(" (= tpd7)\n");
669 } else if (dis->disp2x[j] == tpd8) {
670 wcsprintf(" (= tpd8)\n");
671 } else if (dis->disp2x[j] == tpd9) {
672 wcsprintf(" (= tpd9)\n");
673 } else {
674 wcsprintf("\n");
675 }
676 }
677 WCSPRINTF_PTR(" disx2p: ", dis->disx2p, "\n");
678 for (int j = 0; j < naxis; j++) {
679 wcsprintf(" disx2p[%d]: %s\n", j,
680 wcsutil_fptr2str((void (*)(void))dis->disx2p[j], hext));
681 }
682 WCSPRINTF_PTR(" tmpmem: ", dis->tmpmem, "\n");
683
684 // Memory management.
685 wcsprintf(" m_flag: %d\n", dis->m_flag);
686 wcsprintf(" m_naxis: %d\n", dis->m_naxis);
687 WCSPRINTF_PTR(" m_dtype: ", dis->m_dtype, "");
688 if (dis->m_dtype == dis->dtype) wcsprintf(" (= dtype)");
689 wcsprintf("\n");
690 WCSPRINTF_PTR(" m_dp: ", dis->m_dp, "");
691 if (dis->m_dp == dis->dp) wcsprintf(" (= dp)");
692 wcsprintf("\n");
693 WCSPRINTF_PTR(" m_maxdis: ", dis->m_maxdis, "");
694 if (dis->m_maxdis == dis->maxdis) wcsprintf(" (= maxdis)");
695 wcsprintf("\n");
696
697 return 0;
698 }
699
700 //----------------------------------------------------------------------------
701
disperr(const struct disprm * dis,const char * prefix)702 int disperr(const struct disprm *dis, const char *prefix)
703
704 {
705 if (dis == 0x0) return DISERR_NULL_POINTER;
706
707 if (dis->err) {
708 wcserr_prt(dis->err, prefix);
709 }
710
711 return 0;
712 }
713
714 //----------------------------------------------------------------------------
715
dishdo(struct disprm * dis)716 int dishdo(struct disprm *dis)
717
718 {
719 static const char *function = "dishdo";
720
721 if (dis == 0x0) return DISERR_NULL_POINTER;
722 struct wcserr **err = &(dis->err);
723
724 int status = 0;
725 for (int j = 0; j < dis->naxis; j++) {
726 if (dis->iparm[j][I_DTYPE]) {
727 if (dis->iparm[j][I_DTYPE] == DIS_TPD) {
728 // Implemented as TPD...
729 if (strcmp(dis->dtype[j], "TPD") != 0) {
730 // ... but isn't TPD.
731 dis->iparm[j][I_DTYPE] |= DIS_DOTPD;
732 }
733 } else {
734 // Must be a Polynomial that can't be implemented as TPD.
735 status = wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
736 "Translation of %s to TPD is not possible", dis->dtype[j]);
737 }
738 }
739 }
740
741 return status;
742 }
743
744 //----------------------------------------------------------------------------
745
disset(struct disprm * dis)746 int disset(struct disprm *dis)
747
748 {
749 static const char *function = "disset";
750
751 if (dis == 0x0) return DISERR_NULL_POINTER;
752 struct wcserr **err = &(dis->err);
753
754 int naxis = dis->naxis;
755
756
757 // Do basic checks.
758 if (dis->ndp < 0) {
759 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
760 "disprm::ndp is negative (%d)", dis->ndp);
761 }
762
763 int ndis = 0;
764 for (int j = 0; j < naxis; j++) {
765 if (strlen(dis->dtype[j])) {
766 ndis++;
767 break;
768 }
769 }
770
771 char *dpq;
772 if (dis->ndp) {
773 // Is it prior or sequent?
774 if (dis->dp[0].field[1] == 'P') {
775 dpq = "DPja";
776 } else if (dis->dp[0].field[1] == 'Q') {
777 dpq = "DQia";
778 } else {
779 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
780 "disprm::dp[0].field (%s) is invalid", dis->dp[0].field);
781 }
782
783 } else {
784 if (ndis) {
785 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
786 "No DPja or DQia keywords, NAXES at least is required for each "
787 "distortion");
788 }
789
790 // A Clayton's distortion. Avert compiler warnings about possible use of
791 // uninitialized variables.
792 dpq = 0x0;
793 }
794
795
796 // Free memory allocated separately for each axis.
797 for (int j = 0; j < dis->i_naxis; j++) {
798 if (dis->iparm[j]) free(dis->iparm[j]);
799 if (dis->dparm[j]) free(dis->dparm[j]);
800 dis->iparm[j] = 0x0;
801 dis->dparm[j] = 0x0;
802 }
803
804 // Allocate or reallocate memory, if necessary, for derived parameter and
805 // work arrays sized according to the number of axes.
806 if (dis->i_naxis < naxis) {
807 if (dis->i_naxis) {
808 free(dis->docorr);
809 free(dis->Nhat);
810
811 // Noting that axmap, offset, and scale are allocated in bulk.
812 free(dis->axmap[0]);
813 free(dis->axmap);
814 free(dis->offset[0]);
815 free(dis->offset);
816 free(dis->scale[0]);
817 free(dis->scale);
818
819 free(dis->iparm);
820 free(dis->dparm);
821
822 free(dis->disp2x);
823 free(dis->disx2p);
824
825 free(dis->tmpmem);
826 }
827
828 if ((dis->docorr = calloc(naxis, sizeof(int *))) == 0x0) {
829 disfree(dis);
830 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
831 }
832
833 if ((dis->Nhat = calloc(naxis, sizeof(int *))) == 0x0) {
834 disfree(dis);
835 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
836 }
837
838 // Allocate axmap[][] in bulk and then carve it up.
839 if ((dis->axmap = calloc(naxis, sizeof(int *))) == 0x0) {
840 disfree(dis);
841 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
842 }
843
844 if ((dis->axmap[0] = calloc(naxis*naxis, sizeof(int))) == 0x0) {
845 disfree(dis);
846 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
847 }
848
849 for (int j = 1; j < naxis; j++) {
850 dis->axmap[j] = dis->axmap[j-1] + naxis;
851 }
852
853 // Allocate offset[][] in bulk and then carve it up.
854 if ((dis->offset = calloc(naxis, sizeof(double *))) == 0x0) {
855 disfree(dis);
856 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
857 }
858
859 if ((dis->offset[0] = calloc(naxis*naxis, sizeof(double))) == 0x0) {
860 disfree(dis);
861 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
862 }
863
864 for (int j = 1; j < naxis; j++) {
865 dis->offset[j] = dis->offset[j-1] + naxis;
866 }
867
868 // Allocate scale[][] in bulk and then carve it up.
869 if ((dis->scale = calloc(naxis, sizeof(double *))) == 0x0) {
870 disfree(dis);
871 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
872 }
873
874 if ((dis->scale[0] = calloc(naxis*naxis, sizeof(double))) == 0x0) {
875 disfree(dis);
876 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
877 }
878
879 for (int j = 1; j < naxis; j++) {
880 dis->scale[j] = dis->scale[j-1] + naxis;
881 }
882
883 if ((dis->iparm = calloc(naxis, sizeof(int *))) == 0x0) {
884 disfree(dis);
885 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
886 }
887
888 if ((dis->dparm = calloc(naxis, sizeof(double *))) == 0x0) {
889 disfree(dis);
890 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
891 }
892
893 if ((dis->disp2x = calloc(naxis, sizeof(int (*)(DISP2X_ARGS)))) == 0x0) {
894 disfree(dis);
895 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
896 }
897
898 if ((dis->disx2p = calloc(naxis, sizeof(int (*)(DISX2P_ARGS)))) == 0x0) {
899 disfree(dis);
900 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
901 }
902
903 if ((dis->tmpmem = calloc(5*naxis, sizeof(double))) == 0x0) {
904 disfree(dis);
905 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
906 }
907
908 dis->i_naxis = naxis;
909 }
910
911 // Start with a clean slate.
912 for (int j = 0; j < naxis; j++) {
913 dis->docorr[j] = 1;
914 }
915
916 memset(dis->Nhat, 0, naxis*sizeof(int));
917
918 for (int jhat = 0; jhat < naxis*naxis; jhat++) {
919 dis->axmap[0][jhat] = -1;
920 }
921
922 memset(dis->offset[0], 0, naxis*naxis*sizeof(double));
923
924 for (int jhat = 0; jhat < naxis*naxis; jhat++) {
925 dis->scale[0][jhat] = 1.0;
926 }
927
928 // polyset() etc. must look after iparm[][] and dparm[][].
929
930 dis->i_naxis = naxis;
931 dis->ndis = 0;
932
933 memset(dis->disp2x, 0, naxis*sizeof(int (*)(DISP2X_ARGS)));
934 memset(dis->disx2p, 0, naxis*sizeof(int (*)(DISX2P_ARGS)));
935 memset(dis->tmpmem, 0, naxis*sizeof(double));
936
937
938 // Handle DPja or DQia keywords common to all distortions.
939 struct dpkey *keyp = dis->dp;
940 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
941 // Check that they're all one kind or the other.
942 if (keyp->field[1] != dpq[1]) {
943 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
944 "disprm::dp appears to contain a mix of DPja and DQia keys");
945 }
946
947 int j = keyp->j;
948
949 if (j < 1 || naxis < j) {
950 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
951 "Invalid axis number (%d) in %s", j, keyp->field);
952 }
953
954 char *fp;
955 if ((fp = strchr(keyp->field, '.')) == 0x0) {
956 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
957 "Invalid record field name: %s", j, keyp->field);
958 }
959 fp++;
960
961 // Convert to 0-relative axis number.
962 j--;
963
964 if (strncmp(fp, "DOCORR", 7) == 0) {
965 if (dpkeyi(keyp) == 0) {
966 dis->docorr[j] = 0;
967 }
968
969 } else if (strncmp(fp, "NAXES", 6) == 0) {
970 int Nhat = dpkeyi(keyp);
971 if (Nhat < 0 || naxis < Nhat) {
972 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
973 "Invalid value of Nhat for %s distortion in %s: %d", dis->dtype[j],
974 keyp->field, Nhat);
975 }
976
977 dis->Nhat[j] = Nhat;
978
979 } else if (strncmp(fp, "AXIS.", 5) == 0) {
980 int jhat;
981 sscanf(fp+5, "%d", &jhat);
982 if (jhat < 1 || naxis < jhat) {
983 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
984 "Invalid axis in axis map for %s distortion in %s: %d",
985 dis->dtype[j], keyp->field, jhat);
986 }
987
988 // N.B. axis numbers in the map are 0-relative.
989 dis->axmap[j][jhat-1] = dpkeyi(keyp) - 1;
990
991 } else if (strncmp(fp, "OFFSET.", 7) == 0) {
992 int jhat;
993 sscanf(fp+7, "%d", &jhat);
994 dis->offset[j][jhat-1] = dpkeyd(keyp);
995
996 } else if (strncmp(fp, "SCALE.", 6) == 0) {
997 int jhat;
998 sscanf(fp+6, "%d", &jhat);
999 dis->scale[j][jhat-1] = dpkeyd(keyp);
1000 }
1001 }
1002
1003 // Set defaults and do sanity checks on axmap[][].
1004 for (int j = 0; j < naxis; j++) {
1005 if (strlen(dis->dtype[j]) == 0) {
1006 // No distortion on this axis, check that there are no parameters.
1007 keyp = dis->dp;
1008 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
1009 if (keyp->j == j+1) {
1010 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1011 "No distortion type, yet %s keyvalues are present for axis %d",
1012 dpq, j+1);
1013 }
1014 }
1015
1016 continue;
1017 }
1018
1019 // N.B. NAXES (Nhat) has no default value.
1020 if (dis->Nhat[j] <= 0) {
1021 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1022 "%s.NAXES was not set (or bad) for %s distortion on axis %d",
1023 dpq, dis->dtype[j], j+1);
1024 }
1025
1026 // Set defaults for axmap[][].
1027 int Nhat = dis->Nhat[j];
1028 for (int jhat = 0; jhat < Nhat; jhat++) {
1029 if (dis->axmap[j][jhat] == -1) {
1030 dis->axmap[j][jhat] = jhat;
1031 }
1032 }
1033
1034 // Sanity check on the length of the axis map.
1035 Nhat = 0;
1036 for (int jhat = 0; jhat < naxis; jhat++) {
1037 if (dis->axmap[j][jhat] != -1) Nhat = jhat+1;
1038 }
1039
1040 if (Nhat != dis->Nhat[j]) {
1041 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1042 "Mismatch in length of axis map for %s distortion on axis %d",
1043 dis->dtype[j], j+1);
1044 }
1045
1046 // Check uniqueness of entries in the axis map.
1047 for (int jhat = 0; jhat < Nhat; jhat++) {
1048 for (int k = 0; k < jhat; k++) {
1049 if (dis->axmap[j][jhat] == dis->axmap[j][k]) {
1050 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1051 "Duplicated entry in axis map for %s distortion on axis %d",
1052 dis->dtype[j], j+1);
1053 }
1054 }
1055 }
1056 }
1057
1058
1059 // Identify the distortion functions.
1060 ndis = 0;
1061 for (int j = 0; j < naxis; j++) {
1062 if (strlen(dis->dtype[j]) == 0) {
1063 // No distortion on this axis.
1064 continue;
1065 }
1066
1067 if (dis->Nhat[j] == 0) {
1068 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1069 "Empty axis map for %s distortion on axis %d", dis->dtype[j], j+1);
1070 }
1071
1072 // Invoke the specific setup functions for each distortion.
1073 int status;
1074 if (strcmp(dis->dtype[j], "TPD") == 0) {
1075 // Template Polynomial Distortion.
1076 if ((status = tpdset(j, dis))) {
1077 // (Preserve the error message set by tpdset().)
1078 return status;
1079 }
1080
1081 } else if (strcmp(dis->dtype[j], "TPV") == 0) {
1082 // TPV "projection".
1083 if ((status = tpvset(j, dis))) {
1084 // (Preserve the error message set by tpvset().)
1085 return status;
1086 }
1087
1088 } else if (strcmp(dis->dtype[j], "SIP") == 0) {
1089 // Simple Imaging Polynomial (SIP).
1090 if ((status = sipset(j, dis))) {
1091 // (Preserve the error message set by sipset().)
1092 return status;
1093 }
1094
1095 } else if (strcmp(dis->dtype[j], "DSS") == 0) {
1096 // Digitized Sky Survey (DSS).
1097 if ((status = dssset(j, dis))) {
1098 // (Preserve the error message set by dssset().)
1099 return status;
1100 }
1101
1102 } else if (strncmp(dis->dtype[j], "WAT", 3) == 0) {
1103 // WAT (TNX or ZPX "projections").
1104 if ((status = watset(j, dis))) {
1105 // (Preserve the error message set by watset().)
1106 return status;
1107 }
1108
1109 } else if (strcmp(dis->dtype[j], "Polynomial") == 0 ||
1110 strcmp(dis->dtype[j], "Polynomial*") == 0) {
1111 // General polynomial distortion.
1112 if ((status = polyset(j, dis))) {
1113 // (Preserve the error message set by polyset().)
1114 return status;
1115 }
1116
1117 } else {
1118 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1119 "Unrecognized/unimplemented distortion function: %s", dis->dtype[j]);
1120 }
1121
1122 ndis++;
1123 }
1124
1125 dis->ndis = ndis;
1126 dis->flag = DISSET;
1127
1128 return 0;
1129 }
1130
1131 //----------------------------------------------------------------------------
1132
disp2x(struct disprm * dis,const double rawcrd[],double discrd[])1133 int disp2x(
1134 struct disprm *dis,
1135 const double rawcrd[],
1136 double discrd[])
1137
1138 {
1139 static const char *function = "disp2x";
1140
1141 // Initialize.
1142 if (dis == 0x0) return DISERR_NULL_POINTER;
1143 struct wcserr **err = &(dis->err);
1144
1145 if (dis->flag != DISSET) {
1146 int status;
1147 if ((status = disset(dis))) return status;
1148 }
1149
1150 int naxis = dis->naxis;
1151
1152
1153 // Invoke the distortion functions for each axis.
1154 double *tmpcrd = dis->tmpmem;
1155 for (int j = 0; j < naxis; j++) {
1156 if (dis->disp2x[j]) {
1157 double *offset = dis->offset[j];
1158 double *scale = dis->scale[j];
1159
1160 int Nhat = dis->Nhat[j];
1161 for (int jhat = 0; jhat < Nhat; jhat++) {
1162 int axisj = dis->axmap[j][jhat];
1163 tmpcrd[jhat] = (rawcrd[axisj] - offset[jhat])*scale[jhat];
1164 }
1165
1166 double dtmp;
1167 if ((dis->disp2x[j])(0, dis->iparm[j], dis->dparm[j], Nhat, tmpcrd,
1168 &dtmp)) {
1169 return wcserr_set(DIS_ERRMSG(DISERR_DISTORT));
1170 }
1171
1172 if (dis->docorr[j]) {
1173 // Distortion function computes a correction to be applied.
1174 discrd[j] = rawcrd[j] + dtmp;
1175 } else {
1176 // Distortion function computes corrected coordinates directly.
1177 discrd[j] = dtmp;
1178 }
1179
1180 } else {
1181 discrd[j] = rawcrd[j];
1182 }
1183 }
1184
1185 return 0;
1186 }
1187
1188 //----------------------------------------------------------------------------
1189 // This function is intended for debugging purposes only.
1190 // No documentation or prototype is provided in dis.h.
1191
disitermax(int itermax)1192 int disitermax(int itermax)
1193
1194 {
1195 static int ITERMAX = 30;
1196
1197 if (itermax >= 0) {
1198 ITERMAX = itermax;
1199 }
1200
1201 return ITERMAX;
1202 }
1203
1204 //----------------------------------------------------------------------------
1205
disx2p(struct disprm * dis,const double discrd[],double rawcrd[])1206 int disx2p(
1207 struct disprm *dis,
1208 const double discrd[],
1209 double rawcrd[])
1210
1211 {
1212 static const char *function = "disx2p";
1213
1214 const double TOL = 1.0e-13;
1215
1216 int status;
1217
1218 // Initialize.
1219 if (dis == 0x0) return DISERR_NULL_POINTER;
1220 struct wcserr **err = &(dis->err);
1221
1222 int naxis = dis->naxis;
1223
1224 // Carve up working memory, noting that disp2x() gets to it first.
1225 double *dcrd0 = dis->tmpmem + naxis;
1226 double *dcrd1 = dcrd0 + naxis;
1227 double *rcrd1 = dcrd1 + naxis;
1228 double *delta = rcrd1 + naxis;
1229
1230
1231 // Zeroth approximation. The assumption here and below is that the
1232 // distortion is small so that, to first order in the neighbourhood of
1233 // the solution, discrd[j] ~= a + b*rawcrd[j], i.e. independent of
1234 // rawcrd[i], where i != j. This is effectively equivalent to assuming
1235 // that the distortion functions are separable to first order.
1236 // Furthermore, a is assumed to be small, and b close to unity.
1237 memcpy(rawcrd, discrd, naxis*sizeof(double));
1238
1239 // If available, use disprm::disx2p to improve the zeroth approximation.
1240 for (int j = 0; j < naxis; j++) {
1241 if (dis->disx2p[j]) {
1242 double *offset = dis->offset[j];
1243 double *scale = dis->scale[j];
1244 double *tmpcrd = dis->tmpmem;
1245
1246 int Nhat = dis->Nhat[j];
1247 for (int jhat = 0; jhat < Nhat; jhat++) {
1248 int axisj = dis->axmap[j][jhat];
1249 tmpcrd[jhat] = (discrd[axisj] - offset[jhat])*scale[jhat];
1250 }
1251
1252 double rtmp;
1253 if ((status = (dis->disx2p[j])(1, dis->iparm[j], dis->dparm[j], Nhat,
1254 tmpcrd, &rtmp))) {
1255 return wcserr_set(DIS_ERRMSG(DISERR_DEDISTORT));
1256 }
1257
1258 if (dis->docorr[j]) {
1259 // Inverse distortion function computes a correction to be applied.
1260 rawcrd[j] = discrd[j] + rtmp;
1261 } else {
1262 // Inverse distortion function computes corrected coordinates directly.
1263 rawcrd[j] = rtmp;
1264 }
1265 }
1266 }
1267
1268 // Quick return debugging hook, assumes inverse functions were defined.
1269 int itermax;
1270 if ((itermax = disitermax(-1)) == 0) {
1271 return 0;
1272 }
1273
1274
1275 // Iteratively invert the (well-behaved!) distortion function.
1276 int convergence, iter;
1277 for (iter = 0; iter < itermax; iter++) {
1278 if ((status = disp2x(dis, rawcrd, dcrd0))) {
1279 return wcserr_set(DIS_ERRMSG(status));
1280 }
1281
1282 // Check for convergence.
1283 convergence = 1;
1284 for (int j = 0; j < naxis; j++) {
1285 delta[j] = discrd[j] - dcrd0[j];
1286
1287 double dd;
1288 if (fabs(discrd[j]) < 1.0) {
1289 dd = delta[j];
1290 } else {
1291 // TOL may be below the precision achievable from floating point
1292 // subtraction, so switch to a fractional tolerance.
1293 dd = delta[j] / discrd[j];
1294 }
1295
1296 if (TOL < fabs(dd)) {
1297 // No convergence yet on this axis.
1298 convergence = 0;
1299 }
1300 }
1301
1302 if (convergence) break;
1303
1304 // Determine a suitable test point for computing the gradient.
1305 for (int j = 0; j < naxis; j++) {
1306 // Constrain the displacement.
1307 delta[j] /= 2.0;
1308 if (fabs(delta[j]) < 1.0e-6) {
1309 if (delta[j] < 0.0) {
1310 delta[j] = -1.0e-6;
1311 } else {
1312 delta[j] = 1.0e-6;
1313 }
1314 } else if (1.0 < fabs(delta[j])) {
1315 if (delta[j] < 0.0) {
1316 delta[j] = -1.0;
1317 } else {
1318 delta[j] = 1.0;
1319 }
1320 }
1321 }
1322
1323 if (iter < itermax/2) {
1324 // With the assumption of small distortions (as above), the gradient
1325 // of discrd[j] should be dominated by the partial derivative with
1326 // respect to rawcrd[j], and we can neglect partials with respect
1327 // to rawcrd[i], where i != j. Thus only one test point is needed,
1328 // not one for each axis.
1329 for (int j = 0; j < naxis; j++) {
1330 rcrd1[j] = rawcrd[j] + delta[j];
1331 }
1332
1333 // Compute discrd[] at the test point.
1334 if ((status = disp2x(dis, rcrd1, dcrd1))) {
1335 return wcserr_set(DIS_ERRMSG(status));
1336 }
1337
1338 // Compute the next approximation.
1339 for (int j = 0; j < naxis; j++) {
1340 rawcrd[j] += (discrd[j] - dcrd0[j]) *
1341 (delta[j]/(dcrd1[j] - dcrd0[j]));
1342 }
1343
1344 } else {
1345 // Convergence should not take more than seven or so iterations. As
1346 // it is slow, try computing the gradient in full.
1347 memcpy(rcrd1, rawcrd, naxis*sizeof(double));
1348
1349 for (int j = 0; j < naxis; j++) {
1350 rcrd1[j] += delta[j];
1351
1352 // Compute discrd[] at the test point.
1353 if ((status = disp2x(dis, rcrd1, dcrd1))) {
1354 return wcserr_set(DIS_ERRMSG(status));
1355 }
1356
1357 // Compute the next approximation.
1358 rawcrd[j] += (discrd[j] - dcrd0[j]) *
1359 (delta[j]/(dcrd1[j] - dcrd0[j]));
1360
1361 rcrd1[j] -= delta[j];
1362 }
1363 }
1364 }
1365
1366
1367 if (!convergence) {
1368 double residual = 0.0;
1369 for (int j = 0; j < naxis; j++) {
1370 double dd = discrd[j] - dcrd0[j] ;
1371 residual += dd*dd;
1372 }
1373 residual = sqrt(residual);
1374
1375 return wcserr_set(WCSERR_SET(DISERR_DEDISTORT),
1376 "Convergence not achieved after %d iterations, residual %#7.2g", iter,
1377 residual);
1378 }
1379
1380
1381 return 0;
1382 }
1383
1384 //----------------------------------------------------------------------------
1385
diswarp(struct disprm * dis,const double pixblc[],const double pixtrc[],const double pixsamp[],int * nsamp,double maxdis[],double * maxtot,double avgdis[],double * avgtot,double rmsdis[],double * rmstot)1386 int diswarp(
1387 struct disprm *dis,
1388 const double pixblc[],
1389 const double pixtrc[],
1390 const double pixsamp[],
1391 int *nsamp,
1392 double maxdis[],
1393 double *maxtot,
1394 double avgdis[],
1395 double *avgtot,
1396 double rmsdis[],
1397 double *rmstot)
1398
1399 {
1400 static const char *function = "diswarp";
1401
1402 int status = 0;
1403
1404 // Initialize.
1405 if (dis == 0x0) return DISERR_NULL_POINTER;
1406 struct wcserr **err = &(dis->err);
1407
1408 int naxis = dis->naxis;
1409
1410 if (nsamp) *nsamp = 0;
1411 for (int j = 0; j < naxis; j++) {
1412 if (maxdis) maxdis[j] = 0.0;
1413 if (avgdis) avgdis[j] = 0.0;
1414 if (rmsdis) rmsdis[j] = 0.0;
1415 }
1416 if (maxtot) *maxtot = 0.0;
1417 if (avgtot) *avgtot = 0.0;
1418 if (rmstot) *rmstot = 0.0;
1419
1420 // Quick return if no distortions.
1421 if (dis->ndis == 0) return 0;
1422
1423 // Carve up working memory, noting that disp2x() gets to it first.
1424 double *pixinc = dis->tmpmem + naxis;
1425 double *pixend = pixinc + naxis;
1426 double *sumdis = pixend + naxis;
1427 double *ssqdis = sumdis + naxis;
1428
1429 // Work out increments on each axis.
1430 for (int j = 0; j < naxis; j++) {
1431 double pixspan = pixtrc[j] - (pixblc ? pixblc[j] : 1.0);
1432
1433 if (pixsamp == 0x0) {
1434 pixinc[j] = 1.0;
1435 } else if (pixsamp[j] == 0.0) {
1436 pixinc[j] = 1.0;
1437 } else if (pixsamp[j] > 0.0) {
1438 pixinc[j] = pixsamp[j];
1439 } else if (pixsamp[j] > -1.5) {
1440 pixinc[j] = 2.0*pixspan;
1441 } else {
1442 pixinc[j] = pixspan / ((int)(-pixsamp[j] - 0.5));
1443 }
1444 }
1445
1446 // Get some more memory for coordinate vectors.
1447 double *pix0, *pix1;
1448 if ((pix0 = calloc(2*naxis, sizeof(double))) == 0x0) {
1449 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
1450 }
1451
1452 pix1 = pix0 + naxis;
1453
1454
1455 // Set up the array of pixel coordinates.
1456 for (int j = 0; j < naxis; j++) {
1457 pix0[j] = pixblc ? pixblc[j] : 1.0;
1458 pixend[j] = pixtrc[j] + 0.5*pixinc[j];
1459 }
1460
1461 // Initialize accumulators.
1462 for (int j = 0; j < naxis; j++) {
1463 sumdis[j] = 0.0;
1464 ssqdis[j] = 0.0;
1465 }
1466 double sumtot = 0.0;
1467 double ssqtot = 0.0;
1468
1469
1470 // Loop over N dimensions.
1471 int carry = 0;
1472 while (carry == 0) {
1473 if ((status = disp2x(dis, pix0, pix1))) {
1474 // (Preserve the error message set by disp2x().)
1475 goto cleanup;
1476 }
1477
1478 // Accumulate statistics.
1479 (*nsamp)++;
1480
1481 double dssq = 0.0;
1482 for (int j = 0; j < naxis; j++) {
1483 double dpix = pix1[j] - pix0[j];
1484 double dpx2 = dpix*dpix;
1485
1486 sumdis[j] += dpix;
1487 ssqdis[j] += dpx2;
1488
1489 if (maxdis && (dpix = fabs(dpix)) > maxdis[j]) {
1490 maxdis[j] = dpix;
1491 }
1492
1493 dssq += dpx2;
1494 }
1495
1496 double totdis = sqrt(dssq);
1497 sumtot += totdis;
1498 ssqtot += totdis*totdis;
1499
1500 if (maxtot && *maxtot < totdis) {
1501 *maxtot = totdis;
1502 }
1503
1504 // Next pixel.
1505 for (int j = 0; j < naxis; j++) {
1506 pix0[j] += pixinc[j];
1507 if (pix0[j] < pixend[j]) {
1508 carry = 0;
1509 break;
1510 }
1511
1512 pix0[j] = pixblc ? pixblc[j] : 1.0;
1513 carry = 1;
1514 }
1515 }
1516
1517
1518 // Compute the means and RMSs.
1519 for (int j = 0; j < naxis; j++) {
1520 ssqdis[j] /= *nsamp;
1521 sumdis[j] /= *nsamp;
1522 if (avgdis) avgdis[j] = sumdis[j];
1523 if (rmsdis) rmsdis[j] = sqrt(ssqdis[j] - sumdis[j]*sumdis[j]);
1524 }
1525
1526 ssqtot /= *nsamp;
1527 sumtot /= *nsamp;
1528 if (avgtot) *avgtot = sumtot;
1529 if (rmstot) *rmstot = sqrt(ssqtot - sumtot*sumtot);
1530
1531
1532 cleanup:
1533 free(pix0);
1534
1535 return status;
1536 }
1537
1538 //----------------------------------------------------------------------------
1539
polyset(int j,struct disprm * dis)1540 int polyset(int j, struct disprm *dis)
1541
1542 {
1543 static const char *function = "polyset";
1544
1545 // Initialize.
1546 if (dis == 0x0) return DISERR_NULL_POINTER;
1547 struct wcserr **err = &(dis->err);
1548
1549 int naxis = dis->naxis;
1550
1551 char id[32];
1552 sprintf(id, "Polynomial on axis %d", j+1);
1553
1554
1555 // Find the number of auxiliary variables and terms.
1556 int K = 0;
1557 int M = 0;
1558 struct dpkey *keyp = dis->dp;
1559 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
1560 if (keyp->j-1 != j) continue;
1561
1562 char *fp;
1563 if ((fp = strchr(keyp->field, '.')) == 0x0) {
1564 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1565 "Invalid field name for %s: %s", id, keyp->field);
1566 }
1567 fp++;
1568
1569 if (strcmp(fp, "NAUX") == 0) {
1570 K = dpkeyi(keyp);
1571 } else if (strcmp(fp, "NTERMS") == 0) {
1572 M = dpkeyi(keyp);
1573 }
1574 }
1575
1576 if (K < 0) {
1577 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1578 "Invalid number of auxiliaries (%d) for %s", K, id);
1579 }
1580
1581 if (M <= 0) {
1582 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1583 "Invalid number of terms (%d) for %s", M, id);
1584 }
1585
1586 int Nhat = dis->Nhat[j];
1587 int nKparm = 2*(Nhat + 1);
1588 int nVar = Nhat + K;
1589 int nTparm = 1 + nVar;
1590 int ndparm = K*nKparm + M*nTparm;
1591
1592 // These iparm indices are specific to Polynomial.
1593 #define I_NIDX 3 // No. of indexes in iparm[].
1594 #define I_LENDP 4 // Full (allocated) length of dparm[].
1595 #define I_K 5 // No. of auxiliary variables.
1596 #define I_M 6 // No. of terms in the polynomial.
1597 #define I_NKPARM 7 // No. of parameters used to define each auxiliary.
1598 #define I_NTPARM 8 // No. of parameters used to define each term.
1599 #define I_NVAR 9 // No. of independent + auxiliary variables.
1600 #define I_MNVAR 10 // No. of powers (exponents) in the polynomial.
1601 #define I_DPOLY 11 // dparm offset for polynomial coefficients.
1602 #define I_DAUX 12 // dparm offset for auxiliary coefficients.
1603 #define I_DVPOW 13 // dparm offset for integral powers of variables.
1604 #define I_MAXPOW 14 // iparm offset for max powers.
1605 #define I_DPOFF 15 // iparm offset for dparm offsets.
1606 #define I_FLAGS 16 // iparm offset for flags.
1607 #define I_IPOW 17 // iparm offset for integral powers.
1608 #define I_NPOLY 18
1609
1610 // Add extra for handling integer exponents. See "Optimization" below.
1611 int niparm = I_NPOLY + (2 + 2*M)*nVar;
1612
1613 // Add extra memory for temporaries.
1614 int lendp = ndparm + K;
1615
1616 // Allocate memory for the indexes and parameter array.
1617 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
1618 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
1619 }
1620
1621 if ((dis->dparm[j] = calloc(lendp, sizeof(double))) == 0x0) {
1622 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
1623 }
1624
1625 // These help a bit to stop the code from turning into hieroglyphics.
1626 int *iparm = dis->iparm[j];
1627 double *dparm = dis->dparm[j];
1628
1629
1630 // Record the indexing parameters. The first three are more widely used.
1631 iparm[I_DTYPE] = DIS_POLYNOMIAL;
1632 iparm[I_NIPARM] = niparm;
1633 iparm[I_NDPARM] = ndparm;
1634
1635 iparm[I_NIDX] = I_NPOLY;
1636 iparm[I_LENDP] = lendp;
1637 iparm[I_K] = K;
1638 iparm[I_M] = M;
1639 iparm[I_NKPARM] = nKparm;
1640 iparm[I_NTPARM] = nTparm;
1641 iparm[I_NVAR] = nVar;
1642 iparm[I_MNVAR] = M*nVar;
1643 iparm[I_DPOLY] = K*nKparm;
1644 iparm[I_DAUX] = ndparm;
1645 iparm[I_DVPOW] = ndparm + K;
1646 iparm[I_MAXPOW] = iparm[I_NIDX];
1647 iparm[I_DPOFF] = iparm[I_MAXPOW] + nVar;
1648 iparm[I_FLAGS] = iparm[I_DPOFF] + nVar;
1649 iparm[I_IPOW] = iparm[I_FLAGS] + M*nVar;
1650
1651 // Set default values of POWER for the auxiliary variables.
1652 double *dptr = dparm + (1 + Nhat);
1653 for (int k = 0; k < K; k++) {
1654 for (int jhat = 0; jhat <= Nhat; jhat++) {
1655 dptr[jhat] = 1.0;
1656 }
1657 dptr += nKparm;
1658 }
1659
1660 // Set default values of COEFF for the independent variables.
1661 dptr = dparm + iparm[I_DPOLY];
1662 for (int m = 0; m < M; m++) {
1663 *dptr = 1.0;
1664 dptr += nTparm;
1665 }
1666
1667 // Extract parameter values from DPja or DQia.
1668 int i, k, m;
1669 k = m = 0;
1670 keyp = dis->dp;
1671 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
1672 // N.B. keyp->j is 1-relative, but j is 0-relative.
1673 if (keyp->j-1 != j) continue;
1674
1675 char *fp = strchr(keyp->field, '.') + 1;
1676
1677 if (strncmp(fp, "AUX.", 4) == 0) {
1678 // N.B. k here is 1-relative.
1679 fp += 4;
1680 sscanf(fp, "%d", &k);
1681 if (k < 1 || K < k) {
1682 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1683 "Bad auxiliary variable (%d) for %s: %s", k, id, keyp->field);
1684 }
1685
1686 if ((fp = strchr(fp, '.')) == 0x0) {
1687 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1688 "Invalid field name for %s: %s", id, keyp->field);
1689 }
1690 fp++;
1691
1692 int offset;
1693 if (strncmp(fp, "COEFF.", 6) == 0) {
1694 offset = 0;
1695
1696 } else if (strncmp(fp, "POWER.", 6) == 0) {
1697 offset = 1 + Nhat;
1698
1699 } else {
1700 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1701 "Unrecognized field name for %s: %s", id, keyp->field);
1702 }
1703
1704 fp += 6;
1705 int jhat;
1706 sscanf(fp, "%d", &jhat);
1707 if (jhat < 0 || naxis < jhat) {
1708 // N.B. jhat == 0 is ok.
1709 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1710 "Invalid axis number (%d) for %s: %s", jhat, id, keyp->field);
1711 }
1712
1713 i = (k-1)*nKparm + offset + jhat;
1714 dparm[i] = dpkeyd(keyp);
1715
1716 } else if (strncmp(fp, "TERM.", 5) == 0) {
1717 // N.B. m here is 1-relative.
1718 fp += 5;
1719 sscanf(fp, "%d", &m);
1720 if (m < 1 || M < m) {
1721 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1722 "Bad term (%d) for %s: %s", m, id, keyp->field);
1723 }
1724
1725 if ((fp = strchr(fp, '.')) == 0x0) {
1726 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1727 "Invalid field name for %s: %s", id, keyp->field);
1728 }
1729 fp++;
1730
1731 if (strcmp(fp, "COEFF") == 0) {
1732 i = iparm[I_DPOLY] + (m-1)*nTparm;
1733 dparm[i] = dpkeyd(keyp);
1734
1735 } else if (strncmp(fp, "VAR.", 4) == 0) {
1736 // N.B. jhat here is 1-relative.
1737 fp += 4;
1738 int jhat;
1739 sscanf(fp, "%d", &jhat);
1740 if (jhat < 1 || naxis < jhat) {
1741 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1742 "Invalid axis number (%d) for %s: %s", jhat, id, keyp->field);
1743 }
1744
1745 i = iparm[I_DPOLY] + (m-1)*nTparm + 1 + (jhat-1);
1746 double power = dpkeyd(keyp);
1747 dparm[i] = power;
1748
1749 } else if (strncmp(fp, "AUX.", 4) == 0) {
1750 // N.B. k here is 1-relative.
1751 fp += 4;
1752 sscanf(fp, "%d", &k);
1753 if (k < 1 || K < k) {
1754 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1755 "Bad auxiliary variable (%d) for %s: %s", k, id, keyp->field);
1756 }
1757
1758 i = iparm[I_DPOLY] + (m-1)*nTparm + 1 + Nhat + (k-1);
1759 double power = dpkeyd(keyp);
1760 dparm[i] = power;
1761
1762 } else {
1763 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1764 "Unrecognized field name for %s: %s", id, keyp->field);
1765 }
1766
1767 } else if (strcmp(fp, "DOCORR") &&
1768 strcmp(fp, "NAXES") &&
1769 strncmp(fp, "AXIS.", 5) &&
1770 strncmp(fp, "OFFSET.", 7) &&
1771 strncmp(fp, "SCALE.", 6) &&
1772 strcmp(fp, "NAUX") &&
1773 strcmp(fp, "NTERMS")) {
1774 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1775 "Unrecognized field name for %s: %s", id, keyp->field);
1776 }
1777 }
1778
1779
1780 // Optimization: when the power is integral, it is faster to multiply
1781 // ------------ repeatedly than call pow(). iparm[] is constructed as
1782 // follows:
1783 // I_NPOLY indexing parameters, as above,
1784 // nVar elements record the largest integral power for each variable,
1785 // nVar elements record offsets into dparm for each variable,
1786 // M*nVar flags to signal whether the power is integral,
1787 // M*nVar integral powers.
1788 for (int ivar = 0; ivar < nVar; ivar++) {
1789 // Want at least the first degree power for all variables.
1790 i = iparm[I_MAXPOW] + ivar;
1791 iparm[i] = 1;
1792 }
1793
1794 for (int ivar = 0; ivar < nVar; ivar++) {
1795 for (m = 0; m < M; m++) {
1796 i = iparm[I_DPOLY] + m*nTparm + 1 + ivar;
1797 double power = dparm[i];
1798
1799 // Is it integral? (Positive, negative, or zero.)
1800 int ipow = (int)power;
1801 if (power == (double)ipow) {
1802 // Signal that the power is integral.
1803 i = iparm[I_FLAGS] + m*nVar + ivar;
1804 if (ipow == 0) {
1805 iparm[i] = 3;
1806 } else {
1807 iparm[i] = 1;
1808 }
1809
1810 // The integral power itself.
1811 i = iparm[I_IPOW] + m*nVar + ivar;
1812 iparm[i] = ipow;
1813 }
1814
1815 // Record the largest integral power for each variable.
1816 i = iparm[I_MAXPOW] + ivar;
1817 if (iparm[i] < abs(ipow)) {
1818 iparm[i] = abs(ipow);
1819 }
1820 }
1821 }
1822
1823 // How many of all powers of each variable will there be?
1824 int npow = 0;
1825 for (int ivar = 0; ivar < nVar; ivar++) {
1826 // Offset into dparm.
1827 i = iparm[I_DPOFF] + ivar;
1828 iparm[i] = lendp + npow;
1829
1830 i = iparm[I_MAXPOW] + ivar;
1831 npow += iparm[i];
1832 }
1833
1834 // Expand dparm to store the extra powers.
1835 if (npow) {
1836 lendp += npow;
1837 iparm[I_LENDP] = lendp;
1838 if ((dis->dparm[j] = realloc(dparm, lendp*sizeof(double))) == 0x0) {
1839 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
1840 }
1841 }
1842
1843 // No specialist de-distortions.
1844 dis->disp2x[j] = dispoly;
1845 dis->disx2p[j] = 0x0;
1846
1847 // Translate Polynomial to TPD if possible, it's much faster.
1848 // However don't do it if the name was given as "Polynomial*".
1849 if (strcmp(dis->dtype[j], "Polynomial") == 0) {
1850 pol2tpd(j, dis);
1851 }
1852
1853 return 0;
1854 }
1855
1856 //----------------------------------------------------------------------------
1857
tpdset(int j,struct disprm * dis)1858 int tpdset(int j, struct disprm *dis)
1859
1860 {
1861 static const char *function = "tpdset";
1862
1863 if (dis == 0x0) return DISERR_NULL_POINTER;
1864 struct wcserr **err = &(dis->err);
1865
1866 char id[32];
1867 sprintf(id, "TPD on axis %d", j+1);
1868
1869
1870 // TPD distortion.
1871 if (dis->Nhat[j] < 1 || 2 < dis->Nhat[j]) {
1872 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1873 "Axis map for %s must contain 1 or 2 entries, not %d", id,
1874 dis->Nhat[j]);
1875 }
1876
1877 // Find the number of parameters.
1878 int ncoeff[2] = {0, 0};
1879 int doaux = 0;
1880 int doradial = 0;
1881 struct dpkey *keyp = dis->dp;
1882 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
1883 if (keyp->j-1 != j) continue;
1884
1885 char *fp = strchr(keyp->field, '.') + 1;
1886
1887 if (strncmp(fp, "TPD.", 4) == 0) {
1888 fp += 4;
1889 int idis;
1890 if (strncmp(fp, "FWD.", 4) == 0) {
1891 idis = 0;
1892
1893 } else if (strncmp(fp, "REV.", 4) == 0) {
1894 // TPD may provide a polynomial approximation for the inverse.
1895 idis = 1;
1896
1897 } else {
1898 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1899 "Unrecognized field name for %s: %s", id, keyp->field);
1900 }
1901
1902 int k;
1903 sscanf(fp+4, "%d", &k);
1904 if (0 <= k && k <= 59) {
1905 if (ncoeff[idis] < k+1) ncoeff[idis] = k+1;
1906
1907 // Any radial terms?
1908 if (k == 3 || k == 11 || k == 23 || k == 39 || k == 59) {
1909 doradial = 1;
1910 }
1911
1912 } else {
1913 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1914 "Invalid parameter number (%d) for %s: %s", k, id, keyp->field);
1915 }
1916
1917 } else if (strncmp(fp, "AUX.", 4) == 0) {
1918 // Flag usage of auxiliary variables.
1919 doaux = 1;
1920
1921 } else if (strcmp(fp, "DOCORR") &&
1922 strcmp(fp, "NAXES") &&
1923 strncmp(fp, "AXIS.", 5) &&
1924 strncmp(fp, "OFFSET.", 7) &&
1925 strncmp(fp, "SCALE.", 6)) {
1926 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1927 "Unrecognized field name for %s: %s", id, keyp->field);
1928 }
1929 }
1930
1931 int (*(distpd[2]))(DISP2X_ARGS) = {0x0, 0x0};
1932 for (int idis = 0; idis < 2; idis++) {
1933 if (ncoeff[idis] <= 4) {
1934 if (idis) {
1935 // No inverse polynomial.
1936 break;
1937 }
1938
1939 // First degree.
1940 ncoeff[idis] = 4;
1941 distpd[idis] = tpd1;
1942 } else if (ncoeff[idis] <= 7) {
1943 // Second degree.
1944 ncoeff[idis] = 7;
1945 distpd[idis] = tpd2;
1946 } else if (ncoeff[idis] <= 12) {
1947 // Third degree.
1948 ncoeff[idis] = 12;
1949 distpd[idis] = tpd3;
1950 } else if (ncoeff[idis] <= 17) {
1951 // Fourth degree.
1952 ncoeff[idis] = 17;
1953 distpd[idis] = tpd4;
1954 } else if (ncoeff[idis] <= 24) {
1955 // Fifth degree.
1956 ncoeff[idis] = 24;
1957 distpd[idis] = tpd5;
1958 } else if (ncoeff[idis] <= 31) {
1959 // Sixth degree.
1960 ncoeff[idis] = 31;
1961 distpd[idis] = tpd6;
1962 } else if (ncoeff[idis] <= 40) {
1963 // Seventh degree.
1964 ncoeff[idis] = 40;
1965 distpd[idis] = tpd7;
1966 } else if (ncoeff[idis] <= 49) {
1967 // Eighth degree.
1968 ncoeff[idis] = 49;
1969 distpd[idis] = tpd8;
1970 } else if (ncoeff[idis] <= 60) {
1971 // Ninth degree.
1972 ncoeff[idis] = 60;
1973 distpd[idis] = tpd9;
1974 } else {
1975 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
1976 "Invalid number of parameters (%d) for %s", ncoeff[idis], id);
1977 }
1978 }
1979
1980 // disx2p() only uses the inverse TPD, if present, to provide a better
1981 // zeroth approximation.
1982 dis->disp2x[j] = distpd[0];
1983 dis->disx2p[j] = distpd[1];
1984
1985
1986 // These iparm indices are specific to TPD (matching definitions in wcshdr.c).
1987 #define I_TPDNCO 3 // No. of TPD coefficients, forward...
1988 #define I_TPDINV 4 // ...and inverse.
1989 #define I_TPDAUX 5 // True if auxiliary variables are used.
1990 #define I_TPDRAD 6 // True if the radial variable is used.
1991 #define I_NTPD 7
1992
1993 // Record indexing parameters.
1994 int niparm = I_NTPD;
1995 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
1996 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
1997 }
1998
1999 int ndparm = (doaux?6:0) + ncoeff[0] + ncoeff[1];
2000
2001 // The first three are more widely used.
2002 dis->iparm[j][I_DTYPE] = DIS_TPD;
2003 dis->iparm[j][I_NIPARM] = niparm;
2004 dis->iparm[j][I_NDPARM] = ndparm;
2005
2006 // Number of TPD coefficients.
2007 dis->iparm[j][I_TPDNCO] = ncoeff[0];
2008 dis->iparm[j][I_TPDINV] = ncoeff[1];
2009
2010 // Flag for presence of auxiliary variables.
2011 dis->iparm[j][I_TPDAUX] = doaux;
2012
2013 // Flag for presence of radial terms.
2014 dis->iparm[j][I_TPDRAD] = doradial;
2015
2016
2017 // Allocate memory for the polynomial coefficients and fill it.
2018 if ((dis->dparm[j] = calloc(ndparm, sizeof(double))) == 0x0) {
2019 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2020 }
2021
2022 // Set default auxiliary coefficients.
2023 if (doaux) {
2024 dis->dparm[j][1] = 1.0;
2025 dis->dparm[j][5] = 1.0;
2026 }
2027
2028 keyp = dis->dp;
2029 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2030 if (keyp->j-1 != j) continue;
2031
2032 char *fp = strchr(keyp->field, '.') + 1;
2033
2034 if (strncmp(fp, "AUX.", 4) == 0) {
2035 // Auxiliary variables.
2036 fp += 4;
2037 int k;
2038 sscanf(fp, "%d", &k);
2039 if (k < 1 || 2 < k) {
2040 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2041 "Bad auxiliary variable (%d) for %s: %s", k, id, keyp->field);
2042 }
2043
2044 if ((fp = strchr(fp, '.')) == 0x0) {
2045 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2046 "Invalid field name for %s: %s", id, keyp->field);
2047 }
2048 fp++;
2049
2050 if (strncmp(fp, "COEFF.", 6) != 0) {
2051 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2052 "Unrecognized field name for %s: %s", id, keyp->field);
2053 }
2054
2055 fp += 6;
2056 int m;
2057 sscanf(fp, "%d", &m);
2058 if (m < 0 || 2 < m) {
2059 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2060 "Invalid coefficient number (%d) for %s: %s", m, id, keyp->field);
2061 }
2062
2063 int idis = 3*(k-1) + m;
2064 dis->dparm[j][idis] = dpkeyd(keyp);
2065
2066 } else if (strncmp(fp, "TPD.", 4) == 0) {
2067 fp += 4;
2068 int idis = (doaux?6:0);
2069 if (strncmp(fp, "REV.", 4) == 0) {
2070 idis += ncoeff[0];
2071 }
2072
2073 int k;
2074 sscanf(fp+4, "%d", &k);
2075 dis->dparm[j][idis+k] = dpkeyd(keyp);
2076 }
2077 }
2078
2079 return 0;
2080 }
2081
2082 //----------------------------------------------------------------------------
2083
pol2tpd(int j,struct disprm * dis)2084 int pol2tpd(int j, struct disprm *dis)
2085
2086 {
2087 static const char *function = "pol2tpd";
2088
2089 static const int map[][10] = {{ 0, 2, 6, 10, 16, 22, 30, 38, 48, 58},
2090 { 1, 5, 9, 15, 21, 29, 37, 47, 57, -1},
2091 { 4, 8, 14, 20, 28, 36, 46, 56, -1, -1},
2092 { 7, 13, 19, 27, 35, 45, 55, -1, -1, -1},
2093 {12, 18, 26, 34, 44, 54, -1, -1, -1, -1},
2094 {17, 25, 33, 43, 53, -1, -1, -1, -1, -1},
2095 {24, 32, 42, 52, -1, -1, -1, -1, -1, -1},
2096 {31, 41, 51, -1, -1, -1, -1, -1, -1, -1},
2097 {40, 50, -1, -1, -1, -1, -1, -1, -1, -1},
2098 {49, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
2099
2100 // Initialize.
2101 if (dis == 0x0) return DISERR_NULL_POINTER;
2102 struct wcserr **err = &(dis->err);
2103
2104 int *iparm = dis->iparm[j];
2105 double *dparm = dis->dparm[j];
2106
2107
2108 // Check the number of independent variables, no more than two.
2109 int Nhat = dis->Nhat[j];
2110 if (2 < Nhat) return -1;
2111
2112 // Check auxiliaries: only one is allowed...
2113 int K = iparm[I_K];
2114 if (1 < K) return -1;
2115 if (K) {
2116 // ...and it must be radial.
2117 if (dparm[0] != 0.0) return -1;
2118 if (dparm[1] != 1.0) return -1;
2119 if (dparm[2] != 1.0) return -1;
2120 if (dparm[3] != 0.5) return -1;
2121 if (dparm[4] != 2.0) return -1;
2122 if (dparm[5] != 2.0) return -1;
2123 }
2124
2125 // Check powers...
2126 int *iflgp = iparm + iparm[I_FLAGS];
2127 int *ipowp = iparm + iparm[I_IPOW];
2128 int degree = 0;
2129 for (int m = 0; m < iparm[I_M]; m++) {
2130 int deg = 0;
2131 for (int jhat = 0; jhat < Nhat; jhat++) {
2132 // ...they must be positive integral.
2133 if (*iflgp == 0) return -1;
2134 if (*ipowp < 0) return -1;
2135 deg += *ipowp;
2136 iflgp++;
2137 ipowp++;
2138 }
2139
2140 // The polynomial degree can't be greater than 9.
2141 if (9 < deg) return -1;
2142
2143 if (K) {
2144 // Likewise for the radial variable.
2145 if (*iflgp == 0) return -1;
2146 if (*ipowp) {
2147 if (*ipowp < 0) return -1;
2148 if (9 < *ipowp) return -1;
2149
2150 // Can't mix the radial and other terms.
2151 if (deg) return -1;
2152
2153 // Can't have even powers of the radial variable.
2154 deg = *ipowp;
2155 if (!(deg%2)) return -1;
2156 }
2157 iflgp++;
2158 ipowp++;
2159 }
2160
2161 if (degree < deg) degree = deg;
2162 }
2163
2164
2165 // OK, it ticks all the boxes. Now translate it.
2166 int ndparm = 0;
2167 if (degree == 1) {
2168 ndparm = 4;
2169 dis->disp2x[j] = tpd1;
2170 } else if (degree == 2) {
2171 ndparm = 7;
2172 dis->disp2x[j] = tpd2;
2173 } else if (degree == 3) {
2174 ndparm = 12;
2175 dis->disp2x[j] = tpd3;
2176 } else if (degree == 4) {
2177 ndparm = 17;
2178 dis->disp2x[j] = tpd4;
2179 } else if (degree == 5) {
2180 ndparm = 24;
2181 dis->disp2x[j] = tpd5;
2182 } else if (degree == 6) {
2183 ndparm = 31;
2184 dis->disp2x[j] = tpd6;
2185 } else if (degree == 7) {
2186 ndparm = 40;
2187 dis->disp2x[j] = tpd7;
2188 } else if (degree == 8) {
2189 ndparm = 49;
2190 dis->disp2x[j] = tpd8;
2191 } else if (degree == 9) {
2192 ndparm = 60;
2193 dis->disp2x[j] = tpd9;
2194 }
2195
2196 // No specialist de-distortions.
2197 dis->disx2p[j] = 0x0;
2198
2199 // Record indexing parameters.
2200 int niparm = I_NTPD;
2201 int *tpd_iparm;
2202 if ((tpd_iparm = calloc(niparm, sizeof(int))) == 0x0) {
2203 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2204 }
2205
2206 // The first three are more widely used.
2207 tpd_iparm[I_DTYPE] = DIS_TPD;
2208 tpd_iparm[I_NIPARM] = niparm;
2209 tpd_iparm[I_NDPARM] = ndparm;
2210
2211 // Number of TPD coefficients.
2212 tpd_iparm[I_TPDNCO] = ndparm;
2213 tpd_iparm[I_TPDINV] = 0;
2214
2215 // No auxiliary variables yet.
2216 tpd_iparm[I_TPDAUX] = 0;
2217
2218 // Flag for presence of radial terms.
2219 tpd_iparm[I_TPDRAD] = K;
2220
2221
2222 // Allocate memory for the polynomial coefficients and fill it.
2223 double *tpd_dparm;
2224 if ((tpd_dparm = calloc(ndparm, sizeof(double))) == 0x0) {
2225 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2226 }
2227
2228 ipowp = iparm + iparm[I_IPOW];
2229 double *dpolp = dparm + iparm[I_DPOLY];
2230 for (int m = 0; m < iparm[I_M]; m++) {
2231 if (K && ipowp[Nhat]) {
2232 // The radial variable.
2233 switch (ipowp[Nhat]) {
2234 case 1:
2235 tpd_dparm[3] = *dpolp;
2236 break;
2237 case 3:
2238 tpd_dparm[11] = *dpolp;
2239 break;
2240 case 5:
2241 tpd_dparm[23] = *dpolp;
2242 break;
2243 case 7:
2244 tpd_dparm[39] = *dpolp;
2245 break;
2246 case 9:
2247 tpd_dparm[59] = *dpolp;
2248 break;
2249 }
2250
2251 } else {
2252 // The independent variables.
2253 int p[] = {0, 0};
2254 for (int jhat = 0; jhat < Nhat; jhat++) {
2255 p[jhat] = ipowp[jhat];
2256 }
2257
2258 int n = map[p[0]][p[1]];
2259 tpd_dparm[n] = *dpolp;
2260 }
2261
2262
2263 ipowp += iparm[I_NVAR];
2264 dpolp += iparm[I_NVAR] + 1;
2265 }
2266
2267
2268 // Switch from Polynomial to TPD.
2269 free(iparm);
2270 free(dparm);
2271 dis->iparm[j] = tpd_iparm;
2272 dis->dparm[j] = tpd_dparm;
2273
2274 return 0;
2275 }
2276
2277 //----------------------------------------------------------------------------
2278
tpvset(int j,struct disprm * dis)2279 int tpvset(int j, struct disprm *dis)
2280
2281 {
2282 static const char *function = "tpvset";
2283
2284 // Initialize.
2285 if (dis == 0x0) return DISERR_NULL_POINTER;
2286 struct wcserr **err = &(dis->err);
2287
2288 // TPV "projection".
2289 char id[32];
2290 sprintf(id, "TPV on axis %d", j+1);
2291
2292 // TPV is a sequent distortion, applied to intermediate world coordinates
2293 // (normally used with CDi_ja). It computes corrected coordinates directly.
2294 dis->docorr[j] = 0;
2295
2296 if (dis->Nhat[j] != 2) {
2297 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2298 "Axis map for %s must contain 2 entries, not %d", id, dis->Nhat[j]);
2299 }
2300
2301 // Find the number of parameters.
2302 int ndparm = 0;
2303 int doradial = 0;
2304 struct dpkey *keyp = dis->dp;
2305 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2306 if (keyp->j-1 != j) continue;
2307
2308 char *fp = strchr(keyp->field, '.') + 1;
2309
2310 if (strncmp(fp, "TPV.", 4) == 0) {
2311 int k;
2312 sscanf(fp+4, "%d", &k);
2313 if (0 <= k && k <= 39) {
2314 if (ndparm < k+1) ndparm = k+1;
2315
2316 // Any radial terms?
2317 if (k == 3 || k == 11 || k == 23 || k == 39 || k == 59) {
2318 doradial = 1;
2319 }
2320
2321 } else {
2322 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2323 "Invalid parameter number (%d) for %s: %s", k, id, keyp->field);
2324 }
2325
2326 } else if (strcmp(fp, "NAXES") &&
2327 strncmp(fp, "AXIS.", 5) &&
2328 strncmp(fp, "OFFSET.", 7) &&
2329 strncmp(fp, "SCALE.", 6)) {
2330 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2331 "Unrecognized field name for %s: %s", id, keyp->field);
2332 }
2333 }
2334
2335 // TPD is going to do the dirty work.
2336 if (ndparm <= 4) {
2337 // First degree.
2338 ndparm = 4;
2339 dis->disp2x[j] = tpd1;
2340 } else if (ndparm <= 7) {
2341 // Second degree.
2342 ndparm = 7;
2343 dis->disp2x[j] = tpd2;
2344 } else if (ndparm <= 12) {
2345 // Third degree.
2346 ndparm = 12;
2347 dis->disp2x[j] = tpd3;
2348 } else if (ndparm <= 17) {
2349 // Fourth degree.
2350 ndparm = 17;
2351 dis->disp2x[j] = tpd4;
2352 } else if (ndparm <= 24) {
2353 // Fifth degree.
2354 ndparm = 24;
2355 dis->disp2x[j] = tpd5;
2356 } else if (ndparm <= 31) {
2357 // Sixth degree.
2358 ndparm = 31;
2359 dis->disp2x[j] = tpd6;
2360 } else if (ndparm <= 40) {
2361 // Seventh degree.
2362 ndparm = 40;
2363 dis->disp2x[j] = tpd7;
2364 } else {
2365 // Could go to ninth degree, but that wouldn't be legit.
2366 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2367 "Invalid number of parameters (%d) for %s", ndparm, id);
2368 }
2369
2370 // No specialist de-distortions.
2371 dis->disx2p[j] = 0x0;
2372
2373 // Record indexing parameters.
2374 int niparm = I_NTPD;
2375 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
2376 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2377 }
2378
2379 // The first three are more widely used.
2380 dis->iparm[j][I_DTYPE] = DIS_TPD;
2381 dis->iparm[j][I_NIPARM] = niparm;
2382 dis->iparm[j][I_NDPARM] = ndparm;
2383
2384 // Number of TPD coefficients.
2385 dis->iparm[j][I_TPDNCO] = ndparm;
2386 dis->iparm[j][I_TPDINV] = 0;
2387
2388 // TPV never needs auxiliary variables.
2389 dis->iparm[j][I_TPDAUX] = 0;
2390
2391 // Flag for presence of radial terms.
2392 dis->iparm[j][I_TPDRAD] = doradial;
2393
2394
2395 // Allocate memory for the polynomial coefficients and fill it.
2396 if ((dis->dparm[j] = calloc(ndparm, sizeof(double))) == 0x0) {
2397 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2398 }
2399
2400 keyp = dis->dp;
2401 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2402 if (keyp->j-1 != j) continue;
2403
2404 char *fp = strchr(keyp->field, '.') + 1;
2405
2406 // One-to-one correspondence between TPV and TPD coefficients.
2407 if (strncmp(fp, "TPV.", 4) == 0) {
2408 int k;
2409 sscanf(fp+4, "%d", &k);
2410 dis->dparm[j][k] = dpkeyd(keyp);
2411 }
2412 }
2413
2414 return 0;
2415 }
2416
2417 //----------------------------------------------------------------------------
2418
sipset(int j,struct disprm * dis)2419 int sipset(int j, struct disprm *dis)
2420
2421 {
2422 static const char *function = "sipset";
2423
2424 static const int map[][10] = {{ 0, 2, 6, 10, 16, 22, 30, 38, 48, 58},
2425 { 1, 5, 9, 15, 21, 29, 37, 47, 57, -1},
2426 { 4, 8, 14, 20, 28, 36, 46, 56, -1, -1},
2427 { 7, 13, 19, 27, 35, 45, 55, -1, -1, -1},
2428 {12, 18, 26, 34, 44, 54, -1, -1, -1, -1},
2429 {17, 25, 33, 43, 53, -1, -1, -1, -1, -1},
2430 {24, 32, 42, 52, -1, -1, -1, -1, -1, -1},
2431 {31, 41, 51, -1, -1, -1, -1, -1, -1, -1},
2432 {40, 50, -1, -1, -1, -1, -1, -1, -1, -1},
2433 {49, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
2434
2435 // Initialize.
2436 if (dis == 0x0) return DISERR_NULL_POINTER;
2437 struct wcserr **err = &(dis->err);
2438
2439 // Simple Imaging Polynomial.
2440 char id[32];
2441 sprintf(id, "SIP on axis %d", j+1);
2442
2443
2444 // SIP is a prior distortion that computes an additive correction.
2445 dis->docorr[j] = 1;
2446
2447 if (dis->Nhat[j] != 2) {
2448 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2449 "Axis map for %s must contain 2 entries, not %d", id, dis->Nhat[j]);
2450 }
2451
2452 // Find the polynomial degree, at least 1 for the forward function.
2453 int degree[2] = {1, -1};
2454 struct dpkey *keyp = dis->dp;
2455 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2456 if (keyp->j-1 != j) continue;
2457
2458 char *fp = strchr(keyp->field, '.') + 1;
2459
2460 if (strncmp(fp, "SIP.", 4) == 0) {
2461 fp += 4;
2462 int idis;
2463 if (strncmp(fp, "FWD.", 4) == 0) {
2464 idis = 0;
2465
2466 } else if (strncmp(fp, "REV.", 4) == 0) {
2467 // SIP uses a polynomial approximation for the inverse.
2468 idis = 1;
2469
2470 } else {
2471 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2472 "Unrecognized field name for %s: %s", id, keyp->field);
2473 }
2474
2475 fp += 4;
2476 int p, q;
2477 sscanf(fp, "%d_%d", &p, &q);
2478 int deg = p + q;
2479 if (p < 0 || 9 < p || q < 0 || 9 < q || 9 < deg) {
2480 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2481 "Invalid powers (%d, %d) for %s: %s", p, q, id, keyp->field);
2482 }
2483
2484 if (degree[idis] < deg) degree[idis] = deg;
2485
2486 } else if (strcmp(fp, "NAXES") &&
2487 strncmp(fp, "AXIS.", 5) &&
2488 strncmp(fp, "OFFSET.", 7) &&
2489 strncmp(fp, "SCALE.", 6)) {
2490 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2491 "Unrecognized field name for %s: %s", id, keyp->field);
2492 }
2493 }
2494
2495 if (degree[1] == 0 ) degree[1] = 1;
2496
2497 // TPD is going to do the dirty work.
2498 int (*(distpd[2]))(DISP2X_ARGS) = {0x0, 0x0}, ncoeff[2];
2499 for (int idis = 0; idis < 2; idis++) {
2500 ncoeff[idis] = 0;
2501 if (degree[idis] == 1) {
2502 ncoeff[idis] = 4;
2503 distpd[idis] = tpd1;
2504 } else if (degree[idis] == 2) {
2505 ncoeff[idis] = 7;
2506 distpd[idis] = tpd2;
2507 } else if (degree[idis] == 3) {
2508 ncoeff[idis] = 12;
2509 distpd[idis] = tpd3;
2510 } else if (degree[idis] == 4) {
2511 ncoeff[idis] = 17;
2512 distpd[idis] = tpd4;
2513 } else if (degree[idis] == 5) {
2514 ncoeff[idis] = 24;
2515 distpd[idis] = tpd5;
2516 } else if (degree[idis] == 6) {
2517 ncoeff[idis] = 31;
2518 distpd[idis] = tpd6;
2519 } else if (degree[idis] == 7) {
2520 ncoeff[idis] = 40;
2521 distpd[idis] = tpd7;
2522 } else if (degree[idis] == 8) {
2523 ncoeff[idis] = 49;
2524 distpd[idis] = tpd8;
2525 } else if (degree[idis] == 9) {
2526 ncoeff[idis] = 60;
2527 distpd[idis] = tpd9;
2528 }
2529 }
2530
2531 // SIP uses a polynomial approximation to the inverse. It's not very
2532 // accurate but may provide disx2p() with a better zeroth approximation.
2533 dis->disp2x[j] = distpd[0];
2534 dis->disx2p[j] = distpd[1];
2535
2536
2537 // Record indexing parameters.
2538 int niparm = I_NTPD;
2539 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
2540 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2541 }
2542
2543 int ndparm = ncoeff[0] + ncoeff[1];
2544
2545 // The first three are more widely used.
2546 dis->iparm[j][I_DTYPE] = DIS_TPD;
2547 dis->iparm[j][I_NIPARM] = niparm;
2548 dis->iparm[j][I_NDPARM] = ndparm;
2549
2550 // Number of TPD coefficients.
2551 dis->iparm[j][I_TPDNCO] = ncoeff[0];
2552 dis->iparm[j][I_TPDINV] = ncoeff[1];
2553
2554 // SIP never needs auxiliary variables.
2555 dis->iparm[j][I_TPDAUX] = 0;
2556
2557 // SIP never needs the radial terms.
2558 dis->iparm[j][I_TPDRAD] = 0;
2559
2560
2561 // Allocate memory for the polynomial coefficients and fill it.
2562 if ((dis->dparm[j] = calloc(ndparm, sizeof(double))) == 0x0) {
2563 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2564 }
2565
2566 keyp = dis->dp;
2567 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2568 if (keyp->j-1 != j) continue;
2569
2570 char *fp = strchr(keyp->field, '.') + 1;
2571
2572 if (strncmp(fp, "SIP.", 4) == 0) {
2573 fp += 4;
2574 int idis;
2575 if (strncmp(fp, "FWD.", 4) == 0) {
2576 idis = 0;
2577 } else {
2578 idis = ncoeff[0];
2579 }
2580
2581 int p, q;
2582 sscanf(fp+4, "%d_%d", &p, &q);
2583
2584 // Map to TPD coefficient number.
2585 idis += map[p][q];
2586
2587 dis->dparm[j][idis] = dpkeyd(keyp);
2588 }
2589 }
2590
2591
2592 return 0;
2593 }
2594
2595 //----------------------------------------------------------------------------
2596
dssset(int j,struct disprm * dis)2597 int dssset(int j, struct disprm *dis)
2598
2599 {
2600 static const char *function = "dssset";
2601
2602 // Initialize.
2603 if (dis == 0x0) return DISERR_NULL_POINTER;
2604 struct wcserr **err = &(dis->err);
2605
2606 // Digitized Sky Survey.
2607 char id[32];
2608 sprintf(id, "DSS on axis %d", j+1);
2609
2610
2611 // DSS is translated into a sequent distortion, applied to intermediate
2612 // pixel coordinates. It computes corrected coordinates directly.
2613 dis->docorr[j] = 0;
2614
2615 if (dis->Nhat[j] != 2) {
2616 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2617 "Axis map for %s must contain 2 entries, not %d", id, dis->Nhat[j]);
2618 }
2619
2620 // Safe to assume the polynomial degree is 5 (or less).
2621 int ncoeff = 24;
2622 dis->disp2x[j] = tpd5;
2623
2624 // No specialist de-distortions.
2625 dis->disx2p[j] = 0x0;
2626
2627
2628 // Record indexing parameters.
2629 int niparm = I_NTPD;
2630 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
2631 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2632 }
2633
2634 int ndparm = 6 + ncoeff;
2635
2636 // The first three are more widely used.
2637 dis->iparm[j][I_DTYPE] = DIS_TPD;
2638 dis->iparm[j][I_NIPARM] = niparm;
2639 dis->iparm[j][I_NDPARM] = ndparm;
2640
2641 // Number of TPD coefficients.
2642 dis->iparm[j][I_TPDNCO] = ncoeff;
2643 dis->iparm[j][I_TPDINV] = 0;
2644
2645 // DSS always needs auxiliary variables.
2646 dis->iparm[j][I_TPDAUX] = 1;
2647
2648 // DSS never needs the radial terms.
2649 dis->iparm[j][I_TPDRAD] = 0;
2650
2651
2652 // Allocate memory for the polynomial coefficients and fill it.
2653 if ((dis->dparm[j] = calloc(ndparm, sizeof(double))) == 0x0) {
2654 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2655 }
2656
2657 // This translation follows WCS Paper IV, Sect. 5.2 using the same
2658 // variable names. Find A1, A2, A3, B1, B2, and B3.
2659 double A1, A2, A3, B1, B2, B3;
2660 A1 = A2 = A3 = 0.0;
2661 B1 = B2 = B3 = 0.0;
2662 struct dpkey *keyp = dis->dp;
2663 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2664 char *fp = strchr(keyp->field, '.') + 1;
2665 if (strncmp(fp, "DSS.AMD.", 8) == 0) {
2666 fp += 8;
2667 int m;
2668 sscanf(fp, "%d", &m);
2669
2670 if (m == 1) {
2671 if (keyp->j == 1) {
2672 A1 = dpkeyd(keyp);
2673 } else {
2674 B1 = dpkeyd(keyp);
2675 }
2676 } else if (m == 2) {
2677 if (keyp->j == 1) {
2678 A2 = dpkeyd(keyp);
2679 } else {
2680 B2 = dpkeyd(keyp);
2681 }
2682 } else if (m == 3) {
2683 if (keyp->j == 1) {
2684 A3 = dpkeyd(keyp);
2685 } else {
2686 B3 = dpkeyd(keyp);
2687 }
2688 }
2689 }
2690 }
2691
2692 double X0 = (A2*B3 - A3*B1) / (A1*B1 - A2*B2);
2693 double Y0 = (A3*B2 - A1*B3) / (A1*B1 - A2*B2);
2694
2695 double S = sqrt(fabs(A1*B1 - A2*B2));
2696 if (S == 0.0) {
2697 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2698 "Coefficient scale for %s is zero.", id);
2699 }
2700
2701 // Coefficients for the auxiliary variables.
2702 double *dparm = dis->dparm[j];
2703 if (j == 0) {
2704 dparm[0] = X0;
2705 dparm[1] = -B1/S;
2706 dparm[2] = -A2/S;
2707 dparm[3] = Y0;
2708 dparm[4] = B2/S;
2709 dparm[5] = A1/S;
2710
2711 // Change the sign of S for scaling the A coefficients.
2712 S *= -1.0;
2713
2714 } else {
2715 dparm[0] = Y0;
2716 dparm[1] = B2/S;
2717 dparm[2] = A1/S;
2718 dparm[3] = X0;
2719 dparm[4] = -B1/S;
2720 dparm[5] = -A2/S;
2721 }
2722
2723 // Translate DSS coefficients to TPD.
2724 dparm += 6;
2725 int degree = 3;
2726 keyp = dis->dp;
2727 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2728 if (keyp->j-1 != j) continue;
2729
2730 char *fp = strchr(keyp->field, '.') + 1;
2731
2732 if (strncmp(fp, "DSS.AMD.", 8) == 0) {
2733 // Skip zero coefficients.
2734 double coeff = dpkeyd(keyp);
2735 if (coeff == 0.0) continue;
2736
2737 fp += 8;
2738 int m;
2739 sscanf(fp, "%d", &m);
2740
2741 // Apply the coefficient scale factor.
2742 coeff /= S;
2743
2744 if (m == 1) {
2745 dparm[1] = coeff;
2746 } else if (m == 2) {
2747 dparm[2] = coeff;
2748 } else if (m == 3) {
2749 dparm[0] = coeff;
2750 } else if (m == 4) {
2751 dparm[4] += coeff;
2752 } else if (m == 5) {
2753 dparm[5] = coeff;
2754 } else if (m == 6) {
2755 dparm[6] += coeff;
2756 } else if (m == 7) {
2757 dparm[4] += coeff;
2758 dparm[6] += coeff;
2759 } else if (m == 8) {
2760 dparm[7] += coeff;
2761 } else if (m == 9) {
2762 dparm[8] = coeff;
2763 } else if (m == 10) {
2764 dparm[9] += coeff;
2765 } else if (m == 11) {
2766 dparm[10] = coeff;
2767 } else if (m == 12) {
2768 dparm[7] += coeff;
2769 dparm[9] += coeff;
2770 } else if (m == 13) {
2771 dparm[17] = coeff;
2772 dparm[19] = coeff * 2.0;
2773 dparm[21] = coeff;
2774 degree = 5;
2775 } else if (coeff != 0.0) {
2776 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2777 "Invalid parameter for %s: %s", m, id, keyp->field);
2778 }
2779
2780 } else if (strcmp(fp, "NAXES") &&
2781 strncmp(fp, "AXIS.", 5) &&
2782 strncmp(fp, "OFFSET.", 7) &&
2783 strncmp(fp, "SCALE.", 6)) {
2784 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2785 "Unrecognized field name for %s: %s", id, keyp->field);
2786 }
2787 }
2788
2789 // The DSS polynomial doesn't have 4th degree terms, and the 5th degree
2790 // coefficient is often zero.
2791 if (degree == 3) {
2792 dis->iparm[j][I_TPDNCO] = 12;
2793 dis->disp2x[j] = tpd3;
2794 }
2795
2796 return 0;
2797 }
2798
2799 //----------------------------------------------------------------------------
2800
2801 #define CHEBYSHEV 1
2802 #define LEGENDRE 2
2803 #define MONOMIAL 3
2804
watset(int j,struct disprm * dis)2805 int watset(int j, struct disprm *dis)
2806
2807 {
2808 static const char *function = "watset";
2809
2810 static const int map[][10] = {{ 0, 2, 6, 10, 16, 22, 30, 38, 48, 58},
2811 { 1, 5, 9, 15, 21, 29, 37, 47, 57, -1},
2812 { 4, 8, 14, 20, 28, 36, 46, 56, -1, -1},
2813 { 7, 13, 19, 27, 35, 45, 55, -1, -1, -1},
2814 {12, 18, 26, 34, 44, 54, -1, -1, -1, -1},
2815 {17, 25, 33, 43, 53, -1, -1, -1, -1, -1},
2816 {24, 32, 42, 52, -1, -1, -1, -1, -1, -1},
2817 {31, 41, 51, -1, -1, -1, -1, -1, -1, -1},
2818 {40, 50, -1, -1, -1, -1, -1, -1, -1, -1},
2819 {49, -1, -1, -1, -1, -1, -1, -1, -1, -1}};
2820
2821 // Initialize.
2822 if (dis == 0x0) return DISERR_NULL_POINTER;
2823 struct wcserr **err = &(dis->err);
2824
2825 // WAT (TNX or ZPX) Polynomial.
2826 char id[32];
2827 sprintf(id, "WAT (%s) on axis %d", dis->dtype[0]+4, j+1);
2828
2829
2830 // WAT is a sequent distortion, applied to intermediate world coordinates
2831 // (normally used with CDi_ja). It computes an additive correction.
2832 dis->docorr[j] = 1;
2833
2834 if (dis->Nhat[j] != 2) {
2835 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2836 "Axis map for %s must contain 2 entries, not %d", id, dis->Nhat[j]);
2837 }
2838
2839 // Find the polynomial degree (at least 1), kind, and domain.
2840 int degree = 1;
2841 int kind = 0;
2842 double xmin = 0.0;
2843 double xmax = 0.0;
2844 double ymin = 0.0;
2845 double ymax = 0.0;
2846 struct dpkey *keyp = dis->dp;
2847 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
2848 if (keyp->j-1 != j) continue;
2849
2850 char *fp = strchr(keyp->field, '.') + 1;
2851
2852 if (strncmp(fp, "WAT.", 4) == 0) {
2853 fp += 4;
2854 if (strncmp(fp, "CHBY.", 5) == 0 ||
2855 strncmp(fp, "LEGR.", 5) == 0 ||
2856 strncmp(fp, "MONO.", 5) == 0) {
2857
2858 fp += 5;
2859 int m, n;
2860 sscanf(fp, "%d_%d", &m, &n);
2861 int deg = m + n;
2862 if (m < 0 || 9 < m || n < 0 || 9 < n || 9 < deg) {
2863 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2864 "Invalid powers (%d, %d) for %s: %s", m, n, id, keyp->field);
2865 }
2866
2867 if (degree < deg) degree = deg;
2868
2869 } else if (strcmp(fp, "POLY") == 0) {
2870 kind = dpkeyi(keyp);
2871
2872 } else if (strcmp(fp, "XMIN") == 0) {
2873 xmin = dpkeyd(keyp);
2874
2875 } else if (strcmp(fp, "XMAX") == 0) {
2876 xmax = dpkeyd(keyp);
2877
2878 } else if (strcmp(fp, "YMIN") == 0) {
2879 ymin = dpkeyd(keyp);
2880
2881 } else if (strcmp(fp, "YMAX") == 0) {
2882 ymax = dpkeyd(keyp);
2883 }
2884
2885 } else if (strcmp(fp, "NAXES") &&
2886 strncmp(fp, "AXIS.", 5) &&
2887 strncmp(fp, "OFFSET.", 7) &&
2888 strncmp(fp, "SCALE.", 6)) {
2889 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2890 "Unrecognized field name for %s: %s", id, keyp->field);
2891 }
2892 }
2893
2894 int doaux = (kind == 1 || kind == 2);
2895
2896 // TPD is going to do the dirty work.
2897 int ncoeff = 0;
2898 if (degree == 1) {
2899 // First degree.
2900 ncoeff = 4;
2901 dis->disp2x[j] = tpd1;
2902 } else if (degree == 2) {
2903 // Second degree.
2904 ncoeff = 7;
2905 dis->disp2x[j] = tpd2;
2906 } else if (degree == 3) {
2907 // Third degree.
2908 ncoeff = 12;
2909 dis->disp2x[j] = tpd3;
2910 } else if (degree == 4) {
2911 // Fourth degree.
2912 ncoeff = 17;
2913 dis->disp2x[j] = tpd4;
2914 } else if (degree == 5) {
2915 // Fifth degree.
2916 ncoeff = 24;
2917 dis->disp2x[j] = tpd5;
2918 } else if (degree == 6) {
2919 // Sixth degree.
2920 ncoeff = 31;
2921 dis->disp2x[j] = tpd6;
2922 } else if (degree == 7) {
2923 // Seventh degree.
2924 ncoeff = 40;
2925 dis->disp2x[j] = tpd7;
2926 } else if (degree == 8) {
2927 // Eighth degree.
2928 ncoeff = 49;
2929 dis->disp2x[j] = tpd8;
2930 } else if (degree == 9) {
2931 // Ninth degree.
2932 ncoeff = 60;
2933 dis->disp2x[j] = tpd9;
2934 }
2935
2936 // No specialist de-distortions.
2937 dis->disx2p[j] = 0x0;
2938
2939
2940 // Record indexing parameters.
2941 int niparm = I_NTPD;
2942 if ((dis->iparm[j] = calloc(niparm, sizeof(int))) == 0x0) {
2943 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2944 }
2945
2946 int *iparm = dis->iparm[j];
2947
2948 int ndparm = 6 + ncoeff;
2949
2950 // The first three are more widely used.
2951 iparm[I_DTYPE] = DIS_TPD;
2952 iparm[I_NIPARM] = niparm;
2953 iparm[I_NDPARM] = ndparm;
2954
2955 // Number of TPD coefficients.
2956 iparm[I_TPDNCO] = ncoeff;
2957 iparm[I_TPDINV] = 0;
2958
2959 // The Chebyshev and Legendre polynomials use auxiliary variables.
2960 iparm[I_TPDAUX] = doaux;
2961
2962 // WAT never needs the radial terms.
2963 iparm[I_TPDRAD] = 0;
2964
2965
2966 // Allocate memory for the polynomial coefficients and fill it.
2967 if ((dis->dparm[j] = calloc(ndparm, sizeof(double))) == 0x0) {
2968 return wcserr_set(DIS_ERRMSG(DISERR_MEMORY));
2969 }
2970
2971 double *dparm = dis->dparm[j];
2972
2973
2974 // Coefficients for the auxiliary variables.
2975 if (doaux) {
2976 double x0 = (xmax + xmin)/2.0;
2977 double dx = (xmax - xmin)/2.0;
2978 if (dx == 0.0) {
2979 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2980 "X-span for %s is zero", id);
2981 }
2982
2983 dparm[0] = -x0/dx;
2984 dparm[1] = 1.0/dx;
2985 dparm[2] = 0.0;
2986
2987 double y0 = (ymax + ymin)/2.0;
2988 double dy = (ymax - ymin)/2.0;
2989 if (dy == 0.0) {
2990 return wcserr_set(WCSERR_SET(DISERR_BAD_PARAM),
2991 "Y-span for %s is zero", id);
2992 }
2993
2994 dparm[3] = -y0/dy;
2995 dparm[4] = 0.0;
2996 dparm[5] = 1.0/dy;
2997
2998 dparm += 6;
2999 }
3000
3001
3002 // Unpack the polynomial coefficients.
3003 keyp = dis->dp;
3004 for (int idp = 0; idp < dis->ndp; idp++, keyp++) {
3005 if (keyp->j-1 != j) continue;
3006
3007 char *fp = strchr(keyp->field, '.') + 1;
3008
3009 if ((kind == CHEBYSHEV && strncmp(fp, "WAT.CHBY.", 9) == 0) ||
3010 (kind == LEGENDRE && strncmp(fp, "WAT.LEGR.", 9) == 0) ||
3011 (kind == MONOMIAL && strncmp(fp, "WAT.MONO.", 9) == 0)) {
3012 fp += 9;
3013
3014 int m, n;
3015 sscanf(fp, "%d_%d", &m, &n);
3016
3017 if (kind == MONOMIAL) {
3018 // Monomial coefficient, maps simply to TPD coefficient number.
3019 int idis = map[m][n];
3020 dparm[idis] = dpkeyd(keyp);
3021
3022 } else {
3023 // Coefficient of the product of two Chebyshev or two Legendre
3024 // polynomials. Find the corresponding monomial coefficients.
3025 double coeff = dpkeyd(keyp);
3026
3027 double coeffm[10], coeffn[10];
3028 cheleg(kind, m, n, coeffm, coeffn);
3029 for (int im = 0; im <= m; im++) {
3030 if (coeffm[im] == 0.0) continue;
3031
3032 for (int in = 0; in <= n; in++) {
3033 if (coeffn[in] == 0.0) continue;
3034
3035 int idis = map[im][in];
3036 dparm[idis] += coeff*coeffm[im]*coeffn[in];
3037 }
3038 }
3039 }
3040 }
3041 }
3042
3043 return 0;
3044 }
3045
3046 //----------------------------------------------------------------------------
3047 // Compute the coefficients of Chebyshev or Legendre polynomials of degree
3048 // m and n.
3049
cheleg(int kind,int m,int n,double coeffm[],double coeffn[])3050 int cheleg(int kind, int m, int n, double coeffm[], double coeffn[])
3051
3052 {
3053 int N = (m > n) ? m : n;
3054
3055 // Allocate work arrays.
3056 double *coeff[3];
3057 coeff[0] = calloc(3*(N+1), sizeof(double));
3058 coeff[1] = coeff[0] + (N+1);
3059 coeff[2] = coeff[1] + (N+1);
3060
3061 for (int j = 0; j <= N; j++) {
3062 int j0 = j%3;
3063
3064 if (j == 0) {
3065 coeff[0][0] = 1.0;
3066
3067 } else if (j == 1) {
3068 coeff[1][1] = 1.0;
3069
3070 } else {
3071 // Cyclic buffer indices.
3072 int j1 = (j-1)%3;
3073 int j2 = (j-2)%3;
3074
3075 memset(coeff[j0], 0, (N+1)*sizeof(double));
3076
3077 double d = (double)j;
3078 for (int k = 0; k < N; k++) {
3079 if (kind == CHEBYSHEV) {
3080 coeff[j0][k+1] = 2.0 * coeff[j1][k];
3081 coeff[j0][k] -= coeff[j2][k];
3082 } else if (kind == LEGENDRE) {
3083 coeff[j0][k+1] = ((2.0*d - 1.0) * coeff[j1][k]) / d;
3084 coeff[j0][k] -= ((d - 1.0) * coeff[j2][k]) / d;
3085 }
3086 }
3087 }
3088
3089 if (j == m) memcpy(coeffm, coeff[j0], (m+1)*sizeof(double));
3090 if (j == n) memcpy(coeffn, coeff[j0], (n+1)*sizeof(double));
3091 }
3092
3093 free(coeff[0]);
3094
3095 return 0;
3096 }
3097
3098 //----------------------------------------------------------------------------
3099
dispoly(int dummy,const int iparm[],const double dparm[],int Nhat,const double rawcrd[],double * discrd)3100 int dispoly(
3101 int dummy,
3102 const int iparm[],
3103 const double dparm[],
3104 int Nhat,
3105 const double rawcrd[],
3106 double *discrd)
3107
3108 {
3109 // Avert nuisance compiler warnings about unused parameters.
3110 (void)dummy;
3111
3112 // Check for zeroes.
3113 for (int jhat = 0; jhat < Nhat; jhat++) {
3114 if (rawcrd[jhat] == 0.0) {
3115 *discrd = 0.0;
3116 return 0;
3117 }
3118 }
3119
3120 // Working memory for auxiliaries &c. was allocated at the end of p[].
3121 double *aux = (double *)(dparm + iparm[I_DAUX]);
3122
3123 // Compute the auxiliary variables.
3124 for (int k = 0; k < iparm[I_K]; k++) {
3125 const double *cptr = dparm + k*iparm[I_NKPARM];
3126 const double *pptr = cptr + (1+Nhat);
3127
3128 aux[k] = *(cptr++);
3129 double auxp0 = *(pptr++);
3130
3131 for (int jhat = 0; jhat < Nhat; jhat++) {
3132 aux[k] += *(cptr++)*pow(rawcrd[jhat], *(pptr++));
3133 }
3134
3135 aux[k] = pow(aux[k], auxp0);
3136
3137 // Check for zeroes.
3138 if (aux[k] == 0.0) {
3139 *discrd = 0.0;
3140 return 0;
3141 }
3142 }
3143
3144
3145 // Compute all required integral powers of the variables.
3146 const int *imaxpow = iparm + iparm[I_MAXPOW];
3147 double *dvarpow = (double *)(dparm + iparm[I_DVPOW]);
3148
3149 const int *imaxp = imaxpow;
3150 double *dpowp = dvarpow;
3151 for (int jhat = 0; jhat < Nhat; jhat++, imaxp++) {
3152 double var = 1.0;
3153 for (int ip = 0; ip < *imaxp; ip++, dpowp++) {
3154 var *= rawcrd[jhat];
3155 *dpowp = var;
3156 }
3157 }
3158
3159 for (int k = 0; k < iparm[I_K]; k++, imaxp++) {
3160 double var = 1.0;
3161 for (int ip = 0; ip < *imaxp; ip++, dpowp++) {
3162 var *= aux[k];
3163 *dpowp = var;
3164 }
3165 }
3166
3167 // Loop for each term of the polynomial.
3168 *discrd = 0.0;
3169 const int *iflgp = iparm + iparm[I_FLAGS];
3170 const int *ipowp = iparm + iparm[I_IPOW];
3171 const double *dpolp = dparm + iparm[I_DPOLY];
3172 for (int m = 0; m < iparm[I_M]; m++) {
3173 double term = *(dpolp++);
3174
3175 // Loop over all variables.
3176 imaxp = imaxpow;
3177 dpowp = dvarpow - 1;
3178 for (int ivar = 0; ivar < iparm[I_NVAR]; ivar++) {
3179 if (*iflgp & 2) {
3180 // Nothing (zero power).
3181
3182 } else if (*iflgp) {
3183 // Integral power.
3184 if (*ipowp < 0) {
3185 // Negative.
3186 term /= dpowp[*ipowp];
3187 } else {
3188 // Positive.
3189 term *= dpowp[*ipowp];
3190 }
3191
3192 } else {
3193 // Fractional power.
3194 term *= pow(dpowp[0], *dpolp);
3195 }
3196
3197 iflgp++;
3198 ipowp++;
3199 dpolp++;
3200
3201 dpowp += *imaxp;
3202 imaxp++;
3203 }
3204
3205 *discrd += term;
3206 }
3207
3208 return 0;
3209 }
3210
3211 //----------------------------------------------------------------------------
3212
tpd1(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3213 int tpd1(
3214 int inverse,
3215 const int i[],
3216 const double p[],
3217 int Nhat,
3218 const double rawcrd[],
3219 double *discrd)
3220
3221 {
3222 if (i[I_TPDNCO+inverse] != 4 || 2 < Nhat) {
3223 return 1;
3224 }
3225
3226 double r, s;
3227 double u = rawcrd[0];
3228 double v = rawcrd[1];
3229
3230 // Auxiliary variables?
3231 if (i[I_TPDAUX]) {
3232 r = p[0] + p[1]*u + p[2]*v;
3233 v = p[3] + p[4]*u + p[5]*v;
3234 u = r;
3235 p += 6;
3236 }
3237
3238 if (inverse) p += i[I_TPDNCO];
3239
3240 // First degree.
3241 *discrd = p[0] + u*p[1];
3242
3243 if (Nhat == 1) return 0;
3244
3245 *discrd += v*p[2];
3246
3247 // Radial terms?
3248 if (i[I_TPDRAD]) {
3249 s = u*u + v*v;
3250 r = sqrt(s);
3251
3252 *discrd += r*p[3];
3253 }
3254
3255 return 0;
3256 }
3257
3258 //----------------------------------------------------------------------------
3259
tpd2(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3260 int tpd2(
3261 int inverse,
3262 const int i[],
3263 const double p[],
3264 int Nhat,
3265 const double rawcrd[],
3266 double *discrd)
3267
3268 {
3269 if (i[I_TPDNCO+inverse] != 7 || 2 < Nhat) {
3270 return 1;
3271 }
3272
3273 double r, s;
3274 double u = rawcrd[0];
3275 double v = rawcrd[1];
3276
3277 // Auxiliary variables?
3278 if (i[I_TPDAUX]) {
3279 r = p[0] + p[1]*u + p[2]*v;
3280 v = p[3] + p[4]*u + p[5]*v;
3281 u = r;
3282 p += 6;
3283 }
3284
3285 if (inverse) p += i[I_TPDNCO];
3286
3287 // Second degree.
3288 *discrd = p[0] + u*(p[1] + u*(p[4]));
3289
3290 if (Nhat == 1) return 0;
3291
3292 *discrd +=
3293 v*(p[2] + v*(p[6]))
3294 + u*(p[5])*v;
3295
3296 // Radial terms?
3297 if (i[I_TPDRAD]) {
3298 s = u*u + v*v;
3299 r = sqrt(s);
3300
3301 *discrd += r*p[3];
3302 }
3303
3304 return 0;
3305 }
3306
3307 //----------------------------------------------------------------------------
3308
tpd3(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3309 int tpd3(
3310 int inverse,
3311 const int i[],
3312 const double p[],
3313 int Nhat,
3314 const double rawcrd[],
3315 double *discrd)
3316
3317 {
3318 if (i[I_TPDNCO+inverse] != 12 || 2 < Nhat) {
3319 return 1;
3320 }
3321
3322 double r, s;
3323 double u = rawcrd[0];
3324 double v = rawcrd[1];
3325
3326 // Auxiliary variables?
3327 if (i[I_TPDAUX]) {
3328 r = p[0] + p[1]*u + p[2]*v;
3329 v = p[3] + p[4]*u + p[5]*v;
3330 u = r;
3331 p += 6;
3332 }
3333
3334 if (inverse) p += i[I_TPDNCO];
3335
3336 // Third degree.
3337 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7])));
3338
3339 if (Nhat == 1) return 0;
3340
3341 *discrd +=
3342 v*(p[2] + v*(p[6] + v*(p[10])))
3343 + u*(p[5] + v*(p[9])
3344 + u*(p[8]))*v;
3345
3346 // Radial terms?
3347 if (i[I_TPDRAD]) {
3348 s = u*u + v*v;
3349 r = sqrt(s);
3350
3351 *discrd += r*(p[3] + s*(p[11]));
3352 }
3353
3354 return 0;
3355 }
3356
3357 //----------------------------------------------------------------------------
3358
tpd4(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3359 int tpd4(
3360 int inverse,
3361 const int i[],
3362 const double p[],
3363 int Nhat,
3364 const double rawcrd[],
3365 double *discrd)
3366
3367 {
3368 if (i[I_TPDNCO+inverse] != 17 || 2 < Nhat) {
3369 return 1;
3370 }
3371
3372 double r, s;
3373 double u = rawcrd[0];
3374 double v = rawcrd[1];
3375
3376 // Auxiliary variables?
3377 if (i[I_TPDAUX]) {
3378 r = p[0] + p[1]*u + p[2]*v;
3379 v = p[3] + p[4]*u + p[5]*v;
3380 u = r;
3381 p += 6;
3382 }
3383
3384 if (inverse) p += i[I_TPDNCO];
3385
3386 // Fourth degree.
3387 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12]))));
3388
3389 if (Nhat == 1) return 0;
3390
3391 *discrd +=
3392 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16]))))
3393 + u*(p[5] + v*(p[9] + v*(p[15]))
3394 + u*(p[8] + v*(p[14])
3395 + u*(p[13])))*v;
3396
3397 // Radial terms?
3398 if (i[I_TPDRAD]) {
3399 s = u*u + v*v;
3400 r = sqrt(s);
3401
3402 *discrd += r*(p[3] + s*(p[11]));
3403 }
3404
3405 return 0;
3406 }
3407
3408 //----------------------------------------------------------------------------
3409
tpd5(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3410 int tpd5(
3411 int inverse,
3412 const int i[],
3413 const double p[],
3414 int Nhat,
3415 const double rawcrd[],
3416 double *discrd)
3417
3418 {
3419 if (i[I_TPDNCO+inverse] != 24 || 2 < Nhat) {
3420 return 1;
3421 }
3422
3423 double r, s;
3424 double u = rawcrd[0];
3425 double v = rawcrd[1];
3426
3427 // Auxiliary variables?
3428 if (i[I_TPDAUX]) {
3429 r = p[0] + p[1]*u + p[2]*v;
3430 v = p[3] + p[4]*u + p[5]*v;
3431 u = r;
3432 p += 6;
3433 }
3434
3435 if (inverse) p += i[I_TPDNCO];
3436
3437 // Fifth degree.
3438 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12] + u*(p[17])))));
3439
3440 if (Nhat == 1) return 0;
3441
3442 *discrd +=
3443 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16] + v*(p[22])))))
3444 + u*(p[5] + v*(p[9] + v*(p[15] + v*(p[21])))
3445 + u*(p[8] + v*(p[14] + v*(p[20]))
3446 + u*(p[13] + v*(p[19])
3447 + u*(p[18]))))*v;
3448
3449 // Radial terms?
3450 if (i[I_TPDRAD]) {
3451 s = u*u + v*v;
3452 r = sqrt(s);
3453
3454 *discrd += r*(p[3] + s*(p[11] + s*(p[23])));
3455 }
3456
3457 return 0;
3458 }
3459
3460 //----------------------------------------------------------------------------
3461
tpd6(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3462 int tpd6(
3463 int inverse,
3464 const int i[],
3465 const double p[],
3466 int Nhat,
3467 const double rawcrd[],
3468 double *discrd)
3469
3470 {
3471 if (i[I_TPDNCO+inverse] != 31 || 2 < Nhat) {
3472 return 1;
3473 }
3474
3475 double r, s;
3476 double u = rawcrd[0];
3477 double v = rawcrd[1];
3478
3479 // Auxiliary variables?
3480 if (i[I_TPDAUX]) {
3481 r = p[0] + p[1]*u + p[2]*v;
3482 v = p[3] + p[4]*u + p[5]*v;
3483 u = r;
3484 p += 6;
3485 }
3486
3487 if (inverse) p += i[I_TPDNCO];
3488
3489 // Sixth degree.
3490 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12] + u*(p[17] + u*(p[24]))))));
3491
3492 if (Nhat == 1) return 0;
3493
3494 *discrd +=
3495 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16] + v*(p[22] + v*(p[30]))))))
3496 + u*(p[5] + v*(p[9] + v*(p[15] + v*(p[21] + v*(p[29]))))
3497 + u*(p[8] + v*(p[14] + v*(p[20] + v*(p[28])))
3498 + u*(p[13] + v*(p[19] + v*(p[27]))
3499 + u*(p[18] + v*(p[26])
3500 + u*(p[25])))))*v;
3501
3502 // Radial terms?
3503 if (i[I_TPDRAD]) {
3504 s = u*u + v*v;
3505 r = sqrt(s);
3506
3507 *discrd += r*(p[3] + s*(p[11] + s*(p[23])));
3508 }
3509
3510 return 0;
3511 }
3512
3513 //----------------------------------------------------------------------------
3514
tpd7(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3515 int tpd7(
3516 int inverse,
3517 const int i[],
3518 const double p[],
3519 int Nhat,
3520 const double rawcrd[],
3521 double *discrd)
3522
3523 {
3524 if (i[I_TPDNCO+inverse] != 40 || 2 < Nhat) {
3525 return 1;
3526 }
3527
3528 double r, s;
3529 double u = rawcrd[0];
3530 double v = rawcrd[1];
3531
3532 // Auxiliary variables?
3533 if (i[I_TPDAUX]) {
3534 r = p[0] + p[1]*u + p[2]*v;
3535 v = p[3] + p[4]*u + p[5]*v;
3536 u = r;
3537 p += 6;
3538 }
3539
3540 if (inverse) p += i[I_TPDNCO];
3541
3542 // Seventh degree.
3543 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12] + u*(p[17] + u*(p[24] + u*(p[31])))))));
3544
3545 if (Nhat == 1) return 0;
3546
3547 *discrd +=
3548 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16] + v*(p[22] + v*(p[30] + v*(p[38])))))))
3549 + u*(p[5] + v*(p[9] + v*(p[15] + v*(p[21] + v*(p[29] + v*(p[37])))))
3550 + u*(p[8] + v*(p[14] + v*(p[20] + v*(p[28] + v*(p[36]))))
3551 + u*(p[13] + v*(p[19] + v*(p[27] + v*(p[35])))
3552 + u*(p[18] + v*(p[26] + v*(p[34]))
3553 + u*(p[25] + v*(p[33])
3554 + u*(p[32]))))))*v;
3555
3556 // Radial terms?
3557 if (i[I_TPDRAD]) {
3558 s = u*u + v*v;
3559 r = sqrt(s);
3560
3561 *discrd += r*(p[3] + s*(p[11] + s*(p[23] + s*(p[39]))));
3562 }
3563
3564 return 0;
3565 }
3566
3567 //----------------------------------------------------------------------------
3568
tpd8(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3569 int tpd8(
3570 int inverse,
3571 const int i[],
3572 const double p[],
3573 int Nhat,
3574 const double rawcrd[],
3575 double *discrd)
3576
3577 {
3578 if (i[I_TPDNCO+inverse] != 49 || 2 < Nhat) {
3579 return 1;
3580 }
3581
3582 double r, s;
3583 double u = rawcrd[0];
3584 double v = rawcrd[1];
3585
3586 // Auxiliary variables?
3587 if (i[I_TPDAUX]) {
3588 r = p[0] + p[1]*u + p[2]*v;
3589 v = p[3] + p[4]*u + p[5]*v;
3590 u = r;
3591 p += 6;
3592 }
3593
3594 if (inverse) p += i[I_TPDNCO];
3595
3596 // Eighth degree.
3597 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12] + u*(p[17] + u*(p[24] + u*(p[31] + u*(p[40]))))))));
3598
3599 if (Nhat == 1) return 0;
3600
3601 *discrd +=
3602 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16] + v*(p[22] + v*(p[30] + v*(p[38] + v*(p[48]))))))))
3603 + u*(p[5] + v*(p[9] + v*(p[15] + v*(p[21] + v*(p[29] + v*(p[37] + v*(p[47]))))))
3604 + u*(p[8] + v*(p[14] + v*(p[20] + v*(p[28] + v*(p[36] + v*(p[46])))))
3605 + u*(p[13] + v*(p[19] + v*(p[27] + v*(p[35] + v*(p[45]))))
3606 + u*(p[18] + v*(p[26] + v*(p[34] + v*(p[44])))
3607 + u*(p[25] + v*(p[33] + v*(p[43]))
3608 + u*(p[32] + v*(p[42])
3609 + u*(p[41])))))))*v;
3610
3611 // Radial terms?
3612 if (i[I_TPDRAD]) {
3613 s = u*u + v*v;
3614 r = sqrt(s);
3615
3616 *discrd += r*(p[3] + s*(p[11] + s*(p[23] + s*(p[39]))));
3617 }
3618
3619 return 0;
3620 }
3621
3622 //----------------------------------------------------------------------------
3623
tpd9(int inverse,const int i[],const double p[],int Nhat,const double rawcrd[],double * discrd)3624 int tpd9(
3625 int inverse,
3626 const int i[],
3627 const double p[],
3628 int Nhat,
3629 const double rawcrd[],
3630 double *discrd)
3631
3632 {
3633 if (i[I_TPDNCO+inverse] != 60 || 2 < Nhat) {
3634 return 1;
3635 }
3636
3637 double r, s;
3638 double u = rawcrd[0];
3639 double v = rawcrd[1];
3640
3641 // Auxiliary variables?
3642 if (i[I_TPDAUX]) {
3643 r = p[0] + p[1]*u + p[2]*v;
3644 v = p[3] + p[4]*u + p[5]*v;
3645 u = r;
3646 p += 6;
3647 }
3648
3649 if (inverse) p += i[I_TPDNCO];
3650
3651 // Ninth degree.
3652 *discrd = p[0] + u*(p[1] + u*(p[4] + u*(p[7] + u*(p[12] + u*(p[17] + u*(p[24] + u*(p[31] + u*(p[40] + u*(p[49])))))))));
3653
3654 if (Nhat == 1) return 0;
3655
3656 *discrd +=
3657 v*(p[2] + v*(p[6] + v*(p[10] + v*(p[16] + v*(p[22] + v*(p[30] + v*(p[38] + v*(p[48] + v*(p[58])))))))))
3658 + u*(p[5] + v*(p[9] + v*(p[15] + v*(p[21] + v*(p[29] + v*(p[37] + v*(p[47] + v*(p[57])))))))
3659 + u*(p[8] + v*(p[14] + v*(p[20] + v*(p[28] + v*(p[36] + v*(p[46] + v*(p[56]))))))
3660 + u*(p[13] + v*(p[19] + v*(p[27] + v*(p[35] + v*(p[45] + v*(p[55])))))
3661 + u*(p[18] + v*(p[26] + v*(p[34] + v*(p[44] + v*(p[54]))))
3662 + u*(p[25] + v*(p[33] + v*(p[43] + v*(p[53])))
3663 + u*(p[32] + v*(p[42] + v*(p[52]))
3664 + u*(p[41] + v*(p[51])
3665 + u*(p[50]))))))))*v;
3666
3667 // Radial terms?
3668 if (i[I_TPDRAD]) {
3669 s = u*u + v*v;
3670 r = sqrt(s);
3671
3672 *discrd += r*(p[3] + s*(p[11] + s*(p[23] + s*(p[39] + s*(p[59])))));
3673 }
3674
3675 return 0;
3676 }
3677