Line data Source code
1 : #include "crpropa/massDistribution/Ferriere.h"
2 : #include "crpropa/Common.h"
3 :
4 : #include "kiss/logger.h"
5 :
6 : #include <sstream>
7 :
8 : namespace crpropa {
9 :
10 13 : Vector3d Ferriere::CMZTransformation(const Vector3d &position) const {
11 : // set galactocentric coordinate system with the Sun at (-8.5,0.,0.) instead of (8.5, 0, 0) to be consistand with JF12 implementation
12 13 : double x = -position.x;
13 13 : double y = -position.y;
14 :
15 : double xC = -50*pc; //offset
16 : double yC = 50*pc;
17 : double sinTc = sin(70.*deg);
18 : double cosTc = cos(70.*deg);
19 :
20 : Vector3d pos;
21 13 : pos.x = (x - xC)*cosTc + (y - yC)*sinTc;
22 13 : pos.y = -(x - xC)*sinTc + (y - yC)*cosTc;
23 13 : pos.z = position.z;
24 :
25 13 : return pos;
26 : }
27 :
28 13 : Vector3d Ferriere::DiskTransformation(const Vector3d &position) const {
29 : // set galactocentric coordinate system with the Sun at (-8.5,0.,0.) instead of (8.5, 0, 0) to be consistand with JF12 implementation
30 13 : double x = -position.x;
31 13 : double y = - position.y;
32 13 : double z = position.z;
33 :
34 : double alphaD = 13.5*deg; // rotation arround x-axis
35 : double sinAd = sin(alphaD);
36 : double cosAd = cos(alphaD);
37 : double betaD = 20.*deg; // rotation arround y'-axis
38 : double sinBd = sin(betaD);
39 : double cosBd = cos(betaD);
40 : double TettaD = 48.5*deg; // rotation arround x"-axis
41 : double sinTd = sin(TettaD);
42 : double cosTd = cos(TettaD);
43 :
44 : Vector3d pos;
45 :
46 13 : pos.x = x*cosBd*cosTd - y*(sinAd*sinBd*cosTd -cosAd*sinTd)-z*(cosAd*sinBd*cosTd +sinAd*sinTd);
47 :
48 13 : pos.y = -x*cosBd*sinTd;
49 13 : pos.y += y*(sinAd*sinBd*sinTd +cosAd*cosTd);
50 13 : pos.y += z*(cosAd*sinBd*sinTd -sinAd*cosTd);
51 :
52 13 : pos.z = x*sinBd;
53 13 : pos.z += y*sinAd*cosBd;
54 13 : pos.z += z*cosAd*cosBd;
55 :
56 13 : return pos;
57 : }
58 :
59 12 : double Ferriere::getHIDensity(const Vector3d &position) const {
60 : double n = 0;
61 12 : double R = sqrt(position.x*position.x+position.y*position.y);
62 :
63 12 : if(R<3*kpc)
64 : {
65 : // density at center
66 6 : Vector3d pos = CMZTransformation(position); // coordinate trafo
67 6 : double x = pos.x/pc; // all units in pc
68 6 : double y = pos.y/pc;
69 6 : double z = pos.z/pc;
70 :
71 6 : double A = sqrt(x*x+2.5*2.5*y*y);
72 6 : double nCMZ = 8.8/ccm*exp(-pow_integer<4>((A-125.)/137))*exp(-pow_integer<2>(z/54.));
73 :
74 : // density in disk
75 6 : pos = DiskTransformation(position); // Corrdinate Trafo
76 6 : x = pos.x/pc; // all units in pc
77 6 : y = pos.y/pc;
78 6 : z = pos.z/pc;
79 :
80 6 : A = sqrt(x*x+3.1*3.1*y*y);
81 6 : double nDisk = 0.34/ccm*exp(-pow_integer<4>((A-1200.)/438.))*exp(-pow_integer<2>(z/120));
82 :
83 6 : n = nCMZ + nDisk;
84 : }
85 : else{ // outer region
86 6 : double z = position.z/pc;
87 : double a;
88 6 : if(R<=Rsun){
89 : a = 1;
90 : } else {
91 3 : a = R/Rsun;
92 : }
93 :
94 6 : double nCold = 0.859*exp(-pow_integer<2>(z/(127*a))); // cold HI
95 6 : nCold += 0.047*exp(-pow_integer<2>(z/(318*a)));
96 6 : nCold += 0.094*exp(-fabs(z)/(403*a));
97 6 : nCold *= 0.340/ccm/(a*a);
98 :
99 6 : double nWarm = (1.745 - 1.289/a)*exp(-pow_integer<2>(z/(127*a))); // warm HI
100 6 : nWarm += (0.473 - 0.070/a)*exp(-pow_integer<2>(z/(318*a)));
101 6 : nWarm += (0.283 - 0.142/a)*exp(-fabs(z)/(403*a));
102 6 : nWarm *= 0.226/ccm/a;
103 :
104 6 : n = nWarm + nCold;
105 : }
106 :
107 12 : return n;
108 : }
109 :
110 12 : double Ferriere::getHIIDensity(const Vector3d &position) const {
111 : double n = 0;
112 12 : double R = sqrt(position.x*position.x+position.y*position.y);
113 :
114 12 : if(R< 3*kpc){ // inner
115 6 : double x = position.x/pc;
116 6 : double y = position.y/pc;
117 6 : double z = position.z/pc;
118 :
119 : // warm interstellar matter
120 6 : double H = (x*x + pow_integer<2>(y+10.))/(145*145);
121 6 : double nWIM = exp(-H)* exp(-pow_integer<2>(z+20.)/(26*26));
122 6 : nWIM += 0.009*exp(-pow_integer<2>((R/pc-3700)/(0.5*3700)))*1/pow_integer<2>(cosh(z/140.));
123 6 : nWIM += 0.005*cos(M_PI*R/pc*0.5/17000)*1/pow_integer<2>(cosh(z/950.));
124 6 : nWIM *= 8.0/ccm; // rescaling with density at center
125 :
126 : //very hot interstellar matter
127 : double alphaVH = 21.*deg; // angle for very hot IM in radiant
128 : double cosA = cos(alphaVH);
129 : double sinA = sin(alphaVH);
130 6 : double etta = y*cosA+z*sinA; // coordinate transformation for VHIM along major axis
131 6 : double chi = -y*sinA+z*cosA;
132 :
133 6 : double nVHIM = 0.29/ccm*exp(-((x*x+etta*etta)/(162.*162.)+chi*chi/(90*90)));
134 :
135 6 : n = nWIM + nVHIM;
136 : } else { // outer region
137 6 : double z = position.z;
138 :
139 6 : double nWarm = 0.0237/ccm*exp(-(R*R-Rsun*Rsun)/pow_integer<2>(37*kpc))*exp(-fabs(z)/(1*kpc));
140 6 : nWarm += 0.0013/ccm*exp(-(pow_integer<2>(R-4*kpc)-pow_integer<2>(Rsun-4*kpc))/pow_integer<2>(2*kpc))*exp(-fabs(z)/(150*pc));
141 :
142 6 : double nHot = 0.12*exp(-(R-Rsun)/(4.9*kpc));
143 6 : nHot += 0.88*exp(-(pow_integer<2>(R-4.5*kpc)-pow_integer<2>(Rsun-4.5*kpc))/pow_integer<2>(2.9*kpc));
144 6 : nHot *= pow(R/Rsun, -1.65);
145 6 : nHot *= exp(-fabs(z)/((1500*pc)*pow(R/Rsun,1.65)));
146 6 : nHot *= 4.8e-4/ccm;
147 :
148 6 : n= nWarm + nHot;
149 : }
150 12 : return n;
151 : }
152 :
153 12 : double Ferriere::getH2Density(const Vector3d &position) const{
154 : double n=0;
155 12 : double R=sqrt(position.x*position.x+position.y*position.y);
156 :
157 12 : if(R<3*kpc) {
158 : // density at center
159 6 : Vector3d pos =CMZTransformation(position); // coord trafo for CMZ
160 6 : double x = pos.x/pc; // all units in pc
161 6 : double y = pos.y/pc;
162 6 : double z = pos.z/pc;
163 :
164 6 : double A = sqrt(x*x+pow(2.5*y,2)); // ellipticity
165 6 : double nCMZ = exp(-pow((A-125.)/137.,4))*exp(-pow(z/18.,2));
166 6 : nCMZ *= 150/ccm; // rescaling
167 :
168 : // density in disk
169 6 : pos = DiskTransformation(position); // coord trafo for disk
170 6 : x=pos.x/pc;
171 6 : y=pos.y/pc;
172 6 : z=pos.z/pc;
173 :
174 6 : A = sqrt(x*x+pow_integer<2>(3.1*y));
175 6 : double nDISK = exp(-pow_integer<4>((A-1200)/438))*exp(-pow_integer<2>(z/42));
176 6 : nDISK *= 4.8/ccm; // rescaling
177 :
178 6 : n = nCMZ + nDISK;
179 : } else { // outer region
180 6 : double z = position.z/pc;
181 6 : n = pow(R/Rsun, -0.58);
182 6 : n *= exp(-(pow_integer<2>(R-4.5*kpc)-pow_integer<2>(Rsun-4.5*kpc))/pow_integer<2>(2.9*kpc));
183 6 : n *= exp(-pow_integer<2>(z/(81*pow(R/Rsun,0.58))));
184 6 : n *= 0.58/ccm; // rescaling
185 : }
186 :
187 12 : return n;
188 : }
189 :
190 5 : double Ferriere::getDensity(const Vector3d &position) const{
191 : double n=0;
192 5 : if(isforHI){
193 4 : n += getHIDensity(position);
194 : }
195 5 : if(isforHII){
196 4 : n+=getHIIDensity(position);
197 : }
198 5 : if(isforH2){
199 4 : n+=getH2Density(position);
200 : }
201 : // check if all densities are deactivated and raise warning if so
202 5 : if((isforHI || isforHII || isforH2) == false){
203 2 : KISS_LOG_WARNING
204 1 : << "\nCalled getDensity on fully deactivated Ferriere \n"
205 1 : << "gas density model. The total density is set to 0.";
206 : }
207 :
208 5 : return n;
209 : }
210 :
211 5 : double Ferriere::getNucleonDensity(const Vector3d &position) const{
212 : double n=0;
213 5 : if(isforHI){
214 4 : n += getHIDensity(position);
215 : }
216 5 : if(isforHII){
217 4 : n+=getHIIDensity(position);
218 : }
219 5 : if(isforH2){
220 4 : n+= 2*getH2Density(position);
221 : }
222 :
223 : // check if all densities are deactivated and raise warning if so
224 5 : if((isforHI || isforHII || isforH2) == false){
225 2 : KISS_LOG_WARNING
226 1 : << "\nCalled getNucleonDensity on fully deactivated Ferriere \n"
227 1 : << "gas density model. The total density is set to 0.";
228 : }
229 :
230 5 : return n;
231 : }
232 :
233 1 : void Ferriere::setIsForHI(bool HI){
234 1 : isforHI = HI;
235 1 : }
236 :
237 1 : void Ferriere::setIsForHII(bool HII){
238 1 : isforHII = HII;
239 1 : }
240 :
241 1 : void Ferriere::setIsForH2(bool H2){
242 1 : isforH2 = H2;
243 1 : }
244 :
245 4 : bool Ferriere::getIsForHI(){
246 4 : return isforHI;
247 : }
248 :
249 4 : bool Ferriere::getIsForHII(){
250 4 : return isforHII;
251 : }
252 :
253 4 : bool Ferriere::getIsForH2(){
254 4 : return isforH2;
255 : }
256 :
257 0 : std::string Ferriere::getDescription() {
258 0 : std::stringstream s;
259 0 : s << "Density model Ferriere 2007:\n";
260 0 : s<< "HI component is ";
261 0 : if(!isforHI)
262 0 : s<< "not ";
263 0 : s<< "active.\nHII component is ";
264 0 : if(!isforHII)
265 0 : s<< "not ";
266 0 : s<<"active.\nH2 component is ";
267 0 : if(!isforH2)
268 0 : s<<"not ";
269 0 : s<<"active.";
270 0 : return s.str();
271 0 : }
272 :
273 : } // namespace crpropa
|
