CFITSIO serial implementation

w-stacking-fftw.c

 #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;
@@ -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]);
@@ -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)
             {
+               #ifdef FITSIO
+
 	       printf("%d writing\n",isector);
-               pFilereal = fopen (fftfile2,"ab");
-               pFileimg = fopen (fftfile3,"ab");
+               long * fpixel = (long *) malloc(sizeof(long)*naxis);
+               long * lpixel = (long *) malloc(sizeof(long)*naxis);
                
-	       long global_index = isector*(xaxis*yaxis)*sizeof(double);
+               fpixel[0] = 1;
+               lpixel[0] = xaxis;
+               fpixel[1] = isector*yaxis+1;
+               lpixel[1] = (isector+1)*yaxis;
 
-               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);
+               status = 0;
+	       fits_open_image(&fptreal, testfitsreal, READWRITE, &status);
+               fits_write_subset(fptreal, TDOUBLE, fpixel, lpixel, image_real, &status);
+	       fits_close_file(fptreal, &status);
 
-               fclose(pFilereal);
-               fclose(pFileimg);
+               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