diff --git a/Makefile b/Makefile
index 80ba945444a7f245e0d947789cc977ce4b5d8cc0..0dcde72606c8d3260e7fbde2e13d0b5e9c0ee53e 100644
--- a/Makefile
+++ b/Makefile
@@ -10,7 +10,7 @@ OPT += -DONE_SIDE
# write the final image
OPT += -DWRITE_IMAGE
# perform w-stacking phase correction
-# OPT += PHASE_ON
+OPT += -DPHASE_ON
CC = gcc
CXX = g++
diff --git a/phase_correction.cu b/phase_correction.cu
index 9d1cc7e790b8aec7c5a552745719daf862018cff..77f5ad25a9197dca27a2ae3e6a567f3000fec583 100644
--- a/phase_correction.cu
+++ b/phase_correction.cu
@@ -9,42 +9,41 @@
#ifdef __CUDACC__
#endif
-void phase_correction(double* gridss, double* image_real, double* image_imag, int xaxis, int yaxis, int num_w_planes, int xaxistot, int yaxistot)
+void phase_correction(double* gridss, double* image_real, double* image_imag, int xaxis, int yaxis, int num_w_planes, int xaxistot, int yaxistot,
+ double resolution, double wmin, double wmax)
{
double dnum_w_planes = (double)num_w_planes;
double dxaxistot = (double)xaxistot;
double dyaxistot = (double)yaxistot;
+ double diagonal;
+ double dw = (wmax-wmin)/(double)num_w_planes;
+ double wterm = wmin+0.5*dw;
+ double dwnorm = dw/(wmax-wmin);
for (int iw=0; iw<num_w_planes; iw++)
{
- double wterm = (double)iw/dnum_w_planes;
-
for (int iv=0; iv<yaxis; iv++)
for (int iu=0; iu<xaxis; iu++)
{
- long index = 2*(iu+iv*xaxis+xaxis*yaxis*iw);
- long img_index = iu+iv*xaxis;
+ long index = 2*(iu+iv*xaxis+xaxis*yaxis*iw);
+ long img_index = iu+iv*xaxis;
#ifdef PHASE_ON
+ if (num_w_planes > 1)
+ {
double xcoord = (double)(iu-xaxistot/2);
if(xcoord < 0.0)xcoord = (double)(iu+xaxistot/2);
+ xcoord = sin(xcoord*resolution);
double ycoord = (double)(iv-yaxistot/2);
if(ycoord < 0.0)ycoord = (double)(iv+yaxistot/2);
- xcoord = xcoord/dxaxistot;
- ycoord = ycoord/dyaxistot;
+ ycoord = sin(ycoord*resolution);
double efact;
double preal, pimag;
double radius2 = (xcoord*xcoord+ycoord*ycoord);
- if(xcoord <= 1.0)
- {
- preal = cos(2.0*PI*wterm*(sqrt(1-radius2)-1.0));
- pimag = sin(2.0*PI*wterm*(sqrt(1-radius2)-1.0));
- } else {
- preal = cos(-2.0*PI*wterm*(sqrt(radius2-1.0)-1));
- pimag = 0.0;
- }
+ preal = cos(2.0*PI*wterm*(sqrt(1-radius2)-1.0));
+ pimag = sin(2.0*PI*wterm*(sqrt(1-radius2)-1.0));
double p,q,r,s;
p = gridss[index];
@@ -53,15 +52,21 @@ void phase_correction(double* gridss, double* image_real, double* image_imag, in
s = pimag;
//printf("%d %d %d %ld %ld\n",iu,iv,iw,index,img_index);
- image_real[img_index] += p*r-q*s;
- image_imag[img_index] += p*s+q*r;
-#else
+ image_real[img_index] += (p*r-q*s)*dwnorm*sqrt(1-radius2);
+ image_imag[img_index] += (p*s+q*r)*dwnorm*sqrt(1-radius2);
+ } else {
image_real[img_index] += gridss[index];
image_imag[img_index] += gridss[index+1];
+ }
+#else
+ image_real[img_index] += gridss[index];
+ image_imag[img_index] += gridss[index+1];
#endif
}
+ wterm += dw;
}
+
}
diff --git a/w-stacking-fftw.c b/w-stacking-fftw.c
index 4d8130f87c6b66cef76611158d7b80076a3744ab..55872f61fc2935d8fbcba3e607a696d7b7397042 100644
--- a/w-stacking-fftw.c
+++ b/w-stacking-fftw.c
@@ -15,7 +15,7 @@
#define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
#define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
#define NOVERBOSE
-#define NFILES 10
+#define NFILES 100
// Linked List set-up
struct sectorlist {
@@ -91,15 +91,16 @@ int main(int argc, char * argv[])
double uvmax;
double wmin;
double wmax;
+ double resolution;
// MESH SIZE
int grid_size_x = 2048;
int grid_size_y = 2048;
- int local_grid_size_x;// = 100;
- int local_grid_size_y;// = 100;
+ int local_grid_size_x;// = 8;
+ int local_grid_size_y;// = 8;
int xaxis;
int yaxis;
- int num_w_planes = 1;
+ int num_w_planes = 8;
// DAV: the corresponding KernelLen is calculated within the wstack function. It can be anyway hardcoded for optimization
int w_support = 7;
int num_threads;// = 4;
@@ -176,6 +177,7 @@ int main(int argc, char * argv[])
fscanf(pFile,"%lf",&wmax);
fclose(pFile);
+
// WATCH THIS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
int nsub = 1000;
//int nsub = 10;
@@ -357,13 +359,41 @@ int main(int argc, char * argv[])
reduce_time1 = 0.0;
compose_time = 0.0;
compose_time1 = 0.0;
+
+ // MAIN LOOP OVER FILES
+ //
for (int ifiles=0; ifiles<ndatasets; ifiles++)
{
strcpy(filename,datapath_multi[ifiles]);
printf("Processing %s, %d of %d\n",filename,ifiles+1,ndatasets);
- strcat(filename,weightsfile);
+ // Read metadata
+ strcpy(filename,datapath);
+ strcat(filename,metafile);
+ pFile = fopen (filename,"r");
+ fscanf(pFile,"%ld",&Nmeasures);
+ fscanf(pFile,"%ld",&Nvis);
+ fscanf(pFile,"%ld",&freq_per_chan);
+ fscanf(pFile,"%ld",&polarisations);
+ fscanf(pFile,"%ld",&Ntimes);
+ fscanf(pFile,"%lf",&dt);
+ fscanf(pFile,"%lf",&thours);
+ fscanf(pFile,"%ld",&baselines);
+ fscanf(pFile,"%lf",&uvmin);
+ fscanf(pFile,"%lf",&uvmax);
+ fscanf(pFile,"%lf",&wmin);
+ fscanf(pFile,"%lf",&wmax);
+ fclose(pFile);
+
+ // calculate the resolution in radians
+ resolution = 1.0/MAX(abs(uvmin),abs(uvmax));
+ // calculate the resolution in arcsec
+ double resolution_asec = (3600.0*180.0)/MAX(abs(uvmin),abs(uvmax))/PI;
+ printf("RESOLUTION = %f rad, %f arcsec\n", resolution, resolution_asec);
+
+ strcpy(filename,datapath);
+ strcat(filename,weightsfile);
pFile = fopen (filename,"rb");
fseek (pFile,startrow*polarisations*sizeof(float),SEEK_SET);
fread(weights,(Nweights)*sizeof(float),1,pFile);
@@ -829,7 +859,8 @@ int main(int argc, char * argv[])
if(rank == 0)printf("PHASE CORRECTION\n");
double* image_real = (double*) calloc(xaxis*yaxis,sizeof(double));
double* image_imag = (double*) calloc(xaxis*yaxis,sizeof(double));
- phase_correction(gridss,image_real,image_imag,xaxis,yaxis,num_w_planes,grid_size_x,grid_size_y);
+
+ phase_correction(gridss,image_real,image_imag,xaxis,yaxis,num_w_planes,grid_size_x,grid_size_y,resolution,wmin,wmax);
end = clock();
clock_gettime(CLOCK_MONOTONIC, &finish);
diff --git a/w-stacking.h b/w-stacking.h
index 5c675e474d2174bd47b54b4dd98c2a7faf95fd4b..822b92c4f55e2f7469cfff9cf18fea78103e172d 100644
--- a/w-stacking.h
+++ b/w-stacking.h
@@ -43,6 +43,8 @@ void phase_correction(
int,
int,
int,
- int);
-
+ int,
+ double,
+ double,
+ double);
#endif