Skip to content

GitLab

Commit 5f8198fe authored by Johan Richard's avatar Johan Richard
Browse files

Merged master

parents a2ca462b 1993aedf
Hide whitespace changes
Inline Side-by-side
include/dimension.h
View file @ 5f8198fe
......@@ -17,7 +17,7 @@
#define NLMAX 4000 // maximum number of clumps in the lens[] array
#define NIMAX 8 /* maximum images per family */
#define NFMAX 50000 /* maximum number of families */
#define NPAMAX 35 // number of free parameters (see #define in structure.h)
#define NPAMAX 40 // number of free parameters (see #define in structure.h)
#define NTMAX 1024
#define NPOINT 100000 /* Number of contour points in cleanlens mode*/
#define NPARMAX 50
......@@ -39,7 +39,7 @@
#define IDPARAM1 1 // column of the 1st physical parameter in array (cf readBayesModel.c)
#define LINESIZE 16000 // size of a line in bayes.dat
#define FILENAME_SIZE 50 // size of a filename in .par file
#define FILENAME_SIZE 100 // size of a filename in .par file
#define DISTCOSMO2_ZPREC 1e-3 // precision below which 2 redshifts are considered the same
#define LHOOD_BUFNAME "restart.dat" // name of the restart file
......
include/fonction.h
View file @ 5f8198fe
......@@ -69,6 +69,7 @@ double comp_chi_osv(double *vect);
double comp_chi_osv(double *vect);
void comp_dchi_osv(double *vect, double *dvect);
int comp_desc(double A, double B);
void cl_connectdots(int);
void convertXY( double *x, double *y, int iref, double ra, double dec );
void copyright();
void copyTriplet(struct triplet *in, struct triplet *out);
......@@ -94,6 +95,7 @@ complex dcpx(complex c1, complex c2);
complex dcpxflt(complex c, double f);
double determinant(const struct point *A, const struct point *B, const struct point *C);
double d_gauss(double sig, int *idum);
void d_rbf(double **im, int nx, int ny, double scale);
double deltaVF(double zc);
struct ellipse diag(double a, double b, double c);
double diff_mag(struct galaxie *arclet, struct point *guess);
......@@ -455,6 +457,7 @@ void w_prop(int nprop, char propfile[]);
void wr_mass(char *name, double **map_xx, double **map_yy);
void wr_pot(char *name, double **map);
void w_sicat(struct galaxie ima[NFMAX][NIMAX], struct galaxie ss[NFMAX]);
void wrf_from_header(char *filein, char *filename, int nx, int ny, double **tabpix);
double zero(double c1, double c2, double (*f)(double));
double zero_t(double c1, double c2, double (*f)(double));
void zonemult();
......
include/structure.h
View file @ 5f8198fe
......@@ -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
......@@ -234,6 +239,7 @@ struct g_mode
int iclean;
double zclean;
int flux; //conserve surface brightness ? 1=True 0=False
int output; //writes models and residuals from shapemodel
/* center of multiple images file */
char centerfile[FILENAME_SIZE];
/* bayesys restart flag */
......@@ -583,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;
......@@ -596,6 +604,7 @@ struct galaxie
double dr; // ratio dl0s/dos
double q;
double eps;
double eps2;
double tau;
double dis;
double A;
......
liblt/wfits.c
View file @ 5f8198fe
......@@ -817,3 +817,54 @@ int wrf_cube_fits_abs(char *filename, double ***cube, int nx,int ny, int nz,
return 0;
}
/* write output image based on known header */
void wrf_from_header(char *filein, char *filename, int nx, int ny, double **tabpix)
{
fitsfile *infptr;
fitsfile *outfptr;
int status;
int nkeys;
int fpixel;
int nelements;
float *p_ima;
int ii,jj,k;
int bitpix=FLOAT_IMG;
long naxes[2] = {nx,ny};
status = 0; /* initialize status before calling fitsio routines */
if( fits_open_file(&infptr, filein, READONLY, &status) )
printerror( status );
fpixel = 1;
nelements = nx*ny;
p_ima=(float *)malloc((unsigned) (nelements)*sizeof(float));
k=0;
for (jj=0;jj<ny;jj++)
for (ii=0;ii<nx;ii++)
{
p_ima[k]=(float)tabpix[jj][ii];
k++;
}
remove(filename); /* Delete old file if it already exists */
status = CREATE_DISK_FILE; /* initialize status before calling fitsio routines */
if (fits_create_file(&outfptr, filename, &status)) /* create new FITS file */
printerror( status ); /* call printerror if error occurs */
if ( fits_copy_header(infptr, outfptr, &status) )
printerror( status );
// if ( fits_create_img(outfptr, bitpix, 2, naxes, &status) )
// printerror( status );
if ( fits_write_img(outfptr, TFLOAT, fpixel, nelements, p_ima, &status) )
printerror( status );
if ( fits_close_file(infptr, &status) ) /* close the file */
printerror( status );
if ( fits_close_file(outfptr, &status) )
printerror(status);
free(p_ima);
}
src/bayesapp.c
View file @ 5f8198fe
......@@ -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;
......@@ -787,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;
......@@ -1164,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);
}
......
src/checkpar.c
View file @ 5f8198fe
......@@ -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,53 @@ 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;
sblock[i][STHETA]=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;
}
}
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;
......
src/criticnew.c
View file @ 5f8198fe
......@@ -33,7 +33,6 @@
static void snake(int verbose);
static void marchingSquares(int verbose);
static void cl_connectdots(int verbose);
void criticnew(int verbose)
{
......@@ -47,8 +46,6 @@ void criticnew(int verbose)
else
{
marchingSquares(verbose);
w_critic();
cl_connectdots(verbose);
}
}
......@@ -571,7 +568,7 @@ static void marchingSquares(int verbose)
}
}
static void cl_connectdots(int verbose)
void cl_connectdots(int verbose)
{
extern struct g_cline CL;
extern struct biline radial[], tangent[];
......
src/d_bruiter_omp.c
View file @ 5f8198fe
......@@ -5,6 +5,7 @@
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <math.h>
#include <fonction.h>
// Append a sky background to the pixels of image **z
......@@ -13,7 +14,7 @@ void d_bruiter_omp(double **z, int nx, int ny)
const extern struct g_observ O;
//at this point we mustn't be in parallel
check_not_in_parallel();
check_not_in_parallel("");
int i;
#pragma omp parallel for schedule(static)
......
src/d_profil.c
View file @ 5f8198fe
......@@ -26,11 +26,16 @@
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
//fprintf(stderr,"%f\n",gal->E.a*(1-gal->eps)/(1+gal->eps)-gal->E.b);
if (gal->c == 's')
{
......@@ -115,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);
......
src/d_seeing_omp.c
View file @ 5f8198fe
......@@ -20,6 +20,7 @@
#include <stdlib.h>
#include "structure.h"
#include "lt.h"
#include "fonction.h"
struct cpa2d; //predifiniton of cpa2d (see below)
......
src/do_itos.c
View file @ 5f8198fe
......@@ -73,11 +73,11 @@ void do_itos(double **im, struct pixlist *pl, int npl, double dlsds, double z
if (is >= 0 && js >= 0 && is < ps.ny && js < ps.nx)
{
imult[is][js] += 1;
source[is][js] = ech*((imult[is][js] - 1) * source[is][js] +
pl[k].flux)
source[is][js] = ((imult[is][js] - 1) * source[is][js] +
pl[k].flux*ech)
/ ((double)imult[is][js]);
erreur[is][js] = ech*((imult[is][js] - 1) * erreur[is][js] +
pl[k].flux * pl[k].flux)
erreur[is][js] = ((imult[is][js] - 1) * erreur[is][js] +
pl[k].flux * pl[k].flux*ech*ech)
/ ((double)imult[is][js]);
}
......@@ -113,11 +113,11 @@ void do_itos(double **im, struct pixlist *pl, int npl, double dlsds, double z
if ((is >= 0) && (js >= 0) && (is < ps.ny) && (js < ps.nx))
{
imult[is][js] += 1;
source[is][js] = ech*((imult[is][js] - 1) * source[is][js] +
pl[k].flux)
source[is][js] = ((imult[is][js] - 1) * source[is][js] +
pl[k].flux*ech)
/ ((double)imult[is][js]);
erreur[is][js] = ech*((imult[is][js] - 1) * erreur[is][js] +
pl[k].flux * pl[k].flux)
erreur[is][js] = ((imult[is][js] - 1) * erreur[is][js] +
pl[k].flux * pl[k].flux*ech*ech)
/ ((double)imult[is][js]);
};
A.x += imFrame.pixelx / imFrame.ech * ((double) j);
......
src/interpol.c
View file @ 5f8198fe
......@@ -42,6 +42,9 @@ double interpol(double xx, const double *fxx, const double *fyy, const double *
double bilinear(double **map, double x, double y, int xmax, int ymax)
{
int x1,y1;
//Difference of wcs2pix reference (starts at 1,1) and map reference (starts at pixel 0,0)
x=x-1;
y=y-1;
x1=(int)(x-0.5);
y1=(int)(y-0.5);
if(x1<0) x1=0;
......
src/main.c
View file @ 5f8198fe
......@@ -173,6 +173,8 @@ int main(int argc, char *argv[])
if ( G.nlens > 0 )
{
criticnew(1);
w_critic();
};
// if there is only 1 clump
......@@ -181,7 +183,10 @@ int main(int argc, char *argv[])
/* calcul des zones d'images multiples */
if ( CL.zone != 0 )
{
cl_connectdots(1);
zonemult();
}
}
if ( M.marker != 0 )
......
src/o_chires.c
View file @ 5f8198fe
......@@ -211,6 +211,15 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
M.verbose = verbose_save;
// mod jrichard output shapemodel
if(M.output)
{
fprintf(stderr, "Write output model with same header as input - smodel.fits\n");
wrf_from_header(imFrame.pixfile,"smodel.fits", imFrame.nx, imFrame.ny, ero);
}
if (wo != NULL)
{
// mod bclement shapemodel_contour
......@@ -243,6 +252,20 @@ int o_chi_lhood0(double *chi_ext, double *lhood0_ext, double *np_b0)
}
}
// mod jrichard output shapemodel
if(M.output)
{
fprintf(stderr, "Write output residuals with same header as input - residuals.fits\n");
//invert residuals to produce data - model instead of model - data (within contours)
for( i = 0; i < imFrame.ny; i++ )
for( j = 0; j < imFrame.nx; j++)
ero[i][j]=-ero[i][j];
wrf_from_header(imFrame.pixfile,"residuals.fits", imFrame.nx, imFrame.ny, ero);
}
// restore relative source positions
if( I.n_mult > 0 )
for( i = 0; i < S.ns; i++ )
......@@ -1854,3 +1877,4 @@ static void srcfitLinFit(double *chi2, double *lh0)
#endif
}
src/o_global.c
View file @ 5f8198fe
......@@ -500,6 +500,12 @@ void readConstraints()
O.bin*O.bin*imFrame.nx * imFrame.ny * sizeof(struct pixlist));
keep_cl_extended(imo_pl, imo_npl, imo_pl_ext, &imo_npl_ext);
}
else
{
imo_pl_ext = (struct pixlist*)malloc((unsigned int)
imFrame.nx * imFrame.ny * sizeof(struct pixlist));
keep_cl(imo, imo_pl_ext, &imo_npl_ext);
}
}
// TODO : remove duplicate pixels in list when overlapping contours occur
......@@ -943,7 +949,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++;
......
src/o_print_res.c
View file @ 5f8198fe
......@@ -61,7 +61,7 @@ void o_print_res(double chi0, double evidence)
extern double z_dlsds;
extern double chip, chix, chiy, chis, chil, chi_vel,chi_mass; //I added chi_vel and chi_mass TV
extern double **map_p, **map_axx, **map_ayy;
extern struct g_pixel ps, imFrame;
extern struct g_pixel ps, imFrame, wFrame;
int i, j;
FILE *best, *besto;
char limages[ZMBOUND][IDSIZE];
......@@ -249,6 +249,8 @@ void o_print_res(double chi0, double evidence)
fprintf(best, "\tcleanset %d %f\n", M.iclean, M.zclean);
if(strcmp(imFrame.pixfile, ""))
fprintf(best, "\timframe %d %s\n", imFrame.format, imFrame.pixfile);
if(strcmp(wFrame.pixfile, ""))
fprintf(best, "\twframe %d %s\n", wFrame.format, wFrame.pixfile);
if(strcmp(ps.pixfile, ""))
fprintf(best, "\tsframe %s\n", ps.pixfile);
if(strcmp(M.centerfile, ""))
......@@ -283,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");
}
......@@ -533,7 +544,9 @@ void o_print_res(double chi0, double evidence)
fprintf(besto, "cleanlens\n");
fprintf(besto, "\tcleanset %d %f\n", M.iclean, M.zclean);
if(strcmp(imFrame.pixfile, ""))
fprintf(besto, "\timFrame %d %s\n", imFrame.format, imFrame.pixfile);
fprintf(besto, "\timframe %d %s\n", imFrame.format, imFrame.pixfile);
if(strcmp(wFrame.pixfile, ""))
fprintf(besto, "\twframe %d %s\n", wFrame.format, wFrame.pixfile);
if(strcmp(ps.pixfile, ""))
fprintf(besto, "\tsframe %s\n", ps.pixfile);
if(strcmp(M.centerfile, ""))
......@@ -568,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");
......@@ -585,18 +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][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][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);