Skip to content

Commit

Permalink
changes after review
Browse files Browse the repository at this point in the history
  • Loading branch information
@jemorrison
jemorrison committed Jul 20, 2021
1 parent 62125a3 commit 5e8958e
Show file tree
Hide file tree
Showing 5 changed files with 122 additions and 277 deletions.
19 changes: 5 additions & 14 deletions jwst/cube_build/cube_build_wcs_util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,20 +113,11 @@ def find_corners_MIRI(input, this_channel, instrument_info, coord_system):

if coord_system != 'internal_cal':
# before returning, ra should be between 0 to 360
if a_min < 0:
a_min = a_min + 360
if a_max >= 360.0:
a_max = a_max - 360.0

if a1 < 0:
a1 = a1 + 360
if a1 > 360.0:
a1 = a1 - 360.0

if a2 < 0:
a2 = a2 + 360
if a2 > 360.0:
a2 = a2 - 360.0
a_min %= 360
a_max %= 360

a1 %= 360
a2 %= 360

return a_min, b1, a_max, b2, a1, b_min, a2, b_max, lambda_min, lambda_max
# *****************************************************************************
Expand Down
45 changes: 20 additions & 25 deletions jwst/cube_build/ifu_cube.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -546,11 +546,6 @@ def build_ifucube(self):
self.spaxel_iflux = np.zeros(total_num, dtype=np.float64)
self.spaxel_dq = np.zeros(total_num, dtype=np.uint32)

for i in range(total_num):
self.spaxel_flux[i] = 0.0
self.spaxel_iflux[i] = 0.0
self.spaxel_weight[i] = 0.0
self.spaxel_var[i] = 0.0
# ______________________________________________________________________________
subtract_background = True

Expand Down Expand Up @@ -578,26 +573,23 @@ def build_ifucube(self):
# POINTCLOUD used for skyalign and IFUalign
# --------------------------------------------------------------------------------
if self.interpolation == 'pointcloud':
t0 = time.time()
pixelresult = self.map_detector_to_outputframe(this_par1,
subtract_background,
input_model)

coord1, coord2, wave, flux, err, slice_no, rois_pixel, roiw_pixel, weight_pixel,\
softrad_pixel, scalerad_pixel = pixelresult

# check that there is valid data returned
# If all the data is flagged as DO_NOT_USE - not common then log warning and skip data
# by default flag the dq plane based on the FOV of the detector projected to sky
flag_dq_plane = 1
if wave.size == 0:
log.warning(f'No valid data found on file {input_model.meta.filename}')
else:
flag_dq_plane = 0
if self.skip_dqflagging:
flag_dq_plane = 0

t1 = time.time()
log.debug("Time to transform pixels to output frame = %.1f s" % (t1 - t0,))
# check that there is valid data returned
# If all the data is flagged as DO_NOT_USE - not common- then log warning
if wave.size == 0:
log.warning(f'No valid data found on file {input_model.meta.filename}')
flag_dq_plane = 0

