Commit 2c49b0d7 authored by bepinat's avatar bepinat


parent 3fa2658b
......@@ -219,198 +219,3 @@ class Model2D:
return np.sum(chi2)
class Model3D: # XXX do the same thing: initiate, then produce velocity map, intensity/flux map and velocity dispersion map at high resolution and then produce a cube with the same sampling as data and adjusted with respect to a reference line and wavelengths in the cube. Then convolve with PSF and rebin. Everything needed is in the program with made with Debora
Models in 3D of the datacube around a line of galaxies
Also compute the velocity field and the dispersion of a model using observed cube in input.
This model can be used with any rotational curve models (in '').
def __init__(self, cube, vcube, flux, psf, vel_model, psfz=0, sig0=0, slope=0., header=None, z=0., lbda0=6562.78, cont=0.):
:param Image cube: represents the cube to fit
:param Image vcube: the variance cube
:param Image flux: flux distribution at the same or in high resolution than the velocity
:param PSF psf: the class of the psf
:param func vel_model: velocity function used to create the model
:param float psfz: velocity dispersion corresponding to the spectral resolution
:param float sig0: velocity dispersion of the model
:param float slope: slope of the velocity dispersion
:param hdufits header header: header associated to the datacube
:param float z: systemic redshift of the source
:param float lbda0: restframe wavelength of the line used (could be a list of lines and then need to add line ratios, etc.)
:param float cont: continuum level (assumed constant for now, but could be polynomial)
# Parameters which must be initialized
self.vel_model = vel_model.vel_model
self.model_sig0 = sig0
self.model_slope = slope
self.psfz = psfz
# XXX a reprendre en 3D
self.cube_mod = np.zeros(cube.shape)
#self.vel_map = np.zeros(cube.shape[1:])
self.vel_map_hd = np.zeros(np.array(flux.get_size()))
self.disp_map_hd = np.zeros(np.array(flux.get_size()))
#self.cube_hd = np.zeros((cube.shape[0], np.array(flux.get_size()))) # Check def works (dimension declaration)
self.cube = cube
self.vcube = vcube
self.header = header
self.flux = flux
self.psf = psf
self.oversamp = self.flux.get_oversamp()
self.indices = tools.lowres_indices_highres_frame(flux.get_size(), self.oversamp)
self.z = z
self.lbda0 = lbda0
self.lbda1 = lbda0 * (1 + z)
self.cont = cont # XXX devrait être une image en fait!
# We initiate the wavelength range of the datacube
self.lbda = (np.arange(header['NAXIS3']) - header['CRPIX3']) * header['CDELT3'] + header['CRVAL3'] # XXX Check units!
self.vind = (self.lbda - self.lbda1) / self.lbda1 * ct.c * 1e-3 # km/s ; XXX Ne fonctionne que s'il n'y a qu'une seule raie
# We initiate the geometrical model parameters
self.model_pa = None
self.model_inc = None
self.model_xc = None
self.model_yc = None
self.model_vs = None
# We initiate the associated coordinate maps in the galaxy frame
self.model_radius = None
self.model_theta = None
# We write an ordered list containing the names of the model parameters
self.model_parname = ['xc', 'yc', 'pa', 'inc', 'vs']
self.model_modname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
self.model_compname = ['geometrical', 'geometrical', 'geometrical', 'geometrical', 'geometrical']
# We initiate a dictionary containing the names of rotation curve parameters. The number of arguments adjusts automatically
#ii = 0
for name in vel_model.model_parname:
for name in vel_model.model_modname:
for name in vel_model.model_compname:
#ii += 1
#for name in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount]:
#self.model_RCparams_names = self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount] # can be useful to check the parameters are provided with the correct order
#self.model_RCparams_values = [None for par in self.model_RCparams_names]
#self.model_RCparams = dict((par, None) for par in self.vel_model.__code__.co_varnames[1:self.vel_model.__code__.co_argcount])
self.model_RCparams = None
def set_parameters(self, xc, yc, pa, inc, vs, *args):
Sets the value of the model parameters
:param float xc: abscissa of the center in pixel
:param float yc: ordinate of the center in pixel
:param float pa: position angle of the major axis in degree
:param float inc: inclination of the disk in degree
:param float vs: systemic velocity in km/s
:param list of float args: list that contains other model arguments
self.model_pa = pa
self.model_inc = inc
self.model_xc = xc
self.model_yc = yc
self.model_vs = vs
#self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.model_xc, self.model_yc, self.model_pa, self.model_inc, im_size=np.shape(self.vel_map_hd), oversamp=self.oversamp)
self.model_radius, self.model_theta = tools.sky_coord_to_galactic(self.indices[1], self.indices[0], self.model_xc, self.model_yc, self.model_pa, self.model_inc)
self.model_RCparams = args
def get_parameter(self):
Gets the actual parameters of the model (in low resolution scale)
:return ndarray:
return [self.model_xc, self.model_yc, self.model_pa, self.model_inc, self.model_vs, *self.model_RCparams]
def disk_velocity(self):
Computes the velocity field
:return ndarray:
vr = self.vel_model(self.model_radius, *self.model_RCparams)
# Calculation of the velocity field
v = vr * math.sin(math.radians(self.model_inc)) * self.model_theta + self.model_vs
return v
def linear_velocity_dispersion(self):
Returns the velocity dispersion map needed for the fit process
:return ndarray: velocity dispersion map
# Calculation of the velocity dispersion
sig = self.model_sig0 + self.model_slope * np.abs(self.model_radius)
sig[np.where(sig <= 0)] = 0
return sig
def cube_model(self):
Computes a high resolution cube
: return cube:
self.vel_map_hd = self.disk_velocity()
self.disp_map_hd = self.linear_velocity_dispersion() # besoin du centre du modèle, donc ici! En principe pour la mmodélisation du cube, il faudrait que sigma et slope soient des paramètres libre du modèle
# Il faudrait en principe donner la LSF en input (de la même manière que la PSF). Dans un premier temps, on peut faire l'hypothèse que la LSF est Gaussienne, c'est donc inclus dans sigma
self.cube_hd = self.flux * np.exp(-np.subtract(self.vind.reshape(len(vind), 1, 1), self.vel_map_hd.reshape(1, self.flux.shape[0], self.flux.shape[1])) ** 2 / (2 * self.disp_map_hd ** 2)) + self.cont
self.cube_mod = tools.rebin_data(self.psf.convolution(self.cube_hd), self.oversamp)
def least_square(self, p, fjac=None):
Function minimized by mpfit.
Returns (data-model)^2/err^2 XXX Why no square???
:param ndarray p: array of parameters
:return ndarray:
return [0, np.reshape((self.cube - self.cube_mod)/np.sqrt(self.vcube), -1)] # sqrt, to have uncertainty, i.e. stdev
def log_likelihood(self, cube, ndim, nparams):
Log likelihood function which maximized by multinest
Returns -sum[(data-model)^2/(2*err^2)]
:param ndarray cube: data whith n_params dimension
:param int ndim: number of dimension if different of the number of paramerters
:param int nparams: number of parameters
:return float:
#self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], cube[5], cube[6])
self.set_parameters(cube[0], cube[1], cube[2], cube[3], cube[4], *cube[5:ndim])
chi2 = -(self.cube_mod - self.cube)**2 / (2*self.vcube) # no **2 because we use variance
return np.sum(chi2)
#class Model1D:
This diff is collapsed.
......@@ -10,8 +10,7 @@ from import fits
from Models.velocity_model import combined_velocity_model
import as tools
from Class.Images import Image
from Class.Models import Model2D
#from Class.Model2D import Model3D
from Class.Model2D import Model2D
from Class.PSF import PSF
from SubProcess.use_mpfit import use_mpfit
from SubProcess.use_multinest import use_multinest
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment