FOXSI-R-2 PSF Analysis

In [1]:

    

#this code reads in all the data from the fits files and stores it in an hdf5 file
#only need to do this once
#%run /Users/schriste/Desktop/foxsi2_psf2hdf.py


    

In [2]:

    

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
from astropy.units import Quantity
import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter('ignore', append=True)

import h5py
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
import sunpy.map
from astropy.visualization.mpl_normalize import ImageNormalize
from astropy import visualization
import scipy
from scipy import interpolate
import scipy.optimize as opt
from matplotlib.colors import LogNorm
from astropy.modeling import models, fitting


    

In [3]:

    

warnings.warn('b')


    

In [4]:

    

#%run /Users/schriste/Desktop/FOXSIpsf/foxsi2_psfall.py


    

In [5]:

    

percent_levels = np.arange(1, 100, 10)

filename = '/Users/schriste/Desktop/FOXSIpsf/foxsi2_ccd_corrected.hdf5'
h = h5py.File(filename, 'r+')
images = h['X2/ccd_images']

polar_angles = u.Quantity(h['meta/polar_angle'][...], h['meta/polar_angle'].attrs['units'])
offaxis_angles = u.Quantity(h['meta/offaxis_angle'][...], h['meta/offaxis_angle'].attrs['units'])
nimages = images.shape[0]

def gauss2d((x,y), amplitude, xo, yo, sigma_x, sigma_y, theta):
    """A 2d dimensional Gaussian distribution with an arbitrary rotation, theta"""
    a = (np.cos(theta)**2)/(2*sigma_x**2) + (np.sin(theta)**2)/(2*sigma_y**2)
    b = -(np.sin(2*theta))/(4*sigma_x**2) + (np.sin(2*theta))/(4*sigma_y**2)
    c = (np.sin(theta)**2)/(2*sigma_x**2) + (np.cos(theta)**2)/(2*sigma_y**2)
    return amplitude*np.exp( - (a*((x-xo)**2) + 2*b*(x-xo)*(y-yo) + c*((y-yo)**2)))

def multi_gauss2d((x,y), amplitude, center, sigma_x, sigma_y, theta):
    i = 0
    for amp, sig_x, sig_y in zip(amplitude, sigma_x, sigma_y):
        g = gauss2d((x, y), amp, center[0], center[1], sig_x, sig_y, theta)
        if i == 0:
            result = g
        else:
            result += g
        i += 1
    return result

def three_gaussian((x,y),amp1,amp2,amp3, cen_x1, cen_y1, width_x1, width_x2, width_x3, width_y1, width_y2, width_y3, theta):
    "Returns a 3-gaussian function - all centered at the same x and y"""
    #try to force positive solutions
    if (amp1 < 0):
        amp1 = 1e12
    if (amp2 < 0):
        amp2 = 1e12
    if (amp3 < 0):
        amp3 = 1e12

    amplitude = (amp1, amp2, amp3)
    center = (cen_x1, cen_y1)
    sigma_x = (width_x1, width_x2, width_x3)
    sigma_y = (width_y1, width_y2, width_y3)

    result = multi_gauss2d((x,y), amplitude, center, sigma_x, sigma_y, theta)
    return result.ravel()



def psf(pitch, yaw, params, verbose=False):
    poly_params = []
    offaxis_angle = np.sqrt(pitch ** 2 + yaw ** 2)
    polar_angle = np.arctan2(pitch, yaw) + np.pi/2.
    for g in params:
        f = np.poly1d(g)
        poly_params.append(f(offaxis_angle))

    amplitude = (poly_params[0], poly_params[1], poly_params[2])
    width_x = [poly_params[3], poly_params[5], poly_params[7]]
    width_y = [poly_params[4], poly_params[6], poly_params[8]]
    if verbose:
        print(amplitude)
        print(width_x)
        print(width_y)
        print(polar_angle)
    return lambda (x, y): multi_gauss2d((x, y), amplitude, (0, 0), width_x, width_y, polar_angle)


def psf_map(pitch, yaw, params):
    p = psf(pitch, yaw, params)
    scale = 0.1
    x, y = np.meshgrid(np.arange(-60, 60, scale), np.arange(-60, 60, scale))
    im = p((x,y)).reshape(x.shape)
    iy, ix = np.unravel_index(np.argmax(im), im.shape)
    header =  {'cdelt1': scale,
               'cdelt2': scale,
               'telescop': 'FOXSI-2 Simulate',
               'crpix1': iy, 'crpix2': ix,
               'crval1': 0, 'crval2': 0}
    this_map = sunpy.map.Map(im, header).shift(pitch * u.arcmin, yaw * u.arcmin)
    return this_map



