1 /*
2 SPDX-FileCopyrightText: 2011-2016 Akarsh Simha <akarsh@kde.org>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5 */
6
7 /*
8 * NOTE: I modified nomadbinfile2mysql to do this -- Akarsh
9 */
10
11 /*
12 * TODO: VERY UGLY CODE. Please fix it some day. Preferably now. This
13 * file was created from a modified C file, and needs to be recast
14 * into the C++ way of writing stuff, i.e. with classes etc.
15 */
16
17 #include "nomadbinfile2sqlite.h"
18 #include "binfile.h"
19 #include "angconversion.h"
20 #include "MeshIterator.h"
21 #include <iostream>
22 #include <string.h>
23 #include <stdio.h>
24 #define DEBUG false
25
26 using namespace std;
27
NOMADStarDataWriter(FILE * f,int HTMLevel,sqlite3 * _db,char * _db_tbl)28 NOMADStarDataWriter::NOMADStarDataWriter(FILE *f, int HTMLevel, sqlite3 *_db, char *_db_tbl)
29 {
30 db = _db;
31 DataFile = f;
32 // Create a new shiny HTMesh
33 m_Mesh = new HTMesh(HTMLevel, HTMLevel);
34 strcpy(db_tbl, _db_tbl);
35 m_HeaderRead = false;
36 }
37
~NOMADStarDataWriter()38 NOMADStarDataWriter::~NOMADStarDataWriter()
39 {
40 delete m_Mesh;
41 }
42
bswap_stardata(DeepStarData * stardata)43 void NOMADStarDataWriter::bswap_stardata(DeepStarData *stardata)
44 {
45 stardata->RA = bswap_32(stardata->RA);
46 stardata->Dec = bswap_32(stardata->Dec);
47 stardata->dRA = bswap_16(stardata->dRA);
48 stardata->dDec = bswap_16(stardata->dDec);
49 stardata->B = bswap_16(stardata->B);
50 stardata->V = bswap_16(stardata->V);
51 }
52
53 /**
54 *@short Create the table
55 */
createTable()56 bool NOMADStarDataWriter::createTable()
57 {
58 // TODO: This is not working. Investigate.
59 char create_query[2048];
60 char index_query[2048];
61 char *errorMessage = 0;
62 sprintf(create_query,
63 "CREATE TABLE IF NOT EXISTS `%s` (`Trixel` int(32) NOT NULL , `RA` double NOT NULL , `Dec` double NOT NULL "
64 ", `dRA` double NOT NULL , `dDec` double NOT NULL , `PM` double NOT NULL , `V` float NOT NULL , `B` float "
65 "NOT NULL , `Mag` float NOT NULL , `UID` INTEGER UNIQUE PRIMARY KEY AUTOINCREMENT , `Copies` int(8) NOT "
66 "NULL )",
67 db_tbl);
68
69 if (sqlite3_exec(db, create_query, 0, 0, &errorMessage) != SQLITE_OK)
70 {
71 cerr << "ERROR: Table creation failed: " << errorMessage << endl;
72 sqlite3_free(errorMessage);
73 return false;
74 }
75
76 sprintf(index_query, "CREATE INDEX IF NOT EXISTS `TrixelIndex` ON `%s`(`Trixel`)", db_tbl);
77
78 if (sqlite3_exec(db, index_query, 0, 0, &errorMessage) != SQLITE_OK)
79 {
80 cerr << "ERROR: Trixel index creation failed: " << errorMessage << endl;
81 sqlite3_free(errorMessage);
82 return false;
83 }
84
85 return true;
86 }
87
88 /**
89 *@short Calculate the final destination RA and Dec of a star with the
90 *given initial RA, Dec and proper motion rates after 'years' number
91 *of years
92 */
calculatePMCoords(double startRA,double startDec,double dRA,double dDec,double * endRA,double * endDec,float years)93 void NOMADStarDataWriter::calculatePMCoords(double startRA, double startDec, double dRA, double dDec, double *endRA,
94 double *endDec, float years)
95 {
96 // (Translated from Perl)
97 double theta0 = hour2rad(startRA);
98 double lat0 = deg2rad(startDec);
99
100 double PMperyear = sqrt(dRA * dRA + dDec * dDec);
101
102 double dir0 = (years > 0.0) ? atan2(dRA, dDec) : atan2(-dRA, -dDec); // If years < 0, we're precessing backwards
103 double PM = PMperyear * fabs(years);
104
105 double dst = deg2rad(arcsec2deg(PM / 1000.0)); // Milliarcsecond -> radian
106
107 double phi0 = M_PI / 2.0 - lat0;
108
109 double lat1 = asin(sin(lat0) * cos(dst) + cos(lat0) * sin(dst) * cos(dir0)); // Cosine rule on a sphere
110 double dtheta = atan2(sin(dir0) * sin(dst) * cos(lat0), cos(dst) - sin(lat0) * sin(lat1));
111
112 *endRA = rad2hour(theta0 + dtheta);
113 *endDec = rad2deg(lat1);
114 }
115
116 /**
117 *@short Do some calculations and insert the star data into the database
118 */
insertStarData(unsigned int trixel,const DeepStarData * const data)119 bool NOMADStarDataWriter::insertStarData(unsigned int trixel, const DeepStarData *const data)
120 {
121 char query[2048];
122 char *errorMessage = 0;
123 float mag;
124 float B, V, RA, Dec, dRA, dDec;
125
126 // Rescale the data from the structure
127 B = ((double)data->B) / 1000.0;
128 V = ((double)data->V) / 1000.0;
129 RA = ((double)data->RA) / 1000000.0;
130 Dec = ((double)data->Dec) / 100000.0;
131 dRA = ((double)data->dRA) / 1000.0;
132 dDec = ((double)data->dDec) / 1000.0;
133
134 // Check if the supplied trixel is really the trixel in which the
135 // star is in according to its RA and Dec. If that's not the case,
136 // this star is a duplicate and must be ignored
137 unsigned int originalTrixelID = m_Mesh->index(RA, Dec);
138 if (trixel != originalTrixelID)
139 return true; // Ignore this star.
140
141 // Magnitude for sorting.
142 if (V == 30.0 && B != 30.0)
143 {
144 mag = B - 1.6;
145 }
146 else
147 {
148 mag = V;
149 }
150
151 // Compute the proper motion
152 double RA1, Dec1, RA2, Dec2;
153 double PM; // Magnitude of the proper motion in milliarcseconds per year
154
155 PM = sqrt(dRA * dRA + dDec * dDec);
156
157 calculatePMCoords(RA, Dec, dRA, dDec, &RA1, &Dec1, PM_MILLENIA * -1000.);
158 calculatePMCoords(RA, Dec, dRA, dDec, &RA2, &Dec2, PM_MILLENIA * 1000.);
159
160 unsigned int TrixelList[900];
161 int ntrixels = 0;
162
163 double separation = sqrt(hour2deg(RA1 - RA2) * hour2deg(RA1 - RA2) +
164 (Dec1 - Dec2) * (Dec1 - Dec2)); // Separation in degrees // ugly.
165 if (separation > 50.0 / 60.0) // 50 arcminutes
166 {
167 m_Mesh->intersect(RA1, Dec1, RA2, Dec2);
168 MeshIterator trixels(m_Mesh);
169 while (trixels.hasNext())
170 {
171 TrixelList[ntrixels] = trixels.next();
172 ntrixels++;
173 }
174 }
175 else
176 {
177 TrixelList[0] = originalTrixelID;
178 ntrixels = 1;
179 }
180
181 if (ntrixels == 0)
182 {
183 cerr << "Ntrixels is zero in trixel " << originalTrixelID;
184 return false;
185 }
186
187 sprintf(query,
188 "INSERT INTO `%s` (`Trixel`, `RA`, `Dec`, `dRA`, `dDec`, `B`, `V`, `mag`, `PM`, `Copies`) VALUES (:Trixel, "
189 ":RA, :Dec, :dRA, :dDec, :B, :V, :mag, :PM, :Copies)",
190 db_tbl);
191 sqlite3_stmt *stmt;
192 sqlite3_prepare_v2(db, query, -1, &stmt, 0);
193
194 if (sqlite3_exec(db, "BEGIN TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
195 {
196 cerr << "SQLite INSERT INTO failed! Query was: " << endl << query << endl;
197 cerr << "Error was: " << errorMessage << endl;
198 sqlite3_free(errorMessage);
199 return false;
200 }
201
202 for (int i = 0; i < ntrixels; ++i)
203 {
204 sqlite3_bind_int(stmt, 1, TrixelList[i]);
205 sqlite3_bind_double(stmt, 2, RA);
206 sqlite3_bind_double(stmt, 3, Dec);
207 sqlite3_bind_double(stmt, 4, dRA);
208 sqlite3_bind_double(stmt, 5, dDec);
209 sqlite3_bind_double(stmt, 6, B);
210 sqlite3_bind_double(stmt, 7, V);
211 sqlite3_bind_double(stmt, 8, mag);
212 sqlite3_bind_double(stmt, 9, PM);
213 sqlite3_bind_int(stmt, 10, ((TrixelList[i] == originalTrixelID) ? ntrixels : 0));
214
215 sqlite3_step(stmt);
216 sqlite3_clear_bindings(stmt);
217 sqlite3_reset(stmt);
218 }
219
220 if (sqlite3_exec(db, "END TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
221 {
222 cerr << "SQLite INSERT INTO failed! Query was: " << endl << query << endl;
223 cerr << "Error was: " << errorMessage << endl;
224 sqlite3_free(errorMessage);
225 return false;
226 }
227
228 sqlite3_finalize(stmt);
229
230 return true;
231 }
232
truncateTable()233 bool NOMADStarDataWriter::truncateTable()
234 {
235 // Truncate table. TODO: Issue warning etc
236 char query[60];
237 char *errorMessage = 0;
238 sprintf(query, "DELETE FROM `%s`; VACUUM;", db_tbl);
239 if (sqlite3_exec(db, query, 0, 0, &errorMessage) != SQLITE_OK)
240 {
241 cerr << "Truncate table query \"" << query << "\" failed!" << endl;
242 cerr << "Error was: " << errorMessage << endl;
243 sqlite3_free(errorMessage);
244 return false;
245 }
246 return true;
247 }
248
249 /**
250 *@short Write star data to the database
251 */
writeStarDataToDB()252 bool NOMADStarDataWriter::writeStarDataToDB()
253 {
254 int8_t HTM_Level;
255 u_int16_t MSpT;
256 u_int32_t nstars;
257 u_int32_t offset;
258 unsigned int trixel;
259 DeepStarData data;
260 int16_t mag;
261
262 /*
263 // TODO: FIX THIS // FIXME
264 // We must at least check if the HTM level matches
265 fseek( DataFile, m_IndexOffset - 3, SEEK_SET );
266 fread( &HTM_Level, 1, 1, DataFile );
267 fprintf( stdout, "HTMesh Level: %d\n", HTM_Level );
268 if( HTM_Level != m_Mesh->level() ) {
269 cerr << "ERROR: HTMesh Level in file (" << HTM_Level << ") and HTM_LEVEL in program (" << m_Mesh->level() << ") differ." << endl
270 << "Please set the define directive for HTM_LEVEL in the header file correctly and rebuild."
271 << endl;
272 return false;
273 }
274 */
275 char *errorMessage = 0;
276 char query[2048];
277
278 sprintf(query,
279 "INSERT INTO `%s` (`Trixel`, `RA`, `Dec`, `dRA`, `dDec`, `B`, `V`, `mag`, `PM`, `Copies`) VALUES (:Trixel, "
280 ":RA, :Dec, :dRA, :dDec, :B, :V, :mag, :PM, :Copies)",
281 db_tbl);
282 sqlite3_stmt *stmt;
283 sqlite3_prepare_v2(db, query, -1, &stmt, 0);
284
285 if (sqlite3_exec(db, "BEGIN TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
286 {
287 cerr << "SQLite BEGIN TRANSACTION failed! Query was: " << endl << query << endl;
288 cerr << "Error was: " << errorMessage << endl;
289 sqlite3_free(errorMessage);
290 return false;
291 }
292
293 for (trixel = 0; trixel < ntrixels; ++trixel)
294 {
295 fseek(DataFile, m_IndexOffset + trixel * INDEX_ENTRY_SIZE + 4, SEEK_SET);
296 fread(&offset, 4, 1, DataFile);
297 fread(&nstars, 4, 1, DataFile);
298
299 if (DEBUG)
300 cerr << "Processing trixel " << trixel << " of " << ntrixels << endl;
301
302 /* If offset > 2^31 - 1, do the fseek in two steps */
303 if (offset > unsigned(pow2(31) - 1))
304 {
305 fseek(DataFile, unsigned(pow2(31) - 1), SEEK_SET);
306 fseek(DataFile, offset - pow2(31) + 1, SEEK_CUR);
307 }
308 else
309 fseek(DataFile, offset, SEEK_SET);
310
311 for (int i = 0; i < nstars; ++i)
312 {
313 if (DEBUG)
314 cerr << "Processing star " << i << " of " << nstars << " in trixel " << trixel << endl;
315 fread(&data, sizeof(DeepStarData), 1, DataFile);
316 if (byteswap)
317 bswap_stardata(&data);
318
319 /** CODE FROM INSERTSTARDATA PASTED HERE FOR SPEED */
320 {
321 float mag;
322 float B, V, RA, Dec, dRA, dDec;
323
324 // Rescale the data from the structure
325 B = ((double)data.B) / 1000.0;
326 V = ((double)data.V) / 1000.0;
327 RA = ((double)data.RA) / 1000000.0;
328 Dec = ((double)data.Dec) / 100000.0;
329 dRA = ((double)data.dRA) / 1000.0;
330 dDec = ((double)data.dDec) / 1000.0;
331
332 // Check if the supplied trixel is really the trixel in which the
333 // star is in according to its RA and Dec. If that's not the case,
334 // this star is a duplicate and must be ignored
335 unsigned int originalTrixelID = m_Mesh->index(hour2deg(RA), Dec);
336 if (trixel != originalTrixelID)
337 {
338 cout << "Trixel = " << trixel << ", but this is the original Trixel ID: " << originalTrixelID
339 << ". Skipping" << endl;
340 cout << "Skipped star has (RA, Dec) = " << RA << Dec << "; (dRA, dDec) = " << dRA << dDec
341 << "; and (B, V) = " << B << V << "." << endl;
342 cout << "This suspected duplicate is star " << i << "in trixel " << trixel;
343 continue;
344 }
345
346 // Magnitude for sorting.
347 if (V == 30.0 && B != 30.0)
348 {
349 mag = B - 1.6;
350 }
351 else
352 {
353 mag = V;
354 }
355
356 // Compute the proper motion
357 double RA1, Dec1, RA2, Dec2, RA1deg, RA2deg;
358 double PM; // Magnitude of the proper motion in milliarcseconds per year
359
360 PM = sqrt(dRA * dRA + dDec * dDec);
361
362 calculatePMCoords(RA, Dec, dRA, dDec, &RA1, &Dec1, PM_MILLENIA * -1000.);
363 calculatePMCoords(RA, Dec, dRA, dDec, &RA2, &Dec2, PM_MILLENIA * 1000.);
364 RA1deg = hour2deg(RA1);
365 RA2deg = hour2deg(RA2);
366
367 unsigned int TrixelList[60];
368 int nt = 0;
369
370 double separationsqr = (RA1deg - RA2deg) * (RA1deg - RA2deg) +
371 (Dec1 - Dec2) * (Dec1 - Dec2); // Separation in degrees // ugly.
372 if (separationsqr >
373 0.69) // 50 arcminutes converted to degrees, squared and rounded below = 0.69. (This has nothing to do with sex positions.)
374 {
375 m_Mesh->intersect(RA1deg, Dec1, RA2deg, Dec2);
376 MeshIterator trixels(m_Mesh);
377 while (trixels.hasNext())
378 {
379 TrixelList[nt] = trixels.next();
380 nt++;
381 }
382 }
383 else
384 {
385 TrixelList[0] = originalTrixelID;
386 nt = 1;
387 }
388
389 if (nt == 0)
390 {
391 cerr << "# of trixels is zero in trixel " << originalTrixelID;
392 return false;
393 }
394
395 for (int i = 0; i < nt; ++i)
396 {
397 sqlite3_bind_int(stmt, 1, TrixelList[i]);
398 sqlite3_bind_double(stmt, 2, RA);
399 sqlite3_bind_double(stmt, 3, Dec);
400 sqlite3_bind_double(stmt, 4, dRA);
401 sqlite3_bind_double(stmt, 5, dDec);
402 sqlite3_bind_double(stmt, 6, B);
403 sqlite3_bind_double(stmt, 7, V);
404 sqlite3_bind_double(stmt, 8, mag);
405 sqlite3_bind_double(stmt, 9, PM);
406 sqlite3_bind_int(stmt, 10, ((TrixelList[i] == originalTrixelID) ? ntrixels : 0));
407
408 sqlite3_step(stmt);
409 sqlite3_clear_bindings(stmt);
410 sqlite3_reset(stmt);
411 }
412 }
413 }
414 if (trixel % 100 == 0 && trixel != 0)
415 {
416 if (sqlite3_exec(db, "END TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
417 {
418 cerr << "SQLite END TRANSACTION failed! Query was: " << endl << query << endl;
419 cerr << "Error was: " << errorMessage << endl;
420 sqlite3_free(errorMessage);
421 return false;
422 }
423 sqlite3_finalize(stmt);
424 sqlite3_prepare_v2(db, query, -1, &stmt, 0);
425
426 if (sqlite3_exec(db, "BEGIN TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
427 {
428 cerr << "SQLite BEGIN TRANSACTION failed! Query was: " << endl << query << endl;
429 cerr << "Error was: " << errorMessage << endl;
430 sqlite3_free(errorMessage);
431 return false;
432 }
433
434 cout << "Finished trixel " << trixel << endl;
435 }
436 }
437
438 if (sqlite3_exec(db, "END TRANSACTION", 0, 0, &errorMessage) != SQLITE_OK)
439 {
440 cerr << "SQLite INSERT INTO failed! Query was: " << endl << query << endl;
441 cerr << "Error was: " << errorMessage << endl;
442 sqlite3_free(errorMessage);
443 return false;
444 }
445
446 sqlite3_finalize(stmt);
447
448 return true;
449 }
450
readFileHeader()451 bool NOMADStarDataWriter::readFileHeader()
452 {
453 int i;
454 int16_t endian_id;
455 char ASCII_text[125];
456 u_int8_t version_no;
457 u_int16_t nfields;
458
459 if (!DataFile)
460 return false;
461
462 fread(ASCII_text, 124, 1, DataFile);
463 ASCII_text[124] = '\0';
464 printf("%s", ASCII_text);
465
466 fread(&endian_id, 2, 1, DataFile);
467 if (endian_id != 0x4B53)
468 {
469 fprintf(stdout, "Byteswapping required\n");
470 byteswap = 1;
471 }
472 else
473 {
474 fprintf(stdout, "Byteswapping not required\n");
475 byteswap = 0;
476 }
477
478 fread(&version_no, 1, 1, DataFile);
479 fprintf(stdout, "Version number: %d\n", version_no);
480
481 fread(&nfields, 2, 1, DataFile);
482
483 // Just to read those many bytes
484 // TODO: Don't waste time and memory. fseek.
485 dataElement de;
486 for (i = 0; i < nfields; ++i)
487 fread(&de, sizeof(struct dataElement), 1, DataFile);
488
489 fread(&ntrixels, 4, 1, DataFile);
490 if (byteswap)
491 ntrixels = bswap_32(ntrixels);
492 fprintf(stdout, "Number of trixels reported = %d\n", ntrixels);
493
494 m_IndexOffset = ftell(DataFile);
495
496 m_HeaderRead = true;
497
498 return true;
499 }
500
write()501 bool NOMADStarDataWriter::write()
502 {
503 if (!readFileHeader())
504 return false;
505 if (!createTable())
506 return false;
507 truncateTable();
508 if (!writeStarDataToDB())
509 return false;
510 }
511
main(int argc,char * argv[])512 int main(int argc, char *argv[])
513 {
514 sqlite3 *db;
515 FILE *f;
516 char db_tbl[20];
517 char db_name[20];
518 int rc;
519
520 if (argc <= 3)
521 {
522 fprintf(stderr, "USAGE: %s <NOMAD bin file> <SQLite DB File Name> <Table Name>\n", argv[0]);
523 return 1;
524 }
525
526 strcpy(db_name, argv[2]);
527 strcpy(db_tbl, argv[3]);
528
529 f = fopen(argv[1], "r");
530
531 if (f == NULL)
532 {
533 fprintf(stderr, "ERROR: Could not open file %s for binary read.\n", argv[1]);
534 return 1;
535 }
536
537 /* Open the Database */
538 if (sqlite3_open(db_name, &db))
539 {
540 fprintf(stderr, "ERROR: Failed to open database: %s\n", sqlite3_errmsg(db));
541 sqlite3_close(db);
542 return 1;
543 }
544
545 NOMADStarDataWriter writer(f, HTM_LEVEL, db, db_tbl);
546
547 writer.write();
548
549 fclose(f);
550 sqlite3_close(db);
551 return 0;
552 }
553