Line data Source code
1 : /** \file phoenixSpectrum.hpp
2 : * \author Jared R. Males
3 : * \brief Utilities for working with spectra from phoenix code.
4 : * \ingroup astrophot
5 : *
6 : */
7 :
8 : #ifndef phoenixSpectrum_hpp
9 : #define phoenixSpectrum_hpp
10 :
11 : #include <fstream>
12 : #include <iostream>
13 :
14 : #include "constants.hpp"
15 : #include "../math/vectorUtils.hpp"
16 : #include "../ioutils/fileUtils.hpp"
17 : #include "../ioutils/readColumns.hpp"
18 : #include "../ioutils/stringUtils.hpp"
19 :
20 : namespace mx
21 : {
22 : namespace astro
23 : {
24 :
25 : /// A spectrum from the Phoenix model, for use with ioutils::astro::astroSpectrum.
26 : /** \ingroup astrophot_spectra
27 : */
28 : template <typename _units>
29 : struct phoenixSpectrum
30 : {
31 : typedef _units units;
32 : typedef typename units::realT realT;
33 :
34 : static const bool freq = false;
35 :
36 : /// Convert from A to SI m
37 : static constexpr realT wavelengthUnits = static_cast<realT>( 1e10 );
38 :
39 : /// Convert from erg s-1 cm-2 A-1 to SI W m-3
40 : static constexpr realT fluxUnits =
41 : static_cast<realT>( 1e7 ) / ( static_cast<realT>( 1e4 ) * static_cast<realT>( 1e10 ) );
42 :
43 : static constexpr const char *dataDirEnvVar = "PHOENIX_DATADIR";
44 :
45 : typedef std::string paramsT; ///< The parameter is a string name.
46 :
47 0 : static std::string fileName( const std::string &name )
48 : {
49 0 : return name;
50 : }
51 :
52 : static error_t
53 0 : readSpectrum( std::vector<realT> &rawLambda,
54 : std::vector<realT> &rawSpectrum,
55 : const std::string &path,
56 : const paramsT ¶ms ///< [in] the parameters are passed in case needed to construct the spectrum
57 : )
58 : {
59 0 : error_t errc = ioutils::readColumns( path, rawSpectrum );
60 0 : if( errc != error_t::noerror )
61 0 : return errc;
62 :
63 0 : errc = ioutils::readColumns( ioutils::parentPath( path ) + "/wavelength.dat", rawLambda );
64 0 : if( errc != error_t::noerror )
65 0 : return errc;
66 :
67 0 : return error_t::noerror;
68 : }
69 :
70 0 : static void scaleSpectrum( std::vector<realT> &spectrum, realT radius, realT distance )
71 : {
72 0 : for( int i = 0; i < spectrum.size(); ++i )
73 : {
74 0 : spectrum[i] *=
75 0 : pow( radius / distance * ( constants::radJupiter<units>() / constants::parsec<units>() ), 2 );
76 : }
77 0 : }
78 : };
79 :
80 : /// Read in, crop, scale, and re-write a Phoenix spectrum data file.
81 : /** For working with spectra obtained from: http://perso.ens-lyon.fr/france.allard/
82 : *
83 : * The Phoenix spectra contain many columns of line data which are not often used, are
84 : * not necessarily sorted by wavelength, and are sometimes formated so that points with leading
85 : * minus signs are not in a separate column. We also have to change the fortan D to e.
86 : *
87 : * We also want to apply the dilution factor and take the power of 10.
88 : *
89 : * This function deals with these issues, producing a two-column space-delimited
90 : * file in a specified wavelength range.
91 : *
92 : * References:
93 : * - https://phoenix.ens-lyon.fr/Grids/BT-Settl/README
94 : * - https://phoenix.ens-lyon.fr/Grids/FORMAT
95 : *
96 : * NOTE: This overwrites the input file!
97 : *
98 : *
99 : * \tparam floatT the floating point type in which to work
100 : *
101 : * \ingroup astrophot_spectra
102 : */
103 : template <typename floatT>
104 0 : void rewritePhoenixSpectrum(
105 : const std::string &filename, ///< [in.out] complete name of the file to rewrite
106 : floatT lmin, ///< [in] minimum wavelength to rewrite [microns]
107 : floatT lmax, ///< [in] maximum wavelemngth to rewrite [microns]
108 : int sepWavelength =
109 : 0, ///< [in] [optional] controls how wavelength is handled. 0=> wavelength is included in file as first column.
110 : ///< 1=> wavelength is written to a separate 'wavelength.dat' file. -1=> means don't write wavelength.
111 : floatT DF = -8.0 ///< [in] [optional] the dilution factor. See references.
112 : )
113 : {
114 : float lambda, flambda;
115 :
116 0 : std::ifstream fin;
117 :
118 0 : fin.open( filename );
119 :
120 0 : if( !fin.good() )
121 : {
122 0 : std::cerr << "Error opening file: " << filename << "\n";
123 0 : return;
124 : }
125 :
126 0 : std::vector<floatT> lambdas;
127 0 : std::vector<floatT> flambdas;
128 :
129 0 : int lineSize = 1024;
130 0 : char line[lineSize];
131 0 : std::string sline;
132 :
133 0 : fin.getline( line, lineSize );
134 0 : sline = line;
135 :
136 0 : if( sline.length() < 25 )
137 : {
138 0 : std::cerr << "Error reading file: " << filename << "\n";
139 : }
140 :
141 0 : lambda = ioutils::convertFromString<floatT>( sline.substr( 1, 12 ) );
142 :
143 : // First have to diagnose if this has the column problem
144 : int nst;
145 0 : if( sline[13] == '-' )
146 0 : nst = 13;
147 : else
148 0 : nst = 14;
149 :
150 : // convert the D to e
151 0 : size_t dpos = sline.find( 'D', nst );
152 0 : sline[dpos] = 'e';
153 :
154 0 : flambda = ioutils::convertFromString<floatT>( sline.substr( nst, 12 ) );
155 :
156 0 : if( lambda >= lmin * 1e4 && lambda <= lmax * 1e4 )
157 : {
158 0 : lambdas.push_back( lambda );
159 0 : flambdas.push_back( flambda );
160 : }
161 :
162 0 : fin.getline( line, lineSize );
163 0 : sline = line;
164 :
165 0 : while( fin.good() )
166 : {
167 0 : if( sline.length() < 25 )
168 0 : continue;
169 :
170 0 : dpos = sline.find( 'D', nst );
171 0 : sline[dpos] = 'e';
172 :
173 0 : lambda = ioutils::convertFromString<floatT>( sline.substr( 1, 12 ) );
174 0 : flambda = ioutils::convertFromString<floatT>( sline.substr( nst, 12 ) );
175 :
176 0 : if( lambda >= lmin * 1e4 && lambda <= lmax * 1e4 )
177 : {
178 0 : lambdas.push_back( lambda );
179 0 : flambdas.push_back( flambda );
180 : }
181 :
182 0 : fin.getline( line, lineSize );
183 0 : sline = line;
184 : }
185 :
186 0 : fin.close();
187 :
188 0 : std::vector<size_t> idx = math::vectorSortOrder( lambdas );
189 :
190 0 : std::ofstream fout;
191 :
192 0 : std::string fname = filename; // + ".t";
193 0 : fout.open( filename );
194 0 : fout.precision( 10 );
195 :
196 0 : for( int i = 0; i < lambdas.size(); ++i )
197 : {
198 0 : if( sepWavelength == 0 )
199 : {
200 0 : fout << lambdas[idx[i]] << " ";
201 : }
202 0 : fout << pow( 10, flambdas[idx[i]] + DF ) << "\n";
203 : }
204 0 : fout.close();
205 :
206 0 : if( sepWavelength == 1 )
207 : {
208 0 : fname = ioutils::parentPath( filename ) + "/wavelength.dat";
209 :
210 0 : fout.open( fname );
211 0 : fout.precision( 10 );
212 0 : for( int i = 0; i < lambdas.size(); ++i )
213 : {
214 0 : fout << lambdas[idx[i]] << "\n";
215 : }
216 0 : fout.close();
217 : }
218 :
219 0 : } // rewritePhoenixSpectrum
220 :
221 : /// Call rewritePhoenixSpectrum for all files in a directory.
222 : /** \ingroup astrophot_spectra
223 : */
224 : template <typename floatT>
225 0 : void rewritePhoenixSpectrumBatch( const std::string &dir, floatT lmin, floatT lmax, floatT DF = -8.0 )
226 : {
227 0 : std::vector<std::string> flist;
228 0 : ioutils::getFileNames( flist, dir, "lte", "", ".7" );
229 :
230 0 : int sepWavelength = 1;
231 0 : for( int i = 0; i < flist.size(); ++i )
232 : {
233 0 : rewritePhoenixSpectrum( flist[i], lmin, lmax, sepWavelength, DF );
234 0 : if( i == 0 )
235 0 : sepWavelength = -1;
236 : }
237 0 : }
238 :
239 : } // namespace astro
240 : } // namespace mx
241 :
242 : #endif // phoenixSpectrum_hpp
|
