From c8fcfc40cefe0376a6b6c2fdaf7ac45aabfca463 Mon Sep 17 00:00:00 2001
From: Mihai Cara <mihail.cara@gmail.com>
Date: Fri, 29 Aug 2025 00:22:42 -0400
Subject: [PATCH] Fix PC to CD matrix computation in HSTWCS

---
 CHANGES.rst                |  5 +++++
 stwcs/tests/test_hstwcs.py | 26 ++++++++++++++++++++++++++
 stwcs/wcsutil/hstwcs.py    |  2 +-
 3 files changed, 32 insertions(+), 1 deletion(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 0d64924..4d6f479 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -8,10 +8,15 @@ Unreleased
 
 - Return default response when request to the GSC service failed. [#227]
 
+- Fix a bug in ``stwcs.wcsutil.hstwcs.HSTWCS.pc2cd()`` that would result in
+  incorrect CD matrix when the WCS has a PC matrix. [#228]
+
+
 1.7.5 (2025-07-14)
 ------------------
 - Work around an issue with the AStrometry DB which returns NaNs in certain cases. [#221]
 
+
 1.7.4 (2025-01-07)
 ------------------
 
diff --git a/stwcs/tests/test_hstwcs.py b/stwcs/tests/test_hstwcs.py
index 2f737dd..2c1eff9 100644
--- a/stwcs/tests/test_hstwcs.py
+++ b/stwcs/tests/test_hstwcs.py
@@ -1,5 +1,6 @@
 from astropy.io import fits
 from stwcs.wcsutil import hstwcs, HSTWCS
+import numpy as np
 
 
 def test_radesys():
@@ -22,3 +23,28 @@ def test_prihdu_with_extver_no_extname():
     extname = HSTWCS(hdulist).extname
     assert extname == ('PRIMARY', 7)
     assert hdulist[extname] is hdulist[0]
+
+
+def test_pc2cd():
+    e = np.exp(1)
+    hdulist = fits.HDUList([
+        fits.PrimaryHDU(header=fits.Header([('extver', 1)])),
+        fits.ImageHDU(header=fits.Header([
+            ('pc1_1', 0.1),
+            ('pc1_2', 0.08),
+            ('pc2_1', -0.05),
+            ('pc2_2', 0.7),
+            ('cdelt1', e),
+            ('cdelt2', np.pi),
+            ('crpix1', 512.0),
+            ('crpix2', 512.0),
+            ('crval1', 150.0),
+            ('crval2', 2.0),
+            ('ctype1', 'RA---TAN'),
+            ('ctype2', 'DEC--TAN'),
+        ]))
+    ])
+    wcs = HSTWCS(hdulist, ext=1)
+    assert wcs.wcs.has_cd()
+    cd = np.array([[0.1 * e, 0.08 * e], [-0.05 * np.pi, 0.7 * np.pi]])
+    assert np.allclose(wcs.wcs.cd, cd)
diff --git a/stwcs/wcsutil/hstwcs.py b/stwcs/wcsutil/hstwcs.py
index d19f214..ac19a74 100644
--- a/stwcs/wcsutil/hstwcs.py
+++ b/stwcs/wcsutil/hstwcs.py
@@ -473,7 +473,7 @@ def _idc2hdr(self):
 
     def pc2cd(self):
         if self.wcs.has_pc():
-            self.wcs.cd = self.wcs.pc * self.wcs.cdelt[1]
+            self.wcs.cd = self.pixel_scale_matrix
 
     @deprecated_renamed_argument('accuracy', 'tolerance', '1.6.1', arg_in_kwargs=True)
     def all_world2pix(self, *args, **kwargs):