/* * photom: a package to perform photometric calibration * Copyright (C) 2001 Michael William Richmond * * Contact: Michael William Richmond * Physics Department * Rochester Institute of Technology * 85 Lomb Memorial Drive * Rochester, NY 14623-5603 * E-mail: mwrsps@rit.edu * * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * */ /* * FILE: collate.c * * Combine information on the stars in a field from several different * data files, and from different passbands. Write a single output * file which contains information about stars which appear in all * passbands. * * Usage: * collate passband= coofile astfile passband= coofile astfile [...] * matchrad= jd= exptime= [outfile=] * * where * * passband is name of a passband, such as V or I * The two following args are required: * * coofile is the name of a ".coo" data file for that passband * contains info on each detected star's shape * * astfile is the name of a ".ast" data file for that passband * contains (RA, Dec) and instrumental mag * for each detected star * * (there can be 2-5 sets of "passband= coofile astfile" args, each * describing data in a different passband) * * * matchrad radius (in arcseconds) to use when * matching stars detected in each passband * against those in other passbands * * jd Julian Date when image was taken * * exptime exposure time (in seconds). Is assumed to be * the same for all passbands * * outfile (optional) name of output file * stdout used by default * * MWR 10/2/2000 * * 11/27/2000 Added 'init_merged' function. MWR * 12/10/2000 Added 'exptime=' argument, and code to normalize all * all magnitudes to equiv of 1-second exposure. MWR * 05/28/2001 Added support for a new column in input and output, * a 'quality_flag'. MWR */ #include #include #include #include #include "misc.h" #include "collate.h" #include "mergelist.h" #include "airmass.h" #undef DEBUG /* get some of diagnostic output */ typedef int (*PFI)(); static int read_passband_stuff(char *passband, char *coofile, char *astfile, S_PRE **det_array, int *num_det); #ifdef DEBUG static void sprePrint(S_PRE *star); static void print_pre_array(int num_star, S_PRE *array); #endif static int compare_pre_by_ra(S_PRE *star1, S_PRE *star2); static void smergedPrint(FILE *fp, S_MERGED *star); static void print_merged_array(FILE *fp, int num_star, S_MERGED *array); static int match_and_merge(int num_passbands, char passband[MAX_PASSBANDS][FILTER_LENGTH + 1], S_PRE **pre_array, int *num_pre, double radius, S_MERGED **merged_array, int *num_merged); static void init_merged(int num_stars, S_MERGED *merged_array); static int copy_to_merged(int num_stars, S_PRE *pre_array, S_MERGED *merged_array); static int find_first_index(double start_ra, int nstar, S_PRE *pre_array); static int find_closest_index(S_PRE *target_star, int nstar, S_PRE *array, int start_index, double radius); static int purify_merged_array(S_MERGED **merged_array_ptr, int num_passbands, int *num_merged); static int copy_merged_values(S_MERGED *from, S_MERGED *to); static int set_jd_airmass (S_PRE *star_array, int num_stars, double jd, double latitude, double longitude); static int normalize_exptime(S_PRE *star_array, int num_stars, double exptime); #if 0 static int compare_det_by_ra(S_DET *star1, S_DET *star2); static void print_cat_array(int num_star, S_CAT *array); static int find_matching_stars(int cat_nstar, S_CAT *cat_list, int num_detected_files, int *det_nstar, S_DET *det_list[], double radius); #endif #define USAGE "usage: collate passband= coofile astfile ... "\ "matchrad= jd= lat= long= exptime= [outfile=] " char *progname = "collate"; int main ( int argc, char *argv[] ) { char outfile[MAX_FILENAME_LEN]; char passband[MAX_PASSBANDS][FILTER_LENGTH + 1]; char coofile[MAX_PASSBANDS][MAX_FILENAME_LEN]; char astfile[MAX_PASSBANDS][MAX_FILENAME_LEN]; int i; int num_passbands; int argnum; int num_det[MAX_PASSBANDS]; int num_merged; double match_radius = MATCH_RADIUS; double jd = -1; double exptime = -1; double latitude = -999; double longitude = -999; FILE *fp; S_PRE *pre_array[MAX_PASSBANDS]; S_MERGED *merged_array; /* make special search for --version, --help, and help */ if ((argc < 2) || (strncmp(argv[1], "--help", 6) == 0) || (strncmp(argv[1], "help", 4) == 0)) { printf("%s\n", USAGE); exit(0); } if (strncmp(argv[1], "--version", 9) == 0) { printf("%s: version %s\n", progname, VERSION); exit(0); } if (argc < 9) { fprintf(stderr, "%s\n", USAGE); exit(1); } strcpy(outfile, ""); fp = stdout; /* * parse the arguments * * we start with at least two triplets giving names of passbands * and data files */ for (num_passbands = 0; ; num_passbands++) { if (num_passbands >= MAX_PASSBANDS) { shFatal("%s: num_passbands %d outside range 2-%d ", progname, num_passbands + 1, MAX_PASSBANDS); } argnum = num_passbands*3 + 1; if (strncmp(argv[argnum], "passband=", 9) != 0) { break; } /* okay, this is the start of a triplet. */ for (i = 0; i <= FILTER_LENGTH; i++) { passband[num_passbands][i] = '\0'; } strncpy(passband[num_passbands], argv[argnum] + 9, FILTER_LENGTH); strncpy(coofile[num_passbands], argv[argnum + 1], MAX_FILENAME_LEN); coofile[num_passbands][MAX_FILENAME_LEN - 1] = '\0'; strncpy(astfile[num_passbands], argv[argnum + 2], MAX_FILENAME_LEN); astfile[num_passbands][MAX_FILENAME_LEN - 1] = '\0'; #ifdef DEBUG printf(" passband %d is ..%s.. coo ..%s.. ast ..%s..\n", num_passbands, passband[num_passbands], coofile[num_passbands], astfile[num_passbands]); #endif } if ((num_passbands < 2) || (num_passbands > MAX_PASSBANDS)) { shFatal("%s: num_passbands %d outside range 2-%d ", progname, num_passbands, MAX_PASSBANDS); } /* * the remaining arguments are all preceeded by their names */ for (i = num_passbands*3 + 1; i < argc; i++) { if (strncmp(argv[i], "matchrad=", 9) == 0) { if (sscanf(argv[i] + 9, "%lf", &match_radius) != 1) { shFatal("invalid argument for matchrad argument"); } continue; } if (strncmp(argv[i], "jd=", 3) == 0) { if (sscanf(argv[i] + 3, "%lf", &jd) != 1) { shFatal("invalid argument for jd argument"); } continue; } if (strncmp(argv[i], "lat=", 4) == 0) { if (sscanf(argv[i] + 4, "%lf", &latitude) != 1) { shFatal("invalid argument for latitude argument"); } continue; } if (strncmp(argv[i], "long=", 5) == 0) { if (sscanf(argv[i] + 5, "%lf", &longitude) != 1) { shFatal("invalid argument for longitude argument"); } continue; } if (strncmp(argv[i], "exptime=", 8) == 0) { if (sscanf(argv[i] + 8, "%lf", &exptime) != 1) { shFatal("invalid argument for exptime argument"); } if (exptime <= 0) { shFatal("invalid argument for exptime: must be >= 0 seconds"); } continue; } if (strncmp(argv[i], "outfile=", 8) == 0) { if (sscanf(argv[i] + 8, "%s", outfile) != 1) { shFatal("invalid argument for outfile argument"); } if ((fp = fopen(outfile, "w")) == NULL) { shFatal("can't open file %s for output", outfile); } continue; } } /* verify that the "jd" argument has been provided */ if (jd == -1) { shFatal("jd argument not specified"); } if (jd <= 0) { shFatal("jd value is %lf must be greater than zero", jd); } /* verify that the "exptime" argument has been provided */ if (exptime == -1) { shFatal("exptime argument not specified"); } if (exptime <= 0) { shFatal("exptime value is %lf must be greater than zero", exptime); } /* * verify that the latitude and longitude arguments have been provided * and have valid values */ if (latitude == -999) { shFatal("latitude argument not specified"); } if (longitude == -999) { shFatal("longitude argument not specified"); } if ((latitude < -90) || (latitude > 90)) { shFatal("latitude value %lf is outside valid boundaries", latitude); } if ((longitude < -360) || (longitude > 360)) { shFatal("longitude value %lf is outside valid boundaries", longitude); } #ifdef DEBUG printf(" jd is is %12.7f \n", jd); printf(" latitude is is %7.2f longitude is %7.2f \n", latitude, longitude); printf(" match_radius is %7.2f arcsec, outfile is ..%s..\n", match_radius, outfile); #endif /* * convert "radius" value from arcseconds to degrees, for easier * comparison with calculated distances below */ match_radius /= 3600.0; /* read in the information for each passband, one at a time */ for (i = 0; i < num_passbands; i++) { if (read_passband_stuff(passband[i], coofile[i], astfile[i], &(pre_array[i]), &num_det[i]) != 0) { shFatal("read_passband_stuff fails for passband %s", passband[i]); } } /* * Set the "jd" and "airmass" fields properly for each star, * based on RA, Dec, and JD specified on command line */ for (i = 0; i < num_passbands; i++) { if (set_jd_airmass(pre_array[i], num_det[i], jd, latitude, longitude) != 0) { shFatal("set_jd_airmass fails for passband %s", passband[i]); } } /* * Normalize the raw magnitudes for each star to a 1-second exposure * time equivalent, using the "exptime" value from command line. */ for (i = 0; i < num_passbands; i++) { if (normalize_exptime(pre_array[i], num_det[i], exptime) != 0) { shFatal("normalize_exptime fails for passband %s", passband[i]); } } /* * Okay, now we've finished reading in all the data. * Next, we must match up the stars in each passband. * We pick the first passband specified as the "primary" one, * and match all others against it (and it alone). * * The result of the matching is a single array of S_MERGED structures. */ if (match_and_merge(num_passbands, passband, pre_array, num_det, match_radius, &merged_array, &num_merged) != 0) { shFatal("match_and_merge fails"); } /* * remove all elements from the merged array which have at least * one un-filled passband. In other words, remove all elements * which contain un-merged stars. * * We reset the value of "num_merged" as required. */ if (purify_merged_array(&merged_array, num_passbands, &num_merged) != 0) { shFatal("purify_merged_array fails"); } /* output: print the list of stars with merged data in all passbands */ print_merged_array(fp, num_merged, merged_array); /* free the arrrays we've created */ for (i = 0; i < num_passbands; i++) { shFree(pre_array[i]); } shFree(merged_array); if (fp != stdout) { fclose(fp); } return(0); } /**************************************************************************** * ROUTINE: read_passband_stuff * * Read in data from the ".coo" and ".ast" files in a passband. * Create an array of S_DET structures, and fill it with information * from the files. * * RETURNS: * 0 if all goes well * 1 otherwise */ static int read_passband_stuff ( char *passband, /* I: name of the passband */ char *coofile, /* I: name of the ".coo" data file */ char *astfile, /* I: name of the ".ast" data file */ S_PRE **pre_array, /* O: make this point to the newly allocated */ /* array of S_PRE structs we create */ int *num_stars /* O: make this hold number of stars in */ /* array */ ) { /* sanity checks */ shAssert(passband != NULL); shAssert(coofile != NULL); shAssert(astfile != NULL); #ifdef DEBUG printf("entering read_passband_stuff for passband %s\n", passband); #endif /* step 1: read data from the coofile into the array */ if (read_coo_file(coofile, passband, num_stars, pre_array) != SH_SUCCESS) { shError("read_passband_stuff: read_coo_file fails for %s", coofile); return(1); } #ifdef DEBUG print_pre_array(*num_stars, *pre_array); fflush(stdout); fflush(stderr); #endif /* step 2: read data from the astfile into the array */ if (read_ast_file(astfile, *num_stars, *pre_array) != SH_SUCCESS) { shError("read_passband_stuff: read_ast_file fails for %s", coofile); return(1); } #ifdef DEBUG printf(" after read_ast_file\n"); print_pre_array(*num_stars, *pre_array); fflush(stdout); fflush(stderr); #endif return(0); } /************************************************************************ * * ROUTINE: spreFill * * DESCRIPTION: * Fill the given S_PRE structure with the given information. * * RETURN: * nothing */ void spreFill ( S_PRE *star, /* I: pointer to struct to fill with info */ int id, /* I: ID number for star */ char *filter, /* I: name of filter in which star was observed */ double fwhm, /* I: FWHM */ double round, /* I: roundness parameter */ double sharp /* I: sharpness parameter */ ) { shAssert(star != NULL); star->id = id; shAssert(strlen(filter) <= FILTER_LENGTH); strncpy(star->filter, filter, FILTER_LENGTH); star->filter[FILTER_LENGTH] = '\0'; star->fwhm = fwhm; star->round = round; star->sharp = sharp; } #ifdef DEBUG /************************************************************************ * * ROUTINE: sprePrint * * DESCRIPTION: * Print to stdout an S_PRE structure. * * RETURN: * nothing */ void sprePrint ( S_PRE *star /* I: star to print out */ ) { shAssert(star != NULL); printf("spre %5d ", star->id); printf(" %11.5f %11.5f ", star->ra, star->dec); printf(" %2s ", star->filter); printf(" %6.3f %6.3f ", star->mag, star->magerr); printf(" %5.3f %7.3f %7.3f ", star->fwhm, star->round, star->sharp); printf(" %6u ", star->quality_flag); printf("\n"); } /************************************************************************ * * ROUTINE: print_pre_array * * DESCRIPTION: * Print to stdout the entire contents of an array of pre-merged stars. * This is probably used only in debugging. * * RETURN: * nothing */ static void print_pre_array ( int num_star, /* I: number of stars in the array */ S_PRE *array /* I: array of S_PREs to print */ ) { int i; for (i = 0; i < num_star; i++) { sprePrint(&(array[i])); } } #endif /************************************************************************ * * ROUTINE: smergedPrint * * DESCRIPTION: * Print to the given FILE an S_MERGED structure. * * RETURN: * nothing */ void smergedPrint ( FILE *fp, /* I: print to this stream */ S_MERGED *star /* I: star to print out */ ) { int i; shAssert(star != NULL); fprintf(fp, " %5d ", star->id); fprintf(fp, " %11.5f %11.5f ", star->ra, star->dec); fprintf(fp, " %14.5f %7.3f", star->jd, star->airmass); for (i = 0; i < MAX_PASSBANDS; i++) { if (star->filter[i][0] == '\0') { break; } fprintf(fp, " %3s %6.3f %6.3f ", star->filter[i], star->mag[i], star->magerr[i]); } fprintf(fp, " %8u ", star->quality_flag); fprintf(fp, "\n"); } /************************************************************************ * * ROUTINE: print_merged_array * * DESCRIPTION: * Print to the given FILE the entire contents of an array of merged stars. * This is probably used only in debugging. * * RETURN: * nothing */ static void print_merged_array ( FILE *fp, /* I: print to this stream */ int num_star, /* I: number of stars in the array */ S_MERGED *array /* I: array of S_MERGEDs to print */ ) { int i; for (i = 0; i < num_star; i++) { smergedPrint(fp, &(array[i])); } } /**************************************************************************** * ROUTINE: match_and_merge * * Given a set of arrays, one per passband, with stars detected in each * image, this routine matches objects in each passband against the * "primary" passband, and creates an array of S_MERGED structures * to hold the complete information on stars which appear in ALL passbands. * * In more detail: * * - pick the first passband as the "primary" one * - create an array of S_MERGED structures the same size as the * number of objects in the "primary" passband * (since that's the maximum possible number of merged objects) * - fill this array with info from the "primary" passband * * - for each other passband, * check to see if each detected object matches an object * in the "primary" passband * if so, add its information to the merged object * * - eliminate all objects which didn't have matches in all passbands * * More on matching: * - we consider position only as a criterion, and accept any * match within "radius" degrees. * * The strategy is * - sort the primary array by RA * - for each other array * sort the other array by RA * for each star in the other array * define start, ending points in primary array * check all primary stars between start and end * if more than one primary star matches, pick the one * which is closest to the other position * * * RETURNS: * SH_SUCCESS if all goes well * SH_GENERIC_ERROR if not */ static int match_and_merge ( int num_passbands, /* I: number of passbands of data */ char passband[MAX_PASSBANDS][FILTER_LENGTH + 1], /* I: names of each passband */ S_PRE **pre_array, /* I: array of detections in each passband */ int *num_pre, /* I: number of detections in each passband */ double radius, /* I: max dist for matching stars (degrees) */ S_MERGED **merged_array, /* O: create an array to hold merged objects */ /* and place it here */ int *num_merged /* O: number of objects in merged array */ ) { int i, num_match; int total_match; int i_passband; int start_index, closest_index, other_index, index; double start_ra; double dra, ddec, dist, other_dist; S_PRE *other_star, *match_candidate; S_MERGED *merged_star; total_match = 0; /* * first, sort all arrays by RA */ for (i = 0; i < num_passbands; i++) { qsort((void *) pre_array[i], num_pre[i], sizeof(S_PRE), (PFI) compare_pre_by_ra); #ifdef DEBUG printf(" after sorting %s array \n", passband[i]); print_pre_array(num_pre[i], pre_array[i]); #endif } /* * next, we create an array of S_MERGED structures, one for each * star in the "primary" array. We fill these structures with * info from the primary array */ *merged_array = (S_MERGED *) shMalloc(num_pre[0]*sizeof(S_MERGED)); init_merged(num_pre[0], *merged_array); if (copy_to_merged(num_pre[0], pre_array[0], *merged_array) != 0) { shError("match_and_merge: copy_to_merged fails"); return(1); } #ifdef DEBUG printf("\nHere comes merged array after copying info from %s \n", passband[0]); print_merged_array(stdout, num_pre[0], *merged_array); #endif /* * Now, we must compare each "other" passband in turn to the "primary" * passband. We look for matches between the stars in the two * passbands, and, for each match, we copy info from the "other" * passband into the merged structure */ for (i_passband = 1; i_passband < num_passbands; i_passband++) { #ifdef DEBUG printf("\nAbout to start on merging with passband %s \n", passband[i_passband]); #endif /* * Okay, time to check every star in "other" passband for * counterpart in "primary" passband. * * - pick next star in the array * - make binary search for the primary star which has * an RA just less than (other RA - radius) * - march forward through primary array until we reach * (other RA + radius) */ num_match = 0; for (i = 0; i < num_pre[i_passband]; i++) { other_star = &(pre_array[i_passband][i]); start_ra = other_star->ra - radius; /* * There are three possible cases we have to check here: * A: detected RA very small RA < radius * B: detected RA "normal" RA > radius and RA < 360-radius * C: detected RA very big RA > 360-radius * * In case A, we have to make two searches for a match: one using * the detected RA, and one using RA = (360 - detected RA). * Then we pick the better of the two matches (if there are two). * * In case B (the usual case), we need make only one search. * * In case C, we again have to make two searches: one using the * detected RA, the other using RA = radius. * Then we pick the better of the two matches. */ other_index = -1; /* we make this search in all cases A, B, C */ start_index = find_first_index(start_ra, num_pre[0], pre_array[0]); closest_index = find_closest_index(other_star, num_pre[0], pre_array[0], start_index, radius); #ifdef DEBUG if (closest_index >= 0) { printf(" closest index %5d has RA %8.4f \n", closest_index, pre_array[0][closest_index].ra); } else { printf(" closest index %5d no match \n", closest_index); } #endif /* * now check to see if we have to make a second search; if so, * we'll set "other_index" equal to the index of the star closest * to the other possibility for RA. */ if (other_star->ra < radius) { /* * case A: need to check as if the star's RA * were (360 - det_ra) */ start_ra = (360.0 - other_star->ra) - radius; start_index = find_first_index(start_ra, num_pre[0], pre_array[0]); other_index = find_closest_index(other_star, num_pre[0], pre_array[0], start_index, radius); } else if (other_star->ra > (360.0 - radius)) { /* * case C: we need to check as if the star's RA were "radius", * which means that start_ra = radius - radius = 0 */ start_ra = 0; start_index = find_first_index(start_ra, num_pre[0], pre_array[0]); other_index = find_closest_index(other_star, num_pre[0] , pre_array[0], start_index, radius); } /* * now, check to see if we found any possible matches. If we found * more than 1, pick the star which is closer to the catalog star. */ if ((closest_index != -1) && (other_index == -1)) { index = closest_index; } else if ((closest_index == -1) && (other_index != -1)) { index = other_index; } else if ((closest_index != -1) && (other_index != -1)) { /* * must check to see which is a better match. Calc distance * to each candidate match in arcsec, and pick the one which is * closer. */ match_candidate = &(pre_array[0][closest_index]); dra = (match_candidate->ra - other_star->ra)* cos(match_candidate->dec*DEGTORAD); ddec = (match_candidate->dec - other_star->dec); dist = sqrt(dra*dra + ddec*ddec)*3600.0; match_candidate = &(pre_array[0][other_index]); dra = (match_candidate->ra - other_star->ra)* cos(match_candidate->dec*DEGTORAD); ddec = (match_candidate->dec - other_star->dec); other_dist = sqrt(dra*dra + ddec*ddec)*3600.0; if (dist < other_dist) { index = closest_index; } else { index = other_index; } } else { index = -1; } if (index == -1) { #ifdef DEBUG printf(" no match found for star %4d at (%10.5f, %10.5f)\n", i, other_star->ra, other_star->dec); #endif } else { #ifdef DEBUG /* calculate distance in arcsec */ match_candidate = &(pre_array[0][index]); dra = (match_candidate->ra - other_star->ra)* cos(match_candidate->dec*DEGTORAD); ddec = (match_candidate->dec - other_star->dec); dist = sqrt(dra*dra + ddec*ddec)*3600.0; printf(" MATCH %4d (%9.5f, %8.5f) = %4d (%9.5f, %8.5f) dist %5.1f\n", i, other_star->ra, other_star->dec, index, match_candidate->ra, match_candidate->dec, dist); #endif /* set values in the MERGED structure */ merged_star = &((*merged_array)[index]); merged_star->mag[i_passband] = other_star->mag; merged_star->magerr[i_passband] = other_star->magerr; strncpy(merged_star->filter[i_passband], other_star->filter, FILTER_LENGTH); /* * when we merge, we shift the other star's quality_flag bits * to the left by (8 bits * i_passband). That means we * shouldn't use more than 8 bits for each passband, and * that we shouldn't more more than 4 passbands for 32-bit * integers. */ { int bits_to_shift; bits_to_shift = i_passband*8; merged_star->quality_flag |= (other_star->quality_flag << bits_to_shift); } num_match++; } } #ifdef DEBUG printf("\nHere comes merged array after merging info from %s \n", passband[i_passband]); print_merged_array(stdout, num_pre[0], *merged_array); #endif } /* * the number of merged elements is equal to the number of stars * in the zero'th passband. Some of these don't really contain * merged data, because there was no match. We'll fix that * later by removing all elements which don't contain matches. */ *num_merged = num_pre[0]; #ifdef DEBUG printf(" in passband %d, there are %d matches \n", i_passband, num_match); #endif total_match += num_match; #ifdef DEBUG printf(" grand total of %d matches \n", total_match); #endif return(SH_SUCCESS); } /**************************************************************************** * ROUTINE: compare_pre_by_ra * * Compare two S_PRE structures by RA. * * Return * 1 if first star has larger "RA" * 0 if the two have equal "RA" * -1 if the first has smaller "RA" * */ static int compare_pre_by_ra ( S_PRE *star1, /* I: compare "ra" field of this star .. */ S_PRE *star2 /* I: ... with that of THIS star */ ) { shAssert((star1 != NULL) && (star2 != NULL)); if (star1->ra > star2->ra) { return(1); } if (star1->ra < star2->ra) { return(-1); } return(0); } /**************************************************************************** * ROUTINE: find_first_index * * Given an array of "S_PRE" star structures, which have already * been sorted in ascending order by RA, and given some value "ra", * return the index of the star which occurs JUST before * the given "ra" value. * * We use a binary search to find the desired index. * * If there is no such star, then return the index of the last raw star. * */ static int find_first_index ( double start_ra, /* I: want star just before this RA */ int nstar, /* I: number of stars in array */ S_PRE *pre_array /* I: array of stars, sorted by RA */ ) { int top, bottom, mid; int retval; shAssert(pre_array != NULL); #ifdef DEBUG printf("find_first_index: looking for RA = %.4f\n", start_ra); #endif top = 0; if ((bottom = nstar - 1) < 0) { bottom = 0; } while (bottom - top > 2) { mid = (top + bottom)/2; #ifdef DEBUG printf(" array[%4d] RA=%.4f array[%4d] RA=%.4f array[%4d] RA=%.4f\n", top, pre_array[top].ra, mid, pre_array[mid].ra, bottom, pre_array[bottom].ra); #endif if (pre_array[mid].ra > start_ra) { bottom = mid; } else { top = mid; } } #ifdef DEBUG printf(" array[%4d] RA=%.4f array[%4d] RA=%.4f\n", top, pre_array[top].ra, bottom, pre_array[bottom].ra); #endif /* * if we get here, then the item we seek is either "top" or "bottom" * (which may point to the same item in the array). */ if (pre_array[bottom].ra > start_ra) { #ifdef DEBUG printf(" choosing array[%4d] RA=%.4f \n", top, pre_array[top].ra); #endif retval = top; } else { #ifdef DEBUG printf(" choosing array[%4d] RA=%.4f \n", bottom, pre_array[bottom].ra); #endif retval = bottom; } /* * now, check to make sure that this item has RA < start_ra. If not, * we try to back up in the cat_array until that's the case. */ while ((retval > 0) && (pre_array[retval].ra >= start_ra)) { retval--; #ifdef DEBUG printf(" backing array[%4d] RA=%.4f \n", retval, pre_array[retval].ra); #endif } #ifdef DEBUG printf(" return array[%4d] RA=%.4f \n", retval, pre_array[retval].ra); #endif return(retval); } /**************************************************************************** * ROUTINE: find_closest_index * * Given one particular "target_star", * an array "array" of "nstar" star structures, which have already * been sorted in ascending order by RA, * an index into the "array" at which to start looking, * and a matching radius (in degrees), * try to find the element of the "array" which is closest * to the given detected star. * * If the closest catalog star is within "radius" degrees, then it's a match; * if not, no match. * * RETURN: * index of matching catalog star if there is one * -1 if no match * */ static int find_closest_index ( S_PRE *target_star, /* I: want to find a match to this star */ int nstar, /* I: there are this many stars */ /* in pre_array */ S_PRE *array, /* I: array of stars, sorted by RA */ int start_index, /* I: start looking here in array */ double radius /* I: max distance (degrees) for a match */ ) { int i, closest; double other_ra, other_dec, target_ra, target_dec; double dec_factor, max_ra, min_dist2; double dra, ddec, dist2; S_PRE *other_star; shAssert(target_star != NULL); shAssert(array != NULL); target_ra = target_star->ra; target_dec = target_star->dec; dec_factor = cos(target_dec*DEGTORAD); if (dec_factor <= 0.0) { shError("closest_index: can't handle Dec = 90 degrees"); return(-1); } max_ra = target_ra + (radius/dec_factor); min_dist2 = radius*radius; closest = -1; for (i = start_index; i < nstar; i++) { other_star = &(array[i]); other_ra = other_star->ra; /* check to see if we can stop looking */ if (target_ra > max_ra) { break; } other_dec = other_star->dec; dra = (other_ra - target_ra)*dec_factor; ddec = (other_dec - target_dec); dist2 = (dra*dra + ddec*ddec); if (dist2 < min_dist2) { closest = i; min_dist2 = dist2; } } return(closest); } #if 0 /**************************************************************************** * ROUTINE: compare_det_by_ra * * Compare two S_DET structures by RA. * * Return * 1 if first star has larger "RA" * 0 if the two have equal "RA" * -1 if the first has smaller "RA" * */ static int compare_det_by_ra ( S_DET *star1, /* I: compare "ra" field of this star .. */ S_DET *star2 /* I: ... with that of THIS star */ ) { shAssert((star1 != NULL) && (star2 != NULL)); if (star1->ra > star2->ra) { return(1); } if (star1->ra < star2->ra) { return(-1); } return(0); } #if 0 /**************************************************************************** * ROUTINE: analyze_matched_pairs * * Given an array of "det_nstar" det stars, * an array of "cat_nstar" catalog stars, * the number of matching stars, * and indices into each array for the matching elements, * * do two things with the matched pairs: * * 1. calculate mean and stdev of differences between matched stars * in "identical" passbands (i.e. "V" and "V", etc). * Print out these statistics in lines which start with * COMMENT_CHAR characters * * 2. print a list of the matched stars, showing position, * difference in position (det - catalog), * and magnitudes in all defined passbands. * * However, do NOT include any magnitudes with value "BAD_MAG" in * the calculations -- they denote missing or invalid data. * * RETURN: * nothing */ static void analyze_matched_pairs ( int det_nstar, /* I: number of stars in "det" array */ struct s_mstar *det_array, /* I: array of all "det" star measurements */ int det_nband, /* I: number "det" passbands */ char *det_passband, /* I: names of "det" passbands */ int cat_nstar, /* I: number of stars in "catalog" array */ struct s_mstar *cat_array, /* I: array of all "catalog" stars */ int cat_nband, /* I: number "cat" passbands */ char *cat_passband, /* I: names of "cat" passbands */ int nmatch, /* I: number of matched pairs of stars */ int *det_match_index, /* I: index into "det" array for matches */ int *cat_match_index, /* I: index into "catalog" array for matches */ FILE *fp /* I: print results to this stream */ ) { char matched_passband[MAX_PASSBANDS]; /* names passbands in common */ int i, j; int nband; /* number of passbands in common */ int det_band_index[MAX_PASSBANDS], cat_band_index[MAX_PASSBANDS]; int count[MAX_PASSBANDS]; double dra, ddec, dist; double det_mag, cat_mag, diff; double mean, stdev; double sum[MAX_PASSBANDS], sumsq[MAX_PASSBANDS]; struct s_mstar *det_star, *cat_star; #ifdef DEBUG printf(" found a total of %d matching pairs ... \n", nmatch); #endif /* identify the passbands the "det" and "catalog" stars have in common */ nband = 0; for (i = 0; i < det_nband; i++) { for (j = 0; j < cat_nband; j++) { if (det_passband[i] == cat_passband[j]) { matched_passband[nband] = det_passband[i]; det_band_index[nband] = i; cat_band_index[nband] = j; nband++; break; } } } #ifdef DEBUG printf(" passbands in common are: \n"); for (i = 0; i < nband; i++) { printf(" det band %d = %c equals cat band %d = %c\n", det_band_index[i], det_passband[det_band_index[i]], cat_band_index[i], cat_passband[cat_band_index[i]]); } #endif for (i = 0; i < nband; i++) { sum[i] = 0.0; sumsq[i] = 0.0; count[i] = 0; } /* * walk through the set of matched pairs, calculating statistics * as we go ... */ for (i = 0; i < nmatch; i++) { det_star = &(det_array[det_match_index[i]]); cat_star = &(cat_array[cat_match_index[i]]); #ifdef DEBUG printf("det star %4d vs. cat star %4d\n", det_match_index[i], cat_match_index[i]); #endif /* calculate RA, Dec, and total distance between the stars (arcsec) */ dra = (det_star->ra - cat_star->ra)*cos(cat_star->dec*DEGTORAD)*3600.0; ddec = (det_star->dec - cat_star->dec)*3600.0; dist = sqrt(dra*dra + ddec*ddec); #ifdef DEBUG printf(" dra = %6.2f ddec = %6.2f dist = %6.2f\n", dra, ddec, dist); #endif /* * calc the difference between each corresponding magnitude, and * add to the growing sums */ for (j = 0; j < nband; j++) { det_mag = det_star->mag[det_band_index[j]]; cat_mag = cat_star->mag[cat_band_index[j]]; if ((det_mag == BAD_MAG) || (cat_mag == BAD_MAG)) { #ifdef DEBUG printf(" %c: det %6.3f cat %6.3f skipping \n", matched_passband[j], det_mag, cat_mag); #endif continue; } diff = det_mag - cat_mag; #ifdef DEBUG printf(" %c: det %6.3f cat %6.3f diff = %6.3f\n", matched_passband[j], det_mag, cat_mag, diff); #endif sum[j] += diff; sumsq[j] += diff*diff; count[j]++; } } for (i = 0; i < nband; i++) { if (count[i] == 0) { #ifdef DEBUG printf("no matches in band %c, so no stats\n", matched_passband[i]); #endif } else { if (count[i] == 1) { mean = sum[i]; stdev = 0.0; } else { mean = sum[i]/(double) count[i]; stdev = sqrt((sumsq[i] - nmatch*mean*mean)/(count[i] - 1.0)); } printf("%c passband %c N %3d mean %7.3f stdev %7.3f \n", COMMENT_CHAR, matched_passband[i], count[i], mean, stdev); } } /* * let's make a SECOND pass through the same set of matched pairs, * this time printing out one line per star with a nice format, * plus a nice header line giving labels for the columns. * We make a second pass so that the summary lines, above, will * be printed out first. */ fprintf(fp, "%c \n", COMMENT_CHAR); fprintf(fp, "%c RA, Dec are in degrees, but dRA, dDec, dist in arcsec\n", COMMENT_CHAR); fprintf(fp, "%c %10s %10s %6s %6s %6s", COMMENT_CHAR, "cat RA ", "cat Dec ", "dRA", "dDec", "dist"); for (j = 0; j < det_nband; j++) { fprintf(fp, " det %c", det_passband[j]); } for (j = 0; j < cat_nband; j++) { fprintf(fp, " cat %c", cat_passband[j]); } fprintf(fp, "\n"); for (i = 0; i < nmatch; i++) { det_star = &(det_array[det_match_index[i]]); cat_star = &(cat_array[cat_match_index[i]]); /* calculate RA, Dec, and total distance between the stars (arcsec) */ dra = (det_star->ra - cat_star->ra)*cos(cat_star->dec*DEGTORAD)*3600.0; ddec = (det_star->dec - cat_star->dec)*3600.0; dist = sqrt(dra*dra + ddec*ddec); fprintf(fp, " %10.5f %10.5f %6.2f %6.2f %6.2f ", cat_star->ra, cat_star->dec, dra, ddec, dist); /* * print out all the "det" magnitudes first, then all the "catalog" * magnitudes */ for (j = 0; j < det_nband; j++) { fprintf(fp, " %6.3f", det_star->mag[j]); } for (j = 0; j < cat_nband; j++) { fprintf(fp, " %6.3f", cat_star->mag[j]); } fprintf(fp, "\n"); } } #endif /************************************************************************ * * ROUTINE: scatFill * * DESCRIPTION: * Fill the given S_CAT structure with values ra, dec ... * but leave the 'filter' and 'mag' fields to be filled later. * * RETURN: * nothing */ void scatFill ( S_CAT *star, /* I/O: pointer to struct we fill with info */ double ra, /* I: RA value for new star (decimal degrees) */ double dec, /* I: Dec value for new star */ int num_filt, /* I: number of filters with valid info */ char *filt_array[], /* I: names of the filters */ double *mag_array /* I: magnitudes in each filter */ ) { int i, j; static int id_number = 0; shAssert(star != NULL); shAssert(num_filt > 0); shAssert(filt_array != NULL); shAssert(mag_array != NULL); star->id = id_number++; star->ra = ra; star->dec = dec; for (i = 0; i < num_filt; i++) { shAssert(filt_array[i] != NULL); shAssert(strlen(filt_array[i]) <= FILTER_LENGTH); strncpy(star->filter[i], filt_array[i], FILTER_LENGTH); star->filter[i][FILTER_LENGTH] = '\0'; star->mag[i] = mag_array[i]; } /* now, fill the remaining entries with values indicating "no data" */ for ( ; i < MAX_PASSBANDS; i++) { for (j = 0; j < FILTER_LENGTH; j++) { star->filter[i][j] = '\0'; } star->mag[i] = BAD_MAG; } star->include = 1; } /************************************************************************ * * ROUTINE: scatPrint * * DESCRIPTION: * Print to stdout an SCAT structure. * * RETURN: * nothing */ void scatPrint ( S_CAT *star /* I: star to print out */ ) { int i; shAssert(star != NULL); printf("scat %5d %10.5f %10.5f ", star->id, star->ra, star->dec); for (i = 0; i < MAX_PASSBANDS; i++) { if (star->filter[i] != NULL) { if (strcmp(star->filter[i], "") != 0) { printf(" %2s %5.2f ", star->filter[i], star->mag[i]); } } } printf(" %3d \n", star->include); } /************************************************************************ * * ROUTINE: print_cat_array * * DESCRIPTION: * Print to stdout the entire contents of an array of catalog stars. * This is probably used only in debugging. * * RETURN: * nothing */ static void print_cat_array ( int num_star, /* I: number of stars in the array */ S_CAT *array /* I: array of S_CATs to print */ ) { int i; for (i = 0; i < num_star; i++) { scatPrint(&(array[i])); } } #endif /***************************************************************************** * ROUTINE: init_merged * * Fill all elements of the S_MERGED structs with "reasonable" values * signalling emptiness. That is, strings are full of NULL characters, * numerical values are all set to -1. * * RETURNS: * 0 if all goes well * 1 if not */ static void init_merged ( int num_stars, /* I: number of stars in each array */ S_MERGED *merged_array /* O: ... we put in here */ ) { int i, j, k; S_MERGED *merged_star; shAssert(num_stars >= 0); shAssert(merged_array != NULL); for (i = 0; i < num_stars; i++) { merged_star = &(merged_array[i]); merged_star->id = -1; merged_star->ra = -1; merged_star->dec = -1; merged_star->jd = -1; merged_star->airmass = -1; for (j = 0; j < MAX_PASSBANDS; j++) { for (k = 0; k < FILTER_LENGTH + 1; k++) { merged_star->filter[j][k] = '\0'; } merged_star->mag[j] = -1; merged_star->magerr[j] = -1; } } } /***************************************************************************** * ROUTINE: copy_to_merged * * Copy information from all elements of the given array of S_PRE structs * into the given array of S_MERGED structs. * * RETURNS: * 0 if all goes well * 1 if not */ static int copy_to_merged ( int num_stars, /* I: number of stars in each array */ S_PRE *pre_array, /* I: existing structs with info ... */ S_MERGED *merged_array /* O: ... we put in here */ ) { int i, j; S_PRE *pre_star; S_MERGED *merged_star; shAssert(num_stars >= 0); shAssert(pre_array != NULL); shAssert(merged_array != NULL); for (i = 0; i < num_stars; i++) { pre_star = &(pre_array[i]); merged_star = &(merged_array[i]); merged_star->id = pre_star->id; merged_star->ra = pre_star->ra; merged_star->dec = pre_star->dec; merged_star->jd = pre_star->jd; merged_star->airmass = pre_star->airmass; for (j = 0; j < FILTER_LENGTH + 1; j++) { merged_star->filter[0][j] = '\0'; } strncpy(merged_star->filter[0], pre_star->filter, FILTER_LENGTH); merged_star->mag[0] = pre_star->mag; merged_star->magerr[0] = pre_star->magerr; merged_star->quality_flag = pre_star->quality_flag; } return(0); } /***************************************************************************** * ROUTINE: purify_merged_array * * Given an array of 'num_merged' S_MERGED structs, go through it and * identify all elements which are missing information in at least * one passband. Remove all such elements from the array. * Change the 'num_merged' argument to contain the new (smaller) * number of elements in the purified array. * * RETURNS: * 0 if all goes well * 1 if not */ static int purify_merged_array ( S_MERGED **merged_array_ptr, /* I/O: ptr to array of merged stars .. */ /* replace this with a new array */ int num_passbands, /* I: number of passbands for which a */ /* star ought to have data */ int *num_merged /* I/O: number of stars in the array. */ /* we re-set this after pruning */ /* elements from the array */ ) { int i, j, is_pure, num_pure; S_MERGED *pure_array; S_MERGED *merged_array = *merged_array_ptr; S_MERGED *old_star, *pure_star; /* * count the number of "pure" elements. We use the "filter" * fields of an element as markers. If a "filter" field * is full of NULLs, then that filter has no match. */ num_pure = 0; for (i = 0; i < *num_merged; i++) { old_star = &(merged_array[i]); is_pure = 1; for (j = 0; j < num_passbands; j++) { if (old_star->filter[j][0] == '\0') { /* star has no data in this passband */ is_pure = 0; break; } } if (is_pure == 1) { num_pure++; } } /* * now allocate for a second, "pure" array of S_MERGED structs. * We need only enough to hold all the "pure" stars. */ pure_array = (S_MERGED *) shMalloc(num_pure*sizeof(S_MERGED)); /* * copy over into the new array all "pure" elements of the old array */ num_pure = 0; for (i = 0; i < *num_merged; i++) { old_star = &(merged_array[i]); is_pure = 1; for (j = 0; j < num_passbands; j++) { if (old_star->filter[j][0] == '\0') { /* star has no data in this passband */ is_pure = 0; break; } } if (is_pure == 1) { pure_star = &(pure_array[num_pure]); copy_merged_values(old_star, pure_star); num_pure++; } } /* * free the old array, and set the 'merged_array' argument to point * to the new copy */ shFree(*merged_array_ptr); *merged_array_ptr = pure_array; *num_merged = num_pure; #ifdef DEBUG printf("\nHere comes merged array after purifying \n"); print_merged_array(stdout, num_pure, *merged_array_ptr); #endif return(0); } /***************************************************************************** * ROUTINE: copy_merged_values * * Given pointers to two S_MERGED structures, copy data from 'from' * structure to the 'to' structure. * * RETURNS: * 0 if all goes well * 1 if not */ static int copy_merged_values ( S_MERGED *from, /* I: copy data from this struct ... */ S_MERGED *to /* O: ... into this one */ ) { int i; shAssert(from != NULL); shAssert(to != NULL); to->id = from->id; to->ra = from->ra; to->dec = from->dec; to->jd = from->jd; to->airmass = from->airmass; for (i = 0; i < MAX_PASSBANDS; i++) { strncpy(to->filter[i], from->filter[i], FILTER_LENGTH); to->mag[i] = from->mag[i]; to->magerr[i] = from->magerr[i]; } to->quality_flag = from->quality_flag; return(0); } /************************************************************************* * PROCEDURE: set_jd_airmass * * DESCRIPTION: given an array of S_PRE structures, each of which * has fields for RA and Dec already set, and given * a Julian Date, calculate the airmass for each S_PRE * structure. Set the "airmass" field appropriately. * * RETURNS: * 0 if all goes well * 1 if something goes wrong */ static int set_jd_airmass ( S_PRE *star_array, /* I: array of stars for which we'll calc airmass */ int num_stars, /* I: number of stars in the array */ double jd, /* I: Julian Date of observation of all the stars */ double latitude, /* I: latitude of observatory (degrees) */ double longitude /* I: longitude of observatory */ /* (degrees West of Greenwich) */ ) { int i; double ra, dec, airmass; S_PRE *star; for (i = 0; i < num_stars; i++) { star = &(star_array[i]); ra = star->ra; dec = star->dec; airmass = tass_airmass(ra, dec, jd, latitude, longitude); star->jd = jd; star->airmass = airmass; } return(0); } /************************************************************************* * PROCEDURE: set_jd_airmass * * DESCRIPTION: given an array of S_PRE structures, and given * an exposure time (in seconds), modify the instrumental * magnitudes in each structure: calculate the magnitude * which would have been observed in a 1-second exposure. * Set the "mag" field appropriately. * * * RETURNS: * 0 if all goes well * 1 if something goes wrong */ static int normalize_exptime ( S_PRE *star_array, /* I: array of stars for which we'll calc airmass */ int num_stars, /* I: number of stars in the array */ double exptime /* I: exposure for all stars in array (seconds) */ ) { int i; double old_mag, new_mag; S_PRE *star; shAssert(exptime > 0); for (i = 0; i < num_stars; i++) { star = &(star_array[i]); old_mag = star->mag; if (old_mag != BAD_MAG) { new_mag = old_mag + 2.5*log10(exptime); star->mag = new_mag; } } return(0); }