Commit f18f93e3 authored by Johan Richard's avatar Johan Richard
Browse files

Remove optimisation of sigy and only use sigx and eps in shapemodel

parent 30192305
......@@ -109,7 +109,7 @@ double dist(struct point A, struct point B);
double dlumcosmo1(double z);
void do_itos(double **im, struct pixlist *pl, int npl, double dlsds, double zs, double **source, double **erreur, int **imult);
double d_poisson(double xm, int *idum);
double d_profil(double x, double y, const struct galaxie *gal);
double d_profil(double x, double y, struct galaxie *gal);
double d_random(int *idum);
double dratio(double zl, double zs);
void dratio_gal(struct galaxie *arclet, double zl);
......
......@@ -772,9 +772,6 @@ char* getParName(int ipx, char* name, int type)
case(SA):
strcpy( name, "a (arcsec)" );
break;
case(SB):
strcpy( name, "b (arcsec)" );
break;
case(SEPS):
strcpy( name, "eps" );
break;
......
......@@ -26,6 +26,7 @@ void checkpar()
const extern struct g_mode M;
extern struct g_source S;
extern struct galaxie source[NFMAX];
extern int sblock[NFMAX][NPAMAX];
FILE *OUT;
int i;
......@@ -34,10 +35,13 @@ void checkpar()
//Double check galaxy sources
for (i=0; i<S.ns;i++)
{
// ellipticity has been defined in shapemodel
if(source[i].eps>=0)
{
if(source[i].E.a>=0)
// size has been defined in shapemodel
if(source[i].E.a>0)
{
// three parameters sigx,sigy,eps defined in shapemodel: sigy ignored
if(source[i].E.b>0)
fprintf(stderr, "WARNING: source %d sig_y parameter readjusted according to eps\n", i);
source[i].E.b = source[i].E.a * (1. - source[i].eps) / (1. + source[i].eps);
......@@ -48,6 +52,26 @@ void checkpar()
exit(-1);
}
}
// ellipticity has not been defined in shapemodel
else
{
// sig_x and sig_y defined in shapemodel: compute eps, forbid optimisation
if((source[i].E.a>0)&&(source[i].E.b>0))
{
source[i].eps = (source[i].E.a - source[i].E.b) / (source[i].E.a + source[i].E.b);
//forbid optimisation HERE
fprintf(stderr, "WARNING: ellipticity for source %d not defined. Source shape optimisation forbidden\n", i);
sblock[i][SA]=0;
sblock[i][SB]=0;
sblock[i][SEPS]=0;
}
else if(source[i].E.a>0) // only sig_x defined in shapemodel
{
fprintf(stderr, "WARNING: ellipticity for source %d not defined. Assuming 0\n", i);
source[i].eps = 0;
source[i].E.b=source[i].E.a;
}
}
}
return;
......
......@@ -23,7 +23,7 @@
* alpha=1/(e*l)
*/
double d_profil(double x, double y, const struct galaxie *gal)
double d_profil(double x, double y, struct galaxie *gal)
{
double xx, yy, xxx, yyy;
double res;
......@@ -32,6 +32,10 @@ double d_profil(double x, double y, const struct galaxie *gal)
res = 0;
//jrichard
//Recompute E.b according to E.a and eps
gal->E.b=gal->E.a*(1-gal->eps)/(1+gal->eps);
if (gal->c == 's')
{
xx = x - gal->C.x;
......
......@@ -589,9 +589,7 @@ void o_print_res(double chi0, double evidence)
fprintf(besto, "\ts_angle %d %.6lf %.6lf\n", sblock[i][STHETA], smin[i].E.theta * RTD, smax[i].E.theta * RTD);
if( sblock[i][SA] )
fprintf(besto, "\ts_sigx %d %.6lf %.6lf\n", sblock[i][SA], smin[i].E.a, smax[i].E.a);
if( sblock[i][SB] )
fprintf(besto, "\ts_sigy %d %.6lf %.6lf\n", sblock[i][SB], smin[i].E.b, smax[i].E.b);
if( sblock[i][SB] )
if( sblock[i][SEPS] )
fprintf(besto, "\ts_eps %d %.6lf %.6lf\n", sblock[i][SEPS], smin[i].eps, smax[i].eps);
if( sblock[i][SFLUX] )
fprintf(besto, "\tmag %d %.6lf %.6lf\n", sblock[i][SFLUX], smin[i].mag, smax[i].mag);
......
......@@ -261,9 +261,6 @@ int rescaleParam(int id, int type, int ipx, double *pval)
case(SA):
val = prior(sblock[id][SA], val, smin[id].E.a, smax[id].E.a);
break;
case(SB):
val = prior(sblock[id][SB], val, smin[id].E.b, smax[id].E.b);
break;
case(SEPS):
val = prior(sblock[id][SEPS], val, smin[id].eps, smax[id].eps);
break;
......@@ -508,7 +505,6 @@ int rescaleCube_1Atom(double *cube, int npar)
extern int block[][NPAMAX];
extern int cblock[NPAMAX];
extern int sblock[NFMAX][NPAMAX];
extern struct galaxie source[NFMAX];
extern int vfblock[NPAMAX];
extern double *np_b0; // non parametric accelerator (o_global.c)
......@@ -549,13 +545,9 @@ int rescaleCube_1Atom(double *cube, int npar)
// Rescale the source parameters
if ( M.iclean == 2 )
for ( k = 0; k < S.ns; k++ )
{
for ( ipx = SCX; ipx <= SFLUX; ipx++ )
if ( sblock[k][ipx] != 0 )
valid *= rescaleParam(k, SOURCES, ipx, &cube[ipar++]);
if( source[k].eps>=0 )
source[k].E.b = source[k].E.a * (1. - source[k].eps) / (1. + source[k].eps);
}
// rescale the cosmological parameters
for ( ipx = OMEGAM; ipx <= WA; ipx++ )
......@@ -614,13 +606,6 @@ int rescaleCube_1Atom(double *cube, int npar)
for ( ilens = 0; ilens < G.no_lens; ilens++ )
if ( lens[ilens].b0 != 0. && lens[ilens].rc > lens[ilens].rcut ) valid = 0;
if ( M.iclean == 2 )
{
extern struct galaxie source[NFMAX];
for ( k = 0; k < S.ns; k++ )
if( source[k].E.a < source[k].E.b ) valid = 0;
}
return valid;
}
......
......@@ -90,8 +90,8 @@ void o_set_lmax(int i, int ipx, double x)
case(SA):
smax[i].E.a = x;
break;
case(SB):
smax[i].E.b = x;
case(SEPS):
smax[i].eps = x;
break;
case(STHETA):
smax[i].E.theta = x;
......
......@@ -91,8 +91,8 @@ void o_set_lmin(int i, int ipx, double x)
case(SA):
smin[i].E.a = x;
break;
case(SB):
smin[i].E.b = x;
case(SEPS):
smin[i].eps = x;
break;
case(STHETA):
smin[i].E.theta = x;
......
......@@ -26,9 +26,6 @@ void o_set_source(struct galaxie *source, int ipx, double val)
case(SA):
source->E.a = val;
break;
case(SB):
source->E.b = val;
break;
case(SEPS):
source->eps = val;
break;
......@@ -59,9 +56,6 @@ double o_get_source(struct galaxie *source, int ipx)
case(SA):
val = source->E.a;
break;
case(SB):
val = source->E.b;
break;
case(SEPS):
val = source->eps;
break;
......
......@@ -57,10 +57,10 @@ void r_shapelimit(FILE *IN, FILE *OUT, long int i)
else if (!strcmp(second,"s_sigy") ||
!strcmp(second,"b_arcsec") )
{
sscanf(third, "%d %lf %lf", &sblock[i][SB],
&smin[i].E.b, &smax[i].E.b);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SB],
smin[i].E.b, smax[i].E.b);
//jrichard not allowed any more
//EXIT with ERROR
fprintf(stderr, "ERROR: sig_y parameter cannot be optimised in source %ld. Use ellipticity (eps) instead.\n",i);
exit(-1);
}
else if (!strcmp(second,"s_eps") )
{
......
......@@ -294,8 +294,6 @@ void setBayesModel( long int iVal, long int nVal, double **array)
if( sblock[i][ipx] != 0 )
o_set_source( &source[i], ipx, method(nVal, iParam++, array));
if( sblock[i][SEPS] != 0 )
source[i].E.b = source[i].E.a * (1. - source[i].eps) / (1. + source[i].eps);
}
}
......
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