<jupyter_start><jupyter_text>Here we illustrate how to sample n-body particles with Galaxia.A simple application is to provide a flat starformation history with reasonable age-metallicity relation in order to sample a colour magnitude diagram and look at the occupation of CMD space by the different stellar populations and look at their astrophysical parameters. Like that we can for example sample the age distribution of RGB stars. The technique was first described here: http://adsabs.harvard.edu/abs/2016AN....337..880J<jupyter_code>%pylab inline
import ebf
import shutil
import subprocess
from astropy.io import fits
import os
# Need to specify where the GalaxiaData folder is and how to name the new simulation

nbody_folder = '/home/rybizki/Programme/GalaxiaData/'
folder = 'flat_sfr/'
filename = 'flat_sfr_selection_down_to_5'
folder_cat = '../output/mock_cat'
def create_flat_sfr(nbody_folder,folder,filename):
    
    # 2 input files for Galaxia need to be created
    folder_create = nbody_folder + 'nbody1/' + folder
    if os.path.exists(folder_create):
        shutil.rmtree(folder_create)
        os.mkdir(folder_create)
        print(folder_create, "existed and was recreated")
    else:
        os.mkdir(folder_create)

    # Age metallicity distribution is specified in this input file, using Model A from Just&Jahreiß 2010
    s = '../input/zevolv02.dat'
    data1 = np.genfromtxt(s,skip_header=2, dtype=None, names = True)

    age = data1['t']
    age = age*13.5
    zfe = data1['zfe']
    zox = data1['zox']
    #fehdex = np.log10(zfe)
    # Here we can specify the AMR according to our needs
    age = np.linspace(0,13.5,num = 1001)[1:]
    fehdex = np.linspace(-1.8,-0.7,num = len(age))
    #ohdex = np.log10(zox)
    ohdex = np.ones(len(age))
    enhancement = ohdex - fehdex
    age = age[::-1]
    mass = np.divide(np.ones(len(age)),len(age))
    #################################

    # Here the ebf nbody input file is created
    pos = np.zeros((len(age),3))
    vel = np.zeros((len(age),3))
    for i in range(0,len(pos)):
        fehdex[i] = fehdex[i]
        pos[i] = np.array([ -8.00000000e+00,   0.00000000e+00,   1.50000000e-02])
        vel[i] = np.array([ 1.11000000e+01,   2.39080000e+02,   7.25000000e+00])
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/mass', mass,'w')	
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/feh', fehdex,'a')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/id', 1,'a')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/alpha', enhancement,'a')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/age', age,'a')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/pos3', pos,'a')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '.ebf', '/vel3', vel,'a')

    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '_d3n64_den.ebf', '/h_cubic', np.zeros(len(pos)), 'w')
    ebf.write(nbody_folder + 'nbody1/' + folder + filename + '_d6n64_den.ebf', '/h_cubic', np.zeros((len(pos),2)), 'w')

    # Here the file which tells Galaxia where to find the input file is created
    filedata = 'nbody1/%s\n         1       1\n%s.ebf\n' %(folder, filename)
    file = open(nbody_folder + "nbody1/filenames/" + filename + ".txt", "w")
    file.write(filedata)
    file.close()
    return(fehdex,age)
fehdex,age = create_flat_sfr(nbody_folder,folder,filename)
# Changing the AMR to something like fig.4 of https://ui.adsabs.harvard.edu/#abs/2014ApJ...786...44S/abstract
# Just checking if it looks OK. Change the AMR if needed within the function create_flat_sfr
plt.plot(age,fehdex)
plt.xlabel("age")
plt.ylabel("[Fe/H]")
plt.title("AMR")<jupyter_output>('/home/rybizki/Programme/GalaxiaData/nbody1/flat_sfr/', 'existed and was recreated')<jupyter_text>Preparing a CMD of constant SFR to see what the age distribution of a specific population is<jupyter_code># Here we create the folder where Galaxia will write the mock catalogue

folder_create = folder_cat + '/'
if os.path.exists(folder_create):
    shutil.rmtree(folder_create)
    os.mkdir(folder_create)
    print(folder_create, "existed and was recreated")
else:
    os.mkdir(folder_create)

# Specifying the mag limits of that simulation
magmin = -100
magmax = 5
colmin = -100
colmax = 100
# Oversampling and seed, 
# Try to fix that value to something reasonable such that the age of the stellar population is sampled well
iterations = 1e6
seed = 1

