Index: /branches/eam_branch_20081124/Ohana/src/addstar/Makefile =================================================================== --- /branches/eam_branch_20081124/Ohana/src/addstar/Makefile (revision 20904) +++ /branches/eam_branch_20081124/Ohana/src/addstar/Makefile (revision 20905) @@ -25,4 +25,5 @@ load2mass : $(BIN)/load2mass.$(ARCH) gztest : $(BIN)/gztest.$(ARCH) +mkcmf : $(BIN)/mkcmf.$(ARCH) all: addstar addstarc addstard addstart sedstar load2mass skycells @@ -222,4 +223,7 @@ $(SRC)/SetSignals.$(ARCH).o +MKCMF = \ +$(SRC)/mkcmf.$(ARCH).o + SKYCELLS = \ $(SRC)/skycells.$(ARCH).o \ @@ -230,21 +234,23 @@ $(SRC)/sky_tessalation.$(ARCH).o -$(ADDSTARC) : $(INC)/addstar.h -$(ADDSTARD) : $(INC)/addstar.h -$(ADDSTART) : $(INC)/addstar.h -$(ADDSTAR) : $(INC)/addstar.h -$(SKYCELLS) : $(INC)/addstar.h +$(ADDSTARC) : $(INC)/addstar.h +$(ADDSTARD) : $(INC)/addstar.h +$(ADDSTART) : $(INC)/addstar.h +$(ADDSTAR) : $(INC)/addstar.h +$(SKYCELLS) : $(INC)/addstar.h $(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h - -$(BIN)/addstar.$(ARCH) : $(ADDSTAR) -$(BIN)/addstard.$(ARCH) : $(ADDSTARD) -$(BIN)/addstart.$(ARCH) : $(ADDSTART) -$(BIN)/addstarc.$(ARCH) : $(ADDSTARC) +$(MKCMF) : $(INC)/mkcmf.h + +$(BIN)/addstar.$(ARCH) : $(ADDSTAR) +$(BIN)/addstard.$(ARCH) : $(ADDSTARD) +$(BIN)/addstart.$(ARCH) : $(ADDSTART) +$(BIN)/addstarc.$(ARCH) : $(ADDSTARC) $(BIN)/mkacc-2mass.$(ARCH) : $(MKACC-2MASS) -$(BIN)/sedstar.$(ARCH) : $(SEDSTAR) -$(BIN)/load2mass.$(ARCH) : $(LOAD-2MASS) -$(BIN)/skycells.$(ARCH) : $(SKYCELLS) - -INSTALL = addstar addstarc addstard addstart sedstar load2mass skycells +$(BIN)/sedstar.$(ARCH) : $(SEDSTAR) +$(BIN)/load2mass.$(ARCH) : $(LOAD-2MASS) +$(BIN)/skycells.$(ARCH) : $(SKYCELLS) +$(BIN)/mkcmf.$(ARCH) : $(MKCMF) + +INSTALL = addstar addstarc addstard addstart sedstar load2mass skycells mkcmf # dependancy rules for binary code ######################### Index: /branches/eam_branch_20081124/Ohana/src/addstar/include/mkcmf.h =================================================================== --- /branches/eam_branch_20081124/Ohana/src/addstar/include/mkcmf.h (revision 20905) +++ /branches/eam_branch_20081124/Ohana/src/addstar/include/mkcmf.h (revision 20905) @@ -0,0 +1,17 @@ +# include +# include +# include +# include +# include +# include + +/* solaris requires both of these instead of ip.h: + # include + # include +*/ + +/* linux is happy with this, not solaris */ +# include +# include +# include +# include Index: /branches/eam_branch_20081124/Ohana/src/addstar/src/mkcmf.c =================================================================== --- /branches/eam_branch_20081124/Ohana/src/addstar/src/mkcmf.c (revision 20905) +++ /branches/eam_branch_20081124/Ohana/src/addstar/src/mkcmf.c (revision 20905) @@ -0,0 +1,257 @@ +# include "mkcmf.h" +void gauss_init (int Nbin); +double rnd_gauss (double mean, double sigma); + +# define SKY 100.0 +# define DSKY 2.0 +# define PSFCHI 1.3 +# define CRN 3.0 +# define EXTN 4.0 +# define FX 1.5 +# define FY 1.0 +# define DF 10.0 +# define PSFQUAL 0.98 +# define FLAGS 0x1101 + +int main (int argc, char **argv) { + + // generate a simple, fake cmf file + // load a text table with X,Y,Mag (instrumental?) + + int i, N, Nstars, NSTARS; + double *X, *Y, *M; + double Xmax, Ymax; + float flux, fSN; + + FILE *f, *fits; + Header header, theader; + Matrix matrix; + FTable ftable; + Coords coords; + + PS1_DEV_1 *stars; + + static char *date = "2001-01-01"; + if ((N = get_argument (argc, argv, "-date"))) { + remove_argument (N, &argc, argv); + date = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + static char *time = "00:00:00"; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + time = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + double RA = 10.0; + double DEC = 20.0; + if ((N = get_argument (argc, argv, "-radec"))) { + remove_argument (N, &argc, argv); + RA = atof (argv[N]); + remove_argument (N, &argc, argv); + DEC = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 3) { + fprintf (stderr, "USAGE mkcmf (input) (output) [-date date] [-time time] [-radec ra dec]\n"); + exit (2); + } + + f = fopen (argv[1], "r"); + if (f == NULL) { + fprintf (stderr, "unable to open input file %s\n", argv[1]); + exit (1); + } + + Nstars = 0; + NSTARS = 100; + ALLOCATE (X, double, NSTARS); + ALLOCATE (Y, double, NSTARS); + ALLOCATE (M, double, NSTARS); + + Xmax = Ymax = 0; + + while (fscanf (f, "%lf %lf %lf", &X[Nstars], &Y[Nstars], &M[Nstars]) != EOF) { + if (Nstars == NSTARS - 1) { + NSTARS += 100; + REALLOCATE (X, double, NSTARS); + REALLOCATE (Y, double, NSTARS); + REALLOCATE (M, double, NSTARS); + } + Nstars ++; + } + + // XXX add gaussian-distributed noise based on counts + ALLOCATE (stars, PS1_DEV_1, Nstars); + gauss_init (2048); + for (i = 0; i < Nstars; i++) { + stars[i].detID = i; + + flux = pow (10.0, -0.4*M[i]); + fSN = 1.0 / sqrt(flux); + + stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0); + + Xmax = MAX(Xmax, X[i]); + Ymax = MAX(Ymax, Y[i]); + + stars[i].dX = FX * fSN; + stars[i].dY = FY * fSN; + stars[i].dM = fSN; + + stars[i].Mpeak = M[i] + 1.0; + stars[i].sky = SKY; + stars[i].dSky = DSKY; + stars[i].psfChisq = PSFCHI; + stars[i].crNsigma = CRN; + stars[i].extNsigma = EXTN; + stars[i].fx = FX; + stars[i].fy = FY; + stars[i].df = DF; + stars[i].psfQual = PSFQUAL; + stars[i].nFrames = 1; + stars[i].flags = FLAGS; + } + + // create primary header + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + gfits_print (&header, "NEXTEND", "%d", 1, 1); + + // XXX add minimum needed header fields + gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax + 50)); + gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax + 50)); + + gfits_modify (&header, "NSTARS", "%d", 1, Nstars); + gfits_modify (&header, "PHOTCODE", "%s", 1, "SIMTEST.r.Chip"); + gfits_modify (&header, "DATE-OBS", "%s", 1, date); + gfits_modify (&header, "UTC-OBS", "%s", 1, time); + gfits_modify (&header, "ZERO_PT", "%lf", 1, 25.0); + gfits_modify (&header, "EXPTIME", "%lf", 1, 1.0); + gfits_modify (&header, "AIRMASS", "%lf", 1, 1.0); + gfits_modify (&header, "NASTRO", "%d", 1, 10); + + gfits_modify (&header, "IMAGEID", "%d", 1, (int)(1000*drand48())); + gfits_modify (&header, "SOURCEID", "%d", 1, (int)(100*drand48())); + + /* bore site center guess */ + strcpy (coords.ctype, "DEC--TAN"); + coords.crval1 = RA; + coords.crval2 = DEC; + coords.crpix1 = 0.0; + coords.crpix2 = 0.0; + + coords.cdelt1 = 0.25/3600.0; + coords.cdelt2 = 0.25/3600.0; + + coords.pc1_1 = 1; + coords.pc2_2 = 1; + coords.pc1_2 = 0; + coords.pc2_1 = 0; + coords.Npolyterms = 1; + + PutCoords (&coords, &header); + gfits_modify (&header, "EXTNAME", "%s", 1, "Chip.hdr"); + + ftable.header = &theader; + gfits_table_set_PS1_DEV_1 (&ftable, stars, Nstars); + gfits_modify (&theader, "EXTNAME", "%s", 1, "Chip.psf"); + gfits_modify (&theader, "EXTTYPE", "%s", 1, "PS1_DEV_1"); + gfits_modify (&theader, "EXTHEAD", "%s", 1, "Chip.hdr"); + + fits = fopen (argv[2], "w"); + if (fits == NULL) { + fprintf (stderr, "ERROR: can't open output file %s\n", argv[2]); + exit (1); + } + + gfits_fwrite_header (fits, &header); + gfits_fwrite_matrix (fits, &matrix); + gfits_fwrite_Theader (fits, &theader); + gfits_fwrite_table (fits, &ftable); + fclose (fits); + + exit (0); +} + + +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_branch_20081124/Ohana/src/addstar/test/addstar.sh =================================================================== --- /branches/eam_branch_20081124/Ohana/src/addstar/test/addstar.sh (revision 20905) +++ /branches/eam_branch_20081124/Ohana/src/addstar/test/addstar.sh (revision 20905) @@ -0,0 +1,18 @@ + +macro mkinput + exec rm -f test.in.txt + + output test.in.txt + for i 10 1024 100 + for j 10 1024 100 + fprintf " %4d %4d %6.2f" $i $j {-15.0 + 2.5*($i + $j)/1000.0} + end + end + output stdout +end + +macro showinput + data test.in.txt + read x 1 y 2 m 3 + lim x y; clear; box; zplot -pt 7 x y m -10 -15 +end