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