Skip to content
GitLab
Commit b13bfdc5 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Added ... 

Added  POLAR1ERRMODEL_WORLD,POLAR2ERRMODEL_WORLD,POLARCORRMODEL_WORLD,ELLIP1ERRMODEL_WORLD,ELLIP2ERRMODEL_WORLD and ELLIPCORRMODEL_WORLD measurement parameters.
parent 3959d99d
Hide whitespace changes
Inline Side-by-side
src/astrom.c
View file @ b13bfdc5
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: Astrometrical computations. * Contents: Astrometrical computations.
* *
* Last modify: 03/08/2010 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -854,7 +854,9 @@ Compute profile-fitting shape parameters in WORLD and SKY coordinates. ...@@ -854,7 +854,9 @@ Compute profile-fitting shape parameters in WORLD and SKY coordinates.
void astrom_profshapeparam(picstruct *field, objstruct *obj) void astrom_profshapeparam(picstruct *field, objstruct *obj)
{ {
wcsstruct *wcs; wcsstruct *wcs;
double dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3, ct,st; double mat[9], tempmat[9], mx2wcov[9], dpdmx2[6], cov[4],
dx2,dy2,dxy, xm2,ym2,xym, pm2, lm0,lm1,lm2,lm3, ct,st,
temp, invstemp, den, invden, dval;
int lng,lat, naxis; int lng,lat, naxis;
wcs = field->wcs; wcs = field->wcs;
...@@ -1047,9 +1049,18 @@ void astrom_profshapeparam(picstruct *field, objstruct *obj) ...@@ -1047,9 +1049,18 @@ void astrom_profshapeparam(picstruct *field, objstruct *obj)
dx2 = obj2->prof_mx2; dx2 = obj2->prof_mx2;
dy2 = obj2->prof_my2; dy2 = obj2->prof_my2;
dxy = obj2->prof_mxy; dxy = obj2->prof_mxy;
obj2->prof_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy; mat[0] = lm0*lm0;
obj2->prof_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy; mat[3] = lm1*lm1;
obj2->prof_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy; mat[6] = lm0*lm1;
mat[1] = lm2*lm2;
mat[4] = lm3*lm3;
mat[7] = lm2*lm3;
mat[2] = lm0*lm2;
mat[5] = lm1*lm3;
mat[8] = lm0*lm3+lm1*lm2;
obj2->prof_mx2w = xm2 = mat[0]*dx2 + mat[3]*dy2 + mat[6]*dxy;
obj2->prof_my2w = ym2 = mat[1]*dx2 + mat[4]*dy2 + mat[7]*dxy;
obj2->prof_mxyw = xym = mat[2]*dx2 + mat[5]*dy2 + mat[8]*dxy;
temp=xm2-ym2; temp=xm2-ym2;
if (FLAG(obj2.prof_thetaw)) if (FLAG(obj2.prof_thetaw))
{ {
...@@ -1093,31 +1104,69 @@ void astrom_profshapeparam(picstruct *field, objstruct *obj) ...@@ -1093,31 +1104,69 @@ void astrom_profshapeparam(picstruct *field, objstruct *obj)
obj2->prof_cxyw = (float)(-2*xym/temp); obj2->prof_cxyw = (float)(-2*xym/temp);
} }
/*-- Use the Jacobians to compute the moment covariance matrix */
if (FLAG(obj2.prof_pol1errw) || FLAG(obj2.prof_e1errw))
propagate_covar(obj2->prof_mx2cov, mat, mx2wcov, 3, 3, tempmat);
if (FLAG(obj2.prof_pol1w)) if (FLAG(obj2.prof_pol1w))
{ {
if (xm2+ym2 > 1.0/BIG) if (xm2+ym2 > 1.0/BIG)
{ {
obj2->prof_pol1w = (xm2 - ym2) / (xm2+ym2); obj2->prof_pol1w = (xm2 - ym2) / (xm2+ym2);
obj2->prof_pol2w = 2.0*xym / (xm2 + ym2); obj2->prof_pol2w = 2.0*xym / (xm2 + ym2);
if (FLAG(obj2.prof_pol1errw))
{
/*-------- Compute the Jacobian of polarisation */
invden = 1.0/(xm2+ym2);
dpdmx2[0] = 2.0*ym2*invden*invden;
dpdmx2[1] = -2.0*xm2*invden*invden;
dpdmx2[2] = 0.0;
dpdmx2[3] = -2.0*xym*invden*invden;
dpdmx2[4] = -2.0*xym*invden*invden;
dpdmx2[5] = 2.0*invden;
propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
obj2->prof_pol1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_pol2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_pol12corrw = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
} }
else else
obj2->prof_pol1w = obj2->prof_pol2w = 0.0; obj2->prof_pol1w = obj2->prof_pol2w = obj2->prof_pol1errw
= obj2->prof_pol2errw = obj2->prof_pol12corrw = 0.0;
} }
if (FLAG(obj2.prof_e1w)) if (FLAG(obj2.prof_e1w))
{ {
if (xm2+ym2 > 1.0/BIG) if (xm2+ym2 > 1.0/BIG)
{ {
temp = xm2*ym2-xym*xym; temp = xm2*ym2 - xym*xym;
if (temp>=0.0) den = (temp>=0.0) ? xm2+ym2+2.0*sqrt(temp) : xm2+ym2;
temp = xm2+ym2+2.0*sqrt(temp); invden = 1.0/den;
else obj2->prof_e1w = (float)(invden*(xm2 - ym2));
temp = xm2+ym2; obj2->prof_e2w = (float)(2.0 * invden * xym);
obj2->prof_e1w = (xm2 - ym2) / temp; if (FLAG(obj2.prof_e1errw))
obj2->prof_e2w = 2.0*xym / temp; {
/*------ Compute the Jacobian of ellipticity */
invstemp = (temp>=0.0) ? 1.0/sqrt(temp) : 0.0;
dpdmx2[0] = ( den - (1.0+ym2*invstemp)*(xm2-ym2))*invden*invden;
dpdmx2[1] = (-den - (1.0+xm2*invstemp)*(xm2-ym2))*invden*invden;
dpdmx2[2] = 2.0*xym*invstemp*(xm2-ym2)*invden*invden;
dpdmx2[3] = -2.0*xym*(1.0+ym2*invstemp)*invden*invden;
dpdmx2[4] = -2.0*xym*(1.0+xm2*invstemp)*invden*invden;
dpdmx2[5] = (2.0*den+4.0*xym*xym*invstemp)*invden*invden;
/*------ Use the Jacobian to compute the ellipticity covariance matrix */
propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
obj2->prof_e1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
obj2->prof_e2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
obj2->prof_e12corrw = (dval=cov[0]*cov[3]) > 0.0?
(float)(cov[1]/sqrt(dval)) : 0.0;
}
} }
else else
obj2->prof_e1w = obj2->prof_e2w = 0.0; obj2->prof_e1w = obj2->prof_e2w = obj2->prof_e1errw
= obj2->prof_e2errw = obj2->prof_e12corrw = 0.0;
} }
} }
......
src/catout.c
View file @ b13bfdc5
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: functions for output of catalog data. * Contents: functions for output of catalog data.
* *
* Last modify: 03/08/2010 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -233,8 +233,11 @@ void updateparamflags() ...@@ -233,8 +233,11 @@ void updateparamflags()
FLAG(obj2.prof_bar_aspect) |= FLAG(obj2.prof_bar_aspecterr); FLAG(obj2.prof_bar_aspect) |= FLAG(obj2.prof_bar_aspecterr);
FLAG(obj2.prof_bar_theta) |= FLAG(obj2.prof_bar_thetaerr); FLAG(obj2.prof_bar_theta) |= FLAG(obj2.prof_bar_thetaerr);
FLAG(obj2.prof_arms_mag) |= FLAG(obj2.prof_arms_magerr); FLAG(obj2.prof_arms_mag) |= FLAG(obj2.prof_arms_magerr);
FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w); FLAG(obj2.prof_e1errw) |= FLAG(obj2.prof_e2errw) | FLAG(obj2.prof_e12corrw);
FLAG(obj2.prof_pol1w) |= FLAG(obj2.prof_pol2w); FLAG(obj2.prof_e1w) |= FLAG(obj2.prof_e2w) | FLAG(obj2.prof_e1errw);
FLAG(obj2.prof_pol1errw) |= FLAG(obj2.prof_pol2errw)
| FLAG(obj2.prof_pol12corrw);
FLAG(obj2.prof_pol1w) |= FLAG(obj2.prof_pol2w) | FLAG(obj2.prof_pol1errw);
FLAG(obj2.prof_aw) |= FLAG(obj2.prof_bw); FLAG(obj2.prof_aw) |= FLAG(obj2.prof_bw);
FLAG(obj2.prof_cxxw) |= FLAG(obj2.prof_cyyw) | FLAG(obj2.prof_cxyw); FLAG(obj2.prof_cxxw) |= FLAG(obj2.prof_cyyw) | FLAG(obj2.prof_cxyw);
FLAG(obj2.prof_thetas) |= FLAG(obj2.prof_theta1950) FLAG(obj2.prof_thetas) |= FLAG(obj2.prof_theta1950)
...@@ -325,14 +328,15 @@ void updateparamflags() ...@@ -325,14 +328,15 @@ void updateparamflags()
| FLAG(obj2.prof_arms_scalew) | FLAG(obj2.prof_arms_scalew)
| FLAG(obj2.prof_mx2w); | FLAG(obj2.prof_mx2w);
FLAG(obj2.proferr_e1) |= FLAG(obj2.proferr_e2) | FLAG(obj2.profcorr_e12); FLAG(obj2.prof_e1err) |= FLAG(obj2.prof_e2err) | FLAG(obj2.prof_e12corr)
| FLAG(obj2.prof_e1errw);
FLAG(obj2.prof_e1) |= FLAG(obj2.prof_e2) FLAG(obj2.prof_e1) |= FLAG(obj2.prof_e2)
| FLAG(obj2.proferr_e1) | FLAG(obj2.prof_e1w); | FLAG(obj2.prof_e1err) | FLAG(obj2.prof_e1w);
FLAG(obj2.proferr_pol1) |= FLAG(obj2.proferr_pol2) FLAG(obj2.prof_pol1err) |= FLAG(obj2.prof_pol2err)
| FLAG(obj2.profcorr_pol12); | FLAG(obj2.prof_pol12corr) | FLAG(obj2.prof_pol1errw);
FLAG(obj2.prof_pol1) |= FLAG(obj2.prof_pol2) FLAG(obj2.prof_pol1) |= FLAG(obj2.prof_pol2)
| FLAG(obj2.proferr_pol1) | FLAG(obj2.prof_pol1w); | FLAG(obj2.prof_pol1err) | FLAG(obj2.prof_pol1w);
FLAG(obj2.prof_a) |= FLAG(obj2.prof_b) | FLAG(obj2.prof_theta) FLAG(obj2.prof_a) |= FLAG(obj2.prof_b) | FLAG(obj2.prof_theta)
| FLAG(obj2.prof_aw); | FLAG(obj2.prof_aw);
FLAG(obj2.prof_cxx) |= FLAG(obj2.prof_cyy) FLAG(obj2.prof_cxx) |= FLAG(obj2.prof_cyy)
......
src/globals.h
View file @ b13bfdc5
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: global declarations. * Contents: global declarations.
* *
* Last modify: 01/10/2009 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -54,6 +54,8 @@ extern void alloccatparams(void), ...@@ -54,6 +54,8 @@ extern void alloccatparams(void),
neurclose(void), neurclose(void),
neurresp(double *, double *), neurresp(double *, double *),
preanalyse(int, objliststruct *, int), preanalyse(int, objliststruct *, int),
propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp),
readcatparams(char *), readcatparams(char *),
readdata(picstruct *, PIXTYPE *, int), readdata(picstruct *, PIXTYPE *, int),
readidata(picstruct *, FLAGTYPE *, int), readidata(picstruct *, FLAGTYPE *, int),
......
src/misc.c
View file @ b13bfdc5
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: miscellaneous functions. * Contents: miscellaneous functions.
* *
* Last modify: 24/09/2009 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -73,6 +73,60 @@ float hmedian(float *ra, int n) ...@@ -73,6 +73,60 @@ float hmedian(float *ra, int n)
} }
/****** propagate_covar ******************************************************
PROTO void propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp)
PURPOSE Compute Dt.V.D (propagate covariance matrix errors)
INPUT Pointer to the original covariance matrix,
pointer to the matrix of derivatives,
input number of parameters,
output number of parameters,
pointer to a ni*no work array.
OUTPUT -.
NOTES -.
AUTHOR E. Bertin (IAP)
VERSION 20/08/2010
***/
void propagate_covar(double *vi, double *d, double *vo,
int ni, int no, double *temp)
{
double *vit,*dt,*vot,*tempt,
dval;
int i,j,k;
tempt = temp;
vit = vi;
for (j=0; j<ni; j++)
{
dt = d;
for (i=no; i--;)
{
vit = vi + j*ni;
dval = 0.0;
for (k=ni; k--;)
dval += *(vit++)**(dt++);
*(tempt++) = dval;
}
}
vot = vo;
for (j=0; j<no; j++)
{
for (i=0; i<no; i++)
{
dt = d + j*ni;
tempt = temp + i;
dval = 0.0;
for (k=ni; k--; tempt+=no)
dval += *(dt++)**tempt;
*(vot++) = dval;
}
}
return;
}
/****** counter_seconds ******************************************************* /****** counter_seconds *******************************************************
PROTO double counter_seconds(void) PROTO double counter_seconds(void)
PURPOSE Count the number of seconds (with an arbitrary offset). PURPOSE Count the number of seconds (with an arbitrary offset).
......
src/paramprofit.h
View file @ b13bfdc5
...@@ -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: 03/08/2010 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -198,23 +198,23 @@ ...@@ -198,23 +198,23 @@
"src.ellipticity;stat.fit;instr.det", ""}, "src.ellipticity;stat.fit;instr.det", ""},
{"ELLIP1ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting", {"ELLIP1ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
&outobj2.proferr_e1, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_e1err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""}, "stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"ELLIP2ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting", {"ELLIP2ERRMODEL_IMAGE", "Ellipticity component std.error from model-fitting",
&outobj2.proferr_e2, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_e2err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""}, "stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"ELLIPCORRMODEL_IMAGE", "Corr.coeff between ellip.components from model-fitting", {"ELLIPCORRMODEL_IMAGE", "Corr.coeff between ellip.components from model-fitting",
&outobj2.profcorr_e12, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_e12corr, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit;instr.det", ""}, "stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
{"POLAR1ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting", {"POLAR1ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
&outobj2.proferr_pol1, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_pol1err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""}, "stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"POLAR2ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting", {"POLAR2ERRMODEL_IMAGE", "Polarisation component std.error from model-fitting",
&outobj2.proferr_pol2, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_pol2err, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit;instr.det", ""}, "stat.error;src.ellipticity;stat.fit;instr.det", ""},
{"POLARCORRMODEL_IMAGE", "Corr.coeff between polar. components from fitting", {"POLARCORRMODEL_IMAGE", "Corr.coeff between polar. components from fitting",
&outobj2.profcorr_pol12, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_pol12corr, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit;instr.det", ""}, "stat.correlation;src.ellipticity;stat.fit;instr.det", ""},
{"X2MODEL_WORLD", "Variance along X-WORLD (alpha) from model-fitting", {"X2MODEL_WORLD", "Variance along X-WORLD (alpha) from model-fitting",
...@@ -239,6 +239,26 @@ ...@@ -239,6 +239,26 @@
&outobj2.prof_pol2w, H_FLOAT, T_FLOAT, "%10.6f", "", &outobj2.prof_pol2w, H_FLOAT, T_FLOAT, "%10.6f", "",
"src.ellipticity;stat.fit", ""}, "src.ellipticity;stat.fit", ""},
{"ELLIP1ERRMODEL_WORLD", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e1errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"ELLIP2ERRMODEL_WORLD", "Ellipticity component std.error from model-fitting",
&outobj2.prof_e2errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"ELLIPCORRMODEL_WORLD", "Corr.coeff between ellip.components from model-fitting",
&outobj2.prof_e12corrw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit", ""},
{"POLAR1ERRMODEL_WORLD", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol1errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"POLAR2ERRMODEL_WORLD", "Polarisation component std.error from model-fitting",
&outobj2.prof_pol2errw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.error;src.ellipticity;stat.fit", ""},
{"POLARCORRMODEL_WORLD", "Corr.coeff between polar. components from fitting",
&outobj2.prof_pol12corrw, H_FLOAT, T_FLOAT, "%10.6f", "",
"stat.correlation;src.ellipticity;stat.fit", ""},
{"CXXMODEL_IMAGE", "Cxx ellipse parameter from model-fitting", {"CXXMODEL_IMAGE", "Cxx ellipse parameter from model-fitting",
&outobj2.prof_cxx, H_EXPO, T_FLOAT, "%15.7e", "pixel**(-2)", &outobj2.prof_cxx, H_EXPO, T_FLOAT, "%15.7e", "pixel**(-2)",
"src.impactParam;stat.fit;instr.det", "pix-2"}, "src.impactParam;stat.fit;instr.det", "pix-2"},
......
src/profit.c
View file @ b13bfdc5
...@@ -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: 03/08/2010 * Last modify: 19/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -2027,15 +2027,17 @@ INPUT Profile-fitting structure, ...@@ -2027,15 +2027,17 @@ INPUT Profile-fitting structure,
OUTPUT -. OUTPUT -.
NOTES -. NOTES -.
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 03/08/2010 VERSION 20/08/2010
***/ ***/
void profit_moments(profitstruct *profit, obj2struct *obj2) void profit_moments(profitstruct *profit, obj2struct *obj2)
{ {
profstruct *prof; profstruct *prof;
double *jac,*jact, *pjac,*pjact, double dpdmx2[6], cov[4],
*dmx2,*dmx2t,*dmy2,*dmy2t,*dmxy,*dmxyt, *de1,*de1t,*de2,*de2t, *jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
m0,invm0, mx2,my2,mxy, invden, temp,temp2,invstemp2,temp3, *dmx2,*dmy2,*dmxy,
pmx2,theta, flux, var1,var2,cov, dval; m0,invm0, mx2,my2,mxy, den,invden,
temp, temp2,invtemp2,invstemp2,
pmx2,theta, flux, dval;
float *covart; float *covart;
int findex[PROF_NPROF], int findex[PROF_NPROF],
i,j,p, nparam; i,j,p, nparam;
...@@ -2074,17 +2076,22 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2074,17 +2076,22 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
/* obj2->prof_mxy = mxy = mxy/sum - mx*my;*/ /* obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
nparam = profit->nparam; nparam = profit->nparam;
if (FLAG(obj2.proferr_e1) || FLAG(obj2.proferr_pol1)) if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
{ {
/*-- Set up Jacobian matrices */ /*-- Set up Jacobian matrices */
QCALLOC(jac, double, nparam*3); QCALLOC(jac, double, nparam*3);
QMALLOC(pjac, double, nparam*3); QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
QMALLOC(dcovar, double, nparam*nparam);
dcovart = dcovar;
covart = profit->covar;
for (i=nparam*nparam; i--;)
*(dcovart++) = (double)(*(covart++));
dmx2 = jac; dmx2 = jac;
dmy2 = jac+nparam; dmy2 = jac+nparam;
dmxy = jac+2*nparam; dmxy = jac+2*nparam;
} }
else else
jac = pjac = NULL; jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;
m0 = mx2 = my2 = mxy = 0.0; m0 = mx2 = my2 = mxy = 0.0;
for (p=0; p<profit->nprof; p++) for (p=0; p<profit->nprof; p++)
...@@ -2102,9 +2109,6 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2102,9 +2109,6 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
pjact = pjac; pjact = pjac;
for (j=nparam*3; j--;) for (j=nparam*3; j--;)
*(jact++) += flux * *(pjact++); *(jact++) += flux * *(pjact++);
dmx2[findex[p]] += prof->mx2;
dmy2[findex[p]] += prof->my2;
dmxy[findex[p]] += prof->mxy;
} }
} }
invm0 = 1.0 / m0; invm0 = 1.0 / m0;
...@@ -2116,9 +2120,10 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2116,9 +2120,10 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
{ {
for (p=0; p<profit->nprof; p++) for (p=0; p<profit->nprof; p++)
{ {
dmx2[findex[p]] -= mx2; prof = profit->prof[p];
dmy2[findex[p]] -= my2; dmx2[findex[p]] = prof->mx2 - mx2;
dmxy[findex[p]] -= mxy; dmy2[findex[p]] = prof->my2 - my2;
dmxy[findex[p]] = prof->mxy - mxy;
} }
jact = jac; jact = jac;
for (j=nparam*3; j--;) for (j=nparam*3; j--;)
...@@ -2133,132 +2138,83 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2133,132 +2138,83 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
temp2 = mx2*my2-mxy*mxy; temp2 = mx2*my2-mxy*mxy;
} }
/* Use the Jacobians to compute the moment covariance matrix */
if (jac)
propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
pjac); /* We re-use pjac */
if (FLAG(obj2.prof_pol1)) if (FLAG(obj2.prof_pol1))
{ {
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
if (mx2+my2 > 1.0/BIG) if (mx2+my2 > 1.0/BIG)
{ {
obj2->prof_pol1 = (mx2 - my2) / (mx2+my2); obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
obj2->prof_pol2 = 2.0*mxy / (mx2 + my2); obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
if (jac) if (FLAG(obj2.prof_pol1err))
{ {
/*------ Compute Jacobian of polarisation and ellipticity */ /*------ Compute the Jacobian of polarisation */
if (FLAG(obj2.proferr_pol1)) invden = 1.0/(mx2+my2);
{ dpdmx2[0] = 2.0*my2*invden*invden;
invden = 2.0/((mx2+my2)*(mx2+my2)); dpdmx2[1] = -2.0*mx2*invden*invden;
/*-------- Ellipticity */ dpdmx2[2] = 0.0;
dmx2t = dmx2; dpdmx2[3] = -2.0*mxy*invden*invden;
dmy2t = dmy2; dpdmx2[4] = -2.0*mxy*invden*invden;
dmxyt = dmxy; dpdmx2[5] = 2.0*invden;
de1t = de1 = pjac; /* We re-use pjac */
de2t = de2 = pjac + nparam; /* We re-use pjac */ /*------ Use the Jacobian to compute the polarisation covariance matrix */
for (j=nparam; j--;) propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
{ pjac); /* We re-use pjac */
*(de1t++) = invden*(*dmx2t*my2-*dmy2t*mx2); obj2->prof_pol1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
*(de2t++) = invden*(*(dmxyt++)*(mx2+my2) obj2->prof_pol2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
- mxy*(*(dmx2t++)-*(dmy2t++))); obj2->prof_pol12corr = (dval=cov[0]*cov[3]) > 0.0?
} (float)(cov[1]/sqrt(dval)) : 0.0;
/*-------- Use the Jacobians to compute errors */
covart = profit->covar;
var1 = 0.0;
for (j=0; j<nparam; j++)
{
dval = 0.0;
de1t = de1;
for (i=nparam; i--;)
dval += *(covart++)**(de1t++);
var1 += dval*de1[j];
}
obj2->proferr_pol1 = (float)sqrt(var1<0.0? 0.0: var1);
covart = profit->covar;
cov = var2 = 0.0;
for (j=0; j<nparam; j++)
{
dval = 0.0;
de2t = de2;
for (i=nparam; i--;)
dval += *(covart++)**(de2t++);
cov += dval*de1[j];
var2 += dval*de2[j];
}
obj2->proferr_pol2 = (float)sqrt(var2<0.0? 0.0: var2);
obj2->profcorr_pol12 = (dval=var1*var2) > 0.0? (float)(cov/sqrt(dval))
: 0.0;
}
} }
} }
else else
obj2->prof_pol1 = obj2->prof_pol2 obj2->prof_pol1 = obj2->prof_pol2
= obj2->proferr_pol1 = obj2->proferr_pol2 = obj2->profcorr_pol12 = 0.0; = obj2->prof_pol1err = obj2->prof_pol2err = obj2->prof_pol12corr = 0.0;
} }
if (FLAG(obj2.prof_e1)) if (FLAG(obj2.prof_e1))
{ {
/*--- "Ellipticity", i.e. module = (a-b)/(a+b) */
if (mx2+my2 > 1.0/BIG) if (mx2+my2 > 1.0/BIG)
{ {
if (temp2>=0.0) den = (temp2>=0.0) ? mx2+my2+2.0*sqrt(temp2) : mx2+my2;
invden = 1.0/(mx2+my2+2.0*sqrt(temp2)); invden = 1.0/den;
else
invden = 1.0/(mx2+my2);
obj2->prof_e1 = (float)(invden * (mx2 - my2)); obj2->prof_e1 = (float)(invden * (mx2 - my2));
obj2->prof_e2 = (float)(2.0 * invden * mxy); obj2->prof_e2 = (float)(2.0 * invden * mxy);
if (jac) if (FLAG(obj2.prof_e1err))
{ {
/*------ Compute Jacobian of polarisation and ellipticity */ /*------ Compute the Jacobian of ellipticity */
invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0; invstemp2 = (temp2>=0.0) ? 1.0/sqrt(temp2) : 0.0;
if (FLAG(obj2.proferr_e1)) dpdmx2[0] = ( den - (1.0+my2*invstemp2)*(mx2-my2))*invden*invden;
{ dpdmx2[1] = (-den - (1.0+mx2*invstemp2)*(mx2-my2))*invden*invden;
/*-------- Ellipticity */ dpdmx2[2] = 2.0*mxy*invstemp2*(mx2-my2)*invden*invden;
dmx2t = dmx2; dpdmx2[3] = -2.0*mxy*(1.0+my2*invstemp2)*invden*invden;
dmy2t = dmy2; dpdmx2[4] = -2.0*mxy*(1.0+mx2*invstemp2)*invden*invden;
dmxyt = dmxy; dpdmx2[5] = (2.0*den+4.0*mxy*mxy*invstemp2)*invden*invden;
de1t = de1 = pjac; /* We re-use pjac */
de2t = de2 = pjac + nparam; /* We re-use pjac */ /*------ Use the Jacobian to compute the ellipticity covariance matrix */
for (j=nparam; j--;) propagate_covar(obj2->prof_mx2cov, dpdmx2, cov, 3, 2,
{ pjac); /* We re-use pjac */
temp3 = invden*(*dmx2t+*dmy2t obj2->prof_e1err = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
+ invstemp2 * (*dmx2t*my2+mx2**dmy2t-2.0*mxy**dmxyt)); obj2->prof_e2err = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
*(de1t++) = invden*(*(dmx2t++) - *(dmy2t++) - (mx2-my2)*temp3); obj2->prof_e12corr = (dval=cov[0]*cov[3]) > 0.0?
*(de2t++) = 2.0*invden*(*(dmxyt++) - mxy*temp3); (float)(cov[1]/sqrt(dval)) : 0.0;
}
/*-------- Use the Jacobians to compute errors */
covart = profit->covar;
var1 = 0.0;
for (j=0; j<nparam; j++)
{
dval = 0.0;
de1t = de1;
for (i=nparam; i--;)
dval += *(covart++)**(de1t++);
var1 += dval*de1[j];
}
obj2->proferr_e1 = (float)sqrt(var1<0.0? 0.0: var1);
covart = profit->covar;
cov = var2 = 0.0;
for (j=0; j<nparam; j++)
{
dval = 0.0;
de2t = de2;
for (i=nparam; i--;)
dval += *(covart++)**(de2t++);
cov += dval*de1[j];
var2 += dval*de2[j];
}
obj2->proferr_e2 = (float)sqrt(var2<0.0? 0.0: var2);
obj2->profcorr_e12 = (dval=var1*var2) > 0.0? (float)(cov/sqrt(dval))
: 0.0;
}
} }
} }
else else
obj2->prof_e1 = obj2->prof_e2 obj2->prof_e1 = obj2->prof_e2
= obj2->proferr_e1 = obj2->proferr_e2 = obj2->profcorr_e12 = 0.0; = obj2->prof_e1err = obj2->prof_e2err = obj2->prof_e12corr = 0.0;
} }
if (FLAG(obj2.prof_cxx)) if (FLAG(obj2.prof_cxx))
{ {
obj2->prof_cxx = (float)(my2/temp2); invtemp2 = (temp2>=0.0) ? 1.0/temp2 : 0.0;
obj2->prof_cyy = (float)(mx2/temp2); obj2->prof_cxx = (float)(my2*invtemp2);
obj2->prof_cxy = (float)(-2*mxy/temp2); obj2->prof_cyy = (float)(mx2*invtemp2);
obj2->prof_cxy = (float)(-2*mxy*invtemp2);
} }
if (FLAG(obj2.prof_a)) if (FLAG(obj2.prof_a))
...@@ -2278,6 +2234,7 @@ void profit_moments(profitstruct *profit, obj2struct *obj2) ...@@ -2278,6 +2234,7 @@ void profit_moments(profitstruct *profit, obj2struct *obj2)
/* Free memory used by Jacobians */ /* Free memory used by Jacobians */
free(jac); free(jac);
free(pjac); free(pjac);
free(dcovar);
return; return;
} }
...@@ -3699,7 +3656,7 @@ INPUT Profile-fitting structure, ...@@ -3699,7 +3656,7 @@ INPUT Profile-fitting structure,
OUTPUT Index to the profile flux for further processing. OUTPUT Index to the profile flux for further processing.
NOTES . NOTES .
AUTHOR E. Bertin (IAP) AUTHOR E. Bertin (IAP)
VERSION 21/07/2010 VERSION 20/08/2010
***/ ***/
int prof_moments(profitstruct *profit, profstruct *prof, double *jac) int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
{ {
...@@ -3717,6 +3674,13 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3717,6 +3674,13 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
dmy2 = jac + nparam; dmy2 = jac + nparam;
dmxy = jac + 2*nparam; dmxy = jac + 2*nparam;
} }
else
dmx2 = dmy2 = dmxy = NULL; /* To avoid gcc -Wall warnings */
m20 = 0.0; /* to avoid gcc -Wall warnings */
index = 0; /* to avoid gcc -Wall warnings */
if (prof->posangle) if (prof->posangle)
{ {
a2 = *prof->aspect**prof->aspect; a2 = *prof->aspect**prof->aspect;
...@@ -3731,6 +3695,8 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3731,6 +3695,8 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
ds2 = 2.0*ct*st*DEG; ds2 = 2.0*ct*st*DEG;
dcs = (ct*ct - st*st)*DEG; dcs = (ct*ct - st*st)*DEG;
} }
else
dc2 = ds2 = dcs = 0.0; /* To avoid gcc -Wall warnings */
switch(prof->code) switch(prof->code)
{ {
case PROF_SERSIC: case PROF_SERSIC:
...@@ -3829,8 +3795,6 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac) ...@@ -3829,8 +3795,6 @@ int prof_moments(profitstruct *profit, profstruct *prof, double *jac)
index = profit->paramindex[PARAM_OUTRING_FLUX]; index = profit->paramindex[PARAM_OUTRING_FLUX];
break; break;
default: default:
m20 = 0.0; /* to avoid gcc -Wall warnings */
index = 0; /* to avoid gcc -Wall warnings */
error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ", error(EXIT_FAILURE, "*Internal Error*: Unknown oriented model in ",
"prof_moments()"); "prof_moments()");
break; break;
......
src/types.h
View file @ b13bfdc5
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* *
* Contents: global type definitions. * Contents: global type definitions.
* *
* Last modify: 21/07/2010 * Last modify: 20/08/2010
* *
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/ */
...@@ -350,16 +350,17 @@ typedef struct ...@@ -350,16 +350,17 @@ typedef struct
float poserrcxxw_prof, poserrcyyw_prof, float poserrcxxw_prof, poserrcyyw_prof,
poserrcxyw_prof; /* WORLD error ellipse */ poserrcxyw_prof; /* WORLD error ellipse */
double prof_mx2, prof_my2, prof_mxy; /* Model-fitting moments */ double prof_mx2, prof_my2, prof_mxy; /* Model-fitting moments */
double prof_mx2cov[9]; /* Mod-fit moment cov. matrix */
float prof_a, prof_b, float prof_a, prof_b,
prof_theta; /* Model-fitting ellip. params*/ prof_theta; /* Model-fitting ellip. params*/
float prof_cxx, prof_cyy, float prof_cxx, prof_cyy,
prof_cxy; /* Model-fitting ellip. params*/ prof_cxy; /* Model-fitting ellip. params*/
float prof_pol1, prof_pol2; /* Model-fitting pol. vector*/ float prof_pol1, prof_pol2; /* Model-fitting pol. vector*/
float proferr_pol1, proferr_pol2, float prof_pol1err, prof_pol2err,
profcorr_pol12; /* Model-fitting pol. errors */ prof_pol12corr; /* Model-fitting pol. errors */
float prof_e1, prof_e2; /* Model-fitting ellip.vector*/ float prof_e1, prof_e2; /* Model-fitting ellip.vector*/
float proferr_e1, proferr_e2, float prof_e1err, prof_e2err,
profcorr_e12; /* Model-fitting ellip. errors */ prof_e12corr; /* Model-fitting ellip. errors */
double prof_mx2w, prof_my2w, double prof_mx2w, prof_my2w,
prof_mxyw; /* WORLD model-fitting moments*/ prof_mxyw; /* WORLD model-fitting moments*/
float prof_aw, prof_bw, float prof_aw, prof_bw,
...@@ -370,7 +371,11 @@ typedef struct ...@@ -370,7 +371,11 @@ typedef struct
float prof_cxxw, prof_cyyw, float prof_cxxw, prof_cyyw,
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_pol1errw, prof_pol2errw,
prof_pol12corrw; /* WORLD polarisation errors */
float prof_e1w, prof_e2w; /* WORLD ellipticity vector*/ float prof_e1w, prof_e2w; /* WORLD ellipticity vector*/
float prof_e1errw, prof_e2errw,
prof_e12corrw; /* WORLD ellipticity errors */
float prof_class_star; /* Mod.-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_concentrationerr; /* RMS error */
......
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