Commit 1993aedf authored by Johan Richard's avatar Johan Richard
Browse files

Merge branch 'sigmaplus_ellipse'

parents 3e6583b0 ce2c9637
......@@ -32,14 +32,19 @@
#define STHETA 24
#define SINDEX 25
#define SFLUX 26
#define VFCX 27
#define VFCY 28
#define VFVT 29
#define VFRT 30
#define VFI 31
#define VFTHETA 32
#define VFLCENT 33
#define VFSIGMA 34
#define SA2 27
#define SEPS2 28
#define STHETA2 29
#define SINDEX2 30
#define SFLUX2 31
#define VFCX 32
#define VFCY 33
#define VFVT 34
#define VFRT 35
#define VFI 36
#define VFTHETA 37
#define VFLCENT 38
#define VFSIGMA 39
/*
* structure definition
......@@ -584,10 +589,12 @@ struct galaxie
struct point C;
struct point Grad; // total deflection with all clumps
struct ellipse E;
struct ellipse E2;
char c;
int type;
double magabs;
double mag;
double mag2;
double flux;
double mu;
double I0;
......@@ -597,6 +604,7 @@ struct galaxie
double dr; // ratio dl0s/dos
double q;
double eps;
double eps2;
double tau;
double dis;
double A;
......
......@@ -784,6 +784,21 @@ char* getParName(int ipx, char* name, int type)
case(SFLUX):
strcpy( name, "mag" );
break;
case(SA2):
strcpy( name, "a2 (arcsec)" );
break;
case(SEPS2):
strcpy( name, "eps2" );
break;
case(STHETA2):
strcpy( name, "theta2 (deg)" );
break;
case(SINDEX2):
strcpy( name, "index2" );
break;
case(SFLUX2):
strcpy( name, "mag2" );
break;
case(VFCX):
strcpy( name, "velocity x (arcsec)");
break;
......@@ -1161,11 +1176,12 @@ int UserMonitor( // O 0 = continue, +ve = finish, -ve = abort
extern struct galaxie source[NFMAX];
for( i = 0; i < S.ns; i++ )
for( ipx =SCX; ipx <= SFLUX; ipx++ )
for( ipx =SCX; ipx <= SFLUX2; ipx++ )
if( sblock[i][ipx] != 0 )
{
val = o_get_source(&source[i], ipx);
if( ipx == STHETA ) val *= RTD;
if( ipx == STHETA2 ) val *= RTD;
fprintf(bayes, "%f ", val);
addStat(val, UserCommon, iparam++, samp);
}
......
......@@ -73,6 +73,32 @@ void checkpar()
source[i].E.b=source[i].E.a;
}
}
if((source[i].eps2>0)||(source[i].E2.a>0))
{
if(source[i].type!=6)
{
fprintf(stderr, "ERROR: wrong type (should be 6) for 2nd Sersic profile in source %d\n", i);
exit(-1);
}
if(source[i].eps2>=0)
{
// second size has been defined in shapemodel
if(source[i].E2.a>0)
source[i].E2.b = source[i].E2.a * (1. - source[i].eps2) / (1. + source[i].eps2);
else
{
fprintf(stderr, "ERROR: 2nd size not defined in source %d type 6\n", i);
exit(-1);
}
}
else // only size defined
{
fprintf(stderr, "WARNING: second ellipticity for source %d not defined. Assuming 0\n", i);
source[i].eps2 = 0;
source[i].E2.b=source[i].E2.a;
}
}
}
return;
......
......@@ -26,11 +26,12 @@
double d_profil(double x, double y, const struct galaxie *gal)
{
double xx, yy, xxx, yyy;
double res;
double res,res2;
const extern struct g_observ O;
const extern struct g_large L;
res = 0;
res2 = 0;
//jrichard
//Test E.b according to E.a and eps
......@@ -119,6 +120,19 @@ double d_profil(double x, double y, const struct galaxie *gal)
else
res = 0.;
}
/* Sersic */
else if (L.vitesse == 6 || gal->type == 6)
{
xx = ((x - gal->C.x) * sin(gal->E.theta) - (y - gal->C.y) * cos(gal->E.theta)) / gal->E.b;
yy = ((y - gal->C.y) * sin(gal->E.theta) + (x - gal->C.x) * cos(gal->E.theta)) / gal->E.a;
res = exp(pow(xx * xx + yy * yy, 0.5/gal->var1));
res = pow(10., (26. - gal->mag) / 2.5) / res;
xx = ((x - gal->C.x) * sin(gal->E2.theta) - (y - gal->C.y) * cos(gal->E2.theta)) / gal->E2.b;
yy = ((y - gal->C.y) * sin(gal->E2.theta) + (x - gal->C.x) * cos(gal->E2.theta)) / gal->E2.a;
res2 = exp(pow(xx * xx + yy * yy, 0.5/gal->var2));
res2 = pow(10., (26. - gal->mag2) / 2.5) / res2;
res += res2;
}
else
{
fprintf(stderr, "ERROR: source %s brightness profil type %d unknown\n", gal->n, gal->type);
......
......@@ -947,7 +947,7 @@ int getNParameters()
// source parameters
for ( i = 0; i < S.ns; i++ )
for (ipx = SCX; ipx <= SFLUX; ipx++ )
for (ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[i][ipx] != 0 )
parameters++;
......
......@@ -285,15 +285,24 @@ void o_print_res(double chi0, double evidence)
// SHAPE MODEL
fprintf(best, "shapemodel\n");
fprintf(best, "\tid %s\n", source[i].n);
fprintf(best, "\ttype %d\n", source[i].type);
fprintf(best, "\ts_center_x %.6lf\n", source[i].C.x);
fprintf(best, "\ts_center_y %.6lf\n", source[i].C.y);
fprintf(best, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
fprintf(best, "\ts_sigx %.6lf\n", source[i].E.a);
fprintf(best, "\ts_sigy %.6lf\n", source[i].E.b);
fprintf(best, "\ts_eps %.6lf\n", source[i].eps);
fprintf(best, "\tmag %.6lf\n", source[i].mag);
fprintf(best, "\ttype %d\n", source[i].type);
fprintf(best, "\tindex %.6lf\n", source[i].var1);
fprintf(best, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
fprintf(best, "\ts_mag %.6lf\n", source[i].mag);
fprintf(best, "\ts_index %.6lf\n", source[i].var1);
if(source[i].type==6) // double Sersic
{
fprintf(best, "\ts_sigx2 %.6lf\n", source[i].E2.a);
fprintf(best, "\ts_sigy2 %.6lf\n", source[i].E2.b);
fprintf(best, "\ts_eps2 %.6lf\n", source[i].eps2);
fprintf(best, "\ts_angle2 %.6lf\n", source[i].E2.theta * RTD);
fprintf(best, "\ts_mag2 %.6lf\n", source[i].mag2);
fprintf(best, "\ts_index2 %.6lf\n", source[i].var2);
}
fprintf(best, "\tend\n");
}
......@@ -572,16 +581,25 @@ void o_print_res(double chi0, double evidence)
// SHAPE MODEL
fprintf(besto, "shapemodel\n");
fprintf(besto, "\tid %s\n", source[i].n);
fprintf(besto, "\ttype %d\n", source[i].type);
fprintf(besto, "\ts_center_x %.6lf\n", source[i].C.x);
fprintf(besto, "\ts_center_y %.6lf\n", source[i].C.y);
fprintf(besto, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
fprintf(besto, "\ts_sigx %.6lf\n", source[i].E.a);
fprintf(besto, "\ts_sigy %.6lf\n", source[i].E.b);
fprintf(besto, "\ts_eps %.6lf\n", source[i].eps);
fprintf(besto, "\tmag %.6lf\n", source[i].mag);
fprintf(besto, "\ttype %d\n", source[i].type);
fprintf(besto, "\tindex %.6lf\n", source[i].var1);
fprintf(besto, "\tend\n");
fprintf(besto, "\ts_angle %.6lf\n", source[i].E.theta * RTD);
fprintf(besto, "\ts_mag %.6lf\n", source[i].mag);
fprintf(besto, "\ts_index %.6lf\n", source[i].var1);
if(source[i].type==6) //double Sersic
{
fprintf(besto, "\ts_sigx2 %.6lf\n", source[i].E2.a);
fprintf(besto, "\ts_sigy2 %.6lf\n", source[i].E2.b);
fprintf(besto, "\ts_eps2 %.6lf\n", source[i].eps2);
fprintf(besto, "\ts_angle2 %.6lf\n", source[i].E2.theta * RTD);
fprintf(besto, "\ts_mag2 %.6lf\n", source[i].mag2);
fprintf(besto, "\ts_index2 %.6lf\n", source[i].var2);
fprintf(besto, "\tend\n");
}
// SHAPE LIMITS
fprintf(besto, "shapelimit\n");
......@@ -589,16 +607,29 @@ void o_print_res(double chi0, double evidence)
fprintf(besto, "\ts_center_x %d %.6lf %.6lf\n", sblock[i][SCX], smin[i].C.x, smax[i].C.x);
if( sblock[i][SCY] )
fprintf(besto, "\ts_center_y %d %.6lf %.6lf\n", sblock[i][SCY], smin[i].C.y, smax[i].C.y);
if( sblock[i][STHETA] )
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][SEPS] )
fprintf(besto, "\ts_eps %d %.6lf %.6lf\n", sblock[i][SEPS], smin[i].eps, smax[i].eps);
if( sblock[i][STHETA] )
fprintf(besto, "\ts_angle %d %.6lf %.6lf\n", sblock[i][STHETA], smin[i].E.theta * RTD, smax[i].E.theta * RTD);
if( sblock[i][SFLUX] )
fprintf(besto, "\tmag %d %.6lf %.6lf\n", sblock[i][SFLUX], smin[i].mag, smax[i].mag);
fprintf(besto, "\ts_mag %d %.6lf %.6lf\n", sblock[i][SFLUX], smin[i].mag, smax[i].mag);
if( sblock[i][SINDEX] )
fprintf(besto, "\tindex %d %.6lf %.6lf\n", sblock[i][SINDEX], smin[i].var1, smax[i].var1);
fprintf(besto, "\ts_index %d %.6lf %.6lf\n", sblock[i][SINDEX], smin[i].var1, smax[i].var1);
if(source[i].type==6) //double Sersic
{
if( sblock[i][SA2] )
fprintf(besto, "\ts_sigx2 %d %.6lf %.6lf\n", sblock[i][SA2], smin[i].E2.a, smax[i].E2.a);
if( sblock[i][SEPS2] )
fprintf(besto, "\ts_eps2 %d %.6lf %.6lf\n", sblock[i][SEPS2], smin[i].eps2, smax[i].eps2);
if( sblock[i][STHETA2] )
fprintf(besto, "\ts_angle2 %d %.6lf %.6lf\n", sblock[i][STHETA2], smin[i].E2.theta * RTD, smax[i].E2.theta * RTD);
if( sblock[i][SFLUX2] )
fprintf(besto, "\ts_mag2 %d %.6lf %.6lf\n", sblock[i][SFLUX2], smin[i].mag2, smax[i].mag2);
if( sblock[i][SINDEX2] )
fprintf(besto, "\ts_index2 %d %.6lf %.6lf\n", sblock[i][SINDEX2], smin[i].var2, smax[i].var2);
}
fprintf(besto, "\tend\n");
}
......
......@@ -274,6 +274,21 @@ int rescaleParam(int id, int type, int ipx, double *pval)
case(SFLUX):
val = prior(sblock[id][SFLUX], val, smin[id].mag, smax[id].mag);
break;
case(SA2):
val = prior(sblock[id][SA2], val, smin[id].E2.a, smax[id].E2.a);
break;
case(SEPS2):
val = prior(sblock[id][SEPS2], val, smin[id].eps2, smax[id].eps2);
break;
case(STHETA2):
val = prior(sblock[id][STHETA2], val, smin[id].E2.theta, smax[id].E2.theta);
break;
case(SINDEX2):
val = prior(sblock[id][SINDEX2], val, smin[id].var2, smax[id].var2);
break;
case(SFLUX2):
val = prior(sblock[id][SFLUX2], val, smin[id].mag2, smax[id].mag2);
break;
}
}
......@@ -554,7 +569,7 @@ 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++ )
for ( ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[k][ipx] != 0 )
valid *= rescaleParam(k, SOURCES, ipx, &cube[ipar++]);
......
......@@ -486,7 +486,7 @@ int bayesHeader()
extern struct g_source S;
for ( i = 0; i < S.ns; i++ )
for ( ipx = SCX; ipx <= SFLUX; ipx++ )
for ( ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[i][ipx] != 0 )
{
nDim++;
......
......@@ -147,7 +147,7 @@ void o_set_limit_bayes(int nParam, long int nVal, double **array,
extern struct g_source S;
for ( i = 0; i < S.ns; i++ )
for ( ipx = SCX; ipx <= SFLUX; ipx++ )
for ( ipx = SCX; ipx <= SFLUX2; ipx++ )
if ( sblock[i][ipx] != 0 )
{
sblock[i][ipx] = prior;
......
......@@ -102,6 +102,21 @@ void o_set_lmax(int i, int ipx, double x)
case(SINDEX):
smax[i].var1 = x;
break;
case(SA2):
smax[i].E2.a = x;
break;
case(SEPS2):
smax[i].eps2 = x;
break;
case(STHETA2):
smax[i].E2.theta = x;
break;
case(SFLUX2):
smax[i].mag2 = x;
break;
case(SINDEX2):
smax[i].var2 = x;
break;
case(VFCX):
vfmax.C.x = x;
break;
......
......@@ -103,6 +103,21 @@ void o_set_lmin(int i, int ipx, double x)
case(SINDEX):
smin[i].var1 = x;
break;
case(SA2):
smin[i].E2.a = x;
break;
case(SEPS2):
smin[i].eps2 = x;
break;
case(STHETA2):
smin[i].E2.theta = x;
break;
case(SFLUX2):
smin[i].mag2 = x;
break;
case(SINDEX2):
smin[i].var2 = x;
break;
case(VFCX):
vfmin.C.x = x;
break;
......
......@@ -38,6 +38,21 @@ void o_set_source(struct galaxie *source, int ipx, double val)
case(SFLUX):
source->mag = val;
break;
case(SA2):
source->E2.a = val;
break;
case(SEPS2):
source->eps2 = val;
break;
case(STHETA2):
source->E2.theta = val;
break;
case(SINDEX2):
source->var2 = val;
break;
case(SFLUX2):
source->mag2 = val;
break;
}
}
......@@ -68,6 +83,21 @@ double o_get_source(struct galaxie *source, int ipx)
case(SFLUX):
val = source->mag;
break;
case(SA2):
val = source->E2.a;
break;
case(SEPS2):
val = source->eps2;
break;
case(STHETA2):
val = source->E2.theta;
break;
case(SINDEX2):
val = source->var2;
break;
case(SFLUX2):
val = source->mag2;
break;
}
return val;
......
......@@ -87,13 +87,53 @@ void r_shapelimit(FILE *IN, FILE *OUT, long int i)
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SFLUX],
smin[i].mag, smax[i].mag);
}
else if (!strcmp(second,"index"))
else if (!strcmp(second,"s_index") ||
!strcmp(second,"index"))
{
sscanf(third, "%d %lf %lf", &sblock[i][SINDEX],
&smin[i].var1, &smax[i].var1);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SINDEX],
smin[i].var1, smax[i].var1);
}
else if (!strcmp(second,"s_sigx2"))
{
sscanf(third, "%d %lf %lf", &sblock[i][SA2],
&smin[i].E2.a, &smax[i].E2.a);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SA2],
smin[i].E2.a, smax[i].E2.a);
}
else if (!strcmp(second,"s_eps2"))
{
sscanf(third, "%d %lf %lf", &sblock[i][SEPS2],
&smin[i].eps2, &smax[i].eps2);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SEPS2],
smin[i].eps2, smax[i].eps2);
}
else if (!strcmp(second,"s_angle2"))
{
sscanf(third, "%d %lf %lf", &sblock[i][STHETA2],
&smin[i].E2.theta, &smax[i].E2.theta);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][STHETA2],
smin[i].E2.theta, smax[i].E2.theta);
smin[i].E2.theta *= DTR;
smax[i].E2.theta *= DTR;
}
else if (!strcmp(second,"s_mag2"))
{
sscanf(third, "%d %lf %lf", &sblock[i][SFLUX2],
&smin[i].mag2, &smax[i].mag2);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SFLUX2],
smin[i].mag2, smax[i].mag2);
}
else if (!strcmp(second,"s_index2") ||
!strcmp(second,"index2"))
{
sscanf(third, "%d %lf %lf", &sblock[i][SINDEX2],
&smin[i].var2, &smax[i].var2);
fprintf(OUT, "\t%s\t%d %lf %lf\n", second, sblock[i][SINDEX2],
smin[i].var2, smax[i].var2);
}
// Read the next line
fmot(IN, second);
......
......@@ -26,14 +26,18 @@ void r_shapemodel(FILE *IN,FILE *OUT, long int i)
sshape->C.x = sshape->C.y = 0.;
sshape->E.a = sshape->E.b = sshape->E.theta = 0.;
sshape->E2.a = sshape->E2.b = sshape->E2.theta = 0.;
sshape->eps = -1; // no ellipticity defined originally
sshape->eps2 = -1; // no ellipticity defined originally
sshape->mag = 0;
sshape->mag2 = 0;
sshape->z = 0;
sshape->dl0s = sshape->dos = sshape->dr = -1;
sshape->I0 = 50;
sshape->c = 'g';
sshape->type = 3; // gaussian profile (default)
sshape->var1 = 4; // Sersic index (default: De Vaucouleur)
sshape->var2 = 4; // Sersic index (default: De Vaucouleur)
sprintf(sshape->n, "S%ld", i); // source name (default: S%d)
fmot(IN,second);
......@@ -83,7 +87,8 @@ void r_shapemodel(FILE *IN,FILE *OUT, long int i)
sscanf(third,"%lf",&sshape->mag);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->mag);
}
else if (!strcmp(second,"index"))
else if (!strcmp(second,"s_index") ||
!strcmp(second,"index"))
{
sscanf(third,"%lf",&sshape->var1);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->var1);
......@@ -103,6 +108,35 @@ void r_shapemodel(FILE *IN,FILE *OUT, long int i)
sscanf(third,"%s",sshape->n);
fprintf(OUT,"\t%s\t%s\n",second,sshape->n);
}
else if (!strcmp(second,"s_sigx2") )
{
sscanf(third,"%lf",&sshape->E2.a);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->E2.a);
}
else if (!strcmp(second,"s_eps2") )
{
sscanf(third,"%lf",&sshape->eps2);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->eps2);
}
else if (!strcmp(second,"s_angle2") )
{
sscanf(third,"%lf",&sshape->E2.theta);
sshape->E2.theta*=DTR;
fprintf(OUT,"\t%s\t%lf\n",second,sshape->E.theta);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->E2.theta);
}
else if (!strcmp(second,"s_mag2") ||
!strcmp(second, "mag2"))
{
sscanf(third,"%lf",&sshape->mag2);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->mag2);
}
else if (!strcmp(second,"s_index2") ||
!strcmp(second,"index2"))
{
sscanf(third,"%lf",&sshape->var2);
fprintf(OUT,"\t%s\t%lf\n",second,sshape->var2);
}
// Read the next line
fmot(IN,second);
......
......@@ -290,12 +290,19 @@ void setBayesModel( long int iVal, long int nVal, double **array)
if( pch != NULL )
source[i].C.x = source[i].C.y = 0.;
for( ipx = SCX ; ipx <= SFLUX; ipx++ )
for( ipx = SCX ; ipx <= SFLUX2; ipx++ )
if( sblock[i][ipx] != 0 )
o_set_source( &source[i], ipx, method(nVal, iParam++, array));
if((sblock[i][SA] !=0) & (source[i].E.a<0))
source[i].E.a=-source[i].E.a;
if((sblock[i][SA2] !=0) & (source[i].E2.a<0))
source[i].E2.a=-source[i].E2.a;
if((sblock[i][SEPS] != 0 )|| (sblock[i][SA] !=0))
source[i].E.b = source[i].E.a * (1. - source[i].eps) / (1. + source[i].eps);
if((sblock[i][SEPS2] != 0 )|| (sblock[i][SA2] !=0))
source[i].E2.b = source[i].E2.a * (1. - source[i].eps2) / (1. + source[i].eps2);
}
}
......@@ -843,11 +850,13 @@ static void convertArray(long int nVal, double **array )
// Sourcelimit parameters
if( M.iclean == 2 )
for ( i = 0; i < S.ns; i++ )
for( ipx = SCX; ipx <= SFLUX; ipx ++ )
for( ipx = SCX; ipx <= SFLUX2; ipx ++ )
if( sblock[i][ipx] != 0 )
{
if( ipx == STHETA )
array[nParam][iVal] = array[nParam][iVal] * DTR;
if( ipx == STHETA2 )
array[nParam][iVal] = array[nParam][iVal] * DTR;
nParam ++;
}
}
......
Markdown is supported
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