GRALE
backprojectmatrixnew.h
Go to the documentation of this file.
1 /*
2 
3  This file is a part of GRALE, a library to facilitate the simulation
4  and inversion of gravitational lenses.
5 
6  Copyright (C) 2008-2012 Jori Liesenborgs
7 
8  Contact: jori.liesenborgs@gmail.com
9 
10  This program is free software; you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation; either version 2 of the License, or
13  (at your option) any later version.
14 
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
23 
24 */
25 
26 #ifndef GRALE_BACKPROJECTMATRIXNEW_H
27 
28 #define GRALE_BACKPROJECTMATRIXNEW_H
29 
34 #include "graleconfig.h"
35 #include "deflectionmatrix.h"
36 #include "projectedimagesinterface.h"
37 
38 namespace grale
39 {
40 
41 class ImagesDataExtended;
42 class GravitationalLens;
43 
46 class GRALE_IMPORTEXPORT BackProjectMatrixNew : public ProjectedImagesInterface, public errut::ErrorBase
47 {
48 public:
49  BackProjectMatrixNew();
50  ~BackProjectMatrixNew();
51 
52  bool startInit(double z_d, double D_d, DeflectionMatrix *pDeflectionMatrix,
53  const std::vector<ImagesDataExtended *> &images,
54  const std::vector<bool> &useDeflections,
55  const std::vector<bool> &useDerivatives,
56  const std::vector<bool> &usePotentials,
57  const GravitationalLens *pBaseLens,
58  bool storeOriginalIntensities,
59  bool storeOriginalTimeDelays,
60  bool storeOriginalShearInfo,
61  bool useMassSheet);
62  bool endInit();
63 
64  void storeDeflectionMatrixResults();
65  void calculate(float scaleFactor, float massSheetFactor = 0);
66  void calculateInverseMagnifications() { calculateInverseMagnifications(m_trueFlags); }
67  void calculateShearComponents() { calculateShearComponents(m_trueFlags); }
68  void calculateConvergence() { calculateConvergence(m_trueFlags); }
69  void calculateInverseMagnifications(const std::vector<bool> &sourceMask);
70  void calculateShearComponents(const std::vector<bool> &sourceMask);
71  void calculateConvergence(const std::vector<bool> &sourceMask);
72 
73  double getAngularScale() const { return m_angularScale; }
74  double getMassSheetScale() const { return m_massSheetScale; }
75 
76  const Vector2D<float> *getBetas(int sourceNumber) const { return &(m_betas[sourceNumber][0]); }
77  const Vector2D<float> *getBetas(int sourceNumber, int imageNumber) const { return &(m_betas[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
78  const Vector2D<float> *getThetas(int sourceNumber) const { return &(m_thetas[sourceNumber][0]); }
79  const Vector2D<float> *getThetas(int sourceNumber, int imageNumber) const { return &(m_thetas[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
80  const float *getDerivativesXX(int sourceNumber) const { return &(m_axx[sourceNumber][0]); }
81  const float *getDerivativesXX(int sourceNumber, int imageNumber) const { return &(m_axx[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
82  const float *getDerivativesYY(int sourceNumber) const { return &(m_ayy[sourceNumber][0]); }
83  const float *getDerivativesYY(int sourceNumber, int imageNumber) const { return &(m_ayy[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
84  const float *getDerivativesXY(int sourceNumber) const { return &(m_axy[sourceNumber][0]); }
85  const float *getDerivativesXY(int sourceNumber, int imageNumber) const { return &(m_axy[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
86  const float *getInverseMagnifications(int sourceNumber) const { return &(m_inverseMagnifications[sourceNumber][0]); }
87  const float *getInverseMagnifications(int sourceNumber, int imageNumber) const { return &(m_inverseMagnifications[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
88  const float *getShearComponents1(int sourceNumber) const { return &(m_shearComponent1[sourceNumber][0]); }
89  const float *getShearComponents1(int sourceNumber, int imageNumber) const { return &(m_shearComponent1[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
90  const float *getShearComponents2(int sourceNumber) const { return &(m_axy[sourceNumber][0]); }
91  const float *getShearComponents2(int sourceNumber, int imageNumber) const { return &(m_axy[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
92  const float *getConvergence(int sourceNumber) const { return &(m_convergence[sourceNumber][0]); }
93  const float *getConvergence(int sourceNumber, int imageNumber) const { return &(m_convergence[sourceNumber][m_offsets[sourceNumber][imageNumber]]); }
94  float getTimeDelay(int sourceNumber, int imageNumber, int pointNumber, Vector2D<float> beta) const;
95 private:
96  DeflectionMatrix *m_pDeflectionMatrix;
97  bool m_init, m_initializing;
98 
99  std::vector<std::vector<int> > m_deflectionIndices;
100  std::vector<std::vector<int> > m_derivativeIndices;
101  std::vector<std::vector<int> > m_potentialIndices;
102 
103  std::vector<std::vector<Vector2D<double> > > m_originalPoints;
104 
105  bool m_useBaseLens;
106  std::vector<std::vector<float> > m_baseAxx, m_baseAyy, m_baseAxy;
107  std::vector<std::vector<float> > m_basePotentials;
108 #ifdef AUTODISTFRAC
109  std::vector<std::vector<Vector2D<float> > > m_baseAlphas;
110  std::vector<std::vector<Vector2D<double> > > m_baseAlphasUnscaled;
111 #else
112  std::vector<std::vector<Vector2D<float> > > m_baseBetas;
113  std::vector<std::vector<Vector2D<double> > > m_baseBetasUnscaled;
114 #endif // AUTODISTFRAC
115  std::vector<std::vector<double> > m_basePotentialsUnscaled;
116 
117  double m_Dd, m_zd;
118  double m_angularScale;
119  double m_massSheetScale;
120  float m_timeDelayScale;
121 
122  std::vector<std::vector<Vector2D<float> > > m_thetas;
123  std::vector<std::vector<Vector2D<float> > > m_subDeflectionAngles;
124  std::vector<std::vector<float> > m_subDeflectionDerivatives[3];
125  std::vector<std::vector<float> > m_subPotentialValues;
126 
127  std::vector<std::vector<Vector2D<float> > > m_betas;
128  std::vector<std::vector<float> > m_axx, m_ayy, m_axy;
129  std::vector<std::vector<float> > m_potentials;
130  std::vector<std::vector<float> > m_inverseMagnifications;
131  std::vector<std::vector<float> > m_shearComponent1;
132  std::vector<std::vector<float> > m_convergence;
133 
134  std::vector<std::vector<Vector2D<float> > > m_sheetThetas;
135  std::vector<std::vector<float> > m_sheetPotentials;
136  bool m_useMassSheet;
137 
138  std::vector<float> m_tmpBuffer;
139  std::vector<float> m_oneVector;
140  std::vector<bool> m_trueFlags;
141 };
142 
143 inline float BackProjectMatrixNew::getTimeDelay(int sourceNumber, int imageNumber, int pointNumber, Vector2D<float> beta) const
144 {
145  Vector2D<float> diff = beta - m_thetas[sourceNumber][m_offsets[sourceNumber][imageNumber]+pointNumber];
146 
147  return m_timeDelayScale*(0.5f*diff.getLengthSquared() - m_potentials[sourceNumber][m_offsets[sourceNumber][imageNumber]+pointNumber])/m_distanceFractions[sourceNumber];
148 }
149 
150 } // end namespace
151 
152 #endif // GRALE_BACKPROJECTMATRIXNEW_H
153 
Generated on Wed Dec 19 2012 09:09:23 for GRALE by   doxygen 1.8.1.1