Homec4science
Diffusion ARRAKIHS (master) Browse

arrakihs/Workspace/InjectedDwarfs_conflimit/ark_idealized_dwarfs_generate
d192e801c42bmaster

Last Change
Raw File

ark_idealized_dwarfs_generate
View Options

#!/usr/bin/env python3
import numpy as np
from pNbody import *
from pNbody.mass_models import plummer
from pNbody import ic
from astropy import constants as cte
from astropy import units as u
import argparse
from scipy.stats import skewnorm
from scipy.stats import rv_continuous
from scipy.stats.sampling import NumericalInversePolynomial
from tqdm import tqdm
from pNbody import pychem
from pNbody import Isochrones
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import copy
import pickle
from astropy.io import fits
from pNbody import Mockimgs
from pNbody.Mockimgs.obs_object import Object
from pNbody.Mockimgs import utils
from pNbody.Mockimgs.los import LineOfSights
from pNbody.Mockimgs import filters
from pNbody.Mockimgs import luminosities
####################################################################
# option parser
####################################################################
description="""generate a dwarf galaxy
"""
epilog ="""
Examples:
--------
ark_idealized_dwarfs_generate --Fe -2.5 --age 3 --M0 1e4 -e 0.1 -o 1e4.hdf5
"""
parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument("-o",
action="store",
type=str,
dest="outputfilename",
default=None,
help="Name of the output file")
parser.add_argument("--AgeMean",
action="store",
dest="AgeMean",
metavar='FLOAT',
type=float,
default=1,
help='mean Age')
parser.add_argument("--AgeFlat",
action="store_true",
default=False,
help="flatten age")
parser.add_argument("--Fe",
action="store",
dest="Fe",
metavar='FLOAT',
type=float,
default=-2.5,
help='Fe mode')
parser.add_argument("--M0",
action="store",
dest="M0",
metavar='FLOAT',
type=float,
default=1e5,
help='stellar mass to generate [Msol]')
parser.add_argument("-e",
action="store",
dest="e",
metavar='FLOAT',
type=float,
default=0.5,
help='scale length [kpc]')
parser.add_argument("--rmax_e_ratio",
action="store",
dest="rmax_e_ratio",
metavar='FLOAT',
type=float,
default=10,
help='rmax to scale length ratio')
parser.add_argument("--ftype",
action="store",
type=str,
dest="ftype",
default="swift",
help="type of file")
parser.add_argument("--do_not_compute_rsp",
action="store_true",
default=False,
help="do not compute rsp")
parser.add_argument("--nngb",
action="store",
type=int,
dest="nngb",
default=5,
help="Number of neighbouring particles to consider to compute RSP(==HSML)")
parser.add_argument("-N",
action="store",
dest="N",
metavar='INT',
type=int,
default=10000,
help='number of particles')
class AgesDF():
def __init__(self,AgeMean=2,AgeMax=13.7,flat=False):
self.AgeMax = AgeMax
self.AgeMean = AgeMean
self.s = AgeMean/AgeMax/2
self.flat = flat
self.x = np.linspace(0,1,10000)
self.y = self.DF(self.x)
self.yc = np.add.accumulate(self.y)
self.yc = self.yc/max(self.yc)
def DF(self,x):
'''
max pour x=s
max = 1
'''
n=1
y = x**n * np.exp(-x**n/self.s**n) #/ (self.s**2*np.exp(-1))
if self.flat:
xpeak = self.x[np.argmax(y)]
ypeak = np.max(y)
y = np.where(self.x>=xpeak,ypeak,y)
return y
def Sample(self,n):
X = np.random.random(n)
x = np.interp(X,self.yc,self.x)
return x*self.AgeMax
####################################################################
# main
####################################################################
def Do(opt):
n = opt.N
#############################################
# set particle mass
#############################################
mass = np.ones(n)*opt.M0/n
mass = mass/1e10 # output units are in 1e10 Msol
#############################################
# set ages
#############################################
df = AgesDF(AgeMean=opt.AgeMean,flat=opt.AgeFlat)
Ages = df.AgeMax - df.Sample(n)
Ages.sort()
print('Ages done')
#############################################
# generate metallicities
#############################################
Fes = opt.Fe + np.log10(-np.log(1-np.random.rand(n)))
Fes.sort()
Fes = np.flip(Fes)
#############################################
# generate a plummer model
#############################################
# set maximal radius
rmax = opt.e * opt.rmax_e_ratio
# generate the model
nb = ic.plummer(n,1,1,1,opt.e,rmax,M=1,irand=1,vel='no',name=opt.outputfilename,ftype=opt.ftype)
# compute Hsml
if opt.do_not_compute_rsp is False:
print("Compute Rsp...")
# remove doublets
#u,idx = np.unique(nb.rxyz(),return_index=True)
#nb = nb.selectp(lst=nb.num[idx])
nb.ComputeRsp(opt.nngb)
print("done.")
# set particles to be stars
md = nb.getParticleMatchingDict()
nb.set_tpe(md["stars"])
# set masses (as mass or luminosities)
nb.mass = mass
# store age
if Ages is not None:
nb.age = Ages # store in Gyr
# store metallicity
if Fes is not None:
nb.mh = Fes
# define units
u_Length = 1* u.kpc
u_Mass = 10**10 * u.M_sun
u_Velocity = 1* u.km/u.s
u_Time = u_Length/u_Velocity
toMsol = u_Mass.to(u.M_sun).value
# add units
nb.UnitLength_in_cm = u_Length.to(u.cm).value
nb.UnitMass_in_g = u_Mass.to(u.g).value
nb.UnitVelocity_in_cm_per_s = u_Velocity.to(u.cm/u.s).value
nb.Unit_time_in_cgs = u_Time.to(u.s).value
# non cosmological
nb.setComovingIntegrationOff()
nb.cosmorun=0
nb.time =0
if opt.outputfilename is not None:
nb.rename(opt.outputfilename)
nb.write()
if __name__ == '__main__':
opt = parser.parse_args()
Do(opt)
c4science · Help