Updates for latest xrayvision

TCDSolar · May 15, 2024 · 4d7b12c · 4d7b12c
1 parent c606ef4
commit 4d7b12c
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 19 deletions.
diff --git a/examples/imaging_demo.py b/examples/imaging_demo.py
@@ -114,7 +114,7 @@
 # Set up image parameters
 
 imsize = [512, 512] * u.pixel  # number of pixels of the map to reconstruct
-pixel = [10, 10] * u.arcsec  # pixel size in aresec
+pixel = [10, 10] * u.arcsec/u.pixel  # pixel size in aresec
 
 ###############################################################################
 # Make a full disk back projection (inverse transform) map
@@ -149,18 +149,21 @@
 # Estimate the flare location and plot on top of back projection map.
 
 max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel
-# because WCS axes are reverse order
-max_hpc = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
+# because WCS axes and array are reversed
+max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0])
 
-ax0.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
-ax1.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
+ax0.plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
+ax1.plot_coord(max_stix, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
 
 
 ################################################################################
 # Use estimated flare location to create more accurate visibilities
 
+# because WCS axes and array are reversed and all positions are expected follow array indices
+flare_pos_stix = [max_stix.Ty.to_value('arcsec'), max_stix.Tx.to_value('arcsec')]*u.arcsec
+
 meta_pixels_sci = create_meta_pixels(
-    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=[max_hpc.Tx, max_hpc.Ty], no_shadowing=True
+    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=flare_pos_stix, no_shadowing=True
 )
 
 meta_pixels_bkg_subtracted = {"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
@@ -171,7 +174,7 @@
 uv_data = get_uv_points_data()
 vis.u = uv_data['u']
 vis.v = uv_data['v']
-cal_vis = calibrate_visibility(vis, [max_hpc.Tx, max_hpc.Ty])
+cal_vis = calibrate_visibility(vis, flare_pos_stix)
 
 ###############################################################################
 # Selected detectors 10 to 3
@@ -186,7 +189,7 @@
     vis=cal_vis.obsvis[idx],
     u=cal_vis.u[idx],
     v=cal_vis.v[idx],
-    offset=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec,
+    offset=flare_pos_stix,
 )
 skeys = stix_vis1.__dict__.keys()
 for k, v in cal_vis.__dict__.items():
@@ -197,7 +200,7 @@
 # Set up image parameters
 
 imsize = [129, 129] * u.pixel  # number of pixels of the map to reconstruct
-pixel = [2, 2] * u.arcsec  # pixel size in arcsec
+pixel = [2, 2] * u.arcsec / u.pixel  # pixel size in arcsec
 
 ###############################################################################
 # Create a back projection image with natural weighting
@@ -207,7 +210,7 @@
 ###############################################################################
 # Create a back projection image with uniform weighting
 
-bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, natural=False)
+bp_uni = vis_to_image(stix_vis1, imsize, pixel_size=pixel, scheme='uniform')
 
 ###############################################################################
 # Create a `sunpy.map.Map` with back projection
@@ -221,7 +224,7 @@
 gain = 0.1  # gain used in each clean iteration
 beam_width = 20.0 * u.arcsec
 clean_map, model_map, resid_map = vis_clean(
-    stix_vis1, imsize, pixel=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
+    stix_vis1, imsize, pixel_size=pixel, gain=gain, niter=niter, clean_beam_width=20 * u.arcsec
 )
 
 ###############################################################################
@@ -237,7 +240,7 @@
     vis=cal_vis.obsvis[idx],
     u=cal_vis.u[idx],
     v=cal_vis.v[idx],
-    center=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec,
+    offset=flare_pos_stix,
 )
 skeys = stix_vis1.__dict__.keys()
 for k, v in cal_vis.__dict__.items():
@@ -249,7 +252,7 @@
     stix_vis1,
     shape=imsize,
     pixel_size=pixel,
-    flare_xy=np.array([max_hpc.Tx.value, max_hpc.Ty.value]) * u.arcsec,
+    flare_xy=flare_pos_stix,
     idx=idx,
 )
 

diff --git a/stixpy/imaging/em.py b/stixpy/imaging/em.py
@@ -58,12 +58,9 @@ def get_transmission_matrix(u, v, shape=[64, 64] * apu.pix, pixel_size=[4.0, 4.0
     idx = [6, 28, 0, 24, 4, 22, 5, 29, 1, 14, 26, 30, 23, 7, 27, 20, 25, 3, 15, 13, 31, 2, 19, 21]
     phase_cor = phase_cor[idx]
 
-    m, n = shape.to_value("pix")
+    x = generate_xy(shape[0], pixel_size=pixel_size[0], phase_centre=center[0])
+    y = generate_xy(shape[1], pixel_size=pixel_size[1], phase_centre=center[1])
 
-    # center = center - [26.1, 58.2] * apu.arcsec
-
-    y = generate_xy(m, pixel_size=pixel_size[1], center=center[1])
-    x = generate_xy(n, pixel_size=pixel_size[0], center=center[0])
     x, y = np.meshgrid(x, y)
     # Check apu are correct for exp need to be dimensionless and then remove apu for speed
     if (u[0] * x[0, 0]).unit == apu.dimensionless_unscaled and (v[0] * y[0, 0]).unit == apu.dimensionless_unscaled:
@@ -154,5 +151,5 @@ def em(countrates, vis, shape, pixel_size, maxiter=5000, tolerance=0.001, *, fla
     m, n = shape.to_value("pix").astype(int)
     # Note the order keyword here to match `y.flatten(order='F')` above at line 108
     # TODO investigate why .T is needed probably because of fortran ordering
-    x_im = x.reshape(m, n, order="F").T
+    x_im = x.reshape(m, n, order="F")
     return x_im