-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgraph_densite.m
126 lines (91 loc) · 2.57 KB
/
graph_densite.m
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
% clear *
global dsup dinf Ro elev Rcoro
% Rcoro = 3500;
% elev = 26;
% Ro = 8000;
% dsup = 14000.;
% dinf = 100.; %distance en parsec
global l b
% definition de la fenetre de Baade dans la majeure partie des articles : l = 1 et b = -4
% A priori c'est cette definition qui est juste.
l = 1 *pi/180; % direction d'observation en radian
b = -4 *pi/180;
global sinb cosb cosbl sinl cosl
sinb = abs(sin(b)); cosb = cos(b); cosl = cos(l);
cosbl=cos(b)*cos(l); sinl = sin(l);
x = (0:1e-5:1).*(dsup-dinf)+dinf;
[R, z, th] = toGC(x);
% dens_zhao = rhozhao(R, z, th);
% dens_G2 = rhodwek(R, z, th);
% dens_E2 = rhostanek(R, z, th);
% dens_HetG = rhobuHetG(R, z, th);
dens_d = rhobulbe(R, z, th);
dens_dm = rhodm(R, z, th);
dens_de = rhode(R, z, th);
i0 = find( R <= Rcoro );
i1 = find( R > Rcoro);
% vrot(i1) = vrotdm(R(i1),z(i1),th(i1));
% vrot(i0) = vrotb(R(i0),z(i0),th(i0));
figure(42)
plot(x, (dens_d+dens_dm+dens_de)*4/max(dens_d+dens_dm+dens_de))
figure()
hold on;
% plot(x,vrotdm(x,z,th).*1e-3);
% title('Vitesse de rotation en fonction de R');
% xlabel('distance au centre galactique')
% ylabel('vitesse en km/s')
% plot(x, vrotb(R,z,th).*1e-2);
% plot(x, dens_zhao);
% plot(x, dens_E2);
% plot(x, dens_G2);
% legend('bulbe (Zhao)', 'E2 (Stanek)', 'G2 (dwek)', 'H&G');
plot(x, dens_d)
plot(x, dens_dm)
plot(x, dens_de)
plot(x, (dens_d+dens_dm+dens_de))
legend('bulbe', 'disque mince', 'dsque épais', 'densité totale')
set(gca, 'YScale', 'log')
xlabel('distance au soleil (en pc))');
ylabel('densité de masse en M_{sol}/pc^{3}');
%% Integration pour le calcul de la masse
clear *
maxx = 1e5;
format long
m_bu = integral3(@rho_stanek_simple,0, maxx, 0, maxx, 0, maxx);
disp(['masse bulbe ', num2str(m_bu*8*1e-10)])
function rh = rho_stanek_simple(x,y,z)
x0=890;
y0=x0.*4.3/10;
z0=x0.*2.8/10;
M_b = 2e10;
rho0 = M_b/(x0*y0*z0*8*pi);
a = sqrt((x./x0).^2+(y./y0).^2+(z./z0).^2);
rh = rho0.*exp(-a);
end
function rh = rho_dweck_simple(X,Y,Z)
x0=890;
y0=x0.*4.3/10;
z0=x0.*2.8/10;
M_b = 1e10;
rho0 = M_b/(6.57*pi*x0*y0*z0);
% rho0 = 1/(2*pi*2^(3/2)*2*0.675*x0*y0*z0);
sb2=sqrt(((X/x0).^2+(Y/y0).^2).^2+(Z/z0).^4);
%-----------------
% densite de masse
%-----------------
rh=rho0*(exp(-sb2/2));
end
function rh = rho_zhao_simple(X,Y,Z)
x0=1490;
y0=580;
z0=400;
qa=0.6;
M_b = 2;
rho0 = M_b/0.9221;
sb2=sqrt(((X/x0).^2+(Y/y0).^2).^2+(Z/z0).^4);
sa=sqrt((qa*qa*(X.^2+Y.^2)+Z.^2)/(z0^2));
%------------- ----
% densite de masse
%-----------------
rh=rho0*(exp(-sb2/2)+sa.^(-1.85).*exp(-sa));
end