Commit 6e356d39 authored by Claudio Gheller's avatar Claudio Gheller
Browse files

bug on the UVW normalization fixed

parent 7461a3fc

w-stacking-fftw.c

+20 −22
Original line number Diff line number Diff line
@@ -104,8 +104,8 @@ int main(int argc, char * argv[]) 
	double resolution;



        // Mesh related parameters: global size

	int grid_size_x = 2048;

	int grid_size_y = 2048;

	int grid_size_x = 4096;

	int grid_size_y = 4096;

	// Split Mesh size (auto-calculated)

	int local_grid_size_x;

	int local_grid_size_y;
@@ -219,7 +219,7 @@ if(rank == 0){ 


        // INPUT FILES (only the first ndatasets entries are used)

	int ndatasets = 1;

        strcpy(datapath_multi[0],"/homes/gheller/lofar2/gheller/binMS/ZW2_IFRQ_0-1-5.binMS/");

        strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/SKA3/ZW2_IFRQ_0-1-5.binMS/");



	strcpy(datapath,datapath_multi[0]);

	// Read metadata
@@ -244,8 +244,8 @@ if(rank == 0){ 
	// WATCH THIS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	int nsub = 1000;

	//int nsub = 10;

	printf("Subtracting last %d measurements\n",nsub);

  Nmeasures = Nmeasures-nsub;

	if(rank == 0)printf("Subtracting last %d measurements\n",nsub);

        Nmeasures = Nmeasures-nsub;

	Nvis = Nmeasures*freq_per_chan*polarisations;



  // calculate the coordinates of the center
@@ -302,7 +302,7 @@ if(rank == 0){ 
	strcat(filename,vfile);

	//printf("Reading %s\n",filename);



  pFile = fopen (filename,"rb");

        pFile = fopen (filename,"rb");

	fseek (pFile,startrow*sizeof(double),SEEK_SET);

	fread(vv,Nmeasures*sizeof(double),1,pFile);

	fclose(pFile);
@@ -340,12 +340,12 @@ if(rank == 0){ 
	double maxv_all;

	double maxw_all;

	#ifdef USE_MPI

	MPI_Allreduce(&minu,&minu_all,1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&minv,&minv_all,1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&minw,&minw_all,1, MPI_FLOAT, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&maxu,&maxu_all,1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);

	MPI_Allreduce(&maxv,&maxv_all,1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);

	MPI_Allreduce(&maxw,&maxw_all,1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);

	MPI_Allreduce(&minu,&minu_all,1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&minv,&minv_all,1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&minw,&minw_all,1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);

	MPI_Allreduce(&maxu,&maxu_all,1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

	MPI_Allreduce(&maxv,&maxv_all,1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

	MPI_Allreduce(&maxw,&maxw_all,1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

        #else

	minu_all = minu;

	minv_all = minv;
@@ -361,7 +361,6 @@ if(rank == 0){ 
	maxw = maxw_all;

        ming = ming-offset*ming;

        maxg = maxg+offset*maxg;



        for (long inorm=0; inorm<Nmeasures; inorm++)

	    {

	       uu[inorm] = (uu[inorm]+maxg)/(2.0*maxg);
@@ -483,7 +482,7 @@ if(rank == 0){ 
	for (int ifiles=0; ifiles<ndatasets; ifiles++)

	{

	strcpy(filename,datapath_multi[ifiles]);

        printf("Processing %s, %d of %d\n",filename,ifiles+1,ndatasets);

        if(rank == 0)printf("Processing %s, %d of %d\n",filename,ifiles+1,ndatasets);



        // Read metadata

        strcpy(filename,datapath);
@@ -507,7 +506,7 @@ if(rank == 0){ 
        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);

        if(rank == 0)printf("RESOLUTION = %f rad, %f arcsec\n", resolution, resolution_asec);



        strcpy(filename,datapath);

        strcat(filename,weightsfile);
@@ -543,7 +542,7 @@ if(rank == 0){ 
  float * visreals;

  float * visimgs;

  float * weightss;

	long isector;

  long isector;



  for (long isector_count=0; isector_count<nsectors; isector_count++)

      {
@@ -565,15 +564,14 @@ if(rank == 0){ 


	  // select data for this sector

        long icount = 0;

	      long ip = 0;

	      long inu = 0;

//CLAAAA

	       struct sectorlist * current;

	       current = sectorhead[isector];

	long ip = 0;

	long inu = 0;

	struct sectorlist * current;

	current = sectorhead[isector];







	       while (current->index != -1)

	while (current->index != -1)

          {

             long ilocal = current->index;

      	     //double vvh = vv[ilocal];