'''
Created on 29 nov. 2017
@author: user
'''
import logging
import numpy as np
import copy
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.collections import PatchCollection
from ecpi.common.instru.array_convention import yzegp_to_ijdet_array
from ecpi.common.num.array_square_cell import ArraySquareCell
from ecpi.common import add_path_data_ref_eclairs
from ecpi.common.io import fits_tools
[docs]class InstruXbase(object):
'''
mask convention : 1 for hole, 0 stop photon
'''
def __init__(self, p_mask : ArraySquareCell, p_detec : ArraySquareCell, p_dist_mask_detec):
'''
'''
self.mask = ArraySquareCell()
self.mask = p_mask
self.mask.set_origin_center()
self.detect = ArraySquareCell()
self.detect = p_detec
assert isinstance(self.mask, ArraySquareCell)
assert isinstance(self.detect, ArraySquareCell)
self.detect.set_origin_center()
self.dist_mask_detect = p_dist_mask_detec # [cm]
# declare logger
self.logger = logging.getLogger(__name__)
self.logger.info('creating an instance of InstruXbase')
[docs] def set_pos_source_fov(self, elevation, direction):
"""
set definition of the translation vector linked to the source
vec_trans_source is also position on mask where photon passed to go on origin
"""
self.vec_trans_source = np.array([0,0], dtype=np.float32 )
dd = self.dist_mask_detect/np.tan(np.array(elevation,dtype=np.float32))
self.vec_trans_source[0] = np.cos(direction)
self.vec_trans_source[1] = np.sin(direction)
self.vec_trans_source *= dd
[docs] def get_surface_rect_proj(self):
"""
compute intersection detector and shadow mask then compute its surface
:return: intersection surface in cm2
:rtype: float
"""
self.rect_proj = self.detect.intersection(self.mask, -self.vec_trans_source)
# why -self.vec_trans_source ?
# we want intersection in detector referentiel so we translate position mask to detector
inter_rec = self.rect_proj
surface = (inter_rec[1,0] - inter_rec[0,0])*(inter_rec[1,1] - inter_rec[0,1])
return surface
[docs] def get_full_surface_mask(self):
"""in cm2"""
return self.mask.get_surface()
[docs] def get_open_surface_mask(self):
"""in cm2"""
return self.mask.ar_val.sum()*self.mask.surf_cell
[docs] def get_percent_open_mask(self):
return self.get_open_surface_mask()/self.get_full_surface_mask()
[docs] def get_total_photon(self):
return self.detect.ar_val.sum()
[docs] def get_elevation_limit_fov_side(self):
"""for special case where mask and detector are square and with the same center
return deg
"""
half_size_mask = (self.get_full_surface_mask()**.5)/2
half_size_detect = (self.detect.get_surface()**.5)/2
angle = np.arctan(self.dist_mask_detect/(half_size_mask + half_size_detect))
return np.rad2deg(angle)
[docs] def get_elevation_limit_fov_corner(self):
"""for special case where mask and detector are square and with the same center
return deg
"""
half_size_mask = (self.get_full_surface_mask()**.5)/(2**.5)
half_size_detect = (self.detect.get_surface()**.5)/(2**.5)
angle = np.arctan(self.dist_mask_detect/(half_size_mask + half_size_detect))
return np.rad2deg(angle)
[docs] def skypix_to_elevdir(self, pix_y, pix_z):
"""convert pixel coordinates into elevation and direction
(0,0) is at the center of the sky array
elevation from 0 to 90
direction from -pi to +pi
:param pix_y: y coordinate of the source, in pixels.
:type pix_y: float
:param pix_z: z coordinates of the source, in pixels.
:type pix_z: float
:return: elevation and direction of the source
:rtype: float, float
"""
direction = np.degrees(np.arctan2(pix_z, pix_y))
norm_vec_trans_source = np.sqrt(pix_y**2 + pix_z**2) * self.detect.size_cm
elev = np.degrees(np.arctan2(self.dist_mask_detect, norm_vec_trans_source))
return elev, direction
[docs] def elevdir_to_skypix(self, elev, direction):
"""convert elevation/direction coordinates into sky pixels coordinates
(0,0) is at the center of the sky array
:param elev: elevation in degrees
:type elev: float or array(float)
:param direction: direction in degrees
:type direction: float or array(float)
:return: pixels y,z coordinates
:rtype: float, float or array2D(float)
"""
dd = self.dist_mask_detect/np.tan(np.radians(elev)) / self.detect.size_cm
pix_y = dd*np.cos(np.radians(direction))
pix_z = dd*np.sin(np.radians(direction))
return pix_y, pix_z
[docs]class InstruXMultiMask(InstruXbase):
def __init__(self, p_detec : ArraySquareCell, p_dist_mask_detec):
self.detect = p_detec
self.detect.set_origin_center()
self.dist_mask_detect = p_dist_mask_detec
self.l_mask = []
self.mask = None
[docs] def add_mask(self, p_mask: ArraySquareCell):
p_mask.perct_open = p_mask.ar_val.sum()*1.0/p_mask.ar_val.size
p_mask.name = str(len(self.l_mask)+1)
self.l_mask.append(p_mask)
[docs] def set_mask(self, p_mask):
self.mask = p_mask
[docs] def select_mask(self, p_idx):
self.mask = self.l_mask[p_idx]
[docs] def plot_submasks(self): # pragma: no cover
for mask in self.l_mask:
mask.plot()
plt.figure()
self.logger.debug(f'lenght of l_mask: {len(self.l_mask)}')
plt.title('Pos mask')
for mask in self.l_mask:
plt.plot(mask.to_origin_cm[0], mask.to_origin_cm[1], 'o')
[docs] def plot_mask(self): # pragma: no cover
_, ax = plt.subplots(1)
ax.patch.set_facecolor('red')
amax = 27.01
ax.set_xlim([-amax, amax])
ax.set_ylim([-amax, amax])
l_holes = []
# rect = Rectangle((0, 0), 1, 1)
# l_holes.append(rect)
# rect = Rectangle((3, 3), 1, 1)
# l_holes.append(rect)
for mask in self.l_mask:
assert isinstance(mask, ArraySquareCell)
for idx_x in range(mask.get_nb_cell_x()):
for idx_y in range(mask.get_nb_cell_y()):
idx = np.array([[idx_x, idx_y]])
if mask.value(idx[0]) == 1:
pos = mask.idx2pos_inf(idx)[0]
#print(pos[0], pos[1])
rect = Rectangle((pos[0], pos[1]), mask.size_cm, mask.size_cm)
l_holes.append(rect)
pc = PatchCollection(l_holes, facecolor='white', alpha=1)
# Add collection to axes
ax.add_collection(pc)
plt.show()
[docs] def get_full_surface_mask(self):
raise
[docs] def get_open_surface_mask(self):
"""in cm2"""
s_surf = 0
for r_mask in self.l_mask:
s_surf += r_mask.ar_val.sum()*r_mask.surf_cell
return s_surf
[docs] def get_surface_rect_proj_allmasks(self):
"""compute the surface on the detector of the projection of the masks (the coded detector surface)
It is source dependent
:return: coded detector surface in cm2
:rtype: float
"""
inter_rec = np.array([[99999, 99999], [-99999, -99999]])
for mask in self.l_mask:
self.set_mask(mask)
self.rect_proj = self.detect.intersection(self.mask, -self.vec_trans_source)
inter_rec[0,0] = min(inter_rec[0,0], self.rect_proj[0,0])
inter_rec[0,1] = min(inter_rec[0,1], self.rect_proj[0,1])
inter_rec[1,0] = max(inter_rec[1,0], self.rect_proj[1,0])
inter_rec[1,1] = max(inter_rec[1,1], self.rect_proj[1,1])
coded_surface = (inter_rec[1,0] - inter_rec[0,0])*(inter_rec[1,1] - inter_rec[0,1])
self.select_mask(0)
return coded_surface
[docs]class InstruECLAIRs(InstruXMultiMask):
"""Implementation of the ECLAIRs mask geometry.
"""
def __init__(self, show=None, p_file_mask=None, p_nb_pix=None):
"""**Constructor**
"""
if not p_file_mask:
p_file_mask = add_path_data_ref_eclairs('instru/maskECLAIRs.fits')
# referential detector is center of detector
mask, _, header = fits_tools.set_array_and_header_from_file(p_file_mask,
ugts_standard=1)
# get the characteristics of the instrument from keywords
if not p_nb_pix:
self.p_nb_pix = header['DETPIX']
else:
self.p_nb_pix = p_nb_pix
self.half_cross = header['MASNERV'] / 2
self.det_pix_size = header['DETPSIZE']
self.pitch_pix_size = header['PIXPITCH']
self.dist_mask_detect = header['DETMAS']
self.mask_pix_size = header['MASPSIZE']
self.mask_nb_pix = header['MASPIX']
self.mask_total_size = header['MASSIZE']
#hdul.close()
detect = ArraySquareCell(self.p_nb_pix, self.p_nb_pix, self.pitch_pix_size)
super().__init__(detect, self.dist_mask_detect)
# Now origin is the center of detector
if show == True:
self.logger.debug(f'mask shape: {mask.shape}')
plt.figure()
plt.title("fits mask")
plt.imshow(yzegp_to_ijdet_array(mask))
plt.show()
pm23_hc = self.mask_pix_size * self.mask_nb_pix / 2 + self.half_cross
# left bottom
m1 = ArraySquareCell(self.mask_nb_pix//2, self.mask_nb_pix//2, self.mask_pix_size)
m1.set_array(mask[:self.mask_nb_pix//2, :self.mask_nb_pix//2])
pos_corner_lb = np.array([-pm23_hc, -pm23_hc], dtype=np.float32)
m1.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m1)
# right bottom
m2 = ArraySquareCell(self.mask_nb_pix//2, self.mask_nb_pix//2, self.mask_pix_size)
m2.set_array(mask[self.mask_nb_pix//2:, :self.mask_nb_pix//2])
pos_corner_lb = np.array([self.half_cross, -pm23_hc], dtype=np.float32)
m2.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m2)
# left up
m3 = ArraySquareCell(self.mask_nb_pix//2, self.mask_nb_pix//2, self.mask_pix_size)
m3.set_array(mask[:self.mask_nb_pix//2, self.mask_nb_pix//2:])
pos_corner_lb = np.array([-pm23_hc, self.half_cross], dtype=np.float32)
m3.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m3)
# right up
m4 = ArraySquareCell(self.mask_nb_pix//2, self.mask_nb_pix//2, self.mask_pix_size)
m4.set_array(mask[self.mask_nb_pix//2:, self.mask_nb_pix//2:])
pos_corner_lb = np.array([self.half_cross, self.half_cross], dtype=np.float32)
m4.set_pos_corner_lb(pos_corner_lb)
self.add_mask(m4)
self.select_mask(0)
[docs] def get_full_surface_mask(self):
s_surf = 4*self.mask.get_surface()
# add cross without center
s_surf += (2*self.half_cross)*(4*self.mask.get_sizecm()[0])
# add center cross
s_surf += (2*self.half_cross)**2
return s_surf