Changeset 23924 for trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c

Timestamp:

Apr 18, 2009, 11:11:15 AM (17 years ago)

Author:

eugene

Message:

import updates from camera group (fits table outputs, fixes to persistent trails); modify import to avoid libpsf (fortran)

File:

• trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c (modified) (3 diffs)

trunk/extsrc/gpcsw/gpcsrc/fits/burntool/psfstats.c

@@ -1 @@
-/* psfstats.c - calculate some statistics about a PSF */
+/* dummy version of psfstats.c - disable use of the psf */
 #include <stdio.h>
 #include <stdlib.h>
@@ -12 +12 @@
 #include "psf/psf.h"
 
-//#define PSFSTATS              /* ignore psf stats so we can skip fortan */
-//#define TEST          /* Dump all the PSF info? */
-
-static unsigned short int *psfbuf=NULL;
-static int npsfbuf=0;
-#define PHONYSKY 1000
-
 /****************************************************************/
 /* psf_stats(): Tell us about this PSF star */
@@ -24 +17 @@
                      double *fwhm, double *q)
 {
-
-#ifdef PSFSTATS
-
-   int i, j, r, r2, err, i0, j0;
-   double d, qxx, qxy, qyy, mx, my, sum;
-   PSF_PARAM psfout;      /* Results of fit */
-   PSF_EXTRA psfextra;    /* Extra results from fit */
-   double pi=4*atan(1.0);
-
-   if(nx*ny > npsfbuf) {
-      if(psfbuf != NULL) free(psfbuf);
-      psfbuf = (ushort *)calloc(nx*ny, sizeof(ushort));
-      npsfbuf = nx*ny;
-   }
-   for(i=0; i<nx*ny; i++) {
-      j = ((int)data[i]) - bias + PHONYSKY;
-      if(j >= 0 && j < 65535) psfbuf[i] = j;
-      else psfbuf[i] = 0;
-   }
-
-/* 2D PSF fit */
-   err = psf(nx, ny, psfbuf, PSF_2DIM, NULL, &psfout, &psfextra);
-
-#ifdef TEST
-   printf("err = %d\n", err);
-#endif
-
-   if(!err) {
-#ifdef TEST
-      printf("err = %d\n", err);
-      printf("bias = %d\n", bias);
-      printf("ix, iy = %d %d\n", psfout.ix, psfout.iy);
-      printf("x0, y0 = %.2f %.2f\n", psfout.x0, psfout.y0);
-      printf("fwhm = %.2f\n", psfout.fwhm);
-      printf("peak, bckgnd = %.1f %.1f\n", psfout.peak, psfout.bkgnd);
-      printf("flux, S/N = %.1f %.2f\n", psfout.flux, psfout.sn);
-      printf("weight = %.2f\n", psfout.weight);
-      printf("FWmaj, FWmin = %.2f %.2f\n", psfextra.majfw, psfextra.minfw);
-      printf("theta = %.1f\n", psfextra.thfw*180/pi);
-      printf("wpeak, wbkgnd = %.1f %.1f\n", psfextra.wpeak, psfextra.wbkgnd);
-      printf("dflux, dbckgnd = %.1f %.1f\n", psfextra.dflux, psfextra.dbkgnd);
-      printf("rmsbkgnd = %.1f\n", psfextra.rmsbkgnd);
-#endif
-      fwhm[0] = psfextra.majfw;
-      fwhm[1] = psfextra.minfw;
-      fwhm[2] = psfextra.thfw;
-
-/* Moments */
-      i0 = MIN(psfout.ix, nx-psfout.ix);
-      j0 = MIN(psfout.iy, ny-psfout.iy);
-   } else {
-      fwhm[0] = fwhm[1] = fwhm[2] = 0.0;
-      i0 = nx/2;
-      j0 = ny/2;
-   }
-
-   r = MIN(i0, j0) - 1;
-   sum = mx = my = qxx = qxy = qyy = 0.0;
-   for(j=j0-r; j<j0+r; j++) {
-      for(i=i0-r; i<i0+r; i++) {
-         r2 = (i-i0)*(i-i0) + (j-j0)*(j-j0);
-         if(r2 > r*r) continue;
-         d = (double)data[i+j*nx] - bias;
-         sum += d;
-         mx += (i+0.5) * d;
-         my += (j+0.5) * d;
-         qxx += (i+0.5) * (i+0.5) * d;
-         qxy += (i+0.5) * (j+0.5) * d;
-         qyy += (j+0.5) * (j+0.5) * d;
-      }
-   }
-   if(sum > 0) {
-      mx /= sum;
-      my /= sum;
-      qxx = qxx / sum - mx*mx;
-      qxy = qxy / sum - mx*my;
-      qyy = qyy / sum - my*my;
-   } else {
-      q[0] = q[1] = q[2] = 0.0;
-   }
-
-   q[0] = (qxx + qyy) / 2;
-   q[1] = qxx - qyy;
-   q[2] = 2 * qxy;
-
-/* phi is the angular position of the OTA around the center [rad] */
-//   qtang = (qxx - qyy) * cos(2*phi) + 2 * qxy * sin(2*phi);
-
-#ifdef TEST
-   printf("Sum= %.1f mx= %.2f my= %.2f qxx= %.4f qxy= %.4f qyy= %.4f\n", 
-          sum, mx, my, qxx, qxy, qyy);
-#endif
-
-# else
-   printf("psfstats ignored in this build\n");
-   return (1);
-#endif /* PSFSTATS */
-   return(0);
+   return(1);
 }