/* * GENCAT_CD - Generate CD catalog in CSV format (for processing with XMatch) * Made by Daniel E. Severin (CONICET) in 2015-2017 */ #include "common.h" /**************** * CD FUNCTIONS * ***************/ #define MAXCDSTAR 53000 struct CDstar_struct CDstar[MAXCDSTAR]; int CDstars; /* number of CD stars */ /* * read_cd - read Cordoba Durchmusterung Bytes Format Units Label Explanations 1- 2 A2 --- --- [CD] The catalog prefix 3- 5 I3 deg zone [-22/-89] The declination zone 6- 10 I5 --- num The number of the star within the zone 11 A1 --- suppl *[a-c D] star in corrigenda 12- 15 F4.1 mag mag *Estimated visual magnitude 16- 17 I2 h RAh Hours of right ascension, 1875 18- 19 I2 min RAm Minutes of right ascension, 1875 20- 23 F4.1 s RAs *Seconds of right ascension, 1875 24 A1 --- DE- [-] Sign of declination 25- 26 I2 deg DEd Degree of declination, 1875 27- 30 F4.1 arcmin DEm Minutes of declination, 1875 * * Note1: only stars from -22 to -24 are read * Note2: variable stars and nebulae are discarded */ void read_cd(char *filename) { FILE *stream; char buffer[1024]; char cell[256]; /* open file */ stream = fopen(filename, "rt"); if (stream == NULL) { perror("Error reading CD file"); exit(1); } double min_decl = 90.0, max_decl = -90.0; CDstars = 0; while (fgets(buffer, 1023, stream) != NULL) { /* read declination zone */ readfield(buffer, cell, 3, 3); int decl_ref = atoi(cell); if (decl_ref <= -25) continue; /* read number id and suppl. */ readfield(buffer, cell, 6, 5); int num_ref = atoi(cell); if (num_ref <= 0) { printf("Error in number: %d\n", num_ref); exit(1); } readfield(buffer, cell, 11, 1); char suppl_ref = cell[0]; if (suppl_ref == 'D') continue; /* <-- discard deleted stars */ /* read visual magnitude */ int deleted = 0; readfield(buffer, cell, 12, 4); float mag = atof(cell); if (mag > 12.1) { /* mag is a code */ if (mag > 29.9 && mag < 30.1) { /* discard variable star */ printf("Discarding variable star:\n"); deleted = 1; } else { if ((mag > 19.9 && mag < 20.1) || (mag > 39.9 && mag < 40.1) || (mag > 49.9 && mag < 50.1)) { /* discard nebulae */ printf("Discarding nebulae star code=%.1f:\n", mag); deleted = 2; } else { printf("Code not identified: %.1f, for %d°%d\n", mag, decl_ref, num_ref); exit(1); } } } /* read RA for B1875.0 */ readfield(buffer, cell, 16, 2); float rah = atof(cell); double RA = rah; readfield(buffer, cell, 18, 2); float ramin = atof(cell); RA += ramin / 60.0; readfield(buffer, cell, 20, 4); float raseg = atof(cell); RA += raseg / 3600.0; RA *= 15.0; /* <-- convert hs to degrees */ /* read declination for B1875.0 */ readfield(buffer, cell, 25, 2); float decldeg = atof(cell); double Decl = decldeg; readfield(buffer, cell, 27, 4); float declmin = atof(cell); Decl += declmin / 60.0; Decl = -Decl; /* <-- remember that CD decl. is always negative */ /* compute positions for J2000 using WCS library */ double RA2000 = RA; double DE2000 = Decl; double pm_ra = 0.0, pm_decl = 0.0; wcsconp(WCS_B1950, WCS_J2000, 1875.0, 0.0, 1875.0, 1875.0, &RA2000, &DE2000, &pm_ra, &pm_decl); /* update min and max declination */ if (DE2000 < min_decl) min_decl = DE2000; if (DE2000 > max_decl) max_decl = DE2000; /* store data */ if (CDstars == MAXCDSTAR) bye("Maximum number of CD stars reached!\n"); CDstar[CDstars].decl_ref = decl_ref; CDstar[CDstars].num_ref = num_ref; CDstar[CDstars].suppl_ref = suppl_ref; CDstar[CDstars].rah = rah; CDstar[CDstars].ramin = ramin; CDstar[CDstars].raseg = raseg; CDstar[CDstars].decldeg = decldeg; CDstar[CDstars].declmin = declmin; CDstar[CDstars].RA1875 = RA; CDstar[CDstars].DE1875 = Decl; CDstar[CDstars].RA2000 = RA2000; CDstar[CDstars].DE2000 = DE2000; CDstar[CDstars].mag = mag; CDstar[CDstars].dpl = false; CDstar[CDstars].color = false; CDstar[CDstars].deleted = deleted; readfield(buffer, CDstar[CDstars].reg, 3, 28); if (deleted > 0) printf("%f,%f,\"CD%d°%d%c\"\n", RA2000, DE2000, decl_ref, num_ref, suppl_ref); /* next */ CDstars++; //printf("%d°%d%c RA=%.4f Decl=%.4f mag=%.1f\n", decl_ref, num_ref, suppl_ref, RA2000, DE2000, mag); } printf("Number of Cordoba Durchmusterung stars read: %d\n", CDstars); printf("Declinations J2000: minimum = %lf°, maximum = %lf°\n", min_decl, max_decl); fclose(stream); } /* * read_flag - read color/double flag for CD stars * flag = false -> color, flag = true -> dpl */ void read_flag(char *filename, int decl_ref, bool flag) { FILE *stream; char buffer[1024]; char cell[256]; /* open file */ stream = fopen(filename, "rt"); if (stream == NULL) { perror("Error reading flag file"); exit(1); } int number = 0; while (fgets(buffer, 1023, stream) != NULL) { /* read star number */ readfield(buffer, cell, 1, 5); int num_ref = atoi(cell); /* busca la estrella CD */ int index = -1; for (int i = 0; i < CDstars; i++) { if (CDstar[i].decl_ref == decl_ref && CDstar[i].num_ref == num_ref && CDstar[i].suppl_ref == ' ') { index = i; break; } } if (index == -1) { printf("Star not found CD %d°%d\n", decl_ref, num_ref); exit(1); } /* store data */ number++; if (flag) CDstar[index].dpl = true; else CDstar[index].color = true; } if (flag) printf("Number of double stars in decl. %d: %d\n", decl_ref, number); else printf("Number of color stars in decl. %d: %d\n", decl_ref, number); fclose(stream); } /*************************************** * PPM CATALOG STRUCTURE AND FUNCTIONS * ***************************************/ #define MAXPPMSTAR 10000 struct PPMstar_struct { int ppm_ref; /* PPM ident. */ double RA1875, DE1875, RA2000, DE2000; /* right asc. and declination in B1875 and J2000 */ float vmag; /* magnitude in Johnson V, or 0.0 if it is not present */ float mag; /* magnitude in CD scale */ int cd_index; /* Index to CD */ } PPMstar[MAXPPMSTAR]; int PPMstars; /* number of PPM stars */ /* * read_ppm - read PPM (Sur + Suplemento), for those ones who have indexed in CD * Bytes Format Units Label Explanations 2- 7 I6 --- PPM *[181732/378910]+ Designation of the star 10- 18 A9 --- DM *Durchmusterung, BD or CD (10-12 decl., 13-17 number, 18 component [AB]) 20- 23 F4.1 mag Mag *Magnitude, Visual if Flag5 is 'V' 25- 26 A2 --- Sp *Spectral type 28- 29 I2 h RAh Right Ascension for the Equinox=J2000.0 and Epoch=J2000.0, on the system of FK5 31- 32 I2 min RAm Right Ascension J2000 (minutes) 34- 39 F6.3 s RAs Right Ascension J2000 (seconds) 42 A1 --- DE- Declination J2000 (sign) 43- 44 I2 deg DEd Declination for the equinox and epoch J2000.0, on the system of FK5 46- 47 I2 arcmin DEm Declination J2000 (minutes) 49- 53 F5.2 arcsec DEs Declination J2000 (seconds) 56- 62 F7.4 s/yr pmRA *Proper motion in RA, J2000 64- 69 F6.3 arcsec/yr pmDE *Proper motion in DE, J2000 71- 72 I2 --- Npos *Number of individual positions used 74- 75 I2 10mas e_RA *Mean error of RA 77- 78 I2 10mas e_DE *Mean error of DE 80- 83 F4.1 mas/yr e_pmRA *Mean error of pmRA 85- 88 F4.1 mas/yr e_pmDE *Mean error of pmDE 90- 94 F5.2 yr EpRA-1900 *Weighted mean epoch, RA and pmRA 96-100 F5.2 yr EpDE-1900 *Weighted mean epoch, DE and pmDE 102-107 I6 --- SAO [1/258997]? SAO Designation 109-114 I6 --- HD [1/359083]? Henry Draper Designation 117-125 A9 --- CPD Cape Photographic Durchmusterung Designation 127 A1 --- Flag1 *'P' - problem, 'C' - comment 128 A1 --- Flag2 *'D' - double star 129 A1 --- Flag3 * - not used 130 A1 --- Flag4 *'F' - member of FK5 131 A1 --- Flag5 *'R' - remark, 'V' - V Mag (CPC-2) -------------------------------------------------------------------------------- Note on PPM: Designation of the star in PPM South starting with No. 181732, see chapter 3 of the author's description file "desc.txt". Note on DM: Designation of the star in the Bonner Durchmusterung (for zones from -02 to -22 degrees) or in the Cordoba Durchmusterung (zones -23 to -89 degrees). Note on Mag and Sp: See chapter 3 of description file "desc.txt". Note on pmRA: Proper motion in right ascension for epoch and equinox J2000.0, on the system of FK5 given in seconds of time per Julian year. Note on pmDE: Proper motion in Declination for epoch and equinox J2000.0, on the system of FK5 given in seconds of arc per Julian year. Note on Npos: Number of individual published positions used for the derivation of the position and proper motion given. Note on e_RA: Mean error of right ascension at the mean epoch of right ascension EpRA, multiplied by the cosine of declination, given in units of 0.01 seconds of arc. Note on e_DE: Mean error of declination at the mean epoch of declination EpDE, given in units of 0.01 seconds of arc. Note on e_pmRA: Mean error of the proper motion in right ascension, multiplied by the cosine of declination, given in units of 0.001 seconds of arc per Julian year. Note on e_pmDE: mean error of the proper motion in declination given in units of 0.001 seconds of arc per Julian year Note on EpRA-1900: Weighted mean epoch of the measured positions used for the derivation of RA and pmRA, given in years since 1900.0. Note on EpDE-1900: Weighted mean epoch of the measured positions used for the derivation of DE and pmDE, given in years since 1900.0. Note on Flag1: P - problem case, preferably not to be used as astrometric reference star. C - a critical comment is given in the List of Critical Comments. Not to be used as astrometric reference star. Note on Flag2: D - double star, preferably not to be used as astrometric reference star. Note on Flag3: Not used, void for consistency with PPM (north). Note on Flag4: F - member of FK5, mostly bright stars, original FK5 data are given. Note on Flag5: R - a remark is given in the List of Remarks on Individual Stars. V - the magnitude is a photographic V magnitude copied from CPC-2. * */ void read_ppm(char *filename) { FILE *stream; char buffer[1024]; char cell[256]; stream = fopen(filename, "rt"); if (stream == NULL) { perror("Error reading PPM file"); exit(1); } PPMstars = 0; while (fgets(buffer, 1023, stream) != NULL) { /* read decl. zone of CD star, and check if it is between -23 and -24 * (note: decl. -22 correspond to BD stars) */ readfield(buffer, cell, 10, 3); int decl_ref = atoi(cell); if (decl_ref > -23 || decl_ref < -24) continue; /* read PPM ident. */ readfield(buffer, cell, 2, 6); int ppm_ref = atoi(cell); /* check flags */ readfield(buffer, cell, 127, 1); //if (cell[0] == 'P' || cell[0] == 'C') { // printf("PPM %d not considered. Reason: Star flagged as problematic.\n", ppm_ref); // continue; //} readfield(buffer, cell, 18, 1); if (cell[0] == 'B') continue; // do not consider PPM stars that are secondary components of a CD star /* read mag. V (and convert it to CD scale) */ float vmag = 0.0, mag = 0.0; readfield(buffer, cell, 131, 1); if (cell[0] == 'V') { readfield(buffer, cell, 20, 4); vmag = atof(cell); mag = convert_V(vmag); } /* read CD number and find CD star */ readfield(buffer, cell, 13, 5); int num_ref = atoi(cell); int cd_index = -1; for (int i = 0; i < CDstars; i++) { if (CDstar[i].decl_ref == decl_ref && CDstar[i].num_ref == num_ref) { cd_index = i; break; } } if (cd_index == -1) { printf("Star CD %d°%d (matched with PPM %d) not found.\n", decl_ref, num_ref, ppm_ref); continue; } /* read RA */ readfield(buffer, cell, 28, 2); double RA = atof(cell); readfield(buffer, cell, 31, 2); RA += atof(cell)/60.0; readfield(buffer, cell, 34, 6); RA += atof(cell)/3600.0; RA *= 15.0; /* convet hs to degrees */ /* read declination */ readfield(buffer, cell, 43, 2); double Decl = atof(cell); readfield(buffer, cell, 46, 2); Decl += atof(cell)/60.0; readfield(buffer, cell, 49, 5); Decl += atof(cell)/3600.0; Decl = -Decl; /* read proper motion */ readfield(buffer, cell, 56, 7); double pm_ra = atof(cell); pm_ra /= 240; /* convert from s/yr to deg/yr */ readfield(buffer, cell, 64, 6); double pm_decl = atof(cell); pm_decl /= 3600; /* convert from arcsec/yr to deg/yr */ /* convert coord. to 1875.0: since we need besselian epoch in J2000 format, we use SOFA formula: * B = 1900 + (2451545 - 2415020.31352) / 365.242198781 = 2000.001278 */ double RA1875 = RA; double DE1875 = Decl; wcsconp(WCS_J2000, WCS_B1950, 0.0, 1875.0, 2000.001278, 1875.0, &RA1875, &DE1875, &pm_ra, &pm_decl); /* store data */ if (PPMstars == MAXPPMSTAR) bye("Maximum number of PPM stars reached!\n"); PPMstar[PPMstars].ppm_ref = ppm_ref; PPMstar[PPMstars].RA2000 = RA; PPMstar[PPMstars].DE2000 = Decl; PPMstar[PPMstars].RA1875 = RA1875; PPMstar[PPMstars].DE1875 = DE1875; PPMstar[PPMstars].vmag = vmag; PPMstar[PPMstars].mag = mag; PPMstar[PPMstars].cd_index = cd_index; /* next */ PPMstars++; } printf("Number of PPM stars read: %d\n", PPMstars); fclose(stream); } /* * main */ int main(int argc, char** argv) { printf("GENCAT_CD - Generate CD catalog in CSV format.\n"); printf("Made by Daniel E. Severin (CONICET) in 2015-2017.\n"); /* we read original CD catalog */ read_cd("Cat\\cd.txt"); read_flag("Cat\\color_22.txt", -22, false); read_flag("Cat\\color_23.txt", -23, false); read_flag("Cat\\color_24.txt", -24, false); read_flag("Cat\\dpl_22.txt", -22, true); read_flag("Cat\\dpl_23.txt", -23, true); read_flag("Cat\\dpl_24.txt", -24, true); /* we read PPM catalog */ read_ppm("Cat\\ppm.txt"); /* compare CD and PPM */ for (int i = 0; i < PPMstars; i++) { int cd_index = PPMstar[i].cd_index; if (CDstar[cd_index].deleted > 0) continue; /* magnitude check */ float mag = PPMstar[i].mag; if (mag > 0.00001) { float magdif = fabs(mag - CDstar[cd_index].mag); if (magdif > MAGTOL) { printf("Discarding CD star since it differs from PPM %d (magdif = %.2f):\n", PPMstar[i].ppm_ref, magdif); printf("%f,%f,\"CD%d°%d%c\"\n", PPMstar[i].RA2000, PPMstar[i].DE2000, CDstar[cd_index].decl_ref, CDstar[cd_index].num_ref, CDstar[cd_index].suppl_ref); CDstar[cd_index].deleted = 3; continue; } } /* distance check */ double x1, y1, z1, x2, y2, z2; sph2rec(CDstar[cd_index].RA1875, CDstar[cd_index].DE1875, &x1, &y1, &z1); sph2rec(PPMstar[i].RA1875, PPMstar[i].DE1875, &x2, &y2, &z2); double dist = calcdist(calcscalar(x1, y1, z1, x2, y2, z2)); if (dist > DISTTOL) { printf("Discarding CD star since it differs from PPM %d (dist = %.2f arcmin):\n", PPMstar[i].ppm_ref, (float)(dist*60.0)); printf("%f,%f,\"CD%d°%d%c\"\n", PPMstar[i].RA2000, PPMstar[i].DE2000, CDstar[cd_index].decl_ref, CDstar[cd_index].num_ref, CDstar[cd_index].suppl_ref); CDstar[cd_index].deleted = 3; continue; } } /* save CD positions in csv format */ FILE *stream = fopen("tables\\cd_xmatch.csv", "wt"); if (stream == NULL) { perror("Error creating csv file"); exit(1); } fprintf(stream, "RAJ2000,DEJ2000,id\n"); for (int i = 0; i < CDstars; i++) { if (CDstar[i].deleted > 0) continue; /* process a non-deleted CD star */ fprintf(stream, "%f,%f,%d\n", CDstar[i].RA2000, CDstar[i].DE2000, i); } fclose(stream); /* also save CD catalog in binary format */ stream = fopen("cd_data.bin", "wb"); if (stream == NULL) { perror("Error creating bin file"); exit(1); } fwrite(&CDstars, sizeof(int), 1, stream); fwrite(CDstar, sizeof(struct CDstar_struct), CDstars, stream); fclose(stream); return 0; }