<jupyter_start><jupyter_text>Open data file<jupyter_code>ROOTfile = "./Crab_LZA_excess_uncorr.fits"
HDUlist = fits.open(ROOTfile)
HDUlist.info()
gal_wcs = WCS(HDUlist[0].header)
print(HDUlist[1].header)
print(gal_wcs)
HDUlist.close()

fig_apl = aplpy.FITSFigure(ROOTfile)
fig_apl.show_colorscale(cmap='hot',smooth=5)
col_bar = fig_apl.show_colorbar()
#col_bar

skymap_data, skymap_header = fits.getdata(ROOTfile, header = True)

excess_data_N = numpy.nan_to_num(skymap_data) # rm nan's for fitting<jupyter_output><empty_output><jupyter_text>Create coordinate arrays<jupyter_code>x_bins = skymap_data.shape[0]
y_bins = skymap_data.shape[1]

# Create x and y indices
x = numpy.linspace(0, x_bins-1, x_bins)
y = numpy.linspace(0, y_bins-1, y_bins)
x, y = numpy.meshgrid(x, y)
coords = x, y

x_world, y_world = fig_apl.pixel2world(x,y)
coords_world = x_world, y_world

print(fig_apl.world2pixel(83.6296,22.014)) # Crab<jupyter_output>(129.14929118926517, 128.97964846247766)<jupyter_text>Fitting Central Gaussian<jupyter_code>#Crab_bounds = {'x_mean':(115,117),'y_mean':(115,117)}
Crab_gauss = models.Gaussian2D(amplitude=25,x_mean=129,y_mean=129,x_stddev=2.,y_stddev=2.)

fit_p = fitting.LevMarLSQFitter()
fit_gauss = fit_p(Crab_gauss, x, y, excess_data_N, maxiter=2500)

print(fit_gauss.parameters)

Crab_gauss_1 = models.Gaussian2D(amplitude=20,x_mean=129,y_mean=129,x_stddev=1.,y_stddev=1.)
Crab_gauss_2 = models.Gaussian2D(amplitude=10,x_mean=129,y_mean=129,x_stddev=3.,y_stddev=3.)
Crab_gauss_sum = Crab_gauss_1 + Crab_gauss_2 

sum_fit = fitting.LevMarLSQFitter()
fit_sum = sum_fit(Crab_gauss_sum, x, y, excess_data_N, maxiter=2500)

print(fit_sum.parameters)


print(fig_apl.pixel2world(fit_gauss.x_mean.value,fit_gauss.y_mean.value))

stddevs = fig_apl.pixel2world(fit_gauss.x_stddev.value,fit_gauss.y_stddev.value)

pointsource_filename = "./Crab_gaussian_source.fits"
pointsource_excess = fit_sum(x,y)

if not os.path.exists(pointsource_filename):
    fits.writeto(pointsource_filename, pointsource_excess, skymap_header, clobber=True)

fig_smoothed = pyplot.figure()
pyplot.imshow(numpy.nan_to_num(pointsource_excess))
pyplot.colorbar()
numpy.nanmax(pointsource_excess)

fig_apl = aplpy.FITSFigure(pointsource_filename)
fig_apl.show_colorscale(vmin=-1,vmax=30)<jupyter_output>WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / WCS for this file 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING:astropy:FITSFixedWarning: RADECSYS= 'FK5 ' / WCS for this file 
the RADECSYS keyword is deprecated, use RADESYSa.<jupyter_text>Excess after subtraction of central sources<jupyter_code>residual_excess = skymap_data - pointsource_excess 

subtracted_filename = "./Crab_subtracted.fits"

#if not os.path.exists(pointsource_filename):
fits.writeto(subtracted_filename, numpy.nan_to_num(residual_excess), skymap_header, clobber=True)

fig_smoothed = pyplot.figure()
pyplot.imshow(numpy.nan_to_num(residual_excess))
numpy.nanmax(residual_excess)
pyplot.colorbar()
fig_subtracted = aplpy.FITSFigure(subtracted_filename)
fig_subtracted.show_colorscale(smooth=5)
fig_subtracted.show_colorbar()<jupyter_output>WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / WCS for this file 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING:astropy:FITSFixedWarning: RADECSYS= 'FK5 ' / WCS for this file 
the RADECSYS keyword is deprecated, use RADESYSa.
INFO:astropy:Auto-setting vmin to -1.081e+00
INFO:astropy:Auto-setting vmax to  2.741e+00<jupyter_text>Gaussian Smoothing the Residual<jupyter_code>smooth_stddev = (x_bins/4.0)*0.2 # degrees to bins
smoothed_residual_excess = numpy.array(residual_excess, dtype=numpy.float128)

for i in range(0,x_bins):
    for j in range(0,y_bins):
        #print(residual_excess[i,j],end="")
        if math.isinf(residual_excess[i,j]):
            print("inf: %d,%d" % (i,j))
           
        
print(residual_excess[129,129])
        
