Line data Source code
1 : /** \file cahoyAlbedos.hpp
2 : * \author Jared R. Males
3 : * \brief Type definitions for the Cahoy et al 2010 albedo spectra
4 : * \ingroup astrophot
5 : *
6 : */
7 :
8 : #ifndef mx_astro_cahoyAlbedos_hpp
9 : #define mx_astro_cahoyAlbedos_hpp
10 :
11 : #include "astroSpectrum.hpp"
12 : #include "units.hpp"
13 : #include "../error/mxErrorOld.hpp"
14 : #include "../ioutils/readColumns.hpp"
15 : #include "../math/constants.hpp"
16 :
17 : namespace mx
18 : {
19 : namespace astro
20 : {
21 :
22 : /// An albedo spectrum directly from the Cahoy et al (2010) \cite cahoy_2010 grid.
23 : /**
24 : * \ingroup astrophot_spectra
25 : */
26 : template <typename _units>
27 : struct cahoySpectrumRaw
28 : {
29 : typedef _units units;
30 : typedef typename units::realT realT;
31 :
32 : typedef struct
33 : {
34 : realT sep; ///< The separation, in AU. Can be 0.8, 2.0, 5.0, or 10.0.
35 : int metal; ///< The metallicty. Can be 1, 3, 10, or 30.
36 : realT phase; ///< The phase. Must be one of the values in the grid: 0.0, 0.1745, 0.3491, 0.5236, 0.6981,
37 : ///< 0.8727, 1.0472, 1.2217, 1.3963,1.5708,1.7453,1.9199, 2.0944, 2.2689, 2.4435, 2.6180, 2.7925,2.9671,
38 : ///< 3.139.
39 : } paramsT;
40 :
41 : static const bool freq = false;
42 :
43 : /// Convert from um to SI m
44 : static constexpr realT wavelengthUnits = static_cast<realT>( 1e6 );
45 :
46 : /// Geometric albedos are dimensionless.
47 : static constexpr realT fluxUnits = static_cast<realT>( 1 );
48 :
49 : static constexpr const char *dataDirEnvVar = "CAHOYALBEDO_DATADIR";
50 :
51 : /// The file name is constructed from the paramerters sep, metal, and phase.
52 0 : static std::string fileName( const paramsT ¶ms )
53 : {
54 0 : std::string fname;
55 :
56 0 : if( params.sep == 0.8 )
57 : {
58 0 : fname = "0.8";
59 : }
60 0 : else if( params.sep == 2.0 )
61 : {
62 0 : fname = "2.0";
63 : }
64 0 : else if( params.sep == 5.0 )
65 : {
66 0 : fname = "5";
67 : }
68 0 : else if( params.sep == 10.0 )
69 : {
70 0 : fname = "10";
71 : }
72 : else
73 : {
74 0 : mxError( "cahoyAlbedos::fileName", MXE_INVALIDARG, "invalid separation (params.sep)" );
75 0 : return "";
76 : }
77 :
78 0 : fname += "au_";
79 0 : fname += ioutils::convertToString<int>( params.metal );
80 0 : fname += "x_albedos/00_Spectral_Albedos ";
81 :
82 : char pstr[9];
83 0 : snprintf( pstr, 9, "%0.5f", params.phase );
84 :
85 0 : fname += pstr;
86 :
87 0 : fname += ".txt";
88 :
89 0 : return fname;
90 0 : }
91 :
92 : /// Read the spectrum from the file specified by path. Extra columns are discarded.
93 0 : static error_t readSpectrum( std::vector<realT> &rawLambda,
94 : std::vector<realT> &rawSpectrum,
95 : const std::string &path,
96 : const paramsT ¶ms )
97 : {
98 : ioutils::skipCol sk;
99 :
100 0 : return ioutils::readColumns( path, sk, rawLambda, sk, rawSpectrum );
101 : }
102 : };
103 :
104 : /// A class to manage interpolation in separation and phase on the albedo spectrum grid from Cahoy et al. (2010)
105 : /// \cite cahoy_2010
106 : /** Usage:
107 : \code
108 : std::vector<double> lambda;
109 : math::vectorScale(lambda, 100000, 4e-12, 600e-9); //Construct the wavelength scale
110 :
111 : cahoyGrid<units::si<double>> grid;
112 : grid.loadGrid(3, lambda); //Open the grid for 3x metallicity, pre-interpolating on the wavelength scale.
113 :
114 : grid.setParameters({1.0, 0.1}); //Set the separation to 1.0 AU, and the phase angle to 0.1 radians.
115 : grid.setSpectrum(lambda); //Perform the interpolation. Note that lambda is only used for a size check here.
116 :
117 : \endcode
118 : * After the above steps you now have the Ag spectrum at 1.0 AU and 0.1 radians phase on the lambda scale.
119 : *
120 : * \ingroup astrophot
121 : */
122 : template <typename _units>
123 : struct cahoyGrid : public baseSpectrum<typename _units::realT>
124 : {
125 : typedef _units units;
126 : typedef typename units::realT realT;
127 :
128 : typedef struct
129 : {
130 : realT sep;
131 : realT phase;
132 : } paramsT;
133 :
134 : paramsT _params;
135 :
136 : std::vector<realT> _sep;
137 : std::vector<realT> _phase;
138 :
139 : std::vector<realT> _lambda;
140 : std::vector<std::vector<std::vector<realT>>> _ag;
141 :
142 : /// Load the grid into memory for one of the metallicities and pre-interpolate onto the wavelength scale.
143 0 : int loadGrid( int metal, ///< [in] the metalicity, can be 1, 3, 10 or 30.
144 : std::vector<realT> &lambda ///< [in] the wavelength scale.
145 : )
146 : {
147 0 : _sep = { 0.8, 2.0, 2.0, 5.0, 10.0 }; // 2 is in there twice to prepare for transition distance.
148 0 : _phase = { 0.0,
149 : 0.1745,
150 : 0.3491,
151 : 0.5236,
152 : 0.6981,
153 : 0.8727,
154 : 1.0472,
155 : 1.2217,
156 : 1.3963,
157 : 1.5708,
158 : 1.7453,
159 : 1.9199,
160 : 2.0944,
161 : 2.2689,
162 : 2.4435,
163 : 2.6180,
164 : 2.7925,
165 : 2.9671,
166 : 3.139 };
167 :
168 0 : _ag.resize( _sep.size() );
169 0 : for( int i = 0; i < _ag.size(); ++i )
170 0 : _ag[i].resize( _phase.size() );
171 :
172 0 : astroSpectrum<cahoySpectrumRaw<units>> rawSpect;
173 :
174 0 : for( int i = 0; i < _sep.size(); ++i )
175 : {
176 0 : for( int j = 0; j < _phase.size(); ++j )
177 : {
178 0 : rawSpect.setParameters( { _sep[i], metal, _phase[j] } );
179 0 : if( rawSpect.setSpectrum( lambda ) != error_t::noerror )
180 : {
181 0 : mxError( "cahoGrid::loadGrid", MXE_FILERERR, "Error reading albedo spectrum." );
182 0 : return -1;
183 : }
184 :
185 0 : _ag[i][j] = rawSpect._spectrum;
186 : }
187 : }
188 :
189 0 : _sep[1] = 1.0; // The transition from cloudy to not-cloudy occurs here.
190 :
191 0 : this->_spectrum.resize( lambda.size() );
192 0 : return 0;
193 0 : }
194 :
195 : /// Set the separatio and phase of the spectrum.
196 0 : int setParameters( const paramsT ¶ms /**< [in] Struct containting the separation and phase of the spectrum */ )
197 : {
198 0 : _params = params;
199 :
200 0 : return 0;
201 : }
202 :
203 : /// Get an interpolated spectrum at arbitrary non-grid point using bilinear interpolation.
204 : /** If outside the endpoints of the grid, the grid is simply extended (i.e. constant albedo).
205 : */
206 0 : int setSpectrum(
207 : std::vector<realT>
208 : &lambda /**< [in] the wavelength grid. Must be the same as used for construction and/or openGrid*/ )
209 : {
210 : int i_l, i_u, j_l, j_u;
211 : realT phase, sep;
212 :
213 0 : if( lambda.size() != this->_spectrum.size() )
214 : {
215 0 : mxError(
216 : "cahoyGrid::setSpectrum",
217 : MXE_INVALIDARG,
218 : "wavelength grid (lambda) is not the same size as used for openGrid, or loadGrid not yet called." );
219 0 : return -1;
220 : }
221 :
222 : // Normalize phase
223 0 : phase = fmod( _params.phase, math::pi<realT>() );
224 0 : if( phase < 0 )
225 0 : phase += math::pi<realT>();
226 :
227 0 : sep = _params.sep;
228 :
229 0 : i_l = _sep.size() - 1;
230 0 : while( _sep[i_l] > sep && i_l > 0 )
231 0 : --i_l;
232 :
233 0 : i_u = 0;
234 0 : while( i_u < _sep.size() - 1 )
235 : {
236 0 : if( _sep[i_u] >= sep )
237 0 : break;
238 0 : ++i_u;
239 : }
240 :
241 0 : j_l = _phase.size() - 1;
242 0 : while( j_l > 0 )
243 : {
244 0 : if( _phase[j_l] <= phase )
245 0 : break;
246 0 : --j_l;
247 : }
248 :
249 0 : j_u = 0;
250 0 : while( j_u < _phase.size() - 1 )
251 : {
252 0 : if( _phase[j_u] >= phase )
253 0 : break;
254 0 : ++j_u;
255 : }
256 :
257 0 : realT x = sep - _sep[i_l];
258 0 : if( x < 0 )
259 0 : x = 0;
260 :
261 0 : realT y = ( phase - _phase[j_l] );
262 0 : if( y < 0 )
263 0 : y = 0;
264 :
265 : realT x0, x1;
266 :
267 0 : for( int i = 0; i < this->_spectrum.size(); ++i )
268 : {
269 0 : x0 = _ag[i_l][j_l][i];
270 0 : x1 = _ag[i_u][j_l][i];
271 :
272 0 : if( y != 0 && j_u != j_l )
273 : {
274 0 : x0 += ( _ag[i_l][j_u][i] - _ag[i_l][j_l][i] ) * y / ( _phase[j_u] - _phase[j_l] );
275 0 : x1 += ( _ag[i_u][j_u][i] - _ag[i_u][j_l][i] ) * y / ( _phase[j_u] - _phase[j_l] );
276 : }
277 0 : if( x == 0 || i_u == i_l )
278 0 : this->_spectrum[i] = x0;
279 : else
280 0 : this->_spectrum[i] = x0 + ( x1 - x0 ) * x / ( _sep[i_u] - _sep[i_l] );
281 : }
282 :
283 0 : return 0;
284 : }
285 : };
286 :
287 : } // namespace astro
288 :
289 : } // namespace mx
290 :
291 : #endif // mx_astro_cahoyAlbedos_hpp
|