#def psf2((x, y), amplitude, center_x, center_y, width_x, width_y):
#    g1 = gaussian((x, y), amplitude, center_x, center_y, width_x, width_y)
#    if type(width_x) == 'list':
#        if len(width_x) >= 2:
#            g1 += gaussian((x, y), amplitude, center_x, center_y, width_x[1], width_y[1])
#        if len(width_x) >= 3:
#            g1 += gaussian((x, y), amplitude, center_x, center_y, width_x[2], width_y[2])
#    return g1.ravel()


def get_parabola_vertex(params):
    """Calculate the location of the vertex (the maximum or minimum) of a
        parabola given its parameters
           y = Ax^2 + Bx + C
        RETURNS: [x,y]
        WRITTEN: Steven Christe (25-Apr-2014)
    """
    xv = -params[1] / (2*params[0]);
    yv = params[2] - params[1] * params[1] / (4*params[0]);
    return xv, yv


def get_peak(x, y):
    #max_index = np.argmax(y)

    #xset = [x[max_index-1], x[max_index], x[max_index+1]]
    #yset = [y[max_index-1], y[max_index], y[max_index+1]]
    #params = fit_parabola_to_peak(xset, yset)
    #fit_max_xy = get_parabola_vertex(params)

    f = scipy.interpolate.interp1d(x, y, kind='cubic')
    max_x = opt.fmin(lambda x: -f(x), x[x.shape[0]/2.]+2, disp=False)
    fit_max_xy = [max_x[0], f(max_x[0])]
    return fit_max_xy


def get_fwhm(x, y):

    fit_max_xy = get_peak(x, y)
    half_max = fit_max_xy[1] / 2.0
    #this returns a function called f to be used in the future
    f = scipy.interpolate.interp1d(x, y, kind='cubic', bounds_error=False, fill_value=100)

    xlo = opt.brentq(lambda x: f(x) - half_max, fit_max_xy[0] - 5, fit_max_xy[0])
    xhi = opt.brentq(lambda x: f(x) - half_max, fit_max_xy[0], fit_max_xy[0] + 5)
    return xhi - xlo


def moments(data):
    """Returns (height, x, y, width_x, width_y)
    the gaussian parameters of a 2D distribution by calculating its
    moments """
    total = data.sum()
    X, Y = np.indices(data.shape)
    x = (X*data).sum()/total
    y = (Y*data).sum()/total
    col = data[:, int(y)]
    width_x = np.sqrt(abs((np.arange(col.size)-y)**2*col).sum()/col.sum())
    row = data[int(x), :]
    width_y = np.sqrt(abs((np.arange(row.size)-x)**2*row).sum()/row.sum())
    height = data.max()
    return height, x, y, width_x, width_y


