         printf("Using the same file!\n");
     }
     printf("\n");
-    ////////////////////////////////////////////////////////////////////////////
+
+    if(radec_input==1) int success = doCalcRaDec(infile0, infile1, outfile, silent_on_GPU_testing, scale_factor, nbins, hist_lower_range, hist_upper_range, hist_bin_width, log_binning_flag, two_different_files, conv_factor_angle);
+
+    else  int success = doCalcMpc(infile0, infile1, outfile, silent_on_GPU_testing, scale_factor, nbins, hist_lower_range, hist_upper_range, hist_bin_width, log_binning_flag, two_different_files, conv_factor_angle);
+
+}
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+// Module that runs the calculations for input files in the format ra:dec
+////////////////////////////////////////////////////////////////////////
+
+int doCalcRaDec(FILE *infile0, FILE *infile1, FILE *outfile, bool silent_on_GPU_testing, float scale_factor, int nbins, float hist_lower_range, float hist_upper_range, float hist_bin_width, int log_binning_flag, bool two_different_files, float conv_factor_angle){
+   
+
+
+
+    
+    float *d_alpha0, *d_delta0;
+    float *h_alpha0, *h_delta0;
+
+    float *d_alpha1, *d_delta1;
+    float *h_alpha1, *h_delta1;
+
+    int NUM_GALAXIES0;
+    int NUM_GALAXIES1;
+
+    //////////////////////////////////////////////////////////////////////
+    // Read in the galaxy files.
+    ////////////////////////////////////////////////////////////////////////////
+    // Read in the first file
+    ////////////////////////////////////////////////////////////////////////////
+
+    fscanf(infile0, "%d", &NUM_GALAXIES0);
+
+    int size_of_galaxy_array0 = NUM_GALAXIES0 * sizeof(float);    
+    printf("SIZE 0 # GALAXIES: %d\n",NUM_GALAXIES0);
+
+    h_alpha0 = (float*)malloc(size_of_galaxy_array0);
+    h_delta0 = (float*)malloc(size_of_galaxy_array0);
+    float temp0, temp1;
+
+    for(int i=0; i<NUM_GALAXIES0; i++)
+    {
+        fscanf(infile0, "%f %f", &temp0, &temp1);
+        h_alpha0[i] = temp0/scale_factor;
+        h_delta0[i] = temp1/scale_factor;
+        //if (i<10)
+        //printf("%e %e\n", h_alpha0[i], h_delta0[i]);
+    }
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Read in the second file
+    ////////////////////////////////////////////////////////////////////////////
+
+    fscanf(infile1, "%d", &NUM_GALAXIES1);
+
+    int size_of_galaxy_array1 = NUM_GALAXIES1 * sizeof(float);    
+    printf("SIZE 1 # GALAXIES: %d\n",NUM_GALAXIES1);
+
+    h_alpha1 = (float*)malloc(size_of_galaxy_array1);
+    h_delta1 = (float*)malloc(size_of_galaxy_array1);
+
+    for(int i=0; i<NUM_GALAXIES1; i++)
+    {
+        fscanf(infile1, "%f %f", &temp0, &temp1);
+        h_alpha1[i] = temp0/scale_factor;
+        h_delta1[i] = temp1/scale_factor;
+        //if (i<10)
+        //printf("%e %e\n", h_alpha1[i], h_delta1[i]);
+    }
+
+//get device diagnostics
+ if (!silent_on_GPU_testing) getDeviceDiagnostics(NUM_GALAXIES0+NUM_GALAXIES1, 2);
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Allocation of histogram
+    ///////////////////////////////////////////////////////////////////////////
+
+    int *hist, *dev_hist;
+
+    int size_hist = SUBMATRIX_SIZE * (nbins+2);
+    int size_hist_bytes = size_hist*sizeof(int);
+
+    hist = (int*)malloc(size_hist_bytes);
+    memset(hist, 0, size_hist_bytes);
+
+    printf("Size of histogram: %d bytes\n",size_hist_bytes);
+    cudaMalloc((void **) &dev_hist, (size_hist_bytes));
+    cudaMemset(dev_hist, 0, size_hist_bytes);
+
+    unsigned long  *hist_array;
+
+    int hist_array_size = (nbins+2) * sizeof(unsigned long);
+    hist_array =  (unsigned long*)malloc(hist_array_size);
+    printf("Size of histogram array: %d bytes\n",hist_array_size);
+    memset(hist_array,0,hist_array_size); 
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Define the grid and block size
+    ////////////////////////////////////////////////////////////////////////////
+    dim3 grid, block;
+    // 128*4 = 512, the amount of memory needed for one histogram.
+    // 8192*4 = 32768 is max memory to ask for for the histograms.
+    // 8192/128 = 64, is is the right number of blocks?
+    grid.x = 8192/(DEFAULT_NBINS+2); // Is this the number of blocks?
+    block.x = SUBMATRIX_SIZE/grid.x; // Is this the number of threads per block? NUM_GALAXIES/block.x;
+    // SUBMATRIX is the number of threads per warp? Per kernel call?
+    ////////////////////////////////////////////////////////////////////////////
+
+    cudaMalloc((void **) &d_alpha0, size_of_galaxy_array0 );
+    cudaMalloc((void **) &d_delta0, size_of_galaxy_array0 );
+
+    cudaMalloc((void **) &d_alpha1, size_of_galaxy_array1 );
+    cudaMalloc((void **) &d_delta1, size_of_galaxy_array1 );
+
+    // Check to see if we allocated enough memory.
+    if (0==d_alpha0 || 0==d_delta0 || 0==d_alpha1 || 0==d_delta1 || 0==dev_hist)
+    {
+        printf("couldn't allocate memory\n");
+        return 1;
+    }
+
+    // Initialize array to all 0's
+    cudaMemset(d_alpha0,0,size_of_galaxy_array0);
+    cudaMemset(d_delta0,0,size_of_galaxy_array0);
+    cudaMemset(d_alpha1,0,size_of_galaxy_array1);
+    cudaMemset(d_delta1,0,size_of_galaxy_array1);
+
+    cudaMemcpy(d_alpha0, h_alpha0, size_of_galaxy_array0, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_delta0, h_delta0, size_of_galaxy_array0, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_alpha1, h_alpha1, size_of_galaxy_array1, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_delta1, h_delta1, size_of_galaxy_array1, cudaMemcpyHostToDevice );
+
+    int x, y;
+    int num_submatrices_x = NUM_GALAXIES0 / SUBMATRIX_SIZE;
+    int num_submatrices_y = NUM_GALAXIES1 / SUBMATRIX_SIZE;
+    // Take care of edges of matrix.
+    if (NUM_GALAXIES0%SUBMATRIX_SIZE != 0)
+    {
+        num_submatrices_x += 1;
+    }
+    if (NUM_GALAXIES1%SUBMATRIX_SIZE != 0)
+    {
+        num_submatrices_y += 1;
+    }
+
+
+    printf("Breaking down the calculations.\n");
+    printf("Number of submatrices: %dx%d\n",num_submatrices_x,num_submatrices_y);
+    printf("Number of calculations per submatrices: %dx%d\n",SUBMATRIX_SIZE,SUBMATRIX_SIZE);
+
+    int bin_index = 0;
+    for(int k = 0; k < num_submatrices_y; k++)
+    {
+        y = k*SUBMATRIX_SIZE;
+        //printf("%d %d\n",k,y);
+        for(int j = 0; j < num_submatrices_x; j++)
+        {
+            x = j*SUBMATRIX_SIZE; 
+
+            //printf("----\n");
+            //printf("%d %d\t\t%d %d\n",k,y,j,x);
+            //printf("----\n");
+
+            // Set the histogram to all zeros each time.
+            cudaMemset(dev_hist,0,size_hist_bytes);
+
+            int max_x = NUM_GALAXIES0;
+            int max_y = NUM_GALAXIES1;
+
+            distance<<<grid,block>>>(d_alpha0, d_delta0,d_alpha1, d_delta1, x, y, max_x, max_y, dev_hist, hist_lower_range, hist_upper_range, nbins, hist_bin_width, log_binning_flag, two_different_files,conv_factor_angle);
+            cudaMemcpy(hist, dev_hist, size_hist_bytes, cudaMemcpyDeviceToHost);
+
+            ////////////////////////////////////////////////////////////////////
+            // Sum up the histograms from each thread (hist).
+            ////////////////////////////////////////////////////////////////////
+            for(int m=0; m<size_hist; m++)
+            {
+                bin_index = m%(nbins+2);
+                hist_array[bin_index] += hist[m];
+            }    
+        }  
+    }
+
+    unsigned long total = 0;
+
+    float lo = hist_lower_range;
+    float hi = 0;
+    for(int k=0; k<nbins+1; k++)
+    {
+        if (k==0)
+        {
+            //fprintf(outfile, "Underflow below %.3e %s %lu \n", lo, ",",  hist_array[k]);
+        }
+        else
+        {
+            if (log_binning_flag==0)
+            {
+                hi = lo + hist_bin_width;
+            }
+            else if (log_binning_flag==1)
+            {
+                //printf("lo: %f\t\tlog(lo): %f\n",lo,log(lo));
+                hi = exp(log(lo) + hist_bin_width);
+            }
+            else if (log_binning_flag==2)
+            {
+                //printf("lo: %f\t\tlog10(lo): %f\n",lo,log10(lo));
+                hi = pow(10,(log10(lo) + hist_bin_width));
+            }
+
+            fprintf(outfile, "%.3e %.3e %lu \n",lo,hi,hist_array[k]);
+            total += hist_array[k];
+
+            lo = hi;
+        }
+    }
+    printf("total: %lu \n", total);
+
+    fclose(infile0);
+    fclose(infile1);
+    fclose(outfile);
+
+    free(h_alpha0);
+    free(h_delta0);
+    free(h_alpha1);
+    free(h_delta1);
+    free(hist);
+
+    cudaFree(d_alpha0);
+    cudaFree(d_delta0);  
+    cudaFree(d_alpha1);
+    cudaFree(d_delta1);  
+xs    cudaFree(dev_hist);
+
+    return 0;
+}  
+//////////////////////////////////////////////////////////////////////
+
+
+xs
+
+
+
+
+////////////////////////////////////////////////////////////////////////
+
+////////////////////////////////////////////////////////////////////////
+// Module that runs the calculations for input files in the format x:y:z in Mpc
+////////////////////////////////////////////////////////////////////////
+
+int doCalcMpc(FILE *infile0, FILE *infile1, FILE *outfile, bool silent_on_GPU_testing, float scale_factor, int nbins, float hist_lower_range, float hist_upper_range, float hist_bin_width, int log_binning_flag, bool two_different_files, float conv_factor_angle){
+
+
+
+
+    float *d_x0, *d_y0, *d_z0;
+    float *h_x0, *h_y0, *h_z0;
+
+    float *d_x1, *d_y1, *d_z1;
+    float *h_x1, *h_y1, *h_z1;
+
+    int NUM_GALAXIES0;
+    int NUM_GALAXIES1;
+
+    //////////////////////////////////////////////////////////////////////
+    // Read in the galaxy files.
+    ////////////////////////////////////////////////////////////////////////////
+    // Read in the first file
+    ////////////////////////////////////////////////////////////////////////////
+
+    fscanf(infile0, "%d", &NUM_GALAXIES0);
+
+    int size_of_galaxy_array0 = NUM_GALAXIES0 * sizeof(float);    
+    printf("SIZE 0 # GALAXIES: %d\n",NUM_GALAXIES0);
+
+    h_x0 = (float*)malloc(size_of_galaxy_array0);
+    h_y0 = (float*)malloc(size_of_galaxy_array0);
+    h_z0 = (float*)malloc(size_of_galaxy_array0);
+    float temp0, temp1, temp2;
+
+    for(int i=0; i<NUM_GALAXIES0; i++)
+    {
+        fscanf(infile0, "%f %f", &temp0, &temp1, &temp2);
+        h_x0[i] = temp0/scale_factor;
+        h_y0[i] = temp1/scale_factor;
+        h_z0[i] = temp2/scale_factor;
+        //if (i<10)
+        //printf("%e %e\n", h_x0[i], h_y0[i], h_y0[i],);
+    }
+
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Read in the second file
+    ////////////////////////////////////////////////////////////////////////////
+
+    fscanf(infile1, "%d", &NUM_GALAXIES1);
+
+    int size_of_galaxy_array1 = NUM_GALAXIES1 * sizeof(float);    
+    printf("SIZE 1 # GALAXIES: %d\n",NUM_GALAXIES1);
+
+    h_x1 = (float*)malloc(size_of_galaxy_array1);
+    h_y1 = (float*)malloc(size_of_galaxy_array1);
+    h_z1 = (float*)malloc(size_of_galaxy_array1);
+
+    for(int i=0; i<NUM_GALAXIES1; i++)
+    {
+        fscanf(infile1, "%f %f", &temp0, &temp1, &temp2);
+        h_x1[i] = temp0/scale_factor;
+        h_y1[i] = temp1/scale_factor;
+        h_z1[i] = temp2/scale_factor;
+        //if (i<10)
+        //printf("%e %e\n", h_x1[i], h_y1[i], h_z1[i]);
+    }
+
+// get device diagnostics 
+ if (!silent_on_GPU_testing) getDeviceDiagnostics(NUM_GALAXIES0+NUM_GALAXIES1, 2);
+
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Allocation of histogram
+    ///////////////////////////////////////////////////////////////////////////
+
+    int *hist, *dev_hist;
+
+    int size_hist = SUBMATRIX_SIZE * (nbins+2);
+    int size_hist_bytes = size_hist*sizeof(int);
+
+    hist = (int*)malloc(size_hist_bytes);
+    memset(hist, 0, size_hist_bytes);
+
+    printf("Size of histogram: %d bytes\n",size_hist_bytes);
+    cudaMalloc((void **) &dev_hist, (size_hist_bytes));
+    cudaMemset(dev_hist, 0, size_hist_bytes);
+
+    unsigned long  *hist_array;
+
+    int hist_array_size = (nbins+2) * sizeof(unsigned long);
+    hist_array =  (unsigned long*)malloc(hist_array_size);
+    printf("Size of histogram array: %d bytes\n",hist_array_size);
+    memset(hist_array,0,hist_array_size); 
+
+    ////////////////////////////////////////////////////////////////////////////
+    // Define the grid and block size
+    ////////////////////////////////////////////////////////////////////////////
+    dim3 grid, block;
+    // 128*4 = 512, the amount of memory needed for one histogram.
+    // 8192*4 = 32768 is max memory to ask for for the histograms.
+    // 8192/128 = 64, is is the right number of blocks?
+    grid.x = 8192/(DEFAULT_NBINS+2); // Is this the number of blocks?
+    block.x = SUBMATRIX_SIZE/grid.x; // Is this the number of threads per block? NUM_GALAXIES/block.x;
+    // SUBMATRIX is the number of threads per warp? Per kernel call?
+    ////////////////////////////////////////////////////////////////////////////
+
+    cudaMalloc((void **) &d_x0, size_of_galaxy_array0 );
+    cudaMalloc((void **) &d_y0, size_of_galaxy_array0 );
+    cudaMalloc((void **) &d_z0, size_of_galaxy_array0 );
+
+    cudaMalloc((void **) &d_x1, size_of_galaxy_array1 );
+    cudaMalloc((void **) &d_y1, size_of_galaxy_array1 );
+    cudaMalloc((void **) &d_z1, size_of_galaxy_array1 );
+
+    // Check to see if we allocated enough memory.
+    if (0==d_x0 || 0==d_y0 || 0==d_z0 || 0==d_x1 || 0==d_y1 || 0==d_z1 || 0==dev_hist)
+    {
+        printf("couldn't allocate memory\n");
+        return 1;
+    }
+
+    // Initialize array to all 0's
+    cudaMemset(d_x0,0,size_of_galaxy_array0);
+    cudaMemset(d_y0,0,size_of_galaxy_array0);
+    cudaMemset(d_z0,0,size_of_galaxy_array0);
+    cudaMemset(d_x1,0,size_of_galaxy_array1);
+    cudaMemset(d_y1,0,size_of_galaxy_array1);
+    cudaMemset(d_z1,0,size_of_galaxy_array1);
+
+    cudaMemcpy(d_x0, h_x0, size_of_galaxy_array0, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_y0, h_y0, size_of_galaxy_array0, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_z0, h_z0, size_of_galaxy_array0, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_x1, h_x1, size_of_galaxy_array1, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_y1, h_y1, size_of_galaxy_array1, cudaMemcpyHostToDevice );
+    cudaMemcpy(d_z1, h_z1, size_of_galaxy_array1, cudaMemcpyHostToDevice );
+
+    int x, y;
+    int num_submatrices_x = NUM_GALAXIES0 / SUBMATRIX_SIZE;
+    int num_submatrices_y = NUM_GALAXIES1 / SUBMATRIX_SIZE;
+    // Take care of edges of matrix.
+    if (NUM_GALAXIES0%SUBMATRIX_SIZE != 0)
+    {
+        num_submatrices_x += 1;
+    }
+    if (NUM_GALAXIES1%SUBMATRIX_SIZE != 0)
+    {
+        num_submatrices_y += 1;
+    }
+
+
+    printf("Breaking down the calculations.\n");
+    printf("Number of submatrices: %dx%d\n",num_submatrices_x,num_submatrices_y);
+    printf("Number of calculations per submatrices: %dx%d\n",SUBMATRIX_SIZE,SUBMATRIX_SIZE);
+
+    int bin_index = 0;
+    for(int k = 0; k < num_submatrices_y; k++)
+    {
+        y = k*SUBMATRIX_SIZE;
+        //printf("%d %d\n",k,y);
+        for(int j = 0; j < num_submatrices_x; j++)
+        {
+            x = j*SUBMATRIX_SIZE; 
+
+            //printf("----\n");
+            //printf("%d %d\t\t%d %d\n",k,y,j,x);
+            //printf("----\n");
+
+            // Set the histogram to all zeros each time.
+            cudaMemset(dev_hist,0,size_hist_bytes);
+
+            int max_x = NUM_GALAXIES0;
+            int max_y = NUM_GALAXIES1;
+
+            distanceMpc<<<grid,block>>>(d_x0, d_y0, d_z0,d_x1, d_y1, d_z1, x, y, max_x, max_y, dev_hist, hist_lower_range, hist_upper_range, nbins, hist_bin_width, log_binning_flag, two_different_files,conv_factor_angle);
+            cudaMemcpy(hist, dev_hist, size_hist_bytes, cudaMemcpyDeviceToHost);
+
+            ////////////////////////////////////////////////////////////////////
+            // Sum up the histograms from each thread (hist).
+            ////////////////////////////////////////////////////////////////////
+            for(int m=0; m<size_hist; m++)
+            {
+                bin_index = m%(nbins+2);
+                hist_array[bin_index] += hist[m];
+            }    
+        }  
+    }
+
+    unsigned long total = 0;
+
+    float lo = hist_lower_range;
+    float hi = 0;
+    for(int k=0; k<nbins+1; k++)
+    {
+        if (k==0)
+        {
+            //fprintf(outfile, "Underflow below %.3e %s %lu \n", lo, ",",  hist_array[k]);
+        }
+        else
+        {
+            if (log_binning_flag==0)
+            {
+                hi = lo + hist_bin_width;
+            }
+            else if (log_binning_flag==1)
+            {
+                //printf("lo: %f\t\tlog(lo): %f\n",lo,log(lo));
+                hi = exp(log(lo) + hist_bin_width);
+            }
+            else if (log_binning_flag==2)
+            {
+                //printf("lo: %f\t\tlog10(lo): %f\n",lo,log10(lo));
+                hi = pow(10,(log10(lo) + hist_bin_width));
+            }
+
+            fprintf(outfile, "%.3e %.3e %lu \n",lo,hi,hist_array[k]);
+            total += hist_array[k];
+
+            lo = hi;
+        }
+    }
+    printf("total: %lu \n", total);
+
+    fclose(infile0);
+    fclose(infile1);
+    fclose(outfile);
+
+    free(h_x0);
+    free(h_y0);
+    free(h_z0);
+    free(h_x1);
+    free(h_y1);
+    free(h_z1);
+    free(hist);
+
+    cudaFree(d_x0);
+    cudaFree(d_y0);
+    cudaFree(d_z0);  
+    cudaFree(d_x1);
+    cudaFree(d_y1);
+    cudaFree(d_z1);  
+    cudaFree(dev_hist);
+
+    return 0;
+}  
+//////////////////////////////////////////////////////////////////////
+
+
+
+
+
+void getDeviceDiagnostics(int tot_gals, int n_coords){
+
+ ////////////////////////////////////////////////////////////////////////////
     // Now get the info from the device.
     ////////////////////////////////////////////////////////////////////////////
-    if (!silent_on_GPU_testing)
-    {
+   
         printf("\n------ CUDA device diagnostics ------\n\n");
 
-        int tot_gals = 100000;
+        
         int nx = SUBMATRIX_SIZE;
         int ncalc = nx * nx;
-        int gpu_mem_needed = int(tot_gals * sizeof(float)) * 2; // need to allocate ra, dec.
+        int gpu_mem_needed = int(tot_gals * sizeof(float)) * n_coords; // need to allocate ra, dec.
         printf("Requirements: %d calculations and %d bytes memory on the GPU \n\n", ncalc, gpu_mem_needed);
 
         int deviceCount = 0;


    float scale_factor = 1.0; // For if we need to convert input to arcsec or arcmin
    float conv_factor_angle = 57.2957795; // 180/pi // For if we need to convert arcdistance to arcsec or arcmin

    bool silent_on_GPU_testing = false;

    ////////////////////////////////////////////////////////////////////////////////

    }
    ////////////////////////////////////////////////////////////////////////////


