/* * 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: photom.c * * The basic idea is to calculate a photometric solution for a night, * given a set of catalog stars with known magnitudes and a set * of detected stars with instrumental magnitudes. * * Usage: * photom catalog_file racol deccol cat_nfilt V=4 B=6 * detected_file racol deccol jdcol aircol flagcol * det_nfilt a,b,c d,e,f * [ det_files ... ] * [radius=] [outfile=] * * where * * catalog_file is name of the ASCII data file containing catalog * of standard star data * * racol is the column with the catalog file in which are the * RA values (decimal degrees, J2000) * * deccol is the column with the catalog file in which are the * Dec values (decimal degrees, J2000) * * cat_nfilt is the number of "filter=column" pairs for the * reference catalog, immediately following. * Must be >= 1 * * V=4 is a "filter=column" pair; here, "V" is the name * of a filter, and "4" is the index of the * column containing magnitudes in that filter * in the catalog * * There must be exactly "cat_nfilt" such arguments * denoting passbands in the catalog file. * * * detected_file is the name of the ASCII data file containing * stars detected in some image, with instr. mags * * racol is the column in the detected file in which are the * RA values (decimal degrees, J2000) * * deccol is the column in the detected file in which are the * Dec values (decimal degrees, J2000) * * jdcol is the column in the detected file in which are the * Julian Date values * * aircol is the column in the detected file in which are the * airmass values * * flagcol is the column in the detected file in which are the * 'quality flag' values * * det_nfilt is the number of "filter=magcol,magerrcol" pairs * for the detected files, immediately following. * Must be >= 2 * * a,b,c denotes some passband in the detected file: * The values a,b,c are integers, where * a = column with passband name * b = column with magnitude val * c = column with magerr val * A comma "," separates the integers. * * There must be "det_nfilt" such arguments denoting * the passbands in the detected file. * * * det_files ... are optional additional names of files with data * on stars detected in other images. * They are assumed to have the same arrangement * of data columns and filter names as the * first file * * radius (optional) radius (in arcseconds) to use when * matching stars in catalog against the * detected stars * * outfile (optional) base name of file into which to * write diagnostic output. We create * three types of output file: * * a) a very small file containing summary * statistics: number of stars in solution, * photometric parameters in solution and * their uncertainty. Goes in file with * basename given by "outfile" arg (or * uses "photom" as default), and extension * ".coeff" * * photom.coeff by default * Jun04.coeff if outfile="Jun04" * * b) a large file(s) containing information on * all stars used in the photometric solution; * both their raw input mags, and the catalog * mags, and differences. Used for diagnostics * and debugging. * Goes in file with basename given by * "outfile" arg (or uses "photom" * by default) with extension ".comp". * * photom.comp by default * Jun04.comp if outfile="Jun04" * * c) one big output file * containing the magnitudes of all detected * stars (not just those used to make solution) * after calibration has been applied. * Goes in file with basename given by * "outfile" arg (or uses "photom" by default) * with extension ".cal". * * photom.cal by default * Jun04.cal if outfile="Jun04" * * * Note that all column numbers should be supplied by the user in 0-indexed * form -- this, the first column of data in a file is "0", the second "1", * and so forth. * * 12/10/2000: Added "mark_bad_stars" function, to knock stars with * possible saturation out of photometric solution. * * 12/16/2000: Now don't count stars with PHOTOM_TINY_WEIGHT when * reporting the statistics of photometric solution; * use new "put_star_in_solution" function to decide * whether to keep a star or not. * * 12/21/2000: Fixed bug in "apply_solution", which always applied * the zero'th passband's solution to all stars. * * 1/28/2001: Bug fixes: must add args to "assign_weight" and * "put_star_in_solution" to indicate mag indices * * 2/28/2002: Big change: program now finds zero-point for each pair * of frames individually. This makes much more sense * for a wide-field survey like TASS; the old method * was suitable for all-sky photometry. A future version * ought to have a command-line switch to control the * behavior. */ #include #include #include #include #include "precess.h" #include "misc.h" #include "photom.h" #include "readlist.h" #include "solve_linear.h" #undef DEBUG /* get some of diagnostic output */ /* * if RADEC is defined, coords are in (RA, Dec), decimal degrees. * and we use rough spherical trig to find dist. * * if " " not defined, coords are in arbitrary (x, y) units, * and we use planar trig to find dist. */ #define RADEC typedef int (*PFI)(); static int compare_cat_by_ra(S_CAT *star1, S_CAT *star2); static int compare_det_by_ra(S_DET *star1, S_DET *star2); static void print_cat_array(int num_star, S_CAT *array); static void print_det_array(int num_star, S_DET *array); static int find_first_index(double start_ra, int cat_nstar, S_CAT *cat_array); static int find_closest_index(S_DET *det_star, int cat_nstar, S_CAT *cat_array, int start_index, double radius); 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); static int mark_bad_stars(int num_detected_files, int *det_nstar, S_DET *det_array[]); static int analyze_matched_pairs(int cat_nstar, S_CAT *cat_array, int num_detected_files, int *det_nstar, S_DET *det_list[], S_SOLUTION solution[]); static double calc_dist_between (S_CAT *cat_star, S_DET *det_star, double *dra, double *ddec); static double assign_weight(S_DET *det_star, int pri_index, int sec_index); static int print_coeffs(char *basename, S_SOLUTION solution[], int num_detected_files); static int print_comps (char *basename, int cat_nstar, S_CAT *cat_array, int num_detected_files, int *det_nstar, S_DET *det_array[], S_SOLUTION solution[]); static void apply_solution(S_DET *det_star, S_SOLUTION solution[], int image_index, int pri_index, int sec_index, double *calib_mag, double *calib_magerr); static int calibrate_detected(char *basename, int cat_nstar, S_CAT *cat_array, int num_detected_files, char det_file_name[MAX_DETECTED_FILES][MAX_FILENAME_LEN+1], int *det_nstar, S_DET *det_array[], S_SOLUTION solution[]); #if 0 static int set_calib_name(char *input_name, char *output_name); #endif static int find_passband_indices(int num_cat, S_CAT *cat_array, int num_detected_files, int *num_det, S_DET *det_array[], int pri_index, int *sec_index, int *cat_index); static void initialize_solution(S_SOLUTION solution[], int num_detected_files); static void free_solution(S_SOLUTION solution[]); static void initialize_passbands(S_PASSBANDS *passband_library); static void set_passband_values(int index, int pri_value, int sec_value, int cat_value, char *pri_passband, char *sec_passband, char *cat_passband); static void get_passband_values(int index, int *pri_value, int *sec_value, int *cat_value, char *pri_passband, char *sec_passband, char *cat_passband); static int valid_indices(int pri_index, int sec_index, int cat_index); static int put_star_in_solution(S_DET *det_star, S_CAT *cat_array, int pri_index, int sec_index, int cat_index); #define USAGE "usage: photom catalog_file racol deccol cat_nfilt V=3 ... "\ "detected_file racol deccol jdcol aircol flagcol " \ "det_nfilt a,b,c ... " \ "[det_files ... ] " \ "[radius=] [outfile=] " char *progname = "photom"; S_PASSBANDS passband_library; int main ( int argc, char *argv[] ) { char cat_file[MAX_FILENAME_LEN+1]; char det_file[MAX_DETECTED_FILES][MAX_FILENAME_LEN+1]; char outfile[MAX_FILENAME_LEN+1]; char filter_char; char *cat_filter_array[MAX_PASSBANDS]; int i; int cur_arg_index; int filter_col; int cat_racol, cat_deccol; int cat_numfilt; int cat_magcol_array[MAX_PASSBANDS]; int cat_nstar; int det_racol, det_deccol, det_jdcol, det_aircol, det_flagcol; int det_numfilt, det_magcol, det_magerrcol; int num_detected_files = 0; int det_nstar[MAX_DETECTED_FILES]; int filtercol_array[MAX_PASSBANDS]; int magcol_array[MAX_PASSBANDS], magerrcol_array[MAX_PASSBANDS]; double match_radius = MATCH_RADIUS; S_CAT *cat_array = NULL; S_DET *det_array[MAX_DETECTED_FILES]; S_SOLUTION solution[MAX_PASSBANDS]; /* 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 < 13) { fprintf(stderr, "%s\n", USAGE); exit(1); } strcpy(outfile, progname); initialize_passbands(&passband_library); /* parse the arguments */ /* we start with those pertaining to the "catalog" starlist file */ strncpy(cat_file, argv[1], MAX_FILENAME_LEN); cat_file[MAX_FILENAME_LEN] = '\0'; if (sscanf(argv[2], "%d", &cat_racol) != 1) { shFatal("invalid argument for column for RA values in catalog file"); } if (sscanf(argv[3], "%d", &cat_deccol) != 1) { shFatal("invalid argument for column for Dec values in catalog file"); } if (sscanf(argv[4], "%d", &cat_numfilt) != 1) { shFatal("invalid argument for cat_numfilt"); } cur_arg_index = 5; for (i = 0; i < cat_numfilt; i++) { if (sscanf(argv[cur_arg_index], "%c=%d", &filter_char, &filter_col) != 2) { shFatal("invalid argument %s for catalog primary filter=column", argv[cur_arg_index]); } cat_filter_array[i] = shMalloc(FILTER_LENGTH + 1); sprintf(cat_filter_array[i], "%c", filter_char); cat_magcol_array[i] = filter_col; cur_arg_index++; } #ifdef DEBUG printf(" catfile %20s has RA %2d, Dec %2d ", cat_file, cat_racol, cat_deccol); for (i = 0; i < cat_numfilt; i++) { printf(" %s=%d ", cat_filter_array[i], cat_magcol_array[i]); } printf("\n"); #endif /* * now arguments pertaining to the detected star list file(s) */ strncpy(det_file[num_detected_files], argv[cur_arg_index++], MAX_FILENAME_LEN); det_file[num_detected_files][MAX_FILENAME_LEN] = '\0'; num_detected_files++; if (sscanf(argv[cur_arg_index++], "%d", &det_racol) != 1) { shFatal("invalid argument for column of RA values in detected file"); } if (sscanf(argv[cur_arg_index++], "%d", &det_deccol) != 1) { shFatal("invalid argument for column of Dec values in detected file"); } if (sscanf(argv[cur_arg_index++], "%d", &det_jdcol) != 1) { shFatal("invalid argument for column of JD values in detected file"); } if (sscanf(argv[cur_arg_index++], "%d", &det_aircol) != 1) { shFatal("invalid argument for column of airmass values in detected file"); } if (sscanf(argv[cur_arg_index++], "%d", &det_flagcol) != 1) { shFatal("invalid argument for column of flag values in detected file"); } if (sscanf(argv[cur_arg_index++], "%d", &det_numfilt) != 1) { shFatal("invalid argument for number of passbands in detected file"); } for (i = 0; i < det_numfilt; i++) { if (sscanf(argv[cur_arg_index], "%d,%d,%d", &filter_col, &det_magcol, &det_magerrcol) != 3) { shFatal("invalid argument %s for detected filt,mag,magerr", argv[cur_arg_index]); } filtercol_array[i] = filter_col; magcol_array[i] = det_magcol; magerrcol_array[i] = det_magerrcol; cur_arg_index++; } #ifdef DEBUG printf(" detfiles have RA %2d, Dec %2d, JD %d, airmass %d , flag %d", det_racol, det_deccol, det_jdcol, det_aircol, det_flagcol); for (i = 0; i < det_numfilt; i++) { printf(" %d,%d,%d ", filtercol_array[i], magcol_array[i], magerrcol_array[i]); } #endif /* * now it becomes a bit tricky: the remaining arguments start out * with names of detected star list files, but also may include * "radius=" or "outfile=". So we check each one carefully. */ for ( ; cur_arg_index < argc; cur_arg_index++) { if (strncmp(argv[cur_arg_index], "radius=", 7) == 0) { if (sscanf(argv[cur_arg_index] + 7, "%lf", &match_radius) != 1) { shFatal("invalid argument for radius argument"); } continue; } if (strncmp(argv[cur_arg_index], "outfile=", 8) == 0) { if (sscanf(argv[cur_arg_index] + 8, "%s", outfile) != 1) { shFatal("invalid argument for outfile argument"); } continue; } /* if we get here, interpret the arg as a detected-list file name */ strncpy(det_file[num_detected_files], argv[cur_arg_index], MAX_FILENAME_LEN); num_detected_files++; } #ifdef DEBUG printf(" there are %4d detected files: ", num_detected_files); for (i = 0; i < num_detected_files; i++) { printf(" %s ", det_file[i]); } printf("\n"); 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 */ #ifdef RADEC match_radius /= 3600.0; #endif /* read in the information from the catalog star file */ if (read_catalog_file(cat_file, cat_racol, cat_deccol, cat_numfilt, cat_filter_array, cat_magcol_array, &cat_nstar, &cat_array) != SH_SUCCESS) { shFatal("error reading raw star file %s", cat_file); } #ifdef DEBUG printf(" catalog star list %s contains %d stars\n", cat_file, cat_nstar); print_cat_array(cat_nstar, cat_array); #endif /* read the information from each "detected" star file, one at a time */ for (i = 0; i < num_detected_files; i++) { /* * currently using bogus fixed values of JD and airmass, * same for all files. */ if (read_detected_file(det_file[i], det_racol, det_deccol, det_jdcol, det_aircol, det_flagcol, det_numfilt, filtercol_array, magcol_array, magerrcol_array, &(det_nstar[i]), &(det_array[i])) != SH_SUCCESS) { shFatal("error reading detected star file %s", det_file[i]); } #ifdef DEBUG printf(" detected star list %s contains %d stars\n", det_file[i], det_nstar[i]); print_det_array(det_nstar[i], det_array[i]); #endif } /* * Okay, now we've finished reading in all the data. * * The next step is to go through all the detected-star files. * For each one, we * - look for a catalog match to each star * - if found, set the 'cat_match' field to point to * the corresponding catalog star */ if (find_matching_stars(cat_nstar, cat_array, num_detected_files, det_nstar, det_array, match_radius) != 0) { shFatal("error in find_matching_stars"); } /* * Mark some stars so that we don't include them in the photometric * solution, by setting the 'include' field to zero. */ if (mark_bad_stars(num_detected_files, det_nstar, det_array) != 0) { shFatal("error in mark_bad_stars"); } /* * initialize values in the SOLUTION structures, and allocate * space needed for all the photometric values we're going * to determine. */ initialize_solution(solution, num_detected_files); /* * we analyze the matched pairs, creating photometric solutions * for all possible passbands */ if (analyze_matched_pairs(cat_nstar, cat_array, num_detected_files, det_nstar, det_array, solution) != 0) { shFatal("error in analyze_matched_pairs"); } /* * Now, to print our results to several different sorts of output files */ /* first, files with matched pairs of standard mag and instrumental mag */ if (print_comps(outfile, cat_nstar, cat_array, num_detected_files, det_nstar, det_array, solution) != 0) { shError("print_comps fails"); return(1); } /* * since we've just calculated the RMS of solution, * we can now print the small file with photometric coefficients * and the RMS. */ if (print_coeffs(outfile, solution, num_detected_files) != 0) { shError("print_coeffs fails"); return(1); } /* * now, apply the photometric solution to all measurements of all * detected stars in all files -- not just the ones which had matching * entries in a catalog file. This step creates one ".cal" output * file for ALL the input detected-star files. */ if (calibrate_detected(outfile, cat_nstar, cat_array, num_detected_files, det_file, det_nstar, det_array, solution) != 0) { shError("calibrate_detected fails"); return(1); } /* free the arrrays we've created */ shFree(cat_array); for (i = 0; i < num_detected_files; i++) { shFree(det_array[i]); } free_solution(solution); return(0); } /**************************************************************************** * ROUTINE: find_matching_stars * * Find all catalog matches to detected stars. * * In more detail: * * - for each list of detected stars * for each detected star in the list * look for matching catalog star * * We consider position only as a criterion, and accept any * match within "radius" degrees. * * The strategy is * - sort the catalog list * - for each detected list * sort the detected list * for each star in detected list * define start, ending points in catalog list * check all catalog stars between start and end * if more than one catalog star matches, pick the one * which is closest to the detected position * if a match found, set detected star's "cat_match" field * * * RETURNS: * SH_SUCCESS if all goes well * SH_GENERIC_ERROR if not */ static int find_matching_stars ( int cat_nstar, /* I: number of stars in catalog array */ S_CAT *cat_array, /* I/O: array of catalog stars */ /* we will sort elements in array */ int num_detected_files, /* I: number of detected file lists */ int *det_nstar, /* I: array, number of stars in each */ /* detected array */ S_DET *det_array[], /* I/O: array, each is an array of */ /* detected stars. */ /* we will sort elements in each array */ double radius /* I: max dist for matching stars (degrees) */ ) { int i, num_match, total_match; int file_index; int start_index, closest_index, other_index, index; double start_ra; double dra, ddec, dist, other_dist; S_DET *det_star; S_CAT *cat_star; /* * now, sort the arrays by RA */ qsort((void *) cat_array, cat_nstar, sizeof(S_CAT), (PFI) compare_cat_by_ra); #ifdef DEBUG print_cat_array(cat_nstar, cat_array); #endif for (i = 0; i < num_detected_files; i++) { qsort((void *) det_array[i], det_nstar[i], sizeof(S_DET), (PFI) compare_det_by_ra); #ifdef DEBUG print_det_array(det_nstar[i], det_array[i]); #endif } /* * Okay, time to check every detected star for a catalog counterpart. * We work on one array of detected stars at a time. * * For each detected array, we * - pick next star in the array * - make binary search for the catalog star which has * an RA just less than (detected RA - radius) * - march forward through catalog array until (detected RA + radius)¸ * looking for best match */ total_match = 0; for (file_index = 0; file_index < num_detected_files; file_index++) { #ifdef DEBUG printf("working on detected file %d \n", file_index); #endif num_match = 0; for (i = 0; i < det_nstar[file_index]; i++) { det_star = &(det_array[file_index][i]); start_ra = det_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, cat_nstar, cat_array); closest_index = find_closest_index(det_star, cat_nstar, cat_array, start_index, radius); /* * 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 (det_star->ra < radius) { /* case A: need to check as if the star's RA were (360 - det_ra) */ start_ra = (360.0 - det_star->ra) - radius; start_index = find_first_index(start_ra, cat_nstar, cat_array); other_index = find_closest_index(det_star, cat_nstar, cat_array, start_index, radius); } else if (det_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, cat_nstar, cat_array); other_index = find_closest_index(det_star, cat_nstar, cat_array, 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. */ cat_star = &(cat_array[closest_index]); dist = calc_dist_between(cat_star, det_star, &dra, &ddec); cat_star = &(cat_array[other_index]); other_dist = calc_dist_between(cat_star, det_star, &dra, &ddec); if (dist < other_dist) { index = closest_index; } else { index = other_index; } } else { index = -1; } if (index == -1) { #ifdef DEBUG printf("FAIL no match found for star %4d %6.2f at (%10.5f, %10.5f)\n", i, det_star->mag[0], det_star->ra, det_star->dec); #endif } else { #ifdef DEBUG /* calculate distance in arcsec */ cat_star = &(cat_array[index]); dist = calc_dist_between(cat_star, det_star, &dra, &ddec); printf(" MATCH %4d %6.2f (%9.5f, %8.5f) = %4d (%9.5f, %8.5f) dist %5.1f\n", i, det_star->mag[0], det_star->ra, det_star->dec, index, cat_star->ra, cat_star->dec, dist); #endif det_star->cat_match_id = index; num_match++; } } #ifdef DEBUG printf(" in file %d, there are %d matches \n", file_index, num_match); #endif total_match += num_match; } #ifdef DEBUG printf(" grand total of %d matches \n", total_match); #endif return(SH_SUCCESS); } /**************************************************************************** * ROUTINE: compare_cat_by_ra * * Compare two S_CAT 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_cat_by_ra ( S_CAT *star1, /* I: compare "ra" field of this star .. */ S_CAT *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: 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); } /**************************************************************************** * ROUTINE: find_first_index * * Given an array of "cat_nstar" catalog star structures, which have already * been sorted in ascending order by RA, and given some value "ra", * return the index of the catalog 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 cat_nstar, /* I: number of stars in array */ S_CAT *cat_array /* I: array of catalog stars, sorted by RA */ ) { int top, bottom, mid; int retval; shAssert(cat_array != NULL); #ifdef DEBUG printf("find_first_index: looking for RA = %.4f\n", start_ra); #endif top = 0; if ((bottom = cat_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, cat_array[top].ra, mid, cat_array[mid].ra, bottom, cat_array[bottom].ra); #endif if (cat_array[mid].ra > start_ra) { bottom = mid; } else { top = mid; } } #ifdef DEBUG printf(" array[%4d] RA=%.4f array[%4d] RA=%.4f\n", top, cat_array[top].ra, bottom, cat_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 (cat_array[bottom].ra > start_ra) { #ifdef DEBUG printf(" choosing array[%4d] RA=%.4f \n", top, cat_array[top].ra); #endif retval = top; } else { #ifdef DEBUG printf(" choosing array[%4d] RA=%.4f \n", bottom, cat_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) && (cat_array[retval].ra >= start_ra)) { retval--; #ifdef DEBUG printf(" backing array[%4d] RA=%.4f \n", retval, cat_array[retval].ra); #endif } #ifdef DEBUG printf(" return array[%4d] RA=%.4f \n", retval, cat_array[retval].ra); #endif return(retval); } /**************************************************************************** * ROUTINE: find_closest_index * * Given a detected S_DET structure, * an array of "cat_nstar" catalog star structures, which have already * been sorted in ascending order by RA, * an index into the "cat_array" at which to start looking, * and a matching radius (in degrees), * try to find the catalog star 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_DET *det_star, /* I: want to find a match to this star */ int cat_nstar, /* I: there are this many catalog stars */ /* in cat_array */ S_CAT *cat_array, /* I: array of catalog stars, sorted by RA */ int start_index, /* I: start looking here in catalog array */ double radius /* I: max distance (degrees) for a match */ ) { int i, closest; double cat_ra, cat_dec, det_ra, det_dec; #ifdef RADEC double dec_factor; #endif double max_ra, min_dist2; double dra, ddec, dist2; S_CAT *cat_star; shAssert(det_star != NULL); shAssert(cat_array != NULL); min_dist2 = radius*radius; det_ra = det_star->ra; det_dec = det_star->dec; #ifdef RADEC dec_factor = cos(det_dec*DEGTORAD); if (dec_factor <= 0.0) { shError("closest_index: can't handle Dec = 90 degrees"); return(-1); } max_ra = det_ra + (radius/dec_factor); #else max_ra = det_ra + radius; #endif closest = -1; for (i = start_index; i < cat_nstar; i++) { cat_star = &(cat_array[i]); cat_ra = cat_star->ra; /* check to see if we can stop looking */ if (cat_ra > max_ra) { break; } cat_dec = cat_star->dec; #ifdef RADEC dra = (cat_ra - det_ra)*dec_factor; ddec = (cat_dec - det_dec); #else dra = (cat_ra - det_ra); ddec = (cat_dec - det_dec); #endif dist2 = (dra*dra + ddec*ddec); if (dist2 < min_dist2) { closest = i; min_dist2 = dist2; } } return(closest); } /**************************************************************************** * ROUTINE: mark_bad_stars * * Given 'num_detected_files' arrays of detected stars, * (the number of stars in each is given by corresponding * element of the "det_nstar" array) * * do the following: * * walk through the lists of all detected stars, and set the * "include" field to zero for all stars we don't * want to include in the photometric solution. * * At the moment, this means stars with non-BAD magnitudes, but * with negative values in the "magerr" fields (that means * the CCD may have been saturated). We'll ignore all stars with * "BAD_MAG" magnitudes later on, but this step allows us to * apply other criteria, too. * * * RETURN: * 0 if all goes well * 1 if there's an error */ static int mark_bad_stars ( int num_detected_files, /* I: number of detected file arrays */ int *det_nstar, /* I: array, number of stars in each of */ /* the elements of "det_array" */ S_DET *det_array[] /* I/O: array, each element of which is an */ /* array of "detected" stars */ ) { int i, j, k; int num_marked_stars; S_DET *det_star; #ifdef DEBUG printf(" entering mark_bad_stars ... \n"); #endif num_marked_stars = 0; for (i = 0; i < num_detected_files; i++) { for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); for (k = 0; k < MAX_PASSBANDS; k++) { if ((det_star->mag[k] != BAD_MAG) && (det_star->magerr[k] < 0)) { det_star->include = 0; num_marked_stars++; } } } } #ifdef DEBUG printf(" mark_bad_stars: total number of marked stars is %d \n", num_marked_stars); #endif return(0); } /**************************************************************************** * ROUTINE: analyze_matched_pairs * * Given 'num_detected_files' arrays of detected stars, * (the number of stars in each is given by corresponding * element of the "det_nstar" array) * an array of "cat_nstar" catalog stars, * an output FILE stream, * * do the following: * * 1. verify that the primary catalog passband matches the * the primary detected passband (for ALL stars). * If not, print error message and quit * * 2. set up a photometric solution of the form (for example): * * V = v + a + b(v-i) + kX * j * * where V is the primary catalog magnitude * v is the primary detected instr. mag * i is the secondary detected instr. mag * X is the airmass from detected data file * * a is unknown zero-point term for frame-pair "j" * j * b is unknown color term (same for all frames) * k is unknown extinction term (same for all) * * 3. solve for the unknown coefficients a,b,k by method of least * squares. For each detected star, we calculate a * weight based on the uncertainty in its primary * instrumental magnitude: * * [ 1 ] p * weight = [ ------------ ] where p = 0.5 * [ uncertainty ] * * but we place bounds on this weight: * * if uncertainty < MIN_UNCERT use MIN_UNCERT * uncertainty > MAX_UNCERT set weight to TINY_WEIGHT * * Then, starting with the photometric solution from 2 above, * we can write a big equation including all the stars like so: * * ( w1 (v1-i1)w1 X1w1 ) ( w1(V1­v1) ) * ( w2 (v2-i2)w2 X2w2 ) (a0 ) ( w2(V2­v2) ) * ( w3 (v3-i3)w3 X3w3 ) (a1 ) ( w3(V3­v3) ) * ( . ) ... = ( . ) * ( . ) (aN ) ( . ) * ( . ) ( b ) ( . ) * ( . ) ( k ) ( . ) * ( wN (vN-iN)wN XNwN ) ( wN(VN­vN) ) * * * or, in compact form * * A * x = b * * where the unknown values aj,b,k are the elements of vector "x". * Our job is to calculate the known values of the matrix "A" * and vector "b", then call a routine which will solve the * big matrix equation and give us the values of vector "x". * * 4. after having solved the matrix equation, print out the * values aj,b,k of the photometric solution, together with * their estimated uncertainty. * * * Note that we do NOT include any magnitudes with value "BAD_MAG" in * the calculations -- they denote missing or invalid data. * * RETURN: * 0 if all goes well * 1 if there's an error */ static int analyze_matched_pairs ( int cat_nstar, /* I: number of stars in "cat" array */ S_CAT *cat_array, /* I: array of all "cat" star measurements */ int num_detected_files, /* I: number of detected file arrays */ int *det_nstar, /* I: array, number of stars in each of */ /* the elements of "det_array" */ S_DET *det_array[], /* I: array, each element of which is an */ /* array of "detected" stars */ S_SOLUTION solution[] /* O: place results of solution in here */ ) { int i, j, k; int num_total_stars; int num_tiny_weight; int num_unknowns; /* number of unknowns in photom equations */ int retval; int row; int pri_index, sec_index, cat_index; double *matrix_a, *vector_b, *vector_x, *vector_s; double weight, color; S_CAT *cat_star; S_DET *det_star; #ifdef DEBUG printf(" entering analyze_matched_pairs ... \n"); #endif /* * We loop over all possible detected passbands, checking to see if * we can calculate a solution in each one. If so, we do it. */ for (pri_index = 0; pri_index < MAX_PASSBANDS; pri_index++) { if (find_passband_indices(cat_nstar, cat_array, num_detected_files, det_nstar, det_array, pri_index, &sec_index, &cat_index) != 0) { #ifdef DEBUG printf("analyze_matched_pairs: pri_index %d can't be reduced\n", pri_index); #endif continue; } /* * set up the big matrix equation. First step is to allocate space * for the needed matrix and vectors ... */ num_total_stars = 0; for (i = 0; i < num_detected_files; i++) { for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); if (put_star_in_solution(det_star, cat_array, pri_index, sec_index, cat_index) == 0) { continue; } num_total_stars++; } } #ifdef DEBUG printf(" total number of detected stars is %d \n", num_total_stars); #endif /* * how many unknown values are there? we solve for one zero-point * per frame, plus a single color term, and perhaps a single * extinction coefficient. */ num_unknowns = (num_detected_files - 1) + PHOTOM_COEFFS; /* the matrix of known values */ matrix_a = (double *) shMalloc(num_total_stars*num_unknowns *sizeof(double)); /* the vector of known values on right-hand side of equation */ vector_b = (double *) shMalloc(num_total_stars*sizeof(double)); /* the vector of unknowns: the photometric coefficients */ vector_x = (double *) shMalloc(num_unknowns*sizeof(double)); /* a vector for uncertainty in values of unknowns */ vector_s = (double *) shMalloc(num_unknowns*sizeof(double)); /* * Next step: fill the elements of the matrix. The matrix has * * one row per detected star with matching catalog star * num_unknown columns, for coefficients a0, a1, ..., aN, b, * and possibly k * * the values of elements in each row look like this: * * cols 1-N: weight for star (if current image is number N) * or * zero, otherwise * * col N+1: weight times (instr. primary mag - instr. seconary mag) * * col N+2: weight times airmass * if we are solving for extinction */ row = 0; for (i = 0; i < num_detected_files; i++) { for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); if (put_star_in_solution(det_star, cat_array, pri_index, sec_index, cat_index) == 0) { continue; } weight = assign_weight(det_star, pri_index, sec_index); color = det_star->mag[pri_index] - det_star->mag[sec_index]; /* * the entry for cols 1-N is zero except for the column * corresponding to the current image; the variable "i" * in loop denotes this image */ for (k = 0; k < num_detected_files; k++) { if (k != i) { matrix_a[row*num_unknowns + k] = 0.0; } else { matrix_a[row*num_unknowns + k] = weight; } } /* we always have one column for color term */ matrix_a[row*num_unknowns + num_detected_files] = weight*color; /* if we're solving for extinction, one more column for it */ #if PHOTOM_COEFFS == 3 matrix_a[row*num_unknowns + num_detected_files + 1] = weight*det_star->airmass; #endif row++; } } if (row != num_total_stars) { shFatal("analyze_matched_pairs: mismatch A: num_total_stars %d row %d\n", num_total_stars, row); } #ifdef DEBUG { int i, j; printf("matrix_a follows ... \n"); for (i = 0; i < num_total_stars; i++) { for (j = 0; j < num_unknowns; j++) { printf(" %10.4f ", matrix_a[i*num_unknowns+j]); } printf("\n"); } printf("that was matrix_a \n"); } #endif /* * now fill the elements of "vector_b", which are based on the * measured magnitudes. */ row = 0; num_tiny_weight = 0; for (i = 0; i < num_detected_files; i++) { for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); if (put_star_in_solution(det_star, cat_array, pri_index, sec_index, cat_index) == 0) { continue; } weight = assign_weight(det_star, pri_index, sec_index); if (weight == PHOTOM_TINY_WEIGHT) { num_tiny_weight++; } cat_star = &(cat_array[det_star->cat_match_id]); shAssert(cat_star != NULL); vector_b[row] = weight* (cat_star->mag[cat_index] - det_star->mag[pri_index]); row++; } } #ifdef DEBUG printf("total stars %6d num_tiny_weight %6d num good %6d\n", num_total_stars, num_tiny_weight, num_total_stars - num_tiny_weight); #endif if (row != num_total_stars) { shFatal("analyze_matched_pairs: mismatch B: num_total_stars %d row %d\n", num_total_stars, row); } /* * Okay, we're ready to call the matrix routine to perform a least-squares * solution for our big photometric equation. The routine will * fill the elements of two arrays: * * vector_x element 0 contains "a0", zero-point for image 0 * 1 contains "a1", zero-point for image 1 * ... * N-1 contains zero-point for image N-1 * N "b", color coeff * N+1 (may be) "k", extinction coeff * * vector_s element 0 contains uncertainty in "a0" * 1 uncertainty in "a1" * ... * N-1 uncertainty in "a(N-1)" * N uncertainty in "b" * N+1 (may be) uncertainty in "k" */ retval = solve_matrix_equation(num_total_stars, num_unknowns, matrix_a, vector_b, vector_x, vector_s); if (retval != 0) { shError("analyze_matched_pairs: solve_matrix_equation fails"); return(1); } /* * hurray! We have a decent photometric solution. Now, we need * to tell the user about it. We place the values of the coefficients * into the "solution" structure, together with some auxiliary * information. Note that we do not (yet) set the RMS field. */ solution[pri_index].num_obs = num_total_stars - num_tiny_weight; for (i = 0; i < num_detected_files; i++) { solution[pri_index].a[i] = vector_x[i]; solution[pri_index].asig[i] = vector_s[i]; } solution[pri_index].b = vector_x[num_detected_files]; solution[pri_index].bsig = vector_s[num_detected_files]; #if PHOTOM_COEFFS >= 3 solution[pri_index].k = vector_x[num_detected_files + 1]; solution[pri_index].ksig = vector_s[num_detected_files + 1]; #endif /* de-allocate memory for all the vectors and matrices we used */ shFree(matrix_a); shFree(vector_b); shFree(vector_x); shFree(vector_s); } /* end of loop over all possible detected star primary passbands */ return(0); } /************************************************************************ * * ROUTINE: sdetFill * * DESCRIPTION: * Fill the given S_DET structure with the given information. * * RETURN: * nothing */ void sdetFill ( S_DET *star, /* I: pointer to struct to fill with info */ double ra, /* I: RA value for new star (decimal degrees) */ double dec, /* I: Dec value for new star */ double jd, /* I: Julian Date of observation */ double airmass, /* I: airmass of observation */ unsigned int quality, /* I: "quality flag" -- for all filters */ int num_filt, /* I: number of filters with data */ char *filt_array[], /* I: name of filter(s) in which star was observed */ double *mag, /* I: instrumental magnitude(s) in each filter */ double *magerr /* I: estimated uncertainty(s) in instrumental mag */ /* in each filter */ ) { static int id_number = 0; int i; shAssert(star != NULL); star->id = id_number++; star->ra = ra; star->dec = dec; star->jd = jd; star->airmass = airmass; for (i = 0; i < num_filt; i++) { 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[i]; star->magerr[i] = magerr[i]; star->flag[i] = ((quality & (0xFF << i*8)) >> i*8); } star->cat_match_id = -1; star->include = 1; } /************************************************************************ * * ROUTINE: sdetPrint * * DESCRIPTION: * Print to stdout an SDET structure. * * RETURN: * nothing */ void sdetPrint ( S_DET *star /* I: star to print out */ ) { int i; shAssert(star != NULL); printf("sdet %5d %10.5f %10.5f ", star->id, star->ra, star->dec); printf(" %13.4f %5.2f ", star->jd, star->airmass); for (i = 0; i < MAX_PASSBANDS; i++) { if (star->filter[i][0] != '\0') { printf(" %2s %6.3f %5.3f ", star->filter[i], star->mag[i], star->magerr[i]); } } printf(" %6d %3d\n", (int) star->cat_match_id, star->include); } /************************************************************************ * * ROUTINE: print_det_array * * DESCRIPTION: * Print to stdout the entire contents of an array of detected stars. * This is probably used only in debugging. * * RETURN: * nothing */ static void print_det_array ( int num_star, /* I: number of stars in the array */ S_DET *array /* I: array of S_DETs to print */ ) { int i; for (i = 0; i < num_star; i++) { sdetPrint(&(array[i])); } } /************************************************************************ * * 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])); } } /************************************************************************** * ROUTINE: calc_dist_between * * DESCRIPTION: * Given two stars, calculate the distance between them. * Place the distance in each coordinate into the final two args. * * RETURNS: * Distance between stars */ static double calc_dist_between ( S_CAT *cat_star, /* I: calc distance between this catalog star ... */ S_DET *det_star, /* I: ... and this detected star */ double *ra_dist, /* O: distance between them in the RA (x) dir */ double *dec_dist /* O: distance between them in the Dec (y) dir */ ) { double dist; #ifdef RADEC double dra, ddec; /* calculate distance in arcsec between (RA, Dec) positions */ dra = (cat_star->ra - det_star->ra)*cos(cat_star->dec*DEGTORAD); ddec = (cat_star->dec - det_star->dec); dist = sqrt(dra*dra + ddec*ddec)*3600.0; *ra_dist = dra; *dec_dist = ddec; #else double dx, dy; /* calculate distance in arbitrary units between (x, y) positions */ dx = (cat_star->ra - det_star->ra); dy = (cat_star->dec - det_star->dec); dist = sqrt(dx*dx + dy*dy); *ra_dist = dx; *dec_dist = dy; #endif return(dist); } /************************************************************************** * ROUTINE: assign_weight * * DESCRIPTION: * Based on the given uncertainty in an instrumental magnitude, * assign a weight for the star in the photometric solution. * We watch out for "bad" values, assigning a tiny weight for such stars. * * RETURNS: * weight for star in solution */ static double assign_weight ( S_DET *det_star, /* I: want to calc weight for this star */ int pri_index, /* I: primary magnitude index */ int sec_index /* I: secondary magnitude index */ ) { double weight, uncert; double power = -0.5; /* power to which we raise uncert to calc weight */ /* sanity check */ shAssert(det_star != NULL); /* check that all the instrumental mags and uncertainties are valid */ if ((det_star->mag[pri_index] == BAD_MAG) || (det_star->mag[sec_index] == BAD_MAG)) { weight = PHOTOM_TINY_WEIGHT; return(weight); } if ((det_star->magerr[pri_index] == BAD_MAG) || (det_star->magerr[sec_index] == BAD_MAG)) { weight = PHOTOM_TINY_WEIGHT; return(weight); } /* the mags and magerrs are valid, so use 'em to calculate a weight */ uncert = det_star->magerr[pri_index]; if (uncert <= PHOTOM_MIN_UNCERT) { weight = pow(PHOTOM_MIN_UNCERT, power); } else if (uncert >= PHOTOM_MAX_UNCERT) { weight = PHOTOM_TINY_WEIGHT; } else { weight = pow(uncert, power); } return(weight); } /*************************************************************************** * ROUTINE: print_coeffs * * DESCRIPTION: print the values of the coefficients in the photometric * solution, together with their uncertainties, * to a file with the given "basename", and extension * COEFFS_EXT. * * We check the value of PHOTOM_COEFFS; if it is <= 2, then * we didn't solve for extinction. * * RETURNS: * 0 if all goes well * 1 if there's an error */ static int print_coeffs ( char *basename, /* I: base name of output file */ S_SOLUTION solution[], /* I: holds values of photometric coeffs */ int num_detected_files /* I: number of images included in solution */ ) { char fullname[MAX_FILENAME_LEN+1]; char pri_passband[FILTER_LENGTH + 1]; char sec_passband[FILTER_LENGTH + 1]; char cat_passband[FILTER_LENGTH + 1]; int i, pri_index, sec_index, cat_index; int num_zeropts, zeropt_index; FILE *fp; /* sanity checks */ shAssert(basename != NULL); shAssert(solution != NULL); if ((strlen(basename) + strlen(COEFFS_EXT)) >= MAX_FILENAME_LEN) { shError("print_coeffs: basename %s is too long", basename); return(1); } sprintf(fullname, "%s%s", basename, COEFFS_EXT); if ((fp = fopen(fullname, "w")) == NULL) { shError("print_coeffs: can't open file %s", fullname); return(1); } num_zeropts = num_detected_files; #ifdef DEBUG printf("num_zeropts is %d \n", num_zeropts); #endif /* * loop over all possible detected passbands, and print a line * with coefficients for each one with a photometric solution. */ for (i = 0; i < MAX_PASSBANDS; i++) { get_passband_values(i, &pri_index, &sec_index, &cat_index, pri_passband, sec_passband, cat_passband); if (valid_indices(pri_index, sec_index, cat_index) == 0) { #ifdef DEBUG printf("print_coeffs: no solution in index %d, skipping \n", pri_index); #endif continue; } for (zeropt_index = 0; zeropt_index < num_zeropts; zeropt_index++) { fprintf(fp, "%s=%s,(%s-%s) N %6d a %2d %7.3f %6.3f b %7.3f %6.3f ", cat_passband, pri_passband, pri_passband, sec_passband, solution[pri_index].num_obs, zeropt_index, solution[pri_index].a[zeropt_index], solution[pri_index].asig[zeropt_index], solution[pri_index].b, solution[pri_index].bsig); #if PHOTOM_COEFFS >= 3 fprintf(fp, " k %7.3f %6.3f ", solution[pri_index].k, solution[pri_index].ksig); #endif fprintf(fp, " RMS %6.3f ", solution[pri_index].rms); fprintf(fp, "\n"); } } fclose(fp); return(0); } /***************************************************************************** * PROCEDURE: print_comps * * DESCRIPTION: We have already matched up the detected stars against * the catalog stars. What we do here is to print out one line per * observation, listing the catalog magnitudes, the instrumental magnitudes, * and the differences. This allows us to make diagnostic plots later on. * * We also calculate the RMS of the photometric solution, since * we have all the information we need. * RETURNS: * 0 if all goes well * 1 if an error occurs */ static int print_comps ( char *basename, /* I: base name of output file */ int cat_nstar, /* I: number of stars in "cat" array */ S_CAT *cat_array, /* I: array of all "cat" star measurements */ int num_detected_files, /* I: number of detected file arrays */ int *det_nstar, /* I: array, number of stars in each of */ /* the elements of "det_array" */ S_DET *det_array[], /* I: array, each element of which is an */ /* array of "detected" stars */ S_SOLUTION solution[] /* O: place results of solution in here */ ) { char fullname[MAX_FILENAME_LEN+1]; char pri_passband[FILTER_LENGTH + 1]; char sec_passband[FILTER_LENGTH + 1]; char cat_passband[FILTER_LENGTH + 1]; int i, j, cat_id, num_pairs; int index, pri_index, sec_index, cat_index; double raw_diff, color, predict, predict_err, predict_diff; double sumsq, rms; double weight; FILE *fp; S_DET *det_star; S_CAT *cat_star; /* sanity checks */ shAssert(cat_array != NULL); shAssert(det_nstar != NULL); shAssert(det_array != NULL); shAssert(solution != NULL); /* loop over all possible detected passbands */ for (index = 0; index < MAX_PASSBANDS; index++) { get_passband_values(index, &pri_index, &sec_index, &cat_index, pri_passband, sec_passband, cat_passband); if (valid_indices(pri_index, sec_index, cat_index) == 0) { continue; } /* open the output file */ shAssert(basename != NULL); if ((strlen(basename) + 1 + strlen(pri_passband) + strlen(COMPARE_EXT)) >= MAX_FILENAME_LEN) { shError("print_comps: basename %s is too long", basename); return(1); } #if 1 sprintf(fullname, "%s_%s%s", basename, pri_passband, COMPARE_EXT); #else sprintf(fullname, "%s%s", basename, COMPARE_EXT); #endif if ((fp = fopen(fullname, "w")) == NULL) { shError("print_comps: can't open file %s", fullname); return(1); } /* * we walk through the array of detected files */ sumsq = 0; num_pairs = 0; for (i = 0; i < num_detected_files; i++) { /* * For each detected file, we walk through all observations. * If an observation has a matching catalog star, * we calculate some differences * print them to the output file * accumulate differences for RMS calculation * * After we process all the files, we set the rms field in the * "solution" structure. */ for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); shAssert(det_star != NULL); if (det_star->cat_match_id == -1) { continue; } if (det_star->include == 0) { continue; } cat_id = det_star->cat_match_id; cat_star = &(cat_array[cat_id]); shAssert(cat_star != NULL); if (put_star_in_solution(det_star, cat_array, pri_index, sec_index, cat_index) != 1) { printf("put_star_in_solution says no\n"); } else { printf("put_star_in_solution says yes\n"); } /* calculate difference in raw primary mag, (cat - det) */ /* but ignore stars with bad detected or reference mags */ if ((cat_star->mag[cat_index] == BAD_MAG) || (det_star->mag[pri_index] == BAD_MAG)) { continue; } raw_diff = cat_star->mag[cat_index] - det_star->mag[pri_index]; /* calculate color of star, using instrumental mags */ if (det_star->mag[sec_index] == BAD_MAG) { shError("print_comps: file %d star %d has bad primary color?, i, j"); continue; } color = det_star->mag[pri_index] - det_star->mag[sec_index]; /* skip any star which has PHOTOM_TINY_WEIGHT */ weight = assign_weight(det_star, pri_index, sec_index); if (weight == PHOTOM_TINY_WEIGHT) { continue; } /* * calculate the predicted, calibrated magnitude for this star, * based on its instrumental mags and the photometric solution. */ apply_solution(det_star, solution, i, pri_index, sec_index, &predict, &predict_err); shAssert(predict != BAD_MAG); predict_diff = cat_star->mag[cat_index] - predict; /* now print out a line with lots of info on this observation */ fprintf(fp, "%4d %6d %9.4f %9.4f %6.3f ", i, det_star->id, det_star->ra, det_star->dec, det_star->airmass); fprintf(fp, " %6.3f %6.3f %6.3f ", det_star->mag[pri_index], det_star->mag[sec_index], color); fprintf(fp, " %6.3f ", cat_star->mag[cat_index]); fprintf(fp, " %6.3f ", raw_diff); fprintf(fp, " %6.3f %6.3f ", predict, predict_err); fprintf(fp, " %6.3f ", predict_diff); fprintf(fp, "\n"); /* * accumulate square of differences between catalog mag * and predicted mag, based on photometric solution. * We use this to calculate RMS of solution at end. */ sumsq += predict_diff*predict_diff; num_pairs++; } } #ifdef DEBUG printf("print_comps: num_pairs %6d \n", num_pairs); #endif /* set RMS value here */ if (num_pairs > 0) { rms = sqrt(sumsq/num_pairs); } else { rms = 0.00; } solution[pri_index].rms = rms; fclose(fp); } return(0); } /************************************************************************** * PROCEDURE: apply_solution * * DESCRIPTION: Given a pointer to an array of S_SOLUTION structures, * and a pointer to an S_DET structure containing raw instrumental magnitudes, * apply the photometric solution to the star and derive the * calibrated magnitude. * * We use the 'pri_index' argument to look up the proper element of * the 'solution_array' for the star in question. * * Place the calibrated magnitude into the final argument. * * If the any of the required instrumental magnitudes for the given star * is invalid (marked as BAD_MAG), then set the output calibrated * magnitude to BAD_MAG also. * * RETURNS: * nothing * */ static void apply_solution ( S_DET *det_star, /* I: star with instrumental mags */ S_SOLUTION solution_array[], /* I: coeffs of photometric solution */ int image_index, /* I: index of image */ int pri_index, /* I: index of primary det. passband */ int sec_index, /* I: index of secondary det. passband */ double *calib_mag, /* O: place calibrated magnitude here */ double *calib_magerr /* O: uncertainty in calibrated mag here */ ) { int zeropt_index; double color, color_err, mag, magerr; double contrib_1, contrib_2, contrib_3; double s_pri, s_sec; double stuff; double contrib_4; S_SOLUTION *solution; shAssert(det_star != NULL); shAssert(solution_array != NULL); shAssert((pri_index >= 0) && (pri_index < MAX_PASSBANDS)); shAssert((sec_index >= 0) && (sec_index < MAX_PASSBANDS)); solution = &(solution_array[pri_index]); zeropt_index = image_index; /* * check the input raw magnitudes -- if any are invalid, * then we set the output calibrated magnitude to BAD_MAG * (but don't return with an error). */ if ((det_star->mag[pri_index] == BAD_MAG) || (det_star->mag[sec_index] == BAD_MAG)) { #ifdef DEBUG shError("apply_solution: star ID %d has bad primary mags", det_star->id); #endif *calib_mag = BAD_MAG; return; } /* * calculate the predicted catalog mag, * based on photometric solution */ color = det_star->mag[pri_index] - det_star->mag[sec_index]; mag = det_star->mag[pri_index] + solution->a[zeropt_index] + solution->b*color; #if PHOTOM_COEFFS >= 3 mag += solution->k*det_star->airmass; #endif /* * We show below two ways to calculate the uncertainty in the * final, calibrated magnitude. The first includes the uncertainties * in all photometric coefficients: it yields an honest-to-goodness * uncertainty in the calibrated magnitude, period. The second * ignores some of the systematic uncertainty. */ #if IGNORE_SYSTEMATICS == 0 /* * calculate the uncertainty in the calibrated magnitude, via the * following procedure: * * [ 2 2 2 2 2 2 ] * sigma(cal mag) = sqrt [ s_pri + s_a + b *s_color + X *s_k ] * [ ] * * where * s_pri is uncertainty in raw instrumental primary mag * s_a is uncertainty in zero-point coeff * b is color-term coeff * s_color is uncertainty in raw instrumental color * (see below) * X is airmass * s_k is uncertainty in airmass coeff * * The uncertainty in instrumental color is simply * * 2 2 * sqrt ( s_pri + s_sec ) * * where * s_pri is uncertainty in raw instrumental primary mag * s_sec is uncertainty in raw instrumental secondary mag * * If we are unable to calculate this uncertainty -- because the * uncertainty in some instrumental mag is BAD_MAG -- then set * the output uncertainty to BAD_MAG. */ s_pri = det_star->magerr[pri_index]; s_sec = det_star->magerr[sec_index]; if ((s_pri == BAD_MAG) || (s_sec == BAD_MAG)) { magerr = BAD_MAG; } else { color_err = sqrt(s_pri*s_pri + s_sec*s_sec); contrib_1 = s_pri*s_pri; contrib_2 = solution->asig*solution->asig[zeropt_index]; contrib_3 = solution->b*solution->b*color_err*color_err; stuff = contrib_1 + contrib_2 + contrib_3; #ifdef PHOTOM_COEFFS >= 3 contrib_4 = (solution->ksig*solution->ksig)* (det_star->airmass*det_star->airmass); #else contrib_4 = 0; #endif stuff += contrib_4; magerr = sqrt(stuff); } #else /* * calculate the uncertainty in the calibrated magnitude, via the * following procedure: * * [ 2 2 2 ] * sigma(cal mag) = sqrt [ s_pri + b *s_color ] * [ ] * * where * s_pri is uncertainty in raw instrumental primary mag * b is color-term coeff * s_color is uncertainty in raw instrumental color * (see below) * * The uncertainty in instrumental color is simply * * 2 2 * sqrt ( s_pri + s_sec ) * * where * s_pri is uncertainty in raw instrumental primary mag * s_sec is uncertainty in raw instrumental secondary mag * * If we are unable to calculate this uncertainty -- because the * uncertainty in some instrumental mag is BAD_MAG -- then set * the output uncertainty to BAD_MAG. */ s_pri = det_star->magerr[pri_index]; s_sec = det_star->magerr[sec_index]; if ((s_pri == BAD_MAG) || (s_sec == BAD_MAG)) { magerr = BAD_MAG; } else { color_err = sqrt(s_pri*s_pri + s_sec*s_sec); contrib_1 = s_pri*s_pri; contrib_2 = 0; contrib_3 = solution->b*solution->b*color_err*color_err; contrib_4 = 0; stuff = contrib_1 + contrib_2 + contrib_3 + contrib_4; magerr = sqrt(stuff); } #endif *calib_mag = mag; *calib_magerr = magerr; return; } /***************************************************************************** * PROCEDURE: calibrate_detected * * DESCRIPTION: We have already calculated a photometric solution for * the night. We now walk through all the input files of detected * stars with instrumental magnitudes. We create one giant * output file, with the given output file name base, and * with a suffix CALIBRATED_EXT. For each star in each input file, * we create a line in the output file. The line has a format like this: * * ID RA Dec JD filt mag magerr flag [filt mag magerr flag ...] * * where * ID is some internal ID number for the star, * not a real catalog ID of any sort * RA is the J2000 Right Ascension (decimal degrees) * Dec is the J2000 Declination (decimal degrees) * JD is the Julian Date of midpoint of exposure * * The next set of columns must appear in groups of 4 * * filt name of a passband * mag calibrated magnitude in that passband * magerr uncertainty in magnitude in that passband * flag an integer, the OR of several bits: * 0: star was okay in this passband * 1: star was probably saturated in this passband * 2: star was close to a bad region * 4: star was close to an edge of image * * A value of BAD_MAG indicates invalid or unknown data. * * RETURNS: * 0 if all goes well * 1 if an error occurs */ static int calibrate_detected ( char *basename, /* I: stem of output file name */ int cat_nstar, /* I: number of stars in "cat" array */ S_CAT *cat_array, /* I: array of all "cat" star measurements */ int num_detected_files, /* I: number of detected file arrays */ char det_file_name[MAX_DETECTED_FILES][MAX_FILENAME_LEN+1], /* I: names of the input detected-star files */ int *det_nstar, /* I: array, number of stars in each of */ /* the elements of "det_array" */ S_DET *det_array[], /* I: array, each element of which is an */ /* array of "detected" stars */ S_SOLUTION solution[] /* I: the photometric solution coefficients */ ) { char fullname[MAX_FILENAME_LEN+1]; char pri_passband[FILTER_LENGTH+1]; char sec_passband[FILTER_LENGTH+1]; char cat_passband[FILTER_LENGTH+1]; int i, j; int index, pri_index, sec_index, cat_index; double predict, predict_err; FILE *fp; S_DET *det_star; /* sanity checks */ shAssert(cat_array != NULL); shAssert(det_nstar != NULL); shAssert(det_array != NULL); shAssert(solution != NULL); /* * create the output file, and open for writing */ sprintf(fullname, "%s%s", basename, CALIBRATED_EXT); if ((fp = fopen(fullname, "w")) == NULL) { shError("calibrate_detected: can't open output file %s", fullname); return(1); } /* * we walk through the array of detected files */ for (i = 0; i < num_detected_files; i++) { /* * For each detected file, we walk through all observations. * We try to calculate a calibrated magnitude for each * detected star. We place the calibrated * magnitude into the appropriate 'calib_mag' field, * and the estimated uncertainty in the 'calib_magerr' field. * We may also set the 'flag' field, if necessary. */ for (j = 0; j < det_nstar[i]; j++) { det_star = &(det_array[i][j]); shAssert(det_star != NULL); /* loop over passbands to do the calibration for this star ... */ for (index = 0; index < MAX_PASSBANDS; index++) { get_passband_values(index, &pri_index, &sec_index, &cat_index, pri_passband, sec_passband, cat_passband); if (valid_indices(pri_index, sec_index, cat_index) == 0) { continue; } apply_solution(det_star, solution, i, pri_index, sec_index, &predict, &predict_err); det_star->calib_mag[pri_index] = predict; det_star->calib_magerr[pri_index] = predict_err; } /* end of loop over passbands */ /* now print out a line with lots of info on this observation */ fprintf(fp, " %6d %9.4f %9.4f %13.5f ", det_star->id, det_star->ra, det_star->dec, det_star->jd); /* * loop over passbands again to print information. * We separate these two loops so that have more flexibility * in the printing out -- include all stars, or just those * with valid calibrated values. */ for (index = 0; index < MAX_PASSBANDS; index++) { get_passband_values(index, &pri_index, &sec_index, &cat_index, pri_passband, sec_passband, cat_passband); if (valid_indices(pri_index, sec_index, cat_index) == 0) { continue; } fprintf(fp, " %3s %6.3f %6.3f %-6u ", det_star->filter[pri_index], det_star->calib_mag[pri_index], det_star->calib_magerr[pri_index ], det_star->flag[pri_index]); } fprintf(fp, "\n"); } /* end of loop over all stars detected in this file */ } fclose(fp); return(0); } #if 0 /************************************************************************** * PROCEDURE: set_calib_name * * DESCRIPTION: given the name of an input, detected-star file, * create an associated output file name. Do so by replacing * any existing file extension with CALIBRATED_EXT. * * RETURNS: * 0 if all goes well * 1 if an error occurs */ static int set_calib_name ( char *input_name, /* I: use this as a base for new name */ char *output_name /* O: write new file name here */ ) { int i, len, pos, found; shAssert(input_name != NULL); /* * search for the last occurence of a period "." in the input name. * if found, set 'pos' to that position. */ len = strlen(input_name); found = 0; pos = -1; for (i = len - 1; i > 0; i--) { if (input_name[i] == '.') { found = 1; pos = i; break; } } /* or, if no period, set 'pos' to end of string */ if (found == 0) { pos = len; } if ((pos + 1 + strlen(CALIBRATED_EXT)) >= MAX_FILENAME_LEN) { shError("set_calib_name: input_name %s is too long", input_name); return(1); } /* now, create the output name */ for (i = 0; i < pos; i++) { output_name[i] = input_name[i]; } output_name[i] = '\0'; sprintf(output_name, "%s%s", output_name, CALIBRATED_EXT); #ifdef DEBUG printf(" set_calib_name: input %s output %s \n", input_name, output_name); #endif return(0); } #endif /**************************************************************************** * PROCEDURE: find_passband_indices * * DESCRIPTION: Given an index to the 'filter[]' array in the detected * star file(s), figure out which other passbands we need to reduce * the data. We need * * a primary detected passband the index to which is given * as 'pri_index' arg * * a secondary detected passband we'll figure it out -- and set * the 'sec_index' arg * * a catalog passband we use to compare to primary * passband. We set the * 'cat_index' arg to this * * This is basically a big set of "if-then" statements. * * If the needed passbands aren't present in the detected or catalog * data, then we set the 'sec_index' and 'cat_index' args to -1, * and return with code 1. If we _did_ find the appropriate passbands, * we return with code 0. * * We perform a sanity check: we use the first catalog star, and * the first detected star in the first file, to figure out the * proper passband indices. Then, we check to see if all OTHER * catalog stars have the same passband in the 'cat_index', and * we check that all OTHER detected stars have the same primary and * secondary passbands in 'pri_index, and 'sec_index.' * We also count the number of stars with valid measurements in each * passband. * * RETURNS: * 0 if we could find the needed passbands * 1 if we could not (so we can't reduce the given * detected passband) * */ static int find_passband_indices ( int num_cat, /* I: number of catalog stars we've read */ S_CAT *cat_array, /* I: array of catalog stars we've read */ int num_detected_files, /* I: number of detected data files */ int *num_det, /* I: number of detected stars in each file */ S_DET *det_array[], /* I: array(s) of detected stars */ int pri_index, /* I: we want to reduce the passband which */ /* has this index in S_DET filter[] array */ int *sec_index, /* O: we will use this as index to secondary */ /* filter in S_DET filter[] array */ int *cat_index /* O: we will use this as index to catalog */ /* filter used as comparison */ ) { char pri_passband[FILTER_LENGTH + 1]; char sec_passband[FILTER_LENGTH + 1]; char cat_passband[FILTER_LENGTH + 1]; int i, i_star, i_file; int num_valid_det, num_valid_cat; int retval; S_CAT *one_cat; S_DET *one_det; /* * we don't really need ALL the catalog and detected stars; we are * interested only in the 'filter[]' values, and we'll assume that * the values in the first S_CAT and S_DET are identical to those * in all the structures. */ shAssert(num_cat > 0); shAssert(cat_array != NULL); shAssert(num_detected_files > 0); shAssert(num_det != NULL); shAssert(num_det[0] > 0); shAssert(det_array != NULL); shAssert(det_array[0] != NULL); one_cat = &(cat_array[0]); one_det = &(det_array[0][0]); *sec_index = -1; *cat_index = -1; /* * okay, here are the "rules" governing which passbands are used * to reduce which others. It would be nice to move these rules * into some sort of parameter file, and out of the code itself */ /* pri = "V", sec = "I", cat = "V" */ if (strcmp(one_det->filter[pri_index], "V") == 0) { strcpy(pri_passband, "V"); strcpy(sec_passband, "I"); strcpy(cat_passband, "V"); } /* pri = "I", sec = "V", cat = "I" */ if (strcmp(one_det->filter[pri_index], "I") == 0) { strcpy(pri_passband, "I"); strcpy(sec_passband, "V"); strcpy(cat_passband, "I"); } for (i = 0; i < MAX_PASSBANDS; i++) { if (strcmp(one_det->filter[i], sec_passband) == 0) { *sec_index = i; break; } } for (i = 0; i < MAX_PASSBANDS; i++) { if (strcmp(one_cat->filter[i], cat_passband) == 0) { *cat_index = i; break; } } if (valid_indices(pri_index, *sec_index, *cat_index) == 0) { #ifdef DEBUG printf("find_passband_indices FAIL: pri %d %s sec %d cat %d \n", pri_index, one_det->filter[pri_index], *sec_index, *cat_index); #endif return(1); } #ifdef DEBUG printf("find_passband_indices OKAY: pri %d %s sec %d %s cat %d %s \n", pri_index, one_det->filter[pri_index], *sec_index, one_det->filter[*sec_index], *cat_index, one_cat->filter[*cat_index]); #endif /* set these passband indices and names for future reference */ set_passband_values(pri_index, pri_index, *sec_index, *cat_index, one_det->filter[pri_index], one_det->filter[*sec_index], one_cat->filter[*cat_index]); /* * now, we check to make sure all catalog stars have proper value * in the 'sec_index' filter name, and all detected stars have * proper values in the 'pri_index' and 'sec_index' filter names. * * We also count the number of stars with valid magnitudes. */ retval = 0; num_valid_cat = 0; one_cat = cat_array; for (i_star = 0; i_star < num_cat; i_star++) { if (strncmp(one_cat->filter[*cat_index], cat_passband, FILTER_LENGTH) != 0) { shError("find_passband_indices: bad catalog filter %s in star %d \n", one_cat->filter[*cat_index], i_star); retval = 1; continue; } if (one_cat->mag[*cat_index] != BAD_MAG) { num_valid_cat++; } } num_valid_det = 0; for (i_file = 0; i_file < num_detected_files; i_file++) { one_det = det_array[i_file]; for (i_star = 0; i_star < num_det[i_file]; i_star++) { if (strncmp(one_det->filter[pri_index], pri_passband, FILTER_LENGTH) != 0) { shError("find_passband_indices: bad det pri filter %s in star %d \n", one_det->filter[pri_index], i_star); retval = 1; continue; } if (strncmp(one_det->filter[*sec_index], sec_passband, FILTER_LENGTH) != 0) { shError("find_passband_indices: bad det sec filter %s in star %d \n", one_det->filter[*sec_index], i_star); retval = 1; continue; } if ((one_det->mag[pri_index] != BAD_MAG) && (one_det->mag[*sec_index] != BAD_MAG)) { num_valid_det++; } } } #ifdef DEBUG printf("find_passband_indices: valid det %6d valid cat %6d \n", num_valid_det, num_valid_cat); #endif return(retval); } /**************************************************************************** * PROCEDURE: initialize_solution * * DESCRIPTION: initialize values in the S_SOLUTION structures * to indicate that we have not (yet) calculated any photometric * solution(s). If we later _do_ calculate a solution in some * passband, we'll re-set these values. * */ static void initialize_solution ( S_SOLUTION solution[], /* I/O: init values in these structs */ int num_images /* I: number of images for which we */ /* might solve for zeropoints */ ) { int i, j; for (i = 0; i < MAX_PASSBANDS; i++) { solution[i].num_obs = 0; solution[i].num_zeropts = num_images; for (j = 0; j < num_images; j++) { solution[i].a = (double *) shMalloc(num_images*sizeof(double)); solution[i].asig = (double *) shMalloc(num_images*sizeof(double)); } } } /**************************************************************************** * PROCEDURE: free_solution * * DESCRIPTION: de-allocate space in the "solution" structures */ static void free_solution ( S_SOLUTION solution[] /* I/O: free memory from these structs */ ) { int i; for (i = 0; i < MAX_PASSBANDS; i++) { shFree(solution[i].a); shFree(solution[i].asig); } } /**************************************************************************** * PROCEDURE: initialize_passbands * * DESCRIPTION: initialize values in the given S_PASSBANDS structure * to indicate that we have not (yet) figured out which passbands * were detected, and which can be calibrated. * */ static void initialize_passbands ( S_PASSBANDS *passbands /* I/O: init values in this structs */ ) { int i; for (i = 0; i < MAX_PASSBANDS; i++) { passbands->pri_index[i] = -1; passbands->sec_index[i] = -1; passbands->cat_index[i] = -1; strcpy(passbands->pri_passband[i], ""); strcpy(passbands->sec_passband[i], ""); strcpy(passbands->cat_passband[i], ""); } } /************************************************************************* * PROCEDURE: set_passband_values * * DESCRIPTION: set the values of one element of all the S_PASSBANDS * fields. This is done once, by the 'find_passband_indices' function, * and should not be done again. * * RETURNS: * nothing */ static void set_passband_values ( int index, /* I: set elements with this index */ int pri_value, /* I: value of primary detected filter index */ int sec_value, /* I: value of secondary detected filter index */ int cat_value, /* I: value of catalog filter index */ char *pri_passband, /* I: primary detected filter name */ char *sec_passband, /* I: secondary detected filter name */ char *cat_passband /* I: catalog filter name */ ) { shAssert((index >= 0) && (index < MAX_PASSBANDS)); shAssert((pri_value >= 0) && (pri_value < MAX_PASSBANDS)); shAssert((sec_value >= 0) && (sec_value < MAX_PASSBANDS)); shAssert((cat_value >= 0) && (cat_value < MAX_PASSBANDS)); shAssert((pri_passband != NULL) && (strlen(pri_passband) <= FILTER_LENGTH)); shAssert((sec_passband != NULL) && (strlen(sec_passband) <= FILTER_LENGTH)); shAssert((cat_passband != NULL) && (strlen(cat_passband) <= FILTER_LENGTH)); passband_library.pri_index[index] = pri_value; passband_library.sec_index[index] = sec_value; passband_library.cat_index[index] = cat_value; strncpy(passband_library.pri_passband[index], pri_passband, FILTER_LENGTH); strncpy(passband_library.sec_passband[index], sec_passband, FILTER_LENGTH); strncpy(passband_library.cat_passband[index], cat_passband, FILTER_LENGTH); passband_library.pri_passband[index][FILTER_LENGTH] = '\0'; passband_library.sec_passband[index][FILTER_LENGTH] = '\0'; passband_library.cat_passband[index][FILTER_LENGTH] = '\0'; } /************************************************************************* * PROCEDURE: get_passband_values * * DESCRIPTION: Given a primary filter index, place the values for the * primary, secondary and catalog indices used to reduce that primary * filter into the given arguments. Copy the names of these passbands * into the given arguments. * * This function must be called after 'set_passband_values', of course. * * RETURNS: * nothing */ static void get_passband_values ( int index, /* I: get elements with this index */ int *pri_value, /* I: value of primary detected filter index */ int *sec_value, /* I: value of secondary detected filter index */ int *cat_value, /* I: value of catalog filter index */ char *pri_passband, /* I: primary detected filter name */ char *sec_passband, /* I: secondary detected filter name */ char *cat_passband /* I: catalog filter name */ ) { shAssert((index >= 0) && (index < MAX_PASSBANDS)); shAssert(pri_value != NULL); shAssert(sec_value != NULL); shAssert(cat_value != NULL); shAssert(pri_passband != NULL); shAssert(sec_passband != NULL); shAssert(cat_passband != NULL); *pri_value = passband_library.pri_index[index]; *sec_value = passband_library.sec_index[index]; *cat_value = passband_library.cat_index[index]; strncpy(pri_passband, passband_library.pri_passband[index], FILTER_LENGTH); strncpy(sec_passband, passband_library.sec_passband[index], FILTER_LENGTH); strncpy(cat_passband, passband_library.cat_passband[index], FILTER_LENGTH); pri_passband[FILTER_LENGTH] = '\0'; sec_passband[FILTER_LENGTH] = '\0'; cat_passband[FILTER_LENGTH] = '\0'; } /*************************************************************************** * PROCEDURE: valid_indices * * DESCRIPTION: given a set of indices for the passbands * * pri_index index of primary band in detected files * sec_index index of secondary band in detected files * cat_index index of primary band in catalog files * * decide if the set is "valid"; that is, decide if it's possible to * create a photometric solution using these indices. * * We really just make sure that none are equal to -1, which signals * "there's no data for a passband." * * RETURNS: * 1 if the indices are valid (this is good) * 0 if the indices are invalid (this is bad) */ static int valid_indices ( int pri_index, /* I: index of primary detected measurements */ int sec_index, /* I: index of secondary detected measurements */ int cat_index /* I: index of primary catalog measurements */ ) { if ((pri_index == -1) || (sec_index == -1) || (cat_index == -1)) { return(0); } else { return(1); } } /*************************************************************************** * PROCEDURE: put_star_in_solution * * DESCRIPTION: given a detected star, we check to see if we should * include it in the calculation for the photometric * solution. In order to be included, a star must have * * include field set to 1 * cat_match_id field > 0 (has a matching catalog star) * cat mag not BAD_MAG * weight != PHOTOM_TINY_WEIGHT * * * RETURNS: * 1 if satisfies all conditions * (do include star in solution) * 0 if fails any condition * (don't include star in solution) */ static int put_star_in_solution ( S_DET *det_star, /* I: include this detected star? */ S_CAT *cat_array, /* I: array of catalog star structures */ int pri_index, /* I: index of primary mag in detected star */ int sec_index, /* I: index of secondary mag in detected star */ int cat_index /* I: index of primary mag for catalog star */ ) { S_CAT *cat_star; double weight; /* ignore stars with no matches */ if (det_star->cat_match_id < 0) { return(0); } /* ignore stars the user has marked as bad or saturated */ if (det_star->include <= 0) { return(0); } weight = assign_weight(det_star, pri_index, sec_index); /* ignore stars with excessive uncertainties */ if (weight == PHOTOM_TINY_WEIGHT) { return(0); } cat_star = &(cat_array[det_star->cat_match_id]); shAssert(cat_star != NULL); /* ignore stars with no good reference magnitude */ if (cat_star->mag[cat_index] == BAD_MAG) { return(0); } return(1); }