Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/.mana =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/.mana (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/.mana (revision 37449) @@ -0,0 +1,113 @@ +gaussint +gaussint 10000 0.0 1.0 +gaussint vec 10000 0.0 1.0 +vstat vec +create n 0 vec[] +lim n vec; clear; box; plot n vec +create f 0 1000 +set F = rnd(f) +set r = F^2 +histogram r nr 0 1 0.001 +lim r nr; clear; box; plot r nr -x 1 +lim r nr; clear; box; plot r nr -x 0 +lim r nr; clear; box; plot r nr -x 2 +lim f F; clear; box; plot f F +lim f r; clear; box; plot f r +histogram r nr 0 1 0.001 -range dr +lim dr nr; clear; box; plot dr nr -x 1 +set r = sqrt(F) +histogram r nr 0 1 0.001 -range dr +lim dr nr; clear; box; plot dr nr -x 1 +create f 0 10000 +set F = rnd(f) +set r = sqrt(F) +histogram r nr 0 1 0.001 -range dr +histogram r nr 0 1 0.002 -range dr +lim dr nr; clear; box; plot dr nr -x 1 +create f 0 1000000 +set F = rnd(f) +set r = sqrt(F) +histogram r nr 0 1 0.002 -range dr +lim dr nr; clear; box; plot dr nr -x 1 +set theta = rnd(f)*360.0 +set x = r*dcos(theta); set set y = r*dsin(theta) +set x = r*dcos(theta); set y = r*dsin(theta) +lim x y; clear; box; plot x y +section default -imtool none; parity 0 1; densify im x y -graph; tv im -1 10 +clear -graph +set r = 2000*sqrt(F) +gaussdev +gaussdev z r[] 0.0 50.0 +vectors +set L = theta +set D = sqrt(r^2 + z^2) +set B = (z @ r) +set Balt = atan2 (z , r) +vectors +vstat B +vstat Balt +set dB = B - Balt +vstat dB +set RA = L +set DEC = B +csystem +csystem G C RA DEC +vectors +region 0 0 90 ait +box +tv +tv im +section default -imtool none ; parity 0 1 +cdensify +cdensify img RA DEC +tv img -1 10 +parity 0 1 +cur +cgrid -c red +cur +cur +tvcolor rainbow +tv img -1 100 +cur +stat img +set limg = log(img) +tv limg -0.1 2.0 +pw +input mana.sh +mkgal; showgal +showgal +showgal +head img +echo {3600*2.7691528201E-01} +stat img +keyword img CDELT2 dX +echo $MAX +echo {$MAX / $dX^2} +mkgal; showgal +$Z_GAL = 100 +mkgal; showgal +stat img +echo {$MAX / $dX^2} +input mana.sh +input mana.sh +mkgal; showgal +input mana.sh +mkgal; showgal +buffers +stat img +echo {$MAX / $dX^2} +cur +coords +coords img -c $R1 $D1 +stat +stat img 670 207 20 20 +echo $TOTAL +echo {$TOTAL / $NPIX / $dX^2} +echo {3000 / 5.5} +input mana.sh +mkgal; showgal +stat img 670 207 20 20 +echo {$TOTAL / $NPIX / $dX^2} +stat img +echo {$MAX / $dX^2} +echo {3000 / 34} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile (revision 37449) @@ -0,0 +1,59 @@ +default: fakeastro fakeastro_client +help: + @echo "make options: fakeastro (default)" + +include ../../Makefile.System +HOME = $(ROOT)/src/fakeastro +BIN = $(HOME)/bin +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +include ../../Makefile.Common + +# programs may add their own internal requirements here +FULL_CFLAGS = $(BASE_CFLAGS) +FULL_CPPFLAGS = $(BASE_CPPFLAGS) +FULL_LDFLAGS = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS) + +fakeastro: $(BIN)/fakeastro.$(ARCH) +fakeastro_client: $(BIN)/fakeastro_client.$(ARCH) + +install: $(DESTBIN)/fakeastro $(DESTBIN)/fakeastro_client $(DESTBIN)/testparallax + +FAKEASTRO = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/fakeastro.$(ARCH).o \ +$(SRC)/fakeastro_galaxy.$(ARCH).o \ +$(SRC)/make_fakestars.$(ARCH).o \ +$(SRC)/make_subset.$(ARCH).o \ +$(SRC)/save_fakestars.$(ARCH).o \ +$(SRC)/fakestar_catalog.$(ARCH).o \ +$(SRC)/fakestar_io.$(ARCH).o \ +$(SRC)/insert_fakestar.$(ARCH).o \ +$(SRC)/gaussian.$(ARCH).o \ +$(SRC)/remote_hosts.$(ARCH).o + +#$(SRC)/fakeastro_images.$(ARCH).o \ + +$(FAKEASTRO): $(INC)/fakeastro.h +$(BIN)/fakeastro.$(ARCH): $(FAKEASTRO) + +FAKEASTRO_CLIENT = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/fakeastro_client.$(ARCH).o \ +$(SRC)/fakestar_io.$(ARCH).o \ +$(SRC)/fakestar_catalog.$(ARCH).o \ +$(SRC)/insert_fakestar.$(ARCH).o \ + +$(FAKEASTRO_CLIENT): $(INC)/fakeastro.h +$(BIN)/fakeastro_client.$(ARCH): $(FAKEASTRO_CLIENT) + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/doc/example.sh =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/doc/example.sh (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/doc/example.sh (revision 37449) @@ -0,0 +1,29 @@ + +$Z_GAL = 500.0 +$R_GAL = 2000.0 +$N_GAL = 10000000 + +macro mkgal + create f 0 $N_GAL + set r = $R_GAL*sqrt(rnd(f)) + set L = rnd(f)*360.0 + gaussdev z r[] 0.0 $Z_GAL + set B = atan2 (z , r) + + set RA = L + set DEC = B + + csystem G C RA DEC +end + +macro showgal + region 0 0 90 ait + box + cgrid -c red + section default -imtool none ; parity 0 1 + cdensify img RA DEC + set limg = log(img) + tv limg -0.1 2.0 + tvcolor rainbow +end + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h (revision 37449) @@ -0,0 +1,94 @@ +# include +# include +# include +# include +# include +# include + +typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES} FakeastroOp; + +typedef struct { + double R, D; + StarPar starpar; + int flag; // in a subset? + int found; // assigned to an object? +} FakeAstro_Stars; + +/* used in find_matches, find_matches_refstars */ +# define IN_REGION(R,D) ( \ +((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \ +((R) >= region[0].Rmin) && ((R) < region[0].Rmax)) + +# define myAbortF(FORMAT,...) { fprintf (stderr, FORMAT, __VA_ARGS__); abort(); } + +/**** global variables set in parameter file ***/ + +# define DVO_MAX_PATH 1024 + +char ImageCat[DVO_MAX_PATH]; +char GSCFILE[DVO_MAX_PATH]; +char CATDIR[DVO_MAX_PATH]; +char SKY_TABLE[DVO_MAX_PATH]; + +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +int SKY_DEPTH; /* depth of catalog tables, fix usage */ + +int HOST_ID; +char *HOSTDIR; +char *CPT_FILE; +char *INPUT; + +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +int VERBOSE; +int VERBOSE2; + +SkyRegion UserPatch; + +FakeastroOp FAKEASTRO_OP; + +/*** fakeastro prototypes ***/ + +void ConfigInit PROTO((int *argc, char **argv)); +void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); +int args PROTO((int argc, char **argv)); +int args_client PROTO((int argc, char **argv)); + +void set_db PROTO((FITS_DB *in)); +int Shutdown PROTO((char *format, ...)) OHANA_FORMAT(printf, 1, 2); +void TrapSignal PROTO((int sig)); +void SetProtect PROTO((int mode)); +int SetSignals PROTO((void)); + +void initialize (int argc, char **argv); +void initialize_client (int argc, char **argv); + +int fakeastro_galaxy (SkyList *skylistInput, HostTable *hosts); +FakeAstro_Stars *make_fakestars (int Nstars); +int sortStars (FakeAstro_Stars *stars, int Nstars); + +FakeAstro_Stars *make_subset (FakeAstro_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset); + +int save_fakestars (FakeAstro_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname); + +FakeAstro_Stars *fakestar_load_stars (char *filename, int *nstars); +int fakestar_save_stars (char *filename, FakeAstro_Stars *stars, int Nstars); + +int harvest_host (); +int save_remote_host (HostInfo *host); +int harvest_all (); +int find_empty_slot (); +int init_remote_hosts (); + +int fakestar_catalog (FakeAstro_Stars *stars, int Nstars, SkyRegion *region, char *filename); +int insert_fakestar (SkyRegion *region, FakeAstro_Stars *stars, int Nstars, Catalog *catalog); + +double gaussian (double x, double mean, double sigma); +void gauss_init (int Nbin); +double rnd_gauss (double mean, double sigma); +double int_gauss (int i); + +int strextend (char *input, char *format,...); Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c (revision 37449) @@ -0,0 +1,71 @@ +# include "fakeastro.h" + +void ConfigInit (int *argc, char **argv) { + + char *config, *file; + char CatdirPhotcodeFile[256]; + struct stat filestat; + int status; + + /*** 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); + + // set defaults for all of these if they are not used by parallel / remote clients + // if (!ScanConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS)) ADDSTAR_RADIUS = 1.0; + + // force CATDIR to be absolute (so parallel mode will work) + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + char *tmpcatdir = abspath (CATDIR, DVO_MAX_PATH); + strcpy (CATDIR, tmpcatdir); + free (tmpcatdir); + + GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig(config, "CATMODE", "%s", 0, CATMODE); + ScanConfig(config, "CATFORMAT", "%s", 0, CATFORMAT); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) SKY_DEPTH = 2; + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) SKY_TABLE[0] = 0; + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + // check for existence of CATDIR + status = stat (CATDIR, &filestat); + if (status == -1) { + fprintf (stderr, "directory %s does not exist, giving up\n", CATDIR); + exit (1); + } + + /* update master photcode table if not defined */ + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } + SetZeroPoint (25.0); + + free (config); + free (file); + +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/SetSignals.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/SetSignals.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/SetSignals.c (revision 37449) @@ -0,0 +1,122 @@ +# include "fakeastro.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 SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ +# ifdef SIGPWR + case SIGPWR: /* power failure - why ignore this? (Sys V) */ +# endif +# ifdef SIGWINCH + case SIGWINCH: /* window resized (4.3BSD) */ +# endif + 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: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/Shutdown.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/Shutdown.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/Shutdown.c (revision 37449) @@ -0,0 +1,30 @@ +# include "fakeastro.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: relastro halted\n"); + exit (1); +} + + +/* XXX this is probably not needed anymore : just protect the write statements */ Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/args.c (revision 37449) @@ -0,0 +1,236 @@ +# include "fakeastro.h" +void usage (void); +void usage_client (void); + +int args (int argc, char **argv) { + + int N; + + /* possible operations */ + FAKEASTRO_OP = OP_NONE; + + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + FAKEASTRO_OP = OP_IMAGES; + } + + if ((N = get_argument (argc, argv, "-galaxy"))) { + if (FAKEASTRO_OP != OP_NONE) usage(); + remove_argument (N, &argc, argv); + FAKEASTRO_OP = OP_GALAXY; + } + + if (FAKEASTRO_OP == OP_NONE) usage(); + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + UserPatch.Rmax = UserPatch.Rmin + 0.001; + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + UserPatch.Dmax = UserPatch.Dmin + 0.001; + remove_argument (N, &argc, argv); + } + + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + PARALLEL = FALSE; + if ((N = get_argument (argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // fakeastro will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, &argc, argv); + } + // this is a test mode : rather than launching the fakeastro_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; + remove_argument (N, &argc, argv); + } + + // MaxDensityUse = FALSE; + // if ((N = get_argument (argc, argv, "-max-density"))) { + // remove_argument (N, &argc, argv); + // MaxDensityValue = atof(argv[N]); + // remove_argument (N, &argc, argv); + // MaxDensityUse = TRUE; + // } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE = VERBOSE2 = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 1) usage (); + return TRUE; +} + +int args_client (int argc, char **argv) { + + int N; + + /* possible operations */ + FAKEASTRO_OP = OP_NONE; + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) usage_client(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) usage_client(); + + CPT_FILE = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + CPT_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!CPT_FILE) usage_client(); + + INPUT = NULL; + if ((N = get_argument (argc, argv, "-input"))) { + remove_argument (N, &argc, argv); + INPUT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!INPUT) usage_client(); + + if ((N = get_argument (argc, argv, "-images"))) { + remove_argument (N, &argc, argv); + FAKEASTRO_OP = OP_IMAGES; + } + if ((N = get_argument (argc, argv, "-galaxy"))) { + if (FAKEASTRO_OP != OP_NONE) usage(); + remove_argument (N, &argc, argv); + FAKEASTRO_OP = OP_GALAXY; + } + if (FAKEASTRO_OP == OP_NONE) usage_client(); + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + UserPatch.Rmax = UserPatch.Rmin + 0.001; + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + UserPatch.Dmax = UserPatch.Dmin + 0.001; + remove_argument (N, &argc, argv); + } + + // MaxDensityUse = FALSE; + // if ((N = get_argument (argc, argv, "-max-density"))) { + // remove_argument (N, &argc, argv); + // MaxDensityValue = atof(argv[N]); + // remove_argument (N, &argc, argv); + // MaxDensityUse = TRUE; + // } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE = VERBOSE2 = TRUE; + remove_argument (N, &argc, argv); + } + + if (argc != 1) usage_client (); + return TRUE; +} + +void usage () { + fprintf (stderr, "ERROR: USAGE: fakeastro -images [options]\n"); + fprintf (stderr, " OR: fakeastro -galaxy [options]\n"); + + fprintf (stderr, " specify one of the following modes: \n"); + fprintf (stderr, " -images\n"); + fprintf (stderr, " -galaxy\n"); + + fprintf (stderr, " additional options: \n"); + fprintf (stderr, " -region RA RA DEC DEC\n"); + fprintf (stderr, " -catalog (ra) (dec)\n\n"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " -vv\n"); + fprintf (stderr, " \n"); + exit (2); +} + +void usage_client () { + fprintf (stderr, "ERROR: USAGE: fakeastro_client -images\n"); + fprintf (stderr, " OR: fakeastro_client -galaxy\n"); + + fprintf (stderr, " specify one of the following modes: \n"); + fprintf (stderr, " -images\n"); + fprintf (stderr, " -galaxy\n"); + + fprintf (stderr, " additional options: \n"); + fprintf (stderr, " -region Rmin Rmax Dmin Dmax"); + fprintf (stderr, " -catalog RA DEC"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " -vv\n"); + fprintf (stderr, " \n"); + exit (2); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c (revision 37449) @@ -0,0 +1,57 @@ +# include "fakeastro.h" + +int main (int argc, char **argv) { + + /* get configuration info, args */ + initialize (argc, argv); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + // load the list of hosts + HostTable *hosts = NULL; + if (PARALLEL) { + hosts = HostTableLoad (CATDIR, sky->hosts); + if (!hosts) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + // ensure that the paths are absolute path names + int i; + for (i = 0; i < hosts->Nhosts; i++) { + char *tmppath = abspath (hosts->hosts[i].pathname, DVO_MAX_PATH); + free (hosts->hosts[i].pathname); + hosts->hosts[i].pathname = tmppath; + } + + // set up the array of active hosts + init_remote_hosts (); + } + + switch (FAKEASTRO_OP) { + case OP_GALAXY: + /* the object analysis is a separate process iterating over catalogs */ + fakeastro_galaxy (skylist, hosts); + exit (0); + + case OP_IMAGES: + // fakeastro_images (skylist); + exit (0); + + default: + fprintf (stderr, "impossible!\n"); + abort(); + } + exit (1); +} + +/* fakeastro can be used to generate a DVO database of fake observations. + + top-level modes: + + -galaxy : generate a galaxy model. output is a DVO database with average.d and starpar.d tables populated + -images : generate a set of observations (or just images?) + +*/ Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_client.c (revision 37449) @@ -0,0 +1,24 @@ +# include "fakeastro.h" + +int main (int argc, char **argv) { + + // need to construct these options with args_loadstarpar... + ConfigInit (&argc, argv); + args_client (argc, argv); + + // client is called with a pointer to the file to be loaded + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + SkyList *skylist = SkyRegionByCPT (sky, CPT_FILE); + + int Nstars; + FakeAstro_Stars *stars = fakestar_load_stars (INPUT, &Nstars); + + fakestar_catalog (stars, Nstars, skylist->regions[0], CPT_FILE); + + free (stars); + + exit (0); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_galaxy.c (revision 37449) @@ -0,0 +1,57 @@ +# include "fakeastro.h" + +int fakeastro_galaxy (SkyList *skylistInput, HostTable *hosts) { + + int n, i; + + // INITTIME; + + // XXX check for an existing CATDIR and exit? + + /* outline + + * generate N stars (random draw from galaxy mode) + * assign them to catalogs and distribute + + */ + + gauss_init (50000); + + int Nloop = 10; + int Nstars = 1000000; + + for (n = 0; n < Nloop; n++) { + + FakeAstro_Stars *stars = make_fakestars (Nstars); + + // sort the stars by RA + sortStars (stars, Nstars); + + // scan through the stars, loading the containing catalogs + // skip through table for unsaved stars + for (i = 0; i < Nstars; i++) { + if (stars[i].flag) continue; + + // identify the relevant catalog + SkyList *skylist = SkyRegionByPoint_List (skylistInput, -1, stars[i].R, stars[i].D); + if (skylist[0].Nregions == 0) { + SkyListFree (skylist); + continue; + } + SkyRegion *region = skylist[0].regions[0]; + + // select stars matching this region + int Nsubset; + FakeAstro_Stars *subset = make_subset (stars, Nstars, i, region, &Nsubset); + + // In parallel mode, write out the subset to a disk file. Block until a remote host + // is available. In serial mode, just match against the appropriate region and save + save_fakestars (subset, Nsubset, hosts, region, skylist[0].filename[0]); + } + + // wait for last remote clients to finish + harvest_all (); + } + + return TRUE; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_catalog.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_catalog.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_catalog.c (revision 37449) @@ -0,0 +1,27 @@ +# include "fakeastro.h" + +int fakestar_catalog (FakeAstro_Stars *stars, int Nstars, SkyRegion *region, char *filename) { + + Catalog catalog; + + // now we have all of the loaded stars in this catalog + catalog.filename = filename; + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_STARPAR; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, region, VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + insert_fakestar (region, stars, Nstars, &catalog); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + return (TRUE); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c (revision 37449) @@ -0,0 +1,188 @@ +# include "fakeastro.h" + +int fakestar_save_stars (char *filename, FakeAstro_Stars *stars, int Nstars) { + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + myAbortF ("ERROR: cannot open file for output %s\n", filename); + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + gfits_create_table_header (&theader, "BINTABLE", "STARPAR"); + + gfits_define_bintable_column (&theader, "D", "RA", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "GLON", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "GLAT", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Ebv", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dEbv", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DistMag", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dDistMag", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "M_r", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dM_r", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "FeH", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dFeH", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "averef", "", "", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "objID", "", "", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "catID", "", "", 1.0, FT_BZERO_INT32); + + int i; + double *ra = NULL; ALLOCATE (ra , double , Nstars); for (i = 0; i < Nstars; i++) ra [i] = stars[i].R ; + double *dec = NULL; ALLOCATE (dec , double , Nstars); for (i = 0; i < Nstars; i++) dec [i] = stars[i].D ; + double *glon = NULL; ALLOCATE (glon , double , Nstars); for (i = 0; i < Nstars; i++) glon [i] = stars[i].starpar.galLon ; + double *glat = NULL; ALLOCATE (glat , double , Nstars); for (i = 0; i < Nstars; i++) glat [i] = stars[i].starpar.galLat ; + float *Ebv = NULL; ALLOCATE (Ebv , float , Nstars); for (i = 0; i < Nstars; i++) Ebv [i] = stars[i].starpar.Ebv ; + float *dEbv = NULL; ALLOCATE (dEbv , float , Nstars); for (i = 0; i < Nstars; i++) dEbv [i] = stars[i].starpar.dEbv ; + float *DistMag = NULL; ALLOCATE (DistMag , float , Nstars); for (i = 0; i < Nstars; i++) DistMag [i] = stars[i].starpar.DistMag ; + float *dDistMag = NULL; ALLOCATE (dDistMag, float , Nstars); for (i = 0; i < Nstars; i++) dDistMag[i] = stars[i].starpar.dDistMag; + float *M_r = NULL; ALLOCATE (M_r , float , Nstars); for (i = 0; i < Nstars; i++) M_r [i] = stars[i].starpar.M_r ; + float *dM_r = NULL; ALLOCATE (dM_r , float , Nstars); for (i = 0; i < Nstars; i++) dM_r [i] = stars[i].starpar.dM_r ; + float *FeH = NULL; ALLOCATE (FeH , float , Nstars); for (i = 0; i < Nstars; i++) FeH [i] = stars[i].starpar.FeH ; + float *dFeH = NULL; ALLOCATE (dFeH , float , Nstars); for (i = 0; i < Nstars; i++) dFeH [i] = stars[i].starpar.dFeH ; + unsigned int *averef = NULL; ALLOCATE (averef , unsigned int, Nstars); for (i = 0; i < Nstars; i++) averef [i] = stars[i].starpar.averef ; + unsigned int *objID = NULL; ALLOCATE (objID , unsigned int, Nstars); for (i = 0; i < Nstars; i++) objID [i] = stars[i].starpar.objID ; + unsigned int *catID = NULL; ALLOCATE (catID , unsigned int, Nstars); for (i = 0; i < Nstars; i++) catID [i] = stars[i].starpar.catID ; + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", ra , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DEC", dec , Nstars); + gfits_set_bintable_column (&theader, &ftable, "GLON", glon , Nstars); + gfits_set_bintable_column (&theader, &ftable, "GLAT", glat , Nstars); + gfits_set_bintable_column (&theader, &ftable, "Ebv", Ebv , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dEbv", dEbv , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DistMag", DistMag , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dDistMag", dDistMag, Nstars); + gfits_set_bintable_column (&theader, &ftable, "M_r", M_r , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dM_r", dM_r , Nstars); + gfits_set_bintable_column (&theader, &ftable, "FeH", FeH , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dFeH", dFeH , Nstars); + gfits_set_bintable_column (&theader, &ftable, "averef", averef , Nstars); + gfits_set_bintable_column (&theader, &ftable, "objID", objID , Nstars); + gfits_set_bintable_column (&theader, &ftable, "catID", catID , Nstars); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + return TRUE; +} + + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +FakeAstro_Stars *fakestar_load_stars (char *filename, int *nstars) { + + int i, Ncol; + off_t Nrow; + char type[16]; + + FakeAstro_Stars *stars = NULL; + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read matrix\n"); + goto escape; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + GET_COLUMN(ra , "RA" , double); + GET_COLUMN(dec , "DEC" , double); + GET_COLUMN(glon , "GLON" , double); + GET_COLUMN(glat , "GLAT" , double); + GET_COLUMN(Ebv , "Ebv" , float); + GET_COLUMN(dEbv , "dEbv" , float); + GET_COLUMN(DistMag , "DistMag" , float); + GET_COLUMN(dDistMag, "dDistMag" , float); + GET_COLUMN(M_r , "M_r" , float); + GET_COLUMN(dM_r , "dM_r" , float); + GET_COLUMN(FeH , "FeH" , float); + GET_COLUMN(dFeH , "dFeH" , float); + GET_COLUMN(averef , "averef" , unsigned int); + GET_COLUMN(objID , "objID" , unsigned int); + GET_COLUMN(catID , "catID" , unsigned int); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + ALLOCATE (stars, FakeAstro_Stars, Nrow); + + for (i = 0; i < Nrow; i++) { + stars[i].R = ra [i]; + stars[i].D = dec [i]; + stars[i].starpar.galLon = glon [i]; + stars[i].starpar.galLat = glat [i]; + stars[i].starpar.Ebv = Ebv [i]; + stars[i].starpar.dEbv = dEbv [i]; + stars[i].starpar.DistMag = DistMag [i]; + stars[i].starpar.dDistMag = dDistMag[i]; + stars[i].starpar.M_r = M_r [i]; + stars[i].starpar.dM_r = dM_r [i]; + stars[i].starpar.FeH = FeH [i]; + stars[i].starpar.dFeH = dFeH [i]; + stars[i].starpar.averef = averef [i]; + stars[i].starpar.objID = objID [i]; + stars[i].starpar.catID = catID [i]; + } + + fclose (f); + *nstars = Nrow; + return stars; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + return NULL; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/gaussian.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/gaussian.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/gaussian.c (revision 37449) @@ -0,0 +1,75 @@ +# include "fakeastro.h" + +static int Ngaussint = 0; +static double *gaussint; + +extern double drand48(); + +double gaussian (double x, double mean, double sigma) { + + double f; + + f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); + + return (f); + +} + +/* integrate a gaussian from -5 sigma to +5 sigma */ +void gauss_init (int Nbin) { + + int i; + long A, B; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + /* no need to generate this if it already exists */ + if (Ngaussint == Nbin) return; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + + Ngaussint = Nbin; + ALLOCATE (gaussint, double, Ngaussint + 1); + + val = 0; + dx = 1.0 / Ngaussint; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx) { + df = (3.0*gaussian(x , mean, sigma) + + 9.0*gaussian(x+dx1, mean, sigma) + + 9.0*gaussian(x+dx2, mean, sigma) + + 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i + 0.5) / (double) Ngaussint) { + gaussint[i] = x + dx / 2.0; + i++; + } + } +} + +double rnd_gauss (double mean, double sigma) { + + int i; + double y; + + y = drand48(); + i = Ngaussint*y; + y = gaussint[i]*sigma + mean; + + return (y); + +} + +double int_gauss (int i) { + double y; + y = gaussint[i]; + return (y); +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/initialize.c (revision 37449) @@ -0,0 +1,14 @@ +# include "fakeastro.h" + +void initialize (int argc, char **argv) { + + ConfigInit (&argc, argv); + args (argc, argv); +} + +void initialize_client (int argc, char **argv) { + + ConfigInit (&argc, argv); + args_client (argc, argv); +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/insert_fakestar.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/insert_fakestar.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/insert_fakestar.c (revision 37449) @@ -0,0 +1,75 @@ +# include "fakeastro.h" + +// assume we are working with an empty database? allow for matches? +int insert_fakestar (SkyRegion *region, FakeAstro_Stars *stars, int Nstars, Catalog *catalog) { + + off_t i, j; + off_t Nave, NAVE, Nstarpar, NSTARPAR, Nmatch; + + /* internal counters */ + Nmatch = 0; + NSTARPAR = Nstarpar = catalog[0].Nstarpar; + + // current max obj ID for this catalog + unsigned int objID = catalog[0].objID; + unsigned int catID = catalog[0].catID; + + int Nsecfilt = catalog[0].Nsecfilt; + + /** incorporate unmatched image stars? **/ + for (i = 0; i < Nstars; i++) { + + // skip already matched stars + if (stars[i].found) continue; + if (!IN_REGION (stars[i].R, stars[i].D)) continue; + + /* make sure there is space for next entry */ + if (Nstarpar >= NSTARPAR) { + NSTARPAR = Nstarpar + 1000; + REALLOCATE (catalog[0].starpar, StarPar, NSTARPAR); + } + if (Nave >= NAVE) { + NAVE = Nave + 1000; + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*Nsecfilt); + } + + dvo_average_init (&catalog[0].average[Nave]); + catalog[0].average[Nave].R = stars[i].R; + catalog[0].average[Nave].D = stars[i].D; + + catalog[0].average[Nave].Nstarpar = 1; + catalog[0].average[Nave].starparOffset = Nstarpar; + catalog[0].average[Nave].objID = objID; + catalog[0].average[Nave].catID = catID; + + objID ++; + + for (j = 0; j < Nsecfilt; j++) { + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); + } + + catalog[0].starpar[Nstarpar] = stars[i].starpar; + catalog[0].starpar[Nstarpar].averef = Nave; + catalog[0].starpar[Nstarpar].objID = catalog[0].average[Nave].objID; + catalog[0].starpar[Nstarpar].catID = catID; + + stars[i].found = TRUE; + Nstarpar ++; + Nave ++; + } + + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].starpar, StarPar, Nstarpar); + + catalog[0].sorted = FALSE; + + /* check if the catalog has changed? if no change, no need to write */ + catalog[0].objID = objID; // new max value, save on catalog close + catalog[0].Naverage = Nave; + catalog[0].Nstarpar = Nstarpar; + catalog[0].Nsecf_mem = Nave*Nsecfilt; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nstarpar: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n", Nstars, Nave, Nstarpar, Nmatch); + + return (Nmatch); +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c (revision 37449) @@ -0,0 +1,80 @@ +# include "fakeastro.h" + +# define Z_GAL 50.0 +# define R_GAL 2000.0 + +FakeAstro_Stars *make_fakestars (int Nstars) { + + FakeAstro_Stars *stars; + ALLOCATE (stars, FakeAstro_Stars, Nstars); + + // libdvo uses Liu et al 2011 (A&A 526, A16) for default galactic coords, + // but Greg Green / LSD use the older Reid et al 2004 (ApJ 616, 872) definition + CoordTransform *transform = InitTransform (COORD_GALACTIC, COORD_CELESTIAL); + + int i; + for (i = 0; i < Nstars; i++) { + + // each star has the following value drawn from an appropriate random distribution: + + // ra, dec, distance, (absolute mag, colors) or (mass) or (MK type) + + /* very simple galaxy model: + + * rho(x,y,z) = exp(-0.5*(z/Zpc)^2) (r < 2000) + + */ + + double z = rnd_gauss (0.0, Z_GAL); + double r = sqrt(drand48()) * R_GAL; + double L = drand48() * 2 * M_PI; + + double distance = sqrt (SQ(r) + SQ(z)); + double B = atan2(z,r); + + // double x = r*cos(L); + // double y = r*sin(L); + + // crude Mr distribution from + // http://iopscience.iop.org/1538-3881/139/6/2679/pdf/aj_139_6_2679.pdf + + double Mr = rnd_gauss (11.25, 1.0); + + double R, D; + ApplyTransform (&R, &D, L, B, transform); + + stars[Nstars].R = R; + stars[Nstars].D = D; + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + stars[Nstars].starpar.galLon = L; + stars[Nstars].starpar.galLat = B; + + stars[Nstars].starpar.Ebv = 0.0; + stars[Nstars].starpar.dEbv = 0.0; + stars[Nstars].starpar.DistMag = 5.0*log10(distance) - 5.0; + stars[Nstars].starpar.dDistMag = 0.0; + stars[Nstars].starpar.M_r = Mr; + stars[Nstars].starpar.dM_r = 0.0; + stars[Nstars].starpar.FeH = 0.0; + stars[Nstars].starpar.dFeH = 0.0; + + Nstars ++; + } + return stars; +} + +int sortStars (FakeAstro_Stars *stars, int Nstars) { + +# define SWAPFUNC(A,B){ FakeAstro_Stars temp = stars[A]; stars[A] = stars[B]; stars[B] = temp; } +# define COMPARE(A,B)(stars[A].R < stars[B].R) + + OHANA_SORT (Nstars, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + + return TRUE; +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_subset.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_subset.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_subset.c (revision 37449) @@ -0,0 +1,43 @@ +# include "fakeastro.h" + +// assign stars in the given region to the subset (NOTE: stars are sorted by RA, start is +// first entry in stars array in this region) + +FakeAstro_Stars *make_subset (FakeAstro_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) { + + int i; + + FakeAstro_Stars *subset = NULL; + + // collect array of (Stars *) stars in a new output catalog + int Nsubset = 0; + int NSUBSET = 3000; + ALLOCATE (subset, FakeAstro_Stars, NSUBSET); + + // find the rest of the stars in this output region + for (i = start; i < Nstars; i++) { + if (stars[i].flag) continue; + + // check if in skyregion + if (stars[i].R < region[0].Rmin) continue; + if (stars[i].R > region[0].Rmax) break; + if (stars[i].D < region[0].Dmin) continue; + if (stars[i].D > region[0].Dmax) continue; + + // check if in UserPatch (a GLOBAL) + if (stars[i].R < UserPatch.Rmin) continue; + if (stars[i].R > UserPatch.Rmax) break; + if (stars[i].D < UserPatch.Dmin) continue; + if (stars[i].D > UserPatch.Dmax) continue; + + subset[Nsubset] = stars[i]; + Nsubset ++; + + stars[i].flag = TRUE; + + CHECK_REALLOCATE (subset, FakeAstro_Stars, NSUBSET, Nsubset, 10000); + } + + *nsubset = Nsubset; + return subset; +} Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/remote_hosts.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/remote_hosts.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/remote_hosts.c (revision 37449) @@ -0,0 +1,140 @@ +# include "fakeastro.h" +# include +# include + +# define DEBUG 1 + +// we are running N parallel remote jobs on the remote hosts. In the usual remote +// processing (eg, relphot or relastro), we launch one job per remote host. In this case, +// we are launching one job foreach table (catalog.cpt) being touched. individual hosts +// may have more than one job active at a time. + +int Nhosts = 0; +HostInfo **hosts = NULL; + +int init_remote_hosts () { + + Nhosts = 10; + ALLOCATE (hosts, HostInfo *, Nhosts); + + int i; + for (i = 0; i < Nhosts; i++) { + hosts[i] = NULL; + } + return TRUE; +} + +int find_empty_slot () { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) return i; + } + return -1; +} + +int harvest_all () { + + int slot = -1; + + while (slot != -2) { + slot = harvest_host (); + usleep (50000); + } + return TRUE; +} + +// put the host in the list in a free slot. NOTE: this should never be called if we do not +// already have a free slot. +int save_remote_host (HostInfo *host) { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) { + hosts[i] = host; + return TRUE; + } + } + myAbort ("failed to find an empty slot: this is a programming error"); + return FALSE; +} + +// wait for all children to complete, report output to stdout +// possible states when function is finished: +// * no child was harvested, but children still exist +// * no outstanding children +// * a child was harvested +int harvest_host () { + + // check if any children have finished... + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + + if (!pid) return -1; // there are outstanding children, but none have exited, no slot was opened + if ((pid == -1) && (errno == ECHILD)) return -2; // there are no outstanding children + if ((pid == -1) && (errno != ECHILD)) myAbort ("programming error?"); + + // find the host which has finished + int i; + int found = FALSE; + for (i = 0; (i < Nhosts) && !found; i++) { + if (!hosts[i]) continue; // unassigned slot + if (hosts[i][0].pid != pid) continue; + found = TRUE; + } + myAssert (found, "Programming error: failed to matched finished job to known host!"); + int slot = i; + + HostInfo *host = hosts[slot]; + + // check on the status of this job and report any output + if (DEBUG) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + int Nout, printHead; + + // read stdout + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (DEBUG) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + } + + // close opened connections + close (host->stdio[HOST_STDIN]); + close (host->stdio[HOST_STDOUT]); + close (host->stdio[HOST_STDERR]); + + // free data associated with the host + free (host->hostname); + free (host->pathname); + FreeIOBuffer (&host->stdout); + FreeIOBuffer (&host->stderr); + free (host); + + // free the slot + hosts[slot] = NULL; + return slot; +} + Index: /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fakestars.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fakestars.c (revision 37449) +++ /branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fakestars.c (revision 37449) @@ -0,0 +1,84 @@ +# include "fakeastro.h" +# define DEBUG 1 + +int save_fakestars (FakeAstro_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // if this region is a parallel thing, save and launch remote + if (!PARALLEL) { + fakestar_catalog (stars, Nstars, region, fullname); + } else { + int N = hosts->index[region->hostID]; + HostInfo *hostMach = &hosts->hosts[N]; + + // save to a unique filename + char filename[1024]; // CATDIR/tmpdir/starpar.PID.index.fits + snprintf (filename, 1024, "%s/tmpdir/starpar.%s.%05d.fits", CATDIR, uniquer, region->index); + + // write the data to the given FITS file + fakestar_save_stars (filename, stars, Nstars); + + int slot = -1; + while (slot == -1) { + slot = find_empty_slot (); + if (slot == -1) { + usleep (50000); + slot = harvest_host(); + myAssert (slot != -2, "we should not call harvest_host here if we have open slots"); + } + } + + // allocate a host for this job + HostInfo *host = NULL; + ALLOCATE (host, HostInfo, 1); + + // we want to run this job on the host described by hostMach. copy + // immutable data from hostMach to a locally allocated host: + host->hostID = hostMach->hostID; + host->hostname = strcreate (hostMach->hostname); + host->pathname = strcreate (hostMach->pathname); + InitIOBuffer (&host->stdout, 1000); + InitIOBuffer (&host->stderr, 1000); + + // got a valid slot, so launch a new host + + // need to generate the remote command + char command[1024]; + snprintf (command, 1024, "fakestar_client"); + strextend (command, "-galaxy"); + strextend (command, "-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-hostdir %s", host->pathname); + strextend (command, "-cpt %s", region->name); + strextend (command, "-input %s", filename); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + + save_remote_host (host); + } + return TRUE; +} + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +}