physics77project/dat_generator.py at master · savsky95/physics77project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# this program is adopted from csvoutput.py and modified
# for binary file output instead of csv output.
# 12/4/2018

########### constants ##########
# constants are shared with n_body_simulation.py
# (do not change)
parsec      = 3.086 * 10**16 #[m]
Msun        = 1.989 * 10**30 #[kg]
G           = 6.67408 * 10**(-11) #[m^3 * kg^(-1) * s^(-2)]
pi          = 3.141592653
day         = 3600*24 #[sec]
year        = 3600*24*365 #[sec]
kpc         = 3.086 * 10**19 #[m]

########### parameters ##########
# Data file
filename = 'initial_data.dat'

# Galaxy property
gal_disk_r  = 25 * kpc #[m]
gal_disk_dz = 0.15 * kpc #[m]
gal_bulge_r = 10 * kpc #[m]
gass_mass   = 10**9 * Msun #[kg]
BH_mass     = 8.2 * 10**36 + gass_mass #[kg]
rho_0       = 4 * 10**7 * Msun / (kpc**3) #[kg/m^3] #from arXiv.1304.5127
r_c         = 60 * kpc #[m] #from arXiv:astro-ph/0403206

# Star property
star_v      = 150 * 10**3 #[m/s]
totalstar   = 1.5 * 10  ** 2
bulge_coef  = 0.4
bulgestar   = int(bulge_coef * totalstar)
num_stars   = int(totalstar - bulgestar)
actual_num  = 10 ** 11
mass_coef   = actual_num / num_stars

# simulation parameters
dt          = 5 * 10**5 * year #[sec]
t_max       = 5* 10**9 * year #[sec]
softening   = 15 * gal_disk_r / (num_stars)**(1/2) #[m]: mean distance
####################################


def initial_list_generator():
    import numpy as np
    # position generation
    num = np.arange(num_stars)
    mass = np.random.uniform(1*Msun*mass_coef,20*Msun*mass_coef,num_stars)
    x = np.random.triangular(-1*gal_disk_r,0,gal_disk_r,num_stars)
    y = []
    for i in range(len(x)):
        ylim = np.sqrt(gal_disk_r**2 - x[i]**2)
        y.append(*np.random.uniform(-1*ylim,ylim,1))
    z = np.random.uniform(-1*gal_disk_dz,gal_disk_dz,num_stars)

    # velocity generation
    vx = []
    vy = []
    vz = []
    for i in range(num_stars):
        d = np.sqrt(x[i]**2 + y[i]**2)
        v = np.random.uniform(star_v * 0.8, star_v * 1.2,1)
        vx.append(float(v * -1 * y[i] / d))
        vy.append(float(v * x[i] / d))
        vz.append(float(np.random.uniform(-10*10**3,10*10**3,1)))

    # starlist
    initial_list = []
    for i in range(num_stars):
        stardata = [num[i],mass[i],x[i],y[i],z[i],vx[i],vy[i],vz[i]]
        initial_list.append(stardata)
    return initial_list


def bulge ():
    import numpy as np
    # spherical coordinate
    num_bulge = int(bulgestar / 2)
    radius = np.random.uniform(0.0,gal_bulge_r,num_bulge)
    theta = np.random.uniform(0.,1.,num_bulge)*pi
    phi = np.arccos(1-2*np.random.uniform(0.0,1.,num_bulge))

    # number index
    num = np.arange(num_bulge * 2)

    # position
    x1 = radius * np.sin( theta ) * np.cos( phi )
    x  = np.concatenate((x1,-x1),axis=0)
    y1 = radius * np.sin( theta ) * np.sin( phi )
    y  = np.concatenate((y1,-y1),axis=0)
    z1 = radius * np.cos( theta )
    z  = np.concatenate((z1,z1),axis=0)

    # mass
    massbulge = np.random.uniform(1*Msun*mass_coef,20*Msun*mass_coef,bulgestar)

    # velocity
    vx = []
    vy = []
    vz = []
    for i in range(num_bulge*2):
        d = np.sqrt(x[i]**2 + y[i]**2)
        vel_net = np.sqrt(massbulge[i]*G/d)
        vx.append(float(vel_net * -1 * y[i] / d))
        vy.append(float(vel_net * x[i] / d))
        vz.append(float(np.random.uniform(-10*10**3,10*10**3,1)))

    # starlist_bulge
    starlist_bulge = []
    for i in range(num_bulge * 2):
        each_star = [num[i],massbulge[i],x[i],y[i],z[i],vx[i],vy[i],vz[i]]
        starlist_bulge.append(each_star)
    return starlist_bulge


def condition_data_generator():
    names  = ['disk_r','disk_dz','bulge_r','BH_m','DM_rho0','r_c','star_v',\
            'num_stars','bulgestars','actual_num','mass_coef','dt','t_max','softening']
    values = [gal_disk_r,gal_disk_dz,gal_bulge_r,BH_mass,rho_0,r_c,star_v,\
            num_stars,int(0.3*totalstar),actual_num,mass_coef,dt,t_max,softening]
    condition_data = [names,values]
    print("Initial Condition Parameters:")
    for i in range(len(names)):
        print("\t* {:10}\t= {:1.2}".format(names[i],float(values[i])))

    return condition_data

def output(filename,condition_data,initial_list):
    import pickle
    output = open(filename,'wb')
    outputdata = [condition_data,initial_list]
    pickle.dump(outputdata,output)
    output.close()
    return True


### main ###
initial_list   = initial_list_generator()
condition_data = condition_data_generator()
bulge_list     = bulge()
initial_list = initial_list + bulge_list
if output(filename,condition_data,initial_list):
    print('Data was successfully created and saved to \'{}\''.format(filename))