Line data Source code
1 : #ifndef CRPROPA_ADVECTIONFIELD_H
2 : #define CRPROPA_ADVECTIONFIELD_H
3 :
4 :
5 : #include <string>
6 : #include <iostream>
7 : #include <cmath>
8 : #include <cstdlib>
9 : #include <sstream>
10 :
11 : #include "crpropa/Vector3.h"
12 : #include "crpropa/Referenced.h"
13 : #include "crpropa/Units.h"
14 :
15 : namespace crpropa {
16 :
17 : /**
18 : @class AdvectionField
19 : @brief Abstract base class for advection fields. These are used to model
20 : the deterministic part of the Fokker-Planck equation. The getDivergence()
21 : method is used to model the adibatic cooling/heating.
22 : */
23 : class AdvectionField: public Referenced {
24 : public:
25 : virtual ~AdvectionField() {
26 : }
27 : virtual Vector3d getField(const Vector3d &position, const double &time=0) const = 0;
28 : virtual double getDivergence(const Vector3d &position, const double &time=0) const = 0;
29 : };
30 :
31 :
32 : /**
33 : @class AdvectionFieldList
34 : @brief Advection field decorator implementing a superposition of fields.
35 : */
36 2 : class AdvectionFieldList: public AdvectionField {
37 : std::vector<ref_ptr<AdvectionField> > fields;
38 : public:
39 : void addField(ref_ptr<AdvectionField> field);
40 : Vector3d getField(const Vector3d &position, const double &time=0) const;
41 : double getDivergence(const Vector3d &position, const double &time=0) const;
42 : };
43 :
44 :
45 : /**
46 : @class UniformAdvectionField
47 : @brief Advection field with one velocity/advection-field vector.
48 : */
49 1 : class UniformAdvectionField: public AdvectionField {
50 : Vector3d value;
51 : public:
52 : UniformAdvectionField(const Vector3d &value);
53 : Vector3d getField(const Vector3d &position, const double &time=0) const;
54 : double getDivergence(const Vector3d &position, const double &time=0) const;
55 :
56 : std::string getDescription() const;
57 : };
58 :
59 :
60 : /**
61 : @class ConstantSphericalAdvectionField
62 : @brief Spherical advection field with a constant wind speed
63 : */
64 :
65 : class ConstantSphericalAdvectionField: public AdvectionField {
66 : Vector3d origin; //origin of the advection sphere
67 : double vWind; // wind velocity
68 : public:
69 : /** Constructor
70 : @param origin Origin of the advection field
71 : @param vWind Constant wind velocity
72 :
73 : */
74 :
75 : ConstantSphericalAdvectionField(const Vector3d origin, double vWind);
76 : Vector3d getField(const Vector3d &position, const double &time=0) const;
77 : double getDivergence(const Vector3d &position, const double &time=0) const;
78 :
79 : void setOrigin(const Vector3d origin);
80 : void setVWind(double vMax);
81 :
82 : Vector3d getOrigin() const;
83 : double getVWind() const;
84 :
85 : std::string getDescription() const;
86 :
87 :
88 : };
89 :
90 : /**
91 : @class SphericalAdvectionField
92 : @brief Spherical advection with a exponentially increasing and
93 : exponentially constant velocity.
94 : */
95 :
96 1 : class SphericalAdvectionField: public AdvectionField {
97 : Vector3d origin; //origin of the advection sphere
98 : double radius; //radius of the advection sphere
99 : double vMax; // maximum wind velocity
100 : double tau; // transition distance
101 : double alpha; //tuning parameter
102 : public:
103 : /** Constructor
104 : @param origin Origin of the advection sphere
105 : @param radius Radius of the advection sphere
106 : @param vMax Maximum wind velocity
107 : @param tau Transition distance
108 : @param alpha Tuning parameter
109 : */
110 : SphericalAdvectionField(const Vector3d origin, double radius, double vMax, double tau, double alpha);
111 : Vector3d getField(const Vector3d &position, const double &time=0) const;
112 : double getDivergence(const Vector3d &position, const double &time=0) const;
113 :
114 : double getV(const double &r) const;
115 :
116 : void setOrigin(const Vector3d origin);
117 : void setRadius(double radius);
118 : void setVMax(double vMax);
119 : void setTau(double tau);
120 : void setAlpha(double alpha);
121 :
122 : Vector3d getOrigin() const;
123 : double getRadius() const;
124 : double getVMax() const;
125 : double getTau() const;
126 : double getAlpha() const;
127 :
128 : std::string getDescription() const;
129 : };
130 :
131 : /**
132 : @class OneDimensionalCartesianShock
133 : @brief Advection field in x-direction with shock at x = 0 and width lShock approximated by tanh()
134 : with variable compression ratio vUp/vDown
135 : */
136 : class OneDimensionalCartesianShock: public AdvectionField {
137 : double compressionRatio; //compression ratio of shock
138 : double vUp; //upstream velocity
139 : double lShock; //shock width
140 : public:
141 : /** Constructor
142 : @param compressionRatio //compression ratio of shock
143 : @param vUp //upstream velocity
144 : @param lShock //shock width
145 : */
146 : OneDimensionalCartesianShock(double compressionRatio, double vUp, double lShock);
147 : Vector3d getField(const Vector3d &position, const double &time=0) const;
148 : double getDivergence(const Vector3d &position, const double &time=0) const;
149 :
150 : void setComp(double compressionRatio);
151 : void setVup(double vUp);
152 : void setShockwidth(double lShock);
153 :
154 : double getComp() const;
155 : double getVup() const;
156 : double getShockwidth() const;
157 :
158 : std::string getDescription() const;
159 : };
160 :
161 : /**
162 : @class OneDimensionalSphericalShock
163 : @brief Advection field in x-direction with shock at rShock and width lShock approximated by tanh()
164 : with variable compression ratio ratio vUp/vDown
165 : */
166 : class OneDimensionalSphericalShock: public AdvectionField {
167 : double compressionRatio; //compression ratio of shock
168 : double vUp; //upstream velocity
169 : double lShock; //shock width
170 : double rShock; //shock radius
171 : bool coolUpstream; //flag for upstream cooling
172 : public:
173 : /** Constructor
174 : @param compressionRatio //compression ratio of shock
175 : @param vUp //upstream velocity
176 : @param lShock //shock width
177 : @param rShock //shock radius
178 : @param coolUpstream //flag for upstream cooling
179 : */
180 : OneDimensionalSphericalShock(double rShock, double vUp, double compressionRatio, double lShock, bool coolUpstream);
181 : Vector3d getField(const Vector3d &position, const double &time=0) const;
182 : double getDivergence(const Vector3d &position, const double &time=0) const;
183 :
184 : void setComp(double compressionRatio);
185 : void setVup(double vUp);
186 : void setShockwidth(double lShock);
187 : void setShockRadius(double rShock);
188 : void setCooling(bool coolUpstream);
189 :
190 : double getComp() const;
191 : double getVup() const;
192 : double getShockwidth() const;
193 : double getShockRadius() const;
194 : bool getCooling() const;
195 :
196 : std::string getDescription() const;
197 : };
198 :
199 : /**
200 : @class ObliqueAdvectionShock
201 : @brief Advection field in x-y-direction with shock at x = 0 and width x_sh approximated by tanh()
202 : with variable compression ratio r_comp = vx_up/vx_down. The y component vy is not shocked
203 : and remains constant.
204 : */
205 : class ObliqueAdvectionShock: public AdvectionField {
206 : double compressionRatio; //compression ratio of shock
207 : double vXUp; //upstream velocity x-component
208 : double vY; //constant velocity y-component
209 : double lShock; //shock width
210 :
211 : public:
212 : /** Constructor
213 : @param compressionRatio //compression ratio of shock
214 : @param vXUp //upstream velocity x-component
215 : @param vY //constant velocity y-component
216 : @param lShock //shock width
217 :
218 : */
219 : ObliqueAdvectionShock(double compressionRatio, double vXUp, double vY, double lShock);
220 : Vector3d getField(const Vector3d &position, const double &time=0) const;
221 : double getDivergence(const Vector3d &position, const double &time=0) const;
222 :
223 : void setComp(double compressionRatio);
224 : void setVup(double vXUp);
225 : void setVy(double vY);
226 : void setShockwidth(double lShock);
227 :
228 : double getComp() const;
229 : double getVup() const;
230 : double getVy() const;
231 : double getShockwidth() const;
232 :
233 : std::string getDescription() const;
234 : };
235 :
236 : /**
237 : @class SphericalAdvectionShock
238 : @brief Spherical advection with a constant velocity for r<r_0
239 : at the the shock the velocity drops to v_0/4. followed by
240 : a decrease proportional to 1/r^2.
241 : */
242 :
243 1 : class SphericalAdvectionShock: public AdvectionField {
244 : Vector3d origin; // origin of the advection sphere
245 : double r_0; // position of the shock
246 : double v_0; // constant velocity
247 : double lambda; //transition width
248 : double r_rot; // normalization radius for rotation speed
249 : double v_phi; // rotation speed at r_rot
250 :
251 : public:
252 : /** Constructor
253 : @param origin Origin of the advection sphere
254 : @param r_0 Position of the shock
255 : @param v_0 Constant velocity (r<<r_o)
256 : @param lambda Transition width / width of the shock
257 : */
258 : SphericalAdvectionShock(const Vector3d origin, double r_0, double v_0, double lambda);
259 :
260 : Vector3d getField(const Vector3d &position, const double &time=0) const;
261 : double getDivergence(const Vector3d &position, const double &time=0) const;
262 :
263 : double g(double R) const;
264 : double g_prime(double R) const;
265 :
266 : void setOrigin(const Vector3d Origin);
267 : void setR0(double r);
268 : void setV0(double v);
269 : void setLambda(double l);
270 : void setRRot(double r);
271 : void setAzimuthalSpeed(double vPhi);
272 :
273 : Vector3d getOrigin() const;
274 : double getR0() const;
275 : double getV0() const;
276 : double getLambda() const;
277 : /**
278 : * @param r Normalization radius for rotation speed
279 : */
280 : double getRRot() const;
281 : /**
282 : * @param vPhi Rotation speed at r_rot
283 : */
284 : double getAzimuthalSpeed() const;
285 :
286 : std::string getDescription() const;
287 : };
288 :
289 : } // namespace crpropa
290 :
291 : #endif // CRPROPA_ADVECTIONFIELD_H
|
