/* * match: a package to match lists of stars (or other items) * Copyright (C) 2000 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. * */ /* * little support functions for the matching code * */ #include #include #include #include #include "misc.h" #include "atpmatch.h" /* * this variable holds the value of the TRANS order used in this * instance of the program. It signals whether we're using * linear, quadratic, or cubic terms in the transformation. * * See * atTransOrderSet * atTransOrderGet */ static int at_trans_order = -1; /********************************************************************* * ROUTINE: shMalloc * * Attempt to allocate the given number of bytes. Return the * memory, if we succeeded, or print an error message and * exit with error code if we failed. * * RETURNS: * void * to new memory, if we got it */ void * shMalloc ( int nbytes /* I: allocate a chunk of this many bytes */ ) { void *vptr; if ((vptr = (void *) malloc(nbytes)) == NULL) { shError("shMalloc: failed to allocate for %d bytes", nbytes); exit(1); } return(vptr); } /********************************************************************* * ROUTINE: shFree * * Attempt to free the given piece of memory. * * RETURNS: * nothing */ void shFree ( void *vptr /* I: free this chunk of memory */ ) { free(vptr); } /********************************************************************* * ROUTINE: shError * * Print the given error message to stderr, but continue to execute. * * RETURNS: * nothing */ void shError ( char *format, /* I: format part of printf statement */ ... /* I: optional arguments to printf */ ) { va_list ap; va_start(ap, format); (void) vfprintf(stderr, (const char *)format, ap); fputc('\n', stderr); fflush(stdout); fflush(stderr); va_end(ap); } /********************************************************************* * ROUTINE: shFatal * * Print the given error message to stderr, and halt program execution. * * RETURNS: * nothing */ void shFatal ( char *format, /* I: format part of printf statement */ ... /* I: optional arguments to printf */ ) { va_list ap; va_start(ap, format); (void) vfprintf(stderr, (const char *)format, ap); fputc('\n', stderr); fflush(stdout); fflush(stderr); va_end(ap); exit(1); } /********************************************************************* * ROUTINE: shDebugSet, shDebug * * shDebugSet: sets the static variable 'debug_level' to the * given integer; it starts at 0, but the user may set it to * higher levels. The higher the level, the more messages * may be printed during execution. * * shDebug: If the current 'debug level' is >= the passed 'level', * then print the given message to stdout, and continue execution. * Otherwise, just continue. * * RETURNS: * nothing */ static int debug_level = 0; void shDebugSet ( int level /* I: set debug level to this value */ ) { debug_level = level; } void shDebug ( int level, /* I: debug level at which we print this */ char *format, /* I: format part of printf statement */ ... /* I: optional arguments to printf */ ) { va_list ap; if (level > debug_level) { return; } va_start(ap, format); (void) vfprintf(stdout, (const char *)format, ap); fputc('\n', stdout); fflush(stdout); va_end(ap); } /************************************************************************ * ROUTINE: atTransOrderSet * * DESCRIPTION: * Set the value of the order we'll use for TRANS structures. * Possibilities are: * * AT_TRANS_LINEAR linear transformation * AT_TRANS_QUADRATIC linear plus quadratic terms * AT_TRANS_CUBIC linear plus quadratic plus cubic terms * * RETURN: * nothing * * */ void atTransOrderSet ( int order /* I: order for all TRANS structures */ ) { at_trans_order = order; } /************************************************************************ * ROUTINE: atTransOrderGet * * DESCRIPTION: * Get the value of the order we're using in this instance of the program. * Possibilities are: * * AT_TRANS_LINEAR linear transformation * AT_TRANS_QUADRATIC linear plus quadratic terms * AT_TRANS_CUBIC linear plus quadratic plus cubic terms * * RETURN: * the order value * * */ int atTransOrderGet ( ) { /* sanity check -- make sure it's been set */ if (at_trans_order == -1) { shFatal("atTransOrderGet: at_trans_order not set yet"); } return(at_trans_order); } /************************************************************************ * * * ROUTINE: atTransNew * * DESCRIPTION: * Create a new TRANS structure, and return a pointer to it. * * RETURN: * TRANS * if all goes well * NULL if not * * */ TRANS * atTransNew ( ) { TRANS *new; static int id_number = 0; new = shMalloc(sizeof(TRANS)); new->id = id_number++; new->order = atTransOrderGet(); return(new); } /************************************************************************ * * * ROUTINE: atTransDel * * DESCRIPTION: * Delete the given TRANS structure * * RETURN: * nothing * * */ void atTransDel ( TRANS *trans /* I: structure to delete */ ) { shFree(trans); } /************************************************************************ * * ROUTINE: atStarNew * * DESCRIPTION: * Create a new s_star structure, and return a pointer to it. Fill * it with values 'x', 'y' and 'mag'. * * RETURN: * s_star * if all goes well * NULL if not * * */ struct s_star * atStarNew ( double x, /* I: x value for new star */ double y, /* I: y value for new star */ double mag /* I: mag value for new star */ ) { struct s_star *new; static int id_number = 0; new = (struct s_star *) shMalloc(sizeof(struct s_star)); new->id = id_number++; new->index = -1; new->x = x; new->y = y; new->mag = mag; new->match_id = -1; new->next = (struct s_star *) NULL; return(new); } /********************************************************************* * ROUTINE: read_star_file * * Given the name of a file, and three integers which specify the * columns in which the "x", "y", and "mag" data appear, go through * the file, creating an "s_star" structure for each line. * Place all the stars in a linked list, and return a pointer * to the head of the list, and the number of stars in the list. * * Ignore any line which starts with a COMMENT_CHAR, and any * completely empty line (one with whitespace only). * * If the "idcolumn" is not -1, then read ID numbers for stars from * that column, and use those ID values instead of generating them * internally. * * Note that columns are numbered starting at 0. Thus, * given a file with format like this: * * # x y mag * 234.43 54.33 12.23 * * we have xcolumn = 0, ycolumn = 1, magcolumn = 2. * * If the argument 'ra_hours_col' is >= 0, then it indicates that * the given colunm has Right Ascension values which are in hours * (rather than degrees). In that case, we multiply the column * value by 15.0 to convert from hours -> degrees. * * RETURNS: * SH_SUCCESS if all goes well * SH_GENERIC_ERROR if not */ int read_star_file ( char *filename, /* I: name of file */ int xcolumn, /* I: column in which 'x' data is found */ int ycolumn, /* I: column in which 'y' data is found */ int magcolumn, /* I: column in which 'mag' data is found */ int idcolumn, /* I: column in which 'id' data is found */ int ra_hours_col, /* I: index of column with RA data in hours */ /* if -1, no col has RA in hours */ int *num_stars, /* O: number of stars in new list goes here */ struct s_star **list /* O: new star list will be placed here */ ) { FILE *fp; char line[LINELEN]; char col[MAX_DATA_COL + 1][MAX_COL_LENGTH + 1]; int nline, num, ncol; int last_column = -1; int idval; struct s_star *head, *last, *new; double xval, yval, magval, double_idval; /* sanity checks */ shAssert(xcolumn >= 0); shAssert((ycolumn >= 0) && (ycolumn != xcolumn)); shAssert((magcolumn >= 0) && (magcolumn != xcolumn) && (magcolumn != ycolumn)); if ((fp = fopen(filename, "r")) == NULL) { shError("read_star_file: can't open file %s\n", filename); return(SH_GENERIC_ERROR); } /* find the highest-numbered column we need to read */ last_column = xcolumn; if (ycolumn > last_column) { last_column = ycolumn; } if (magcolumn > last_column) { last_column = magcolumn; } if (idcolumn > last_column) { last_column = idcolumn; } shAssert(last_column >= 2); if (last_column > MAX_DATA_COL) { shError("read_star_file: only %d columns allowed", MAX_DATA_COL); return(SH_GENERIC_ERROR); } nline = 0; head = (struct s_star *) NULL; last = head; num = 0; while (fgets(line, LINELEN, fp) != NULL) { if (line[0] == COMMENT_CHAR) { nline++; continue; } if (is_blank(line)) { nline++; continue; } ncol = sscanf(line, "%s %s %s %s %s %s %s %s %s %s", &(col[0][0]), &(col[1][0]), &(col[2][0]), &(col[3][0]), &(col[4][0]), &(col[5][0]), &(col[6][0]), &(col[7][0]), &(col[8][0]), &(col[9][0])); if (last_column > ncol) { shError("read_data_file: not enough entries in following line; skipping"); shError(" %s", line); nline++; continue; } /* now read values from each column */ if (get_value(col[xcolumn], &xval) != SH_SUCCESS) { shError("read_data_file: can't read X value from %s; skipping", col[xcolumn]); nline++; continue; } if (ra_hours_col == xcolumn) { xval *= 15.0; } if (get_value(col[ycolumn], &yval) != SH_SUCCESS) { shError("read_data_file: can't read Y value from %s; skipping", col[ycolumn]); nline++; continue; } if (ra_hours_col == ycolumn) { yval *= 15.0; } if (get_value(col[magcolumn], &magval) != SH_SUCCESS) { shError("read_data_file: can't read mag value from %s; skipping", col[magcolumn]); nline++; continue; } if (idcolumn != -1) { if (get_value(col[idcolumn], &double_idval) != SH_SUCCESS) { shError("read_data_file: can't read id value from %s; skipping", col[idcolumn]); nline++; continue; } else { idval = (int) double_idval; } } /* okay, it's safe to create a new s_star for this line */ nline++; num++; new = atStarNew(xval, yval, magval); if (idcolumn != -1) { new->id = idval; } if (head == NULL) { head = new; last = new; } else { last->next = new; last = new; } } *num_stars = num; *list = head; return(SH_SUCCESS); } /********************************************************************* * ROUTINE: read_matched_file * * Given the name of a file -- which has been created by the * 'match' program just a bit earlier -- read in data for * stars, creating an "s_star" structure for each line. * Place all the stars in a linked list, and return a pointer * to the head of the list, and the number of stars in the list. * * The format of the file is exactly like this: * * 61 2422.000000 -1175.000000 7.40 * * where the first column is an ID number * second x position (in coord system B) * third y position * fourth magnitude * * RETURNS: * SH_SUCCESS if all goes well * SH_GENERIC_ERROR if not */ int read_matched_file ( char *filename, /* I: name of file */ int *num_stars, /* O: number of stars in new list goes here */ struct s_star **list /* O: new star list will be placed here */ ) { FILE *fp; char line[LINELEN]; int nline, num, ncol; int idval; struct s_star *head, *last, *new; double xval, yval, magval; if ((fp = fopen(filename, "r")) == NULL) { shError("read_matched_file: can't open file %s\n", filename); return(SH_GENERIC_ERROR); } nline = 0; head = (struct s_star *) NULL; last = head; num = 0; while (fgets(line, LINELEN, fp) != NULL) { if (line[0] == COMMENT_CHAR) { nline++; continue; } if (is_blank(line)) { nline++; continue; } ncol = sscanf(line, "%d %lf %lf %lf", &idval, &xval, &yval, &magval); if (ncol != 4) { shError("read_matched_file: bad line; skipping"); shError(" %s", line); nline++; continue; } /* okay, it's safe to create a new s_star for this line */ nline++; num++; new = atStarNew(xval, yval, magval); new->id = idval; if (head == NULL) { head = new; last = new; } else { last->next = new; last = new; } } *num_stars = num; *list = head; return(SH_SUCCESS); } /********************************************************************** * ROUTINE: is_blank * * If the given string consists only of whitespace, return 1. * Otherwise, return 0. */ int is_blank ( char *line /* I: string to be checked */ ) { char *p; for (p = line; (*p != '\0') && (isspace(*p)); p++) { ; } if (p == '\0') { return(1); } else { return(0); } } /********************************************************************** * ROUTINE: get_value * * Given a string containing a numerical value, read the numerical * value and place it into the given double argument. * * Return 0 if all goes well. * Return 1 if there is an error. */ int get_value ( char *str, /* I: string to be converted to double */ double *val /* O: place value here */ ) { if (sscanf(str, "%lf", val) != 1) { return(1); } else { return(0); } } /************************************************************************ * * * ROUTINE: print_trans * * DESCRIPTION: * Print the elements of a TRANS structure. * * RETURNS: * nothing * * */ void print_trans ( TRANS *trans /* I: TRANS to print out */ ) { switch (trans->order) { case 1: /* linear transformation */ printf("TRANS: a=%-15.9f b=%-15.9f c=%-15.9f d=%-15.9f e=%-15.9f f=%-15.9f\n", trans->a, trans->b, trans->c, trans->d, trans->e, trans->f); break; case 2: /* quadratic terms */ printf("TRANS: a=%-15.9f b=%-15.9f c=%-15.9f d=%-15.9f e=%-15.9f f=%-15.9f ", trans->a, trans->b, trans->c, trans->d, trans->e, trans->f); printf(" g=%-15.9f h=%-15.9f i=%-15.9f j=%-15.9f k=%-15.9f l=%-15.9f\n", trans->g, trans->h, trans->i, trans->j, trans->k, trans->l); break; case 3: /* cubic terms */ printf("TRANS: a=%-15.9f b=%-15.9f c=%-15.9f d=%-15.9f e=%-15.9f f=%-15.9f g=%-15.9f h=%-15.9f", trans->a, trans->b, trans->c, trans->d, trans->e, trans->f, trans->g, trans->h); printf(" i=%-15.9f j=%-15.9f k=%-15.9f l=%-15.9f m=%-15.9f n=%-15.9f o=%-15.9f p=%-15.9f\n", trans->i, trans->j, trans->k, trans->l, trans->m, trans->n, trans->o, trans->p); break; default: shFatal("print_trans: invalid trans->order %d \n", trans->order); exit(1); } }