t0 = time.time()
roiw_ave = np.mean(roiw_pixel)
Expand Down Expand Up @@ -1391,10 +1383,10 @@ def map_detector_to_outputframe(self, this_par1,

if self.linear_wavelength:
min_wave_tolerance = self.zcoord[0] - self.roiw
max_wave_tolerance = (self.zcoord[-1] + self.roiw)
max_wave_tolerance = self.zcoord[-1] + self.roiw
else:
min_wave_tolerance = self.zcoord[0] - np.max(self.roiw_table)
max_wave_tolerance = (self.zcoord[-1] + np.max(self.roiw_table))
max_wave_tolerance = self.zcoord[-1] + np.max(self.roiw_table)

valid_min = np.where(wave_all >= min_wave_tolerance)
not_mapped_low = wave_all.size - len(valid_min[0])
Expand All @@ -1414,8 +1406,8 @@ def map_detector_to_outputframe(self, this_par1,
# all_flags = (dqflags.pixel['DO_NOT_USE'] +
# dqflags.pixel['NON_SCIENCE'])

valid3 = np.logical_and((wave_all >= min_wave_tolerance),
(wave_all <= max_wave_tolerance))
valid3 = np.logical_and(wave_all >= min_wave_tolerance,
wave_all <= max_wave_tolerance)

# find the location of good data
# good_data = np.where((np.bitwise_and(dq_all, all_flags) == 0) &
Expand All @@ -1427,13 +1419,16 @@ def map_detector_to_outputframe(self, this_par1,

# good data holds the location of pixels we want to map to cube
# define variables as numpy arrays (numba needs this defined)
flux = np.zeros(flux_all[good_data].shape, dtype=np.float64)
err = np.zeros(flux_all[good_data].shape, dtype=np.float64)
coord1 = np.zeros(flux_all[good_data].shape, dtype=np.float64)
coord2 = np.zeros(flux_all[good_data].shape, dtype=np.float64)
wave = np.zeros(flux_all[good_data].shape, dtype=np.float64)
slice_no = np.zeros(flux_all[good_data].shape)
flux[:] = flux_all[good_data]
flux_all_good = flux_all[good_data]
good_shape = flux_all_good.shape
flux = np.zeros(good_shape, dtype=np.float64)
err = np.zeros(good_shape, dtype=np.float64)
coord1 = np.zeros(good_shape, dtype=np.float64)
coord2 = np.zeros(good_shape, dtype=np.float64)
wave = np.zeros(good_shape, dtype=np.float64)
slice_no = np.zeros(good_shape)

flux[:] = flux_all_good
err[:] = err_all[good_data]
wave[:] = wave_all[good_data]
slice_no[:] = slice_no_all[good_data]
Expand Down
74 changes: 12 additions & 62 deletions jwst/cube_build/src/cube_match_internal.c
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,13 @@
/*
The detector pixels are represented by a 'point could' on the sky. The IFU cube is
represented by a 3-D regular grid. This module finds the point cloud members contained
in a region centered on the center of the cube spaxel. The size of the spaxel is spatial
coordinates is cdetl1 and cdelt2, while the wavelength size is zcdelt3.
This module uses the e modified shephard weighting method to determine how to weight each point clold member
in the spaxel.
Thie module forms the IFU cube in the "detector plane". This method of creating cubes is
for an engineering mode. The detector pixels are mapped to the along slice coordinate
and wavelength. The slice number represents the across slice coordinate. The pixels from
each slice are assumed to not overlap on the output plane with pixels from other slices. This
breaks the 3-D representation of the pixel's coordinates in the output plane to 2D. Standard
drizzling type routines are used to calculate the percentage of each pixel's overlap with
the output spaxels. The percentage of the area of the detector pixel overlapped onto an output
spaxel is used to weight the pixel flux for determining the spaxel flux.
Main function for Python: cube_wrapper_internal
Expand Down Expand Up @@ -91,52 +94,8 @@ extern double sh_find_overlap(const double xcenter, const double ycenter,
extern double find_area_quad(double MinX, double MinY, double Xcorner[], double Ycorner[]);


//_______________________________________________________________________
// Allocate the memory to for the spaxel arrays in the cube
//_______________________________________________________________________

int mem_alloc(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv) {

double *f, *w, *v, *i;
const char *msg = "Couldn't allocate memory for output arrays.";

// flux:
f = (double*)malloc(nelem* sizeof(double));
if (f) {
*fluxv = f;
} else {
PyErr_SetString(PyExc_MemoryError, msg);
return 1;
}

// weight:
w = (double*)malloc(nelem* sizeof(double));
if (w) {
*weightv = w;
} else {
PyErr_SetString(PyExc_MemoryError, msg);
return 1;
}

// variance:
v = (double*)malloc(nelem* sizeof(double));
if (v) {
*varv = v;
} else {
PyErr_SetString(PyExc_MemoryError, msg);
return 1;
}
// iflux
i = (double*)malloc(nelem* sizeof(double));
if (i) {
*ifluxv = i;
} else {
PyErr_SetString(PyExc_MemoryError, msg);
return 1;
}
extern alloc_flux_arrays(int nelem, double **fluxv, double **weightv, double **varv, double **ifluxv);

return 0;
}

//_______________________________________________________________________
// based on a 2-D drizzling type of algorithm find the contribution of
Expand All @@ -151,20 +110,11 @@ int match_detector_cube(int instrument, int naxis1, int naxis2, int nz, int npt,
double *pixel_flux, double *pixel_err,
double **spaxel_flux, double **spaxel_weight, double **spaxel_var,double **spaxel_iflux){

double *fluxv = NULL, *weightv=NULL, *varv=NULL ; // vectors for spaxel
double *ifluxv = NULL; // vector for spaxel
double *fluxv, *weightv, *varv, *ifluxv ; // vectors for spaxel

// allocate memory to hold output
if (mem_alloc(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1;
if (alloc_flux_arrays(ncube, &fluxv, &weightv, &varv, &ifluxv)) return 1;

double set_zero=0.0;
// Set all data to zero
for (int i = 0; i < ncube; i++){
varv[i] = set_zero;
fluxv[i] = set_zero;
ifluxv[i] = set_zero;
weightv[i] = set_zero;
}

// loop over each valid point on detector and find match to IFU cube based
// on along slice coordinate and wavelength
Expand Down
Loading

0 comments on commit 5e8958e

Please sign in to comment.