Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/Makefile
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/Makefile (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/Makefile (revision 10360)
@@ -0,0 +1,37 @@
+default: uniphot
+help:
+ @echo "make options: uniphot (default)"
+
+include ../../Configure
+HOME = $(ROOT)/src/uniphot
+BIN = $(HOME)/bin
+LIB = $(HOME)/lib
+SRC = $(HOME)/src
+MAN = $(HOME)/doc
+INC = $(HOME)/include
+include ../../Makefile.Common
+
+LDFLAGS := -ldvo -lFITS -lohana $(LDFLAGS)
+
+uniphot: $(BIN)/uniphot.$(ARCH)
+install: $(DESTBIN)/uniphot
+
+UNIPHOT = \
+$(SRC)/uniphot.$(ARCH).o \
+$(SRC)/initialize.$(ARCH).o \
+$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/args.$(ARCH).o \
+$(SRC)/liststats.$(ARCH).o \
+$(SRC)/load_images.$(ARCH).o \
+$(SRC)/subset_images.$(ARCH).o \
+$(SRC)/find_image_tgroups.$(ARCH).o \
+$(SRC)/find_image_sgroups.$(ARCH).o \
+$(SRC)/fit_groups.$(ARCH).o \
+$(SRC)/update.$(ARCH).o \
+$(SRC)/update_catalog.$(ARCH).o \
+$(SRC)/SetSignals.$(ARCH).o \
+$(SRC)/Shutdown.$(ARCH).o \
+$(SRC)/sort.$(ARCH).o
+
+$(UNIPHOT): $(INC)/uniphot.h
+$(BIN)/uniphot.$(ARCH): $(UNIPHOT)
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/bin/.cvsignore
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/bin/.cvsignore (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/bin/.cvsignore (revision 10360)
@@ -0,0 +1,2 @@
+*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix
+*.linrh
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/ChangeLog.txt
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/ChangeLog.txt (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/ChangeLog.txt (revision 10360)
@@ -0,0 +1,24 @@
+
+- uniphot-1-4:
+ * converted to gfits APIs (forces libfits 1.6)
+ * converted to new DVO APIs (forces libdvo 1.3)
+ * removed old, unused code
+
+- uniphot-1-3:
+ * dropped IMAGE_CATALOG from config
+ * convert to dvo_image_lock,unlock
+ * dropped _PS from average.R,D,M,dM
+
+- uniphot-1-2:
+
+ * various changes to use mode/format concepts
+ * mods to use new dvo load functions / libdvo
+ * minor changes to use new libohana (v1.5) / libfits (v1.4)
+
+- uniphot-1-1:
+
+ * changes from loneos.h to dvo.h (libohana)
+
+- uniphot-1-0:
+
+ * initial release
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/notes.txt
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/notes.txt (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/doc/notes.txt (revision 10360)
@@ -0,0 +1,14 @@
+
+ uniphot
+
+
+- load images (select on photcode, time range)
+
- create time groups
+
- create space groups
+
- solve for:
+ tgroup[i].M = ave(image[tgroup].Mcal - image[tgroup].sgroup.Mcal)
+ sgroup[i].M = ave(image[sgroup].Mcal - image[tgroup].tgroup.Mcal)
+
+ image.tgroup -> tgroup
+ group.Nimage, image
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/include/uniphot.h
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/include/uniphot.h (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/include/uniphot.h (revision 10360)
@@ -0,0 +1,105 @@
+# include
+# include
+# include
+
+# ifdef ANSI
+# define F_SETFL 4
+# define O_NONBLOCK 0200000
+# define AF_UNIX 1
+# define SOCK_STREAM 1
+# define ENOENT 2 /* No such file or directory */
+# endif /* ANSI */
+
+typedef struct {
+ double median;
+ double mean;
+ double sigma;
+ double error;
+ double chisq;
+ double min;
+ double max;
+ double total;
+ int Nmeas;
+} StatType;
+
+typedef struct {
+ void *tgroup;
+ void *sgroup;
+} ImageLink;
+
+typedef struct {
+ char label[64];
+ float M;
+ float dM;
+ float dMsub;
+ double v1, v2;
+ Image **image;
+ ImageLink **imlink;
+ int Nimage, Ngood;
+} Group;
+
+/* global variables set in parameter file */
+char ImageCat[256];
+char CATDIR[256];
+char CATMODE[16]; /* raw, mef, split, mysql */
+char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */
+char STATMODE[64];
+int VERBOSE;
+int NLOOP;
+int TimeSelect;
+int VERBOSE;
+int UPDATE;
+int IMAGE_BAD;
+double RADIUS;
+double TRANGE;
+time_t TSTART;
+time_t TSTOP;
+PhotCode *photcode;
+
+enum {black, white, red, orange, yellow, green, blue, indigo, violet};
+
+typedef struct {
+ double xmin, xmax, ymin, ymax;
+ int style, ptype, ltype, etype, color;
+ double lweight, size;
+} Graphdata;
+
+/***** prototypes ****/
+void ConfigInit PROTO((int *argc, char **argv));
+void DonePlotting PROTO((Graphdata *graphmode, int N));
+void JpegPlot PROTO((Graphdata *graphmode, int N, char *filename));
+void PSPlot PROTO((Graphdata *graphmode, int N, char *filename));
+void PlotLabel PROTO((char *string, int N));
+void PlotVector PROTO((int Npts, double *vect, int mode, int N));
+void PrepPlotting PROTO((int Npts, Graphdata *graphmode, int N));
+void XDead PROTO(());
+int args PROTO((int argc, char **argv));
+void dumpresult PROTO(());
+Group *find_image_sgroups PROTO((FITS_DB *db, ImageLink **imlink, int *Nsgroup));
+Group *find_image_tgroups PROTO((FITS_DB *db, ImageLink **imlink, int *Ntgroup));
+void fit_sgroup PROTO((Group *sgroup, int Nsgroup));
+void fit_tgroup PROTO((Group *tgroup, int Ntgroup));
+int gcatalog PROTO((Catalog *catalog));
+void initialize PROTO((int argc, char **argv));
+void initstats PROTO((char *mode));
+int liststats PROTO((double *value, double *dvalue, int N, StatType *stats));
+int load_images PROTO((FITS_DB *db));
+int main PROTO((int argc, char **argv));
+int open_graph PROTO((int N));
+void sort PROTO((unsigned int *X, int N));
+void sortB PROTO((double *X, double *Y, int N));
+void sortD PROTO((double *X, double *Y, double *Z, int N));
+void update_dvo_catalog PROTO((Catalog *catalog, Group *sgroup, int warn));
+void wcatalog PROTO((Catalog *catalog));
+void wimages PROTO((Image *image, int Nimage));
+void check_permissions (char *basefile);
+void lock_image_db (FITS_DB *db, char *filename);
+void unlock_image_db (FITS_DB *db);
+void create_image_db (FITS_DB *db);
+void set_db (FITS_DB *in);
+int Shutdown (char *format, ...);
+void TrapSignal (int sig);
+void SetProtect (int mode);
+int SetSignals ();
+int subset_images (FITS_DB *db);
+void update (FITS_DB *db, Group *sgroup, int Nsgroup);
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/ConfigInit.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/ConfigInit.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/ConfigInit.c (revision 10360)
@@ -0,0 +1,39 @@
+# include "uniphot.h"
+
+void ConfigInit (int *argc, char **argv) {
+
+ double ZERO_POINT;
+ char *config, *file;
+ char PhotCodeFile[256];
+
+ /*** load configuration info ***/
+ file = SelectConfigFile (argc, argv, "ptolemy");
+ config = LoadConfigFile (file);
+ if (config == (char *) NULL) {
+ fprintf (stderr, "ERROR: can't find configuration file %s\n", file);
+ if (file != (char *) NULL) free (file);
+ exit (0);
+ }
+ if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
+
+ ScanConfig (config, "CATDIR", "%s", 0, CATDIR);
+ ScanConfig (config, "CATMODE", "%s", 0, CATMODE);
+ ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT);
+ sprintf (ImageCat, "%s/Images.dat", CATDIR);
+
+ ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile);
+ ScanConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT);
+ SetZeroPoint (ZERO_POINT);
+
+ free (config);
+ free (file);
+
+ if (*CATMODE == 0) strcpy (CATMODE, "RAW");
+ if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
+
+ if (!LoadPhotcodes (PhotCodeFile)) {
+ fprintf (stderr, "error loading photcodes\n");
+ exit (0);
+ }
+
+}
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/SetSignals.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/SetSignals.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/SetSignals.c (revision 10360)
@@ -0,0 +1,118 @@
+# include "uniphot.h"
+
+static int Protect = FALSE;
+static int Trapped = FALSE;
+
+void TrapSignal (int sig) {
+ fprintf (stderr, "trapped signal %d\n", sig);
+ if (sig == 11) {
+ fprintf (stderr, "seg fault\n");
+ exit (1);
+ }
+ if (Protect) {
+ Trapped = TRUE;
+ fprintf (stderr, "blocking until protected sections are clear\n");
+ return;
+ }
+ Shutdown ("halted by signal (trapped)");
+}
+
+void SetProtect (int mode) {
+ Protect = mode;
+ if (Trapped && !Protect) Shutdown ("halted by signal (protect)");
+}
+
+int SetSignals () {
+
+ int i;
+
+ /* disable almost all signal interrupts */
+ for (i = 0; i < 36; i++) {
+ switch (i) {
+ /* can't redirect this signals */
+ case SIGKILL: /* kill -9: cannot be caught or ignored */
+ case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */
+ /* ignore these signals */
+ case SIGCHLD: /* child halted: ignore */
+ case SIGPWR: /* power failure - why ignore this? */
+ case SIGWINCH: /* window resized */
+ case SIGCONT: /* continue - maintain this action */
+ case SIGTSTP: /* stop signal sent from tty - why ignore? */
+ case SIGURG: /* socket signal, ignore this */
+ break;
+
+ default:
+ signal (i, TrapSignal);
+ }
+ }
+ return (TRUE);
+}
+/*
+
+ Signal Value Action Comment
+ -------------------------------------------------------------------------
+ SIGHUP 1 A Hangup detected on controlling terminal
+ or death of controlling process
+ SIGINT 2 A Interrupt from keyboard
+ SIGQUIT 3 A Quit from keyboard
+ SIGILL 4 A Illegal Instruction
+ SIGABRT 6 C Abort signal from abort(3)
+ SIGFPE 8 C Floating point exception
+ SIGKILL 9 AEF Kill signal
+ SIGSEGV 11 C Invalid memory reference
+ SIGPIPE 13 A Broken pipe: write to pipe with no readers
+ SIGALRM 14 A Timer signal from alarm(2)
+ SIGTERM 15 A Termination signal
+ SIGUSR1 30,10,16 A User-defined signal 1
+ SIGUSR2 31,12,17 A User-defined signal 2
+ SIGCHLD 20,17,18 B Child stopped or terminated
+ SIGCONT 19,18,25 Continue if stopped
+ SIGSTOP 17,19,23 DEF Stop process
+ SIGTSTP 18,20,24 D Stop typed at tty
+ SIGTTIN 21,21,26 D tty input for background process
+ SIGTTOU 22,22,27 D tty output for background process
+
+ Next various other signals.
+
+ Signal Value Action Comment
+ ---------------------------------------------------------------------
+ SIGTRAP 5 CG Trace/breakpoint trap
+ SIGIOT 6 CG IOT trap. A synonym for SIGABRT
+ SIGEMT 7,-,7 G
+ SIGBUS 10,7,10 AG Bus error
+ SIGSYS 12,-,12 G Bad argument to routine (SVID)
+ SIGSTKFLT -,16,- AG Stack fault on coprocessor
+ SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD)
+ SIGIO 23,29,22 AG I/O now possible (4.2 BSD)
+ SIGPOLL AG A synonym for SIGIO (System V)
+ SIGCLD -,-,18 G A synonym for SIGCHLD
+ SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD)
+ SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD)
+ SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD)
+ SIGPROF 27,27,29 AG Profile alarm clock
+ SIGPWR 29,30,19 AG Power failure (System V)
+ SIGINFO 29,-,- G A synonym for SIGPWR
+ SIGLOST -,-,- AG File lock lost
+ SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun)
+ SIGUNUSED -,31,- AG Unused signal
+ (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and
+ 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
+ sparc.)
+
+ The letters in the "Action" column have the following meanings:
+
+ A Default action is to terminate the process.
+
+ B Default action is to ignore the signal.
+
+ C Default action is to dump core.
+
+ D Default action is to stop the process.
+
+ E Signal cannot be caught.
+
+ F Signal cannot be ignored.
+
+ G Not a POSIX.1 conformant signal.
+
+*/
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/Shutdown.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/Shutdown.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/Shutdown.c (revision 10360)
@@ -0,0 +1,28 @@
+# include "uniphot.h"
+
+static FITS_DB *db;
+
+void set_db (FITS_DB *in) {
+ db = in;
+}
+
+/* clean up open / locked ImageCat before shutting down */
+int Shutdown (char *format, ...) {
+ va_list argp;
+ char *formatplus;
+
+ ALLOCATE (formatplus, char, strlen(format));
+ strcpy (formatplus, format);
+ strcat (formatplus, "\n");
+
+ va_start (argp, format);
+ vfprintf (stderr, formatplus, argp);
+ free (formatplus);
+ va_end (argp);
+
+ SetProtect (TRUE);
+ gfits_db_close (db);
+ fprintf (stderr, "ERROR: addstar halted\n");
+ exit (1);
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/args.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/args.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/args.c (revision 10360)
@@ -0,0 +1,73 @@
+# include "uniphot.h"
+
+int args (int argc, char **argv) {
+
+ int N;
+
+ /* define time */
+ TimeSelect = FALSE;
+ if ((N = get_argument (argc, argv, "-time"))) {
+ remove_argument (N, &argc, argv);
+ if (!str_to_time (argv[N], &TSTART)) {
+ fprintf (stderr, "ERROR: syntax error\n");
+ exit (1);
+ }
+ remove_argument (N, &argc, argv);
+ if (!str_to_time (argv[N], &TSTOP)) {
+ fprintf (stderr, "ERROR: syntax error\n");
+ exit (1);
+ }
+ remove_argument (N, &argc, argv);
+ TimeSelect = TRUE;
+ }
+
+ VERBOSE = FALSE;
+ if ((N = get_argument (argc, argv, "-v"))) {
+ VERBOSE = TRUE;
+ remove_argument (N, &argc, argv);
+ }
+
+ NLOOP = 8;
+ if ((N = get_argument (argc, argv, "-n"))) {
+ remove_argument (N, &argc, argv);
+ NLOOP = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ strcpy (STATMODE, "MEAN");
+ if ((N = get_argument (argc, argv, "-statmode"))) {
+ remove_argument (N, &argc, argv);
+ strcpy (STATMODE, argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ UPDATE = FALSE;
+ if ((N = get_argument (argc, argv, "-update"))) {
+ remove_argument (N, &argc, argv);
+ UPDATE = TRUE;
+ }
+
+ /* max separation for unique space group, in degrees */
+ RADIUS = 2.0;
+ if ((N = get_argument (argc, argv, "-radius"))) {
+ remove_argument (N, &argc, argv);
+ RADIUS = atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ /* max separation for unique time group, in days -> seconds */
+ TRANGE = 86400*7.0;
+ if ((N = get_argument (argc, argv, "-trange"))) {
+ remove_argument (N, &argc, argv);
+ TRANGE = 86400*atof (argv[N]);
+ remove_argument (N, &argc, argv);
+ }
+
+ if (argc != 2) {
+ fprintf (stderr, "ERROR: USAGE: uniphot (photcode) [options]\n");
+ exit (2);
+ }
+
+ return (TRUE);
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/dumpresult.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/dumpresult.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/dumpresult.c (revision 10360)
@@ -0,0 +1,25 @@
+# include "uniphot.h"
+
+void dumpresult () {
+ short Mcal, Mgrp, Mset;
+ FILE *f;
+ char outfile[64];
+ Group *tgrp;
+
+ for (i = 0; i < Nsgroup; i++) {
+ sprintf (outfile, "test.%02d.dat", i);
+ f = fopen (outfile, "w");
+ for (j = 0; j < sgroup[i].Nimage; j++) {
+ if (sgroup[i].image[j][0].code & IMAGE_BAD) continue;
+ tgrp = (Group *) sgroup[i].imlink[j][0].tgroup;
+ Mcal = sgroup[i].image[j][0].Mcal;
+ Mset = sgroup[i].M;
+ Mgrp = tgrp[0].M;
+ fprintf (f, "%7.4f %7.4f %7.4f %7.4f %10.6f %10.6f %f %s\n",
+ 0.001*Mcal, 0.001*Mgrp, 0.001*Mset, 0.001*sgroup[i].image[j][0].dMcal,
+ sgroup[i].image[j][0].coords.crval1, sgroup[i].image[j][0].coords.crval2, (sgroup[i].image[j][0].tzero-915148800)/86400.0, tgrp[0].label);
+ }
+ fclose (f);
+ }
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_sgroups.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_sgroups.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_sgroups.c (revision 10360)
@@ -0,0 +1,85 @@
+# include "uniphot.h"
+
+Group *find_image_sgroups (FITS_DB *db, ImageLink **Imlink, int *Nsgroup) {
+
+ int i, j, Nimage, Ngroup, Nentry, NENTRY;
+ double r, d, x, y, radius;
+ Group *group;
+ Coords coords;
+ Image *image;
+ ImageLink *imlink;
+
+ imlink = *Imlink;
+ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+
+ coords.crpix1 = coords.crpix2 = 0.0;
+ coords.cdelt1 = coords.cdelt2 = 1.0; /* scale is degrees, radius in degrees */
+ coords.pc1_1 = coords.pc2_2 = 1.0;
+ coords.pc1_2 = coords.pc2_1 = 0.0;
+ strcpy (coords.ctype, "RA---TAN");
+
+ Ngroup = 0;
+ ALLOCATE (group, Group, Nimage);
+
+ /* set imlink.sgroups = NULL as a marker */
+ for (i = 0; i < Nimage; i++) imlink[i].sgroup = NULL;
+
+ for (i = 0; i < Nimage; i++) {
+ if (imlink[i].sgroup != NULL) continue;
+ if (image[i].code & ID_IMAGE_NOCAL) continue;
+
+ /* new sgroup, set ref coords */
+ XY_to_RD (&r, &d, 0.5*image[i].NX, 0.5*image[i].NX, &image[i].coords);
+ coords.crval1 = r;
+ coords.crval2 = d;
+
+ /* init sgroup structure */
+ Nentry = 0;
+ NENTRY = 100;
+ ALLOCATE (group[Ngroup].image, Image *, NENTRY);
+ ALLOCATE (group[Ngroup].imlink, ImageLink *, NENTRY);
+ group[Ngroup].M = 0;
+ group[Ngroup].dM = 0;
+ snprintf (group[Ngroup].label, 64, "%10.6f - %10.6f", r, d);
+ group[Ngroup].v1 = r;
+ group[Ngroup].v2 = d;
+
+ /* link this image to sgroup */
+ group[Ngroup].image[Nentry] = &image[i];
+ group[Ngroup].imlink[Nentry] = &imlink[i];
+ imlink[i].sgroup = &group[Ngroup];
+ Nentry ++;
+
+ for (j = 0; j < Nimage; j++) {
+ if (image[j].code & ID_IMAGE_NOCAL) continue;
+ if (imlink[j].sgroup != NULL) continue;
+
+ /* project image center to local coords, check radius */
+ XY_to_RD (&r, &d, 0.5*image[j].NX, 0.5*image[j].NX, &image[j].coords);
+ if (!RD_to_XY (&x, &y, r, d, &coords)) continue;
+ /* RD_to_XY returns FALSE if opposite hemispheres */
+ radius = hypot (x, y);
+
+ if (radius > RADIUS) continue;
+
+ /* image in sgroup, add to entry */
+ group[Ngroup].image[Nentry] = &image[j];
+ group[Ngroup].imlink[Nentry] = &imlink[j];
+ imlink[j].sgroup = &group[Ngroup];
+ Nentry ++;
+ if (Nentry == NENTRY) {
+ NENTRY += 100;
+ REALLOCATE (group[Ngroup].image, Image *, NENTRY);
+ REALLOCATE (group[Ngroup].imlink, ImageLink *, NENTRY);
+ }
+ }
+ group[Ngroup].Nimage = Nentry;
+ Ngroup ++;
+ }
+ *Nsgroup = Ngroup;
+ return (group);
+}
+
+ /* this is a bit weak: since we use pointers, we can't
+ reallocate group after the pointers are assigned.
+ therefore, we allocate the max possible groups */
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_tgroups.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_tgroups.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/find_image_tgroups.c (revision 10360)
@@ -0,0 +1,84 @@
+# include "uniphot.h"
+
+Group *find_image_tgroups (FITS_DB *db, ImageLink **Imlink, int *Ntgroup) {
+
+ char *start, *stop;
+ int i, j, Nimage, Ngroup, NGROUP, Nentry, NENTRY;
+ unsigned int *time, *tmin, *tmax;
+ Group *group;
+ Image *image;
+ ImageLink *imlink;
+
+ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+ ALLOCATE (imlink, ImageLink, Nimage);
+
+ /* sort time list (use only valid images?) */
+ ALLOCATE (time, unsigned int, Nimage);
+ for (i = 0; i < Nimage; i++) {
+ time[i] = image[i].tzero;
+ }
+ sort (time, Nimage);
+
+ /* find groups with dt < TRANGE */
+ Ngroup = 0;
+ NGROUP = 100;
+ ALLOCATE (tmin, unsigned int, NGROUP);
+ ALLOCATE (tmax, unsigned int, NGROUP);
+ tmin[Ngroup] = time[0];
+
+ /* generate tgroups */
+ for (i = 0; i < Nimage - 1; i++) {
+ if (time[i+1] - time[i] < TRANGE) continue;
+
+ tmax[Ngroup] = time[i];
+
+ Ngroup ++;
+ if (Ngroup == NGROUP) {
+ NGROUP += 100;
+ REALLOCATE (tmin, unsigned int, NGROUP);
+ REALLOCATE (tmax, unsigned int, NGROUP);
+ }
+ tmin[Ngroup] = time[i + 1];
+ }
+ tmax[Ngroup] = time[Nimage - 1];
+ Ngroup ++;
+ ALLOCATE (group, Group, Ngroup);
+
+ /* assign images to groups */
+ for (i = 0; i < Ngroup; i++) {
+ Nentry = 0;
+ NENTRY = 100;
+ ALLOCATE (group[i].image, Image *, NENTRY);
+ ALLOCATE (group[i].imlink, ImageLink *, NENTRY);
+ group[i].M = 0;
+ group[i].dM = 0;
+
+ start = sec_to_date (tmin[i]);
+ stop = sec_to_date (tmax[i]);
+ snprintf (group[i].label, 64, "%s - %s", start, stop);
+ free (start);
+ free (stop);
+
+ for (j = 0; j < Nimage; j++) {
+ if (image[j].tzero < tmin[i]) continue;
+ if (image[j].tzero > tmax[i]) continue;
+ if (image[j].code & ID_IMAGE_NOCAL) continue;
+
+ group[i].image[Nentry] = &image[j];
+ group[i].imlink[Nentry] = &imlink[j];
+ imlink[j].tgroup = &group[i];
+ Nentry ++;
+ if (Nentry == NENTRY) {
+ NENTRY += 100;
+ REALLOCATE (group[i].image, Image *, NENTRY);
+ REALLOCATE (group[i].imlink, ImageLink *, NENTRY);
+ }
+ }
+ group[i].Nimage = Nentry;
+ }
+ *Imlink = imlink;
+ *Ntgroup = Ngroup;
+ return (group);
+}
+
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/fit_groups.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/fit_groups.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/fit_groups.c (revision 10360)
@@ -0,0 +1,76 @@
+# include "uniphot.h"
+
+void fit_tgroup (Group *tgroup, int Ntgroup) {
+
+ int i, j, Nlist;
+ float Mcal, Mgrp;
+ double *mlist, *dlist;
+ StatType stats;
+ Group *sgroup;
+
+ Nlist = tgroup[0].Nimage;
+ for (i = 0; i < Ntgroup; i++) { Nlist = MAX (Nlist, tgroup[i].Nimage); }
+ ALLOCATE (mlist, double, Nlist);
+ ALLOCATE (dlist, double, Nlist);
+
+ for (i = 0; i < Ntgroup; i++) {
+ for (j = Nlist = 0; j < tgroup[i].Nimage; j++) {
+ if (tgroup[i].image[j][0].code & IMAGE_BAD) continue;
+ sgroup = (Group *) tgroup[i].imlink[j][0].sgroup;
+ Mcal = tgroup[i].image[j][0].Mcal_PS;
+ Mgrp = sgroup[0].M;
+ mlist[Nlist] = (Mcal - Mgrp);
+ dlist[Nlist] = tgroup[i].image[j][0].dMcal_PS;
+ Nlist ++;
+ }
+ liststats (mlist, dlist, Nlist, &stats);
+ tgroup[i].M = stats.mean;
+ tgroup[i].dMsub = stats.sigma;
+ tgroup[i].Ngood = stats.Nmeas;
+
+ initstats ("MEAN");
+ liststats (mlist, dlist, Nlist, &stats);
+ tgroup[i].dM = stats.sigma;
+ initstats (STATMODE);
+ }
+ free (mlist);
+ free (dlist);
+}
+
+void fit_sgroup (Group *sgroup, int Nsgroup) {
+
+ int i, j, Nlist;
+ float Mcal, Mgrp;
+ double *mlist, *dlist;
+ StatType stats;
+ Group *tgroup;
+
+ Nlist = sgroup[0].Nimage;
+ for (i = 0; i < Nsgroup; i++) { Nlist = MAX (Nlist, sgroup[i].Nimage); }
+ ALLOCATE (mlist, double, Nlist);
+ ALLOCATE (dlist, double, Nlist);
+ Nlist = 0;
+
+ for (i = 0; i < Nsgroup; i++) {
+ for (j = Nlist = 0; j < sgroup[i].Nimage; j++) {
+ if (sgroup[i].image[j][0].code & IMAGE_BAD) continue;
+ tgroup = (Group *) sgroup[i].imlink[j][0].tgroup;
+ Mcal = sgroup[i].image[j][0].Mcal_PS;
+ Mgrp = tgroup[0].M;
+ mlist[Nlist] = (Mcal - Mgrp);
+ dlist[Nlist] = sgroup[i].image[j][0].dMcal_PS;
+ Nlist ++;
+ }
+ liststats (mlist, dlist, Nlist, &stats);
+ sgroup[i].M = stats.mean;
+ sgroup[i].dMsub = stats.sigma;
+ sgroup[i].Ngood = stats.Nmeas;
+
+ initstats ("MEAN");
+ liststats (mlist, dlist, Nlist, &stats);
+ sgroup[i].dM = stats.sigma;
+ initstats (STATMODE);
+ }
+ free (mlist);
+ free (dlist);
+}
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/initialize.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/initialize.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/initialize.c (revision 10360)
@@ -0,0 +1,22 @@
+# include "uniphot.h"
+
+void initialize (int argc, char **argv) {
+
+ /* are these set correctly? */
+ ConfigInit (&argc, argv);
+ args (argc, argv);
+
+ if ((photcode = GetPhotcodebyName (argv[1])) == NULL) {
+ fprintf (stderr, "ERROR: photcode not found in photcode table\n");
+ exit (1);
+ }
+ if ((photcode[0].type == PHOT_DEP) && (photcode[0].type == PHOT_REF)) {
+ fprintf (stderr, "photcode must be primary or secondary type\n");
+ exit (1);
+ }
+
+ IMAGE_BAD = ID_IMAGE_NOCAL | ID_IMAGE_POOR | ID_IMAGE_SKIP | ID_IMAGE_FEW;
+
+ initstats (STATMODE);
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/liststats.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/liststats.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/liststats.c (revision 10360)
@@ -0,0 +1,117 @@
+# include "uniphot.h"
+
+enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
+ M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
+
+static int statmode;
+
+void initstats (char *mode) {
+
+ statmode = -1;
+ if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
+ if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
+ if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN;
+ if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
+ if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
+ if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
+ if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
+
+ if (statmode == -1) {
+ fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
+ exit (1);
+ }
+}
+
+int liststats (double *value, double *dvalue, int N, StatType *stats) {
+
+ int i, ks, ke, Nm;
+ double Mo, dMo, M, dM, X2, dS, *chi;
+
+ ke = ks = dMo = 0;
+
+ stats[0].Nmeas = N;
+ stats[0].mean = 0;
+ stats[0].sigma = 0;
+ if (N < 2) return (FALSE);
+
+ sortB (value, dvalue, N);
+ stats[0].median = value[(int)(0.5*N)];
+ stats[0].min = value[0];
+ stats[0].max = value[N-1];
+
+ switch (statmode) {
+ case M_MEDIAN:
+ ks = 0;
+ ke = N;
+ Mo = stats[0].median;
+ Nm = N;
+ goto chisq;
+ break;
+ case M_MEAN:
+ case M_WT_MEAN:
+ ks = 0;
+ ke = N;
+ break;
+ case M_INNER_MEAN:
+ case M_INNER_WTMEAN:
+ case M_CHI_INNER_MEAN:
+ case M_CHI_INNER_WTMEAN:
+ ks = 0.25*N + 0.50;
+ ke = 0.75*N + 0.25;
+ if (N <= 3) {
+ ks = 0;
+ ke = N;
+ }
+ break;
+ }
+
+ if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
+ ALLOCATE (chi, double, N);
+ for (i = 0; i < N; i++) {
+ chi[i] = (value[i] - stats[0].median) / dvalue[i];
+ }
+ sortD (chi, value, dvalue, N);
+ free (chi);
+ }
+
+ /* calculating the per-star offset based on the weighted average */
+ M = dM = Nm = 0;
+ if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
+ for (i = ks; i < ke; i++) {
+ M += value[i] / SQ (dvalue[i]);
+ dM += 1.0 / SQ (dvalue[i]);
+ Nm ++;
+ }
+ Mo = M / dM;
+ dMo = sqrt (1.0 / dM);
+ } else {
+ for (i = ks; i < ke; i++) {
+ M += value[i];
+ dM += SQ (dvalue[i]);
+ Nm ++;
+ }
+ Mo = M / (double) Nm;
+ dMo = sqrt (dM / (double) Nm);
+ }
+
+ chisq:
+ /* find sigma and chisq */
+ X2 = dS = 0;
+ for (i = ks; i < ke; i++) {
+ M = SQ (value[i] - Mo);
+ dM = SQ (dvalue[i]);
+ X2 += M / dM;
+ dS += M;
+ }
+ X2 = X2 / Nm;
+ dS = sqrt (dS / Nm);
+
+ stats[0].mean = Mo;
+ stats[0].Nmeas = Nm;
+ stats[0].chisq = X2;
+ stats[0].sigma = dS;
+ stats[0].error = dMo;
+
+ return (TRUE);
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/load_images.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/load_images.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/load_images.c (revision 10360)
@@ -0,0 +1,11 @@
+# include "uniphot.h"
+
+int load_images (FITS_DB *db) {
+
+ if (VERBOSE) fprintf (stderr, "finding images\n");
+
+ /* read entire db table */
+ if (!dvo_image_load (db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db[0].filename);
+
+ return (TRUE);
+}
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/plotstuff.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/plotstuff.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/plotstuff.c (revision 10360)
@@ -0,0 +1,190 @@
+# include "uniphot.h"
+# include
+
+static int Xgraph[5] = {0,0,0,0,0};
+static int active;
+
+void XDead () {
+ signal (SIGPIPE, XDead);
+ fprintf (stderr, "kapa is dead, must restart\n");
+ Xgraph[active] = -1;
+}
+
+int open_graph (int N) {
+
+ int i, InitSocket, status, addreslen;
+ struct sockaddr_un Address;
+ char temp[100], *display_name;
+ char socket_name[100];
+ FILE *f;
+
+ active = N;
+ sprintf (socket_name, "/tmp/kapa%d", N);
+ sprintf (temp, "rm -f %s", socket_name);
+ system (temp);
+
+ strcpy (Address.sun_path, socket_name);
+ Address.sun_family = AF_UNIX;
+ InitSocket = socket (AF_UNIX, SOCK_STREAM, 0);
+ status = bind (InitSocket, &Address, sizeof (Address));
+ status = listen (InitSocket, 1);
+
+ sprintf (temp, "kapa %s &", socket_name);
+# ifndef DEBUG
+ system (temp);
+# else
+ fprintf (stderr, "start kapa, press return: %s\n", temp);
+ fscanf (stdin, "%d", &i);
+# endif
+
+ addreslen = sizeof (Address);
+ Xgraph[N] = accept (InitSocket, &Address, &addreslen);
+ if (Xgraph[N] < 0) {
+ fprintf (stderr, "error starting kapa\n");
+ return (FALSE);
+ }
+ else {
+ return (TRUE);
+ }
+
+}
+
+void DonePlotting (Graphdata *graphmode, int N) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "DBOX", 4);
+ sprintf (buffer, "%f %f %f %f", graphmode[0].xmin, graphmode[0].xmax, graphmode[0].ymin, graphmode[0].ymax);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+ sprintf (buffer, "%s %s %s", "2222", "2222", "2222");
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+}
+
+void JpegPlot (Graphdata *graphmode, int N, char *filename) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "PNGF", 4);
+ sprintf (buffer, "LEN: %11d", strlen (filename));
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], filename, strlen (filename));
+ read (Xgraph[N], buffer, 4);
+}
+
+void PSPlot (Graphdata *graphmode, int N, char *filename) {
+ char buffer[65], buffer2[65];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "PSIT", 4);
+ sprintf (buffer, "LEN: %11d", strlen (filename));
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], filename, strlen (filename));
+ read (Xgraph[N], buffer, 4);
+}
+
+void PrepPlotting (int Npts, Graphdata *graphmode, int N) {
+
+ int i, status;
+ char buffer[128], buffer2[128];
+
+ active = N;
+ if (Npts < 1) return;
+
+ if (Xgraph[N] < 1) if (!open_graph(N)) return;
+
+ /* test Xgraph[N], flush junk from pipe */
+ signal (SIGPIPE, XDead);
+ fcntl (Xgraph[N], F_SETFL, O_NONBLOCK);
+ for (i = 0; (read (Xgraph[N], buffer, 64) > 0) && (i < 20); i++);
+ fcntl (Xgraph[N], F_SETFL, !O_NONBLOCK);
+
+ write (Xgraph[N], "ERAS", 4);
+
+ /* tell kapa to look for the incoming image */
+ status = write (Xgraph[N], "PLOT", 4);
+
+ /* send Xgraph[N] the plot details */
+ sprintf (buffer, "%8d %8d %d %d %d %d %f %f",
+ Npts, graphmode[0].style,
+ graphmode[0].ptype, graphmode[0].ltype,
+ graphmode[0].etype, graphmode[0].color,
+ graphmode[0].lweight, graphmode[0].size);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+ sprintf (buffer, "%f %f %f %f",
+ graphmode[0].xmin, graphmode[0].xmax,
+ graphmode[0].ymin, graphmode[0].ymax);
+ sprintf (buffer2, "NBYTES: %6d", strlen (buffer));
+ write (Xgraph[N], buffer2, 16);
+ write (Xgraph[N], buffer, strlen (buffer));
+
+}
+
+void PlotLabel (char *string, int N) {
+
+ char buffer[32];
+
+ if (Xgraph[N] == 0) return;
+
+ write (Xgraph[N], "LABL", 4);
+ sprintf (buffer, " %6d %6d", strlen (string), 2);
+ write (Xgraph[N], buffer, 16);
+ write (Xgraph[N], string, strlen (string));
+}
+
+void PlotVector (int Npts, double *vect, int mode, int N) {
+
+ float *values;
+ int i, Nbytes;
+
+ active = N;
+ if (Npts < 1) return;
+
+ ALLOCATE (values, float, Npts);
+ for (i = 0; i < Npts; i++) {
+ values[i] = vect[i];
+ }
+
+ Nbytes = Npts * sizeof (float);
+ write (Xgraph[N], values, Nbytes);
+ free (values);
+
+}
+
+/* include these lines to plot a pair of vectors:
+
+ typedef struct {
+ double xmin, xmax, ymin, ymax;
+ int style, ptype, ltype, etype, color;
+ double lweight, size;
+ } Graphdata;
+ Graphdata graphdata;
+
+ graphdata.xmin = -200;
+ graphdata.xmax = 4200;
+ graphdata.ymin = -500;
+ graphdata.ymax = 500;
+ graphdata.style = 2;
+ graphdata.ptype = 2;
+ graphdata.ltype = 0;
+ graphdata.etype = 0;
+ graphdata.color = 0;
+ graphdata.lweight = 0;
+ graphdata.size = 0.5;
+
+ PrepPlotting (N, &graphdata, n);
+ PlotVector (N, Y, 0, n);
+ PlotVector (N, dM, 1, n);
+ DonePlotting (&graphdata, n);
+
+ */
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/sort.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/sort.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/sort.c (revision 10360)
@@ -0,0 +1,128 @@
+# include "uniphot.h"
+
+void sort (unsigned int *X, int N) {
+
+ int l,j,ir,i;
+ unsigned tX;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ }
+}
+
+void sortB (double *X, double *Y, int N) {
+
+ int l,j,ir,i;
+ double tX, tY;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ tY = Y[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ tY = Y[ir];
+ Y[ir] = Y[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ Y[0] = tY;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ Y[i] = Y[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ Y[i] = tY;
+ }
+}
+
+void sortD (double *X, double *Y, double *Z, int N) {
+
+ int l,j,ir,i;
+ double tX, tY, tZ;
+
+ if (N < 2) return;
+
+ l = N >> 1;
+ ir = N - 1;
+ for (;;) {
+ if (l > 0) {
+ l--;
+ tX = X[l];
+ tY = Y[l];
+ tZ = Z[l];
+ }
+ else {
+ tX = X[ir];
+ X[ir] = X[0];
+ tY = Y[ir];
+ Y[ir] = Y[0];
+ tZ = Z[ir];
+ Z[ir] = Z[0];
+ if (--ir == 0) {
+ X[0] = tX;
+ Y[0] = tY;
+ Z[0] = tZ;
+ return;
+ }
+ }
+ i = l;
+ j = (l << 1) + 1;
+ while (j <= ir) {
+ if (j < ir && X[j] < X[j+1]) j++;
+ if (tX < X[j]) {
+ X[i] = X[j];
+ Y[i] = Y[j];
+ Z[i] = Z[j];
+ j += (i=j) + 1;
+ }
+ else j = ir + 1;
+ }
+ X[i] = tX;
+ Y[i] = tY;
+ Z[i] = tZ;
+ }
+}
+
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/subset_images.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/subset_images.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/subset_images.c (revision 10360)
@@ -0,0 +1,36 @@
+# include "uniphot.h"
+
+int subset_images (FITS_DB *db) {
+
+ int i, Nimage, equiv;
+ int Nkeep, *keep;
+ Image *image;
+
+ /* use a vtable to keep the images to be calibrated */
+ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+
+ Nkeep = 0;
+ ALLOCATE (keep, int, Nimage);
+
+ /* mark images to be calibrated */
+ for (i = 0; i < Nimage; i++) {
+
+ image[i].code |= ID_IMAGE_NOCAL;
+
+ /* select images by photcode */
+ equiv = GetPhotcodeEquivCodebyCode (image[i].source);
+ if (equiv != photcode[0].code) continue;
+
+ /* select images by time */
+ if (TimeSelect) {
+ if (image[i].tzero < TSTART) continue;
+ if (image[i].tzero > TSTOP) continue;
+ }
+ image[i].code &= ~ID_IMAGE_NOCAL;
+ keep[Nkeep] = i;
+ Nkeep ++;
+ }
+
+ gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, keep, Nkeep);
+ return (TRUE);
+}
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/uniphot.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/uniphot.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/uniphot.c (revision 10360)
@@ -0,0 +1,63 @@
+# include "uniphot.h"
+
+int main (int argc, char **argv) {
+
+ int i, Ntgroup, Nsgroup, status;
+ Group *tgroup, *sgroup;
+ ImageLink *imlinks;
+ FITS_DB db;
+
+ /* get configuration info, args, lockfile */
+ initialize (argc, argv);
+
+ set_db (&db);
+ status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
+ if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+ if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
+ if (!UPDATE) dvo_image_unlock (&db);
+
+ /* load images */
+ load_images (&db);
+
+ /* filter image list by selection */
+ subset_images (&db);
+
+ tgroup = find_image_tgroups (&db, &imlinks, &Ntgroup);
+ sgroup = find_image_sgroups (&db, &imlinks, &Nsgroup);
+
+ /* determine fit values */
+ for (i = 0; i < NLOOP; i++) {
+
+ fit_tgroup (tgroup, Ntgroup);
+ fit_sgroup (sgroup, Nsgroup);
+
+ }
+
+ fprintf (stdout, "# uniphot results for filter %s\n", photcode[0].name);
+ fprintf (stdout, "# STATMODE: %s\n", STATMODE);
+ fprintf (stdout, "# NLOOP: %d\n", NLOOP);
+ fprintf (stdout, "# time groups : %d\n", Ntgroup);
+ for (i = 0; i < Ntgroup; i++) {
+ fprintf (stdout, "%s %5d %5d %7.4f %7.4f %7.4f\n", tgroup[i].label,
+ tgroup[i].Nimage, tgroup[i].Ngood, 0.001*tgroup[i].M, 0.001*tgroup[i].dM, 0.001*tgroup[i].dMsub);
+ }
+ fprintf (stdout, "\n");
+
+ fprintf (stdout, "# space groups : %d\n", Nsgroup);
+ for (i = 0; i < Nsgroup; i++) {
+ fprintf (stdout, "%s %5d %5d %7.4f %7.4f %7.4f\n", sgroup[i].label,
+ sgroup[i].Nimage, sgroup[i].Ngood, 0.001*sgroup[i].M, 0.001*sgroup[i].dM, 0.001*sgroup[i].dMsub);
+ }
+ if (!UPDATE) exit (0);
+
+ update (&db, sgroup, Nsgroup);
+ dvo_image_unlock (&db);
+
+ exit (0);
+}
+
+/* add a mode to check / force consistency between image.Mcal and
+ measure.Mcal.
+
+ use the method in delstar to match image with measure (gregions (image), etc)
+*/
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update.c (revision 10360)
@@ -0,0 +1,122 @@
+# include "uniphot.h"
+# include
+
+void update (FITS_DB *db, Group *sgroup, int Nsgroup) {
+
+ int i, j, status, Nmin, Nimage;
+ char line[256];
+ glob_t pglob;
+ double Rmin, Rmax, Dmin, Dmax, x, y, radius;
+ Image *image;
+ Catalog catalog;
+ Coords coords;
+
+ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+
+ /* clear the NOCAL flags */
+ for (i = 0; i < Nimage; i++) {
+ image[i].code &= ~ID_IMAGE_NOCAL;
+ }
+
+ /* apply calculated space-group offset to image Mcal values */
+ for (i = 0; i < Nsgroup; i++) {
+ for (j = 0; j < sgroup[i].Nimage; j++) {
+ sgroup[i].image[j][0].Mcal_PS -= sgroup[i].M;
+ }
+ }
+
+ /** write image table **/
+ dvo_image_update (db, VERBOSE);
+
+ // XXX this process uses the existence of the file to perform the update
+ // XXX convert this to an examination of the SkyTable
+ /** update catalog tables **/
+ pglob.gl_offs = 0;
+ sprintf (line, "%s/*/*.cpt", CATDIR);
+ status = glob (line, 0, NULL, &pglob);
+
+ coords.crpix1 = coords.crpix2 = 0.0;
+ coords.cdelt1 = coords.cdelt2 = 1.0; /* scale is degrees, radius in degrees */
+ coords.pc1_1 = coords.pc2_2 = 1.0;
+ coords.pc1_2 = coords.pc2_1 = 0.0;
+ strcpy (coords.ctype, "RA---TAN");
+
+ for (i = 0; i < pglob.gl_pathc; i++) {
+ catalog.filename = pglob.gl_pathv[i];
+ catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+ catalog.Nsecfilt = GetPhotcodeNsecfilt ();
+
+ // an error exit status here is a significant error
+ if (!dvo_catalog_open (&catalog, NULL, VERBOSE, "a")) {
+ fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename);
+ exit (2);
+ }
+ if (!catalog.Nave_disk) {
+ dvo_catalog_unlock (&catalog);
+ dvo_catalog_free (&catalog);
+ continue;
+ }
+
+ gfits_scan (&catalog.header, "RA0", "%lf", 1, &Rmin);
+ gfits_scan (&catalog.header, "RA1", "%lf", 1, &Rmax);
+ gfits_scan (&catalog.header, "DEC0", "%lf", 1, &Dmin);
+ gfits_scan (&catalog.header, "DEC1", "%lf", 1, &Dmax);
+ while (Rmin < 0.0) { Rmin += 360.0; }
+ while (Rmin > 360.0) { Rmin -= 360.0; }
+ while (Rmax < 0.0) { Rmax += 360.0; }
+ while (Rmax > 360.0) { Rmax -= 360.0; }
+
+ coords.crval1 = 0.5*(Rmin + Rmax);
+ coords.crval2 = 0.5*(Dmin + Dmax);
+
+ Nmin = 0;
+ Rmin = 1000;
+ /* primitive version: match catalog with closest sgroup */
+ for (j = 0; j < Nsgroup; j++) {
+ if (!RD_to_XY (&x, &y, sgroup[j].v1, sgroup[j].v2, &coords)) continue;
+ radius = hypot (x, y);
+ if ((j == 0) || (radius < Rmin)) {
+ Rmin = radius;
+ Nmin = j;
+ }
+ }
+
+ fprintf (stderr, "catalog: %s sgroup: %d %s %f\n", catalog.filename, Nmin, sgroup[Nmin].label, Rmin);
+ update_dvo_catalog (&catalog, &sgroup[Nmin], (Rmin > 2*RADIUS));
+ dvo_catalog_save (&catalog, VERBOSE);
+ dvo_catalog_unlock (&catalog);
+ dvo_catalog_free (&catalog);
+ }
+}
+
+/* loop over all average
+
+ if (source equiv photcode)
+
+ loop over sgroups
+
+ if (radius < RADIUS)
+
+ average.M (+/-) source.M
+*/
+
+/* loop over all measure
+
+ if (source equiv photcode)
+
+ loop over sgroups
+
+ if (radius < RADIUS)
+
+ measure.Mcal (+/-) source.M
+*/
+
+/* alternative (slower, more robust)
+
+ loop over all average
+ loop over all average.Nm
+ select measure with source equiv photcode
+ match measure to image, get Mcal & offset
+ if any offset for average[i] is different, give error
+ apply offset to average.M
+*/
Index: /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update_catalog.c
===================================================================
--- /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update_catalog.c (revision 10360)
+++ /tags/uniphot-1-5-1/Ohana/src/uniphot/src/update_catalog.c (revision 10360)
@@ -0,0 +1,35 @@
+# include "uniphot.h"
+
+void update_dvo_catalog (Catalog *catalog, Group *sgroup, int warn) {
+
+ int i, j, m, found;
+ int Primary, Nsec, Nsecfilt;
+ float *Mp;
+ PhotCode *code;
+
+ Primary = (photcode[0].type == PHOT_PRI);
+ Nsec = GetPhotcodeNsec (photcode[0].code);
+ Nsecfilt = GetPhotcodeNsecfilt ();
+
+ found = 0;
+ for (i = 0; i < catalog[0].Naverage; i++) {
+
+ Mp = Primary ? &catalog[0].average[i].M : &catalog[0].secfilt[i*Nsecfilt+Nsec].M_PS;
+ if (*Mp != NO_MAG) *Mp += sgroup[0].M;
+
+ m = catalog[0].average[i].offset;
+ for (j = 0; j < catalog[0].average[i].Nm; j++, m++) {
+ code = GetPhotcodebyCode (catalog[0].measure[m].source);
+ if (code == NULL) continue;
+ if (code[0].type != PHOT_DEP) continue;
+ if (code[0].equiv != photcode[0].code) continue;
+ catalog[0].measure[m].Mcal_PS -= sgroup[0].M;
+ found ++;
+ }
+ }
+
+ if (found) {
+ fprintf (stderr, "found %d matches\n", found);
+ if (warn) fprintf (stderr, "warning: updated values, perhaps out of range?\n");
+ }
+}