# Creating the parameterfile
filedata = 'outputFile                          %s\noutputDir                           %s\nphotoSys                            UBV\nmagcolorNames                       V,B-V\nappMagLimits[0]                     -100\nappMagLimits[1]                     1000\nabsMagLimits[0]                     %f\nabsMagLimits[1]                     %f\ncolorLimits[0]                      %f\ncolorLimits[1]                      %f\ngeometryOption                      0\nlongitude                           0\nlatitude                            90\nsurveyArea                          1000\nfSample                             %d\npopID                               -1\nwarpFlareOn                         1\nseed                                %d\nr_max                               1000\nstarType                            0\nphotoError                          0\n' %(filename, folder_cat,magmin,magmax,colmin,colmax,iterations,seed)
myparameterfile = folder_cat + '/myparameterfile'
file = open(myparameterfile, "w")
file.write(filedata)
file.close()

# Starting Galaxia
p = subprocess.Popen(['galaxia', '-r', '--nfile=%s' %(filename), '--hdim=3', myparameterfile], stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
print out
print err
# Adding other photometric bands
args = ['galaxia', '-a', '--psys=Gaia' ,folder_cat + '/' + filename + '.ebf']
p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)	
print("Galaxia adds Gaia bands")
(output, err) = p.communicate()
#p_status = p.wait()
print(output)
# Adding other photometric bands
args = ['galaxia', '-a', '--psys=SDSS' ,folder_cat + '/' + filename + '.ebf']
p = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)	
print("Galaxia adds SDSS bands")
(output, err) = p.communicate()
#p_status = p.wait()
print(output)
# Adding the library
path = os.path.abspath('../library/')
if path not in sys.path:
    sys.path.append(path)
from convert_to_recarray import convert

# Reading in the Catalogue and converting it to npy file. (If you added photometric bands then you will have to edit the convert routine)
bes = ebf.read(folder_cat + '/' + filename + '.ebf','/')
x = convert(bes)
x = x.view(np.recarray)
x.age = np.divide(np.power(10,x.age),1e9)
x.teff = np.power(10,x.teff)
print(len(x))
print(x.dtype.names)
# Writing into a fits file so that we can inspect the result in topcat
fits_create = folder_cat + '/' + filename + '.fits'
if os.path.isfile(fits_create):
    os.remove(fits_create)
    print(fits_create, "existed and was deleted and respawned")
else:
    pass
fits.writeto(fits_create,x)
os.remove(folder_cat + '/' + filename + '.ebf')
#Load mock catalogue
x = fits.getdata(fits_create)
print(len(x),x.dtype)
x = x.view(np.recarray)
plt.plot(x.ubv_b-x.ubv_v,x.ubv_v,".")
plt.ylim(6,-7)

rgb = x[(x.ubv_v < -0.5) & (x.ubv_b - x.ubv_v > 0.8)]
plt.plot(rgb.ubv_b-rgb.ubv_v,rgb.ubv_v,".")
plt.ylim(-0.4,-7)
n,bins,patch =plt.hist(rgb.age, bins = 30, normed=True)
plt.xlabel("age")
plt.ylabel('sampling of red clump star population')
new_bins= np.divide(bins[:-1] + bins [1:],2)
rgb_time = new_bins[::-1]
rgb_dist = n[::-1]
plt.plot(rgb_time,rgb_dist)
plt.xlabel("age")
plt.ylabel('sampling of RGB star population')
np.save('../output/selection/rgb_time.npy',rgb_time)
np.save('../output/selection/rgb.npy',rgb_dist)
# This format now looks similar to the Chempy input, 
# except that the time axis does not exactly start and end at 0 and 13.5Gyr
# Might, but should not have an impact
rc = np.load('../input/chempy_selection/red_clump_new.npy')
rc_time = np.load('../input/chempy_selection/time_red_clump_new.npy')
plt.plot(rc_time,rc)
plt.xlabel("age")
plt.ylabel('sampling of red clump star population')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14053272267084488, 'fraction_numerical': 0.02615403157767868, 'mean_word_length': 3.1740442655935612, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '254388', 'license': 'MIT', 'n_tokens_mistral': 3064, 'n_tokens_neox': 2877, 'n_words': 710, 'path': 'notebook/notebook_aida/[2]Providing n-body particle input for Galaxia.ipynb', 'repo_name': 'jan-rybizki/Galaxia_wrap', 'size': '8297'}