class FoxsiCCDImage(object):

    def __init__(self, im, polar_angle, offaxis_angle, range=None, fit=True, shift=False):
        #self.background = np.average(im[0:100, 0:100])
        #im = im - self.background
        #im = im/im.max()
        # find the maximum
        self.scale = u.Quantity([1.3, 1.3], 'arcsec')
        self.map_fit = None
        iy, ix = np.unravel_index(np.argmax(im), im.shape)
        self.header = {'cdelt1': self.scale[0].value,
                       'cdelt2': self.scale[0].value,
                       'telescop': 'FOXSI-2',
                       'crpix1': iy, 'crpix2': ix,
                       'crval1': 0, 'crval2': 0}

        #yaw = np.sin(polar_angle*!PI/180.0)* offaxis_angle * 60
        #pitch = np.cos(polar_angle*!PI/180.0)* offaxis_angle * 60
        # the pitch must be reversed to correct for the yaw stage
        pitch, yaw = (-np.sin(polar_angle) * offaxis_angle, np.cos(polar_angle) * offaxis_angle)
        warnings
        self.map = self._create_map(im, range=range)
        if shift:
            self.map = self.map.shift(pitch, yaw)
        self.xaxis = self.map.pixel_to_data(np.arange(0, self.map.data.shape[0])*u.pix, 0 * u.pix)[0].value
        self.yaxis =  self.map.pixel_to_data(0*u.pix, np.arange(0, self.map.data.shape[1]) * u.pix)[1].value
        self.polar_angle = polar_angle
        self.offaxis_angle = offaxis_angle
        self.fit_parameters = {}
        if fit:
            self._fit()

    def _create_map(self, im, range=None):
        if range is None:
            range = [-35,35] * u.arcsec
        map = sunpy.map.Map(im, self.header).submap(range, range)
        map.plot_settings['norm'] = ImageNormalize(stretch=visualization.AsinhStretch(0.01))
        #m = moments(map.data)
        #map.meta['crpix1'] = m[1]
        #map.meta['crpix2'] = m[2]
        #map = sunpy.map.Map(map.data, map.meta)
        #print(m)
        return map

    def _fit(self):

        x, y = np.meshgrid(self.xaxis, self.yaxis)

        # develop some guesses for the fit
        moments_of_data = moments(self.map.data)
        #get the moment info into the right coordinates (data coords not pixel coords)
        guess_xcen = self.xaxis[0] + (self.map.scale[0].value * moments_of_data[1])
        guess_ycen = self.yaxis[0] + (self.map.scale[1].value * moments_of_data[2])
        guess_xwid = self.scale[0].value * moments_of_data[3]
        guess_ywid = self.scale[1].value * moments_of_data[4]
        m = [moments_of_data[0], guess_xcen, guess_ycen, guess_xwid, guess_ywid]
        angle = 0.0
        guess = [m[0],m[0]/2,m[0]/4,m[1],m[2],m[3],m[3]*1.4,m[3]*2,m[4],m[4]*1.4,m[4]*2, angle]

        self.fit_parameters.update({'guess': guess})

        # now fit the data
        popt, pcov = opt.curve_fit(three_gaussian, (x,y), self.map.data.ravel(), p0 = guess, full_output=False)

        self.fit_function = lambda x, y: three_gaussian((x,y), *popt)

        data_fitted = three_gaussian((x,y), *popt).reshape(self.map.data.shape)
        self.fit_parameters.update({'fit_function':lambda (x,y):three_gaussian((x,y), *popt)})
        self.fit_parameters.update({'fit_data':data_fitted})
        self.fit_parameters.update({'fit_residuals':(self.map.data - data_fitted) / self.map.data})
        #data_fitted = psf((x, y), *popt).reshape(self.map.data.shape)
        self.fit_parameters.update({'popt':popt})
        self.fit_parameters.update({'pcov':popt})

        # get the error on the parameters
        errorp = np.sqrt(np.diag(pcov))
        meta_fit = self.map.meta.copy()
        meta_fit.update({'telescop':'FOXSI-2 Fit'})
        self.map_fit = sunpy.map.Map(data_fitted, meta_fit)
        self.map_fit.plot_settings['norm'] = ImageNormalize(stretch=visualization.AsinhStretch(0.01))
        meta_resid = self.map.meta.copy()
        meta_resid.update({'telescop':'FOXSI-2 Residuals'})
        self.map_fit_residuals = sunpy.map.Map(self.fit_parameters['fit_residuals'], meta_resid)
        self.map_fit_residuals.plot_settings['norm'] = ImageNormalize(stretch=visualization.AsinhStretch(0.01))
        self.fit_parameters.update({'reduced_chi_squared':np.sum((self.map.data - self.map_fit.data) ** 2)/(self.map.data.size - len(popt))})
        self.fit_parameters.update({'fit_params':{'amp1':(popt[0], errorp[0]),
                           'amp2':(popt[1], errorp[1]),
                           'amp3':(popt[2], errorp[2]),
                           'center_x':(popt[3], errorp[3]),
                           'center_y':(popt[4], errorp[4]),
                           'width_x1':(popt[5], errorp[5]),
                           'width_x2':(popt[6], errorp[6]),
                           'width_x3':(popt[7], errorp[7]),
                           'width_y1':(popt[8], errorp[8]),
                           'width_y2':(popt[9], errorp[9]),
                           'width_y3':(popt[10],  errorp[10])
                           }})
        self.fit_parameters.update({'popt': popt})

    def plot(self, title=None):
        if self.map_fit:
            cmap = sunpy.map.Map([self.map, self.map, self.map_fit], composite=True)
            cmap.set_levels(1, percent_levels, percent=True)
            cmap.set_levels(2, percent_levels, percent=True)
        else:
            cmap = sunpy.map.Map([self.map, self.map], composite=True)
            cmap.set_levels(1, percent_levels, percent=True)
        cmap.plot()

        if title:
            plt.title(title)
        else:
            plt.title(str(self.polar_angle) + ' ' + str(self.offaxis_angle))

    def plot_fit(self, title=None):
        cmap = sunpy.map.Map([self.map, self.map_fit, self.map_fit], composite=True)
        cmap.set_alpha(0, 0)
        cmap.set_levels(2, percent_levels, percent=True)
        cmap.plot()

        if title:
            plt.title(title)
        else:
            plt.title(str(self.polar_angle) + ' ' + str(self.offaxis_angle))

    def cut_data(self, direction='x'):
        # create an interpolation function of the data
        x, y = np.meshgrid(self.xaxis, self.yaxis)
        fim = interpolate.interp2d(x, y, self.map.data)
        xmin, ymin = opt.fmin(lambda (x,y):-fim(x,y), (0,0))
        min = [xmin, ymin]
        maxval = self.fim(self.min[0], self.min[1])

        if direction == 'x':
            x = self.xaxis
            y = self.fim(x, self.min[1] ) / np.max(self.fim(x, self.min[1] ))
        else:
            x = self.yaxis
            y = self.fim(self.min[0], x ) / np.max(self.fim(self.min[0], x ))

        return x, y

    def cut_plot(self, direction='x'):
        if direction == 'x':
            x = self.xaxis
            y = self.fim(x, self.min[1] ) / np.max(self.fim(x, self.min[1] ))
            #x = self.map.pixel_to_data(np.arange(0, self.map.data.shape[0])*u.pix, 0*u.pix)[0].value
            #y = self.map.data[np.int(self.fit_params['center_y']),:]
            #yfit = psf((x, self.min[1]), *popt)
            yfit = three_gaussian((x,self.min[1]),*popt)
            yfit = yfit / np.max(yfit)
            label = 'X [arcsec]'
        else:
            x = self.yaxis
            y = (self.fim(self.min[0], x ) / np.max(self.fim(self.min[0], x ))).ravel()
            #yfit = psf((self.min[0], y), *popt)
            yfit = three_gaussian((self.min[0],x),*popt)
            yfit = yfit / np.max(yfit)
            label = 'Y [arcsec]'
        peak_xy = get_peak(x, y)
        fwhm = get_fwhm(x, y)

        plt.plot(x, y, label='data')
        plt.plot(x, yfit, 'b', label='fit',linewidth=2,linestyle='dashed')
        plt.xlabel(label)
        plt.axvline(peak_xy[0], color='g', label='x=%0.2f' % peak_xy[0])
        plt.axhline(peak_xy[1], color='g', label='y=%0.2f' % peak_xy[1])
        f = scipy.interpolate.interp1d(x, y, kind='cubic')
        xinter = np.arange(-30, 30, 0.1)
        plt.plot(xinter, f(xinter), label='inter')

        plt.axvspan(peak_xy[0] - fwhm/2., peak_xy[0] + fwhm/2., alpha=0.3, color='g', label='fwhm=%0.2f' % fwhm)
        plt.ylim(0, 1.1)
        plt.xlim(-30,30)
        leg = plt.legend(prop={'size':10})

        for legobj in leg.legendHandles:
            legobj.set_linewidth(2.0)

    def plot_all(self):
        fig = plt.subplots(figsize=(12, 4))

        plt.subplot(1,4,1)
        self.plot()

        plt.subplot(1,4,2)
        self.map_fit.plot()

        plt.subplot(1,4,3)
        self.map.plot()

        plt.subplot(1,4,4)
        self.map_fit_residuals.plot()

        plt.show()


    

