IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3339


Ignore:
Timestamp:
Feb 28, 2005, 11:07:26 AM (21 years ago)
Author:
eugene
Message:

working towards addstar-2-0

Location:
trunk/Ohana/src/addstar
Files:
22 added
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/addstar/Makefile

    r2510 r3339  
    1919
    2020ADDSTAR = \
    21 $(SRC)/addstar.$(ARCH).o        $(SRC)/gcatalog.$(ARCH).o   \
    22 $(SRC)/gregions.$(ARCH).o       $(SRC)/gimages.$(ARCH).o    \
    23 $(SRC)/wimage.$(ARCH).o         $(SRC)/gstars.$(ARCH).o     \
    24 $(SRC)/find_matches.$(ARCH).o   $(SRC)/wcatalog.$(ARCH).o   \
    25 $(SRC)/sort_lists.$(ARCH).o     $(SRC)/ConfigInit.$(ARCH).o     \
    26 $(SRC)/aregion.$(ARCH).o        \
    27 $(SRC)/check_permissions.$(ARCH).o $(SRC)/calibrate.$(ARCH).o \
    28 $(SRC)/args.$(ARCH).o           $(SRC)/parse_time.$(ARCH).o \
    29 $(SRC)/mkcatalog.$(ARCH).o
     21$(SRC)/addstar.$(ARCH).o        $(SRC)/gcatalog.$(ARCH).o     \
     22$(SRC)/gregion_image.$(ARCH).o  $(SRC)/gregion_star.$(ARCH).o \
     23$(SRC)/gimages.$(ARCH).o        $(SRC)/wimage.$(ARCH).o       \
     24$(SRC)/gstars.$(ARCH).o         $(SRC)/grefstars.$(ARCH).o    \
     25$(SRC)/find_matches.$(ARCH).o   $(SRC)/find_matches_refstars.$(ARCH).o \
     26$(SRC)/wcatalog.$(ARCH).o       $(SRC)/conversions.$(ARCH).o     \
     27$(SRC)/sort_lists.$(ARCH).o     $(SRC)/ConfigInit.$(ARCH).o   \
     28$(SRC)/aregion.$(ARCH).o        $(SRC)/find_subset.$(ARCH).o  \
     29$(SRC)/load_subpix.$(ARCH).o    $(SRC)/airmass.$(ARCH).o      \
     30$(SRC)/mkcatalog.$(ARCH).o      $(SRC)/calibrate.$(ARCH).o    \
     31$(SRC)/args.$(ARCH).o           $(SRC)/parse_time.$(ARCH).o   \
     32$(SRC)/match_refstars.$(ARCH).o         \
     33$(SRC)/check_permissions.$(ARCH).o
    3034
    3135OBJ = $(ADDSTAR)

  • trunk/Ohana/src/addstar/include/addstar.h

    r2457 r3339  
    2424double DEFAULT_RADIUS, NSIGMA, SNLIMIT;
    2525double ZeroPt;
     26double Latitude, SiderealTime;
    2627int    VERBOSE;
    2728int    SKIP_MISSED;
    … …  
    4041char AirmassKeyword[64];
    4142char CCDNumKeyword[64];
     43char STKeyword[64];
     44char SubpixDatafile[256];
    4245
    4346int CalReference;
    … …  
    4750int ONLY_IMAGES;
    4851int FORCE_READ;
     52int SKYPROBE;
     53int ADDREFS;
     54int REPLACE;
     55int ONLY_MATCH;
    4956
     57time_t TIMEREF;
     58
     59int main (int argc, char **argv);
    5060int Shutdown (); 
    5161void TrapSignal (int sig);
    52 int main (int argc, char **argv);
    5362int SetSignals ();
     63void help ();
     64int args (int argc, char **argv);
     65void SetProtect (int mode);
     66
    5467int gcatalog (Catalog *catalog);
    55 GSCRegion *gregions (Image *image, int *Nregions);
     68void gregion_star (Stars *star, GSCRegion *region);
     69GSCRegion *gregion_image (Image *image, int *Nregions);
    5670Image *gimages (FILE *f, Image *image, int *Npimage);
    5771int edge_check (double *x1, double *y1, double *x2, double *y2);
    … …  
    5973void wimage (FILE *f, int dbstate, Image *image, int Nstars);
    6074Stars *gstars (char *file, int *NSTARS, Image *image);
     75void find_matches_refstars (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog);
    6176void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage);
    6277int in_image (double r, double d, Image *image);
    6378void wcatalog (Catalog *catalog);
     79void ConfigInit (int *argc, char **argv);
     80void aregion (GSCRegion *region, FILE *f, double ra, double dec);
     81void wimage (FILE *f, int dbstate, Image *image, int Nstars);
     82double airmass (double ra, double dec, double st, double latitude);
     83void InitCalibration ();
     84void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y);
     85void FindCalibration (Image *image);
     86double get_subpix (double x, double y);
     87void load_subpix ();
     88double scat_subpix (double x, double y);
     89void mkcatalog (GSCRegion *region, Catalog *catalog);
     90Stars *grefstars (char *file, int *Nstars);
     91Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, int *NSTARS);
     92int match_refstars (Stars *stars, int Nstars);
     93int parse_time (Header *header);
     94
     95void check_permissions (char *basefile);
     96void make_backup (char *filename);
     97void uppercase (char *string);
     98
    6499void sort_lists (float *X, float *Y, int *S, int N);
    65100void fsort (float *X, int N);
    66 void ConfigInit (int *argc, char **argv);
     101void fsort2 (float *X, float *Y, int N);
     102
     103char *sec_to_date (time_t second);
     104double sec_to_jd (time_t second);
     105int chk_time (char *line);
     106int dms_to_ddd (double *Value, char *string);
    67107int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string);
    68 void aregion (GSCRegion *region, FILE *f, double ra, double dec);
    69 void wimage (FILE *f, int dbstate, Image *image, int Nstars);
    70 void check_permissions (char *basefile);
    71 void make_backup (char *filename);
    72 void InitCalibration ();
    73 void SaveCalibration (float M, float dM, float Mr, float Mc, float A, int N, float ra, float dec, float x, float y);
    74 void FindCalibration (Image *image);
    75 void help ();
    76 int args (int argc, char **argv);
    77 int parse_time (Header *header);
    78 void uppercase (char *string);
    79 void mkcatalog (GSCRegion *region, Catalog *catalog);
     108int str_to_dtime (char *line, double *second);
     109int str_to_radec (double *ra, double *dec, char *str1, char *str2);
     110int str_to_time (char *line, time_t *second);
     111time_t jd_to_sec (double jd);
     112time_t date_to_sec (char *date);
     113
     114/**
     115  there is an inconsistency to be resolved: fixed structures (like Image)
     116  need a fixed bit-length time (e_time), but these functions all use the
     117  UNIX time_t types, which may be 32 or 64 bits, depending on the machine.
     118  This can be resolved by using time_t with these functions, but casting
     119  between e_time and time_t when necessary (ie, cannot return data to an
     120  e_time pointer from one of these functions)
     121**/

  • trunk/Ohana/src/addstar/src/ConfigInit.c

    r2457 r3339  
    4141  ScanConfig (config, "PHOTCODE_FILE",          "%s", 0, PhotCodeFile);
    4242
     43  /* used by parse_time to find time-related keywords */
    4344  ScanConfig (config, "DATE-KEYWORD",           "%s", 0, DateKeyword);
    4445  ScanConfig (config, "DATE-MODE",              "%s", 0, DateMode);
    … …  
    5152  ScanConfig (config, "CCDNUM-KEYWORD",         "%s", 0, CCDNumKeyword);
    5253
     54  ScanConfig (config, "ST-KEYWORD",             "%s",  0, STKeyword);
     55  ScanConfig (config, "OBSERVATORY-LATITUDE",   "%lf", 0, &Latitude);
     56  ScanConfig (config, "SUBPIX_DATAFILE",        "%s",  0, SubpixDatafile);
     57
    5358  if (!LoadPhotcodes (PhotCodeFile)) {
    5459    fprintf (stderr, "error loading photcodes\n");

  • trunk/Ohana/src/addstar/src/addstar.c

    r2457 r3339  
    11# include "addstar.h"
    2 
    3 /* these variables are needed by Shutdown */
    4 static FILE *f = (FILE *) NULL;
    5 static int Protect = FALSE;
    6 static int Trapped = FALSE;
    7 
    8 /* clean up open / locked ImageCat before shutting down */
    9 int Shutdown () { 
    10   int dbstate;
    11 
    12   Protect = TRUE;
    13   fclearlockfile (ImageCat, f, LCK_HARD, &dbstate);
    14   fprintf (stderr, "ERROR: addstar halted\n");
    15   exit (1);
    16 }
    17 
    18 void TrapSignal (int sig) {
    19     fprintf (stderr, "trapped signal %d, exiting gracefully\n", sig);
    20     if (sig == 11) {
    21       fprintf (stderr, "seg fault\n");
    22       exit (1);
    23     }
    24     if (Protect) {
    25       Trapped = TRUE;
    26       fprintf (stderr, "blocking until protected sections are clear\n");
    27       return;
    28     }
    29     Shutdown ();
    30 }   
    312
    323int main (int argc, char **argv) {
    … …  
    345  int i, mode, dbstate;
    356  int Nstars, Nimage, Nregions, Nmissed;
    36   Stars *stars, *gstars();
    37   Image image, *pimage, *gimages();
    38   GSCRegion *region, *gregions();
     7  Stars *stars;
     8  Image image, *image_db;
     9  GSCRegion *region;
    3910  Catalog catalog;
    40  
     11
    4112  SetSignals ();
     13  ConfigInit (&argc, argv);
    4214  args (argc, argv);
    4315
    44   /* load the input datafile */
    45   stars = gstars (argv[1], &Nstars, &image);
    46   if (Nstars == 0) {
    47     if (VERBOSE) fprintf (stderr, "no stars in data file, skipping\n");
    48     fprintf (stderr, "SUCCESS\n");
    49     exit (0);
     16  lock_image_db ();
     17
     18  switch (mode) {
     19  case IMAGE:
     20    stars = gstars (argv[1], &Nstars, &image);
     21    region = gregion_image (&image, &Nregions);
     22    image_db = gimages (&image, &Nimage_db);
     23    break;
     24  case REFLIST:
     25    stars = grefstars (argv[1], &Nstars);
     26    region = gregion_stars (stars, &Nregion);
     27    break;
     28  case PATCH :
     29    regions = gregion_patch (patch, &Nregions);
     30    break;
     31  case MATCH:
     32    regions = gregion_match (&Nregions);
     33    break;
    5034  }
    5135
    52   /* lock the image catalog */
    53   check_permissions (ImageCat);
    54   f = fsetlockfile (ImageCat, 3600.0, LCK_HARD, &dbstate);
    55   if ((f == (FILE *) NULL) && (dbstate != LCK_EMPTY)) {
    56     fprintf (stderr, "ERROR: can't lock image catalog\n");
    57     exit (1);
    58   }
    59   fseek (f, 0, SEEK_SET);
     36  for (i = 0; i < Nregions; i++) {
     37    load_catalog (&catalog, &region[i]);
    6038
    61   /* load images */
    62   if (dbstate == LCK_EMPTY) {
    63     Nimage = 0;
    64   } else {
    65     pimage = gimages (f, &image, &Nimage);
    66   }
    67   region = gregions (&image, &Nregions);
    68   fprintf (stderr, "region %f %f %f %f\n", region[0].RA[0], region[0].RA[1], region[0].DEC[0], region[0].DEC[1]);
    69 
    70   for (i = 0; i < Nregions; i++) {
    71     check_permissions (region[i].filename);
    72   }
    73 
    74   if (CALIBRATE) { InitCalibration (); }
    75 
    76   for (i = 0; i < Nregions; i++) {
    77     catalog.filename = region[i].filename;  /* don't free region before catalog! */
    78     fprintf (stderr, "adding to %s\n", region[i].filename);
    79    
    80     switch (lock_catalog (&catalog, LCK_XCLD)) {
    81     case 0:
    82       /* this is a serious error because ImageCat was NOT locked */
    83       fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename);
    84       exit (1);
    85     case 1:
    86       gcatalog (&catalog); /* load from disk */
     39    switch (mode) {
     40    case IMAGE:
     41      find_matches (&region[i], stars, Nstars, &catalog, &image, pimage, Nimage);
    8742      break;
    88     case 2:
    89       mkcatalog (&region[i], &catalog); /* fills in new header info */
     43    case PATCH:
     44    case MATCH:
     45      stars = grefcat (region, &Nstars);
     46    case REFLIST:
     47      subset = find_subset (&region, stars, Nstars, &subref, &Nsubset);
     48      find_matches_refstars (&region, subset, Nsubset, &catalog);
    9049      break;
    9150    }
    92     find_matches (&region[i], stars, Nstars, &catalog, &image, pimage, Nimage);
    93 
    94     /* protect wcatalog from interrupt signals */
    95     Protect = TRUE;
    96     if (!DUMP_MATCHES && !ONLY_IMAGES) wcatalog (&catalog);
    97     if (Trapped) Shutdown ();
    98     Protect = FALSE;
    99     unlock_catalog (&catalog);
     51    wcatalog (&catalog);
    10052  }
    10153
    … …  
    10355  if (DUMP_MATCHES) Shutdown ();
    10456
    105   for (Nmissed = i = 0; i < Nstars; i++) {
    106     if (stars[i].found == -1) {
    107       fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM);
    108       Nmissed ++;
    109     }
    110   }
    111   if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed);
    112 
    113   /* protect wimage from interrupt signals */
    114   Protect = TRUE;
    115   wimage (f, dbstate, &image, Nstars);
    116   if (Trapped) Shutdown ();
    117   Protect = FALSE;
    118 
    119   fclearlockfile (ImageCat, f, LCK_HARD, &dbstate);
    120 
    121   mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
    122   chmod (ImageCat, mode);
    123 
    124   fprintf (stderr, "SUCCESS\n");
     57  if (mode == IMAGE) wimage (f, dbstate, &image, Nstars);
     58  unlock_image_db ();
    12559  exit (0);
    126 
    12760}
    12861
    129 int SetSignals () {
     62/* names:
    13063
    131   int i;
    132 
    133   /* disable almost all signal interrupts */
    134   for (i = 0; i < 36; i++) {
    135     switch (i) {
    136       /* can't redirect this signals */
    137     case SIGKILL:    /* kill -9: cannot be caught or ignored */
    138     case SIGSTOP:    /* SIGSTOP: cannot be caught or ignored */
    139       /* ignore these signals */
    140     case SIGCHLD:    /* child halted: ignore */
    141     case SIGPWR:     /* power failure - why ignore this? */
    142     case SIGWINCH:   /* window resized */
    143     case SIGCONT:    /* continue - maintain this action */
    144     case SIGTSTP:    /* stop signal sent from tty - why ignore? */
    145     case SIGURG:     /* socket signal, ignore this */
    146       break;
    147      
    148     default:
    149       signal (i, TrapSignal);
    150     }
    151   }
    152   return (TRUE);
    153 }
    154 /*
    155 
    156        Signal     Value     Action   Comment
    157        -------------------------------------------------------------------------
    158        SIGHUP        1        A      Hangup detected on controlling terminal
    159                                      or death of controlling process
    160        SIGINT        2        A      Interrupt from keyboard
    161        SIGQUIT       3        A      Quit from keyboard
    162        SIGILL        4        A      Illegal Instruction
    163        SIGABRT       6        C      Abort signal from abort(3)
    164        SIGFPE        8        C      Floating point exception
    165        SIGKILL       9       AEF     Kill signal
    166        SIGSEGV      11        C      Invalid memory reference
    167        SIGPIPE      13        A      Broken pipe: write to pipe with no readers
    168        SIGALRM      14        A      Timer signal from alarm(2)
    169        SIGTERM      15        A      Termination signal
    170        SIGUSR1   30,10,16     A      User-defined signal 1
    171        SIGUSR2   31,12,17     A      User-defined signal 2
    172        SIGCHLD   20,17,18     B      Child stopped or terminated
    173        SIGCONT   19,18,25            Continue if stopped
    174        SIGSTOP   17,19,23    DEF     Stop process
    175        SIGTSTP   18,20,24     D      Stop typed at tty
    176        SIGTTIN   21,21,26     D      tty input for background process
    177        SIGTTOU   22,22,27     D      tty output for background process
    178 
    179        Next various other signals.
    180 
    181        Signal       Value     Action   Comment
    182        ---------------------------------------------------------------------
    183        SIGTRAP        5         CG     Trace/breakpoint trap
    184        SIGIOT         6         CG     IOT trap. A synonym for SIGABRT
    185        SIGEMT       7,-,7       G
    186        SIGBUS      10,7,10      AG     Bus error
    187        SIGSYS      12,-,12      G      Bad argument to routine (SVID)
    188        SIGSTKFLT    -,16,-      AG     Stack fault on coprocessor
    189        SIGURG      16,23,21     BG     Urgent condition on socket (4.2 BSD)
    190        SIGIO       23,29,22     AG     I/O now possible (4.2 BSD)
    191        SIGPOLL                  AG     A synonym for SIGIO (System V)
    192        SIGCLD       -,-,18      G      A synonym for SIGCHLD
    193        SIGXCPU     24,24,30     AG     CPU time limit exceeded (4.2 BSD)
    194        SIGXFSZ     25,25,31     AG     File size limit exceeded (4.2 BSD)
    195        SIGVTALRM   26,26,28     AG     Virtual alarm clock (4.2 BSD)
    196        SIGPROF     27,27,29     AG     Profile alarm clock
    197        SIGPWR      29,30,19     AG     Power failure (System V)
    198        SIGINFO      29,-,-      G      A synonym for SIGPWR
    199        SIGLOST      -,-,-       AG     File lock lost
    200        SIGWINCH    28,28,20     BG     Window resize signal (4.3 BSD, Sun)
    201        SIGUNUSED    -,31,-      AG     Unused signal
    202        (Here - denotes that a signal is absent; there where three values are given, the first one is usually  valid  for  alpha  and
    203        sparc,  the  middle  one  for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a
    204        sparc.)
    205 
    206        The letters in the "Action" column have the following meanings:
    207 
    208        A      Default action is to terminate the process.
    209 
    210        B      Default action is to ignore the signal.
    211 
    212        C      Default action is to dump core.
    213 
    214        D      Default action is to stop the process.
    215 
    216        E      Signal cannot be caught.
    217 
    218        F      Signal cannot be ignored.
    219 
    220        G      Not a POSIX.1 conformant signal.
    221 
     64   catalog - existing object db table
     65   region  - sky area which may or may contain data
     66   patch   - RA,DEC bounded portion of sky
     67   
    22268*/

  • trunk/Ohana/src/addstar/src/aregion.c

    r2457 r3339  
    8888
    8989  if (!done) {
    90     fprintf (stderr, "error in search: %f %f\n", ra, dec);
     90    fprintf (stderr, "ERROR in search: %f %f\n", ra, dec);
    9191    exit (1);
    9292  }

  • trunk/Ohana/src/addstar/src/args.c

    r2457 r3339  
    1515int args (int argc, char **argv) {
    1616 
    17   int N;
     17  int i, N;
    1818
    1919  /* check for help request */
    … …  
    2222    help ();
    2323  }
    24 
    25   /* configuration info */
    26   ConfigInit (&argc, argv);
    2724
    2825  /* check for command line options */
    … …  
    4340    remove_argument (N, &argc, argv);
    4441  }
     42  ONLY_IMAGES = FALSE;
     43  if ((N = get_argument (argc, argv, "-image"))) {
     44    ONLY_IMAGES = TRUE;
     45    remove_argument (N, &argc, argv);
     46  }
    4547  DUMP_MATCHES = FALSE;
    4648  if ((N = get_argument (argc, argv, "-dump"))) {
    … …  
    5052  FORCE_READ = FALSE;
    5153  if ((N = get_argument (argc, argv, "-force"))) {
     54    FORCE_READ = TRUE;
    5255    remove_argument (N, &argc, argv);
    53     FORCE_READ = TRUE;
    5456  }
    5557  SKIP_MISSED = FALSE;
    … …  
    5860    remove_argument (N, &argc, argv);
    5961  }
    60   ONLY_IMAGES = FALSE;
    61   if ((N = get_argument (argc, argv, "-image"))) {
    62     ONLY_IMAGES = TRUE;
     62  SKYPROBE = FALSE;
     63  if ((N = get_argument (argc, argv, "-skyprobe"))) {
     64    SKYPROBE = TRUE;
    6365    remove_argument (N, &argc, argv);
     66  }
     67  ADDREFS = FALSE;
     68  if ((N = get_argument (argc, argv, "-addrefs"))) {
     69    ADDREFS = TRUE;
     70    remove_argument (N, &argc, argv);
     71  }
     72  REPLACE = FALSE;
     73  if ((i = get_argument (argc, argv, "-replace"))) {
     74    REPLACE = TRUE;
     75    remove_argument (i, &argc, argv);
     76  }
     77  TIMEREF = 0;
     78  if ((i = get_argument (argc, argv, "-time"))) {
     79    remove_argument (i, &argc, argv);
     80    if (!str_to_time (argv[i], &TIMEREF)) {
     81      fprintf (stderr, "syntax error in time\n");
     82      exit (1);
     83    }
     84    remove_argument (i, &argc, argv);
     85  }
     86  ONLY_MATCH = FALSE;
     87  if ((i = get_argument (argc, argv, "-only-match"))) {
     88    ONLY_MATCH = TRUE;
     89    remove_argument (i, &argc, argv);
    6490  }
    6591  CALIBRATE = FALSE;
    … …  
    85111  }
    86112
     113  if (SKYPROBE) load_subpix ();
     114
    87115  if (argc != 2) {
    88116    fprintf (stderr, "USAGE: addstar (filename) [-dump] [-image] [-v] [-p photcode] [-accept] [-cal code1 code2]\n");
    … …  
    92120}
    93121
     122/** addstar modes:
     123 
     124    addstar (image.smp)  - add cmp/smp image data to db
     125    addstar -ref (file.dat) (photcode)
     126    addstar -cat (USNO/2MASS/GSC) -region (ra dec - ra dec)
     127
     128    -replace : ref/cat - replace existing match (photcode/time)
     129    -match   : ref/cat - only add measures to existing averages
     130
     131    ref types:
     132    ASCII - RA,DEC,M,dM in a table
     133
     134**/
     135

  • trunk/Ohana/src/addstar/src/calibrate.c

    r2457 r3339  
    22
    33static int Ncal, NCAL;
    4 static float *Mobs, *dMobs, *Mref, *Cref, *Airm;
     4static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm;
    55static float *RA, *DEC, *X, *Y;
    66int *Nstar;
    … …  
    1414    ALLOCATE (dMobs, float, NCAL);
    1515    ALLOCATE (Mref,  float,  NCAL);
     16    ALLOCATE (dMref,  float,  NCAL);
    1617    ALLOCATE (Cref,  float,  NCAL);
    1718    ALLOCATE (Airm,  float,  NCAL);
    … …  
    2324}
    2425 
    25 void SaveCalibration (float M, float dM, float Mr, float Mc, float A, int N, float ra, float dec, float x, float y) {
     26void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y) {
    2627
    2728  Mobs[Ncal] = M;
    2829  dMobs[Ncal] = dM;
    2930  Mref[Ncal] = Mr;
     31  dMref[Ncal] = dMr;
    3032  Cref[Ncal] = Mc;
    3133  Airm[Ncal] = A;
    … …  
    4345    REALLOCATE (Mobs,  float, NCAL);
    4446    REALLOCATE (Mref,  float, NCAL);
     47    REALLOCATE (dMref,  float, NCAL);
    4548    REALLOCATE (Cref,  float, NCAL);
    4649    REALLOCATE (Airm,  float, NCAL);
    … …  
    5659
    5760  int i, MaxN, *Nlist, Nkeep;
    58   float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag;
    59   float dAs, dMs, dYs;
     61  float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag, *dDmag;
     62  float dMo, dMr, Mw, Dmed, W1, W2, NSigma;
    6063 
    6164  /* reject multiple matched-stars */
    … …  
    8083  /* accumulate delta mags */
    8184  ALLOCATE (Dmag, float, Ncal);
     85  ALLOCATE (dDmag, float, Ncal);
    8286  Nkeep = 0;
    8387  Klam = thiscode[0].K;
    … …  
    8993    Mabs = 0.001*(Mobs[i] + Klam*(Airm[i] - 1000) + Clam + Xlam*(Mref[i] - Cref[i])) - ZeroPt;
    9094
    91     /* subpix correction : make parameters global! */
    92     dAs = 0;
    93     if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i];
    94     dYs = Y[i] - (int) (Y[i]) + 0.55;
    95     dMs = dAs * sin (6.28*dYs);
    96     /* Mabs -= dMs; */
     95    /* note: subpix correction is applied in gstars */
    9796
    9897    if (DUMP_MATCHES)
    99       fprintf (stdout, "%d  %6.3f %6.3f %6.3f  %10.6f %10.6f  %7.2f %7.2f\n",
    100                i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i]);
     98      fprintf (stdout, "%d  %6.3f %6.3f %6.3f  %10.6f %10.6f  %7.2f %7.2f  %7.2f  %7.2f\n",
     99               i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i], Airm[i], 0.001*dMobs[i]);
    101100
    102101    /* skip stars brighter than 8.0 */
    103     /* if (Mabs < 8.0) continue; */
     102    if (SKYPROBE) {
     103      if (Mabs > 9.0) continue;
     104      if (Mabs < 5.0) continue;
     105    }
    104106   
     107    dMr = MAX (0.005, 0.001*dMref[i]);
     108    dMo = MAX (0.005, 0.001*dMobs[i]);
     109
    105110    Dmag[Nkeep] = (Mabs - 0.001*Mref[i]);
     111    dDmag[Nkeep] = (dMr*dMr + dMo*dMo);
    106112    Nkeep ++;
    107113  }
    108114
     115  if (Nkeep < 5) {
     116    fprintf (stderr, "too few stars\n");
     117    image[0].Mcal = 10000;
     118    image[0].dMcal = 10000;
     119    return;
     120  }
     121  fsort2 (Dmag, dDmag, Nkeep);
     122
     123  /* take sort list of Dmag, find median */
     124  Dmed = Dmag[(int)(0.5*Nkeep)];
     125
     126  /* exclude points with abs(Dmag - Dmed) / dDmag > 2.5 */
     127
     128  /* accumulate delta mags (25% - 75% of clipped range) */
     129  W1 = 0.0;
     130  W2 = 0.0;
     131  M1 = 0.0;
     132  M2 = 0.0;
     133  N  = 0.0;
     134  for (i = 0; i < Nkeep; i++) {
     135    NSigma = fabs (Dmag[i] - Dmed) / sqrt (dDmag[i]);
     136    if (NSigma > 2.5) continue;
     137    W1 += Dmag[i] / dDmag[i];
     138    W2 += 1 / dDmag[i];
     139    M1 += Dmag[i];
     140    M2 += SQ (Dmag[i]);
     141    N  += 1.0;
     142  }
     143
     144  if (N > 1) {
     145    M1 = M1 / N;
     146    M2 = sqrt (fabs(M2/N - M1*M1));
     147    Mw = W1 / W2;
     148    fprintf (stderr, "N: %.0f, mean: %f, wt mean: %f, stdev: %f, precision: %f\n", N, M1, Mw, M2, M2 / sqrt (N));
     149    image[0].Mcal = 1000 * M1;
     150    image[0].dMcal = 1000 * M2 / sqrt (N);
     151    image[0].Mxxxx = N;
     152  } else {
     153    fprintf (stderr, "too few stars\n");
     154    image[0].Mcal = 10000;
     155    image[0].dMcal = 10000;
     156    image[0].Mxxxx = 0;
     157  }
     158}
     159
     160# if (0)
     161
     162/* old correction - superceeded by subpix grid */
     163
     164  /* find lo (Dmed - 0.1) + hi (Dmed + 0.1)
     165  lo = -1; hi = -1;
     166  for (i = 0; (lo == -1) && (i < Nkeep); i++) { if (Dmed - Dmag[i] < 0.1) lo = i; }
     167  for (i = Nkeep - 1; (hi == -1) && (i >= 0); i--) { if (Dmag[i] - Dmed < 0.1) hi = i; }
     168  if (lo == -1) lo = 0;
     169  if (hi == -1) hi = Nkeep - 1;
     170  */
     171
    109172  if (Nkeep < 3) {
    110173    fprintf (stderr, "too few stars\n");
    … …  
    113176    return;
    114177  }
    115   fsort (Dmag, Nkeep);
     178  fsort2 (Dmag, dDmag, Nkeep);
    116179
    117180  /* accumulate delta mags (25% - 75% of list) */
    … …  
    139202    image[0].dMcal = 10000;
    140203  }
    141 }
    142 
    143 
     204
     205
     206    /* subpix correction : make parameters global! */
     207    dAs = 0;
     208    if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i];
     209    dYs = Y[i] - (int) (Y[i]) + 0.55;
     210    dMs = dAs * sin (6.28*dYs);
     211    /* Mabs -= dMs; */
     212
     213# endif

  • trunk/Ohana/src/addstar/src/find_matches.c

    r2457 r3339  
    77void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage) {
    88
     9  double secz;
    910  int i, j, k, n, m, N, first_j, found, found0, found1;
    1011  double X, Y, RADIUS, RADIUS2;
    11   float *X1, *Y1, *X2, *Y2, CalM, CalC, *Xs, *Ys;
     12  float *X1, *Y1, *X2, *Y2, CalM, CalC, dCalM, *Xs, *Ys;
    1213  float dX, dY, dR;
    1314  int *N1, *N2,  *next, *next_miss, last, last_miss;
    … …  
    1920  short int Mrel, *Mval;
    2021
    21   /** allocate local arrays **/
     22  /* photcode data */
     23  CalM = CalC = dCalM = 0;
     24  Nsecfilt = GetPhotcodeNsecfilt ();
     25  Nsec = GetPhotcodeNsec (thiscode[0].code);
     26
     27  /** allocate local arrays (stars) **/
    2228  ALLOCATE (X1, float, Nstars);
    2329  ALLOCATE (Y1, float, Nstars);
    24   ALLOCATE (N1, int, Nstars);
    25 
    26   /** allocate local arrays **/
    27 
    28   Nsecfilt = GetPhotcodeNsecfilt ();
    29 
     30  ALLOCATE (N1, int,   Nstars);
     31
     32  /** allocate local arrays (catalog) **/
    3033  Nave = catalog[0].Naverage;
    3134  NAVE = Nave + 1000;
    3235  ALLOCATE (X2, float, NAVE);
    3336  ALLOCATE (Y2, float, NAVE);
    34   ALLOCATE (N2, int, NAVE);
     37  ALLOCATE (N2, int,   NAVE);
    3538  ALLOCATE (Xs, float, NAVE);
    3639  ALLOCATE (Ys, float, NAVE);
    … …  
    5154  REALLOCATE (catalog[0].missing, Missing, NMISS);
    5255 
    53   /* check on photcode */
    54   Nsec = GetPhotcodeNsec (thiscode[0].code);
    55 
    5656  /* project onto rectilinear grid with 1 arcsec pixels, sort by X */
    5757  /* reference for coords is this image */
    … …  
    6363 
    6464  for (i = 0; i < Nstars; i++) {
    65     fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords);
     65    fRD_to_XY (&X1[i], &Y1[i], stars[0][i].R, stars[0][i].D, &tcoords);
    6666    N1[i] = i;
    6767  }
    … …  
    9090
    9191  /* choose a radius for matches */
    92   /* 0.02 corrects cerror to arcsec from storage units */
    9392  if (DEFAULT_RADIUS == 0) {
    94     RADIUS = NSIGMA * 0.02 * image[0].cerror;
     93    RADIUS = NSIGMA * 0.02 * image[0].cerror;  /* 0.02 corrects cerror to arcsec from storage units */
    9594  } else {
    9695    RADIUS = DEFAULT_RADIUS;
    9796  }
    9897  RADIUS2 = RADIUS*RADIUS;
     98
    9999  /* correct instrumental mags for exposure time */
    100100  if (image[0].exptime > 0) {
    … …  
    108108   
    109109    dX = X1[i] - X2[j];
    110     /* fprintf (stderr, "%f   %f %f  %f %f    %f %f  %f %f\n", dX, X1[i], Y1[i], X2[j], Y2[j], stars[i].R, stars[i].D, catalog[0].average[j].R, catalog[0].average[j].D); */
    111110
    112111    if (dX <= -2*RADIUS) {
    … …  
    127126      if (dR < RADIUS2) {  /* new measurement of this star */
    128127        Nmatch ++;
    129         /** insert star in measurement list */
    130128        n = N2[j];
    131129        N = N1[i];
    132         m = catalog[0].average[n].offset;  /* first measurement of this star */
    133         for (k = 0; k < catalog[0].average[n].Nm - 1; k++)
    134           m = next[m];
     130
     131        /** insert star in measurement list */
     132        /* find last measurement of this star */
     133        m = catalog[0].average[n].offset; 
     134        for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m];
     135        /* set up references */
    135136        next[Nmeas] = next[m];
    136137        next[m] = Nmeas;
    137         if (next[Nmeas] == -1) { /* last just was moved */
    138           last = Nmeas;
    139         }
     138        /* last just was moved */
     139        if (next[Nmeas] == -1) last = Nmeas;
    140140       
     141        /* calculate accurate per-star airmass */
     142        secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude);
     143
     144        /* find ref measurements & add to calibration table */
     145        /* is this search right?  it does not use the next[] list */
    141146        if (CALIBRATE) {
    142147          found = FALSE;
    … …  
    145150          m = catalog[0].average[n].offset;  /* first measurement of this star */
    146151          for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) {
    147             if (catalog[0].measure[m+k].source == CalReference) { found0 = TRUE; CalM = catalog[0].measure[m+k].M; }
    148             if (catalog[0].measure[m+k].source == CalColor)     { found1 = TRUE; CalC = catalog[0].measure[m+k].M; }
     152            if (catalog[0].measure[m+k].source == CalReference) {
     153              found0 = TRUE;
     154              CalM = catalog[0].measure[m+k].M;
     155              dCalM = catalog[0].measure[m+k].dM;
     156            }
     157            if (catalog[0].measure[m+k].source == CalColor) {
     158              found1 = TRUE;
     159              CalC = catalog[0].measure[m+k].M;
     160            }
    149161            if (found0 && found1) {
    150162              found = TRUE;
    151               SaveCalibration (MIN (NO_MAG, 1000*stars[N].M + MTIME),
    152                                MIN (255, stars[N].dM),
    153                                CalM, CalC,
    154                                image[0].secz, N,
    155                                stars[N].R, stars[N].D, Xs[n], Ys[n]);
     163              SaveCalibration (MIN (NO_MAG, 1000*stars[0][N].M + MTIME),
     164                               MIN (255, stars[0][N].dM),
     165                               CalM, dCalM, CalC,
     166                               secz, N,
     167                               stars[0][N].R, stars[0][N].D, Xs[n], Ys[n]);
    156168            }
    157169          }
    … …  
    159171             
    160172        /** add measurements for this star **/
    161         catalog[0].measure[Nmeas].dR     = 360000.0*(catalog[0].average[n].R - stars[N].R);
    162         catalog[0].measure[Nmeas].dD     = 360000.0*(catalog[0].average[n].D - stars[N].D);
    163         catalog[0].measure[Nmeas].M      = MIN (NO_MAG, 1000*stars[N].M + MTIME);
     173        catalog[0].measure[Nmeas].dR     = 360000.0*(catalog[0].average[n].R - stars[0][N].R);
     174        catalog[0].measure[Nmeas].dD     = 360000.0*(catalog[0].average[n].D - stars[0][N].D);
     175        catalog[0].measure[Nmeas].M      = MIN (NO_MAG, 1000*stars[0][N].M + MTIME);
    164176        catalog[0].measure[Nmeas].Mcal   = image[0].Mcal;
    165         catalog[0].measure[Nmeas].dM     = MIN (255, stars[N].dM);  /* error in input files stored in thousandths of mag */
    166         catalog[0].measure[Nmeas].t      = image[0].tzero + 1e-4*stars[N].Y*image[0].trate;  /* trate is in 0.1 msec / row */
     177        catalog[0].measure[Nmeas].dM     = MIN (255, stars[0][N].dM);  /* error in input files stored in thousandths of mag */
     178        catalog[0].measure[Nmeas].t      = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate;  /* trate is in 0.1 msec / row */
    167179        catalog[0].measure[Nmeas].averef = n;
    168180        catalog[0].measure[Nmeas].source = image[0].source;  /* photometry source */
    169         catalog[0].measure[Nmeas].dophot = stars[N].dophot; 
     181        catalog[0].measure[Nmeas].dophot = stars[0][N].dophot; 
    170182        catalog[0].measure[Nmeas].flags  = 0;
    171183        catalog[0].measure[Nmeas].dt     = MTIME;
    172184
    173         catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000*stars[N].Mgal + MTIME);
    174         catalog[0].measure[Nmeas].airmass  = image[0].secz;
    175         catalog[0].measure[Nmeas].FWx      = MIN (NO_MAG, 100*stars[N].fx);
    176         catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[N].fy / stars[N].fx));
    177         catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[N].df);
     185        catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000*stars[0][N].Mgal + MTIME);
     186        catalog[0].measure[Nmeas].airmass  = secz;
     187        catalog[0].measure[Nmeas].FWx      = MIN (NO_MAG, 100*stars[0][N].fx);
     188        catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx));
     189        catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[0][N].df);
    178190        /* refers to same number as first measurement */
    179191       
    … …  
    185197        /*** handle multiple stars */
    186198        /* this image star matches more than one catalog star */
    187         if (stars[N].found > -1) {
    188           catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE;
     199        if (stars[0][N].found > -1) {
     200          catalog[0].measure[stars[0][N].found].flags |= BLEND_IMAGE;
    189201          catalog[0].measure[Nmeas].flags |= BLEND_IMAGE;
    190202        }
    191         if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
     203        if (stars[0][N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
    192204          catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR;
    193205        }
    194         if (stars[N].found == -1) { /* this image star matches only this catalog star */
    195           stars[N].found = Nmeas;  /* save first match, in case coincidences are found */
     206        if (stars[0][N].found == -1) { /* this image star matches only this catalog star */
     207          stars[0][N].found = Nmeas;  /* save first match, in case coincidences are found */
    196208        }
    197209        /* this catalog star matches more than one image star */
    … …  
    291303  for (i = 0; i < Nstars; i++) {
    292304    N = N1[i];
    293     if ((stars[N].found < 0) && IN_CATALOG (stars[N].R, stars[N].D)) {
    294       catalog[0].average[Nave].R = stars[N].R;
    295       catalog[0].average[Nave].D = stars[N].D;
     305    if ((stars[0][N].found < 0) && IN_CATALOG (stars[0][N].R, stars[0][N].D)) {
     306      catalog[0].average[Nave].R = stars[0][N].R;
     307      catalog[0].average[Nave].D = stars[0][N].D;
    296308      catalog[0].average[Nave].M = NO_MAG;
    297309      for (j = 0; j < Nsecfilt; j++) {
    … …  
    300312      }
    301313     
    302       Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K;
    303       Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M;
    304       if (*Mval == NO_MAG) *Mval = Mrel;
     314      if (SKYPROBE)
     315        secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude);
     316      else
     317        secz = image[0].secz;
    305318
    306319      catalog[0].average[Nave].Nm        = 1;
    … …  
    315328      catalog[0].measure[Nmeas].dR       = 0.0;
    316329      catalog[0].measure[Nmeas].dD       = 0.0;
    317       catalog[0].measure[Nmeas].M        = MIN (NO_MAG, 1000.0*stars[N].M + MTIME);
     330      catalog[0].measure[Nmeas].M        = MIN (NO_MAG, 1000.0*stars[0][N].M + MTIME);
    318331      catalog[0].measure[Nmeas].Mcal     = image[0].Mcal;
    319       catalog[0].measure[Nmeas].dM       = MIN (255, stars[N].dM);
    320       catalog[0].measure[Nmeas].t        = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */
     332      catalog[0].measure[Nmeas].dM       = MIN (255, stars[0][N].dM);
     333      catalog[0].measure[Nmeas].t        = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate; /* trate is in 0.1 msec / row */
    321334      catalog[0].measure[Nmeas].averef   = Nave;
    322335      catalog[0].measure[Nmeas].source   = image[0].source;  /* photometry source */
    323       catalog[0].measure[Nmeas].dophot   = stars[N].dophot; 
     336      catalog[0].measure[Nmeas].dophot   = stars[0][N].dophot; 
    324337      catalog[0].measure[Nmeas].flags    = 0;
    325338      catalog[0].measure[Nmeas].dt       = MTIME;
    326339
    327       catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000.0*stars[N].Mgal + MTIME);
    328       catalog[0].measure[Nmeas].airmass  = image[0].secz;
    329       catalog[0].measure[Nmeas].FWx      = MIN (0x7fff, 100*stars[N].fx);
    330       catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[N].fy / stars[N].fx));
    331       catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[N].df);
    332       stars[N].found = Nmeas;
     340      catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000.0*stars[0][N].Mgal + MTIME);
     341      catalog[0].measure[Nmeas].airmass  = secz;
     342      catalog[0].measure[Nmeas].FWx      = MIN (0x7fff, 100*stars[0][N].fx);
     343      catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx));
     344      catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[0][N].df);
     345
     346      Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K;
     347      Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M;
     348      if (*Mval == NO_MAG) *Mval = Mrel;
     349
     350      stars[0][N].found = Nmeas;
    333351      next[last] = Nmeas;
    334352      next[Nmeas] = -1;
    … …  
    340358        REALLOCATE (catalog[0].measure, Measure, NMEAS);
    341359      }
     360
    342361      /** now add references from all previous non-detection observations of this spot on the sky */
    343362      for (j = 0; (j < Nimage) && !SKIP_MISSED; j++) {
    … …  
    408427  /* note stars which have been found in this catalog */
    409428  for (i = 0; i < Nstars; i++) {
    410     if (stars[i].found > -1) {
    411       stars[i].found = -2;
     429    if (stars[0][i].found > -1) {
     430      stars[0][i].found = -2;
    412431    }
    413432  }

  • trunk/Ohana/src/addstar/src/gimages.c

    r2457 r3339  
    11# include "addstar.h"
    2 double opening_angle ();
    32
    4 Image *gimages (FILE *f, Image *image, int *Npimage) {
     3/* given image, find catalog images which overlap it */
     4Image *gimages (Image *image, int *Npimage) {
    55 
    66  int i, j, k, Nimage, addtolist, size;
    … …  
    1515  struct stat filestatus;
    1616 
     17  f = GetDB ();
    1718  fseek (f, 0, SEEK_SET);
    1819
    1920  /* read header */
    2021  if (!fits_load_header (f, &header)) {
    21     if (VERBOSE) fprintf (stderr, "ERROR: can't read image catalog %s\n", ImageCat);
    22     exit (1);
     22    if (VERBOSE) fprintf (stderr, "error: can't read image catalog %s\n", ImageCat);
     23    Shutdown ();
    2324  }
    2425
    2526  fits_scan (&header, "ZERO_PT", "%lf", 1, &zeropt);
    2627  if (fabs (ZeroPt - zeropt) > 1e-4) {
    27     fprintf (stderr, "ERROR: zero point in image (%f:%s) inconsistent with zero point in catalog (%f)\n",
     28    fprintf (stderr, "error: zero point in image (%f:%s) inconsistent with zero point in catalog (%f)\n",
    2829             ZeroPt, image[0].name, zeropt);
    29     exit (1);
     30    Shutdown ();
    3031  }
    3132
    … …  
    5657  /* check that file size makes sense */
    5758  if (stat (ImageCat, &filestatus) == -1) {
    58     fprintf (stderr, "ERROR: failed to get status of image catalog\n");
    59     exit (1);
     59    fprintf (stderr, "error: failed to get status of image catalog\n");
     60    Shutdown ();
    6061  }
    6162  size = Nimage*sizeof(Image) + header.size;
    6263  if (size != filestatus.st_size) {
    6364    int Ndata;
     65
    6466    Ndata = (filestatus.st_size - header.size) / sizeof (Image);
    65     fprintf (stderr, "ERROR: image catalog has inconsistent size\n");
    66     fprintf (stderr, "header: %d, data: %d\n", Nimage, Ndata);
     67    fprintf (stderr, "error: image catalog has inconsistent size\n");
     68    fprintf (stderr, "header: %d images, data: %d images\n", Nimage, Ndata);
    6769    if (!FORCE_READ) Shutdown ();
    68     Nimage = Ndata;
     70
     71    Nimage = MIN (Nimage, Ndata);
     72    fprintf (stderr, "using %d images\n", Nimage);
    6973  }
    7074
    … …  
    7680  ALLOCATE (pimage, Image, NPIMAGE);
    7781
     82  /* run through image table, reading images in blocks */
    7883  for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) {
    7984    ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image");
    … …  
    9297      }
    9398      /* check if any edges cross */
     99      /* not every robust - images must be close in size */
    94100      for (j = 0; (j < 4) && !addtolist; j++) {
    95101        for (k = 0; (k < 4) && !addtolist; k++) {
    … …  
    108114
    109115    if (ntimage == 0) {
    110       fprintf (stderr, "ERROR: expected %d images, only found %d\n", Nimage, Ntimage);
    111       exit (1);
     116      fprintf (stderr, "error: expected %d images, only found %d\n", Nimage, Ntimage);
     117      Shutdown ();
    112118    }
    113119  }

  • trunk/Ohana/src/addstar/src/gstars.c

    r2457 r3339  
    1616  double tmp;
    1717  int done, doneread, Nskip, Nextra, itmp;
     18  int hour, min;
     19  double sec, dMs;
     20  char line[80];
    1821
    1922  if (!fits_read_header (file, &header)) {
    … …  
    5255  }
    5356   
     57  /* get ST (used for airmass calculation) */
     58  fits_scan (&header, STKeyword, "%s", 1, line);
     59  /* remove ':' characters */
     60  for (c = strchr (line, ':'); c != (char *) NULL; c = strchr (line, ':')) { *c = ' '; }
     61  sscanf (line, "%d %d %lf", &hour, &min, &sec);
     62  SiderealTime = hour + min/60.0 + sec/3600.0;
     63
    5464  /* CERROR in data file is in arcsec */
    5565  if (!fits_scan (&header, "CERROR",   "%lf", 1, &tmp)) tmp = 1.0;
    … …  
    120130  /* secz is in units milli-airmass */
    121131  image[0].Mcal = 0.0;
     132  image[0].Xm   = NO_MAG;
    122133  image[0].code = 0;
    123   image[0].Xm   = NO_MAG;
    124134  bzero (image[0].dummy, sizeof(image[0].dummy));
    125135
    … …  
    131141  }
    132142  ALLOCATE (stars, Stars, image[0].nstar);
     143
     144  /* temporary:  get sky background from Flips data */
     145  if (SKYPROBE) {
     146    char *p;
     147    int sky;
     148   
     149    sky = 0;
     150    fits_scan (&header, "HISTORY", "%S", 10, line);
     151    p = strstr (line, "|I=");
     152    if (p != (char *) NULL) {
     153      p += 3;
     154      sky = atoi (p);
     155    }
     156    image[0].Myyyy = MIN (sky - 0x8000, 0xffff);
     157  }
    133158
    134159  /* load in stars by blocks of 1000 */
    … …  
    200225      while (stars[N].R >= 360.0) stars[N].R -= 360.0;
    201226      stars[N].found = -1;
     227
     228      if (SKYPROBE) {
     229        dMs = get_subpix (stars[N].X, stars[N].Y);
     230        stars[N].M    -= dMs;
     231        stars[N].Mgal -= dMs;
     232        stars[N].Map  -= dMs;
     233        dMs = 1000.0*scat_subpix (stars[N].X, stars[N].Y);
     234        stars[N].dM = hypot (stars[N].dM, dMs);
     235      }
    202236    }
    203237  }
    … …  
    210244  return (stars);
    211245}
    212 
    213 
    214 
    215246
    216247# if (0)

  • trunk/Ohana/src/addstar/src/parse_time.c

    r2457 r3339  
    7575  for (c = strchr (line, 0x2d); c != (char *) NULL; c = strchr (line, 0x2d)) { *c = ' '; }
    7676
     77  Nf = 0;
    7778  switch (mode) {
    7879  case 1:
    … …  
    106107  }   
    107108
     109  /* this should probably use localtime */
     110
    108111  /* convert yy.mm.dd hh.mm.ss to Nsec since 1970 (jd = 2440587.5) */
    109112  /* note that in this section, tm_mon has range 1-12, unlike for gmtime () */

  • trunk/Ohana/src/addstar/src/sort_lists.c

    r2457 r3339  
     1
     2void sort_GSCregions (GSCRegion *region, int Nregion) {
     3
     4  int l,j,ir,i;
     5  GSCRegion tmp;
     6 
     7  l = N >> 1;
     8  ir = N - 1;
     9  for (;;) {
     10    if (l > 0) {
     11      l--;
     12      tmp = region[l];
     13    }
     14    else {
     15      tmp = region[ir];
     16      region[ir] = region[0];
     17      if (--ir == 0) {
     18        region[0] = tmp;
     19        return;
     20      }
     21    }
     22    i = l;
     23    j = (l << 1) + 1;
     24    while (j <= ir) {
     25      if (j < ir && region[j].DEC < region[j+1].DEC) j++;
     26      if (tmp.DEC < region[j].DEC) {
     27        region[i] = region[j];
     28        j += (i=j) + 1;
     29      }
     30      else j = ir + 1;
     31    }
     32    region[i] = tmp;
     33  }
     34}
    135
    236void sort_lists (float *X, float *Y, int *S, int N) {
    … …  
    79113  }
    80114}
     115
     116void fsort2 (float *X, float *Y, int N) {
     117
     118  int l,j,ir,i;
     119  float tX, tY;
     120 
     121  l = N >> 1;
     122  ir = N - 1;
     123  for (;;) {
     124    if (l > 0) {
     125      l--;
     126      tX = X[l];
     127      tY = Y[l];
     128    }
     129    else {
     130      tX = X[ir];
     131      tY = Y[ir];
     132      X[ir] = X[0];
     133      Y[ir] = Y[0];
     134      if (--ir == 0) {
     135        X[0] = tX;
     136        Y[0] = tY;
     137        return;
     138      }
     139    }
     140    i = l;
     141    j = (l << 1) + 1;
     142    while (j <= ir) {
     143      if (j < ir && X[j] < X[j+1]) j++;
     144      if (tX < X[j]) {
     145        X[i] = X[j];
     146        Y[i] = Y[j];
     147        j += (i=j) + 1;
     148      }
     149      else j = ir + 1;
     150    }
     151    X[i] = tX;
     152    Y[i] = tY;
     153  }
     154}

  • trunk/Ohana/src/addstar/src/wimage.c

    r2457 r3339  
    33void wimage (FILE *f, int dbstate, Image *image, int Nstars) {
    44 
    5   int Nimages, status;
     5  int Nimages, status, offset;
    66  Header header;
    77
    … …  
    3737  }
    3838
    39   /* position to end of file for new image data */
    40   fseek (f, 0, SEEK_END);
     39  /* position to end of data array */
     40  offset = (Nimages - 1)*sizeof(Image) + header.size;
     41  fseek (f, offset, SEEK_SET);
     42
    4143  status = Fwrite (image, sizeof(Image), 1, f, "image");
    4244  if (status != 1) {
Note: See TracChangeset for help on using the changeset viewer.