Commit 96c6e08a authored by Claudio Gheller's avatar Claudio Gheller
Browse files

w-stacking completed

parent f43f86ed

Makefile

+1 −1
Original line number Diff line number Diff line
@@ -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++

phase_correction.cu

+23 −18
Original line number Diff line number Diff line
@@ -9,16 +9,19 @@ 
#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++)

            {
@@ -26,25 +29,21 @@ void phase_correction(double* gridss, double* image_real, double* image_imag, in 
		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;

		    } 





		    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;

		    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;

	}







}

w-stacking-fftw.c

+37 −6
Original line number Diff line number Diff line
@@ -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);

w-stacking.h

+4 −2
Original line number Diff line number Diff line
@@ -43,6 +43,8 @@ void phase_correction( 
     int,

     int,

     int,

     int);



     int,

     double,

     double,

     double);

#endif