From 952f32754f23e78436702c52c14fe26728317ee4 Mon Sep 17 00:00:00 2001 From: Emmanuel Bertin Date: Mon, 20 Jul 2009 19:41:29 +0000 Subject: [PATCH] Added FOCAL coordinates for centering parameters. Re-organized astrometric calls. --- man/sex.1 | 2 +- src/analyse.c | 73 +++++++++++++++++-- src/astrom.c | 180 ++++++++++++++++++++++++++++++---------------- src/astrom.h | 7 +- src/catout.c | 15 ++-- src/param.h | 21 +++++- src/paramprofit.h | 9 ++- src/profit.c | 90 ++++++++++++++++++----- src/profit.h | 3 +- src/types.h | 15 ++-- 10 files changed, 316 insertions(+), 99 deletions(-) diff --git a/man/sex.1 b/man/sex.1 index 05dd7dc..34db77a 100644 --- a/man/sex.1 +++ b/man/sex.1 @@ -1,4 +1,4 @@ -.TH SEXTRACTOR "1" "June 2009" "SWarp 2.8.8" "User Commands" +.TH SEXTRACTOR "1" "July 2009" "SWarp 2.8.8" "User Commands" .SH NAME sex \- extract a source catalog from an astronomical FITS image .SH SYNOPSIS diff --git a/src/analyse.c b/src/analyse.c index d982c43..3085d7d 100644 --- a/src/analyse.c +++ b/src/analyse.c @@ -9,7 +9,7 @@ * * Contents: analyse(), endobject()...: measurements on detections. * -* Last modify: 18/11/2008 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -385,9 +385,11 @@ Final processing of object data, just before saving it to the catalog. void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, picstruct *dwfield, int n, objliststruct *objlist) { + objstruct *obj; checkstruct *check; + double rawpos[NAXIS], + pixscale2; int i,j, ix,iy,selecflag, newnumber,nsub; - objstruct *obj; obj = &objlist->obj[n]; @@ -468,6 +470,43 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, obj2->polar = (obj->a*obj->a - obj->b*obj->b) / (obj->a*obj->a + obj->b*obj->b); +/*-- Express positions in FOCAL or WORLD coordinates */ + if (FLAG(obj2.mxf) || FLAG(obj2.mxw)) + astrom_pos(field, obj); + + if (FLAG(obj2.mx2w) + || FLAG(obj2.win_mx2w) + || FLAG(obj2.poserr_mx2w) + || FLAG(obj2.winposerr_mx2w) + || FLAG(obj2.poserrmx2w_prof) + || FLAG(obj2.prof_flagw) + || ((!prefs.pixel_scale) && (FLAG(obj2.npixw) + || FLAG(obj2.fdnpixw) + || FLAG(obj2.fwhmw)))) + { + rawpos[0] = obj2->posx; + rawpos[1] = obj2->posy; + pixscale2 = wcs_jacobian(field->wcs, rawpos, obj2->jacob); + } + +/*-- Express shape parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.mx2w)) + astrom_shapeparam(field, obj); +/*-- Express position error parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.poserr_mx2w)) + astrom_errparam(field, obj); + + if (FLAG(obj2.npixw)) + obj2->npixw = obj->npix * (prefs.pixel_scale? + field->pixscale/3600.0*field->pixscale/3600.0 : pixscale2); + if (FLAG(obj2.fdnpixw)) + obj2->fdnpixw = obj->fdnpix * (prefs.pixel_scale? + field->pixscale/3600.0*field->pixscale/3600.0 : pixscale2); + + if (FLAG(obj2.fwhmw)) + obj2->fwhmw = obj->fwhm * (prefs.pixel_scale? + field->pixscale/3600.0 : sqrt(pixscale2)); + /*------------------------------- Photometry -------------------------------*/ /*-- Convert the father of photom. error estimates from variance to RMS */ @@ -493,11 +532,25 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, /*--------------------------- Windowed barycenter --------------------------*/ if (FLAG(obj2.winpos_x)) + { compute_winpos(field, wfield, obj); +/*---- Express positions in FOCAL or WORLD coordinates */ + if (FLAG(obj2.winpos_xf) || FLAG(obj2.winpos_xw)) + astrom_winpos(field, obj); +/*---- Express shape parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.mx2w)) + astrom_winshapeparam(field, obj); +/*---- Express position error parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.poserr_mx2w)) + astrom_winerrparam(field, obj); + } -/*-- What about the peak of the profile? */ +/*---------------------------- Peak information ----------------------------*/ if (obj->peak+obj->bkg >= field->satur_level) obj->flag |= OBJ_SATUR; +/*-- Express positions in FOCAL or WORLD coordinates */ + if (FLAG(obj2.peakxf) || FLAG(obj2.peakxw)) + astrom_peakpos(field, obj); /*-- Check-image CHECK_APERTURES option */ @@ -639,14 +692,24 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, /*----------------------------- Profile fitting -----------------------------*/ nsub = 1; if (prefs.prof_flag) + { profit_fit(theprofit, field, wfield, obj, obj2); +/*---- Express positions in FOCAL or WORLD coordinates */ + if (FLAG(obj2.winpos_xf) || FLAG(obj2.winpos_xw)) + astrom_profpos(field, obj); +/*---- Express shape parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.mx2w)) + astrom_profshapeparam(field, obj); +/*---- Express position error parameters in the FOCAL or WORLD frame */ + if (FLAG(obj2.poserr_mx2w)) + astrom_proferrparam(field, obj); + } /*--- Express everything in magnitude units */ computemags(field, obj); /*-------------------------------- Astrometry ------------------------------*/ - if (prefs.world_flag) - computeastrom(field, obj); + /*-- Edit min and max coordinates to follow the FITS conventions */ obj->xmin += 1; obj->ymin += 1; diff --git a/src/astrom.c b/src/astrom.c index 86ff9c6..c5c9663 100644 --- a/src/astrom.c +++ b/src/astrom.c @@ -9,7 +9,7 @@ * * Contents: Astrometrical computations. * -* Last modify: 19/05/2008 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -79,11 +79,11 @@ void initastrom(picstruct *field) } -/**************************** computeastrom *********************************/ +/***************************** astrom_pos **********************************/ /* -Compute real WORLD coordinates and dimensions according to FITS info. +Compute real FOCAL and WORLD coordinates according to FITS info. */ -void computeastrom(picstruct *field, objstruct *obj) +void astrom_pos(picstruct *field, objstruct *obj) { wcsstruct *wcs; @@ -94,7 +94,16 @@ void computeastrom(picstruct *field, objstruct *obj) wcs = field->wcs; lng = wcs->lng; lat = wcs->lat; - pixscale2 = 0.0; /* To avoid gcc -Wall warnings */ + +/* If working with WCS, compute FOCAL coordinates and local matrix */ + if (FLAG(obj2.mxf)) + { + rawpos[0] = obj2->posx; + rawpos[1] = obj2->posy; + raw_to_red(wcs, rawpos, wcspos); + obj2->mxf = wcspos[0]; + obj2->myf = wcspos[1]; + } /* If working with WCS, compute WORLD coordinates and local matrix */ if (FLAG(obj2.mxw)) @@ -145,7 +154,58 @@ void computeastrom(picstruct *field, objstruct *obj) } } -/* Idem for peak-flux positions */ +/* Custom coordinate system for the MAMA machine */ + if (FLAG(obj2.mamaposx)) + { + rawpos[0] = obj2->posx - 0.5; + rawpos[1] = obj2->posy - 0.5; + raw_to_wcs(wcs, rawpos, wcspos); + obj2->mamaposx = wcspos[1]*(MAMA_CORFLEX+1.0); + obj2->mamaposy = wcspos[0]*(MAMA_CORFLEX+1.0); + } + + if (FLAG(obj2.mx2w) + || FLAG(obj2.win_mx2w) + || FLAG(obj2.poserr_mx2w) + || FLAG(obj2.winposerr_mx2w) + || FLAG(obj2.poserrmx2w_prof) + || FLAG(obj2.prof_flagw) + || ((!prefs.pixel_scale) && (FLAG(obj2.npixw) + || FLAG(obj2.fdnpixw) + || FLAG(obj2.fwhmw)))) + { + rawpos[0] = obj2->posx; + rawpos[1] = obj2->posy; + pixscale2 = wcs_jacobian(wcs, rawpos, obj2->jacob); + } + + return; + } + + +/***************************** astrom_peakpos *******************************/ +/* +Compute real FOCAL and WORLD peak coordinates according to FITS info. +*/ +void astrom_peakpos(picstruct *field, objstruct *obj) + + { + wcsstruct *wcs; + double rawpos[NAXIS], wcspos[NAXIS]; + int lng,lat; + + wcs = field->wcs; + lng = wcs->lng; + lat = wcs->lat; + + if (FLAG(obj2.peakxf)) + { + rawpos[0] = obj->peakx; + rawpos[1] = obj->peaky; + raw_to_red(wcs, rawpos, wcspos); + obj2->peakxf = wcspos[0]; + obj2->peakyf = wcspos[1]; + } if (FLAG(obj2.peakxw)) { rawpos[0] = obj->peakx; @@ -174,7 +234,33 @@ void computeastrom(picstruct *field, objstruct *obj) } } -/* Idem for Windowed positions */ + return; + } + + +/****************************** astrom_winpos *******************************/ +/* +Compute real FOCAL and WORLD windowed coordinates according to FITS info. +*/ +void astrom_winpos(picstruct *field, objstruct *obj) + + { + wcsstruct *wcs; + double rawpos[NAXIS], wcspos[NAXIS]; + int lng,lat; + + wcs = field->wcs; + lng = wcs->lng; + lat = wcs->lat; + + if (FLAG(obj2.winpos_xf)) + { + rawpos[0] = obj2->winpos_x; + rawpos[1] = obj2->winpos_y; + raw_to_red(wcs, rawpos, wcspos); + obj2->winpos_xf = wcspos[0]; + obj2->winpos_yf = wcspos[1]; + } if (FLAG(obj2.winpos_xw)) { rawpos[0] = obj2->winpos_x; @@ -203,7 +289,33 @@ void computeastrom(picstruct *field, objstruct *obj) } } -/* Idem for Model-fitted positions */ + return; + } + + +/****************************** astrom_profpos *******************************/ +/* +Compute real FOCAL and WORLD profit coordinates according to FITS info. +*/ +void astrom_profpos(picstruct *field, objstruct *obj) + + { + wcsstruct *wcs; + double rawpos[NAXIS], wcspos[NAXIS]; + int lng,lat; + + wcs = field->wcs; + lng = wcs->lng; + lat = wcs->lat; + + if (FLAG(obj2.xf_prof)) + { + rawpos[0] = obj2->x_prof; + rawpos[1] = obj2->y_prof; + raw_to_red(wcs, rawpos, wcspos); + obj2->xf_prof = wcspos[0]; + obj2->yf_prof = wcspos[1]; + } if (FLAG(obj2.xw_prof)) { rawpos[0] = obj2->x_prof; @@ -232,58 +344,6 @@ void computeastrom(picstruct *field, objstruct *obj) } } -/* Custom coordinate system for the MAMA machine */ - if (FLAG(obj2.mamaposx)) - { - rawpos[0] = obj2->posx - 0.5; - rawpos[1] = obj2->posy - 0.5; - raw_to_wcs(wcs, rawpos, wcspos); - obj2->mamaposx = wcspos[1]*(MAMA_CORFLEX+1.0); - obj2->mamaposy = wcspos[0]*(MAMA_CORFLEX+1.0); - } - - if (FLAG(obj2.mx2w) - || FLAG(obj2.win_mx2w) - || FLAG(obj2.poserr_mx2w) - || FLAG(obj2.winposerr_mx2w) - || FLAG(obj2.poserrmx2w_prof) - || FLAG(obj2.prof_flagw) - || ((!prefs.pixel_scale) && (FLAG(obj2.npixw) - || FLAG(obj2.fdnpixw) - || FLAG(obj2.fwhmw)))) - { - rawpos[0] = obj2->posx; - rawpos[1] = obj2->posy; - pixscale2 = wcs_jacobian(wcs, rawpos, obj2->jacob); - } - -/* Express shape parameters in WORLD frame */ - if (FLAG(obj2.mx2w)) - astrom_shapeparam(field, obj); - if (FLAG(obj2.win_mx2w)) - astrom_winshapeparam(field, obj); - if (FLAG(obj2.prof_flagw)) - astrom_profshapeparam(field, obj); - -/* Express position error parameters in WORLD frame */ - if (FLAG(obj2.poserr_mx2w)) - astrom_errparam(field, obj); - if (FLAG(obj2.winposerr_mx2w)) - astrom_winerrparam(field, obj); - if (FLAG(obj2.poserrmx2w_prof)) - astrom_proferrparam(field, obj); - - if (FLAG(obj2.npixw)) - obj2->npixw = obj->npix * (prefs.pixel_scale? - field->pixscale/3600.0*field->pixscale/3600.0 : pixscale2); - if (FLAG(obj2.fdnpixw)) - obj2->fdnpixw = obj->fdnpix * (prefs.pixel_scale? - field->pixscale/3600.0*field->pixscale/3600.0 : pixscale2); - - if (FLAG(obj2.fwhmw)) - obj2->fwhmw = obj->fwhm * (prefs.pixel_scale? - field->pixscale/3600.0 : sqrt(pixscale2)); - return; } diff --git a/src/astrom.h b/src/astrom.h index 06f3039..871be94 100644 --- a/src/astrom.h +++ b/src/astrom.h @@ -9,7 +9,7 @@ * * Contents: Astrometrical stuff. * -* Last modify: 18/05/2008 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -30,12 +30,15 @@ /*------------------------------- structures --------------------------------*/ /*------------------------------- functions ---------------------------------*/ extern void astrom_errparam(picstruct *, objstruct *), + astrom_peakpos(picstruct *, objstruct *), + astrom_pos(picstruct *, objstruct *), astrom_proferrparam(picstruct *, objstruct *), + astrom_profpos(picstruct *, objstruct *), astrom_profshapeparam(picstruct *, objstruct *), astrom_shapeparam(picstruct *, objstruct *), astrom_winerrparam(picstruct *, objstruct *), + astrom_winpos(picstruct *, objstruct *), astrom_winshapeparam(picstruct *, objstruct *), - computeastrom(picstruct *, objstruct *), initastrom(picstruct *), j2b(double, double, double, double *, double *), precess(double,double,double,double,double *,double *); diff --git a/src/catout.c b/src/catout.c index f7afd78..267ec50 100644 --- a/src/catout.c +++ b/src/catout.c @@ -9,7 +9,7 @@ * * Contents: functions for output of catalog data. * -* Last modify: 29/05/2009 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -201,13 +201,15 @@ void updateparamflags() | FLAG(obj2.alpha2000_prof); FLAG(obj2.xw_prof) |= FLAG(obj2.yw_prof) | FLAG(obj2.alphas_prof); - + FLAG(obj2.xf_prof) |= FLAG(obj2.yf_prof); FLAG(obj2.x_prof) |= FLAG(obj2.y_prof) | FLAG(obj2.xw_prof) + | FLAG(obj2.xf_prof) | FLAG(obj2.poserra_prof) | FLAG(obj2.poserrcxx_prof) | FLAG(obj2.prof_concentration) | FLAG(obj2.prof_class_star); + FLAG(obj2.mag_prof) |= FLAG(obj2.magerr_prof); FLAG(obj2.magerr_prof) |= FLAG(obj2.fluxerr_prof); FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof) | FLAG(obj2.fluxerr_prof); @@ -392,6 +394,7 @@ void updateparamflags() | FLAG(obj2.winposerr_theta2000); FLAG(obj2.winpos_alphas) |= FLAG(obj2.winpos_deltas) | FLAG(obj2.winpos_alpha2000); + FLAG(obj2.winpos_xf) |= FLAG(obj2.winpos_yf); FLAG(obj2.winpos_xw) |= FLAG(obj2.winpos_yw) | FLAG(obj2.winpos_alphas); @@ -432,8 +435,11 @@ void updateparamflags() | FLAG(obj2.npixw) | FLAG(obj2.fdnpixw) | FLAG(obj2.fwhmw); + FLAG(obj2.peakxw) |= FLAG(obj2.peakyf); FLAG(obj2.peakxw) |= FLAG(obj2.peakyw) | FLAG(obj2.peakalphas); - FLAG(obj.peakx) |= FLAG(obj.peaky) | FLAG(obj2.peakxw); + FLAG(obj.peakx) |= FLAG(obj.peaky) | FLAG(obj2.peakxw) | FLAG(obj2.peakxf); + + FLAG(obj2.mxf) |= FLAG(obj2.myf); FLAG(obj2.mxw) |= FLAG(obj2.myw) | FLAG(obj2.mx2w) | FLAG(obj2.alphas) | FLAG(obj2.poserr_mx2w); @@ -442,7 +448,8 @@ void updateparamflags() | FLAG(obj2.flux_win) | FLAG(obj2.fluxerr_win); FLAG(obj2.winpos_x) |= FLAG(obj2.winpos_y) | FLAG(obj2.winposerr_mx2) | FLAG(obj2.win_mx2) - | FLAG(obj2.winpos_xw) | FLAG(obj2.win_flag) + | FLAG(obj2.winpos_xw) | FLAG(obj2.winpos_xf) + | FLAG(obj2.win_flag) | FLAG(obj2.flux_win) |FLAG(obj2.winpos_niter); /*------------------------------ Photometry ---------------------------------*/ diff --git a/src/param.h b/src/param.h index f2ac9a4..330646f 100644 --- a/src/param.h +++ b/src/param.h @@ -9,7 +9,7 @@ * * Contents: parameter list for catalog data. * -* Last modify: 19/05/2009 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -222,6 +222,12 @@ keystruct objkey[] = { {"YPEAK_IMAGE", "y-coordinate of the brightest pixel", &outobj.peaky, H_INT, T_LONG, "%10d", "pixel", "pos.cartesian.y", "pix"}, + {"XPEAK_FOCAL", "Focal-plane x coordinate of the brightest pixel", + &outobj2.peakxf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.x", ""}, + {"YPEAK_FOCAL", "Focal-plane y coordinate of the brightest pixel", + &outobj2.peakyf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.y", ""}, {"XPEAK_WORLD", "World-x coordinate of the brightest pixel", &outobj2.peakxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", "pos.eq.ra", "deg"}, @@ -262,6 +268,12 @@ keystruct objkey[] = { {"Y_IMAGE_DBL", "Object position along y (double precision)", &outobj2.posy, H_FLOAT, T_DOUBLE, "%10.3f", "pixel", "pos.cartesian.x;pos.barycenter;instr.det", "pix"}, + {"X_FOCAL", "Barycenter position along focal-plane x axis", + &outobj2.mxf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.x", ""}, + {"Y_FOCAL", "Barycenter position along focal-plane y axis", + &outobj2.myf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.y", ""}, {"X_WORLD", "Barycenter position along world x axis", &outobj2.mxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", "pos.eq.ra", "deg"}, @@ -435,6 +447,13 @@ keystruct objkey[] = { &outobj2.winpos_y, H_FLOAT, T_DOUBLE, "%10.3f", "pixel", "pos.cartesian.y;instr.det", "pix"}, + {"XWIN_FOCAL", "Windowed position along focal-plane x axis", + &outobj2.winpos_xf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.x", ""}, + {"YWIN_FOCAL", "Windowed position along focal-plane y axis", + &outobj2.winpos_yf, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.y", ""}, + {"XWIN_WORLD", "Windowed position along world x axis", &outobj2.winpos_xw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", "pos.eq.ra", "deg"}, diff --git a/src/paramprofit.h b/src/paramprofit.h index 46ae118..2413d99 100644 --- a/src/paramprofit.h +++ b/src/paramprofit.h @@ -9,7 +9,7 @@ * * Contents: Model-fitting parameter list for catalog data. * -* Last modify: 27/05/2009 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -50,6 +50,13 @@ &outobj2.y_prof, H_FLOAT, T_FLOAT, "%10.3f", "pixel", "pos.cartesian.y;stat.fit.param;instr.det;meta.main", "pix"}, + {"XMODEL_WORLD", "Fitted position along focal-plane x axis", + &outobj2.xf_prof, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.x;stat.fit.param", ""}, + {"YMODEL_WORLD", "Fitted position along focal-plane y axis", + &outobj2.yf_prof, H_FLOAT, T_DOUBLE, "%18.10e", "", + "pos.cartesian.y;stat.fit.param", ""}, + {"XMODEL_WORLD", "Fitted position along world x axis", &outobj2.xw_prof, H_FLOAT, T_DOUBLE, "%18.10e", "deg", "pos.eq.ra;stat.fit.param", "deg"}, diff --git a/src/profit.c b/src/profit.c index daeed35..c076b53 100644 --- a/src/profit.c +++ b/src/profit.c @@ -9,7 +9,7 @@ * * Contents: Fit an arbitrary profile combination to a detection. * -* Last modify: 29/05/2009 +* Last modify: 13/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -168,7 +168,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the fit). NOTES It is a modified version of the lm_minimize() of lmfit. AUTHOR E. Bertin (IAP) -VERSION 29/05/2009 +VERSION 13/07/2009 ***/ void profit_fit(profitstruct *profit, picstruct *field, picstruct *wfield, @@ -178,8 +178,10 @@ void profit_fit(profitstruct *profit, patternstruct *pattern; psfstruct *psf; checkstruct *check; - double psf_fwhm, a , cp,sp, emx2,emy2,emxy, dchi2, err; - int i,j,p, nparam, ncomp, nprof; + double *oldparaminit, + psf_fwhm, a , cp,sp, emx2,emy2,emxy, dchi2, err, + oldchi2; + int i,j,p, nparam, ncomp, nprof, oldniter; nparam = profit->nparam; nprof = profit->nprof; @@ -257,7 +259,21 @@ the_gal++; /* Actual minimisation */ profit->niter = profit_minimize(profit, PROFIT_MAXITER); profit_residuals(profit,field,wfield, 10.0, profit->param,profit->resi); - +/* + oldchi2 = profit->chi2; + oldniter = profit->niter; + oldparaminit = profit_reresetparams(profit); + profit->niter = profit_minimize(profit, PROFIT_MAXITER); + if (profit->chi2 > oldchi2) + { + memcpy(profit->paraminit, oldparaminit, nparam*sizeof(double)); + profit->chi2 = oldchi2; + profit->niter = oldniter; + } + else + obj2->prof_flag |= PROFIT_FLIPPED; + free(oldparaminit); +*/ /* QMEMCPY(profit->paraminit, oldparaminit, double, nparam); if (profit_setparam(profit, PARAM_ARMS_PITCH, 160.0, 130.0, 175.0)==RETURN_OK) @@ -593,7 +609,7 @@ INPUT Profile-fitting structure. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 22/10/2008 +VERSION 30/06/2009 ***/ void profit_psf(profitstruct *profit) { @@ -644,6 +660,7 @@ void profit_psf(profitstruct *profit) *(pixout++) *= norm; } + return; } @@ -868,8 +885,8 @@ double *profit_compresi(profitstruct *profit, double dynparam, double *resi) val = *(objpix++); if ((wval=*(objweight++))>0.0) { - val2 = (val - *lmodpix)*wval*invsig; - val2 = val2>0.0? LOGF(1.0+val2) : -LOGF(1.0-val2); + val2 = (val - *lmodpix)*invsig; + // val2 = val2>0.0? sqrt(1.0+val2) - 1.0 : 1.0 - sqrt(1.0-val2); *(resit++) = val2*dynparam; error += val2*val2; } @@ -1341,7 +1358,7 @@ INPUT Profile-fitting structure. OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 27/05/2009 +VERSION 29/06/2009 ***/ void profit_moments(profitstruct *profit) { @@ -1349,13 +1366,14 @@ void profit_moments(profitstruct *profit) obj2struct *obj2; double *pix, hw,hh, x,y, xstart, val, - mx,my, sum, mx2,my2,mxy, den; + mx,my, sum, mx2,my2,mxy, den, r2max; int ix,iy; obj = profit->obj; obj2 = profit->obj2; hw = (double)(profit->modnaxisn[0]/2); hh = (double)(profit->modnaxisn[1]/2); + r2max = hwmodpix; @@ -1364,15 +1382,18 @@ void profit_moments(profitstruct *profit) { x = xstart; for (ix=profit->modnaxisn[0]; ix--; x+=1.0) - { - val = *(pix++); - sum += val; - mx += val*x; - my += val*y; - mx2 += val*x*x; - mxy += val*x*y; - my2 += val*y*y; - } + if (y*y+x*x <= r2max) + { + val = *(pix++); + sum += val; + mx += val*x; + my += val*y; + mx2 += val*x*x; + mxy += val*x*y; + my2 += val*y*y; + } + else + pix++; } if (sum <= 1.0/BIG) @@ -1450,6 +1471,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) obj = profit->obj; obj2 = profit->obj2; param = parammin = parammax = 0.0; /* Avoid gcc -Wall warnings*/ + switch(paramtype) { case PARAM_BACK: @@ -1619,6 +1641,35 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) } +/****** profit_reresetparams **************************************************** +PROTO double *profit_reresetparams(profitstruct *profit) +PURPOSE Reset profile parameters once more according to previous findings. +INPUT Pointer to the profit structure. +OUTPUT -. +NOTES -. +AUTHOR E. Bertin (IAP) +VERSION 13/07/2009 + ***/ +double *profit_reresetparams(profitstruct *profit) + { + double *param, *parami; + int p, nparam; + + QMEMCPY(profit->paraminit, param, double, PARAM_NPARAM); + nparam = profit->nparam; + profit_resetparams(profit); + profit_boundtounbound(profit, profit->paraminit); + profit_boundtounbound(profit, param); + parami = profit->paraminit; + for (p=0; pparaminit); + profit_unboundtobound(profit, param); + + return param; + } + + /****** profit_resetparams **************************************************** PROTO void profit_resetparams(profitstruct *profit) PURPOSE Set the initial, lower and upper boundary values of profile parameters. @@ -1640,6 +1691,7 @@ void profit_resetparams(profitstruct *profit) } + /****** profit_setparam **************************************************** PROTO void profit_setparam(profitstruct *profit, paramenum paramtype, double param, double parammin, double parammax) diff --git a/src/profit.h b/src/profit.h index 325f2cd..bc31ffe 100644 --- a/src/profit.h +++ b/src/profit.h @@ -9,7 +9,7 @@ * * Contents: Include file for profit.c. * -* Last modify: 20/03/2009 +* Last modify: 13/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -144,6 +144,7 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode); double *profit_compresi(profitstruct *profit, double dynparam, double *resi), + *profit_reresetparams(profitstruct *profit), *profit_residuals(profitstruct *profit, picstruct *field, picstruct *wfield, double dynparam, double *param, double *resi), diff --git a/src/types.h b/src/types.h index 682461c..cf47b46 100644 --- a/src/types.h +++ b/src/types.h @@ -9,7 +9,7 @@ * * Contents: global type definitions. * -* Last modify: 18/03/2009 +* Last modify: 20/07/2009 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -153,6 +153,7 @@ typedef struct /* ---- astrometric data */ double posx,posy; /* "FITS" pos. in pixels */ double jacob[NAXIS*NAXIS]; /* Local deproject. Jacobian */ + double mamaposx,mamaposy; /* "MAMA" pos. in pixels */ float sposx,sposy; /* single precision pos. */ float poserr_a, poserr_b, @@ -170,7 +171,9 @@ typedef struct poserr_cxyw; /* WORLD error ellipse */ double mx2w,my2w,mxyw; /* WORLD var. and covar. */ double peakxw, peakyw; /* WORLD of brightest pix */ + double peakxf, peakyf; /* FOCAL of brightest pix */ double mxw, myw; /* WORLD barycenters */ + double mxf, myf; /* FOCAL barycenters */ double alphas, deltas; /* native alpha, delta */ float thetas; /* native position angle E/N*/ double peakalphas, peakdeltas; /* native for brightest pix */ @@ -228,12 +231,13 @@ typedef struct float win_aw, win_bw, win_thetaw; /* WORLD ellipse parameters */ float win_polarw; /* WORLD WIN "polarization" */ - float win_thetas; /* native error pos. angle */ - float win_theta2000; /* J2000 error pos. angle */ - float win_theta1950; /* B1950 error pos. angle */ + float win_thetas; /* native error pos. angle */ + float win_theta2000; /* J2000 error pos. angle */ + float win_theta1950; /* B1950 error pos. angle */ float win_cxxw, win_cyyw, win_cxyw; /* WORLD ellipse parameters */ double winpos_xw, winpos_yw; /* WORLD coordinates */ + double winpos_xf, winpos_yf; /* FOCAL coordinates */ double winpos_alphas, winpos_deltas; /* native alpha, delta */ double winpos_alpha2000, winpos_delta2000; /* J2000 alpha, delta */ double winpos_alpha1950, winpos_delta1950; /* B1950 alpha, delta */ @@ -311,7 +315,8 @@ typedef struct float mag_prof; /* Mag from model-fitting */ float magerr_prof; /* RMS mag from model-fitting */ float x_prof, y_prof; /* Coords from model-fitting*/ - double xw_prof, yw_prof; /* WORLD coords */ + double xw_prof, yw_prof; /* WORLD coordinates */ + double xf_prof, yf_prof; /* FOCAL coordinates */ double alphas_prof, deltas_prof; /* native alpha, delta */ double alpha2000_prof, delta2000_prof; /* J2000 alpha, delta */ double alpha1950_prof, delta1950_prof; /* B1950 alpha, delta */ -- GitLab