Unfortunately SWARP doesn't provide the version number of WCSLIB that it is using. (just as a side-note, in Gnuastro, version information is taken very seriously and written to all output headers as described in the book.)

So I guess the only way would be to either dig into its code or ask Emmanuel. But from the formatting of the WCS related keywords, you can guess that SWARP doesn't use WCSLIB's FITS keyword output. In Gnuastro we directly write WCSLIB's produced headers into the FITS file. You can see that both the ordering and the comments on the keywords differs between SWARP and WCSLIB (Gnuastro's outputs). I think SWARP writes the keywords itsself directly from the structure.

Maybe we can make a feature request to WSCLIB. Currently `wcshdo' does take some flags for telling it how to print the header information. So we might be able to convince Mark Calabretta to add another flag for writing the `CDi_j' elements.

I just changed the category of this issue to "All Gnuastro" because all Gnuastro's programs ultimately use a single function to read and write FITS files, thus this discussion applies to all of them that deal with FITS WCS information (the bug #50073 and bug #50073 also applied to all Gnuastro's programs).

Thanks for this very intriguing suggestion. Indeed the `CDi_j' formalism is the most common in most imaging archives.

This is what the WCSLIB manual says under `wcsprm::altin':

"These alternate specifications of the linear transformation matrix are translated immediately to PCi_ja by wcsset() and are invisible to the lower-level WCSLIB routines. In particular, wcsset() resets wcsprm::cdelt to unity if CDi_ja is present (and no PCi_ja )."

I couldn't find any better explanation. I would guess that since WCSLIB is intended for any kind of data (not just 2D images), for example spectra, 3D data cubes, wide full-sky surveys, and etc, it has chosen this formalism since the scale is separate and more clearly defined. I should just add that with the fix to bug #50073, the `CDELTi's in the `PCi_j' formalism are corrected.

Ultimately, if its not possible to tell WCSLIB to print the CD matrix, it is not hard to manually change the `PCi_j' and `CDELTi' key words to `CDi_j' immediately before saving then in the FITS. I am just slightly concerned because WCSLIB has already been written to be the most general library for the job, taking the most general applications in mind. But here, we are only focused on 2D imaging surveys. This is also the type of data I deal with most, but we don't want to limit Gnuastro to this type of data.

Just for reference the conversion from `wcsprm' to text is done by the call to `wcshdo' in `lib/fits.c' and `gal_fits_add_wcs_to_header' is in charge of writing that text into the output FITS header. Some corrections can be done in either of these two steps (before writing into the FITS or after it).

So, let's keep this discussion going on here (not in the task tracker, so everyone can easily post comments). Please share your thoughts on what the best way forward can be.