CFITSIO and parallel images writing implementation

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

w-stacking-fftw.c

@@ -985,28 +985,65 @@ if(rank == 0){
         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 PARALLEL_FITS
+        #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
+
+        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
 	    MPI_Barrier(MPI_COMM_WORLD);
             #endif
 	    if(isector == rank)
-            //if(rank == 0)
 	    {
                #ifdef FITSIO
 
 	       printf("%d writing\n",isector);
                //long * fpixel = (long *) malloc(sizeof(long)*naxis);
                //long * lpixel = (long *) malloc(sizeof(long)*naxis);
-               
+
                fpixel[0] = 1;
                fpixel[1] = isector*yaxis+1;
                lpixel[0] = xaxis;
                lpixel[1] = (isector+1)*yaxis;
-               //printf("fpixel %d, %d\n", fpixel[0], fpixel[1]);
-	       //printf("lpixel %d, %d\n", lpixel[0], lpixel[1]);
-	       
+               
                status = 0;
 	       fits_open_image(&fptreal, testfitsreal, READWRITE, &status);
                fits_write_subset(fptreal, TDOUBLE, fpixel, lpixel, image_real, &status);
@@ -1033,10 +1070,12 @@ if(rank == 0){
                fclose(pFileimg);
 	    }
 	}
+	#endif
 	#ifdef USE_MPI
 	MPI_Barrier(MPI_COMM_WORLD);
         #endif
 
+
 #endif //WRITE_IMAGE