Line data Source code
1 : //----------------------------------------------------------------------
2 : // This file is part of PARSEC (http://physik.rwth-aachen.de/parsec)
3 : // a parametrized simulation engine for cosmic rays.
4 : //
5 : // Copyright (C) 2011 Martin Erdmann, Peter Schiffer, Tobias Winchen
6 : // RWTH Aachen University, Germany
7 : // Contact: winchen@physik.rwth-aachen.de
8 : //
9 : // This program is free software: you can redistribute it and/or
10 : // modify it under the terms of the GNU General Public License as
11 : // published by the Free Software Foundation, either version 3 of
12 : // the License, or (at your option) any later version.
13 : //
14 : // This program is distributed in the hope that it will be useful,
15 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : // GNU General Public License for more details.
18 : //
19 : // You should have received a copy of the GNU General Public License
20 : // along with this program. If not, see <http://www.gnu.org/licenses/>.
21 : //----------------------------------------------------------------------
22 :
23 : #include "crpropa/magneticLens/MagneticLens.h"
24 :
25 : #include "crpropa/Random.h"
26 : #include "crpropa/Units.h"
27 :
28 : // needed for memcpy in gcc 4.3.2
29 : #include <cstring>
30 :
31 : namespace crpropa
32 : {
33 :
34 0 : void MagneticLens::loadLens(const string &filename)
35 : {
36 0 : ifstream infile(filename.c_str());
37 0 : if (!infile)
38 : {
39 0 : throw std::runtime_error("Can't read file: " + filename);
40 : }
41 : string line;
42 :
43 : string prefix;
44 0 : int sp = filename.find_last_of("/");
45 0 : if (sp >= 0)
46 : {
47 0 : prefix = filename.substr(0, sp);
48 0 : prefix.append("/");
49 : }
50 : string mfdatfile;
51 : double emin, emax;
52 :
53 : int lineCounter = 0;
54 0 : while (!infile.eof())
55 : {
56 0 : getline(infile, line);
57 0 : lineCounter++;
58 0 : if (line.find('#') == string::npos)
59 : {
60 0 : stringstream ss;
61 : ss << line;
62 0 : ss >> mfdatfile >> emin >> emax;
63 0 : if (ss.fail())
64 : {
65 0 : cerr << "WARNING! Cannot read line " << lineCounter << " of file " << filename << " :\n [" << lineCounter << " ]: " << line << " \n";
66 0 : cerr << " ... Skipping line and continue\n";
67 : }
68 : else
69 : {
70 0 : this->loadLensPart(prefix + mfdatfile, pow(10, emin) * eV, pow(10, emax) * eV);
71 : }
72 0 : }
73 : }
74 0 : }
75 :
76 12293 : bool MagneticLens::transformCosmicRay(double rigidity, double& phi,
77 : double& theta)
78 : {
79 12293 : uint32_t c = _pixelization->direction2Pix(phi, theta);
80 12293 : LensPart *lenspart = getLensPart(rigidity);
81 12293 : if (!lenspart)
82 : {
83 1 : std::cerr << "Warning. Trying to transform cosmic ray with rigidity " << rigidity / eV << " eV which is not covered by this lens!.\n";
84 2 : std::cerr << " This lens covers the range " << _minimumRigidity /eV << " eV - " << _maximumRigidity << " eV.\n";
85 1 : return false;
86 : }
87 :
88 12292 : ModelVectorType v = lenspart->getMatrix().col(c);
89 :
90 : uint32_t r;
91 :
92 : // the random number to compare with
93 12292 : double rn = Random::instance().rand();
94 :
95 : ModelVectorType::InnerIterator i(v);
96 : double cpv = 0;
97 12292 : while (i)
98 : {
99 12292 : cpv += i.value();
100 12292 : if (rn < cpv)
101 : {
102 12292 : _pixelization->pix2Direction(i.index(), phi, theta);
103 : return true;
104 : }
105 : else
106 : {
107 : ++i;
108 : }
109 : }
110 : return false;
111 : }
112 :
113 4 : bool MagneticLens::transformCosmicRay(double rigidity, Vector3d &p){
114 :
115 4 : double galacticLongitude = atan2(-p.y, -p.x);
116 4 : double galacticLatitude = M_PI / 2 - acos(-p.z/ sqrt(p.x*p.x + p.y*p.y + p.z*p.z));
117 4 : bool result = transformCosmicRay(rigidity, galacticLongitude, galacticLatitude);
118 :
119 4 : p.x = -1 * cos(galacticLongitude) * sin(M_PI / 2 - galacticLatitude);
120 4 : p.y = -1 * sin(galacticLongitude) * sin(M_PI / 2 - galacticLatitude);
121 4 : p.z = -1. * cos(M_PI / 2 - galacticLatitude);
122 :
123 4 : return result;
124 : }
125 :
126 0 : void MagneticLens::loadLensPart(const string &filename, double rigidityMin,
127 : double rigidityMax)
128 : {
129 0 : updateRigidityBounds(rigidityMin, rigidityMax);
130 :
131 0 : LensPart *p = new LensPart(filename, rigidityMin, rigidityMax);
132 : p->loadMatrixFromFile();
133 0 : _checkMatrix(p->getMatrix());
134 :
135 0 : _lensParts.push_back(p);
136 0 : }
137 :
138 3 : void MagneticLens::_checkMatrix(const ModelMatrixType &M)
139 : {
140 3 : if (M.rows() != M.cols())
141 : {
142 0 : throw std::runtime_error("Not a square Matrix!");
143 : }
144 :
145 3 : if (_pixelization)
146 : {
147 3 : if (_pixelization->nPix() != M.cols())
148 : {
149 : std::cerr << "*** ERROR ***" << endl;
150 0 : std::cerr << " Pixelization: " << _pixelization->nPix() << endl;
151 : std::cerr << " Matrix Size : " << M.cols() << endl;
152 0 : throw std::runtime_error("Matrix doesn't fit into Lense");
153 : }
154 : }
155 : else
156 : {
157 0 : uint32_t morder = Pixelization::pix2Order(M.cols());
158 0 : if (morder == 0)
159 : {
160 0 : throw std::runtime_error(
161 0 : "Matrix size doesn't match healpix scheme!");
162 : }
163 0 : _pixelization = new Pixelization(morder);
164 : }
165 3 : }
166 :
167 3 : void MagneticLens::updateRigidityBounds(double rigidityMin, double rigidityMax)
168 : {
169 3 : if (rigidityMin >= rigidityMax)
170 : {
171 0 : throw std::runtime_error("rigidityMin >= rigidityMax");
172 : }
173 3 : if (rigidityMin < _minimumRigidity)
174 : {
175 3 : _minimumRigidity = rigidityMin;
176 : }
177 :
178 3 : if (_maximumRigidity < rigidityMin)
179 : {
180 3 : _maximumRigidity = rigidityMax;
181 : }
182 3 : }
183 :
184 3 : void MagneticLens::setLensPart(const ModelMatrixType &M, double rigidityMin,
185 : double rigidityMax)
186 : {
187 3 : updateRigidityBounds(rigidityMin, rigidityMax);
188 3 : LensPart *p = new LensPart("Direct Input", rigidityMin, rigidityMax);
189 :
190 : p->setMatrix(M);
191 :
192 3 : _checkMatrix(p->getMatrix());
193 3 : _lensParts.push_back(p);
194 3 : }
195 :
196 12293 : LensPart* MagneticLens::getLensPart(double rigidity) const
197 : {
198 : const_LensPartIter i = _lensParts.begin();
199 12294 : while (i != _lensParts.end())
200 : {
201 12293 : if (((*i)->getMinimumRigidity() < rigidity / eV)
202 12293 : && ((*i)->getMaximumRigidity() >= rigidity / eV))
203 : {
204 : return (*i);
205 : }
206 : ++i;
207 : }
208 : return NULL;
209 : }
210 :
211 0 : bool MagneticLens::rigidityCovered(double rigidity) const
212 : {
213 0 : if (getLensPart(rigidity))
214 : return true;
215 : else
216 0 : return false;
217 : }
218 :
219 :
220 0 : void MagneticLens::normalizeMatrixColumns()
221 : {
222 0 : for (LensPartIter iter = _lensParts.begin(); iter != _lensParts.end();
223 : ++iter)
224 : {
225 0 : normalizeColumns((*iter)->getMatrix());
226 : }
227 0 : }
228 :
229 :
230 0 : void MagneticLens::normalizeLens()
231 : {
232 : // get maximum of sums of columns, and normalize each matrix to that
233 : double norm = 0;
234 0 : for (LensPartIter iter = _lensParts.begin(); iter != _lensParts.end();
235 : ++iter)
236 : {
237 0 : if ((*iter)->getMaximumOfSumsOfColumns() > norm)
238 : {
239 0 : norm = (*iter)->getMaximumOfSumsOfColumns();
240 : }
241 : }
242 0 : for (LensPartIter iter = _lensParts.begin(); iter != _lensParts.end();
243 : ++iter)
244 : {
245 0 : normalizeMatrix((*iter)->getMatrix(), norm);
246 : }
247 0 : _norm = norm;
248 0 : }
249 :
250 0 : void MagneticLens::normalizeLensparts()
251 : {
252 0 : for (LensPartIter iter = _lensParts.begin(); iter != _lensParts.end();
253 : ++iter)
254 : {
255 0 : double norm = (*iter)->getMaximumOfSumsOfColumns();
256 0 : normalizeMatrix((*iter)->getMatrix(), norm);
257 : }
258 0 : }
259 :
260 0 : void MagneticLens::transformModelVector(double* model, double rigidity) const
261 : {
262 0 : LensPart* lenspart = getLensPart(rigidity);
263 :
264 0 : if (!lenspart)
265 : {
266 0 : std::cerr << "Warning. Trying to transform vector with rigidity " << rigidity / eV << "eV which is not covered by this lens!.\n" << std::endl;
267 0 : return;
268 : }
269 :
270 0 : prod_up(lenspart->getMatrix(), model);
271 :
272 : }
273 :
274 :
275 :
276 : } // namespace parsec
277 :
|
