CFITSIO and parallel images writing implementation

Use CFITSIO to write .fits images as output other than .bin files

w-stacking-fftw.c

 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#ifdef FITSIO
+#include "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;
@@ -123,6 +127,19 @@ int main(int argc, char * argv[])
 	// Half support size
 	double w_supporth = (double)((w_support-1)/2)*dx;
 
+
+	// Initialize FITS image parameters
+	fitsfile *fptreal;
+	fitsfile *fptrimg;
+	int status;
+	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 };
+
+
+
 	// 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;
@@ -195,15 +212,14 @@ if(rank == 0){
 	// LOCAL grid size
 	xaxis = local_grid_size_x;
 	yaxis = local_grid_size_y;
-
 	clock_gettime(CLOCK_MONOTONIC, &begin);
 	start = clock();
 
         // 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]);
@@ -931,6 +947,26 @@ if(rank == 0){
 
         if(rank == 0)
         {
+            #ifdef FITSIO
+            printf("REMOVING RESIDUAL FITS FILE\n");
+            remove(testfitsreal);
+            remove(testfitsimag);
+
+
+	    printf("FITS CREATION\n");
+            status = 0;
+
+            fits_create_file(&fptrimg, testfitsimag, &status);
+            fits_create_img(fptrimg, DOUBLE_IMG, naxis, naxes, &status);
+            fits_close_file(fptrimg, &status);
+
+	    status = 0;
+
+            fits_create_file(&fptreal, testfitsreal, &status);
+	    fits_create_img(fptreal, DOUBLE_IMG, naxis, naxes, &status);
+	    fits_close_file(fptreal, &status);
+	    #endif
+
             pFilereal = fopen (fftfile2, "wb");
 	    pFileimg = fopen (fftfile3, "wb");
 	    fclose(pFilereal);
@@ -940,6 +976,49 @@ if(rank == 0){
 	MPI_Barrier(MPI_COMM_WORLD);
         #endif
         if(rank == 0)printf("WRITING IMAGE\n");
+        long * fpixel = (long *) malloc(sizeof(long)*naxis);
+        long * lpixel = (long *) malloc(sizeof(long)*naxis);
+
+        #ifdef USE_MPI
+        MPI_Barrier(MPI_COMM_WORLD);
+        #endif
+	#ifdef PARALLELIO
+        #ifdef FITSIO
+
+        fpixel[0] = 1;
+        fpixel[1] = rank*yaxis+1;
+        lpixel[0] = xaxis;
+        lpixel[1] = (rank+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 //FITSIO
+
+        pFilereal = fopen (fftfile2,"ab");
+        pFileimg = fopen (fftfile3,"ab");
+
+	long global_index = rank*(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
+	MPI_Barrier(MPI_COMM_WORLD);
+        #endif
+        #else
 	for (int isector=0; isector<size; isector++)
 	{
 	    #ifdef USE_MPI
@@ -947,7 +1026,27 @@ if(rank == 0){
             #endif
 	    if(isector == rank)
 	    {
+               #ifdef FITSIO
+
 	       printf("%d writing\n",isector);
+               
+               fpixel[0] = 1;
+               fpixel[1] = isector*yaxis+1;
+               lpixel[0] = xaxis;
+               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 //FITSIO
+
                pFilereal = fopen (fftfile2,"ab");
                pFileimg = fopen (fftfile3,"ab");
 
@@ -962,10 +1061,12 @@ if(rank == 0){
                fclose(pFileimg);
 	    }
 	}
+	#endif //PARALLELIO
 	#ifdef USE_MPI
 	MPI_Barrier(MPI_COMM_WORLD);
         #endif
 
+
 #endif //WRITE_IMAGE