Commit 952f3275 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Added FOCAL coordinates for centering parameters.

Re-organized astrometric calls.
parent f823a212
.TH SEXTRACTOR "1" "June 2009" "SWarp 2.8.8" "User Commands" .TH SEXTRACTOR "1" "July 2009" "SWarp 2.8.8" "User Commands"
.SH NAME .SH NAME
sex \- extract a source catalog from an astronomical FITS image sex \- extract a source catalog from an astronomical FITS image
.SH SYNOPSIS .SH SYNOPSIS
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: analyse(), endobject()...: measurements on detections. * 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. ...@@ -385,9 +385,11 @@ Final processing of object data, just before saving it to the catalog.
void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
picstruct *dwfield, int n, objliststruct *objlist) picstruct *dwfield, int n, objliststruct *objlist)
{ {
objstruct *obj;
checkstruct *check; checkstruct *check;
double rawpos[NAXIS],
pixscale2;
int i,j, ix,iy,selecflag, newnumber,nsub; int i,j, ix,iy,selecflag, newnumber,nsub;
objstruct *obj;
obj = &objlist->obj[n]; obj = &objlist->obj[n];
...@@ -468,6 +470,43 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, ...@@ -468,6 +470,43 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
obj2->polar = (obj->a*obj->a - obj->b*obj->b) obj2->polar = (obj->a*obj->a - obj->b*obj->b)
/ (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 -------------------------------*/ /*------------------------------- Photometry -------------------------------*/
/*-- Convert the father of photom. error estimates from variance to RMS */ /*-- Convert the father of photom. error estimates from variance to RMS */
...@@ -493,11 +532,25 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, ...@@ -493,11 +532,25 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
/*--------------------------- Windowed barycenter --------------------------*/ /*--------------------------- Windowed barycenter --------------------------*/
if (FLAG(obj2.winpos_x)) if (FLAG(obj2.winpos_x))
{
compute_winpos(field, wfield, obj); 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) if (obj->peak+obj->bkg >= field->satur_level)
obj->flag |= OBJ_SATUR; 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 */ /*-- Check-image CHECK_APERTURES option */
...@@ -639,14 +692,24 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield, ...@@ -639,14 +692,24 @@ void endobject(picstruct *field, picstruct *dfield, picstruct *wfield,
/*----------------------------- Profile fitting -----------------------------*/ /*----------------------------- Profile fitting -----------------------------*/
nsub = 1; nsub = 1;
if (prefs.prof_flag) if (prefs.prof_flag)
{
profit_fit(theprofit, field, wfield, obj, obj2); 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 */ /*--- Express everything in magnitude units */
computemags(field, obj); computemags(field, obj);
/*-------------------------------- Astrometry ------------------------------*/ /*-------------------------------- Astrometry ------------------------------*/
if (prefs.world_flag)
computeastrom(field, obj);
/*-- Edit min and max coordinates to follow the FITS conventions */ /*-- Edit min and max coordinates to follow the FITS conventions */
obj->xmin += 1; obj->xmin += 1;
obj->ymin += 1; obj->ymin += 1;
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Astrometrical computations. * Contents: Astrometrical computations.
* *
* Last modify: 19/05/2008 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -79,11 +79,11 @@ void initastrom(picstruct *field) ...@@ -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; wcsstruct *wcs;
...@@ -94,7 +94,16 @@ void computeastrom(picstruct *field, objstruct *obj) ...@@ -94,7 +94,16 @@ void computeastrom(picstruct *field, objstruct *obj)
wcs = field->wcs; wcs = field->wcs;
lng = wcs->lng; lng = wcs->lng;
lat = wcs->lat; 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 working with WCS, compute WORLD coordinates and local matrix */
if (FLAG(obj2.mxw)) if (FLAG(obj2.mxw))
...@@ -145,7 +154,58 @@ void computeastrom(picstruct *field, objstruct *obj) ...@@ -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)) if (FLAG(obj2.peakxw))
{ {
rawpos[0] = obj->peakx; rawpos[0] = obj->peakx;
...@@ -174,7 +234,33 @@ void computeastrom(picstruct *field, objstruct *obj) ...@@ -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)) if (FLAG(obj2.winpos_xw))
{ {
rawpos[0] = obj2->winpos_x; rawpos[0] = obj2->winpos_x;
...@@ -203,7 +289,33 @@ void computeastrom(picstruct *field, objstruct *obj) ...@@ -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)) if (FLAG(obj2.xw_prof))
{ {
rawpos[0] = obj2->x_prof; rawpos[0] = obj2->x_prof;
...@@ -232,58 +344,6 @@ void computeastrom(picstruct *field, objstruct *obj) ...@@ -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; return;
} }
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Astrometrical stuff. * Contents: Astrometrical stuff.
* *
* Last modify: 18/05/2008 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -30,12 +30,15 @@ ...@@ -30,12 +30,15 @@
/*------------------------------- structures --------------------------------*/ /*------------------------------- structures --------------------------------*/
/*------------------------------- functions ---------------------------------*/ /*------------------------------- functions ---------------------------------*/
extern void astrom_errparam(picstruct *, objstruct *), extern void astrom_errparam(picstruct *, objstruct *),
astrom_peakpos(picstruct *, objstruct *),
astrom_pos(picstruct *, objstruct *),
astrom_proferrparam(picstruct *, objstruct *), astrom_proferrparam(picstruct *, objstruct *),
astrom_profpos(picstruct *, objstruct *),
astrom_profshapeparam(picstruct *, objstruct *), astrom_profshapeparam(picstruct *, objstruct *),
astrom_shapeparam(picstruct *, objstruct *), astrom_shapeparam(picstruct *, objstruct *),
astrom_winerrparam(picstruct *, objstruct *), astrom_winerrparam(picstruct *, objstruct *),
astrom_winpos(picstruct *, objstruct *),
astrom_winshapeparam(picstruct *, objstruct *), astrom_winshapeparam(picstruct *, objstruct *),
computeastrom(picstruct *, objstruct *),
initastrom(picstruct *), initastrom(picstruct *),
j2b(double, double, double, double *, double *), j2b(double, double, double, double *, double *),
precess(double,double,double,double,double *,double *); precess(double,double,double,double,double *,double *);
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: functions for output of catalog data. * Contents: functions for output of catalog data.
* *
* Last modify: 29/05/2009 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -201,13 +201,15 @@ void updateparamflags() ...@@ -201,13 +201,15 @@ void updateparamflags()
| FLAG(obj2.alpha2000_prof); | FLAG(obj2.alpha2000_prof);
FLAG(obj2.xw_prof) |= FLAG(obj2.yw_prof) FLAG(obj2.xw_prof) |= FLAG(obj2.yw_prof)
| FLAG(obj2.alphas_prof); | FLAG(obj2.alphas_prof);
FLAG(obj2.xf_prof) |= FLAG(obj2.yf_prof);
FLAG(obj2.x_prof) |= FLAG(obj2.y_prof) FLAG(obj2.x_prof) |= FLAG(obj2.y_prof)
| FLAG(obj2.xw_prof) | FLAG(obj2.xw_prof)
| FLAG(obj2.xf_prof)
| FLAG(obj2.poserra_prof) | FLAG(obj2.poserra_prof)
| FLAG(obj2.poserrcxx_prof) | FLAG(obj2.poserrcxx_prof)
| 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.mag_prof) |= FLAG(obj2.magerr_prof);
FLAG(obj2.magerr_prof) |= FLAG(obj2.fluxerr_prof); FLAG(obj2.magerr_prof) |= FLAG(obj2.fluxerr_prof);
FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof) | FLAG(obj2.fluxerr_prof); FLAG(obj2.flux_prof) |= FLAG(obj2.mag_prof) | FLAG(obj2.fluxerr_prof);
...@@ -392,6 +394,7 @@ void updateparamflags() ...@@ -392,6 +394,7 @@ void updateparamflags()
| FLAG(obj2.winposerr_theta2000); | FLAG(obj2.winposerr_theta2000);
FLAG(obj2.winpos_alphas) |= FLAG(obj2.winpos_deltas) FLAG(obj2.winpos_alphas) |= FLAG(obj2.winpos_deltas)
| FLAG(obj2.winpos_alpha2000); | FLAG(obj2.winpos_alpha2000);
FLAG(obj2.winpos_xf) |= FLAG(obj2.winpos_yf);
FLAG(obj2.winpos_xw) |= FLAG(obj2.winpos_yw) FLAG(obj2.winpos_xw) |= FLAG(obj2.winpos_yw)
| FLAG(obj2.winpos_alphas); | FLAG(obj2.winpos_alphas);
...@@ -432,8 +435,11 @@ void updateparamflags() ...@@ -432,8 +435,11 @@ void updateparamflags()
| FLAG(obj2.npixw) | FLAG(obj2.fdnpixw) | FLAG(obj2.npixw) | FLAG(obj2.fdnpixw)
| FLAG(obj2.fwhmw); | FLAG(obj2.fwhmw);
FLAG(obj2.peakxw) |= FLAG(obj2.peakyf);
FLAG(obj2.peakxw) |= FLAG(obj2.peakyw) | FLAG(obj2.peakalphas); 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.mxw) |= FLAG(obj2.myw) | FLAG(obj2.mx2w) | FLAG(obj2.alphas)
| FLAG(obj2.poserr_mx2w); | FLAG(obj2.poserr_mx2w);
...@@ -442,7 +448,8 @@ void updateparamflags() ...@@ -442,7 +448,8 @@ void updateparamflags()
| FLAG(obj2.flux_win) | FLAG(obj2.fluxerr_win); | FLAG(obj2.flux_win) | FLAG(obj2.fluxerr_win);
FLAG(obj2.winpos_x) |= FLAG(obj2.winpos_y) FLAG(obj2.winpos_x) |= FLAG(obj2.winpos_y)
| FLAG(obj2.winposerr_mx2) | FLAG(obj2.win_mx2) | 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); | FLAG(obj2.flux_win) |FLAG(obj2.winpos_niter);
/*------------------------------ Photometry ---------------------------------*/ /*------------------------------ Photometry ---------------------------------*/
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: parameter list for catalog data. * Contents: parameter list for catalog data.
* *
* Last modify: 19/05/2009 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -222,6 +222,12 @@ keystruct objkey[] = { ...@@ -222,6 +222,12 @@ keystruct objkey[] = {
{"YPEAK_IMAGE", "y-coordinate of the brightest pixel", {"YPEAK_IMAGE", "y-coordinate of the brightest pixel",
&outobj.peaky, H_INT, T_LONG, "%10d", "pixel", &outobj.peaky, H_INT, T_LONG, "%10d", "pixel",
"pos.cartesian.y", "pix"}, "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", {"XPEAK_WORLD", "World-x coordinate of the brightest pixel",
&outobj2.peakxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", &outobj2.peakxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
"pos.eq.ra", "deg"}, "pos.eq.ra", "deg"},
...@@ -262,6 +268,12 @@ keystruct objkey[] = { ...@@ -262,6 +268,12 @@ keystruct objkey[] = {
{"Y_IMAGE_DBL", "Object position along y (double precision)", {"Y_IMAGE_DBL", "Object position along y (double precision)",
&outobj2.posy, H_FLOAT, T_DOUBLE, "%10.3f", "pixel", &outobj2.posy, H_FLOAT, T_DOUBLE, "%10.3f", "pixel",
"pos.cartesian.x;pos.barycenter;instr.det", "pix"}, "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", {"X_WORLD", "Barycenter position along world x axis",
&outobj2.mxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", &outobj2.mxw, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
"pos.eq.ra", "deg"}, "pos.eq.ra", "deg"},
...@@ -435,6 +447,13 @@ keystruct objkey[] = { ...@@ -435,6 +447,13 @@ keystruct objkey[] = {
&outobj2.winpos_y, H_FLOAT, T_DOUBLE, "%10.3f", "pixel", &outobj2.winpos_y, H_FLOAT, T_DOUBLE, "%10.3f", "pixel",
"pos.cartesian.y;instr.det", "pix"}, "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", {"XWIN_WORLD", "Windowed position along world x axis",
&outobj2.winpos_xw, H_FLOAT, T_DOUBLE, "%18.10e", "deg", &outobj2.winpos_xw, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
"pos.eq.ra", "deg"}, "pos.eq.ra", "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: 27/05/2009 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -50,6 +50,13 @@ ...@@ -50,6 +50,13 @@
&outobj2.y_prof, H_FLOAT, T_FLOAT, "%10.3f", "pixel", &outobj2.y_prof, H_FLOAT, T_FLOAT, "%10.3f", "pixel",
"pos.cartesian.y;stat.fit.param;instr.det;meta.main", "pix"}, "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", {"XMODEL_WORLD", "Fitted position along world x axis",
&outobj2.xw_prof, H_FLOAT, T_DOUBLE, "%18.10e", "deg", &outobj2.xw_prof, H_FLOAT, T_DOUBLE, "%18.10e", "deg",
"pos.eq.ra;stat.fit.param", "deg"}, "pos.eq.ra;stat.fit.param", "deg"},
......
...@@ -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: 29/05/2009 * Last modify: 13/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -168,7 +168,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the ...@@ -168,7 +168,7 @@ OUTPUT Pointer to an allocated fit structure (containing details about the
fit). fit).
NOTES It is a modified version of the lm_minimize() of lmfit. NOTES It is a modified version of the lm_minimize() of lmfit.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 29/05/2009 VERSION 13/07/2009
***/ ***/
void profit_fit(profitstruct *profit, void profit_fit(profitstruct *profit,
picstruct *field, picstruct *wfield, picstruct *field, picstruct *wfield,
...@@ -178,8 +178,10 @@ void profit_fit(profitstruct *profit, ...@@ -178,8 +178,10 @@ void profit_fit(profitstruct *profit,
patternstruct *pattern; patternstruct *pattern;
psfstruct *psf; psfstruct *psf;
checkstruct *check; checkstruct *check;
double psf_fwhm, a , cp,sp, emx2,emy2,emxy, dchi2, err; double *oldparaminit,
int i,j,p, nparam, ncomp, nprof; psf_fwhm, a , cp,sp, emx2,emy2,emxy, dchi2, err,
oldchi2;
int i,j,p, nparam, ncomp, nprof, oldniter;
nparam = profit->nparam; nparam = profit->nparam;
nprof = profit->nprof; nprof = profit->nprof;
...@@ -257,7 +259,21 @@ the_gal++; ...@@ -257,7 +259,21 @@ the_gal++;
/* Actual minimisation */ /* Actual minimisation */
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);
/*
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); QMEMCPY(profit->paraminit, oldparaminit, double, nparam);
if (profit_setparam(profit, PARAM_ARMS_PITCH, 160.0, 130.0, 175.0)==RETURN_OK) if (profit_setparam(profit, PARAM_ARMS_PITCH, 160.0, 130.0, 175.0)==RETURN_OK)
...@@ -593,7 +609,7 @@ INPUT Profile-fitting structure. ...@@ -593,7 +609,7 @@ INPUT Profile-fitting structure.
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 22/10/2008 VERSION 30/06/2009
***/ ***/
void profit_psf(profitstruct *profit) void profit_psf(profitstruct *profit)
{ {
...@@ -644,6 +660,7 @@ void profit_psf(profitstruct *profit) ...@@ -644,6 +660,7 @@ void profit_psf(profitstruct *profit)
*(pixout++) *= norm; *(pixout++) *= norm;
} }
return; return;
} }
...@@ -868,8 +885,8 @@ double *profit_compresi(profitstruct *profit, double dynparam, double *resi) ...@@ -868,8 +885,8 @@ double *profit_compresi(profitstruct *profit, double dynparam, double *resi)
val = *(objpix++); val = *(objpix++);
if ((wval=*(objweight++))>0.0) if ((wval=*(objweight++))>0.0)
{ {
val2 = (val - *lmodpix)*wval*invsig; val2 = (val - *lmodpix)*invsig;
val2 = val2>0.0? LOGF(1.0+val2) : -LOGF(1.0-val2); // val2 = val2>0.0? sqrt(1.0+val2) - 1.0 : 1.0 - sqrt(1.0-val2);
*(resit++) = val2*dynparam; *(resit++) = val2*dynparam;
error += val2*val2; error += val2*val2;
} }
...@@ -1341,7 +1358,7 @@ INPUT Profile-fitting structure. ...@@ -1341,7 +1358,7 @@ INPUT Profile-fitting structure.
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 27/05/2009 VERSION 29/06/2009
***/ ***/
void profit_moments(profitstruct *profit) void profit_moments(profitstruct *profit)
{ {
...@@ -1349,13 +1366,14 @@ void profit_moments(profitstruct *profit) ...@@ -1349,13 +1366,14 @@ void profit_moments(profitstruct *profit)
obj2struct *obj2; obj2struct *obj2;
double *pix, double *pix,
hw,hh, x,y, xstart, val, hw,hh, x,y, xstart, val,
mx,my, sum, mx2,my2,mxy, den; mx,my, sum, mx2,my2,mxy, den, r2max;
int ix,iy; int ix,iy;
obj = profit->obj; obj = profit->obj;
obj2 = profit->obj2; obj2 = profit->obj2;
hw = (double)(profit->modnaxisn[0]/2); hw = (double)(profit->modnaxisn[0]/2);
hh = (double)(profit->modnaxisn[1]/2); hh = (double)(profit->modnaxisn[1]/2);
r2max = hw<hh? hw*hw : hh*hh;
xstart = -hw; xstart = -hw;
y = -hh; y = -hh;
pix = profit->modpix; pix = profit->modpix;
...@@ -1364,15 +1382,18 @@ void profit_moments(profitstruct *profit) ...@@ -1364,15 +1382,18 @@ void profit_moments(profitstruct *profit)
{ {
x = xstart; x = xstart;
for (ix=profit->modnaxisn[0]; ix--; x+=1.0) for (ix=profit->modnaxisn[0]; ix--; x+=1.0)
{ if (y*y+x*x <= r2max)
val = *(pix++); {
sum += val; val = *(pix++);
mx += val*x; sum += val;
my += val*y; mx += val*x;
mx2 += val*x*x; my += val*y;
mxy += val*x*y; mx2 += val*x*x;
my2 += val*y*y; mxy += val*x*y;
} my2 += val*y*y;
}
else
pix++;
} }
if (sum <= 1.0/BIG) if (sum <= 1.0/BIG)
...@@ -1450,6 +1471,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) ...@@ -1450,6 +1471,7 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype)
obj = profit->obj; obj = profit->obj;
obj2 = profit->obj2; obj2 = profit->obj2;
param = parammin = parammax = 0.0; /* Avoid gcc -Wall warnings*/ param = parammin = parammax = 0.0; /* Avoid gcc -Wall warnings*/
switch(paramtype) switch(paramtype)
{ {
case PARAM_BACK: case PARAM_BACK:
...@@ -1619,6 +1641,35 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) ...@@ -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; p<nparam; p++)
parami[p] = 3.0*param[p] - 2.0*parami[p];
profit_unboundtobound(profit, profit->paraminit);
profit_unboundtobound(profit, param);
return param;
}
/****** profit_resetparams **************************************************** /****** profit_resetparams ****************************************************
PROTO void profit_resetparams(profitstruct *profit) PROTO void profit_resetparams(profitstruct *profit)
PURPOSE Set the initial, lower and upper boundary values of profile parameters. PURPOSE Set the initial, lower and upper boundary values of profile parameters.
...@@ -1640,6 +1691,7 @@ void profit_resetparams(profitstruct *profit) ...@@ -1640,6 +1691,7 @@ void profit_resetparams(profitstruct *profit)
} }
/****** profit_setparam **************************************************** /****** profit_setparam ****************************************************
PROTO void profit_setparam(profitstruct *profit, paramenum paramtype, PROTO void profit_setparam(profitstruct *profit, paramenum paramtype,
double param, double parammin, double parammax) double param, double parammin, double parammax)
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Include file for profit.c. * 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); ...@@ -144,6 +144,7 @@ profstruct *prof_init(profitstruct *profit, proftypenum profcode);
double *profit_compresi(profitstruct *profit, double dynparam, double *profit_compresi(profitstruct *profit, double dynparam,
double *resi), double *resi),
*profit_reresetparams(profitstruct *profit),
*profit_residuals(profitstruct *profit, picstruct *field, *profit_residuals(profitstruct *profit, picstruct *field,
picstruct *wfield, double dynparam, picstruct *wfield, double dynparam,
double *param, double *resi), double *param, double *resi),
......
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: global type definitions. * Contents: global type definitions.
* *
* Last modify: 18/03/2009 * Last modify: 20/07/2009
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -153,6 +153,7 @@ typedef struct ...@@ -153,6 +153,7 @@ typedef struct
/* ---- astrometric data */ /* ---- astrometric data */
double posx,posy; /* "FITS" pos. in pixels */ double posx,posy; /* "FITS" pos. in pixels */
double jacob[NAXIS*NAXIS]; /* Local deproject. Jacobian */ double jacob[NAXIS*NAXIS]; /* Local deproject. Jacobian */
double mamaposx,mamaposy; /* "MAMA" pos. in pixels */ double mamaposx,mamaposy; /* "MAMA" pos. in pixels */
float sposx,sposy; /* single precision pos. */ float sposx,sposy; /* single precision pos. */
float poserr_a, poserr_b, float poserr_a, poserr_b,
...@@ -170,7 +171,9 @@ typedef struct ...@@ -170,7 +171,9 @@ typedef struct
poserr_cxyw; /* WORLD error ellipse */ poserr_cxyw; /* WORLD error ellipse */
double mx2w,my2w,mxyw; /* WORLD var. and covar. */ double mx2w,my2w,mxyw; /* WORLD var. and covar. */
double peakxw, peakyw; /* WORLD of brightest pix */ double peakxw, peakyw; /* WORLD of brightest pix */
double peakxf, peakyf; /* FOCAL of brightest pix */
double mxw, myw; /* WORLD barycenters */ double mxw, myw; /* WORLD barycenters */
double mxf, myf; /* FOCAL barycenters */
double alphas, deltas; /* native alpha, delta */ double alphas, deltas; /* native alpha, delta */
float thetas; /* native position angle E/N*/ float thetas; /* native position angle E/N*/
double peakalphas, peakdeltas; /* native for brightest pix */ double peakalphas, peakdeltas; /* native for brightest pix */
...@@ -228,12 +231,13 @@ typedef struct ...@@ -228,12 +231,13 @@ typedef struct
float win_aw, win_bw, float win_aw, win_bw,
win_thetaw; /* WORLD ellipse parameters */ win_thetaw; /* WORLD ellipse parameters */
float win_polarw; /* WORLD WIN "polarization" */ float win_polarw; /* WORLD WIN "polarization" */
float win_thetas; /* native error pos. angle */ float win_thetas; /* native error pos. angle */
float win_theta2000; /* J2000 error pos. angle */ float win_theta2000; /* J2000 error pos. angle */
float win_theta1950; /* B1950 error pos. angle */ float win_theta1950; /* B1950 error pos. angle */
float win_cxxw, win_cyyw, float win_cxxw, win_cyyw,
win_cxyw; /* WORLD ellipse parameters */ win_cxyw; /* WORLD ellipse parameters */
double winpos_xw, winpos_yw; /* WORLD coordinates */ double winpos_xw, winpos_yw; /* WORLD coordinates */
double winpos_xf, winpos_yf; /* FOCAL coordinates */
double winpos_alphas, winpos_deltas; /* native alpha, delta */ double winpos_alphas, winpos_deltas; /* native alpha, delta */
double winpos_alpha2000, winpos_delta2000; /* J2000 alpha, delta */ double winpos_alpha2000, winpos_delta2000; /* J2000 alpha, delta */
double winpos_alpha1950, winpos_delta1950; /* B1950 alpha, delta */ double winpos_alpha1950, winpos_delta1950; /* B1950 alpha, delta */
...@@ -311,7 +315,8 @@ typedef struct ...@@ -311,7 +315,8 @@ typedef struct
float mag_prof; /* Mag from model-fitting */ float mag_prof; /* Mag from model-fitting */
float magerr_prof; /* RMS mag from model-fitting */ float magerr_prof; /* RMS mag from model-fitting */
float x_prof, y_prof; /* Coords 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 alphas_prof, deltas_prof; /* native alpha, delta */
double alpha2000_prof, delta2000_prof; /* J2000 alpha, delta */ double alpha2000_prof, delta2000_prof; /* J2000 alpha, delta */
double alpha1950_prof, delta1950_prof; /* B1950 alpha, delta */ double alpha1950_prof, delta1950_prof; /* B1950 alpha, delta */
......
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