Page MenuHomec4science
Files F97581409

mockimgs_pkl_compute_area
No OneTemporary
Actions

File Metadata

Created
Sun, Jan 5, 11:03

mockimgs_pkl_compute_area
View Options

#!/usr/bin/env python3
import numpy as np
from astropy import constants as cte
from astropy import units as u
import argparse
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from pNbody import *
from pNbody import ic
from pNbody.Mockimgs import lib as libMockimgs
import pickle
####################################################################
# option parser
####################################################################
description=""
epilog =""""""
parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument(action="store",
dest="files",
metavar='FILE',
type=str,
default=None,
nargs='*',
help='a file name')
parser.add_argument("-o",
action="store",
type=str,
dest="outputfilename",
default=None,
help="Name of the output file")
parser.add_argument("--magmin",
action="store",
type=float,
dest="magmin",
default=26,
help="minimum surface brightness")
parser.add_argument("--magmax",
action="store",
type=float,
dest="magmax",
default=31,
help="maximal surface brightness")
parser.add_argument("--nbins",
action="store",
type=int,
dest="nbins",
default=50,
help="number of bins")
parser.add_argument("--diff",
action="store_true",
dest="diff",
default=False,
help="substract the area of at the minimum magnitude")
####################################################################
# main
####################################################################
cs = ['r','g','b','y']
if __name__ == '__main__':
opt = parser.parse_args()
nmagbins = opt.nbins
magmin = opt.magmin
magmax = opt.magmax
magbins = np.linspace(magmin,magmax,nmagbins)
for j,f in enumerate(opt.files):
print(f)
fd = open(f,"rb")
images = pickle.load(fd)
fd.close()
# convert to a surface density
for i in range(len(images)):
FluxMap,params = images[i]
# apply some filters
FluxMap = libMockimgs.FilterFluxMap(FluxMap,None)
# compute the surface brightness map from the luminosity/flux map
images[i] = libMockimgs.FluxToSurfaceBrightness(FluxMap,params["object_distance"],params["ccd_pixel_area"])
# loop over images of the same model
Npixs = np.zeros((len(images),nmagbins))
for k in range(len(images)):
image = images[k]
image = np.ravel(image)
image = np.compress(image<100,image)
Npix = np.zeros(nmagbins)
for i in tqdm(range(nmagbins)):
img1 = np.where(image<magbins[i],1,0)
N = np.ravel(img1).sum()
if i==0:
N0 = N
Npix[i] = N
if opt.diff:
Npix[i] = Npix[i] - N0
#plt.plot(magbins,Npix,color=cs[j])
# sum
Npixs[k] = Npix
mean = Npixs.mean(axis=0)
std = Npixs.std(axis=0)
label = os.path.basename(f)
label = label[:label.find("keV")+3]
#label = label[:-4]
plt.plot(magbins,mean,color=cs[j],label=r"$%s$"%label,lw=2)
#plt.plot(magbins,mean+1*std,color=cs[j],lw=5)
#plt.plot(magbins,mean-1*std,color=cs[j],lw=5)
ax = plt.gca()
ax.fill_between(magbins,mean+1*std,mean-1*std,facecolor=cs[j],alpha=0.5)
plt.legend()
plt.xlabel(r"$\rm{VIS}\,\,\rm{magnitude}$")
plt.ylabel(r"$\rm{Area\,\,\rm{pixels}}$")
if opt.outputfilename is None:
plt.show()
else:
plt.savefig(opt.outputfilename)

Event Timeline

c4science · Help