Index: /tags/gophot-1-0/Ohana/src/gophot/Makefile =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/Makefile (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/Makefile (revision 2556) @@ -0,0 +1,70 @@ +include ../../Configure +HOME = $(ROOT)/src/gophot +PROGRAM = gophot + +CC = gcc -g +BIN = $(HOME)/bin +INC = $(HOME)/include +SRC = $(HOME)/src +MAN = $(HOME)/doc +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +# +INCS = -I$(INC) -I$(LINC) -I$(XINC) +LIBS = -L$(LLIB) -lFITS -lohana -lm +CFLAGS = -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +default: $(PROGRAM) + +OBJ = \ + $(SRC)/gophot.$(ARCH).o $(SRC)/ConfigInit.$(ARCH).o \ + $(SRC)/message.$(ARCH).o $(SRC)/dophot.$(ARCH).o \ + $(SRC)/makenoise.$(ARCH).o $(SRC)/findsky.$(ARCH).o \ + $(SRC)/makemask.$(ARCH).o $(SRC)/paravg.$(ARCH).o \ + $(SRC)/objout.$(ARCH).o $(SRC)/addstar.$(ARCH).o \ + $(SRC)/ellipse.$(ARCH).o $(SRC)/improve.$(ARCH).o \ + $(SRC)/fillerup.$(ARCH).o $(SRC)/parinterp.$(ARCH).o \ + $(SRC)/guess.$(ARCH).o $(SRC)/offpic.$(ARCH).o \ + $(SRC)/addlims.$(ARCH).o $(SRC)/parupd.$(ARCH).o \ + $(SRC)/toofaint.$(ARCH).o $(SRC)/isearch.$(ARCH).o \ + $(SRC)/shape.$(ARCH).o $(SRC)/galaxy.$(ARCH).o \ + $(SRC)/transmask.$(ARCH).o $(SRC)/cosmic.$(ARCH).o \ + $(SRC)/impaper2.$(ARCH).o $(SRC)/sort.$(ARCH).o \ + $(SRC)/oblit.$(ARCH).o $(SRC)/toobright.$(ARCH).o \ + $(SRC)/Chisq.$(ARCH).o $(SRC)/oblims.$(ARCH).o \ + $(SRC)/ludcmp.$(ARCH).o $(SRC)/lubksb.$(ARCH).o \ + $(SRC)/pseud2d.$(ARCH).o $(SRC)/twofit.$(ARCH).o \ + $(SRC)/skyfun_plane.$(ARCH).o $(SRC)/pseud4d.$(ARCH).o \ + $(SRC)/large_features.$(ARCH).o $(SRC)/set_thresholds.$(ARCH).o \ + $(SRC)/outline.$(ARCH).o $(SRC)/delete_ellipse.$(ARCH).o \ + $(SRC)/feature_fluxes.$(ARCH).o $(SRC)/mediansky.$(ARCH).o + +# dependancy rules for binary code ########################## +$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + +$(BIN)/$(PROGRAM).$(ARCH): $(OBJ) + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi + $(CC) $(OBJ) -o $(BIN)/$(PROGRAM).$(ARCH) $(CCFLAGS) + +install: $(DESTBIN)/$(PROGRAM) + +$(DESTBIN)/$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$(PROGRAM) + cp $(BIN)/$(PROGRAM).$(ARCH) $(DESTBIN)/$(PROGRAM) + +# utilities ################################################# +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $< Index: /tags/gophot-1-0/Ohana/src/gophot/bin/.cvsignore =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/bin/.cvsignore (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/bin/.cvsignore (revision 2556) @@ -0,0 +1,1 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix Index: /tags/gophot-1-0/Ohana/src/gophot/doc/Notes =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/doc/Notes (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/doc/Notes (revision 2556) @@ -0,0 +1,144 @@ + +Converting dophot to C: (I use dophot.f and dophot.c to refer to the +fortran and C versions of dophot). + +some rules and goals for the first pass: + +1) I am maintaining the fortran dophot variable names as much as + possible, at least for now. + +2) I am avoiding making proceedural changes / method changes, but I am +willing to make minor modifications of things like looping strategies +and if-statement order to clean the code a bit. + +3) I have placed all of the common-block variables in a single include +file, gophot.h, with some attempt to define the names. I've also placed the +image data and noise matricies in gophot.h, along with nfast, nslow, +so these are not passed endlessly back and forth among functions. + +4) My (sort-of) working fortran f2c dophot used an external loop, +getfits.c which enclosed dophot.f and performed the image loading, +big matrix allocation stuff. I have renamed that gophot.c and cleaned +it somewhat. + +5) The biggest structural change so far is the infamous tuneup.f. I +have replaced it with my own ConfigInit function, based on my ohana +config functions. These have a much cleaner syntax for reading +configuration values. I have also placed the ConfigInit function much +earlier, in gophot.c at the very beginning. this means that config +values can be overridden by image header values, if the fits_scan +lines are coded. currently, itop is read from the header (SATVALUE). + +For the moment, I am only making those changes needed to get a working +dophot version, not a complete version. Thus, I am skipping: +warmstart, autothresh, median sky, hubble sky. + +some function with significant adjustments, other than getfits.c and +tuneup.f: + +message.c: actually this didn't exist in dophot.f - I have created a +function mprint, which conditionally prints a message based on the +verbosity level and an associated threshold. the verbosity is +registered with set_verbosity (level) and is stored static in +message.c (not visible to the rest of dophot.c). the prototype is: +mprint (int level, char format, ...), with a format and argument +syntax identical to printf. all messages go to stderr, for now. + +makenoise.c: I am using pointers to do the loops. this should speed +life alot. + +findsky.c: a new addition of mine to determine skyvalue from median of +image. this is implemented in the most recent dophot.f, but only +since March 2000. + +There are some typical coding structures in the fortran that I am +changing to make more legible: + +1) a = a + 1, a = a + b --> I'm using C inc args: a++, a+=b + +2) nested logic: + +{ /* some block */ + + code, code, code; + + if (condition) then + + big block; + + endif + +} /* end block */ + +there are many of these structures which conditionally perform a large +operation, and otherwise do nothing or almost nothing. I'm inverting +these like this: + +{ /* some block */ + + code, code, code; + + if (!condition) continue + + big block; + +} /* end block */ + +data hard.soft.mat +read m1 4 m2 15 +set dm = m1 - m2 +lim -17 -7 -1 1; clear; box; plot m1 dm + +I have been making some significant changes to the dophot strategy for +finding stars. The goal is to handle saturated stars more +effectively. The strategy should not really change for faint stars. +There are a few problems with the old strategy. + +First, we were using an image array and a noise array. To inhibit +finding fake stars in the wings of bright stars, the noise array was +enhanced when the stars were subtracted. But for each loop, we were +fitting and re-fitting each star several times. As a result, we were +adding, subtracting, adding, and subtracting the noise, which was of +the order Nstar^2, many times. Round off errors eventually made the +noise array useless in the vicinity of bright stars. Second, bright +stars likely have many bad pixels, either low from sagging counts at +or near saturation or high from bleeding. Since these all have quite +high flux, they highly discrepant points drive the fits +significantly. Finally, the initial guess location used the peak of +the flux in a box around the trigger pixel. This was problematic for +bright stars because bleeding and saturation meant the peak was +frequently far from the centroid. I am addressing each of these +problems separately. + +First, for the fits, we are now softening the errors by the residual +from the initial guess fit. This is a fitting method called 'robust +fitting', and it allows significantly discrepant points to be part of +the data without dominating or driving the fit to funny places. This +means the saturated stars are less biased by the few pixels in the +core which are saturated and discrepant from the fit for the other +pixels. + +Second, for the noise problem, I have made several related changes. +First, I am now converting the image from DN to electrons (and +re-scaling all relevant parameters to match). This means the Poisson +noise in the original image is now simply sqrt(Ncts). I am putting +only the square of the read noise in the noise array. When I subtract +a star, I move the counts from the image array to the noise array. +Thus, the noise per pixel is always (Icts + Ncts) [where Icts = image +array counts and Ncts = noise array counts] before or after the star +subtraction. I now have two tests for bump significance: 1) the +signal-to-noise must be above a threshold: + +sq(sum[image-sky]) > sum[noise + image] * ratio + +and the image counts must be significant: + +sub[image-sky] > sum[noise - readnoise^2] * ratio + +Finally, the initial guess will use a centroid in fillerup if the flux +is high enough. + + + + + Index: /tags/gophot-1-0/Ohana/src/gophot/include/fit.h =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/include/fit.h (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/include/fit.h (revision 2556) @@ -0,0 +1,27 @@ + +# ifndef ALLOCATE +# define ALLOCATE(X,T,S) \ + X=(T *)malloc((unsigned) (MAX(((S)*sizeof(T)),1)));\ + if(X==NULL) \ + { \ + fprintf(stderr,"failed to malloc X\n");\ + exit (10);\ + } +# define REALLOCATE(X,T,S) \ + X=(T *)realloc(X,(unsigned) (MAX(((S)*sizeof(T)),1))); \ + if(X==NULL) \ + { \ + fprintf(stderr,"failed to realloc X\n"); \ + exit (10); \ + } +# endif /* ALLOCATE */ + +float mrq2dinit (int *, int *, float *, float *, int, float *, int, float (funcs)(int *, int *, float *, float *)); +float mrq2dmin (int *, int *, float *, float *, int, float *, int, float (funcs)(int *, int *, float *, float *)); +float **mrq2dcovar (int); + +# define NPARS 8 +# define MIN(X,Y) ((X) < (Y) ? (X) : (Y)) +# define MAX(X,Y) ((X) > (Y) ? (X) : (Y)) +# define SQ(X) (double) (((double)(X))*((double)(X))) +# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} Index: /tags/gophot-1-0/Ohana/src/gophot/include/gophot.h =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/include/gophot.h (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/include/gophot.h (revision 2556) @@ -0,0 +1,128 @@ +/* includes */ + +# include + +typedef char bool; + +# define TRUE (1) +# define FALSE (0) +# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) +# define ROUND(X) ((int) ((X) + 0.5*SIGN(X))) +# define SQ(X) (double) (((double)(X))*((double)(X))) +# define MIN(X,Y) ((X) < (Y) ? (X) : (Y)) +# define MAX(X,Y) ((X) > (Y) ? (X) : (Y)) +# define SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} + +/* constants */ + +# define NRMAX 1024 /* max Ny, deprecated */ +# define NCMAX 1024 /* max Nx, deprecated */ + +# define CHIPAR 0.9 /* used by chisq.f */ +# define NSMAX 100000 /* max number of stars */ +# define NPMAX 8 /* max number of parameters */ +# define NPAR 8 /* max number of parameters used */ +# define NPSKY 8 /* max sky fit parameters */ +# define NSKYFIT 3 /* max sky fit parameters used */ +# define NFF 20 /* max number of files / flags */ +# define NAPPLE 5 /* number of aper data somethings? */ +# define NAPMAX 30 /* number of correction apertures in file? */ +# define NFIT0 2 /* N par in fit 0 */ +# define NFIT1 4 /* N par in fit 1 */ +# define NFIT2 7 /* N par in fit 2 */ +# define NFIT3 8 /* N par in fit 3 */ + +# define MAXFIL 5000 /* max size of subraster vector */ +# define NMASK 17 /* max mask size */ +# define MAGIC HUGE_VAL /* sentinel for bad pixels */ +/* # define MAGICSET 2e30 old value for sentinel */ + +# define ADD +1 +# define SUB -1 + +/* global variables */ + +/* int lverb; * verbosity - no longer global */ +float chipar; /* unknown chisq scale factor */ +float ufactor; /* star scaling factor */ + +/* float b[2*NPMAX]; * two-star fit array */ +/* float fb[2*NPMAX]; * two-star fit error array */ + +bool test7; +bool needit; /* deprecated? */ + +int xs[MAXFIL], ys[MAXFIL], nrect[3]; /* subraster vectors */ +float zs[MAXFIL], dzs[MAXFIL]; /* subraster vectors */ +float ts[MAXFIL]; + +float starmask[NMASK][NMASK]; +float parms[NPMAX]; + +/* float a[NPMAX], fa[NPMAX], c[NPMAX][NPMAX]; fit param arrays */ +float chiimp, apertime, filltime, addtime; /* deprecated? */ + +float sum0, sum1, sum2, maxval, xmax, ymax, xmax2, ymax2; /* crude star statistics */ +int npt; + +float chi[5]; /* almost deprecated, but still in galaxy & shape */ + +/* tuneup parameters */ + +enum {NONE1, PGAUSS}; +enum {NONE2, PLANE, HUBBLE, MEDIAN}; +enum {NONE3, COMPLETE, INCOMPLETE, INTERNAL, OLDSTYLE}; + +char flags[NFF]; +char files[64][NFF]; +bool fixpos; + +float skyguess, tmin, tmax, tfac; +float fac, xpnd, ctpersat, widobl, cmax; +float stograt, discrim, sig[4], arect[3]; +float chicrit, xtra, crit7, snlim, bumpcrit, sn2cos; +float enuff4, enuff7; +float eperdn, rnoise; +float acc[NPMAX], parlim[NPMAX], ava[NPMAX]; +float beta4, beta64; +float pixthresh; +float apmagmaxerr; +float nphsub, nphob, apmax, apskymin, apskymax, aperrmax; + +int irect[3], krect[3], ibot, itop, nit, grect[3]; +int icrit, ixby2, iyby2; +int n0left, n0right, nthpix, nbadleft, nbadright, nbadtop, nbadbot; +int jhxwid, jhywid, mprec, napertures; + +/* image data */ +float *big, *noise; +int nfast, nslow; /* NAXIS1, NAXIS2 of image */ + +/* star data */ +float starpar[NSMAX][NPMAX]; +float galpar[NSMAX][NPMAX]; +float shadow[NSMAX][NPMAX]; +float shaderr[NSMAX][NPMAX]; +float apple[NSMAX][NAPPLE]; +int imtype[NSMAX]; +int nstot; +float thresh; +float probgal[NSMAX]; /* deprecated */ +float rchisq[NSMAX]; /* deprecated */ +bool fixxy; + +int nregion; +float region[100][8]; + +/* sky data */ +float skypar[NPSKY]; + +/* image data */ +Header header; +Matrix matrix; + +float (*onestar)(int, int, float *, float *); +float (*twostar)(int, int, float *, float *); +float (*skyfun)(int, int, float *, float *); + +# include "prototypes.h" Index: /tags/gophot-1-0/Ohana/src/gophot/include/prototypes.h =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/include/prototypes.h (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/include/prototypes.h (revision 2556) @@ -0,0 +1,46 @@ +float pseud2d (int ix, int iy, float *a, float *fa); +float pseud4d (int ix, int iy, float *a, float *fa); +float skyfun_plane (int ix, int iy, float *a, float *fa); +float chisq (float (function)(int, int, float *, float *), int *, int *, float *, float *, int, float *, float*, int, float *, float *, int); +int ConfigInit (int *argc, char **argv); +int addlims (float *star, int *jrect); +int addstar (float *instar, int iadd, int type); +char *SelectConfigFile (int *argc, char **argv, char *progname); +char *LoadConfigFile (char *filename); +char *ScanConfig (char *config, char *field, char *mode, int Nentry, ...); +char *expandline (char *line, char *config); +bool cosmic (float *star); +int dophot (); +int ellipse (float, float, float, float *, float *, float *, float *); +int errupd (float *, int); +int fillerup (int, int, int); +float findsky (); +bool galaxy (float *, float *, float *); +float guess1 (float *, float *, int, int); +float guess2 (float *, float *, int *, int *); +float guess3 (float *, float *, int *, int *); +int impaper2 (int); +int improve (int); +int isearch (); +void lubksb (float **, int, int *, float *); +void ludcmp (float **, int, int *, float *); +int makemask (); +int makenoise (); +int mprint (int, char *, ...); +int set_verbosity (int); +int completeout (); +int oblims (float *, int *); +bool oblit (float *); +bool oblit (float *); +int paravg (); +float parinterp (float, float, float *); +int parupd (float *, float *, int, int); +int shape (); +int sort (float *, int); +bool toobright (float *); +bool toofaint (float *, float *); +bool transmask (int, int, float); +float twofit (float *, float *, float *, float *); +int variparplane (); +bool offpic (float *, int, int, float *, float *); +float newguess (float *, float *, int, int); Index: /tags/gophot-1-0/Ohana/src/gophot/src/Chisq.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/Chisq.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/Chisq.c (revision 2556) @@ -0,0 +1,202 @@ +# include "gophot.h" + +float chisq (function, ix, iy, z, dz, npts, pars, dpars, npars, Acc, Parlim, niter) + float (function)(int, int, float *, float *); + int *ix; + int *iy; + float *z, *dz; + float *pars, *dpars, *Acc, *Parlim; + int npts, npars, niter; +{ + + float **covmatr, **tmpmatr; + float *v, *tmpvec, *fpars; + int *indx; + bool conv, marq, limit, islimit; + int i, j, k, jj, kk; + float ifact, fact, chinew, chiold, f, dz1, fakk, d, tmpval, value; + float perdeg, save; + + if (npars >= NPMAX) return (MAGIC); + + ALLOCATE (v, float, npars); + ALLOCATE (indx, int, npars); + ALLOCATE (fpars, float, NPMAX); + ALLOCATE (tmpvec, float, npars); + ALLOCATE (covmatr, float *, npars); + ALLOCATE (tmpmatr, float *, npars); + for (i = 0; i < npars; i++) { + ALLOCATE (covmatr[i], float, npars); + ALLOCATE (tmpmatr[i], float, npars + 1); + } + + conv = FALSE; + limit = FALSE; + + /* check if parameters exceeds limits at the start. */ + for (j = 0; j < npars; j++) { + if (Parlim[j] < 0.0) { + limit = fabs(pars[j]) > fabs(Parlim[j]); + if (limit) mprint (4, "self-deception has occured: initial limits\n"); + } + } + + /* + for (j = 0; j < npars; j++) { + mprint (3, "%d %f %f\n", j, pars[j], fpars[j]); + } + */ + + ifact = 0; + + for (i = 0; (i < niter) && !conv && !limit; i++) { + chinew = 0.0; + for (j = 0; j < npars; j++) { + for (kk = 0; kk < npars + 1; kk++) { + tmpmatr[j][kk] = 0; + } + } + + for (j = 0; j < npts; j++) { + f = function (ix[j], iy[j], pars, fpars) - z[j]; + /* fprintf (stderr, "%d %d %d %f %f %f\n", j, ix[j], iy[j], f, z[j], dz[j]); */ + dz1 = 1.0 / dz[j]; + chinew += SQ(f)*dz1; + for (kk = 0; kk < npars; kk++) { + if (fabs(fpars[kk]) > 1e-12) { + fakk = fpars[kk]*dz1; + tmpmatr[kk][npars] += fakk*f; + for (jj = 0; jj <= kk; jj++) { + if (fabs(fpars[kk]) > 1e-12) tmpmatr[kk][jj] += fakk*fpars[jj]; + } + } + } + } + + chiold = chinew; + marq = FALSE; + for (k = 1; (k <= 10) && !marq && !limit; k++) { + conv = (k == 1); + fact = (k == 1) ? 0.0 : pow (2.0, ifact); + for (j = 0; j < npars; j++) { + for (jj = 0; jj < j; jj++) { + covmatr[j][jj] = tmpmatr[j][jj]; + covmatr[jj][j] = tmpmatr[j][jj]; + } + covmatr[j][j] = (1+fact)*tmpmatr[j][j]; + v[j] = tmpmatr[j][npars]; + } + ludcmp (covmatr, npars, indx, &d); + + /* if d = 0, the matrix was singular; no convergence. */ + if (d == 0) { + mprint (4, "singular matrix!\n"); + /* need to free arrays */ + return (MAGIC); + } + lubksb (covmatr, npars, indx, v); + + /* + check if change in parameters exceeds limits. if Parlim(j) > 0, then + consider fractional changes. if Parlim(j) < 0, consider absolute + changes. if Parlim(j) = 0, ignore this test. + */ + for (j = 0; j < npars; j++) { + pars[j] -= chipar*v[j]; + if (Parlim[j] > 0.0) { + tmpval = fabs (v[j]/pars[j]); + islimit = (tmpval > Parlim[j]); + limit = limit || islimit; + if (islimit) mprint (4, "self-deception has occured: frac limits: %d %f %f %f\n", j, pars[j], v[j], Parlim[j]); + } + if (Parlim[j] < 0.0) { + islimit = (fabs (pars[j]) > fabs (Parlim[j])); + limit = limit || islimit; + if (islimit) mprint (4, "self-deception has occured: abs limits: %d %f %f %f\n", j, pars[j], v[j], Parlim[j]); + } + /* check convergence */ + if (Acc[j] > 0) { + tmpval = fabs (v[j]/pars[j]); + conv = conv && (tmpval <= Acc[j]); + } else { + conv = conv && (fabs (v[j]) <= fabs (Acc[j])); + } + } + + if (conv) { + marq = TRUE; + } else { + if (!limit) { + chinew = 0.0; + for (j = 0; j < npts; j++) { + f = function(ix[j], iy[j], pars, (float *) NULL) - z[j]; + chinew += SQ(f)/dz[j]; + } + if (k == 2) ifact --; + if (chinew < 1.0001*chiold) { + marq = TRUE; + } else { + if (k == 2) ifact += 2; + if (k > 2) ifact ++; + if (ifact > 10) goto escape; + for (j = 0; j < npars; j++) pars[j] += chipar*v[j]; + } + } + } + + mprint (4, "%d, %d, ", i, k); + for (kk = 0; kk < npars; kk++) mprint (4, "%f ", pars[kk]); + mprint (4, " %f\n", chinew); + + } + } + escape: + + if (!limit) { + for (j = 0; j < npars; j++) { + for (i = 0; i < npars; i++) tmpvec[i] = 0; + tmpvec[j] = 1; + lubksb (covmatr, npars, indx, tmpvec); + for (i = 0; i < npars; i++) tmpmatr[i][j] = tmpvec[i]; + } + } + + if (conv && !limit) { + perdeg = sqrt (chinew / MAX (npts - npars, 1)); + for (i = 0; i < npars; i++) { + if (tmpmatr[i][i] > 0) { + save = sqrt (tmpmatr[i][i]); + } else { + mprint (4, "trouble: negative autovariance for tmpmatr[%d][%d] = %f\n", i, i, tmpmatr[i][i]); + save = 1e10; + } + for (j = 0; j < npars; j++) { + covmatr[i][j] = tmpmatr[i][j]/save; + covmatr[j][i] = tmpmatr[j][i]/save; + } + covmatr[i][i] = save*perdeg; + } + value = chiold; + } else { + value = MAGIC; + } + + for (i = 0; i < npars; i++) dpars[i] = SQ(covmatr[i][i]); + + for (kk = 0; kk < npars; kk++) mprint (4, "%f ", pars[kk]); + mprint (4, " %f\n", chiold); + + for (i = 0; i < npars; i++) { + free (covmatr[i]); + free (tmpmatr[i]); + } + free (v); + free (indx); + free (fpars); + free (tmpvec); + free (covmatr); + free (tmpmatr); + + return (value); + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/ConfigInit.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/ConfigInit.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/ConfigInit.c (revision 2556) @@ -0,0 +1,264 @@ +# include "gophot.h" + +# define TestConfig(A,B,C,D,E) { if (!ScanConfig (A,B,C,D,E)) { fprintf (stderr, B); exit (1); }} + +ConfigInit (int *argc, char **argv) { + + char *config, *file; + float fwhm, ar, tilt; + char line[128]; + int level; + float gmajwid, gxwid, gywid, fwhmx, fwhmy; + float scalefb, fbmin, scaleab, abmin, scalemb, ambmin; + + /*** load configuration info ***/ + file = argv[3]; + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + exit (0); + } + + set_verbosity (0); /* initializes function data */ + + /* flags -- OK for now, not really used much */ + flags[1] = PGAUSS; /* psf type */ + flags[2] = NONE2; /* sky type */ + flags[3] = NONE3; /* objtype_out */ + flags[4] = FALSE; /* output shadow file? */ + flags[5] = FALSE; /* image out? */ + flags[6] = FALSE; /* warmstart? */ + flags[7] = FALSE; /* load input shadow file? */ + flags[8] = FALSE; /* objtype_in */ + flags[10] = FALSE; /* aperture correction file? */ + + onestar = pseud2d; + twostar = pseud4d; + skyfun = skyfun_plane; + + /* files: + files[2] - image_out + files[3] - objects_in - UNUSED + files[5] - shadowfile_out - UNUSED + files[7] - shadowfile_in - UNUSED + files[6] - logfile - UNUSED + files[8] - apcorrfile - UNUSED + */ + + /** legacy - delete? */ + n0left = 0; + n0right = 0; + nthpix = 0; + + TestConfig (config, "FWHM", "%f", 0, &fwhm); /* Approx FWHM of objects (pixels) along major axis. */ + TestConfig (config, "AXIS_RATIO", "%f", 0, &ar); /* For star objects. AR=b/a; b=minor axis. */ + TestConfig (config, "TILT", "%f", 0, &tilt); /* Angle of major axis in degrees; +x=0; +y=90. */ + TestConfig (config, "SKY", "%f", 0, &skyguess); /* Approximate mean sky value in data numbers. */ + + /* convert to internal values */ + fwhm *= 1.2; /* is this really necessary? */ + tilt = tilt/57.29578; + gmajwid = SQ(fwhm/2.3548); + gxwid = gmajwid*(SQ(cos(tilt)) + SQ(ar*sin(tilt))); + gywid = gmajwid*(SQ(ar*cos(tilt)) + SQ(sin(tilt))); + fwhmx = 2.3548*sqrt(gxwid); + fwhmy = 2.3548*sqrt(gywid); + mprint (1, "fwhm x & y: %f, %f\n", fwhmx, fwhmy); + + /* no longer use ava[0-3] */ + ava[4] = gxwid; + ava[6] = gywid; + ava[5] = 0.01/sqrt(gxwid*gywid); + + TestConfig (config, "NFITBOX_X", "%d", 0, &irect[1]); /* Size of fit box in the x-direction. */ + TestConfig (config, "NFITBOX_Y", "%d", 0, &irect[2]); /* Size of fit box in the y-direction. */ + TestConfig (config, "MASKBOX_X", "%d", 0, &ixby2); /* Size of mask box size in x. */ + TestConfig (config, "MASKBOX_Y", "%d", 0, &iyby2); /* Size of mask box size in y. */ + TestConfig (config, "APBOX_X", "%f", 0, &arect[1]); /* Size of aperture photometry box in x. */ + TestConfig (config, "APBOX_Y", "%f", 0, &arect[2]); /* Size of aperture photometry box in y. */ + TestConfig (config, "NGALBOX_X", "%d", 0, &grect[1]); /* Size of fit box in the x-direction. */ + TestConfig (config, "NGALBOX_Y", "%d", 0, &grect[2]); /* Size of fit box in the y-direction. */ + + TestConfig (config, "IBOTTOM", "%d", 0, &ibot); /* Lowest allowed data value in data numbers. */ + TestConfig (config, "ITOP", "%d", 0, &itop); /* Level where Saturation begins. */ + TestConfig (config, "THRESHMIN", "%f", 0, &tmin); /* Sigmas above sky for min threshold */ + TestConfig (config, "THRESHMAX", "%f", 0, &tmax); /* Value of maximum threshold. */ + TestConfig (config, "THRESHDEC", "%f", 0, &tfac); /* Threshold decrement in powers-of-2. */ + TestConfig (config, "EPERDN", "%f", 0, &eperdn); /* Electrons per data number. */ + TestConfig (config, "RDNOISE", "%f", 0, &rnoise); /* Readout noise in electrons. */ + rnoise = SQ (rnoise); /* we will store SQ(rnoise) since this is always needed */ + + /* I don't like AUTOSCALE because it is not dynamic - + it sets the values once here for the run, so FWHM better be right */ + TestConfig (config, "AUTOSCALE", "%s", 0, line); /* Auto-scaling of sizes by FWHM. */ + if (!strncasecmp (line, "y", 1)) { + ScanConfig (config, "SCALEFITBOX", "%f", 0, &scalefb); /* Size of fit box in units of FWHM. */ + ScanConfig (config, "FITBOXMIN", "%f", 0, &fbmin); /* Smallest allowed fit box size. */ + ScanConfig (config, "SCALEAPBOX", "%f", 0, &scaleab); /* Size of aperture phot box in units of FWHM. */ + ScanConfig (config, "APBOXMIN", "%f", 0, &abmin); /* Smallest allowed aperture phot box size. */ + ScanConfig (config, "SCALEMASKBOX", "%f", 0, &scalemb); /* Size of mask box in units of FWHM. */ + ScanConfig (config, "AMASKBOXMIN", "%f", 0, &ambmin); /* Smallest allowed mask box size. */ + irect[1] = MAX (fwhmx*scalefb, fbmin); + irect[2] = MAX (fwhmy*scalefb, fbmin); + arect[1] = MAX (fwhmx*scaleab, abmin); + arect[2] = MAX (fwhmy*scaleab, abmin); + ixby2 = MAX(fwhmx*scalemb, ambmin); + ixby2 = MAX(fwhmy*scalemb, ambmin); + } + /* force boxes to have odd sizes */ + if (((int)arect[1]) % 2 == 0) arect[1]++; + if (((int)arect[2]) % 2 == 0) arect[2]++; + if (irect[1] % 2 == 0) irect[1]++; + if (irect[2] % 2 == 0) irect[2]++; + if (ixby2 % 2 == 0) ixby2 ++; + if (iyby2 % 2 == 0) iyby2 ++; + ixby2 = (ixby2 - 1)/2; + iyby2 = (iyby2 - 1)/2; + + fixpos = FALSE; + TestConfig (config, "FIXPOS", "%s", 0, line); /* Fix star positions? */ + if (!strncasecmp (line, "y", 1)) fixpos = TRUE; + + if (ScanConfig (config, "IMAGE_OUT", "%s", 0, files[2])) flags[5] = TRUE; /* Output image name. */ + + TestConfig (config, "OBJTYPE_OUT", "%s", 0, &line); /* Output format: (COMPLETE, INCOMPLETE, INTERNAL) */ + if (!strcasecmp (line, "complete")) flags[3] = COMPLETE; + if (!strcasecmp (line, "incomplete")) flags[3] = INCOMPLETE; + if (!strcasecmp (line, "internal")) flags[3] = INTERNAL; + if (!strcasecmp (line, "oldstyle")) flags[3] = OLDSTYLE; + if (flags[3] == NONE3) { + fprintf (stderr, "invalid OBJTYPE_OUT: %s\n", line); + exit (1); + } + + TestConfig (config, "LOGVERBOSITY", "%d", 0, &level); /* Verbosity of log file; (0-4). */ + set_verbosity (level); + TestConfig (config, "RESIDNOISE", "%f", 0, &fac); /* Fraction of noise to ADD to noise file. */ + TestConfig (config, "FOOTPRINT_NOISE", "%f", 0, &xpnd); /* Expand stars in noise file by this amount. */ + TestConfig (config, "NPHSUB", "%f", 0, &nphsub); /* Limiting surface brightness for subtractions. */ + TestConfig (config, "NPHOB", "%f", 0, &nphob); /* Limiting surface brightness for obliterations. */ + TestConfig (config, "ICRIT", "%d", 0, &icrit); /* Obliterate if # of pixels > ITOP exceeds this. */ + TestConfig (config, "CENTINTMAX", "%f", 0, &cmax); /* Obliterate if central intensity exceeds this. */ + TestConfig (config, "CTPERSAT", "%f", 0, &ctpersat); /* Assumed intensity for saturated pixels. */ + + TestConfig (config, "STARGALKNOB", "%f", 0, &stograt); /* Star/galaxy discriminator: bigger number, more stars */ + TestConfig (config, "STARCOSKNOB", "%f", 0, &discrim); /* Object/cosmic-ray discriminator: bigger number, more cosmics */ + TestConfig (config, "SNLIM7", "%f", 0, &crit7); /* Minimum S/N for 7-parameter fit. */ + crit7 = SQ(crit7); + TestConfig (config, "SNLIM", "%f", 0, &snlim); /* Minimum S/N for a pixel to be in fit subraster. */ + TestConfig (config, "SNLIMMASK", "%f", 0, &bumpcrit); /* Minimum S/N through mask to identify an object. */ + TestConfig (config, "SNLIMCOS", "%f", 0, &sn2cos); /* Minimum S/N to be called a cosmic ray. */ + sn2cos = SQ(sn2cos); + TestConfig (config, "NBADLEFT", "%d", 0, &nbadleft); /* Ignore pixels closer to the left edge than this. */ + TestConfig (config, "NBADRIGHT", "%d", 0, &nbadright); /* Ignore pixels closer to the right edge than this. */ + TestConfig (config, "NBADTOP", "%d", 0, &nbadtop); /* Ignore pixels closer to the top edge than this. */ + TestConfig (config, "NBADBOT", "%d", 0, &nbadbot); /* Ignore pixels closer to the bottom edge than this. */ + + TestConfig (config, "SKYTYPE", "%s", 0, line); /* SKY type: (PLANE, HUBBLE, MEDIAN) */ + if (!strcasecmp (line, "plane")) flags[2] = PLANE; + /* if (!strcasecmp (line, "hubble")) flags[2] = HUBBLE; */ + /* if (!strcasecmp (line, "median")) flags[2] = MEDIAN; */ + if (flags[2] == NONE2) { + fprintf (stderr, "invalid SKYTYPE: %s\n", line); + exit (1); + } + + TestConfig (config, "NFITITER", "%d", 0, &nit); /* Maximum number of iterations. */ + TestConfig (config, "NFITBOXFIRST_X", "%d", 0, &krect[1]); /* Size of fit box in x for first pass. */ + TestConfig (config, "NFITBOXFIRST_Y", "%d", 0, &krect[2]); /* Size of fit box in y for first pass. */ + TestConfig (config, "CHI2MINBIG", "%f", 0, &chicrit); /* Critical CHI-squared for a large object. */ + TestConfig (config, "XTRA", "%f", 0, &xtra); /* We need more S/N if some pixels are missing. */ + TestConfig (config, "SIGMA1", "%f", 0, &sig[1]); /* Max. frac. scatter in sigma_x for stars. */ + TestConfig (config, "SIGMA2", "%f", 0, &sig[2]); /* Max. scatter in xy cross term for stars. */ + TestConfig (config, "SIGMA3", "%f", 0, &sig[3]); /* Max. frac. scatter in sigma_y for stars. */ + TestConfig (config, "ENUFF4", "%f", 0, &enuff4); /* Fraction of pixels needed for 4-param fit. */ + TestConfig (config, "ENUFF7", "%f", 0, &enuff7); /* Fraction of pixels needed for 7-param fit. */ + TestConfig (config, "COSOBLSIZE", "%f", 0, &widobl); /* Size of obliteration box for a cosmic ray. */ + TestConfig (config, "APMAG_MAXERR", "%f", 0, &apmagmaxerr); /* Max anticipated error for aperture phot report. */ + TestConfig (config, "PIXTHRESH", "%f", 0, &pixthresh); /* Trigger on pixels higher than noise*PIXTHRESH. */ + TestConfig (config, "BETA4", "%f", 0, &beta4); /* R**4 coefficient modifier. */ + TestConfig (config, "BETA6", "%f", 0, &beta64); /* R**6 coefficient modifier. */ + beta4 = 1.0; + beta64 = 1.0; + + TestConfig (config, "RELACC1", "%f", 0, &acc[0]); /* Convergence criterion for sky. */ + TestConfig (config, "RELACC2", "%f", 0, &acc[1]); /* Convergence criterion for for central intensity. */ + TestConfig (config, "RELACC3", "%f", 0, &acc[2]); /* Convergence criterion for x-position. */ + TestConfig (config, "RELACC4", "%f", 0, &acc[3]); /* Convergence criterion for y-position. */ + TestConfig (config, "RELACC5", "%f", 0, &acc[4]); /* Convergence criterion for sigma-x. */ + TestConfig (config, "RELACC6", "%f", 0, &acc[5]); /* Convergence criterion for sigma-xy. */ + TestConfig (config, "RELACC7", "%f", 0, &acc[6]); /* Convergence criterion for sigma-y. */ + TestConfig (config, "PARLIM1", "%f", 0, &parlim[0]); /* Allowed change for sky value. */ + TestConfig (config, "PARLIM2", "%f", 0, &parlim[1]); /* Allowed change for central intensity. */ + TestConfig (config, "PARLIM3", "%f", 0, &parlim[2]); /* Allowed change for x-position. */ + TestConfig (config, "PARLIM4", "%f", 0, &parlim[3]); /* Allowed change for y-position. */ + TestConfig (config, "PARLIM5", "%f", 0, &parlim[4]); /* Allowed change for sigma-x. */ + TestConfig (config, "PARLIM6", "%f", 0, &parlim[5]); /* Allowed change for sigma-xy. */ + TestConfig (config, "PARLIM7", "%f", 0, &parlim[6]); /* Allowed change for sigma-y. */ + + /* other initial values for parameters */ + needit = TRUE; + fixxy = FALSE; + test7 = FALSE; + ufactor = 100; + chipar = 0.9; + + free (config); + /* free (file); */ + +} + + +# if (0) /* things used by MEDIAN sky, disabled */ + /* only if 'median' */ + ScanConfig (config, "JHXWID", "%f", 0, &jhxwid); /* X Half-size of median box (.le. 0 -> autoscale) */ + ScanConfig (config, "JHYWID", "%f", 0, &jhywid); /* Y (same as above) */ + ScanConfig (config, "MPREC", "%f", 0, &mprec); /* Median precision in DN (use .le. 0 for autocalc) */ + ScanConfig (config, "NTHPIX", "%f", 0, &nthpix); /* Frequency of sky updates in pixels for 1st pass */ + + /* log is now going to stderr, up to user to redirect */ + ScanConfig (config, "LOGFILE", "%s", 0, &A); /* Log file name. TERM for screen. */ + ScanConfig (config, "IMAGE_IN", "%s", 0, &A); /* Input image name. */ + ScanConfig (config, "OBJECTS_OUT", "%s", 0, &A); /* Output object list file name. */ + ScanConfig (config, "PARAMS_DEFAULT", "%s", 0, &A); /* Default parameters file name. */ + ScanConfig (config, "PARAMS_OUT", "%s", 0, &A); /* Output parameters file name. */ + ScanConfig (config, "SHADOWFILE_OUT", "%s", 0, &A); /* Output shadow file name. */ + ScanConfig (config, "SHADOWFILE_IN", "%s", 0, &A); /* Input shadow file name. */ + + ScanConfig (config, "ABSLIM1", "%f", 0, &A); /* Allowed range for sky value. */ + ScanConfig (config, "ABSLIM2", "%f", 0, &A); /* Allowed range for central intensity. */ + ScanConfig (config, "ABSLIM3", "%f", 0, &A); /* Allowed range for x-position. */ + ScanConfig (config, "ABSLIM4", "%f", 0, &A); /* Allowed range for y-position. */ + ScanConfig (config, "ABSLIM5", "%f", 0, &A); /* Allowed range for sigma-x. */ + ScanConfig (config, "ABSLIM6", "%f", 0, &A); /* Allowed range for sigma-xy. */ + ScanConfig (config, "ABSLIM7", "%f", 0, &A); /* Allowed range for sigma-y. */ + ScanConfig (config, "ABSLIM8", "%f", 0, &A); /* Allowed range for sigma-y. */ + + ScanConfig (config, "NPARAM", "%f", 0, &A); /* Maximum number of PSF fit parameters. */ + ScanConfig (config, "NFITMAG", "%f", 0, &A); /* No. of PSF parameters to get magnitudes. */ + ScanConfig (config, "NFITSHAPE", "%f", 0, &A); /* No. of PSF parameters to get shape and mags. */ + + ScanConfig (config, "MAXSTARS", "%f", 0, &A); /* Ignore pixels closer to the bottom edge than this. */ + ScanConfig (config, "PSFTYPE", "%s", 0, &A); /* PSF type: (PGAUSS) */ +# endif + +# if (0) /* autothresh eliminated: we are using threshmin = Nsigma above sky always */ + Autothresh = FALSE; + ScanConfig (config, "AUTOTHRESH", "%s", 0, line); /* Auto-scaling of thresholds. */ + if (!strncasecmp (line, "y")) Autothresh = TRUE; + /* only if 'autothresh' */ + ScanConfig (config, "SIGMAIBOTTOM", "%f", 0, &sigbot); /* Level of IBOTTOM below sky in units of noise. */ + ScanConfig (config, "SIGMATHRESHMIN", "%f", 0, &sigthresh); /* Level of THRESHMIN above sky in units of noise. */ +# endif + +# if (0) /* warmstart code */ + status = ScanConfig (config, "OBJECTS_IN", "%s", 0, &A); /* Input object list file name. */ + if (status) { + flags[6] = TRUE; + ScanConfig (config, "OBJTYPE_IN", "%s", 0, &A); /* Input format: (COMPLETE, INTERNAL) */ + } +# endif + +# if (0) /* aperture correction file - there are several other parameters, see tuneup.f */ + ScanConfig (config, "APCORRFILE", "%s", 0, &A); /* Aperture correction file name. */ +# endif Index: /tags/gophot-1-0/Ohana/src/gophot/src/addlims.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/addlims.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/addlims.c (revision 2556) @@ -0,0 +1,34 @@ +# include "gophot.h" + +addlims (float *star, int *jrect) { + + float temp, fudgex, fudgey; + + if (star[1] > 0) { + if (beta4 < 0.1) { + temp = star[1]/nphsub - 1.0; + } else { + temp = pow (6*star[1]/nphsub, 0.33333); + } + if (star[4] > 0) { + fudgex = sqrt(temp*star[4]*2); + } else { + fudgex = 1.5*irect[1]/2; + } + if (star[6] > 0) { + fudgey = sqrt(temp*star[6]*2); + } else { + fudgey = 1.5*irect[2]/2; + } + } else { + fudgex = 1.5*irect[1]/2; + fudgey = 1.5*irect[2]/2; + } + /* don't extrapolate beyond 2*fit box */ + fudgex = MIN (irect[1], fudgex); + fudgey = MIN (irect[2], fudgex); + jrect[1] = star[2] - fudgex; + jrect[2] = star[2] + fudgex; + jrect[3] = star[3] - fudgey; + jrect[4] = star[3] + fudgey; +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/addstar.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/addstar.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/addstar.c (revision 2556) @@ -0,0 +1,60 @@ +# include "gophot.h" + +addstar (float *instar, int iadd, int type) { + + float star0[NPMAX], star1[NPMAX], *bigval, *noiseval; + float sky, bfactor, bsky, cfactor, csky, val; + int ix, iy, jrect[5]; + int i, j, ihi, ilo, jhi, jlo, ixin, iyin; + + if (type == 12) beta4 = 0.01; + + needit = FALSE; + + sky = guess2 (star0, instar, &ixin, &iyin); + sky = guess2 (star1, instar, &ixin, &iyin) / ufactor; + + addlims (instar, jrect); + star1[4] = star0[4]*SQ(xpnd); + star1[5] = star0[5]/SQ(xpnd); + star1[6] = star0[6]*SQ(xpnd); + + bfactor = ufactor*iadd; + bsky = iadd*(0.5 - sky); /* why the 0.5? */ + cfactor = fac*ufactor; + csky = 0.5 - fac*sky; + + ilo = MAX (jrect[1], 0); + ihi = MIN (jrect[2], nfast-1); + jlo = MAX (jrect[3], 0); + jhi = MIN (jrect[4], nslow-1); + + for (j = jlo; j <= jhi; j++) { + iy = j - iyin; + bigval = &big[ilo + j*nfast]; + noiseval = &noise[ilo + j*nfast]; + for (i = ilo; i <= ihi; i++, bigval++, noiseval++) { + if (!finite (*noiseval)) continue; + ix = i - ixin; + /* shouldn't bfactor multiply both onestar and bsky? */ + val = ufactor*(onestar (ix, iy, star0, (float *) NULL) - sky); + *bigval += iadd*val; + + val = fac*ufactor*fabs(onestar (ix, iy, star1, (float *) NULL) - sky); + *noiseval -= iadd*val; + + if (*noiseval <= -1000) { + mprint (2, "i,j,noise = %d, %d, %f, %f %f %f %d, negative noise! obliterating\n", i, j, *noiseval, val, *bigval, cfactor, iadd); + /* + *bigval = MAGIC; + *noiseval = MAGIC; + */ + } + } + } + needit = TRUE; + beta4 = 1.0; +} + +/* star0 is the nominial star region, star1 is the expanded are for the noise array */ + Index: /tags/gophot-1-0/Ohana/src/gophot/src/cosmic.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/cosmic.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/cosmic.c (revision 2556) @@ -0,0 +1,84 @@ +# include "gophot.h" + +bool cosmic (float *star) { + + float dummy[NPMAX]; + float obs, maxob, sky, chistar, chicos, pred, temp; + float sn2, tnoise, *Bval, *Nval; + int ix, iy, I, J, ii, jj, npix, imax, jmax; + bool pointy; + + pointy = FALSE; + + sky = guess2 (dummy, star, &ix, &iy); + temp = big[ix + iy*nfast] / ufactor - dummy[0]; + if (finite (temp)) dummy[1] = temp; + maxob = -HUGE_VAL; + chistar = 0; + chicos = 0; + npix = 0; + + for (J = -1; J <= 1; J++) { + jj = iy + J; + if (jj < 0) continue; + if (jj > nslow - 1) continue; + + Bval = &big[ix - 1 + jj*nfast]; + Nval = &noise[ix - 1 + jj*nfast]; + + for (I = -1; I <= 1; I++, Bval++, Nval++) { + ii = ix + I; + if (ii < 0) continue; + if (ii > nfast - 1) continue; + if (!finite (*Nval)) continue; + npix ++; + pred = ufactor*onestar (I, J, dummy, (float *) NULL); + obs = *Bval; + temp = 1.0 / (*Nval + *Bval); + chistar += SQ (obs - pred) * temp; + sn2 = SQ (obs - sky) * temp; + chicos += sn2; + if ((obs > maxob) && (sn2 >= sn2cos)) { + imax = ii; + jmax = jj; + maxob = obs; + tnoise = temp; + } + } + } + + /* this is meant to test if the object is mostly a single pixel event. + doesn't do a good job of distinguishing a single pixel event from a + generally poor fit to the top of the star. */ + if ((npix >= 7) && (maxob > -HUGE_VAL)) { + chicos -= SQ (maxob - sky) * tnoise; + pointy = (chicos/chistar < discrim); + mprint (3, "location %d, %d, chi-star & chi-cosmic = %f, %f\n", ix, iy, chistar, chicos); + } + + if (pointy) { + mprint (2, "cosmic ray intensity, x, y: %f %d %d\n", maxob, imax, jmax); + star[0] = sky; + star[1] = maxob; + star[2] = imax; + star[3] = jmax; + star[4] = widobl; + star[5] = -1; + star[6] = widobl; + } + return (pointy); +} + +/* We are basically comparing if the object is more like a stellar + object, or if it is compatible with a single high point over a + general sky-like region. + + we have a danger here that a wide object will boost Io + significantly, which in turn makes the fit to the object + particularly poor in the center areas. This in turn makes the + object seem comparable to the sky, and therefore a cosmic. Since + we really only care if the inner 9 pixels are comparable to the + star shape, vaguely, we can fudge this by temporarily setting Io to + the central pixel value. We don't bother if that pixel is bad (NaN). + +*/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/delete_ellipse.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/delete_ellipse.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/delete_ellipse.c (revision 2556) @@ -0,0 +1,81 @@ +# include "gophot.h" + +float delete_ellipse (float *par, float sky) { + + int i, j, npix; + float theta, phi; + float xp, yp, x, y, dx, dy, Dx, Dy; + float Chi, dv, Dv, R2, F2, r1, r2; + float flux, mean; + + npix = flux = 0; + + dx = MAX (par[2], par[3]) + 2; + dy = MAX (par[2], par[3]) + 2; + + for (j = par[1] - dy; j < par[1] + dy; j++) { + if (j < 0) continue; + if (j >= nslow) continue; + for (i = par[0] - dx; i < par[0] + dx; i++) { + if (i < 0) continue; + if (i >= nfast) continue; + + Dx = i - par[0]; + Dy = j - par[1]; + phi = atan2 (Dy, Dx) - RAD_DEG * par[4]; + theta = atan2 (par[2]*sin(phi), par[3]*cos(phi)); + + /* this is the point on the ellipse at the same angle as ref point */ + xp = par[2] * cos (theta); + yp = par[3] * sin (theta); + + x = xp * cos (par[4] * RAD_DEG) - yp * sin (par[4] * RAD_DEG); + y = xp * sin (par[4] * RAD_DEG) + yp * cos (par[4] * RAD_DEG); + + r1 = hypot (Dx, Dy); + r2 = hypot (x, y); + + if (r1 < r2) { + flux += big[i + nfast*j] - sky; + npix ++; + } + + } + } + mean = flux / npix; + fprintf (stderr, "flux: %f, mag: %f, mean: %f\n", flux, -2.5*log10(flux), flux / npix); + + for (j = par[1] - dy; j < par[1] + dy; j++) { + if (j < 0) continue; + if (j >= nslow) continue; + for (i = par[0] - dx; i < par[0] + dx; i++) { + if (i < 0) continue; + if (i >= nfast) continue; + + Dx = i - par[0]; + Dy = j - par[1]; + phi = atan2 (Dy, Dx) - RAD_DEG * par[4]; + theta = atan2 (par[2]*sin(phi), par[3]*cos(phi)); + + /* this is the point on the ellipse at the same angle as ref point */ + xp = par[2] * cos (theta); + yp = par[3] * sin (theta); + + x = xp * cos (par[4] * RAD_DEG) - yp * sin (par[4] * RAD_DEG); + y = xp * sin (par[4] * RAD_DEG) + yp * cos (par[4] * RAD_DEG); + + r1 = hypot (Dx, Dy); + r2 = hypot (x, y); + + if (r1 < r2) { + big[i + nfast*j] -= mean; + } + if (r1 < 1.2*r2) { + noise[i + nfast*j] += fac*(mean + sky); + } + + } + } + fprintf (stderr, "flux: %f, mag: %f, mean: %f\n", flux, -2.5*log10(flux), flux / npix); + return (flux); +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/dophot.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/dophot.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/dophot.c (revision 2556) @@ -0,0 +1,94 @@ +# include "gophot.h" +float set_thresholds (float, float, float *, int *); +void large_features (float, float); +int fix_mediansky (float); +float get_mediansky (int, int); + +dophot () { + + bool first, lastround; + float factor, sky, dsky; + int nstar, Nit, i, j; + char c; + int xtest, ytest; + struct timeval now, then; + + makenoise (); + make_mediansky (); + get_skystats (&sky, &dsky); + fix_mediansky (sky); + + /* + gettimeofday (&then, (void *) NULL); + for (i = 0; i < nfast; i++) { + for (j = 0; j < nslow; j++) { + sky = get_mediansky (i, j); + big[j*nfast + i] -= sky; + } + } + gettimeofday (&now, (void *) NULL); + fprintf (stderr, "elapsed time = %.2f sec\n", + (now.tv_sec - then.tv_sec) + 1e-6*(now.tv_usec - then.tv_usec)); + + fits_write_header ("test.sub", &header); + fits_write_matrix ("test.sub", &matrix); + exit (0); + */ + + thresh = set_thresholds (sky, dsky, &factor, &Nit); + + /* large_features (sky, dsky); */ + + first = TRUE; + lastround = FALSE; + + for (i = 0; i < Nit + 1; i++) { + mprint (0, "starting loop at threshold level %f\n", thresh); + + if (i == Nit) lastround = TRUE; + + makemask (); + + /* fix pos needs to be defined correctly */ + if (fixpos && first) improve (FALSE); + + nstar = isearch (first); + shape (); + paravg (); + + improve (lastround); + + mprint (1, " ending loop at threshold level %f\n", thresh); + mprint (1, " number of new objects found on this threshold = %f\n", nstar); + mprint (1, " total number of objects found so far = %d\n", nstot); + + thresh /= factor; + + first = FALSE; + +# if (0) + fits_write_header ("test.sub", &header); + fits_write_matrix ("test.sub", &matrix); + + completeout (nstot); + fprintf (stderr, "type return to continue: "); + fscanf (stdin, "%c", &c); +# endif + + } + + feature_fluxes (); + + if (lastround) { + if (flags[3] == COMPLETE) completeout (nstot); + } + + /* + fits_write_header ("test.sub", &header); + fits_write_matrix ("test.sub", &matrix); + completeout (nstot); + fprintf (stderr, "type return to continue: "); + fscanf (stdin, "%c", &c); + */ + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/ellipse.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/ellipse.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/ellipse.c (revision 2556) @@ -0,0 +1,57 @@ +# include +# include +# define SQ(X) (double) (((double)(X))*((double)(X))) + +int ellipse (float Sx, float Sxy, float Sy, + float *area, float *amaj, float *amin, float *angle) { + + float A1, A2, A3, R, root1, root2; + + if (fabs(Sxy) >= 1.0 / sqrt(fabs(Sx*Sy))) { + *area = 0.0; + *amaj = 0.0; + *amin = 0.0; + *angle = 0.0; + /* this is a poor fit - not an ellipse but a hyperbola */ + return (0); + } + + A1 = 1/(2*Sx); + A2 = 1/(2*Sy); + A3 = Sxy; + + *angle = atan2(-A3, A2 - A1) / 2.0; + R = sqrt( SQ(A2 - A1) + SQ(A3)); + root1 = (A1 + A2 + R); + root2 = (A1 + A2 - R); + + *area = 2.0*M_PI/sqrt(root1*root2); + *amaj = 2.35482*sqrt(1.0/root2); + *amin = 2.35482*sqrt(1.0/root1); + + return (1); + +} + + +/* In this function, Sx & Sy represent sx^2, sy^2 + + given an elliptical Gaussian of the form: + + exp (-z); z = x^2 / 2 sx^2 + y^2 / 2 sy^2 + xy Sxy + + Sa, Sb, and angle can be found by: + + A1 = 1/(2 Sx) + A2 = 1/(2 Sy) + A3 = Sxy + + R = sqrt ((A2 - A1)^2 + A3^2) + root1 = (A1 + A2 + R) + root2 = (A1 + A2 - R) + + angle = atan2 (-A3, A2 - A1) / 2.0 + Sa = sqrt (1.0/root2) + Sb = sqrt (1.0/root1) + +*/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/feature_fluxes.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/feature_fluxes.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/feature_fluxes.c (revision 2556) @@ -0,0 +1,56 @@ +# include "gophot.h" + +float feature_fluxes () { + + int i, j, npix, n; + float theta, phi, cs, sn; + float xp, yp, x, y, dx, dy, Dx, Dy; + float r1, r2; + float flux, mean, sky; + + for (n = 0; n < nregion; n++) { + + npix = flux = 0; + + dx = MAX (region[n][4], region[n][5]) + 2; + dy = MAX (region[n][4], region[n][5]) + 2; + cs = cos (region[n][6] * RAD_DEG); + sn = sin (region[n][6] * RAD_DEG); + sky = region[n][0]; + + for (j = region[n][3] - dy; j < region[n][3] + dy; j++) { + if (j < 0) continue; + if (j >= nslow) continue; + for (i = region[n][2] - dx; i < region[n][2] + dx; i++) { + if (i < 0) continue; + if (i >= nfast) continue; + + Dx = i - region[n][2]; + Dy = j - region[n][3]; + phi = atan2 (Dy, Dx) - RAD_DEG * region[n][6]; + theta = atan2 (region[n][4]*sin(phi), region[n][5]*cos(phi)); + + /* this is the point on the ellipse at the same angle as ref point */ + xp = region[n][4] * cos (theta); + yp = region[n][5] * sin (theta); + + x = xp * cs - yp * sn; + y = xp * sn + yp * cs; + + r1 = hypot (Dx, Dy); + r2 = hypot (x, y); + + if (r1 < r2) { + flux += big[i + nfast*j] - sky; + npix ++; + } + + } + } + /* existing image has mean subtracted = flux*Npix */ + flux += region[n][1]; + mean = flux / npix; + fprintf (stderr, "sub flux: %f, mag: %f, mean: %f\n", flux, -2.5*log10(flux), flux / npix); + region[n][7] = flux; + } +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/fillerup.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/fillerup.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/fillerup.c (revision 2556) @@ -0,0 +1,169 @@ +# include "gophot.h" + +fillerup (int xo, int yo, int findstats) { + + float *bigval, *noiseval, ufactor2, snlim2; + int i, j, ilo, ihi, jlo, jhi, ix, iy, idist, nsky, maxnpt, outer; + int *x, *y; + float *z, *d; + float N, S, F, Sx, Sy, Sx2, Sy2, xv, yv; + float value, SNR, sky; + + ufactor2 = SQ(ufactor); + snlim2 = SQ(snlim); + npt = 0; + nsky = 0; + + ilo = MAX (xo - nrect[1]/2, 0); + ihi = MIN (xo + nrect[1]/2, nfast-1); + jlo = MAX (yo - nrect[2]/2, 0); + jhi = MIN (yo + nrect[2]/2, nslow-1); + + outer = MIN (nrect[1]/3, nrect[2]/3); + + for (j = jlo; j <= jhi; j++) { + bigval = &big[ilo + j*nfast]; + noiseval = &noise[ilo + j*nfast]; + for (i = ilo; i <= ihi; i++, bigval++, noiseval++) { + if (!finite (*bigval)) continue; + if (!finite (*noiseval)) continue; + if (SQ(*bigval) < (*bigval + *noiseval) * snlim2) continue; + if (*bigval < (*noiseval - rnoise)) continue; + ix = i - xo; + iy = j - yo; + xs[npt] = ix; + ys[npt] = iy; + zs[npt] = *bigval/ufactor; + dzs[npt] = (*bigval + *noiseval) / ufactor2; + idist = MAX (abs(ix), abs(iy)); + if (idist > outer) { + ts[nsky] = zs[npt]; + nsky ++; + } + npt ++; + } + } + + if (npt == 0) return (FALSE); + if (!findstats) return (TRUE); + + if (nsky == 0) return (FALSE); + + sort (ts, nsky); + sky = ts[(int)(0.5*nsky)]; + sum2 = sky; + + z = zs; + d = dzs; + S = N = 0; + for (i = 0; i < npt; i++, z++, d++) { + if (*z < sky) continue; + N += *d; + S += *z - sky; + } + + if (S < 0) return (FALSE); + SNR = S*S / N; + if (SNR < SQ(bumpcrit)) return (FALSE); + + if (SNR > 49) { + x = xs; + y = ys; + z = zs; + d = dzs; + N = S = Sx = Sy = Sx2 = Sy2 = 0; + for (i = 0; i < npt; i++, x++, y++, z++, d++) { + idist = MAX (abs(*x), abs(*y)); + if (idist > outer) continue; + if (*z < sky) continue; + value = *z - sky; + xv = *x * value; + yv = *y * value; + N += *d; + S += value; + Sx += xv; + Sx2 += *x * xv; + Sy += yv; + Sy2 += *y * yv; + } + xmax = Sx / S; + ymax = Sy / S; + xmax2 = fabs(Sx2 / S - xmax*xmax); + ymax2 = fabs(Sy2 / S - ymax*ymax); + maxval = S / sqrt (xmax2*ymax2); + } else { + z = zs; + maxval = *z; + maxnpt = 0; + for (i = 0; i < npt; i++, z++) { + if (*z > maxval) { + maxnpt = i; + maxval = *z; + } + } + xmax = xs[maxnpt]; + ymax = ys[maxnpt]; + xmax2 = ava[4]; + ymax2 = ava[6]; + maxval -= sum2; + } + + maxval *= ufactor; + sum2 *= ufactor; + + if (!finite (maxval)) return (FALSE); + return (TRUE); + +} + +/* adjust the errors by the deviation from the guess fit */ +filladjust (float *star) { + + int i; + float f, df; + + return (0); + if (star[1] < cmax/(ufactor*10)) return (0); + + for (i = 0; i < npt; i++) { + + f = onestar ((int)xs[i], (int)ys[i], star, (float *) NULL); + + df = 10 * fabs(f - zs[i]); + dzs[i] += df/ufactor; + } + +} + +/* check if weighted center is in OK location + weighted center is based on useable pixels, not + on entire box */ +centertest (float *star, int *x, int *y) { + + int i; + float f, df, Sf, Sx, Sy, dx, dy; + + Sf = Sx = Sy = 0; + + for (i = 0; i < npt; i++) { + + f = onestar ((int)xs[i], (int)ys[i], star, (float *) NULL); + + Sf += f; + Sx += f*xs[i]; + Sy += f*ys[i]; + + } + + dx = star[2] - Sx/Sf; + dy = star[3] - Sy/Sf; + + if ((fabs(dx) > 0.25*nrect[1]) || (fabs(dy) > 0.25*nrect[2])) { + mprint (3, "center %f %f %f\n", dx, dy, Sf); + return (TRUE); + } + + return (FALSE); + +} + Index: /tags/gophot-1-0/Ohana/src/gophot/src/findsky.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/findsky.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/findsky.c (revision 2556) @@ -0,0 +1,36 @@ +# include "gophot.h" + +float findsky () { + + int i, skip, nsky, bin; + int zvalue[0x10000], nvalue; + float median, *valB; + + skip = 100; + valB = big; + + for (i = 0; i < nfast*nslow; i+=skip, valB+=skip) { + bin = *valB; + if ((bin >= 0) && (bin < 0x10000)) { + nvalue ++; + zvalue[bin] ++; + } + } + + nsky = 0; + for (i = 0; (i < 0x10000) && (nsky < nvalue / 2); i++) { + nsky += zvalue[i]; + } + if (nsky < nvalue / 2) { + mprint (1, "weird situation: no median found\n"); + exit (1); + } + median = i; + mprint (1, "median = %f\n", median); + + return (median); +} + +/* arguments OK */ + +/* this median method assumes 64k data range */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/galaxy.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/galaxy.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/galaxy.c (revision 2556) @@ -0,0 +1,58 @@ +# include "gophot.h" + +/* this is kind of stupid: this function takes two + lists of the parameters for the object because + star1 has x, y, relative to 0, */ + +bool galaxy (float *star1, float *err, float *star2) { + + float dummy[NPMAX], tot[4], temp, nsigma; + int i; + bool value; + + value = FALSE; + + for (i = 0; i < 4; i++) chi[i] = 0; + + parinterp (star2[2], star2[3], dummy); + + /* if errors on sigma x or y are so large, call it a galaxy */ + if (dummy[4] < 3*sqrt(err[4])) return (TRUE); + if (dummy[6] < 3*sqrt(err[6])) return (TRUE); + + /* too weak a measurement to be called galaxy */ + if (star1[4] < 2*sqrt(err[4])) return (FALSE); + if (star1[6] < 2*sqrt(err[6])) return (FALSE); + if (star1[1] < 2*sqrt(err[1])) return (FALSE); + + /************************************************************* + I'm rather concerned about 'parsm' where did it come from? + parms[] is set in varipar_plane.c */ + temp = SQ (sig[1]*dummy[4]); + tot[1] = MAX (parms[4], temp); + temp = SQ (sig[2])/(dummy[4]*dummy[6]); + tot[2] = MAX (parms[5], temp); + temp = SQ (sig[3]*dummy[6]); + tot[3] = MAX (parms[6], temp); + chi[1] = SQ (star1[4] - dummy[4]) / (tot[1] + err[4]); + chi[2] = SQ (star1[5] - dummy[5]) / (tot[2] + err[5]); + chi[3] = SQ (star1[6] - dummy[6]) / (tot[3] + err[6]); + + if (star1[4] < dummy[4]) chi[1] = 0; + if (star1[6] < dummy[6]) chi[3] = 0; + chi[4] = chi[1] + chi[2] + chi[3]; + + mprint (2, "galaxy test, object at %f %f %f\n", star2[2], star2[3], star2[1]); + mprint (2, "chisqs: %f %f %f %f\n", chi[1], chi[2], chi[3], chi[4]); + + nsigma = MIN (sqrt(chi[4]), 1.0e8); + mprint (2, "nsigma = %f\n", nsigma); + + value = chi[4] >= chicrit; + + return (value); + +} + +/* this function uses C 0,N-1 for a[], fa[] */ + Index: /tags/gophot-1-0/Ohana/src/gophot/src/gophot.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/gophot.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/gophot.c (revision 2556) @@ -0,0 +1,54 @@ +# include "gophot.h" + +main (int argc, char **argv) { + + float tmp; + + if (argc < 4) { + fprintf (stderr, "ERROR: usage: dophot imagename outfile paramfile\n"); + exit (0); + } + ConfigInit (&argc, argv); + + fprintf (stderr, "reading from %s, writing to %s, %s param file\n", argv[1], argv[2], argv[3]); + + FT_UNSIGN_MODE = FALSE; + strcpy (files[4], argv[2]); + + /* load image header and data */ + if (!fits_read_header (argv[1], &header)) { + fprintf (stderr, "ERROR: can't open FITS header %s\n", argv[1]); + exit (1); + } + if (!fits_read_matrix (argv[1], &matrix)) { + fprintf (stderr, "ERROR: can't open FITS matrix %s\n", argv[1]); + exit (1); + } + /* convert to float, set up noise array and axes */ + fits_convert_format (&header, &matrix, -32, 1.0, 0.0, 0); + ALLOCATE (noise, float, matrix.size); + big = (float *) matrix.buffer; + nfast = matrix.Naxis[0]; + nslow = matrix.Naxis[1]; + + /* override config values with header values */ + if (fits_scan (&header, "SATVALUE", "%f", 1, &tmp)) itop = tmp; + if (fits_scan (&header, "GAIN", "%f", 1, &tmp)) eperdn = tmp; + if (fits_scan (&header, "NEWGAIN", "%f", 1, &tmp)) eperdn = tmp; + + /* we are scaling the image to electrons, must also scale parameters */ + itop *= eperdn; + ibot *= eperdn; + tmax *= eperdn; + cmax *= eperdn; + ctpersat *= eperdn; + + dophot (); + + fprintf (stderr, "SUCCESS\n"); +} + +/* in the rest of the program: + NFAST = matrix.Naxis[0] -- # of X-axis pixels + NSLOW = matrix.Naxis[1] -- # of Y-axis pixels +*/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/guess.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/guess.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/guess.c (revision 2556) @@ -0,0 +1,72 @@ +# include "gophot.h" + +/* fill in values for star[NPMAX] based on rough stats and sky model */ +float guess1 (float *star, float *dummy, int ix, int iy) { + + /* parinterp fills in star[4,5,6] from averages */ + parinterp (ix, iy, star); + + star[0] = sum2/ufactor; + star[1] = maxval/ufactor; + star[2] = xmax; + star[3] = ymax; + + return (0.0); + +} + +/* fill in values for star[NPMAX] based on rough stats and sky model */ +/* use sigma_x, sigma_y from fillerup */ +float newguess (float *star, float *dummy, int ix, int iy) { + + star[0] = sum2 / ufactor; + star[1] = maxval / ufactor; + star[2] = xmax; + star[3] = ymax; + star[4] = xmax2; + star[5] = 0; + star[6] = ymax2; + + return (0.0); +} + +/* fill in values for star from instar, rescaling */ +float guess2 (float *star, float *instar, int *ix, int *iy) { + + /* i'm concerned that the ix, iy values need to be passed back. */ + float value; + + *ix = (int) (instar[2] + 0.5); + *iy = (int) (instar[3] + 0.5); + + star[0] = instar[0]/ufactor; + star[1] = instar[1]/ufactor; + star[2] = instar[2] - *ix; + star[3] = instar[3] - *iy; + star[4] = instar[4]; + star[5] = instar[5]; + star[6] = instar[6]; + + value = instar[0]; + return (value); +} + + +/* fill in values for star from instar, rescaling */ +float guess3 (float *star, float *instar, int *ix, int *iy) { + + float value; + + *ix = (int) (instar[2] + 0.5); + *iy = (int) (instar[3] + 0.5); + + parinterp (instar[2], instar[3], star); + + star[0] = instar[0]/ufactor; + star[1] = instar[1]/ufactor; + star[2] = instar[2] - *ix; + star[3] = instar[3] - *iy; + + value = instar[0]; + return (value); +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/impaper2.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/impaper2.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/impaper2.c (revision 2556) @@ -0,0 +1,42 @@ +# include "gophot.h" + +impaper2 (int k) { + + int i, nsky, nsort; + float sky, sum, r, dr, zsort[MAXFIL]; + + /* find median sky outside inner radius */ + + dr = 0.5*(SQ (0.5*arect[1]) + SQ (0.5*arect[2])); + nsort = 0; + for (i = 0; i < npt; i++) { + r = SQ (xs[i]) + SQ (ys[i]); + if (r >= dr) { + zsort[nsort] = zs[i]; + nsort ++; + } + } + + if (nsort < 1) return (0); + + sort (zsort, nsort); + + sky = 0; + nsky = 0; + for (i = 0.25*nsort; i < 0.75*nsort; i++) { + sky += zsort[i]; + nsky ++; + } + sky /= nsky; + + sum = 0; + for (i = 0; i < npt; i++) { + sum += (zs[i] - sky); + } + + apple[k][1] = sum*ufactor; + + return (1); + +} + Index: /tags/gophot-1-0/Ohana/src/gophot/src/improve.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/improve.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/improve.c (revision 2556) @@ -0,0 +1,117 @@ +# include "gophot.h" + +# define ADD +1 +# define SUB -1 + +improve (int last) { + + float err[NPMAX], star[NPMAX], sky, dx, dy, starchi; + int i, k, ix, iy, jmtype, niter, nfit; + bool skip, snok; + + /* imtype may be 1-9 for regular objects, 101-109 for fixed objects */ + + for (i = 0; i < nstot; i++) { + + /* set up flags for fixed objects */ + fixxy = fixpos && (imtype[i] >= 100); + jmtype = fixxy ? imtype[i] - 100 : imtype[i]; + /* inverse: imtype[i] = fixxy ? jmtype + 100 : jmtype; */ + + if (jmtype == 0) continue; + if (jmtype == 6) continue; + if (jmtype == 16) continue; + if (jmtype == 8) continue; + if (jmtype == 18) continue; + + /* if (jmtype == 3) continue; */ + if (jmtype == 2) continue; + if (jmtype == 12) continue; + if (jmtype == 10) continue; + + /* add star back to frame, get sky, get subraster */ + addstar (starpar[i], ADD, jmtype); + sky = guess3 (star, starpar[i], &ix, &iy); + mprint (3, " improving star %d at %d, %d\n", i, ix, iy); + nrect[1] = irect[1]; + nrect[2] = irect[2]; + fillerup (ix, iy, FALSE); /* we don't use the center information here */ + + /* skip star if off picture, has negative flux, or if s/n too low for non-fixed objects */ + skip = offpic (star, ix, iy, &dx, &dy); + if (star[1] <= 0) { + skip = TRUE; + starpar[i][1] = 0; + } + if (!skip && !fixxy) { + snok = transmask (ix,iy,sky); + skip = skip && !snok; + } + if (skip) { + jmtype = 6; + mprint (3, " deactivating star %d at %d, %d\n", i, ix, iy); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + /* fixed objects only vary x,y */ + if (fixxy) { + nfit = NFIT0; + niter = 2; + } else { + nfit = NFIT1; + niter = nit; + } + + filladjust (star); + + /* fit for magnitudes -- 4 param fit (or 2 for fixed obj) (initial guesses in are the previous fits) */ + starchi = chisq (onestar, xs, ys, zs, dzs, npt, star, err, nfit, acc, parlim, niter); + + /* set non-converge objects to type 4, don't bother with other calcs */ + if (!finite(starchi)) { + if (jmtype != 3) jmtype = 4; + for (k = 0; k < NPMAX; k++) galpar[i][k] = starpar[i][k]; + for (k = 0; k < NPMAX; k++) shadow[i][k] = starpar[i][k]; + apple[i][4] = MIN (1.086*sqrt(err[1])/star[1], 1); + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + /* things which get negative flux, turn in to type 6 and skip */ + if (star[1] <= 0) { + starpar[i][1] = 0; + jmtype = 6; + mprint (3, " deactivating star %d at %d, %d\n", i, ix, iy); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + /* this section applies to type 1 3 5 7 */ + if ((jmtype != 2) && (jmtype != 10)) { + /* store the fitted values */ + parupd (star, starpar[i], ix, iy); + } + + if ((jmtype != 2) && (jmtype != 3) && toofaint(starpar[i], err)) jmtype = 7; + + apple[i][4] = MIN (1, 1.086*sqrt(err[1])/star[1]); + + if (last) { + nrect[1] = arect[1]; + nrect[2] = arect[2]; + fillerup (ix, iy, FALSE); + impaper2 (i); + } + + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + + } + + return (0); + +} + +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/isearch.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/isearch.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/isearch.c (revision 2556) @@ -0,0 +1,170 @@ +# include "gophot.h" + +int isearch (int first) { + + int jrect[3]; + int i, j, nsprev, nfound, nnew, ngoodpts, ngoodmask, ngoodfit, nfit; + float dummy[2*NPMAX], star[2*NPMAX], err[NPMAX], pixt2, bestsky; + float tthresh, thresh2, highsky, highthresh, tmp; + float sky, Chi, dx, dy, rdnoise2; + float *Bval, *Nval; + bool iscosmic, isbright, isfaint, hole; + float get_mediansky (); + FILE *f; + + bzero (dummy, NPMAX*sizeof(float)); + + pixt2 = SQ(pixthresh); + + /* use jrect as temp in isearch.c */ + for (i = 1; i < 3; i++) { + jrect[i] = (first) ? krect[i] : irect[i]; + } + + nsprev = nstot; + nnew = 0; + ngoodpts = 0; + ngoodmask = 0; + ngoodfit = 0; + thresh2 = SQ(thresh); + nfit = (first) ? NFIT2 : NFIT1; + /* on first pass or so, we use 7 par fit, not 4 par fit */ + + Bval = big; + Nval = noise; + for (i = 0; i < nslow; i++) { + for (j = 0; j < nfast; j++, Bval++, Nval++) { + + /* occasionally update sky guess based on median sky image */ + if (!(j % nthpix)) { + sky = get_mediansky (j, i); + tthresh = sky + thresh; + } + + imtype[nstot] = 1; + + /* skip bad pixels */ + if (*Bval < tthresh) continue; + if (!finite (*Nval)) continue; + if (*Bval < *Nval - rnoise) continue; + + /* test for significant peak */ + if (!transmask (j, i, sky)) continue; + nnew ++; + + /* nrect is passed globally to fillerup or through - rather obscure */ + nrect[1] = jrect[1]; + nrect[2] = jrect[2]; + + /* fills the vectors xs, ys, zs, dzs */ + if (!fillerup (j, i, TRUE)) continue; + if ((fabs(xmax) > 0.33*nrect[1]) || (fabs(ymax) > 0.33*nrect[2])) { + mprint (3, "centroid moved: %d %d %f %f\n", j, i, xmax, ymax); + } + + tmp = enuff4*jrect[1]*jrect[2]; + if (npt < tmp) { + mprint (3, "skipping: npt = %d\n", npt); + continue; + } + ngoodpts ++; + + /* if sky guess has changed significantly, retest for significant peak */ + if (fabs(sum2 - sky) / maxval > 0.2) { + if (!transmask (j, i, sum2)) { + mprint (3, "failed transmask on local sky, %f\n", sum2); + continue; + } + } + ngoodmask ++; + + /* fill in guess for this object */ + if (first) { + newguess (star, dummy, j, i); + } else { + guess1 (star, dummy, j, i); + } + filladjust (star); + + /* fit star using 4 parameterfit */ + Chi = chisq (onestar, xs, ys, zs, dzs, npt, star, err, nfit, acc, parlim, nit); + + /* if (finite(Chi) && ((err[2] > 0.1) || (err[3] > 0.1))) fprintf (stderr, "errors: %d %d %f %f %f %f %f\n", j, i, err[2], err[3], err[4], err[6], Chi/npt); */ + if (!finite(Chi)) { + mprint (3, "failed to converge: no entry in starlist\n"); + continue; + } + if ((fabs(star[2]) > 0.33*nrect[1]) || (fabs(star[3]) > 0.33*nrect[2])) { + mprint (3, "fit center moved: %d %d %f %f\n", j, i, star[2], star[3]); + continue; + } + if (centertest (star, xs, ys, npt)) { + mprint (3, "star %d at %d %d moved\n", nstot, j, i); + continue; + } + if (offpic (star, j, i, &dx, &dy)) { + mprint (3, "fitted star off image\n"); + continue; + } + + /* save fit parameters */ + parupd (star, starpar[nstot], j, i); + parupd (star, shadow[nstot], j, i); + rchisq[nstot] = Chi / npt; + apple[nstot][4] = MIN (1, 1.086*sqrt(err[1])/star[1]); + mprint (2, "star %d: %f, %f peak: %f chisq: %f\n", nstot, starpar[nstot][2], starpar[nstot][3], maxval, rchisq[nstot]); + ngoodfit ++; + + if (verybright (starpar[nstot])) { + imtype[nstot] = 10; + goto finish; + } + if (cosmic (starpar[nstot])) { + imtype[nstot] = 8; + hole = oblit (starpar[nstot]); + goto finish_noadd; + } + if (toofaint (starpar[nstot], err)) { + imtype[nstot] = 7; + mprint (3, "faint imtype %d\n", imtype[nstot]); + goto finish; + } + + finish: + addstar (starpar[nstot], SUB, imtype[nstot]); + + finish_noadd: + nstot ++; + if (nstot >= NSMAX) { + mprint (0, "too many stars!\n"); + return (0); + } + + } + } + nfound = nstot - nsprev; + mprint (0, "stars found %d, stars tested %d\n", nfound, nnew); + mprint (0, "ngoodpts: %d, ngoodfit: %d, ngoodmask: %d\n", ngoodpts, ngoodfit, ngoodmask); + + return (nfound); +} + + +/* this function uses C 0,N-1 for a[], fa[] */ + + /* before star subtraction: + + Nval contains just sq(RN) (e) + Bval contains counts (e) + dBval = sqrt(Bval + Nval) + + after star subtraction: + + Nval contains fit + sq(RN) + Bval contains (obs-fit) + dBval = sqrt(Bval + Nval) [Bval' + Nval' = Bval + Nval] + */ + + /* note that the coordinates used in the fit are relative to j, i + and are adjusted back in parupd */ + Index: /tags/gophot-1-0/Ohana/src/gophot/src/large_features.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/large_features.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/large_features.c (revision 2556) @@ -0,0 +1,161 @@ +# include "gophot.h" + +/* these don't need to be 'static' unless we intend to call these + functions outside of this file... */ + +static int *x; +static int *y; +static int Npts = 0; +static int NPTS = 0; + +static int dx[] = {-1, 1, 0, 0}; +static int dy[] = {0, 0, -1, 1}; +static int Ntry = 4; + +static float *mediansky; +static int Nx, Ny; +static float fx, fy; + +float *copy_mediansky (int *, int *, float *, float *); +int set_value (int, int, float); +void large_features (float, float); +void get_neighbors (int, int, float); +float get_value (int, int); +int outline (float, float, float, float, float, float, float *); +float delete_ellipse (float *, float); + +void large_features (float sky, float dsky) { + + int i, j, k; + float min, flux, dSky; + float Xo, Yo, dX, dY; + float Xmin, Xmax, Ymin, Ymax; + float fitpars[5]; + + nregion = 0; + + /* use a copy to protect original */ + mediansky = copy_mediansky (&Nx, &Ny, &fx, &fy); + + dSky = sqrt(sky + rnoise); + dsky = MAX (dsky, dSky); + + min = sky + 3*dsky; + fprintf (stderr, "sky %f, %f, %f\n", sky, dsky, min); + /* find pixels which stand above threshold */ + + for (j = 0; j < Ny; j++) { + for (i = 0; i < Nx; i++) { + if (mediansky [i + j*Nx] > min) { + /* this is always the first point in the group */ + /* init the storage arrays */ + Npts = 0; + NPTS = 100; + ALLOCATE (x, int, NPTS); + ALLOCATE (y, int, NPTS); + x[Npts] = i; + y[Npts] = j; + Npts ++; + get_neighbors (i, j, min); + /* we now have a list x, y, Npts */ + Xmin = Xmax = x[0]; + Ymin = Ymax = y[0]; + for (k = 0; k < Npts; k++) { + Xmin = MIN (Xmin, x[k]); + Ymin = MIN (Ymin, y[k]); + Xmax = MAX (Xmax, x[k]); + Ymax = MAX (Ymax, y[k]); + } + Xo = 0.5*(Xmax + 1 + Xmin); + Yo = 0.5*(Ymax + 1 + Ymin); + dX = (Xmax + 1 - Xmin); + dY = (Ymax + 1 - Ymin); + convert_coords (&Xo, &Yo, &dX, &dY); + fprintf (stderr, "large feature: %f %f %f %f\n", Xo, Yo, dX, dY); + outline (Xo, Yo, 0.4*dX, 0.4*dY, sky + 7*dsky, 2*dsky, fitpars); + flux = delete_ellipse (fitpars, sky); + fprintf (stderr, " flux: %f %f\n", flux, sky); + region[nregion][0] = sky; + region[nregion][1] = flux; + region[nregion][2] = fitpars[0]; + region[nregion][3] = fitpars[1]; + region[nregion][4] = fitpars[2]; + region[nregion][5] = fitpars[3]; + region[nregion][6] = fitpars[4]; + region[nregion][7] = flux; + nregion ++; + if (nregion == 100) { + fprintf (stderr, "too many regions!\n"); + exit (0); + } + } + } + } +} + +void get_neighbors (int ix, int iy, float min) { + + int i, Ix, Iy; + + for (i = 0; i < Ntry; i++) { + + Ix = ix + dx[i]; + Iy = iy + dy[i]; + + if (get_value(Ix, Iy) > min) { + + x[Npts] = Ix; + y[Npts] = Iy; + Npts ++; + if (Npts == NPTS) { + NPTS += 100; + REALLOCATE (x, int, NPTS); + REALLOCATE (y, int, NPTS); + } + + set_value (Ix, Iy, 0); + + get_neighbors (Ix, Iy, min); + + } + } +} + +float get_value (int i, int j) { + + float value; + + if (i < 0) return (0); + if (i >= Nx) return (0); + + if (j < 0) return (0); + if (j >= Ny) return (0); + + value = mediansky[i + Nx*j]; + + return (value); + +} + +int set_value (int i, int j, float value) { + + if (i < 0) return (0); + if (i >= Nx) return (0); + + if (j < 0) return (0); + if (j >= Ny) return (0); + + mediansky[i + Nx*j] = value;; + + return (1); + +} + +convert_coords (float *X, float *Y, float *dX, float *dY) { + + *X /= fx; + *dX /= fx; + *Y /= fy; + *dY /= fy; + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/lubksb.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/lubksb.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/lubksb.c (revision 2556) @@ -0,0 +1,31 @@ +void lubksb (float **a, int n, int *indx, float *b) { + + int i, j, ii, kk; + float sum; + + ii = -1; + for (i = 0; i < n; i++) { + kk = indx[i]; + sum = b[kk]; + b[kk] = b[i]; + if (ii != -1) { + for (j = ii; j < i; j++) { + sum -= a[i][j]*b[j]; + } + } else { + if (sum != 0) ii = i; + } + b[i] = sum; + } + + for (i = n-1; i >= 0; i--) { + sum=b[i]; + if (i < n-1) { + for (j = i+1; j < n; j++) { + sum -= a[i][j]*b[j]; + } + } + b[i] = sum/a[i][i]; + } + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/ludcmp.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/ludcmp.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/ludcmp.c (revision 2556) @@ -0,0 +1,71 @@ +# include + +# define TINY 1.0e-20 + +void ludcmp (float **a, int n, int *indx, float *D) { + + float *vv; + int i, j, k, imax; + float aamax, sum, dum, d; + + ALLOCATE (vv, float, n); + + d = 1.0; + for (i = 0; i < n; i++) { + aamax = 0.0; + for (j = 0; j < n; j++) { + if (fabs (a[i][j]) > aamax) aamax = fabs (a[i][j]); + } + if (aamax == 0.0) { + *D = 0; + free (vv); + return; + } + vv[i] = 1.0/aamax; + } + + for (j = 0; j < n; j++) { + for (i = 0; i < j; i++) { + sum = a[i][j]; + for (k = 0; k < i; k++) { + sum -= a[i][k]*a[k][j]; + } + a[i][j] = sum; + } + aamax = 0.0; + for (i = j; i < n; i++) { + sum = a[i][j]; + for (k = 0; k < j; k++) { + sum -= a[i][k]*a[k][j]; + } + a[i][j] = sum; + dum = vv[i]*fabs(sum); + if (dum >= aamax) { + imax = i; + aamax = dum; + } + } + if (j != imax) { + for (k = 0; k < n; k++) { + dum = a[imax][k]; + a[imax][k] = a[j][k]; + a[j][k] = dum; + } + d = -d; + vv[imax] = vv[j]; + } + indx[j] = imax; + if (j != n-1) { + if (a[j][j] == 0.0) a[j][j] = TINY; + dum = 1.0 / a[j][j]; + for (i = j+1; i < n; i++) { + a[i][j] = a[i][j]*dum; + } + } + } + if (a[n-1][n-1] == 0.0) a[n-1][n-1] = TINY; + free (vv); + + *D = d; + return; +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/makemask.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/makemask.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/makemask.c (revision 2556) @@ -0,0 +1,26 @@ +# include "gophot.h" + +makemask () { + + int i, j, ii, jj; + float value; + float dumx, dumy, dummy, star[NPMAX]; + + dumx = 0; + dumy = 0; + dummy = parinterp (dumx, dumy, star); + + star[0] = 0; + star[1] = 1; + star[2] = 0; + star[3] = 0; + + for (j = 0, jj = -iyby2; j < 2*iyby2 + 1; j++, jj++) { + for (i = 0, ii = -ixby2; i < 2*ixby2 + 1; i++, ii++) { + starmask[i][j] = value = onestar (ii, jj, star, (float *) NULL); + mprint (2, "starmask[%d][%d]: %f\n", ii, jj, value); + } + } +} + +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/makenoise.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/makenoise.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/makenoise.c (revision 2556) @@ -0,0 +1,47 @@ +# include "gophot.h" + +makenoise () { + + int i, j; + float *valN, *valB; + + valN = noise; + valB = big; + + for (i = 0; i < nslow; i++) { + if ((i < nbadbot) || (i > nslow - nbadtop - 1)) { + for (j = 0; j < nfast; j++, valN++, valB++) *valN = MAGIC; + continue; + } + for (j = 0; j < nfast; j++, valN++, valB++) { + if (j < nbadleft) { + *valN = MAGIC; + continue; + } + if (j > nfast - nbadright - 1) { + *valN = MAGIC; + continue; + } + *valB *= eperdn; + if (*valB > itop) { + *valN = MAGIC; + continue; + } + if (*valB <= ibot) { + *valN = MAGIC; + continue; + } + *valN = rnoise; + } + } +} + +/* arguments OK */ + + +/* NEW CHOICES FOR NOISE AND BIG: + + 'big' is converted to electrons - this is the Poisson component of the error + 'noise' will contain just the SQ(read-noise) - this is added to when a star is subtracted + +*/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/mediansky.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/mediansky.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/mediansky.c (revision 2556) @@ -0,0 +1,197 @@ +# include "gophot.h" + +static float *mediansky; +static int Nx, Ny; +static float fx, fy; +sort (float *, int); + +int make_mediansky () { + + /* we have an image of nfast x nslow pix. + make a new image of size sqrt(nfast) x sqrt(nslow) + (2k x 4k -> 45 x 63) + */ + + float *temp; + int i, j, I0, I1, J0, J1, I, J, n; + int nx, ny; + float Mv, Nv, Mv2, value; + + Nx = sqrt (nfast); + Ny = sqrt (nslow); + + /* nthpix is used in isearch.c to update sky guess */ + nthpix = Nx; + + fx = (float) Nx / nfast; + fy = (float) Ny / nslow; + + ALLOCATE (mediansky, float, Nx*Ny); + + nx = 1 + 1/fx; + ny = 1 + 1/fy; + + ALLOCATE (temp, float, 2*nx*ny); + + Nv = Mv = Mv2 = 0.0; + + for (j = 0; j < Ny; j++) { + for (i = 0; i < Nx; i++) { + + I0 = i / fx; + J0 = j / fy; + + I1 = (i + 1) / fx; + J1 = (j + 1) / fy; + + n = 0; + temp[0] = 0; + + for (J = J0; J < J1; J++) { + if (J < nbadbot) continue; + if (J > nslow - nbadtop - 1) continue; + for (I = I0; I < I1; I++) { + if (I < nbadleft) continue; + if (I > nfast - nbadright - 1) continue; + temp[n] = big[J*nfast + I]; + n++; + } + } + + sort (temp, n); + value = temp[(int)(0.5*n)]; + + if (n < 2) { + mediansky[j*Nx + i] = MAGIC; + } else { + mediansky[j*Nx + i] = value; + } + } + } + return (0); +} + +/* take array which is the median sky, copy to a separate vector, + sort, find median and sigma in 50% interval around median */ + +get_skystats (float *sky, float *dsky) { + + float *temp, sky2; + int i, Npix, Nsky; + + Npix = Nx*Ny; + + ALLOCATE (temp, float, Npix); + + memcpy (temp, mediansky, Npix*sizeof(float)); + + sort (temp, Npix); + + *sky = temp[(int)(0.5*Npix)]; + + Nsky = 0; + sky2 = 0; + for (i = 0.25*Npix; i < 0.75*Npix; i++) { + if (!finite(temp[i])) continue; + sky2 += SQ(temp[i] - *sky); + Nsky ++; + } + + *dsky = sqrt (sky2/Nsky); + + free (temp); + + mprint (0, "median sky: %f %f\n", *sky, *dsky); + + return (1); + +} + +fix_mediansky (float sky) { + + int i; + + for (i = 0; i < Nx*Ny; i++) { + if (!finite (mediansky[i])) { + mediansky[i] = sky; + } + } +} + +/* this version is fast and has no edge problems. */ + +# if (1) +float get_mediansky (int i, int j) { + + float value; + int I, J; + + I = i * fx; + J = j * fy; + + value = mediansky[J*Nx + I]; + + return (value); + +} +# endif + +/* this version is slow, but more accurate. + on the other hand, it needs to be fixed for + edge problems. */ + +# if (0) +float get_mediansky (int i, int j) { + + float Fx, Fy, Vo, Vm, Vx, Vy; + float value; + int I, J, dx, dy; + + I = i * fx; + J = j * fy; + + Vo = mediansky[J*Nx + I]; + + Fx = i * fx - I; + if (Fx < 0.5) { + dx = -1; + } else { + dx = +1; + Fx = 1.0 - Fx; + } + Vm = mediansky[J*Nx + I + dx]; + Vx = (Vo*(0.5+Fx) + Vm*(0.5-Fx)); + + Fy = j * fy - J; + if (Fy < 0.5) { + dy = -1; + } else { + dy = +1; + Fy = 1.0 - Fy; + } + Vm = mediansky[J*Nx + I + dy*Nx]; + Vy = (Vo*(0.5+Fy) + Vm*(0.5-Fy)); + + value = 0.5*(Vx + Vy); + + return (value); + +} +# endif + +float *copy_mediansky (int *nx, int *ny, float *Fx, float *Fy) { + + float *temp; + + ALLOCATE (temp, float, Nx*Ny); + memcpy (temp, mediansky, Nx*Ny*sizeof(float)); + + *nx = Nx; + *ny = Ny; + *Fx = fx; + *Fy = fy; + + return (temp); + +} + Index: /tags/gophot-1-0/Ohana/src/gophot/src/message.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/message.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/message.c (revision 2556) @@ -0,0 +1,29 @@ +# include +# include + +static FILE *f; +static int verbosity = 0; + +int mprint (int level, char *mode, ...) { + + int status; + va_list argp; + + if (level > verbosity) return (0); + + va_start (argp, mode); + + status = vfprintf (f, mode, argp); + + va_end (argp); + + return (status); + +} + +set_verbosity (int level) { + + f = stderr; + verbosity = level; + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/objout.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/objout.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/objout.c (revision 2556) @@ -0,0 +1,88 @@ +# include "gophot.h" +# define NCHAR 104 + +completeout () { + + FILE *f; + int i, Nchar; + float gcorr, gmag; + float area, amajor, aminor, tilt, fmag, xc, yc, apmag, tmp; + char line[NCHAR]; + + f = fopen (files[4], "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "error, can't save data in file %s\n", files[4]); + exit (1); + } + + ellipse (ava[4], ava[5], ava[6], &area, &amajor, &aminor, &tilt); + fprintf (f, "# Average Star: %f %f %f\n", amajor, aminor, tilt); + + for (i = 0; i < nregion; i++) { + + fmag = -2.5*log10 (region[i][1]); + gmag = -2.5*log10 (region[i][7]); + apmag = 99.999; + + Nchar = snprintf (line, NCHAR, "%3d %8.2f %8.2f %8.3f %6.3f %9.2f %9.3f %9.3f %7.2f %8.3f %8.3f %8.2f", + 20, region[i][2], region[i][3], fmag, 0.01, region[i][0], region[i][4], region[i][5], region[i][6], + gmag, apmag, 10.0); + fprintf (f, "%s\n", line); + } + + for (i = 0; i < nstot; i++) { + + fmag = 99.999; + gmag = 99.999; + apmag = 99.999; + + if (imtype[i] != 8) { + /* pure gaussian fit mags */ + ellipse (starpar[i][4], starpar[i][5], starpar[i][6], &area, &amajor, &aminor, &tilt); + tmp = area*starpar[i][1]/eperdn; + if (tmp > 0.0) fmag = -2.5 * log10 (tmp); + /* galaxy-non-gauss fit mags */ + ellipse (shadow[i][4], shadow[i][5], shadow[i][6], &area, &amajor, &aminor, &tilt); + gcorr = 1.0; + tmp = area*gcorr*shadow[i][1]/eperdn; + if (tmp > 0.0) gmag = -2.5 * log10 (tmp); + tilt = 57.29578 * tilt; + } else { + /* get correct orientation for oblit boxes */ + if (starpar[i][5] != -1) fmag = -99.999; + if (starpar[i][4] >= starpar[i][6]) { + amajor = starpar[i][4]; + aminor = starpar[i][6]; + tilt = 0.0; + } else { + amajor = starpar[i][6]; + aminor = starpar[i][4]; + tilt = 90.0; + } + } + + /* what is the value of the center of a pixel? */ + /* this assumes the pixel center is at 0,0, not 0.5, 0.5 */ + xc = starpar[i][2]; + yc = starpar[i][3]; + + if (apple[i][1] > 0.0) apmag = -2.5 * log10 (apple[i][1]/eperdn); + + Nchar = snprintf (line, NCHAR, "%3d %8.2f %8.2f %8.3f %6.3f %9.2f %9.3f %9.3f %7.2f %8.3f %8.3f %8.2f", + imtype[i], xc, yc, fmag, apple[i][4], shadow[i][0], amajor, aminor, tilt, + gmag, apmag, rchisq[i]); + if (Nchar != NCHAR - 1) { + mprint (1, "funny line %d\n", i); + } + fprintf (f, "%s\n", line); + } + + /* + 1 37.51 1193.91 -15.246 .034 594.48 6.657 5.654 -7.06 -15.611 -15.822 + */ + fclose (f); + + return (0); +} + +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/oblims.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/oblims.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/oblims.c (revision 2556) @@ -0,0 +1,31 @@ +# include "gophot.h" + +oblims (float *star, int *jrect) { + + float temp, fudgex, fudgey; + + if (star[1] > 0) { + /* think about this line: */ + temp = pow ((6*star[1]/nphob), 0.33333333); + if (star[4] > 0) { + fudgex = sqrt(temp*star[4]*2); + } else { + fudgex = 10; + } + if (star[6] > 0) { + fudgey = sqrt(temp*star[6]*2); + } else { + fudgey = 10; + } + } else { + fudgex = 10; + fudgey = 10; + } + + mprint (3, "fudge size: %f %f\n", fudgex, fudgey); + jrect[1] = star[2] - fudgex; + jrect[2] = star[2] + fudgex; + jrect[3] = star[3] - fudgey; + jrect[4] = star[3] + fudgey; + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/oblit.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/oblit.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/oblit.c (revision 2556) @@ -0,0 +1,34 @@ +# include "gophot.h" + +bool oblit (float *star) { + + bool value; + + float wx, wy, dum, duma[NPMAX]; + int ix, iy, ixhi, ixlo, iyhi, iylo, jx, jy; + + dum = guess2 (duma, star, &ix, &iy); + wx = star[4]; + wy = star[6]; + + ixhi = MIN ((int)(ix + 0.5*wx + 0.5), nfast-1); + ixlo = MAX ((int)(ix - 0.5*wx + 0.5), 0); + iyhi = MIN ((int)(iy + 0.5*wy + 0.5), nslow-1); + iylo = MAX ((int)(iy - 0.5*wy + 0.5), 0); + + mprint (3, "obliterating following region : %d - %d, %d - %d\n", ixlo, ixhi, iylo, iyhi); + for (jy = iylo; jy <= iyhi; jy++) { + for (jx = ixlo; jx <= ixhi; jx++) { + big[jx+jy*nfast] = MAGIC; + noise[jx+jy*nfast] = MAGIC; + } + } + + mprint (2, "obliteration: %d, %d (%f x %f)\n", ix, iy, wx, wy); + + return (TRUE); + +} + + + Index: /tags/gophot-1-0/Ohana/src/gophot/src/offpic.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/offpic.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/offpic.c (revision 2556) @@ -0,0 +1,24 @@ +# include "gophot.h" + +bool offpic (float *star, int ix, int iy, float *dx, float *dy) { + + bool nogood; + float x, y; + + x = ix + star[2]; + y = iy + star[3]; + + *dx = 0; + if (x < 0) *dx = -x; + if (x > nfast) *dx = x - nfast; + + *dy = 0; + if (y < 0) *dy = -y; + if (y > nslow) *dy = y - nslow; + + nogood = (*dx != 0) || (*dy != 0); + + if (!fixpos) nogood = nogood || (star[1] < 0); + + return (nogood); +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/outline.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/outline.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/outline.c (revision 2556) @@ -0,0 +1,188 @@ +# include "gophot.h" + +int Npts; +float *xv, *yv, *zv; +float par[5]; +float dpar[5]; +float outline_chi (float); + +/* fit pars[5] to ellipse at Xo, Yo, dX, dY */ +int outline (float Xo, float Yo, float dX, float dY, float Io, float dIo, float *pars) { + + int i, j, k, Nx, Ny, NPTS, BigChange, ABigChange; + float oChi, dchi, Chi, Chi_p, Chi_m, dp, tmp_par, nChi, tmp1, tmp2; + float curve, frac, value; + + par[0] = Xo; + par[1] = Yo; + par[2] = dX; + par[3] = dY; + par[4] = 0.0; + + dpar[0] = 10; + dpar[1] = 10; + dpar[2] = 10; + dpar[3] = 10; + dpar[4] = 10; + + /* find all pixels within range Io-dIo : Io+dIo, in region about center guess */ + + Npts = 0; + NPTS = 1000; + ALLOCATE (xv, float, NPTS); + ALLOCATE (yv, float, NPTS); + ALLOCATE (zv, float, NPTS); + for (j = par[1]-2*par[3]; j < par[1]+2*par[3]; j++) { + if (j < 0) continue; + if (j >= nslow) continue; + for (i = par[0]-2*par[2]; i < par[0]+2*par[2]; i++) { + if (i < 0) continue; + if (i >= nfast) continue; + value = big[i + nfast*j]; + if (fabs (value - Io) < dIo) { + xv[Npts] = i; + yv[Npts] = j; + zv[Npts] = value; + Npts ++; + if (Npts == NPTS) { + NPTS += 1000; + REALLOCATE (xv, float, NPTS); + REALLOCATE (yv, float, NPTS); + REALLOCATE (zv, float, NPTS); + } + } + } + } + + if (Npts == 0) { + fprintf (stderr, "no valid points in box, try again\n"); + free (xv); + free (yv); + free (zv); + return (FALSE); + } + + Chi = outline_chi (Io); + for (j = 0; j < 15; j++) { + + oChi = Chi; + for (i = 0; i < 5; i++) { + /* find +Chi, -Chi for this par & adjust par as needed */ + + for (k = 0, BigChange = TRUE; (k < 3) && BigChange; k++) { + tmp_par = par[i]; + par[i] = tmp_par + dpar[i]; + Chi_p = outline_chi (Io); + par[i] = tmp_par - dpar[i]; + Chi_m = outline_chi (Io); + + /* have we braketted a minimum? (curve < 0) */ + curve = (Chi_p - Chi) * (Chi - Chi_m); + if (curve > 0) { + dp = 2*dpar[i]; + } else { + dp = 0.5 * dpar[i] * (Chi_m - Chi_p) / (Chi_m + Chi_p - 2*Chi); + } + if (Chi_m + Chi_p - 2*Chi == 0) dp = 0; + /* don't let extrapolation go too far */ + if (fabs (dp) > 2*fabs(dpar[i])) { dp = SIGN(dp) * fabs (2*dpar[i]); } + + par[i] = tmp_par + dp; + Chi = outline_chi (Io); + + BigChange = FALSE; + if (Chi <= 1.001*oChi) { + /* got better */ + dchi = (oChi - Chi) / oChi; + if ((dchi > 0.03) || (curve > 0)) BigChange = TRUE; + } else { + par[i] = tmp_par; + Chi = oChi; + if (Chi_m < Chi) { + Chi = Chi_m; + par[i] = tmp_par - dpar[i]; + } + if (Chi_p < Chi) { + Chi = Chi_p; + par[i] = tmp_par + dpar[i]; + } + } + oChi = Chi; + } + if (!BigChange) dpar[i] *= 0.8; + } + + mprint (0, "try: %d %f ", j, Chi); + for (i = 0; i < 5; i++) { + mprint (0, "%f ", par[i]); + } + mprint (0, "\n"); + for (i = 0; i < 5; i++) { + mprint (2, "%f ", dpar[i]); + } + mprint (2, "\n"); + dchi -= Chi; + + } + + free (xv); + free (yv); + free (zv); + + pars[0] = par[0]; + pars[1] = par[1]; + pars[2] = par[2]; + pars[3] = par[3]; + pars[4] = par[4]; + return (TRUE); + +} + +/* par[0] = x + par[1] = y + par[2] = dx + par[3] = dy + par[4] = dxy + + xp = par[2] * cos (t); + yp = par[3] * sin (t); + + x = xp * cos (par[4] * RAD_DEG) - yp * sin (par[4] * RAD_DEG) + par[0]; + y = xp * sin (par[4] * RAD_DEG) + yp * cos (par[4] * RAD_DEG) + par[1]; + +*/ + +float outline_chi (float Io) { + + int i; + float theta, phi; + float xp, yp, x, y; + float Chi, dv, Dv, R2, F2; + float cs, sn; + + Chi = 0; + + cs = cos (par[4] * RAD_DEG); + sn = sin (par[4] * RAD_DEG); + + for (i = 0; i < Npts; i++) { + + phi = atan2 (yv[i] - par[1], xv[i] - par[0]) - RAD_DEG * par[4]; + theta = atan2 (par[2]*sin(phi), par[3]*cos(phi)); + + /* this is the point on the ellipse at the same angle as ref point */ + xp = par[2] * cos (theta); + yp = par[3] * sin (theta); + + x = xp * cs - yp * sn + par[0]; + y = xp * sn + yp * cs + par[1]; + + R2 = sqrt (SQ (x - xv[i]) + SQ (y - yv[i])); + Chi += R2; + + } + + Chi = Chi / Npts; + return (Chi); + +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/paravg.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/paravg.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/paravg.c (revision 2556) @@ -0,0 +1,98 @@ +# include "gophot.h" + +paravg () { + + float sum[NPMAX][2]; + int i, j, n; + FILE *f; + char c; + float *ave4, *ave5, *ave6, *key; + int Ngood, NGOOD, *index, *idx; + int Nave; + float a4, a5, a6; + + for (j = 0; j < NPAR; j++) { + sum[j][0] = 0; + sum[j][1] = 0; + } + + if (nstot < 1) return (0); + + NGOOD = 100; + Ngood = 0; + ALLOCATE (index, int, NGOOD); + + /* find good stars, store index */ + for (i = 0; i < nstot; i++) { + if (starpar[i][1] < 100) continue; + if ((imtype[i] == 1) || (imtype[i] == 101)) { + index[Ngood] = i; + Ngood ++; + if (Ngood == NGOOD) { + NGOOD += 100; + REALLOCATE (index, int, NGOOD); + } + } + } + for (i = 0; i < Ngood; i++) { + fprintf (stderr, "S: %d %d %d %f %f %f %f\n", i, index[i], imtype[index[i]], starpar[index[i]][1], shaderr[index[i]][1], starpar[index[i]][4], shadow[index[i]][4]); + } + + /* find good stars, store index */ + if (Ngood < 10) { /* accept type 2 as well... */ + for (i = 0; i < nstot; i++) { + if ((imtype[i] == 2) || (imtype[i] == 102)) { + index[Ngood] = i; + Ngood ++; + if (Ngood == NGOOD) { + NGOOD += 100; + REALLOCATE (index, int, NGOOD); + } + } + } + } + + if (Ngood > 10) { /* don't change if not enough 'stars' */ + + ALLOCATE (ave4, float, Ngood); + ALLOCATE (ave5, float, Ngood); + ALLOCATE (ave6, float, Ngood); + ALLOCATE (key, float, Ngood); + ALLOCATE (idx, int, Ngood); + + for (i = 0; i < Ngood; i++) { + idx[i] = i; + key[i] = shadow[index[i]][4]; + ave4[i] = shadow[index[i]][4]; + ave5[i] = shadow[index[i]][5]; + ave6[i] = shadow[index[i]][6]; + } + sortkey (idx, key, Ngood); + + a4 = a5 = a6 = Nave = 0; + for (i = 0.4*Ngood; i < 0.6*Ngood; i++) { + a4 += ave4[idx[i]]; + a5 += ave5[idx[i]]; + a6 += ave6[idx[i]]; + Nave ++; + } + ava[4] = a4/Nave; + ava[5] = a5/Nave; + ava[6] = a6/Nave; + + free (key); + free (idx); + free (ave4); + free (ave5); + free (ave6); + } + free (index); + + + mprint (1, "%d stars used to find parameter averages\n", Ngood); + mprint (1, "average values so far of shape parameters for stars: \n"); + for (j = 4; j < NPAR; j++) mprint (1, "%f ", ava[j]); + mprint (1, "\n"); + +} +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/parinterp.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/parinterp.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/parinterp.c (revision 2556) @@ -0,0 +1,10 @@ +# include "gophot.h" + +float parinterp (float x, float y, float *star) { + + star[4] = ava[4]; + star[5] = ava[5]; + star[6] = ava[6]; + + return (star[0]); +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/parupd.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/parupd.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/parupd.c (revision 2556) @@ -0,0 +1,25 @@ +# include "gophot.h" + +parupd (float *instar, float *outstar, int ix, int iy) { + + outstar[0] = instar[0]*ufactor; + outstar[1] = instar[1]*ufactor; + outstar[2] = instar[2] + ix; + outstar[3] = instar[3] + iy; + outstar[4] = instar[4]; + outstar[5] = instar[5]; + outstar[6] = instar[6]; +} + +/* a legacy function +twoupd (float *instar, float *outstar, int ix, int iy) { + + outstar[0] = instar[0]; + outstar[1] = instar[1]; + outstar[2] = instar[2] + ix; + outstar[3] = instar[3] + iy; + outstar[4] = instar[4]; + outstar[5] = instar[5] + ix; + outstar[6] = instar[6] + iy; +} +*/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/pseud2d.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/pseud2d.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/pseud2d.c (revision 2556) @@ -0,0 +1,84 @@ +# include "gophot.h" + +# define HALF 0.5 +# define THIRD 0.333 +# define EXPMIN -15.0 + +float pseud2d (int ix, int iy, float *a, float *fa) { + + double x, y, t5, t6, t7, T; + double denom, dt, pexp, R; + double Q; + float value; + int i; + + x = ix - a[2]; + y = iy - a[3]; + + t6 = a[5]*y; + t5 = x/a[4]; + t7 = y/a[6]; + + T = 0.5*((t5 + 2.0*t6)*x + t7*y); + + if (!finite (a[1])) fprintf (stderr, "bad star!\n"); + + if (T >= 0) { + denom = 1.0 + T*(beta4 + HALF*beta64*T*(1.0 + THIRD*T)); + pexp = 1.0 / denom; + dt = beta4 + beta64*T*(1.0 + HALF*T); + if (!finite (pexp)) { + fprintf (stderr, "error in pseud2d: %f %f %f %d %d\n", T, denom, pexp, ix, iy); + denom = 1.0 + T + HALF*beta4*T*T + THIRD*beta64*T*T*T; + pexp = 1.0 / denom; + if (!finite (pexp)) { + fprintf (stderr, "error in pseud2d: %f %f %f %d %d\n", T, denom, pexp, ix, iy); + fprintf (stderr, "%f %f %f %f %f\n", x, y, t5, t6, t7); + fprintf (stderr, "%f %f %f\n", T, T*HALF*beta4*T, T*HALF*beta4*T*THIRD*beta64*T); + for (i = 0; i < 7; i++) fprintf (stderr, "%d %f\n", i, a[i]); + fprintf (stderr, "nstot: %d\n", nstot); + fprintf (stderr, "invalid value!\n"); + exit (1); + } + } + } else { + denom = MAX (T, EXPMIN); + denom = fabs (denom); + pexp = exp (denom); + denom = 1.0; + dt = 1.0; + } + + value = a[1]*pexp; + + if (fa != (float *) NULL) { + R = value*dt/denom; + fa[1] = pexp; + fa[2] = R*(t5 + t6); + fa[3] = R*(a[5]*x + t7); + fa[4] = R*t5*t5*0.5; + fa[6] = R*t7*t7*0.5; + + fa[5] = -R*x*y; + fa[0] = 1.0; + } + + value += + a[0]; + + if (!finite(value)) { + fprintf (stderr, "error in pseud2d: %f %f %f %f\n", T, denom, pexp, a[1]); + exit (1); + } + + return (value); + +} + +/* this function uses C 0,N-1 for a[], fa[] */ + +/* a[0] - sky + a[1] - Io + a[2], a[3] - X, Y + a[4], a[5], a[6] - shape +*/ + Index: /tags/gophot-1-0/Ohana/src/gophot/src/pseud4d.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/pseud4d.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/pseud4d.c (revision 2556) @@ -0,0 +1,73 @@ +# include "gophot.h" + +# define HALF 0.5 +# define THIRD 0.33333333333333 +# define EXPMIN -15 + +float pseud4d (int ix, int iy, float *a, float *fa) { + + int i, ioff; + float x, y; + float denom, T, dt, pexp[2]; + float a2, a5, a6, a7, t5, t6, t7, value; + + a5 = 1.0/a[7]; + a7 = 1.0/a[9]; + a6 = a[8]; + + for (i = 0; i < 2; i++) { + + ioff = 3*i; + x = ix - a[2 + ioff]; + y = iy - a[3 + ioff]; + + t5 = a5*x; + t6 = a6*y; + t7 = a7*y; + T = HALF*((t5 + 2*t6)*x + t7*y); + + if (T > 0) { + denom = 1.0 + T*(1.0 + HALF*beta4*T*(1.0 + THIRD*beta64*T)); + pexp[i] = 1./denom; + } else { + denom = MAX (T, EXPMIN); + denom = fabs (denom); + pexp[i] = exp (denom); + denom = 1.0; + } + + a2 = exp(a[1 + ioff]); + value = a2*pexp[i]; + pexp[i] = value; + + if (fa != (float *) NULL) { + if (T > 0) { + dt = 1.0 + beta4*T*(1.0 + HALF*beta64*T); + } else { + dt = 1.0; + } + fa[1 + ioff] = value; + T = pexp[i]*dt/denom; + fa[2 + ioff] = (a5*x + a6*y)*T; + fa[3 + ioff] = (a6*x + a7*y)*T; + + fa[0] = 1; + } + } + + value = pexp[0] + pexp[1] + a[0]; + + return (value); + +} + +/* this function uses C 0,N-1 for a[], fa[] */ + +/* a[0] - sky + a[1] - I1 + a[2], a[3] - X1, Y1 + a[4] - I2 + a[5], a[6] - X2, Y2 + a[7], a[8], a[9] - shape +*/ + Index: /tags/gophot-1-0/Ohana/src/gophot/src/set_thresholds.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/set_thresholds.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/set_thresholds.c (revision 2556) @@ -0,0 +1,33 @@ +# include "gophot.h" + +float set_thresholds (float sky, float dSky, float *factor, int *Nit) { + + float dsky, tmp; + + /* tmin starts with Nsigma above sky */ + + dsky = sqrt(sky + rnoise); + fprintf (stderr, "dsky: %f %f %f %f\n", rnoise, sky, dsky, dSky); + dsky = MAX (dsky, dSky); + tmin = tmin * dsky; + + mprint (0, "gain is %f, median sky is %f, min threshold is %f\n", + eperdn, sky, tmin); + + /* set starting threshold */ + tmp = pow (2.0, tfac); + tmax = 0.5*itop; + *Nit = log (tmax / tmin) / log (tmp); + *factor = pow ((tmax / tmin), (1.0 / *Nit)); + + return (tmax); + +} + +/* we are going to evenly divide the range from + sky + tmin*dsky to 0.5*itop in slices with ratios + close to tfac */ + +/* images with very messy backgrounds should not be pushed to the + absolute minimum. We let the sky sigma be the max of the + formal and the measured sigmas */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/shape.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/shape.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/shape.c (revision 2556) @@ -0,0 +1,267 @@ +# include "gophot.h" + +# define ZERO_SHADOW {shadow[i][0] = 0; \ + shadow[i][1] = 0; \ + shadow[i][4] = 1; \ + shadow[i][5] = 0; \ + shadow[i][6] = 1; } + +/* try 7 par fit for each object, decide on object types */ +shape () { + + bool notnuff, gotfaint, verybig, offp, converge; + + /* watch twofpar - i think we only need 2 entries, for storage in this loop */ + float star[NPMAX], err[NPMAX], star1[NPMAX], star2[NPMAX]; + float sky, galchi, dx, dy, starchi; + int i, j, n, k, ix, iy, jmtype, nsprev; + + test7 = TRUE; + nsprev = nstot; + + for (i = 0; i < nsprev; i++) { + mprint (3, " determining shape for object no. %d\n", i); + + gotfaint = FALSE; + fixxy = imtype[i] >= 100; + jmtype = (fixxy) ? imtype[i] - 100 : imtype[i]; + + switch (jmtype) { + case 7: + mprint (3, "too faint, skipping object %d at %f, %f\n", i, starpar[i][2], starpar[i][3]); + continue; + break; + case 4: + case 9: + mprint (3, "nonconverge, skipping object %d at %f, %f\n", i, starpar[i][2], starpar[i][3]); + continue; + break; + case 8: + case 18: + case 6: + case 16: + mprint (3, "poor obj, skipping object %d at %f, %f\n", i, starpar[i][2], starpar[i][3]); + ZERO_SHADOW; + continue; + break; + } + + if (shadow[i][0] == 0) { + sky = guess2 (star, starpar[i], &ix, &iy); + } else { + sky = guess2 (star, shadow[i], &ix, &iy); + } + + addstar (starpar[i], ADD, jmtype); + + nrect[1] = irect[1]; + nrect[2] = irect[2]; + fillerup (ix, iy, FALSE); + notnuff = (npt < enuff7*irect[1]*irect[2]); + + if (notnuff) { + mprint (3, "obj %d, npts %d, %d & %d - skipping star: not enough pixels for 7-parm fit \n", i, npt, ix, iy); + if (jmtype != 2) jmtype = 5; + if (jmtype == 12) beta4 = 0.01; + addstar (starpar[i], SUB, jmtype); + beta4 = 1.0; + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + if (fixxy) { + gotfaint = !transmask (ix, iy, sky); + if (gotfaint) { + if (jmtype != 3) jmtype = 7; + mprint (3, "obj %d, npts %d, %d & %d - skipping star: too faint for 7 par fit\n", i, npt, ix, iy); + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + } + + filladjust (star); + mprint (2, "obj %d, npts %d, %d & %d\n", i, npt, ix, iy); + + /* fit 'extended' to star (7 parameter fit) */ + galchi = chisq (onestar, xs, ys, zs, dzs, npt, star, err, NFIT2, acc, parlim, nit); + /* + if (finite(galchi) && ((err[4] > 0.25) || (err[6] > 0.25))) { + fprintf (stderr, "7 errors: %d %d %f %f %f %f %f\n", ix, iy, err[2], err[3], err[4], err[6], galchi/npt); + fprintf (stderr, " %f %f %f %f\n", star[2], star[3], star[4], star[6]); + } + */ + if (fabs(star[5]) >= 1.0 / sqrt(fabs(star[4]*star[6]))) { + mprint (3, "hyperbola! problem with %d at %f %f\n", i, starpar[i][2], starpar[i][3]); + /* just a warning... */ + } + if ((star[4] > SQ(0.25*nrect[1])) || (star[6] > SQ(0.25*nrect[1]))) { + mprint (3, "warning - dubious fit for %d at %f %f\n", i, starpar[i][2], starpar[i][3]); + /* just a warning... */ + } + if (!finite (galchi)) { + /* shadow, starpar keep 4-par fit */ + mprint (3, "non-converge object %d at %f %f\n", i, starpar[i][2], starpar[i][3]); + if (jmtype != 3) jmtype = 9; + mprint (3, "obj %d, npts %d, %d & %d .... failed to converge\n", i, npt, ix, iy); + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + /* is new fit position reasonable? */ + if (centertest (star, xs, ys, npt) || (fabs(star[2]) > 0.33*nrect[1]) || (fabs(star[3]) > 0.33*nrect[2])) { + mprint (3, "fit center moved 2: %f %f %f %f\n", starpar[i][2], starpar[i][3], star[2], star[3]); + jmtype = 16; + /* shadow, starpar keep 4-par fit */ + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + /* calculate error here for type 2, 10. redone in improve for 1, 3, 5, 7 */ + apple[i][4] = MIN (1, 1.086*sqrt(err[1])/star[1]); + rchisq[i] = galchi / npt; + + /* save 7-par fit in shadow */ + parupd (star, shadow[i], ix, iy); + for (j = 0; j < NPAR; j++) shaderr[i][j] = err[j]; + + /* allow type 10 to go back to 1 and vice-versa */ + if ((jmtype == 10) && !verybright (starpar[i])) jmtype = 1; + if (verybright (starpar[i])) jmtype = 10; + + if (jmtype == 10) { + mprint (3, " very bright, keep 7 par fit only\n"); + + if ((thresh < tmax / 4) && (thresh > tmax/16)) { + /* test if type 10 can be safely split in two */ + starchi = twofit (star, starpar[i], star1, star2); + /* has to be a really significant improvement */ + if (finite(starchi) && (starchi/galchi < 10*stograt) && (starchi < 10*npt)) { + mprint (3, " result-> a split star: gal-chi: %f, star-chi: %f\n", galchi, starchi); + jmtype = 3; + parupd (star1, shadow[i], ix, iy); + parupd (star1, starpar[i], ix, iy); + addstar (starpar[i], SUB, jmtype); + imtype[nstot] = 3; + parupd (star2, starpar[nstot], ix, iy); + parupd (star2, shadow[nstot], ix, iy); + addstar (starpar[nstot], SUB, jmtype); + nstot ++; + continue; + } + } + parupd (star, starpar[i], ix, iy); + + if ((thresh < tmax / 4) && (thresh > tmax/16)) { + + float err2[NPMAX], galchi2; + + beta4 = 0.01; + for (n = 0; n < NPMAX; n++) star2[n] = star[n]; + galchi2 = chisq (onestar, xs, ys, zs, dzs, npt, star2, err2, NFIT2, acc, parlim, nit); + if (1.2*galchi2 < galchi) { + jmtype = 12; + parupd (star2, starpar[i], ix, iy); + mprint (3, "extended: %f vs %f %f %f %f - %d\n", galchi/npt, galchi2/npt, starpar[i][2], starpar[i][3], star[1], jmtype); + } + beta4 = 1.0; + } + + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + verybig = galaxy (star, shaderr[i], starpar[i]); + if (jmtype == 3) verybig = verybig && (chi[4] > xtra); + + /* shouldn't this test come earlier? */ + converge = finite(galchi); + offp = offpic (star, ix, iy, &dx, &dy); + verybig = verybig && !offp && converge;; + + if (!verybig || fixxy) { + if (jmtype != 3) jmtype = 1; + if (!converge) { + if (jmtype != 3) jmtype = 9; + mprint (3, "obj %d, npts %d, %d & %d .... failed to converge\n", i, npt, ix, iy); + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + if (offp) { + jmtype = 9; + imtype[i] = fixxy ? jmtype + 100 : jmtype; + mprint (3, "obj %d, npts %d, %d & %d .... fit center outside fit subraster\n", i, npt, ix, iy); + } + addstar (starpar[i], SUB, jmtype); + imtype[i] = fixxy ? jmtype + 100 : jmtype; + continue; + } + + mprint (3, "obj %d, npts %d, %d & %d .... is very big....\n", i, npt, ix, iy); + mprint (3, ".... testing galaxy vs. double-star \n"); + + starchi = twofit (star, starpar[i], star1, star2); + + if (finite(starchi) && (starchi/galchi < stograt)) { + mprint (3, " result-> a split star: gal-chi: %f, star-chi: %f\n", galchi, starchi); + jmtype = 3; + parupd (star1, shadow[i], ix, iy); + parupd (star1, starpar[i], ix, iy); + addstar (starpar[i], SUB, jmtype); + imtype[nstot] = 3; + parupd (star2, starpar[nstot], ix, iy); + parupd (star2, shadow[nstot], ix, iy); + addstar (starpar[nstot], SUB, jmtype); + nstot ++; + } else { + mprint (3, " result-> a galaxy: gal-chi: %f, star-chi: %f\n", galchi, starchi); + jmtype = 2; + parupd (star, starpar[i], ix, iy); + if (star[1]*star[4]*star[6]*ufactor > 20000) { + float err2[NPMAX], galchi2; + + beta4 = 0.01; + for (n = 0; n < NPMAX; n++) star2[n] = star[n]; + galchi2 = chisq (onestar, xs, ys, zs, dzs, npt, star2, err2, NFIT2, acc, parlim, nit); + if (1.2*galchi2 < galchi) { + jmtype = 12; + parupd (star2, starpar[i], ix, iy); + } + /* fprintf (stderr, "extended: %f vs %f %f %f %f - %d\n", galchi/npt, galchi2/npt, starpar[i][2], starpar[i][3], star[1], jmtype); + */ + beta4 = 1.0; + } + addstar (starpar[i], SUB, jmtype); + } + imtype[i] = fixxy ? jmtype + 100 : jmtype; + } + test7 = FALSE; +} + + + +/***************** + + This routine seems really poorly written. there is too much obfuscation. + + we just need to do a few things: + + 1) do we try this object with 7-par fit? (if no, continue) + + 2) fit the object with 7-par fit + + 3) did it succeed? (if no, set some types) + + 4) should it be a type 10? (ie, very bright) + + 5) should it be a type 2? (ie, significantly extended) + + 6) should if be a type 3? (ie, split in two) + + 7) should it be a type 12? (ie, hubble fit) + +******************/ Index: /tags/gophot-1-0/Ohana/src/gophot/src/skyfun_plane.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/skyfun_plane.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/skyfun_plane.c (revision 2556) @@ -0,0 +1,22 @@ +# include "gophot.h" + +# define HALF 0.5 + +float skyfun_plane (int ix, int iy, float *a, float *fa) { + + int i; + float value; + + fa[0] = 1.0; + fa[1] = HALF*(ix - HALF*nfast)/(HALF*nfast); + fa[2] = HALF*(iy - HALF*nslow)/(HALF*nslow); + + value = 0; + for (i = 0; i < 1; i++) { + value += a[i]*fa[i]; + } + + return (value); +} + +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/sort.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/sort.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/sort.c (revision 2556) @@ -0,0 +1,77 @@ + +sort (float *value, int N) { + + int l,j,ir,i; + float temp; + + if (N < 2) return; + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + temp = value[--l]; + } + else { + temp = value[ir]; + value[ir] = value[0]; + if (--ir == 0) { + value[0] = temp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && value[j] < value[j+1]) ++j; + if (temp < value[j]) { + value[i]=value[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + value[i] = temp; + } +} + +sortkey (int *index, float *key, int N) { + + int l,j,ir,i; + float ktmp; + int itmp; + + if (N < 2) return; + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + ktmp = key[l]; + itmp = index[l]; + } + else { + ktmp = key[ir]; + key[ir] = key[0]; + itmp = index[ir]; + index[ir] = index[0]; + ir--; + if (ir == 0) { + key[0] = ktmp; + index[0] = itmp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && key[j] < key[j+1]) ++j; + if (ktmp < key[j]) { + key[i] = key[j]; + index[i] = index[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + key[i] = ktmp; + index[i] = itmp; + } +} Index: /tags/gophot-1-0/Ohana/src/gophot/src/toobright.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/toobright.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/toobright.c (revision 2556) @@ -0,0 +1,64 @@ +# include "gophot.h" + +bool toobright (float *star) { + + bool value; + bool flag; + int jrect[5]; + int nsat, nbad, ix, iy, ixhi, ixlo, iyhi, iylo, jx, jy; + float dum, duma[NPMAX]; + + value = FALSE; + + if (star[1] < itop/4) return (FALSE); + + nsat = 0; + nbad = 0; + dum = guess2 (duma, star, &ix, &iy); + ixhi = MIN ((int)(ix + krect[1]/2 + 0.5), nfast - 1); + ixlo = MAX ((int)(ix - krect[1]/2 + 0.5), 0); + iyhi = MIN ((int)(iy + krect[2]/2 + 0.5), nslow - 1); + iylo = MAX ((int)(iy - krect[2]/2 + 0.5), 0); + + for (jy = iylo; jy <= iyhi; jy++) { + for (jx = ixlo; jx <= ixhi; jx++) { + if (!finite (noise[jx+jy*nfast])) nbad ++; + if (big[jx+jy*nfast] >= itop) nsat ++; + } + } + + value = (star[1] > cmax) || (nsat >= icrit) || (nbad >= icrit); + + mprint (3, "object at: %d, %d - %d: %d %d %f %f %f\n", ix,iy, value, nsat, nbad, star[1], big[ix+iy*nfast], noise[ix+iy*nfast]); + + if (value) { + if (nsat >= icrit) star[1] = ctpersat*nsat; + + oblims(star,jrect); + + star[4] = jrect[2] - jrect[1]; + star[5] = nsat; + star[6] = jrect[4] - jrect[3]; + + mprint (0, "%d, %d funny pixels in object at %d, %d (%f peak)\n", nsat, nbad, ix, iy, star[1]); + } + + return (value); +} + +bool verybright (float *star) { + + float flux; + + /* rough scaling to compare fits with different scale */ + + flux = star[1]*sqrt((star[4]*star[6])/(ava[4]*ava[6]));; + + if (flux > cmax) { + return (TRUE); + } else { + return (FALSE); + } + +} + Index: /tags/gophot-1-0/Ohana/src/gophot/src/toofaint.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/toofaint.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/toofaint.c (revision 2556) @@ -0,0 +1,22 @@ +# include "gophot.h" + +bool toofaint (float *star, float *err) { + + bool value; + float cmin, sig2; + + cmin = 1.0; + value = (star[1] < cmin); + + if (err[1] > 0) { + sig2 = SQ(star[1]/ufactor) / err[1]; + value = value || (sig2 < crit7); + mprint (3, "sig2 = %f (%f %f %f)\n", sig2, star[1], ufactor, err[1]); + mprint (3, "value, crit7: %d, %f\n", value, crit7); + } + + return (value); + +} + +/* this function uses C 0,N-1 for a[], fa[] */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/transmask.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/transmask.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/transmask.c (revision 2556) @@ -0,0 +1,94 @@ +# include "gophot.h" + +bool transmask (int ix, int iy, float sky) { + + bool value; + int i, j, ii, jj; + float tmp, tmp0, tmp1, tmp2, ratio, hump2; + float tmpB, tmpN, df; + float *Bval, *Nval, cval, tmp3; + + value = FALSE; + + if (test7) + hump2 = crit7; + else + hump2 = SQ (bumpcrit); + + tmp0 = tmp1 = tmp2 = 0; + tmpB = tmpN = tmp3 = 0; + + cval = big[ix + iy*nfast]; + + jj = iy - iyby2; + for (j = 0; j < 2*iyby2 + 1; j++, jj++) { + if (jj < 0) continue; + if (jj > nslow - 1) continue; + + ii = ix - ixby2; + Bval = &big[ii + jj*nfast]; + Nval = &noise[ii + jj*nfast]; + for (i = 0; i < 2*ixby2 + 1; i++, ii++, Bval++, Nval++) { + if (ii < 0) continue; + if (ii > nfast - 1) continue; + if (!finite (*Nval)) continue; + +# if (1) + tmp0 ++; + tmp = starmask[i][j] / (*Nval + *Bval); + tmp1 += tmp*starmask[i][j]; + tmp2 += tmp*(*Bval-sky); + +# else + + tmp0 ++; + tmp1 += (*Nval + *Bval); + tmp2 += (*Bval - sky); + + df = (cval - *Bval); + tmp3 += df * fabs(df) / (*Nval + *Bval); + +# endif + + tmpB += (*Bval - sky); + tmpN += (*Nval - rnoise); + } + } + + if (tmp0 == 0) return (FALSE); + + if (tmpB < tmpN) return (FALSE); + + if (tmp2 > 0) { + ratio = tmp2*tmp2/tmp1; + +# if (1) + if (ratio > hump2) { + value = TRUE; + mprint (3, " trigger 1 on new object: %d, %d, %f %f %f %f %f\n", ix, iy, ratio, tmp3, sky, tmp2, tmp1); + if (ratio < 1.1*hump2) mprint (3, "marginal: (s/n)**2 through mask = %f\n", ratio); + } +# else + if ((ratio > hump2) && (tmp3 > 0.25*ratio)) { + value = TRUE; + mprint (3, " trigger 1 on new object: %d, %d, %f %f %f %f %f\n", ix, iy, ratio, tmp3, sky, tmp2, tmp1); + if (ratio < 1.1*hump2) mprint (3, "marginal: (s/n)**2 through mask = %f\n", ratio); + } +# endif + + } + + if (value) mprint (3, " trigger 2 on new object: %d, %d\n", ix, iy); + + return (value); + +} + + + /* t1 = 1/N + t2 = d/N + rat = N d^2 / N^2 = d^2 / N + */ + +/* does not check if bump is significant in region, + only if mask region is significant over sky */ Index: /tags/gophot-1-0/Ohana/src/gophot/src/twofit.c =================================================================== --- /tags/gophot-1-0/Ohana/src/gophot/src/twofit.c (revision 2556) +++ /tags/gophot-1-0/Ohana/src/gophot/src/twofit.c (revision 2556) @@ -0,0 +1,127 @@ +# include "gophot.h" + +/**************************************** + there are big problems here with the values of b and probably + other things. fix please! */ + +/* we are fitting two stars at the location of the input star */ +/* starraw has the real X,Y coords, star is relative to centroid */ + +float twofit (float *star, float *starraw, float *star1, float *star2) { + + int i; + float b[2*NPMAX], err[2*NPMAX], bacc[2*NPMAX], bparlim[2*NPMAX], tstar[NPMAX]; + bool badfit, conv; + float a5, a7, angle, root, root1, root2, dx2, dy2, dx, dy, value; + float garea, sarea, dx74, dy74, dot, fac1, fac2, dxmax, dymax; + + badfit = FALSE; + conv = TRUE; + + /* fills in values 1, 5, 6, 7, 8 */ + parinterp (starraw[2], starraw[3], tstar); + + /* this all seems way too obtuse. */ + a5 = 1.0/star[4]; + a7 = 1.0/star[6]; + angle = atan2 (-2*star[5], a7 - a5)/2.0; + root = sqrt (SQ(a5 - a7) + 4*SQ(star[5])); + root1 = (a5 + a7 + root)/2.0; + root2 = root1 - root; + dx2 = star[4] - tstar[4]; + dy2 = star[6] - tstar[6]; + dx = sqrt (MAX (dx2, 0.0)) / 2.0; + dy = sqrt (MAX (dy2, 0.0)) / 2.0; + badfit = (MAX (dx, dy) == 0); + /* we can use these since dx & dy are always >= 0 */ + dx = (cos(angle) < 0) ? -dx : dx; + dy = (sin(angle) < 0) ? -dy : dy; + garea = 1.0 / sqrt (fabs (root1*root2)); + sarea = sqrt (fabs (tstar[4]*tstar[6])); + dx74 = starraw[2] - star[2]; + dy74 = starraw[3] - star[3]; + dot = dx74*dx + dy74*dy; + if (dot > 0) { + fac1 = 0.6666666; + fac2 = 1.3333333; + } else { + fac1 = 1.3333333; + fac2 = 0.6666666; + } + + /* + b[0] - sky + b[1] - I1 + b[2], b[3] - X1, Y1 + b[4] - I2 + b[5], b[6] - X2, Y2 + b[7], b[8], b[9] - shape + */ + + b[0] = star[0]; + b[1] = (star[1]*garea*fac2 / (sarea*2)); + b[2] = star[2] - dx*fac1; + b[3] = star[3] - dy*fac1; + b[4] = (star[1]*garea*fac1 / (sarea*2)); + b[5] = star[2] + dx*fac2; + b[6] = star[3] + dy*fac2; + + /* for reasons unclear, the fit is done in ln(Io), not Io */ + b[1] = log(b[1]); + b[4] = log(b[4]); + + b[7] = tstar[4]; + b[8] = tstar[5]; + b[9] = tstar[6]; + + bacc[0] = acc[0]; + bacc[1] = bacc[4] = acc[1]; + bacc[2] = bacc[5] = acc[2]; + bacc[3] = bacc[6] = acc[3]; + + bparlim[0] = parlim[0]; + bparlim[1] = bparlim[4] = parlim[1]; + bparlim[2] = bparlim[5] = parlim[2]; + bparlim[3] = bparlim[6] = parlim[3]; + + dxmax = MAX (fabs(b[2]), fabs(b[5])); + dymax = MAX (fabs(b[3]), fabs(b[6])); + badfit = badfit || (dxmax > irect[1]/2.0); + badfit = badfit || (dymax > irect[2]/2.0); + + if (!badfit) { + value = chisq (twostar, xs, ys, zs, dzs, npt, b, err, NFIT2, bacc, bparlim, 2*nit); + conv = finite (value); + } + + dxmax = MAX (fabs(b[2]), fabs(b[5])); + dymax = MAX (fabs(b[3]), fabs(b[6])); + badfit = badfit || (dxmax > 0.4*irect[1]); + badfit = badfit || (dymax > 0.4*irect[2]); + + if (conv) { + if (b[1] < b[4]) { + SWAP (b[1], b[4]); + SWAP (b[2], b[5]); + SWAP (b[3], b[6]); + } + + star1[0] = star2[0] = b[0]; + star1[4] = star2[4] = b[7]; + star1[5] = star2[5] = b[8]; + star1[6] = star2[6] = b[9]; + star1[1] = exp(b[1]); + star1[2] = b[2]; + star1[3] = b[3]; + star2[1] = exp(b[4]); + star2[2] = b[5]; + star2[3] = b[6]; + } + + badfit = badfit || centertest (star1, xs, ys); + badfit = badfit || centertest (star2, xs, ys); + if (!conv || badfit) value = MAGIC; + + return (value); +} +/* this function uses C 0,N-1 for a[], fa[] */