import copy
import os
import numpy as np
from astropy.io import fits
from scipy.optimize import bisect
from scipy.special import gammainc, gammaincinv
from dataclasses import dataclass, field
from galfitools.galin.getStar import GetFits
from typing import Dict
[docs]
@dataclass
class GalHead:
"""
Data class to store the header info of the galfit file
"""
inputimage = "none.fits"
outimage = "none-out.fits"
sigimage = "none"
psfimage = "none"
psfsamp = 1
maskimage = "none"
constraints = "none"
xmin = 0
xmax = 1
ymin = 0
ymax = 1
convx = 1
convy = 1
mgzpt = 25.0
scale = 1.0
scaley = 1.0
display = "regular"
P = 0
# extra info
imgidx = "sci"
flagidx = False
num = 1
flagnum = False
exptime = 1
stampmask = "stampmask.fits"
[docs]
@dataclass
class GalComps:
"""
Data class to store the galfit components of the galfit file
"""
N: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
NameComp = np.array([]) # 0)
PosX = np.array([]) # 1)
PosY = np.array([]) # 2)
Mag = np.array([]) # 3)
Rad = np.array([]) # 4)
Exp = np.array([]) # 5)
Exp2 = np.array([]) # 6) for moffat or Ferrer
Exp3 = np.array([]) # 7) for moffat
# 8) There is No 8 in any galfit model
AxRat = np.array([]) # 9) AxisRatio
PosAng = np.array([]) # 10) position angle
skip: np.array = field(
default_factory=lambda: np.array([], dtype=np.int64)
) # z) skip model
Active = np.array([]) # activate component for galaxy
# store the flags related to parameters
PosXFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
PosYFree: np.ndarray = field(default_factory=lambda: np.array([], dtype=np.int64))
PosYFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
MagFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
RadFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
ExpFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
Exp2Free: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
Exp3Free: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
AxRatFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
PosAngFree: np.array = field(default_factory=lambda: np.array([], dtype=np.int64))
# computed parameters:
Rad50 = np.array([])
SerInd = np.array([])
Rad50kpc = np.array([])
Rad50sec = np.array([])
Rad90 = np.array([])
AbsMag = np.array([])
Lum = np.array([])
Flux = np.array([])
PerLight = np.array([])
me = np.array([])
mme = np.array([])
kser = np.array([])
KronRad = np.array([])
PetRad = np.array([])
[docs]
@dataclass
class GalSky:
"""
Data class to store the sky component of the GALFIT file
"""
sky = 0
dskyx = 0
dskyy = 0
skip = 0
skyfree = 1
dskyxfree = 0
dskyyfree = 0
[docs]
@dataclass
class DataImg:
"""
Data class to store the galaxy, model and residual images
of the output file of GALFIT
"""
img = np.array([[1, 1], [1, 1]])
model = np.array([[1, 1], [1, 1]])
imres = np.array([[1, 1], [1, 1]])
mask = np.array([[1, 1], [1, 1]])
sigma = np.array([[1, 1], [1, 1]])
imsnr = np.array([[1, 1], [1, 1]])
imchi = np.array([[1, 1], [1, 1]])
impsf = np.array([[1, 1], [1, 1]])
[docs]
def readDataImg(galhead, sigma=None):
"""
reads galfit output cube data image
"""
dataimg = DataImg()
# reading galaxy and model images from file
errmsg = "file {} does not exist".format(galhead.outimage)
assert os.path.isfile(galhead.outimage), errmsg
# hdu 1 => image hdu 2 => model
hdu = fits.open(galhead.outimage)
dataimg.img = (hdu[1].data.copy()).astype(float)
dataimg.model = (hdu[2].data.copy()).astype(float)
dataimg.imres = (hdu[3].data.copy()).astype(float)
hdu.close()
### reading mask image from file
if os.path.isfile(galhead.maskimage):
GetFits(
galhead.maskimage,
galhead.stampmask,
0,
galhead.xmin,
galhead.xmax,
galhead.ymin,
galhead.ymax,
)
errmsg = "file {} does not exist".format(galhead.stampmask)
assert os.path.isfile(galhead.stampmask), errmsg
hdu = fits.open(galhead.stampmask)
mask = hdu[0].data
dataimg.mask = np.array(mask, dtype=bool)
hdu.close()
else:
dataimg.mask = None
if sigma is not None:
hdu = fits.open(sigma)
dataimg.sigma = hdu[0].data
hdu.close()
return dataimg
[docs]
class Galfit:
"""
Class to read GALFIT parameters file
Attributes
----------
File: str
name of the GALFIT file
Methods
-------
ReadHead : GalHead
reads the header of the GALFIT file
ReadComps : GalComps
reads the galfit components of the GALFIT file
ReadSky : GalSky
reads the galfit sky components of the GALFIT file
"""
def __init__(self, File: str):
self.File = File
[docs]
def ReadHead(self) -> GalHead:
inputf = self.File
galhead = GalHead() # class for header
GalfitFile = open(inputf, "r")
# All lines including the blank ones
lines = (line.rstrip() for line in GalfitFile)
lines = (line.split("#", 1)[0] for line in lines) # remove comments
# remove lines containing only comments
lines = (line.rstrip() for line in lines)
lines = (line for line in lines if line) # Non-blank lines
lines = list(lines)
index = 0
while index < len(lines):
# ================================================================
# IMAGE and GALFIT CONTROL PARAMETERS
# A) tempfits/A2399-3-2.fits # Input data image (FITS file)
# B) A2399-215615.96-m073822.7-337-out.fits # Output data image block
# C) tempfits/none-3-2.fits # Sigma image name
# D) psfs/PSF-1309-721.fits # Input PSF image
# E) 1 # PSF fine sampling factor relative to data
# F) mask-337 # Bad pixel mask (FITS image or ASCII coord list)
# G) constraints # File with parameter constraints (ASCII file)
# H) 129 809 265 809 # Image region to fit (xmin xmax ymin ymax)
# I) 60 60 # Size of the convolution box (x y)
# J) 21.672 # Magnitude photometric zeropoint
# K) 0.680 0.680 # Plate scale (dx dy) [arcsec per pixel]
# O) regular # Display type (regular, curses, both)
# P) 0 # Choose: 0=optimize,1=model,2=imgblock,3=subcomps
line = lines[index]
(tmp) = line.split()
if len(tmp) > 1: # avoids empty options
if tmp[0] == "A)": # input image
galhead.inputimage = tmp[1]
if tmp[0] == "B)": # out image
galhead.outimage = tmp[1]
if tmp[0] == "C)": # sigma
galhead.sigimage = tmp[1]
if tmp[0] == "D)": # psf file
galhead.psfimage = tmp[1]
if tmp[0] == "E)": # psf sampling
galhead.psfsamp = int(tmp[1])
if tmp[0] == "F)": # mask image
try:
galhead.maskimage = tmp[1]
except IndexError:
galhead.maskimage = "None"
if tmp[0] == "G)": # psf file
galhead.constraints = tmp[1]
if tmp[0] == "H)": # region fit box
galhead.xmin = int(tmp[1])
galhead.xmax = int(tmp[2])
galhead.ymin = int(tmp[3])
galhead.ymax = int(tmp[4])
if tmp[0] == "I)": # convolution size
galhead.convx = int(tmp[1])
try:
galhead.convy = int(tmp[2])
except Exception:
galhead.convy = galhead.convx
if tmp[0] == "J)": # mgzpt
galhead.mgzpt = float(tmp[1])
if tmp[0] == "K)": # plate scale
galhead.scale = float(tmp[1])
try:
galhead.scaley = float(tmp[2])
except Exception:
galhead.scaley = galhead.scale
if tmp[0] == "O)": # display
galhead.display = tmp[1]
if tmp[0] == "P)": # optional output
galhead.P = int(tmp[1])
index += 1
# check for extension in name
chars = set("[]")
numbers = set("1234567890")
if any((c in chars) for c in galhead.inputimage):
print("Ext Found")
galhead.flagidx = True
(filename, imgidxc) = galhead.inputimage.split("[")
(imgidx, trash) = imgidxc.split("]")
if any((n in numbers) for n in imgidx):
galhead.flagnum = True
(imgidx, num) = imgidx.split(",")
num = int(num)
galhead.inputimage = filename
galhead.imgidx = imgidx
galhead.num = num
####################
return galhead
[docs]
def ReadComps(self) -> GalComps:
File = self.File
galcomps = GalComps()
GalfitFile = open(File, "r")
# All lines including the blank ones
lines = (line.rstrip() for line in GalfitFile)
lines = (line.split("#", 1)[0] for line in lines) # remove comments
# remove lines containing only comments
lines = (line.rstrip() for line in lines)
lines = (line for line in lines if line) # Non-blank lines
lines = list(lines)
index = 0
N = 0
while index < len(lines):
line = lines[index]
(tmp) = line.split()
# init values
NameComp = "none"
PosX = 0
PosY = 0
Mag = 99
Rad = 0
Exp = 0
Exp2 = 0
Exp3 = 0
AxRat = 1
PosAng = 0
skip = 0
flagcomp = False
Xfree = 1
Yfree = 1
Magfree = 1
Radfree = 1
Expfree = 1
Exp2free = 1
Exp3free = 1
AxRatfree = 1
PosAngfree = 1
if (tmp[0] == "0)") and (tmp[1] != "sky"):
namec = tmp[1]
N = N + 1
NameComp = namec
while tmp[0] != "Z)":
index += 1
line = lines[index]
(tmp) = line.split()
if tmp[0] == "1)": # center
PosX = float(tmp[1])
PosY = float(tmp[2])
Xfree = int(tmp[3])
Yfree = int(tmp[4])
if tmp[0] == "3)": # mag
Mag = float(tmp[1])
Magfree = int(tmp[2])
if tmp[0] == "4)":
Rad = float(tmp[1])
Radfree = int(tmp[2])
if tmp[0] == "5)":
Exp = float(tmp[1])
Expfree = int(tmp[2])
if tmp[0] == "6)":
Exp2 = float(tmp[1])
try:
Exp2free = int(tmp[2])
except Exception:
Exp2free = 0
if tmp[0] == "7)":
Exp3 = float(tmp[1])
try:
Exp3free = int(tmp[2])
except Exception:
Exp3free = 0
if tmp[0] == "9)":
AxRat = float(tmp[1])
try:
AxRatfree = int(tmp[2])
except Exception:
AxRatfree = 0
if tmp[0] == "10)":
PosAng = float(tmp[1])
PosAngfree = int(tmp[2])
if tmp[0] == "Z)":
skip = int(tmp[1])
galcomps.PosX = np.append(galcomps.PosX, PosX)
galcomps.PosY = np.append(galcomps.PosY, PosY)
galcomps.NameComp = np.append(galcomps.NameComp, NameComp)
galcomps.N = np.append(galcomps.N, N)
galcomps.Mag = np.append(galcomps.Mag, Mag)
galcomps.Rad = np.append(galcomps.Rad, Rad)
galcomps.Exp = np.append(galcomps.Exp, Exp)
galcomps.Exp2 = np.append(galcomps.Exp2, Exp2)
galcomps.Exp3 = np.append(galcomps.Exp3, Exp3)
galcomps.AxRat = np.append(galcomps.AxRat, AxRat)
galcomps.PosAng = np.append(galcomps.PosAng, PosAng)
galcomps.skip = np.append(galcomps.skip, skip)
galcomps.Active = np.append(galcomps.Active, flagcomp)
galcomps.PosXFree = np.append(galcomps.PosXFree, Xfree)
galcomps.PosYFree = np.append(galcomps.PosYFree, Yfree)
galcomps.MagFree = np.append(galcomps.MagFree, Magfree)
galcomps.RadFree = np.append(galcomps.RadFree, Radfree)
galcomps.ExpFree = np.append(galcomps.ExpFree, Expfree)
galcomps.Exp2Free = np.append(galcomps.Exp2Free, Exp2free)
galcomps.Exp3Free = np.append(galcomps.Exp3Free, Exp3free)
galcomps.AxRatFree = np.append(galcomps.AxRatFree, AxRatfree)
galcomps.PosAngFree = np.append(galcomps.PosAngFree, PosAngfree)
index += 1
GalfitFile.close()
# filling the rest of arrays with zeros:
arsiz = len(galcomps.N)
galcomps.Rad50 = np.zeros(arsiz)
galcomps.SerInd = np.zeros(arsiz)
galcomps.Rad50kpc = np.zeros(arsiz)
galcomps.Rad50sec = np.zeros(arsiz)
galcomps.Rad90 = np.zeros(arsiz)
galcomps.AbsMag = np.zeros(arsiz)
galcomps.Lum = np.zeros(arsiz)
galcomps.Flux = np.zeros(arsiz)
galcomps.PerLight = np.zeros(arsiz)
galcomps.me = np.zeros(arsiz)
galcomps.mme = np.zeros(arsiz)
galcomps.kser = np.zeros(arsiz)
galcomps.KronRad = np.zeros(arsiz)
galcomps.PetRad = np.zeros(arsiz)
return galcomps
[docs]
def ReadSky(self) -> GalSky:
File = self.File
galsky = GalSky()
GalfitFile = open(File, "r")
# All lines including the blank ones
lines = (line.rstrip() for line in GalfitFile)
lines = (line.split("#", 1)[0] for line in lines) # remove comments
# remove lines containing only comments
lines = (line.rstrip() for line in lines)
lines = (line for line in lines if line) # Non-blank lines
lines = list(lines)
index = 0
while index < len(lines):
line = lines[index]
(tmp) = line.split()
# init values
NameComp = "sky"
sky = 0
dskyx = 0
dskyy = 0
skip = 0
skyfree = 1
dskyxfree = 0
dskyyfree = 0
if (tmp[0] == "0)") and (tmp[1] == NameComp):
while tmp[0] != "Z)":
index += 1
line = lines[index]
(tmp) = line.split()
if tmp[0] == "1)": # center
sky = float(tmp[1])
skyfree = int(tmp[2])
if tmp[0] == "2)": # mag
dskyx = float(tmp[1])
dskyxfree = int(tmp[2])
if tmp[0] == "3)":
dskyy = float(tmp[1])
dskyyfree = int(tmp[2])
if tmp[0] == "Z)":
skip = int(tmp[1])
# note if this work in this way or using append
galsky.NameComp = NameComp
galsky.sky = sky
galsky.dskyx = dskyx
galsky.dskyy = dskyy
galsky.skip = skip
galsky.skyfree = skyfree
galsky.dskyxfree = dskyxfree
galsky.dskyyfree = dskyyfree
index += 1
GalfitFile.close()
return galsky
[docs]
def numParFree(galcomps: GalComps) -> int:
"""Obtains the number of free parameters
Count the number of free parameters of the surface brightness model
of the GALFIT file.
Parameters
----------
galcomps : GalComps Data Class defined above
Returns
-------
pt : int
Number of free parameters
See Also
--------
numParSkyFree: Obtains the number of free parameters of the sky
component
Notes
-----
It only counts the components with galcomps.Active == True
This prevents to count components of other galaxies that
were simultanously fitted.
Sky component is ignored.
"""
p1 = 0
p2 = 0
p3 = 0
p4 = 0
p5 = 0
p6 = 0
p7 = 0
p8 = 0
p9 = 0
parmask1 = (galcomps.Active == 1) & (galcomps.PosXFree == 1)
parmask2 = (galcomps.Active == 1) & (galcomps.PosYFree == 1)
parmask3 = (galcomps.Active == 1) & (galcomps.MagFree == 1)
parmask4 = (galcomps.Active == 1) & (galcomps.RadFree == 1)
parmask5 = (galcomps.Active == 1) & (galcomps.ExpFree == 1)
parmask6 = (galcomps.Active == 1) & (galcomps.Exp2Free == 1)
parmask7 = (galcomps.Active == 1) & (galcomps.Exp3Free == 1)
parmask8 = (galcomps.Active == 1) & (galcomps.AxRatFree == 1)
parmask9 = (galcomps.Active == 1) & (galcomps.PosAngFree == 1)
if parmask1.any():
p1 = np.sum(galcomps.PosXFree[parmask1])
if parmask2.any():
p2 = np.sum(galcomps.PosYFree[parmask2])
if parmask3.any():
p3 = np.sum(galcomps.MagFree[parmask3])
if parmask4.any():
p4 = np.sum(galcomps.RadFree[parmask4])
if parmask5.any():
p5 = np.sum(galcomps.ExpFree[parmask5])
if parmask6.any():
p6 = np.sum(galcomps.Exp2Free[parmask6])
if parmask7.any():
p7 = np.sum(galcomps.Exp3Free[parmask7])
if parmask8.any():
p8 = np.sum(galcomps.AxRatFree[parmask8])
if parmask9.any():
p9 = np.sum(galcomps.PosAngFree[parmask9])
pt = p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
return int(pt)
[docs]
def numParSkyFree(galsky: GalSky) -> int:
"""Obtains the number of free parameters of the sky component
Count the number of free parameters of the sky component
of the GALFIT file.
Parameters
----------
galsky: GalSky data Class defined above
Returns
-------
pt : int
Number of free parameters
See Also
--------
numParFree : Obtains the number of free parameters
"""
p1 = 0
p2 = 0
p3 = 0
parmask1 = galsky.skyfree == 1
parmask2 = galsky.dskyxfree == 1
parmask3 = galsky.dskyyfree == 1
if parmask1:
p1 = 1
if parmask2:
p2 = 1
if parmask3:
p3 = 1
pt = p1 + p2 + p3
return int(pt)
[docs]
def numComps(galcomps: GalComps, name: str) -> int:
"""gets the number of components for a galaxy.
Given the number of components of a GALFIT file,
it discerns which components belong to the galaxy. This
is useful when simultaneous fitting was used.
Parameters
----------
galcomps : GalComps data class defined above
name : str
indicates which model components to count:
all, sersic, expdisk, gaussian, devauc
Returns
-------
N : Number of components of the galaxy
"""
if name == "all":
nummask = (galcomps.Active == 1) & (
(galcomps.NameComp == "sersic")
| (galcomps.NameComp == "expdisk")
| (galcomps.NameComp == "gaussian")
| (galcomps.NameComp == "devauc")
| (galcomps.NameComp == "ferrer")
)
else:
nummask = (galcomps.Active == 1) & (galcomps.NameComp == name)
N = galcomps.Active[nummask].size
return N
[docs]
def SelectGal(galcomps: GalComps, distmax: float, n_comp: int) -> GalComps:
"""Selects which components belong to the galaxy.
From the data class GalComps, it changes the flag Active to
True for those components where distance is less than distmax
from component n_comp.
This is useful when there are simultaneous
fitting of galaxies in a single GALFIT file.
Parameters
----------
galcomps : GalComps data class defined above
distmax : float
maximum distance among components
n_comp : int
The center of the galaxy will be determined from this
component (n_comp). The distances of the remaining
components will be calculated relative to this center.
Returns
-------
galcomps : GalComps data class with flag Active = True for
those components which belong to the galaxy selected
by n_comp.
Notes
------
If a galaxy is made up of three components, it does not matter which
of the three you select as a center. As long these belong to the
same galaxy (and center).
"""
galcomps.Active.fill(False)
idx = np.where(galcomps.N == n_comp)
assert idx[0].size != 0, "component not found"
n_idx = idx[0].item(0)
if distmax > 10:
print("Warning: maximum distance among components is greater than 10")
print(
"Equations' solutions only apply for components that share the same center"
)
galcomps.Active[n_idx] = True # main component
posx = galcomps.PosX[n_idx]
posy = galcomps.PosY[n_idx]
dx = galcomps.PosX - posx
dy = galcomps.PosY - posy
dist = np.sqrt((dx) ** 2 + (dy) ** 2)
maskdist = dist <= distmax
galcomps.Active[maskdist] = True
return galcomps
[docs]
def SelectComp(galcomps: GalComps, n_comp: int) -> GalComps:
"""Selects only a component for galcomps.
From the data class GalComps, it changes the flag Active to
True for one components indicated in n_comp.
useful for selecting for one component
Parameters
----------
galcomps : GalComps data class defined above
n_comp : int
the component to select to be active
Returns
-------
galcomps : GalComps data class with flag Active = True for
one component selected
by n_comp.
"""
comps = copy.deepcopy(galcomps)
comps.Active.fill(False)
idx = np.where(comps.N == n_comp)
assert idx[0].size != 0, "component not found"
n_idx = idx[0].item(0)
comps.Active[n_idx] = True
return comps
[docs]
def conver2Sersic(galcomps: GalComps) -> GalComps:
"""Converts the exponential, gaussian or De Vaucouleurs to Sersic
Using the format of GALFIT with the GalComps data class,
it converts the parameters of the functions exponential,
gaussian or De Vaucouleurs to parameters of Sersic function.
Parameters
-----------
galcomps : GalComps data class defined above
Returns
-------
galcomps : GalComps data class with the functions converted to Sersic
parameters
"""
comps = copy.deepcopy(galcomps)
maskdev = comps.NameComp == "devauc"
maskexp = comps.NameComp == "expdisk"
maskgas = comps.NameComp == "gaussian"
K_GAUSS = 0.6931471805599455 # constant k for gaussian
K_EXP = 1.6783469900166612 # constant k for expdisk
SQ2 = np.sqrt(2)
# for gaussian functions
if maskgas.any():
comps.NameComp[maskgas] = "sersic"
comps.Exp[maskgas] = 0.5
comps.Rad[maskgas] = comps.Rad[maskgas] / 2.354 # converting to sigma
comps.Rad[maskgas] = (
SQ2 * (K_GAUSS**0.5) * comps.Rad[maskgas]
) # converting to Re
# for de vaucouleurs
if maskdev.any():
comps.NameComp[maskdev] = "sersic"
comps.Exp[maskdev] = 4
# for exponential disks
if maskexp.any():
comps.NameComp[maskexp] = "sersic"
comps.Exp[maskexp] = 1
comps.Rad[maskexp] = K_EXP * comps.Rad[maskexp] # converting to Re
return comps
[docs]
def conver2Ferrer(galcomps: GalComps, scale, N) -> GalComps:
"""Converts the Sersic to Ferrer
Using the format of GALFIT with the GalComps data class,
it converts the parameters of the functions Sersic n = 0.5,
to parameters of Ferrer function.
Parameters
-----------
galcomps : GalComps data class defined above
scale : Plate Scale ''/pixel
N : number of component to be converted
Returns
-------
galcomps : GalComps data class with the functions converted to Ferrer
parameters
"""
comps = copy.deepcopy(galcomps)
K_GAUSS = 0.6931471805599455 # constant k for gaussian
SQ2 = np.sqrt(2)
masksec = (comps.NameComp == "sersic") * (comps.N == N)
rout = comps.Rad[masksec] / np.sqrt(K_GAUSS)
rout_arcsec = rout * scale
I0 = comps.Flux[masksec] / (np.pi * rout_arcsec**2)
comps.NameComp[masksec] = "ferrer"
comps.Mag[masksec] = -2.5 * np.log10(I0)
comps.Rad[masksec] = rout
comps.Exp[masksec] = 1
comps.Exp2[masksec] = 0
return comps
[docs]
def expdisk_to_edgedisk(
mag_total: float,
rs_pix: float,
axis_ratio: float,
pa_deg: float,
pixscale: float,
) -> Dict[str, float]:
"""
Convert approximate GALFIT `expdisk` parameters into initial `edgedisk`
parameters.
This conversion is based on two assumptions:
1. The radial scale length is preserved:
rs_edge = rs_exp
2. The vertical scale height is approximated from the exponential disk
axis ratio:
hs ~= q * rs
The central surface brightness of the edge-on disk is then estimated by
conserving the total flux:
mu0 = mag_total + 2.5 * log10(2 * pi * rs * hs * pixscale^2)
where:
- rs and hs are in pixels
- pixscale is in arcsec/pixel
- mu0 is returned in mag/arcsec^2
Parameters
----------
mag_total : float
Total magnitude of the `expdisk` model.
rs_pix : float
Exponential scale length in pixels.
axis_ratio : float
Axis ratio b/a of the `expdisk` model.
pa_deg : float
Position angle in degrees.
pixscale : float
Pixel scale in arcsec/pixel.
Returns
-------
Dict[str, float]
Dictionary with approximate `edgedisk` parameters:
- ``mu0_mag_arcsec2`` : central surface brightness [mag/arcsec^2]
- ``hs_pix`` : disk scale height [pixels]
- ``rs_pix`` : disk scale length [pixels]
- ``pa_deg`` : position angle [degrees]
Raises
------
ValueError
If any input value is not physically valid.
Notes
-----
This is an approximation for initial guesses only. It is most useful when
the galaxy is close to edge-on. For moderately inclined systems, the
relation hs ~= q * rs is only a rough estimate.
"""
if rs_pix <= 0:
raise ValueError("rs_pix must be > 0.")
if not (0 < axis_ratio <= 1):
raise ValueError("axis_ratio must satisfy 0 < axis_ratio <= 1.")
if pixscale <= 0:
raise ValueError("pixscale must be > 0.")
hs_pix = axis_ratio * rs_pix
if hs_pix <= 0:
raise ValueError("Computed hs_pix must be > 0.")
area_arcsec2 = 2.0 * np.pi * rs_pix * hs_pix * pixscale**2
mu0_mag_arcsec2 = mag_total + 2.5 * np.log10(area_arcsec2)
return {
"mu0_mag_arcsec2": mu0_mag_arcsec2,
"hs_pix": hs_pix,
"rs_pix": rs_pix,
"pa_deg": pa_deg,
}
[docs]
def conver2Edge(galcomps: GalComps, scale, N) -> GalComps:
"""Converts the expdisk to edgedisk
Using the format of GALFIT with the GalComps data class,
it converts the parameters of the functions Sersic n = 1 ,
to parameters of edgedisk function.
Parameters
-----------
galcomps : GalComps data class defined above
scale : Plate Scale ''/pixel
N : number of component to be converted
Returns
-------
galcomps : GalComps data class with the functions converted to Ferrer
parameters
"""
comps = copy.deepcopy(galcomps)
R_EXP = 1.6783469900166612 # constant for scale length convertion
masksec = (comps.NameComp == "sersic") * (comps.N == N)
mag_total = comps.Mag[masksec]
rs_pix = comps.Rad[masksec] / R_EXP
axis_ratio = comps.AxRat[masksec]
pa_deg = comps.PosAng[masksec]
pixscale = scale # arcsec/pixel, example only
params = expdisk_to_edgedisk(
mag_total=mag_total,
rs_pix=rs_pix,
axis_ratio=axis_ratio,
pa_deg=pa_deg,
pixscale=pixscale,
)
comps.NameComp[masksec] = "edgedisk"
comps.Mag[masksec] = params["mu0_mag_arcsec2"]
comps.Rad[masksec] = params["hs_pix"]
comps.Exp[masksec] = params["rs_pix"]
comps.PosAng[masksec] = params["pa_deg"]
comps.AxRat[masksec] = 1
comps.AxRatFree[masksec] = -1
comps.ExpFree[masksec] = 1
return comps
[docs]
def GetRadAng(R: float, q: list, pa: list, theta: float) -> float:
"""Obtains the radius along the specified angular direction.
Given the angular position and axis ratio of the ellipse
it gives the radius of the ellipse in the theta direction.
If theta is in the major axis direction it returns the
major axis. On the other hand, if theta is in the direction
of the minor axis it returns the minor axis.
Parameters
----------
R : float
Radius along the major axis
q : list
Axis ratio of ellipse
pa : list
Angular position of ellipse measured from Y-axis (same as GALFIT)
theta : float
angle that defines the direction
Returns
-------
aell : float
radius along the angle specified by theta
"""
# changing measured angle from y-axis to x-axis
# and changing to rads:
newpa = (pa + 90) * np.pi / 180 # angle of every component
theta = (theta + 90) * np.pi / 180 # angle of direction of R
# bim = q * R
ecc = np.sqrt(1 - q**2)
alpha = theta - newpa # this is the direction
bell = R * np.sqrt(1 - (ecc * np.cos(alpha)) ** 2)
aell = bell / q # rad to evalue for every component
return aell
[docs]
def galfitLastFit(directory: str) -> str:
"""determine the last fit completed by GALFIT
Every time GALFIT produces a successful fit it produces
a file called galfit.XX, where the XX represent a number
that increases every time where GALFIT finds a solution
Parameters
----------
directory : str
directory containing the galfit files
Returns
-------
max_file : str
last fitting model file produced by GALFIT
"""
# Directory containing the files
# directory = "." #actual directory
# List all files in the directory
files = os.listdir(directory)
# Filter files that match the prefix 'galfit.'
galfit_files = [f for f in files if f.startswith("galfit.") and f[7:].isdigit()]
# Extract the numerical suffix and find the maximum
max_file = max(galfit_files, key=lambda x: int(x[7:]))
return max_file
[docs]
def galPrintComp(hdl: str, ncomp: int, idx: int, galcomps: GalComps) -> bool:
"""prints GALFIT component parameters to a file
Given a file handler, prints the GALFIT component parameters of
one model the surface brightness models to a GALFIT file
Parameters
----------
hdl : str
file handler where the component parameter information will be print
ncomp : int
prints the number of component in the GALFIT file
idx : int
index to indicate which component of the galcomps data class will be print
galcomps : GalComps data class defined above
Returns
-------
bool
"""
# print to filehandle
line00 = "# Object number: {} \n".format(ncomp)
line01 = " 0) {} # Object type \n".format(galcomps.NameComp[idx])
line02 = " 1) {:.2f} {:.2f} {} {} # position x, y [pixel] \n".format(
galcomps.PosX[idx],
galcomps.PosY[idx],
galcomps.PosXFree[idx],
galcomps.PosYFree[idx],
)
line03 = " 3) {:.4f} {} # total magnitude \n".format(
galcomps.Mag[idx], galcomps.MagFree[idx]
)
line04 = " 4) {:.2f} {} # R_e [Pixels] \n".format(
galcomps.Rad[idx], galcomps.RadFree[idx]
)
line05 = " 5) {:.2f} {} # Sersic exponent (deVauc=4, expdisk=1) \n".format(
galcomps.Exp[idx], galcomps.ExpFree[idx]
)
line06 = " 6) {} {} # ---------------- \n".format(
galcomps.Exp2[idx], galcomps.Exp2Free[idx]
)
line07 = " 7) {} {} # ---------------- \n".format(
galcomps.Exp3[idx], galcomps.Exp3Free[idx]
)
line08 = (
" 8) 0.0000 0 # ---------------- \n"
)
line09 = " 9) {:.3f} {} # axis ratio (b/a) \n".format(
galcomps.AxRat[idx], galcomps.AxRatFree[idx]
)
line10 = "10) {:.2f} {} # position angle (PA) \n".format(
galcomps.PosAng[idx], galcomps.PosAngFree[idx]
)
lineZ = " Z) {} # Skip this model in output image? \n".format(
galcomps.skip[idx]
)
line11 = "\n"
hdl.write(line00)
hdl.write(line01)
hdl.write(line02)
hdl.write(line03)
hdl.write(line04)
hdl.write(line05)
hdl.write(line06)
hdl.write(line07)
hdl.write(line08)
hdl.write(line09)
hdl.write(line10)
hdl.write(lineZ)
hdl.write(line11)
return True
[docs]
def galPrintSky(hdl: str, ncomp: int, galsky: GalSky) -> bool:
"""prints GALFIT sky component parameters to a file
Given a file handler, prints the GALFIT sky component parameters
to a GALFIT file
Parameters
----------
hdl : str
file handler where the sky component parameter information will be print
ncomp : int
prints the number of component in the GALFIT file
galsky : GalSky data class defined above
Returns
-------
bool
"""
line00 = "# Object number: {} \n".format(ncomp)
line01 = " 0) sky # Object type \n"
line02 = " 1) {:.3f} {} # sky background [ADU counts] \n".format(
galsky.sky, galsky.skyfree
)
line03 = " 2) {:.2f} {} # dsky/dx (sky gradient in x) \n".format(
galsky.dskyx, galsky.dskyxfree
)
line04 = " 3) {:.2f} {} # dsky/dy (sky gradient in y) \n".format(
galsky.dskyy, galsky.dskyyfree
)
line05 = " Z) {} # Skip this model in output image? (yes=1, no=0) \n".format(
galsky.skip
)
line06 = "\n"
line07a = "================================================"
line07b = "================================\n"
line07 = line07a + line07b
hdl.write(line00)
hdl.write(line01)
hdl.write(line02)
hdl.write(line03)
hdl.write(line04)
hdl.write(line05)
hdl.write(line06)
hdl.write(line07)
return True