Commit 1d1d4318 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Fixed BITPIX display info.

Changed CONCENTRATION_MODEL to SPREAD_MODEL.
Added SPREADERR_MODEL measurement (currently a crude estimate).
Optimized FFT calls.
Set PSF_NMAX to 2 by default.
parent ff13ac00
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: functions for output of catalog data. * Contents: functions for output of catalog data.
* *
* Last modify: 24/09/2009 * Last modify: 29/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -175,6 +175,7 @@ void updateparamflags() ...@@ -175,6 +175,7 @@ void updateparamflags()
| FLAG(obj2.prof_disk_patternspiral); | FLAG(obj2.prof_disk_patternspiral);
FLAG(obj2.prof_disk_scale) |= prefs.pattern_flag; FLAG(obj2.prof_disk_scale) |= prefs.pattern_flag;
FLAG(obj2.prof_concentration) |= FLAG(obj2.prof_concentrationerr);
FLAG(obj2.poserraw_prof) |= FLAG(obj2.poserrbw_prof); FLAG(obj2.poserraw_prof) |= FLAG(obj2.poserrbw_prof);
FLAG(obj2.poserrcxxw_prof) |= FLAG(obj2.poserrcyyw_prof) FLAG(obj2.poserrcxxw_prof) |= FLAG(obj2.poserrcyyw_prof)
| FLAG(obj2.poserrcxyw_prof); | FLAG(obj2.poserrcxyw_prof);
...@@ -210,8 +211,8 @@ void updateparamflags() ...@@ -210,8 +211,8 @@ void updateparamflags()
| FLAG(obj2.prof_concentration) | FLAG(obj2.prof_concentration)
| FLAG(obj2.prof_class_star); | FLAG(obj2.prof_class_star);
FLAG(obj2.mag_prof) |= FLAG(obj2.magerr_prof); FLAG(obj2.fluxerr_prof) |= FLAG(obj2.magerr_prof)
FLAG(obj2.magerr_prof) |= FLAG(obj2.fluxerr_prof); | FLAG(obj2.prof_concentrationerr);
FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof) | FLAG(obj2.fluxerr_prof); FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof) | FLAG(obj2.fluxerr_prof);
FLAG(obj2.prof_spheroid_mag) |= FLAG(obj2.prof_spheroid_magerr); FLAG(obj2.prof_spheroid_mag) |= FLAG(obj2.prof_spheroid_magerr);
FLAG(obj2.prof_spheroid_reff) |= FLAG(obj2.prof_spheroid_refferr); FLAG(obj2.prof_spheroid_reff) |= FLAG(obj2.prof_spheroid_refferr);
...@@ -612,10 +613,10 @@ void dumpparams(void) ...@@ -612,10 +613,10 @@ void dumpparams(void)
for (i=0; *objkey[i].name ; i++) for (i=0; *objkey[i].name ; i++)
if (*objkey[i].unit) if (*objkey[i].unit)
printf("#%-22.22s %-58.58s [%s]\n", printf("#%-24.24s %-57.57s [%s]\n",
objkey[i].name, objkey[i].comment, objkey[i].unit); objkey[i].name, objkey[i].comment, objkey[i].unit);
else else
printf("#%-22.22s %-58.58s\n", printf("#%-24.24s %-57.57s\n",
objkey[i].name, objkey[i].comment); objkey[i].name, objkey[i].comment);
return; return;
......
...@@ -32,11 +32,12 @@ ...@@ -32,11 +32,12 @@
#include "threads.h" #include "threads.h"
#endif #endif
fftwf_plan fplan, bplan;
int firsttimeflag; int firsttimeflag;
#ifdef USE_THREADS #ifdef USE_THREADS
pthread_mutex_t fftmutex; pthread_mutex_t fftmutex;
#endif #endif
fftwf_complex *fdata1;
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
/****** fft_init ************************************************************ /****** fft_init ************************************************************
...@@ -98,6 +99,30 @@ void fft_end(void) ...@@ -98,6 +99,30 @@ void fft_end(void)
} }
/****** fft_reset ***********************************************************
PROTO void fft_reset(void)
PURPOSE Reset the FFT plans
INPUT -.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 26/06/2009
***/
void fft_reset(void)
{
if (fplan)
{
QFFTWFREE(fdata1);
fftwf_destroy_plan(fplan);
}
if (bplan)
fftwf_destroy_plan(bplan);
fplan = bplan = NULL;
return;
}
/****** fft_conv ************************************************************ /****** fft_conv ************************************************************
PROTO void fft_conv(float *data1, float *fdata2, int *size) PROTO void fft_conv(float *data1, float *fdata2, int *size)
PURPOSE Optimized 2-dimensional FFT convolution using the FFTW library. PURPOSE Optimized 2-dimensional FFT convolution using the FFTW library.
...@@ -112,8 +137,7 @@ VERSION 26/03/2007 ...@@ -112,8 +137,7 @@ VERSION 26/03/2007
***/ ***/
void fft_conv(float *data1, float *fdata2, int *size) void fft_conv(float *data1, float *fdata2, int *size)
{ {
fftwf_plan plan; float *fdata1p,*fdata2p,
float *fdata1,*fdata1p,*fdata2p,
real,imag, fac; real,imag, fac;
int i, npix,npix2; int i, npix,npix2;
...@@ -125,18 +149,23 @@ void fft_conv(float *data1, float *fdata2, int *size) ...@@ -125,18 +149,23 @@ void fft_conv(float *data1, float *fdata2, int *size)
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex); QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif #endif
QFFTWMALLOC(fdata1, float, npix2); if (!fplan)
plan = fftwf_plan_dft_r2c_2d(size[1], size[0], data1, {
(fftwf_complex *)fdata1, FFTW_ESTIMATE|FFTW_DESTROY_INPUT); QFFTWMALLOC(fdata1, fftwf_complex, npix2);
fplan = fftwf_plan_dft_r2c_2d(size[1], size[0], data1,
(fftwf_complex *)fdata1, FFTW_MEASURE|FFTW_DESTROY_INPUT);
}
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
fftwf_execute(plan); fftwf_execute_dft_r2c(fplan, data1, fdata1);
// fftwf_execute(plan);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex); QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif #endif
fftwf_destroy_plan(plan); // fftwf_destroy_plan(plan);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
...@@ -157,20 +186,21 @@ void fft_conv(float *data1, float *fdata2, int *size) ...@@ -157,20 +186,21 @@ void fft_conv(float *data1, float *fdata2, int *size)
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex); QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif #endif
plan = fftwf_plan_dft_c2r_2d(size[1], size[0], (fftwf_complex *)fdata1, if (!bplan)
data1, FFTW_ESTIMATE|FFTW_DESTROY_INPUT); bplan = fftwf_plan_dft_c2r_2d(size[1], size[0], (fftwf_complex *)fdata1,
data1, FFTW_MEASURE|FFTW_DESTROY_INPUT);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
fftwf_execute_dft_c2r(bplan, fdata1, data1);
fftwf_execute(plan); // fftwf_execute(plan);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_LOCK(&fftmutex); QPTHREAD_MUTEX_LOCK(&fftmutex);
#endif #endif
fftwf_destroy_plan(plan); // fftwf_destroy_plan(plan);
/* Free the fdata1 scratch array */ /* Free the fdata1 scratch array */
QFFTWFREE(fdata1);
#ifdef USE_THREADS #ifdef USE_THREADS
QPTHREAD_MUTEX_UNLOCK(&fftmutex); QPTHREAD_MUTEX_UNLOCK(&fftmutex);
#endif #endif
......
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
/*---------------------------------- protos --------------------------------*/ /*---------------------------------- protos --------------------------------*/
extern void fft_conv(float *data1, float *fdata2, int *size), extern void fft_conv(float *data1, float *fdata2, int *size),
fft_end(), fft_end(),
fft_init(int nthreads); fft_init(int nthreads),
fft_reset(void);
extern float *fft_rtf(float *data, int *size); extern float *fft_rtf(float *data, int *size);
...@@ -113,7 +113,7 @@ void makeit() ...@@ -113,7 +113,7 @@ void makeit()
changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1); changecatparamarrays("VECTOR_MODEL", &theprofit->nparam, 1);
changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1); changecatparamarrays("VECTOR_MODELERR", &theprofit->nparam, 1);
nparam2[0] = nparam2[1] = theprofit->nparam; nparam2[0] = nparam2[1] = theprofit->nparam;
changecatparamarrays("MATRIX_MODELERR", &nparam2, 2); changecatparamarrays("MATRIX_MODELERR", nparam2, 2);
if (prefs.pattern_flag) if (prefs.pattern_flag)
{ {
npat = prefs.prof_disk_patternvectorsize; npat = prefs.prof_disk_patternvectorsize;
......
...@@ -829,7 +829,7 @@ keystruct objkey[] = { ...@@ -829,7 +829,7 @@ keystruct objkey[] = {
{"ERRBPSF_IMAGE", "PSF RMS position error along minor axis", {"ERRBPSF_IMAGE", "PSF RMS position error along minor axis",
&outobj2.poserrb_psf, H_FLOAT, T_FLOAT, "%8.4f", "pixel", &outobj2.poserrb_psf, H_FLOAT, T_FLOAT, "%8.4f", "pixel",
"stat.stdev;stat.min;pos.errorEllipse;instr.det", "pix"}, "stat.stdev;stat.min;pos.errorEllipse;instr.det", "pix"},
{"ERRTHTPSF_IMAGE", "PSF error ellipse position angle (CCW/x)", {"ERRTHETAPSF_IMAGE", "PSF error ellipse position angle (CCW/x)",
&outobj2.poserrtheta_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg", &outobj2.poserrtheta_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
"pos.posAng;pos.errorEllipse;instr.det", "deg"}, "pos.posAng;pos.errorEllipse;instr.det", "deg"},
{"ERRAPSF_WORLD", "World PSF RMS position error along major axis", {"ERRAPSF_WORLD", "World PSF RMS position error along major axis",
...@@ -838,16 +838,16 @@ keystruct objkey[] = { ...@@ -838,16 +838,16 @@ keystruct objkey[] = {
{"ERRBPSF_WORLD", "World PSF RMS position error along minor axis", {"ERRBPSF_WORLD", "World PSF RMS position error along minor axis",
&outobj2.poserrbw_psf, H_FLOAT, T_FLOAT, "%12.7g", "pixel", &outobj2.poserrbw_psf, H_FLOAT, T_FLOAT, "%12.7g", "pixel",
"stat.stdev;stat.min;pos.errorEllipse", "deg"}, "stat.stdev;stat.min;pos.errorEllipse", "deg"},
{"ERRTHTPSF_WORLD", "PSF error ellipse pos. angle (CCW/world-x)", {"ERRTHETAPSF_WORLD", "PSF error ellipse pos. angle (CCW/world-x)",
&outobj2.poserrthetaw_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg", &outobj2.poserrthetaw_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
"pos.posAng;pos.errorEllipse", "deg"}, "pos.posAng;pos.errorEllipse", "deg"},
{"ERRTHTPSF_SKY", "Native PSF error ellipse pos. angle (east of north)", {"ERRTHETAPSF_SKY", "Native PSF error ellipse pos. angle (east of north)",
&outobj2.poserrthetas_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg", &outobj2.poserrthetas_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
"pos.posAng;pos.errorEllipse", "deg"}, "pos.posAng;pos.errorEllipse", "deg"},
{"ERRTHTPSF_J2000", "J2000 PSF error ellipse pos. angle (east of north)", {"ERRTHETAPSF_J2000", "J2000 PSF error ellipse pos. angle (east of north)",
&outobj2.poserrtheta2000_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg", &outobj2.poserrtheta2000_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
"pos.posAng;pos.errorEllipse", "deg"}, "pos.posAng;pos.errorEllipse", "deg"},
{"ERRTHTPSF_B1950", "B1950 PSF error ellipse pos. angle (east of north)", {"ERRTHETAPSF_B1950", "B1950 PSF error ellipse pos. angle (east of north)",
&outobj2.poserrtheta1950_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg", &outobj2.poserrtheta1950_psf, H_FLOAT, T_FLOAT, "%6.2f", "deg",
"pos.posAng;pos.errorEllipse", "deg"}, "pos.posAng;pos.errorEllipse", "deg"},
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Model-fitting parameter list for catalog data. * Contents: Model-fitting parameter list for catalog data.
* *
* Last modify: 24/09/2009 * Last modify: 29/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -266,8 +266,11 @@ ...@@ -266,8 +266,11 @@
&outobj2.prof_theta1950, H_FLOAT, T_FLOAT, "%+6.2f", "deg", &outobj2.prof_theta1950, H_FLOAT, T_FLOAT, "%+6.2f", "deg",
"pos.posAng;stat.fit", "deg"}, "pos.posAng;stat.fit", "deg"},
{"CONCENTRATION_MODEL", "Concentration parameter from model-fitting", {"SPREAD_MODEL", "Spread parameter from model-fitting",
&outobj2.prof_concentration, H_FLOAT, T_FLOAT, "%7.4f", "", &outobj2.prof_concentration, H_FLOAT, T_FLOAT, "%8.5f", "",
"src.morph.param", ""},
{"SPREADERR_MODEL", "Spread parameter error from model-fitting",
&outobj2.prof_concentrationerr, H_FLOAT, T_FLOAT, "%8.5f", "",
"src.morph.param", ""}, "src.morph.param", ""},
{"CLASS_STAR_MODEL", "S/G classifier from model-fitting", {"CLASS_STAR_MODEL", "S/G classifier from model-fitting",
&outobj2.prof_class_star, H_FLOAT, T_FLOAT, "%7.4f", "", &outobj2.prof_class_star, H_FLOAT, T_FLOAT, "%7.4f", "",
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Keywords for the configuration file. * Contents: Keywords for the configuration file.
* *
* Last modify: 20/11/2008 * Last modify: 24/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -287,7 +287,7 @@ char *default_prefs[] = ...@@ -287,7 +287,7 @@ char *default_prefs[] =
"*#--------------------------- Experimental Stuff -----------------------------", "*#--------------------------- Experimental Stuff -----------------------------",
"*", "*",
"*PSF_NAME default.psf # File containing the PSF model", "*PSF_NAME default.psf # File containing the PSF model",
"*PSF_NMAX 9 # Max.number of PSFs fitted simultaneously", "*PSF_NMAX 2 # Max.number of PSFs fitted simultaneously",
"*PSFDISPLAY_TYPE SPLIT # Catalog type for PSF-fitting: SPLIT or VECTOR", "*PSFDISPLAY_TYPE SPLIT # Catalog type for PSF-fitting: SPLIT or VECTOR",
"*PATTERN_TYPE RINGS-HARMONIC # can RINGS-QUADPOLE, RINGS-OCTOPOLE,", "*PATTERN_TYPE RINGS-HARMONIC # can RINGS-QUADPOLE, RINGS-OCTOPOLE,",
"* # RINGS-HARMONICS or GAUSS-LAGUERRE", "* # RINGS-HARMONICS or GAUSS-LAGUERRE",
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Fit an arbitrary profile combination to a detection. * Contents: Fit an arbitrary profile combination to a detection.
* *
* Last modify: 24/09/2009 * Last modify: 29/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -49,12 +49,22 @@ static void make_kernel(float pos, float *kernel, interpenum interptype); ...@@ -49,12 +49,22 @@ static void make_kernel(float pos, float *kernel, interpenum interptype);
/*------------------------------- variables ---------------------------------*/ /*------------------------------- variables ---------------------------------*/
char profname[][32]={"background offset", "Sersic spheroid", const char profname[][32]={"background offset", "Sersic spheroid",
"De Vaucouleurs spheroid", "exponential disk", "spiral arms", "De Vaucouleurs spheroid", "exponential disk", "spiral arms",
"bar", "inner ring", "outer ring", "tabulated model", "bar", "inner ring", "outer ring", "tabulated model",
""}; ""};
int interp_kernwidth[5]={1,2,4,6,8}; const int interp_kernwidth[5]={1,2,4,6,8};
const int flux_flag[PARAM_NPARAM] = {0,0,0,
1,0,0,0,0,
1,0,0,0,
1,0,0,0,0,0,0,0,
1,0,0,
1,0,0,
1,0,0
};
int theniter, the_gal; int theniter, the_gal;
/* "Local" global variables; it seems dirty but it simplifies a lot */ /* "Local" global variables; it seems dirty but it simplifies a lot */
/* interfacing to the LM routines */ /* interfacing to the LM routines */
...@@ -175,8 +185,11 @@ void profit_fit(profitstruct *profit, ...@@ -175,8 +185,11 @@ void profit_fit(profitstruct *profit,
psfstruct *psf; psfstruct *psf;
checkstruct *check; checkstruct *check;
double emx2,emy2,emxy, a , cp,sp, cn, bn, n; double emx2,emy2,emxy, a , cp,sp, cn, bn, n;
float psf_fwhm, dchi2, err; float **list,
int i,j,p, nparam, nparam2, ncomp, nprof; *cov,
psf_fwhm, dchi2, err;
int *index,
i,j,p, nparam, nparam2, ncomp, nprof;
nparam = profit->nparam; nparam = profit->nparam;
nparam2 = nparam*nparam; nparam2 = nparam*nparam;
...@@ -242,7 +255,8 @@ void profit_fit(profitstruct *profit, ...@@ -242,7 +255,8 @@ void profit_fit(profitstruct *profit,
profit->modnaxisn[0] = profit->modnaxisn[1]; profit->modnaxisn[0] = profit->modnaxisn[1];
/* Allocate memory for the complete model */ /* Allocate memory for the complete model */
QCALLOC(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QFFTWMALLOC(profit->modpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
memset(profit->modpix, 0, profit->modnaxisn[0]*profit->modnaxisn[1]*sizeof(float));
QMALLOC(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QMALLOC(profit->psfpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
/* Allocate memory for the partial model */ /* Allocate memory for the partial model */
QMALLOC(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]); QMALLOC(profit->pmodpix, float, profit->modnaxisn[0]*profit->modnaxisn[1]);
...@@ -259,6 +273,7 @@ void profit_fit(profitstruct *profit, ...@@ -259,6 +273,7 @@ void profit_fit(profitstruct *profit,
//the_gal++; //the_gal++;
/* Actual minimisation */ /* Actual minimisation */
fft_reset();
profit->niter = profit_minimize(profit, PROFIT_MAXITER); profit->niter = profit_minimize(profit, PROFIT_MAXITER);
profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi); profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi);
...@@ -307,20 +322,17 @@ void profit_fit(profitstruct *profit, ...@@ -307,20 +322,17 @@ void profit_fit(profitstruct *profit,
if (FLAG(obj2.fluxerr_prof)) if (FLAG(obj2.fluxerr_prof))
{ {
err = 0.0; err = 0.0;
i = j = 0; /* avoid gcc -Wall warning */ cov = profit->covar;
if (profit->paramlist[PARAM_DISK_FLUX]) index = profit->paramindex;
{ list = profit->paramlist;
i = profit->paramindex[PARAM_DISK_FLUX]; for (i=0; i<PARAM_NPARAM; i++)
err += profit->covar[i*(nparam+1)]; if (flux_flag[i] && list[i])
}
if (profit->paramlist[PARAM_SPHEROID_FLUX])
{ {
j = profit->paramindex[PARAM_SPHEROID_FLUX]; cov = profit->covar + nparam*index[i];
err += profit->covar[j*(nparam+1)]; for (j=0; j<PARAM_NPARAM; j++)
if (flux_flag[j] && list[j])
err += cov[index[j]];
} }
if (profit->paramlist[PARAM_DISK_FLUX]
&& profit->paramlist[PARAM_SPHEROID_FLUX])
err += profit->covar[i+j*nparam]+profit->covar[j+i*nparam];
obj2->fluxerr_prof = err>0.0? sqrt(err): 0.0; obj2->fluxerr_prof = err>0.0? sqrt(err): 0.0;
} }
...@@ -626,6 +638,9 @@ void profit_fit(profitstruct *profit, ...@@ -626,6 +638,9 @@ void profit_fit(profitstruct *profit,
obj2->prof_concentration = 99.0; obj2->prof_concentration = 99.0;
else if (pprofit.flux > 0.0) else if (pprofit.flux > 0.0)
obj2->prof_concentration = -99.0; obj2->prof_concentration = -99.0;
if (FLAG(obj2.prof_concentrationerr))
obj2->prof_concentrationerr = (obj2->flux_prof > 0.0?
1.086*(obj2->fluxerr_prof / obj2->flux_prof) : 99.0);
} }
prof_end(pprofit.prof[0]); prof_end(pprofit.prof[0]);
free(pprofit.prof); free(pprofit.prof);
...@@ -861,11 +876,11 @@ int profit_minimize(profitstruct *profit, int niter) ...@@ -861,11 +876,11 @@ int profit_minimize(profitstruct *profit, int niter)
profit_boundtounbound(profit, profit->paraminit); profit_boundtounbound(profit, profit->paraminit);
/* Perform fit */ /* Perform fit */
lm_opts[0] = 1.0e-4; lm_opts[0] = 1.0e-3;
lm_opts[1] = 1.0e-18; lm_opts[1] = 1.0e-17;
lm_opts[2] = 1.0e-18; lm_opts[2] = 1.0e-17;
lm_opts[3] = 1.0e-18; lm_opts[3] = 1.0e-17;
lm_opts[4] = 1.0e-4; lm_opts[4] = 1.0e-6;
niter = slevmar_dif(profit_evaluate, profit->paraminit, profit->resi, niter = slevmar_dif(profit_evaluate, profit->paraminit, profit->resi,
n, m, niter, lm_opts, info, NULL, profit->covar, profit); n, m, niter, lm_opts, info, NULL, profit->covar, profit);
...@@ -2259,13 +2274,13 @@ INPUT Profile structure, ...@@ -2259,13 +2274,13 @@ INPUT Profile structure,
OUTPUT Corrected flux contribution. OUTPUT Corrected flux contribution.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 21/09/2009 VERSION 29/09/2009
***/ ***/
float prof_add(profstruct *prof, profitstruct *profit) float prof_add(profstruct *prof, profitstruct *profit)
{ {
double xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n; double xscale, yscale, saspect, ctheta,stheta, flux, scaling, bn, n;
float posin[PROFIT_MAXEXTRA], posout[2], dnaxisn[2], float posin[PROFIT_MAXEXTRA], posout[2], dnaxisn[2],
*pixin, *pixout, *pixin, *pixin2, *pixout,
fluxfac, amp,cd11,cd12,cd21,cd22, dcd11,dcd21, dx1,dx2, fluxfac, amp,cd11,cd12,cd21,cd22, dcd11,dcd21, dx1,dx2,
x1,x10,x2, x1cin,x2cin, x1cout,x2cout, x1max,x2max, x1,x10,x2, x1cin,x2cin, x1cout,x2cout, x1max,x2max,
x1in,x2in, odx, ostep, x1in,x2in, odx, ostep,
...@@ -2276,7 +2291,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2276,7 +2291,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
r2max1, r2max2, r2min, invr2xdif, r2max1, r2max2, r2min, invr2xdif,
val, theta, thresh, ra,rb,rao, num; val, theta, thresh, ra,rb,rao, num;
int npix, noversamp, threshflag, int npix, noversamp, threshflag,
d,e,i, ix1,ix2, idx1,idx2; d,e,i, ix1,ix2, idx1,idx2, nx2, npix2;
npix = profit->modnaxisn[0]*profit->modnaxisn[1]; npix = profit->modnaxisn[0]*profit->modnaxisn[1];
...@@ -2311,6 +2326,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2311,6 +2326,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
x1cout = (float)(profit->modnaxisn[0]/2); x1cout = (float)(profit->modnaxisn[0]/2);
x2cout = (float)(profit->modnaxisn[1]/2); x2cout = (float)(profit->modnaxisn[1]/2);
nx2 = profit->modnaxisn[1]/2 + 1;
/*-- Compute the largest r^2 that fits in the frame */ /*-- Compute the largest r^2 that fits in the frame */
num = cd11*cd22-cd12*cd21; num = cd11*cd22-cd12*cd21;
...@@ -2333,7 +2349,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2333,7 +2349,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
x10 = -x1cout - dx1; x10 = -x1cout - dx1;
x2 = -x2cout - dx2; x2 = -x2cout - dx2;
pixin = profit->pmodpix; pixin = profit->pmodpix;
for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) for (ix2=nx2; ix2--; x2+=1.0)
{ {
x1 = x10; x1 = x10;
for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0) for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
...@@ -2374,6 +2390,13 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2374,6 +2390,13 @@ float prof_add(profstruct *prof, profitstruct *profit)
} }
} }
} }
/*---- Copy the symmetric part */
if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
{
pixin2 = pixin - profit->modnaxisn[0] - 1;
for (i=npix2; i--;)
*(pixin++) = *(pixin2--);
}
prof->lostfluxfrac = 1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn)); prof->lostfluxfrac = 1.0 - prof_gammainc(2.0*n, bn*pow(r2max, hinvn));
threshflag = 0; threshflag = 0;
break; break;
...@@ -2382,7 +2405,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2382,7 +2405,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
x10 = -x1cout - dx1; x10 = -x1cout - dx1;
x2 = -x2cout - dx2; x2 = -x2cout - dx2;
pixin = profit->pmodpix; pixin = profit->pmodpix;
for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) for (ix2=nx2; ix2--; x2+=1.0)
{ {
x1 = x10; x1 = x10;
for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0) for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
...@@ -2423,6 +2446,13 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2423,6 +2446,13 @@ float prof_add(profstruct *prof, profitstruct *profit)
} }
} }
} }
/*---- Copy the symmetric part */
if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
{
pixin2 = pixin - profit->modnaxisn[0] - 1;
for (i=npix2; i--;)
*(pixin++) = *(pixin2--);
}
prof->lostfluxfrac = 1.0-prof_gammainc(8.0, 7.66924944*pow(r2max, 0.125)); prof->lostfluxfrac = 1.0-prof_gammainc(8.0, 7.66924944*pow(r2max, 0.125));
threshflag = 0; threshflag = 0;
break; break;
...@@ -2430,7 +2460,7 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2430,7 +2460,7 @@ float prof_add(profstruct *prof, profitstruct *profit)
x10 = -x1cout - dx1; x10 = -x1cout - dx1;
x2 = -x2cout - dx2; x2 = -x2cout - dx2;
pixin = profit->pmodpix; pixin = profit->pmodpix;
for (ix2=profit->modnaxisn[1]; ix2--; x2+=1.0) for (ix2=nx2; ix2--; x2+=1.0)
{ {
x1 = x10; x1 = x10;
for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0) for (ix1=profit->modnaxisn[0]; ix1--; x1+=1.0)
...@@ -2471,6 +2501,13 @@ float prof_add(profstruct *prof, profitstruct *profit) ...@@ -2471,6 +2501,13 @@ float prof_add(profstruct *prof, profitstruct *profit)
} }
} }
} }
/*---- Copy the symmetric part */
if ((npix2=(profit->modnaxisn[1]-nx2)*profit->modnaxisn[0]) > 0)
{
pixin2 = pixin - profit->modnaxisn[0] - 1;
for (i=npix2; i--;)
*(pixin++) = *(pixin2--);
}
rmax = sqrt(r2max); rmax = sqrt(r2max);
prof->lostfluxfrac = (1.0 + rmax)*exp(-rmax); prof->lostfluxfrac = (1.0 + rmax)*exp(-rmax);
threshflag = 0; threshflag = 0;
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Include file for profit.c. * Contents: Include file for profit.c.
* *
* Last modify: 21/09/2009 * Last modify: 24/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
#define PROFIT_MAXITER 1000 /* Max. nb of iterations in profile fitting */ #define PROFIT_MAXITER 1000 /* Max. nb of iterations in profile fitting */
#define PROFIT_MAXPROF 8 /* Max. nb of profile components */ #define PROFIT_MAXPROF 8 /* Max. nb of profile components */
#define PROFIT_OVERSAMP 15 /* Max. profile oversamp. factor on each axis */ #define PROFIT_OVERSAMP 5 /* Max. profile oversamp. factor on each axis */
#define PROFIT_HIDEFRES 201 /* Resolution of the high def. model raster */ #define PROFIT_HIDEFRES 201 /* Resolution of the high def. model raster */
#define PROFIT_REFFFAC 6.0 /* Factor in r_eff for measurement radius*/ #define PROFIT_REFFFAC 6.0 /* Factor in r_eff for measurement radius*/
#define PROFIT_DYNPARAM 10.0 /* Dynamic compression param. in sigma units */ #define PROFIT_DYNPARAM 10.0 /* Dynamic compression param. in sigma units */
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: functions for input of image data. * Contents: functions for input of image data.
* *
* Last modify: 11/10/2007 * Last modify: 29/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -245,6 +245,8 @@ void readimagehead(picstruct *field) ...@@ -245,6 +245,8 @@ void readimagehead(picstruct *field)
field->width = tab->naxisn[0]; field->width = tab->naxisn[0];
field->height = tab->naxisn[1]; field->height = tab->naxisn[1];
field->npix = (KINGSIZE_T)field->width*field->height; field->npix = (KINGSIZE_T)field->width*field->height;
field->bitpix = tab->bitpix;
field->bytepix = tab->bytepix;
if (tab->bitsgn && prefs.fitsunsigned_flag) if (tab->bitsgn && prefs.fitsunsigned_flag)
tab->bitsgn = 0; tab->bitsgn = 0;
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: global type definitions. * Contents: global type definitions.
* *
* Last modify: 24/09/2009 * Last modify: 29/09/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -361,8 +361,9 @@ typedef struct ...@@ -361,8 +361,9 @@ typedef struct
prof_cxyw; /* WORLD ellipse parameters */ prof_cxyw; /* WORLD ellipse parameters */
float prof_pol1w, prof_pol2w; /* WORLD polarisation vector*/ float prof_pol1w, prof_pol2w; /* WORLD polarisation vector*/
float prof_e1w, prof_e2w; /* WORLD ellipticity vector*/ float prof_e1w, prof_e2w; /* WORLD ellipticity vector*/
float prof_class_star; /* Model-fitting star/gal class*/ float prof_class_star; /* Mod.-fitting star/gal class*/
float prof_concentration; /* Model-fitting concentration*/ float prof_concentration; /* Model-fitting concentration*/
float prof_concentrationerr; /* RMS error */
float prof_offset_flux; /* Background offset */ float prof_offset_flux; /* Background offset */
float prof_offset_fluxerr; /* RMS error */ float prof_offset_fluxerr; /* RMS error */
float prof_spheroid_flux; /* Spheroid total flux */ float prof_spheroid_flux; /* Spheroid total flux */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment