DP3
DynSpecFitsWriter.h
Go to the documentation of this file.
1 // Copyright (C) 2025 ASTRON (Netherlands Institute for Radio Astronomy)
2 // SPDX-License-Identifier: GPL-3.0-or-later
3 
4 #ifndef DP3_COMMON_DYNSPEC_WRITER_H_
5 #define DP3_COMMON_DYNSPEC_WRITER_H_
6 
7 #include <aocommon/fits/fitsbase.h>
8 #include <aocommon/fits/fitswriter.h>
9 #include <aocommon/uvector.h>
10 
11 #include <array>
12 #include <cmath>
13 #include <limits>
14 #include <map>
15 #include <string>
16 #include <vector>
17 
18 #include <fitsio.h>
19 
20 namespace dp3::common {
21 
24 class DynSpecFitsWriter : private aocommon::FitsBase {
25  public:
26  DynSpecFitsWriter()
27  : n_times_(0),
28  n_channels_(0),
29  n_stokes_(4),
30  object_ra_(0.0),
31  object_dec_(0.0),
32  obs_date_(0.0),
33  time_resolution_(0.0),
34  frequency_(0.0),
35  frequency_resolution_(0.0),
36  telescope_name_(),
37  observer_(),
38  object_name_(),
39  origin_("DP3"),
40  origin_comment_("DP3/DynSpec - dynamic spectrum extracter"),
41  obs_id_() {}
42 
43  ~DynSpecFitsWriter() {
44  if (file_pointer != nullptr) {
45  int status = 0;
46  fits_close_file(file_pointer, &status);
47  }
48  }
49 
50  template <typename NumType>
51  void Write(const std::string& filename, const NumType* image) {
52  int status = 0;
53  if (file_pointer != nullptr) {
54  fits_close_file(file_pointer, &status);
55  file_pointer = nullptr;
56  }
57 
58  fits_create_file(&file_pointer, (std::string("!") + filename).c_str(),
59  &status);
60  checkStatus(status, filename);
61 
62  writeHeaders(file_pointer, filename);
63 
64  long first_pixel[3] = {1, 1, 1};
65  writeImage(file_pointer, filename, image, first_pixel);
66 
67  fits_close_file(file_pointer, &status);
68  checkStatus(status, filename);
69  file_pointer = nullptr;
70  }
71 
72  void InitializeTimeAxis(size_t n_times, double time_resolution,
73  double obs_start_time) {
74  n_times_ = n_times;
75  time_resolution_ = time_resolution;
76  obs_date_ = obs_start_time;
77  }
78 
79  void InitializeFrequencyAxis(size_t n_channels, double frequency_resolution,
80  double reference_frequency) {
81  n_channels_ = n_channels;
82  frequency_resolution_ = frequency_resolution;
83  frequency_ = reference_frequency;
84  }
85 
86  void InitializeStokesAxis(size_t n_stokes_parameters) {
87  n_stokes_ = n_stokes_parameters;
88  }
89 
90  void SetObjectName(const std::string& object_name) {
91  object_name_ = object_name;
92  }
93 
94  void SetObjectCoordinates(double ra, double dec) {
95  object_ra_ = ra;
96  object_dec_ = dec;
97  }
98 
99  void SetOrigin(const std::string& origin, const std::string& comment) {
100  origin_ = origin;
101  origin_comment_ = comment;
102  }
103 
104  void SetObsId(const std::string& obs_id) { obs_id_ = obs_id; }
105 
106  private:
107  void writeHeaders(fitsfile*& fptr, const std::string& filename) const {
108  // append image HDU
109  int bitPixInt = FLOAT_IMG;
110  std::vector<long> naxes(3);
111  naxes[0] = n_times_;
112  naxes[1] = n_channels_;
113  naxes[2] = n_stokes_;
114 
115  int status = 0;
116  fits_create_img(fptr, bitPixInt, naxes.size(), naxes.data(), &status);
117  checkStatus(status, filename);
118  const double zero = 0.0;
119  const double one = 1.0;
120  fits_write_key(fptr, TDOUBLE, "BSCALE", (void*)&one, "", &status);
121  checkStatus(status, filename);
122  fits_write_key(fptr, TDOUBLE, "BZERO", (void*)&zero, "", &status);
123  checkStatus(status, filename);
124 
125  fits_write_key(fptr, TSTRING, "BUNIT", (void*)"JY", "Units are in Jansky",
126  &status);
127  checkStatus(status, filename);
128 
129  if (!telescope_name_.empty()) {
130  fits_write_key(fptr, TSTRING, "TELESCOP", (void*)telescope_name_.c_str(),
131  "", &status);
132  checkStatus(status, filename);
133  }
134  if (!object_name_.empty()) {
135  fits_write_key(fptr, TSTRING, "OBJECT", (void*)object_name_.c_str(), "",
136  &status);
137  checkStatus(status, filename);
138  }
139  if (!obs_id_.empty()) {
140  fits_write_key(fptr, TSTRING, "OBSID", (void*)obs_id_.c_str(),
141  "Observation ID", &status);
142  checkStatus(status, filename);
143  }
144  fits_write_key(fptr, TSTRING, "ORIGIN", (void*)origin_.c_str(),
145  origin_comment_.c_str(), &status);
146  checkStatus(status, filename);
147 
148  // First axis is TIME
149  fits_write_key(fptr, TSTRING, "CTYPE1", (void*)"TIME", "", &status);
150  checkStatus(status, filename);
151  fits_write_key(fptr, TDOUBLE, "CRPIX1", (void*)&one, "", &status);
152  checkStatus(status, filename);
153  fits_write_key(fptr, TDOUBLE, "CRVAL1", (void*)&obs_date_, "", &status);
154  checkStatus(status, filename);
155  fits_write_key(fptr, TDOUBLE, "CDELT1", (void*)&time_resolution_, "",
156  &status);
157  checkStatus(status, filename);
158  fits_write_key(fptr, TSTRING, "CUNIT1", (void*)"s", "", &status);
159  checkStatus(status, filename);
160 
161  // Second axis is FREQUENCY
162  fits_write_key(fptr, TSTRING, "CTYPE2", (void*)"FREQ", "", &status);
163  checkStatus(status, filename);
164  fits_write_key(fptr, TDOUBLE, "CRPIX2", (void*)&one, "", &status);
165  checkStatus(status, filename);
166  fits_write_key(fptr, TDOUBLE, "CRVAL2", (void*)&frequency_, "", &status);
167  checkStatus(status, filename);
168  fits_write_key(fptr, TDOUBLE, "CDELT2", (void*)&frequency_resolution_, "",
169  &status);
170  checkStatus(status, filename);
171  fits_write_key(fptr, TSTRING, "CUNIT2", (void*)"Hz", "", &status);
172  checkStatus(status, filename);
173 
174  // Third axis is POLARIZATION: I, Q, U, V
175  fits_write_key(fptr, TSTRING, "CTYPE3", (void*)"STOKES", "", &status);
176  checkStatus(status, filename);
177  fits_write_key(fptr, TDOUBLE, "CRPIX3", (void*)&one, "", &status);
178  checkStatus(status, filename);
179  fits_write_key(fptr, TDOUBLE, "CRVAL3", (void*)&one, "", &status);
180  checkStatus(status, filename);
181  fits_write_key(fptr, TDOUBLE, "CDELT3", (void*)&one, "", &status);
182  checkStatus(status, filename);
183  fits_write_key(fptr, TSTRING, "CUNIT3", (void*)"", "", &status);
184  checkStatus(status, filename);
185 
186  // Store the observation date
187  int year, month, day, hour, min, sec, deciSec;
188  aocommon::FitsWriter::ModifiedJulianDateToYMD(obs_date_, year, month, day);
189  aocommon::FitsWriter::MJDToHMS(obs_date_, hour, min, sec, deciSec);
190  char dateStr[40];
191  std::sprintf(dateStr, "%d-%02d-%02dT%02d:%02d:%02d.%01d", year, month, day,
192  hour, min, sec, deciSec);
193  fits_write_key(fptr, TSTRING, "DATE-OBS", (void*)dateStr, "", &status);
194  checkStatus(status, filename);
195 
196  // Store the target's right ascension and declination
197  double ra_deg = (object_ra_ / M_PI) * 180.0,
198  dec_deg = (object_dec_ / M_PI) * 180.0;
199  fits_write_key(fptr, TDOUBLE, "OBJ-RA", (void*)&ra_deg,
200  "Target right ascension in degrees", &status);
201  checkStatus(status, filename);
202  fits_write_key(fptr, TDOUBLE, "OBJ-DEC", (void*)&dec_deg,
203  "Target declination in degrees", &status);
204  checkStatus(status, filename);
205  }
206 
207  void writeImage(fitsfile* fptr, const std::string& filename,
208  const double* image, long* currentPixel) const {
209  double null_value = std::numeric_limits<double>::max();
210  int status = 0;
211  const size_t total_size = n_times_ * n_channels_ * n_stokes_;
212  fits_write_pixnull(fptr, TDOUBLE, currentPixel, total_size,
213  const_cast<double*>(image), &null_value, &status);
214  checkStatus(status, filename);
215  }
216 
217  void writeImage(fitsfile* fptr, const std::string& filename,
218  const float* image, long* currentPixel) const {
219  float null_value = std::numeric_limits<float>::max();
220  int status = 0;
221  const size_t total_size = n_times_ * n_channels_ * n_stokes_;
222  fits_write_pixnull(fptr, TFLOAT, currentPixel, total_size,
223  const_cast<float*>(image), &null_value, &status);
224  checkStatus(status, filename);
225  }
226 
227  template <typename NumType>
228  void writeImage(fitsfile* fptr, const std::string& filename,
229  const NumType* image, long* currentPixel) const {
230  double null_value = std::numeric_limits<double>::max();
231  int status = 0;
232  const size_t total_size = n_times_ * n_channels_ * n_stokes_;
233  std::vector<double> copy(total_size);
234  for (size_t i = 0; i != total_size; ++i) copy[i] = image[i];
235  fits_write_pixnull(fptr, TDOUBLE, currentPixel, total_size, &copy[0],
236  &null_value, &status);
237  checkStatus(status, filename);
238  }
239 
240  fitsfile* file_pointer = nullptr;
241  std::size_t n_times_;
242  std::size_t n_channels_;
243  std::size_t n_stokes_;
244  double object_ra_;
245  double object_dec_;
246  double obs_date_;
247  double time_resolution_;
248  double frequency_;
249  double frequency_resolution_;
250  std::string telescope_name_;
251  std::string observer_;
252  std::string object_name_;
253  std::string origin_;
254  std::string origin_comment_;
255  std::string obs_id_;
256 };
257 } // namespace dp3::common
258 
259 #endif
dp3::common::DynSpecFitsWriter
Definition: DynSpecFitsWriter.h:24
dp3::common::DynSpecFitsWriter::SetObsId
void SetObsId(const std::string &obs_id)
Definition: DynSpecFitsWriter.h:104
dp3::common::DynSpecFitsWriter::SetOrigin
void SetOrigin(const std::string &origin, const std::string &comment)
Definition: DynSpecFitsWriter.h:99
dp3::common::DynSpecFitsWriter::SetObjectName
void SetObjectName(const std::string &object_name)
Definition: DynSpecFitsWriter.h:90
dp3::common::DynSpecFitsWriter::DynSpecFitsWriter
DynSpecFitsWriter()
Definition: DynSpecFitsWriter.h:26
dp3::common::DynSpecFitsWriter::InitializeTimeAxis
void InitializeTimeAxis(size_t n_times, double time_resolution, double obs_start_time)
Definition: DynSpecFitsWriter.h:72
dp3::common::DynSpecFitsWriter::InitializeFrequencyAxis
void InitializeFrequencyAxis(size_t n_channels, double frequency_resolution, double reference_frequency)
Definition: DynSpecFitsWriter.h:79
dp3::common::DynSpecFitsWriter::InitializeStokesAxis
void InitializeStokesAxis(size_t n_stokes_parameters)
Definition: DynSpecFitsWriter.h:86
dp3::common::DynSpecFitsWriter::Write
void Write(const std::string &filename, const NumType *image)
Definition: DynSpecFitsWriter.h:51
dp3::common::DynSpecFitsWriter::~DynSpecFitsWriter
~DynSpecFitsWriter()
Definition: DynSpecFitsWriter.h:43
dp3::common::DynSpecFitsWriter::SetObjectCoordinates
void SetObjectCoordinates(double ra, double dec)
Definition: DynSpecFitsWriter.h:94
dp3::common
Definition: BaselineSelection.h:20
Generated on Thu May 7 2026 07:44:42 for DP3 by doxygen 1.9.1