smooth_stddev = (x_bins/4.0)*0.2 # degrees to bins
smoothed_residual_excess = numpy.array(residual_excess, dtype=numpy.float128)

for i in range(0,x_bins):
    for j in range(0,y_bins):
        #print(residual_excess[i,j],end="")
        if not math.isnan(residual_excess[i,j]) and not math.isinf(residual_excess[i,j]):
            #print("%d,%d" % (i,j))
            point_gauss = models.Gaussian2D(amplitude=residual_excess[i,j],x_mean=i,y_mean=j,x_stddev=smooth_stddev,y_stddev=smooth_stddev)
            smoothed_residual_excess += point_gauss(x,y)
            del point_gauss
        
smoothed_residual_excess_N = numpy.nan_to_num(smoothed_residual_excess)
smoothed_residual_excess_F = numpy.array(smoothed_residual_excess_N, dtype=numpy.float)
print(smoothed_residual_excess_F)
print(smoothed_residual_excess_F[129,129])

fig_smoothed = pyplot.figure()
pyplot.imshow(smoothed_residual_excess_F)


def symmetric_2d_gaussian(x, y, amp, x0, y0, sigma):
    g = amp*numpy.exp(-((x-x0)**2+(y-y0)**2)/(2*sigma**2))
    return g

x_max = fig_apl.pixel2world(0,y_bins/2)[0]
x_min = fig_apl.pixel2world(x_bins,y_bins/2)[0]
delta_x = x_max - x_min
print(delta_x)

y_min = fig_apl.pixel2world(x_bins/2,0)[1]
y_max = fig_apl.pixel2world(x_bins/2,y_bins)[1]
delta_y = y_max - y_min
print(y_max-y_min)

sigma = 0.2
sigma_bins = sigma/delta_y * x_bins 
print(sigma_bins)

amplitude = 1
x_lower = -0.5
x_upper = 0.5
area = integrate.dblquad(symmetric_2d_gaussian, x_lower, x_upper, lambda y_lower: -0.5, lambda y_upper: 0.5, args=(amplitude,0,0,sigma_bins))

area2 = integrate.dblquad(symmetric_2d_gaussian, x_lower, x_upper, lambda y_lower: -0.5, lambda y_upper: 0.5, args=(amplitude,0,0,1))
area3 = integrate.dblquad(symmetric_2d_gaussian, x_lower, x_upper, lambda y_lower: -0.5, lambda y_upper: 0.5, args=(amplitude,0,0,8))

x = lambda y: 5
print(x)
print(sigma_bins)
print(area)
print(area2)
print(area3)<jupyter_output><empty_output><jupyter_text>create symmetric gaussian kernel<jupyter_code>sig = 8
n_dim = 35 # should be odd and greater than 1 
gauss_kernel = numpy.zeros((n_dim, n_dim))

#coord_range = range(math.ceil(-n_dim/2),math.ceil(n_dim/2))
kernel_range = range(0, n_dim)

for i in kernel_range:
    for j in kernel_range:
        area = integrate.dblquad(symmetric_2d_gaussian,i,i+1,lambda yl:j,lambda yu:j+1,(1,n_dim/2,n_dim/2,sig))
        gauss_kernel[i, j] = area[0]
        
print(gauss_kernel)
d = math.floor(n_dim/2)

residual_excess_n = numpy.nan_to_num(residual_excess)
smoothed_residual_excess_conv = numpy.zeros((residual_excess_n.shape))

for i in range(d,residual_excess_n.shape[0]-d):
    for j in range(d,residual_excess_n.shape[1]-d):
        # sum over the other cells 
        for ip in kernel_range:
            for jp in kernel_range: 
                smoothed_residual_excess_conv[i,j] += gauss_kernel[ip,jp] * residual_excess_n[i+d-ip,j+d-jp]
                
                

fig_custom = pyplot.figure()
pyplot.imshow(smoothed_residual_excess_conv)
pyplot.colorbar()<jupyter_output><empty_output><jupyter_text>fix header<jupyter_code>print(type(HDUlist[0].header))
fixed_header = HDUlist[1].header.copy()
fixed_header['RADECSYS']
verify_header(fixed_header)<jupyter_output><empty_output><jupyter_text>King Function<jupyter_code>models.custom_model?
def KingFunction(x, sigma, gamma):
    K = 1/(2*math.pi*sigma**2)*(1-1/gamma)*math.pow((1+1/(2*gamma)*x**2/sigma**2),-gamma)
    return K

king_custom = models.custom_model(func=KingFunction)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1271598272138229, 'fraction_numerical': 0.034152267818574514, 'mean_word_length': 6.124038461538461, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '343582', 'license': 'MIT', 'n_tokens_mistral': 3328, 'n_tokens_neox': 3086, 'n_words': 440, 'path': 'python/Crab_Gaussian.ipynb', 'repo_name': 'mbuchove/VEGAS_scripts', 'size': '7408'}