diff --git a/w-stacking-fftw.c b/w-stacking-fftw.c
index 5469014518e9e03a3d65ee9461bccb9c58cebfc1..6e00e73eb4ede8a460fa899779e706a2a6ec180b 100644
--- a/w-stacking-fftw.c
+++ b/w-stacking-fftw.c
@@ -1,6 +1,9 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#ifdef FITSIO
+#include "/m100/home/userexternal/ederubei/cfitsio-3.49/fitsio.h"
+#endif
#ifdef USE_MPI
#include <mpi.h>
#ifdef USE_FFTW
@@ -44,6 +47,7 @@ int main(int argc, char * argv[])
int rank;
int size;
+
// Define main filenames
FILE * pFile;
FILE * pFile1;
@@ -51,7 +55,7 @@ int main(int argc, char * argv[])
FILE * pFileimg;
// Global filename to be composed
char filename[1000];
-
+
// MS paths
char datapath[900];
char datapath_multi[NFILES][900];
@@ -118,11 +122,29 @@ int main(int argc, char * argv[])
int num_threads;
// Resolution
- double dx = 1.0/(double)grid_size_x;
- double dw = 1.0/(double)num_w_planes;
+ double dx = 0.5/(double)grid_size_x;
+ double dw = 0.5/(double)num_w_planes;
// Half support size
double w_supporth = (double)((w_support-1)/2)*dx;
+
+ // Initialize FITS image parameters
+ fitsfile *fptreal;
+ fitsfile *fptrimg;
+ int status;
+ long nelements;
+ long fpixel, lpixel;
+ char testfitsreal[FILENAMELENGTH] = "parallel_np2_real.fits";
+ char testfitsimag[FILENAMELENGTH] = "parallel_np2_img.fits";
+
+ long naxis = 2;
+ long naxes[2] = { grid_size_x, grid_size_y };
+
+ nelements = naxes[0] * naxes[1];
+
+
+
+
// Internal profiling parameters
clock_t start, end, start0, startk, endk;
double setup_time, process_time, mpi_time, fftw_time, tot_time, kernel_time, reduce_time, compose_time, phase_time;
@@ -201,9 +223,9 @@ if(rank == 0){
// INPUT FILES (only the first ndatasets entries are used)
int ndatasets = 1;
- //strcpy(datapath_multi[0],"data/newgauss2noconj_t201806301100_SBL180.binMS/");
+ strcpy(datapath_multi[0],"/m100_scratch/userexternal/ederubei/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/");
//strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/newgauss4_t201806301100_SBL180.binMS/");
- strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/");
+ //strcpy(datapath_multi[0],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_146MHz.pre-cal.binMS/");
//strcpy(datapath_multi[1],"/m100_scratch/userexternal/cgheller/gridding/Lofar/L798046_SB244_uv.uncorr_130B27932t_134MHz.pre-cal.binMS/");
strcpy(datapath,datapath_multi[0]);
@@ -917,8 +939,8 @@ if(rank == 0){
clock_gettime(CLOCK_MONOTONIC, &begin);
start = clock();
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));
+ 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,resolution,wmin,wmax,num_threads);
@@ -931,10 +953,33 @@ if(rank == 0){
if(rank == 0)
{
- pFilereal = fopen (fftfile2,"wb");
- pFileimg = fopen (fftfile3,"wb");
- fclose(pFilereal);
- fclose(pFileimg);
+ #ifdef FITSIO
+ printf("REMOVING RESIDUAL FITS FILE\n");
+ remove(testfitsreal);
+ remove(testfitsimag);
+
+
+ printf("FITS CREATING\n");
+ status = 0;
+
+ fits_create_file(&fptrimg, testfitsimag, &status);
+ fits_create_img(fptrimg, DOUBLE_IMG, naxis, naxes, &status);
+// fits_write_img(fptrimg, TDOUBLE, fpixel, nelements, image_imag, &status);
+ fits_close_file(fptrimg, &status);
+
+ status = 0;
+
+ fits_create_file(&fptreal, testfitsreal, &status);
+ fits_create_img(fptreal, DOUBLE_IMG, naxis, naxes, &status);
+// fits_write_img(fptreal, TDOUBLE, fpixel, nelements, image_real, &status);
+ fits_close_file(fptreal, &status);
+ #endif
+
+// pFilereal = fopen (fftfile2,"wb");
+// pFileimg = fopen (fftfile3,"wb");
+
+// fclose(pFilereal);
+// fclose(pFileimg);
}
#ifdef USE_MPI
MPI_Barrier(MPI_COMM_WORLD);
@@ -947,19 +992,41 @@ if(rank == 0){
#endif
if(isector == rank)
{
- printf("%d writing\n",isector);
- pFilereal = fopen (fftfile2,"ab");
- pFileimg = fopen (fftfile3,"ab");
-
- long global_index = isector*(xaxis*yaxis)*sizeof(double);
+ #ifdef FITSIO
- fseek(pFilereal, global_index, SEEK_SET);
- fwrite(image_real, xaxis*yaxis, sizeof(double), pFilereal);
- fseek(pFileimg, global_index, SEEK_SET);
- fwrite(image_imag, xaxis*yaxis, sizeof(double), pFileimg);
-
- fclose(pFilereal);
- fclose(pFileimg);
+ printf("%d writing\n",isector);
+ long * fpixel = (long *) malloc(sizeof(long)*naxis);
+ long * lpixel = (long *) malloc(sizeof(long)*naxis);
+
+ fpixel[0] = 1;
+ lpixel[0] = xaxis;
+ fpixel[1] = isector*yaxis+1;
+ lpixel[1] = (isector+1)*yaxis;
+
+ status = 0;
+ fits_open_image(&fptreal, testfitsreal, READWRITE, &status);
+ fits_write_subset(fptreal, TDOUBLE, fpixel, lpixel, image_real, &status);
+ fits_close_file(fptreal, &status);
+
+ status = 0;
+ fits_open_image(&fptrimg, testfitsimag, READWRITE, &status);
+ fits_write_subset(fptrimg, TDOUBLE, fpixel, lpixel, image_imag, &status);
+ fits_close_file(fptrimg, &status);
+
+ #endif
+
+// pFilereal = fopen (fftfile2,"ab");
+// pFileimg = fopen (fftfile3,"ab");
+
+// long global_index = isector*(xaxis*yaxis)*sizeof(double);
+
+// fseek(pFilereal, global_index, SEEK_SET);
+// fwrite(image_real, xaxis*yaxis, sizeof(double), pFilereal);
+// fseek(pFileimg, global_index, SEEK_SET);
+// fwrite(image_imag, xaxis*yaxis, sizeof(double), pFileimg);
+
+// fclose(pFilereal);
+// fclose(pFileimg);
}
}
#ifdef USE_MPI