In [6]:

    

c = FoxsiCCDImage(images[0], 90 * u.deg, 0.3 * u.arcmin, shift=True)


    

In [7]:

    

c.plot_all()


    

In [8]:

    

polar_angles


    

    
Out[8]:





$[0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~0,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~45,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~90,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135,~135] \; \mathrm{{}^{\circ}}$

In [9]:

    

offaxis_angles


    

    
Out[9]:





$[0,~0.5,~1,~1.5,~2,~3,~5,~7,~9,~-0.5,~-1,~-1.5,~-2,~-3,~-5,~-7,~-9,~0,~0.5,~1,~1.5,~2,~3,~5,~7,~9,~-0.5,~-1,~-1.5,~-2,~-3,~-5,~-7,~-9,~0,~0.5,~1,~1.5,~2,~3,~5,~7,~9,~-0.5,~-1,~-1.5,~-2,~-3,~-5,~-7,~-9,~0,~0.5,~1,~1.5,~2,~3,~5,~7,~9,~-0.5,~-1,~-1.5,~-2,~-3,~-5,~-7,~-9] \; \mathrm{{}^{\prime}}$

In [10]:

    

fit_params = []
ccds = []
for im, polar_angle, offaxis_angle in zip(images, polar_angles, offaxis_angles):
    c = FoxsiCCDImage(im, polar_angle, offaxis_angle, fit=False, shift=True)
    ccds.append(c)


    

In [11]:

    

fig = plt.figure(figsize=(8, 8))
for ccd in ccds:
    ccd.plot()

circle_rad = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9]) * 60
for rad in circle_rad:
    circle=plt.Circle((0,0),rad, fill=False, color='grey', linestyle='dashed')
    fig.gca().add_artist(circle)
plt.title("Data")
#fov = plt.Rectangle((-3.2 * 160 / 2., -3.2 * 160 / 2.), 3.2 * 160, 3.2 * 160., fill=False)
#fig.gca().add_artist(fov)
#plt.savefig('foxsi2_allpsf_data.pdf')


    

    
Out[11]:





<matplotlib.text.Text at 0x11c4de890>