text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "hdf5.h"
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_bessel.h>
#include "mclib.h"

#define PROP_DIM1 1
#define PROP_DIM2 8
#define PROP_DIM3 8
#define COORD_DIM1 2

//define constants
const double A_RAD=7.56e-15, C_LIGHT=2.99792458e10, PL_CONST=6.6260755e-27;
const double K_B=1.380658e-16, M_P=1.6726231e-24, THOMP_X_SECT=6.65246e-25, M_EL=9.1093879e-28  ;

void printPhotons(struct photon *ph, int num_ph, int frame,char dir[200] )
{
    //function to save the photons' positions and 4 momentum
    int i=0;
    char mc_file_p0[200], mc_file_p1[200],mc_file_p2[200], mc_file_p3[200];
    char mc_file_r0[200], mc_file_r1[200], mc_file_r2[200], mc_file_ns[200];
    FILE *fPtr=NULL, *fPtr1=NULL,*fPtr2=NULL,*fPtr3=NULL,*fPtr4=NULL,*fPtr5=NULL,*fPtr6=NULL,*fPtr7=NULL;

    //make strings for proper files
    snprintf(mc_file_p0,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_P0.dat" );
    snprintf(mc_file_p1,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_P1.dat" );
    snprintf(mc_file_p2,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_P2.dat" );
    snprintf(mc_file_p3,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_P3.dat" );
    snprintf(mc_file_r0,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_R0.dat" );
    snprintf(mc_file_r1,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_R1.dat" );
    snprintf(mc_file_r2,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_R2.dat" );
    snprintf(mc_file_ns,sizeof(mc_file_p0),"%s%s%d%s",dir,"mcdata_", frame,"_NS.dat" ); //for number of scatterings each photon went through

    //save the energy
    fPtr=fopen(mc_file_p0, "a");
    fPtr1=fopen(mc_file_p1, "a");
    fPtr2=fopen(mc_file_p2, "a");
    fPtr3=fopen(mc_file_p3, "a");
    fPtr4=fopen(mc_file_r0, "a");
    fPtr5=fopen(mc_file_r1, "a");
    fPtr6=fopen(mc_file_r2, "a");
    fPtr7=fopen(mc_file_ns, "a");

    //printf("Writing P0\n");
    for (i=0;i<num_ph;i++)
    {
        fprintf(fPtr,"%0.13e\t",  (ph+i)->p0);
        //printf("%d: %0.13e \n", i, (ph+i)->p0);

        fprintf(fPtr1,"%0.13e\t",  (ph+i)->p1);
        //printf("%d: %0.13e \n", i, (ph+i)->p1);

        fprintf(fPtr2,"%0.13e\t",  (ph+i)->p2);
        //printf("%d: %0.13e \n", i, (ph+i)->p2);

        fprintf(fPtr3,"%0.13e\t",  (ph+i)->p3);
        //printf("%d: %0.13e \n", i, (ph+i)->p3);

        fprintf(fPtr4,"%0.13e\t",  (ph+i)->r0);
        //printf("%d: %0.13e \n", i, (ph+i)->r0);

        fprintf(fPtr5,"%0.13e\t",  (ph+i)->r1);
        //printf("%d: %0.13e \n", i, (ph+i)->r1);

        fprintf(fPtr6,"%0.13e\t",  (ph+i)->r2);
        //printf("%d: %0.13e \n", i, (ph+i)->r2);

        //fprintf(fPtr7,"%0.13e\t",  *(ph_num_scatt+i));
        fprintf(fPtr7,"%e\t",  (ph+i)->num_scatt);
        //printf("%d: %0.13e \n", i, (ph+i)->num_scatt);

    }
    fclose(fPtr);
    fclose(fPtr1);
    fclose(fPtr2);
    fclose(fPtr3);
    fclose(fPtr4);
    fclose(fPtr5);
    fclose(fPtr6);
    fclose(fPtr7);

    //printf("%s\n%s\n%s\n", mc_file_p0, mc_file_r0, mc_file_ns);
}

void saveCheckpoint(char dir[200], int frame, int scatt_frame, int ph_num,double time_now, struct photon *ph, int last_frame )
{
    //function to save data necessary to restart simulation if it ends
    //need to save all photon data
    FILE *fPtr=NULL;
    char checkptfile[200];
    char restart;
    int i=0;

    snprintf(checkptfile,sizeof(checkptfile),"%s%s",dir,"mc_chkpt.dat" );

    fPtr=fopen(checkptfile, "wb");

    if (scatt_frame!=last_frame)
    {
        restart='c';
        fwrite(&restart, sizeof(char), 1, fPtr);
        fwrite(&frame, sizeof(int), 1, fPtr);
        fwrite(&scatt_frame, sizeof(int), 1, fPtr);
        fwrite(&time_now, sizeof(double), 1, fPtr);
        fwrite(&ph_num, sizeof(int), 1, fPtr);

        //for(i=0;i<ph_num;i++)
        //{
            fwrite((ph), sizeof (struct photon )*ph_num, ph_num, fPtr);
        //}

    }
    else
    {
        //just finished last iteration of scatt_frame
        restart='r';
        fwrite(&restart, sizeof(char), 1, fPtr);
        fwrite(&frame, sizeof(int), 1, fPtr);
    }
    fclose(fPtr);

}

void readCheckpoint(char dir[200], struct photon **ph, int *framestart, int *scatt_framestart, int *ph_num, char *restart, double *time )
{
    //function to read in data from checkpoint file
    FILE *fPtr=NULL;
    char checkptfile[200];
    int i=0;
    //int frame, scatt_frame, ph_num, i=0;
    struct photon *phHolder=NULL; //pointer to struct to hold data read in from checkpoint file

    snprintf(checkptfile,sizeof(checkptfile),"%s%s",dir,"mc_chkpt.dat" );

    fPtr=fopen(checkptfile, "rb");

    fread(restart, sizeof(char), 1, fPtr);
    printf("%c\n", *restart);
    fread(framestart, sizeof(int), 1, fPtr);
    printf("%d\n", *framestart);

    if((*restart)=='c')
    {
        fread(scatt_framestart, sizeof(int), 1, fPtr);
        *scatt_framestart+=1; //add one to start at the next frame after the siomulation was interrrupted
        printf("%d\n", *scatt_framestart);
        fread(time, sizeof(double), 1, fPtr);
        printf("%e\n", *time);
        fread(ph_num, sizeof(int), 1, fPtr);
        printf("%d\n", *ph_num);

        phHolder=malloc(sizeof(struct photon));
        (*ph)=malloc(sizeof(struct photon)*(*ph_num)); //allocate memory to hold photon data

        for (i=0;i<(*ph_num);i++)
        {
            fread(phHolder, sizeof(struct photon), 1, fPtr);
            //printf("%e,%e,%e, %e,%e,%e, %e, %e\n",(ph)->p0, (ph)->p1, (ph)->p2, ph->p3, (ph)->r0, (ph)->r1, (ph)->r2, ph->num_scatt );

            (*ph)[i].p0=phHolder->p0;
            (*ph)[i].p1=phHolder->p1;
            (*ph)[i].p2=phHolder->p2;
            (*ph)[i].p3=phHolder->p3;
            (*ph)[i].r0= phHolder->r0;
            (*ph)[i].r1=phHolder->r1 ;
            (*ph)[i].r2=phHolder->r2;
            (*ph)[i].num_scatt=phHolder->num_scatt;

        }

        free(phHolder);
    }
    else
    {
        *framestart+=1; //if the  checkpoint file saved and the program was inturrupted before the frame variable had just increased and before the scatt_frame iteration was saved, add one to the frame start
    }

    fclose(fPtr);
}

void readMcPar(char file[200], double *fps, double *theta_jmin, double *theta_j, double *inj_radius, int *frm0,int *last_frm, int *frm2, int *photon_num, double *ph_weight, char *spect, char *restart)
{
    //function to read mc.par file
	FILE *fptr=NULL;
	char buf[100];
	double theta_deg;

	//open file
	fptr=fopen(file,"r");
	//read in frames per sec and other variables outlined in main()
	fscanf(fptr, "%lf",fps);
	//printf("%f\n", *fps );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%d",frm0);
	//printf("%d\n", *frm0 );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%d",last_frm);
	//printf("%d\n", *last_frm );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%d",frm2);
    *frm2+=*frm0; //frame to go to is what is given in the file plus the starting frame
	//printf("%d\n", *frm2 );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%d",photon_num);
	//printf("%d\n", *photon_num );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%lf",inj_radius);
	//printf("%lf\n", *inj_radius );

	fgets(buf, 100,fptr);
	//theta jmin
	fscanf(fptr, "%lf",&theta_deg);
	*theta_jmin=theta_deg*M_PI/180;
	//printf("%f\n", *theta_jmin );

	fgets(buf, 100,fptr);

	fscanf(fptr, "%lf",&theta_deg);
    *theta_j=theta_deg*M_PI/180;
	//printf("%f\n", *theta_j );

	fgets(buf, 100,fptr);

    fscanf(fptr, "%lf",ph_weight);
    fgets(buf, 100,fptr);

    *spect=getc(fptr);
    fgets(buf, 100,fptr);
    //printf("%c\n",*spect);

    *restart=getc(fptr);

	//close file
	fclose(fptr);
}

void readAndDecimate(char flash_file[200], double r_inj, double **x, double **y, double **szx, double **szy, double **r,\
 double **theta, double **velx, double **vely, double **dens, double **pres, double **gamma, double **dens_lab, double **temp, int *number)
{
    //function to read in data from FLASH file
    hid_t  file,dset, space;
    herr_t status;
    hsize_t dims[2]={0,0}; //hold dimension size for coordinate data set (mostly interested in dims[0])
    double **vel_x_buffer=NULL, **vel_y_buffer=NULL, **dens_buffer=NULL, **pres_buffer=NULL, **coord_buffer=NULL, **block_sz_buffer=NULL;
    double *velx_unprc=NULL, *vely_unprc=NULL, *dens_unprc=NULL, *pres_unprc=NULL, *x_unprc=NULL, *y_unprc=NULL, *r_unprc=NULL, *szx_unprc=NULL, *szy_unprc=NULL;
    int  i,j,count,x1_count, y1_count, r_count, **node_buffer=NULL, num_nodes=0;
    double x1[8]={-7.0/16,-5.0/16,-3.0/16,-1.0/16,1.0/16,3.0/16,5.0/16,7.0/16};

    file = H5Fopen (flash_file, H5F_ACC_RDONLY, H5P_DEFAULT);
    printf(">> mc.py: Reading positional, density, pressure, and velocity information...\n");

    //printf("Reading coord\n");
    dset = H5Dopen (file, "coordinates", H5P_DEFAULT);

    //get dimensions of array and save it
    space = H5Dget_space (dset);

    H5Sget_simple_extent_dims(space, dims, NULL); //save dimesnions in dims

    /*
     * Allocate array of pointers to rows. OPTIMIZE HERE: INITALIZE ALL THE BUFFERS AT ONCE IN 1 FOR LOOP
     */
    coord_buffer = (double **) malloc (dims[0] * sizeof (double *));

    coord_buffer[0] = (double *) malloc (dims[0] * dims[1] * sizeof (double));

    block_sz_buffer= (double **) malloc (dims[0] * sizeof (double *));

    block_sz_buffer[0] = (double *) malloc (dims[0] * COORD_DIM1 * sizeof (double));

    node_buffer= (int **) malloc (dims[0] * sizeof (int *));
    node_buffer[0] = (int *) malloc (dims[0] * sizeof (int));

    vel_x_buffer= (double **) malloc (dims[0] * sizeof (double *));
    vel_x_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));

    vel_y_buffer= (double **) malloc (dims[0] * sizeof (double *));
    vel_y_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));

    dens_buffer= (double **) malloc (dims[0] * sizeof (double *));
    dens_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));

    pres_buffer= (double **) malloc (dims[0] * sizeof (double *));
    pres_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));

    /*
     * Set the rest of the pointers to rows to the correct addresses.
     */
     for (i=1; i<dims[0]; i++)
     {
         coord_buffer[i] = coord_buffer[0] + i * dims[1];
         block_sz_buffer[i] = block_sz_buffer[0] + i * COORD_DIM1;
         node_buffer[i] = node_buffer[0] + i ;
         vel_x_buffer[i] = vel_x_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
         vel_y_buffer[i] = vel_y_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
         dens_buffer[i] = dens_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
         pres_buffer[i] = pres_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
     }

//    for (i=1; i<dims[0]; i++)
//        coord_buffer[i] = coord_buffer[0] + i * dims[1];

    //read data such that first column is x and second column is y
    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,coord_buffer[0]);

    //close dataset
    status = H5Sclose (space);
    status = H5Dclose (dset);

    //printf("Reading block size\n");

    dset = H5Dopen (file, "block size", H5P_DEFAULT);

//    block_sz_buffer= (double **) malloc (dims[0] * sizeof (double *));
//
//    block_sz_buffer[0] = (double *) malloc (dims[0] * COORD_DIM1 * sizeof (double));
//
//    for (i=1; i<dims[0]; i++){
//        block_sz_buffer[i] = block_sz_buffer[0] + i * COORD_DIM1;
//    }
    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,block_sz_buffer[0]);

    // first column of buffer is x and second column is y
    status = H5Dclose (dset);

    //printf("Reading node type\n");
    dset = H5Dopen (file, "node type", H5P_DEFAULT);

//    node_buffer= (int **) malloc (dims[0] * sizeof (int *));
//    node_buffer[0] = (int *) malloc (dims[0] * sizeof (int));
//
//    for (i=1; i<dims[0]; i++){
//        node_buffer[i] = node_buffer[0] + i ;
//    }
    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,node_buffer[0]);
    status = H5Dclose (dset);

    //printf("Reading velx\n");
    dset = H5Dopen (file, "velx", H5P_DEFAULT);

//   vel_x_buffer= (double **) malloc (dims[0] * sizeof (double *));
//   vel_x_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));
//
//   for (i=1; i<dims[0]; i++)
//   {
//        vel_x_buffer[i] = vel_x_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
//   }
   //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,vel_x_buffer[0]);
    status = H5Dclose (dset);

    //printf("Reading vely\n");
    dset = H5Dopen (file, "vely", H5P_DEFAULT);

//     vel_y_buffer= (double **) malloc (dims[0] * sizeof (double *));
//     vel_y_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));
//
//     for (i=1; i<dims[0]; i++)
//     {
//        vel_y_buffer[i] = vel_y_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
//     }
    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,vel_y_buffer[0]);
    status = H5Dclose (dset);

    //printf("Reading dens\n");
    dset = H5Dopen (file, "dens", H5P_DEFAULT);

//    dens_buffer= (double **) malloc (dims[0] * sizeof (double *));
//    dens_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));
//     for (i=1; i<dims[0]; i++)
//     {
//        dens_buffer[i] = dens_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
//     }

    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,dens_buffer[0]);
    status = H5Dclose (dset);

    //printf("Reading pres\n");

    dset = H5Dopen (file, "pres", H5P_DEFAULT);

//    pres_buffer= (double **) malloc (dims[0] * sizeof (double *));
//    pres_buffer[0]= (double *) malloc (dims[0] * PROP_DIM1  *PROP_DIM2*PROP_DIM3* sizeof (double));
//     for (i=1; i<dims[0]; i++)
//     {
//        pres_buffer[i] = pres_buffer[0] + i * PROP_DIM1*PROP_DIM2*PROP_DIM3;
//     }
    //printf("Reading Dataset\n");
    status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,pres_buffer[0]);
    status = H5Dclose (dset);

    status = H5Fclose (file);

    printf(">> Selecting good node types (=1)\n");
    //find out how many good nodes there are
    for (i=0;i<dims[0];i++)
    {
        if (node_buffer[i][0]==1 ){
            num_nodes++;
        }
    }

    //allocate memory for arrays to hold unprocessed data
    pres_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    dens_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    velx_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    vely_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    x_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    y_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    r_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    szx_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    szy_unprc=malloc (num_nodes* PROP_DIM1  *PROP_DIM2*PROP_DIM3 * sizeof (double ));

    //find where the good values corresponding to the good gones (=1) and save them to the previously allocated pointers which are 1D arrays
    //also create proper x and y arrays and block size arrays
    //and then free up the buffer memory space
    printf(">> Creating and reshaping arrays\n");
    count=0;
    for (i=0;i<dims[0];i++)
    {
        if (node_buffer[i][0]==1 )
        {
            x1_count=0;
            y1_count=0;
            for (j=0;j<(PROP_DIM1*PROP_DIM2*PROP_DIM3);j++)
            {
                *(pres_unprc+count)=pres_buffer[i][j];
                *(dens_unprc+count)=dens_buffer[i][j];
                *(velx_unprc+count)=vel_x_buffer[i][j];
                *(vely_unprc+count)=vel_y_buffer[i][j];
                *(szx_unprc+count)=((block_sz_buffer[i][0])/8)*1e9; //divide by 8 for resolution, multiply by 1e9 to scale properly?
                *(szy_unprc+count)=((block_sz_buffer[i][1])/8)*1e9;
                if (j%8==0)
                {
                    x1_count=0;
                }
                if ((j%8==0) &&  (j!=0))
                {
                    y1_count++;
                }
                *(x_unprc+count)=(coord_buffer[i][0]+block_sz_buffer[i][0]*x1[x1_count])*1e9;
                *(y_unprc+count)=(coord_buffer[i][1]+block_sz_buffer[i][1]*x1[y1_count])*1e9;

                //printf("%d,%d,%d,%d\n",count,j,x1_count,y1_count);
                x1_count++;
                count++;
            }
        }
    }
    free (pres_buffer[0]); free (dens_buffer[0]);free (vel_x_buffer[0]);free (vel_y_buffer[0]); free(coord_buffer[0]);free(block_sz_buffer[0]);free(node_buffer[0]);
    free (pres_buffer);free(dens_buffer);free(vel_x_buffer);free(vel_y_buffer);free(coord_buffer);free(block_sz_buffer);free(node_buffer);

    //fill in radius array and find in how many places r > injection radius
    r_count=0;
    for (i=0;i<count;i++)
    {
        *(r_unprc+i)=pow((pow(*(x_unprc+i),2)+pow(*(y_unprc+i),2)),0.5);
        if (*(r_unprc+i)> (0.95*r_inj) )
        {
            r_count++;
        }
    }
        /*
    //find in how many places r > injection radius
    r_count=0;
    for (i=0;i<count;i++)
    {
        if (*(r_unprc+i)> (0.95*r_inj) )
        {
            r_count++;
        }
    }
    */
    //allocate memory to hold processed data
    (*pres)=malloc (r_count * sizeof (double ));
    (*velx)=malloc (r_count * sizeof (double ));
    (*vely)=malloc (r_count * sizeof (double ));
    (*dens)=malloc (r_count * sizeof (double ));
    (*x)=malloc (r_count * sizeof (double ));
    (*y)=malloc (r_count * sizeof (double ));
    (*r)=malloc (r_count * sizeof (double ));
    (*theta)=malloc (r_count * sizeof (double ));
    (*gamma)=malloc (r_count * sizeof (double ));
    (*dens_lab)=malloc (r_count * sizeof (double ));
    (*szx)=malloc (r_count * sizeof (double ));
    (*szy)=malloc (r_count * sizeof (double ));
    (*temp)=malloc (r_count * sizeof (double ));

    //assign values based on r> 0.95*r_inj
    j=0;
    for (i=0;i<count;i++)
    {
        if (*(r_unprc+i)> (0.95*r_inj) )
        {
            (*pres)[j]=*(pres_unprc+i);
            (*velx)[j]=*(velx_unprc+i);
            (*vely)[j]=*(vely_unprc+i);
            (*dens)[j]=*(dens_unprc+i);
            (*x)[j]=*(x_unprc+i);
            (*y)[j]=*(y_unprc+i);
            (*r)[j]=*(r_unprc+i);
            (*szx)[j]=*(szx_unprc+i);
            (*szy)[j]=*(szy_unprc+i);
            (*theta)[j]=atan2( *(x_unprc+i) , *(y_unprc+i) );//theta in radians in relation to jet axis
            (*gamma)[j]=pow(pow(1.0-(pow(*(velx_unprc+i),2)+pow(*(vely_unprc+i),2)),0.5),-1); //v is in units of c
            (*dens_lab)[j]= (*(dens_unprc+i)) * (pow(pow(1.0-(pow(*(velx_unprc+i),2)+pow(*(vely_unprc+i),2)),0.5),-1));
            (*temp)[j]=pow(3*(*(pres_unprc+i))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
            j++;
        }
    }
    *number=j;

    free(pres_unprc); free(velx_unprc);free(vely_unprc);free(dens_unprc);free(x_unprc); free(y_unprc);free(r_unprc);free(szx_unprc);free(szy_unprc);

}

void photonInjection( struct photon **ph, int *ph_num, double r_inj, double ph_weight, char spect, int array_length, double fps, double theta_min, double theta_max,\
double *x, double *y, double *szx, double *szy, double *r, double *theta, double *temps, double *vx, double *vy, gsl_rng * rand)
{
    int i=0, block_cnt=0, *ph_dens=NULL, ph_tot=0, j=0,k=0;
    double ph_dens_calc=0.0, fr_dum=0.0, y_dum=0.0, yfr_dum=0.0, fr_max=0, bb_norm=0, position_phi;
    double com_v_phi, com_v_theta, *p_comv=NULL, *boost=NULL; //comoving phi, theta, comoving 4 momentum for a photon, and boost for photon(to go to lab frame)
    double *l_boost=NULL; //pointer to hold array of lorentz boost, to lab frame, values
    float num_dens_coeff;

    if (spect=='w') //from MCRAT paper, w for wien spectrum
    {
        num_dens_coeff=8.44;
        printf("in wien spectrum\n");
    }
    else
    {
        num_dens_coeff=20.29; //this is for black body spectrum
        //printf("in BB spectrum");
    }

    //find how many blocks are near the injection radius within the angles defined in mc.par, get temperatures and calculate number of photons to allocate memory for
    //and then rcord which blocks have to have "x" amount of photons injected there

    for(i=0;i<array_length;i++)
    {
        //look at all boxes in width delta r=c/fps and within angles we are interested in NEED TO IMPLEMENT
            if ((*(r+i) > (r_inj - C_LIGHT/fps))  &&   (*(r+i)  < (r_inj + C_LIGHT/fps)  ) && (*(theta+i)< theta_max) && (*(theta+i) > theta_min) )
            {
                block_cnt++;
            }
    }
    //printf("Blocks: %d\n", block_cnt);

    //allocate memory to record density of photons for each block
    ph_dens=malloc(block_cnt * sizeof(int));

    //calculate the photon density for each block and save it to the array
    j=0;
    ph_tot=0;
    for (i=0;i<array_length;i++)
    {
        //printf("%d\n",i);
        //printf("%e, %e, %e, %e, %e, %e\n", *(r+i),(r_inj - C_LIGHT/fps), (r_inj + C_LIGHT/fps), *(theta+i) , theta_max, theta_min);
            if ((*(r+i) > (r_inj - C_LIGHT/fps))  &&   (*(r+i)  < (r_inj + C_LIGHT/fps)  ) && (*(theta+i)< theta_max) && (*(theta+i) > theta_min) )
            {
                ph_dens_calc=num_dens_coeff*2.0*M_PI*(*(x+i))*pow(*(temps+i),3.0)*pow(*(szx+i),2.0) /(ph_weight) ; //a*T^3/(weight) dV, dV=2*PI*x*dx^2,

                 *(ph_dens+j)=gsl_ran_poisson(rand,ph_dens_calc) ; //choose from poission distribution with mean of ph_dens_calc

                //printf("%d, %lf \n",*(ph_dens+j), ph_dens_calc);

                 //sum up all the densities to get total number of photons
                 ph_tot+=(*(ph_dens+j));

                 j++;
            }
    }

    //printf("%d\n", ph_tot);

    //allocate memory for that many photons and also allocate memory to hold comoving 4 momentum of each photon and the velocity of the fluid
    (*ph)=malloc (ph_tot * sizeof (struct photon ));

    p_comv=malloc(4*sizeof(double));
    boost=malloc(3*sizeof(double));
    l_boost=malloc(4*sizeof(double));

    //go through blocks and assign random energies/locations to proper number of photons
    ph_tot=0;
    k=0;
     for (i=0;i<array_length;i++)
    {
       if ((*(r+i) > (r_inj - C_LIGHT/fps))  &&   (*(r+i)  < (r_inj + C_LIGHT/fps)  ) && (*(theta+i)< theta_max) && (*(theta+i) > theta_min) )
        {
            //*(temps+i)=0.76*(*(temps+i));
            for(j=0;j<(*(ph_dens+k));j++ )
            {
                    //have to get random frequency for the photon comoving frequency
                    y_dum=1; //initalize loop
                    yfr_dum=0;
                    while (y_dum>yfr_dum)
                    {
                        fr_dum=gsl_rng_uniform_pos(rand)*6.3e11*(*(temps+i)); //in Hz
                        //printf("%lf, %lf ",gsl_rng_uniform_pos(rand), (*(temps+i)));
                        y_dum=gsl_rng_uniform_pos(rand);
                        //printf("%lf ",fr_dum);

                        if (spect=='w')
                        {
                            yfr_dum=(1.0/(1.29e31))*pow((fr_dum/(*(temps+i))),3.0)/(exp((PL_CONST*fr_dum)/(K_B*(*(temps+i)) ))-1); //curve is normalized to maximum
                        }
                        else
                        {
                            fr_max=(C_LIGHT*(*(temps+i)))/(0.29); //max frequency of bb
                            bb_norm=pow((fr_max/(*(temps+i))),2.0)/(exp(PL_CONST*fr_max/K_B/(*(temps+i)))-1); //find value of bb at fr_max
                            yfr_dum=(1.0/bb_norm)*pow((fr_dum/(*(temps+i))),2.0)/(exp((PL_CONST*fr_dum)/(K_B*(*(temps+i)) ))-1);//curve is normalized to vaue of bb @ max frequency
                        }
                        //printf("%lf, %lf,%lf,%e \n",(*(temps+i)),fr_dum, y_dum, yfr_dum);

                    }
                   //printf("%lf\n ",fr_dum);
                   position_phi=gsl_rng_uniform(rand)*2*M_PI;
                   com_v_phi=gsl_rng_uniform(rand)*2*M_PI;
                   com_v_theta=acos((gsl_rng_uniform(rand)*2)-1);
                   //printf("%lf, %lf, %lf\n", position_phi, com_v_phi, com_v_theta);

                   //populate 4 momentum comoving array
                   *(p_comv+0)=PL_CONST*fr_dum/C_LIGHT;
                   *(p_comv+1)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*cos(com_v_phi);
                   *(p_comv+2)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*sin(com_v_phi);
                   *(p_comv+3)=(PL_CONST*fr_dum/C_LIGHT)*cos(com_v_theta);

                    //populate boost matrix, not sure why multiplying by -1, seems to give correct answer in old python code...
                    *(boost+0)=-1*(*(vx+i))*cos(position_phi);
                    *(boost+1)=-1*(*(vx+i))*sin(position_phi);
                    *(boost+2)=-1*(*(vy+i));
                    //printf("%lf, %lf, %lf\n", *(boost+0), *(boost+1), *(boost+2));

                    //boost to lab frame
                    lorentzBoost(boost, p_comv, l_boost, 'p');
                    //printf("Assignemnt: %e, %e, %e, %e\n", *(l_boost+0), *(l_boost+1), *(l_boost+2),*(l_boost+3));

                (*ph)[ph_tot].p0=(*(l_boost+0));
                (*ph)[ph_tot].p1=(*(l_boost+1));
                (*ph)[ph_tot].p2=(*(l_boost+2));
                (*ph)[ph_tot].p3=(*(l_boost+3));
                (*ph)[ph_tot].r0= (*(x+i))*cos(position_phi); //put photons @ center of box that they are supposed to be in with random phi
                (*ph)[ph_tot].r1=(*(x+i))*sin(position_phi) ;
                (*ph)[ph_tot].r2=(*(y+i)); //y coordinate in flash becomes z coordinate in MCRaT
                (*ph)[ph_tot].num_scatt=0;
                //printf("%d\n",ph_tot);
                ph_tot++;
            }
            k++;
        }
    }

    *ph_num=ph_tot; //save number of photons
    free(ph_dens); free(p_comv);free(boost); free(l_boost);

}

void lorentzBoost(double *boost, double *p_ph, double *result, char object)
{
    //function to perform lorentz boost
    //if doing boost for an electron last argument is 'e' and there wont be a check for zero norm
    //if doing boost for a photon  last argument is 'p' and there will be a check for zero norm
    double beta=0, gamma=0, *boosted_p=NULL;

    gsl_vector_view b=gsl_vector_view_array(boost, 3); //make boost pointer into vector
    gsl_vector_view p=gsl_vector_view_array(p_ph, 4); //make boost pointer into vector
    gsl_matrix *lambda1= gsl_matrix_calloc (4, 4); //create matrix thats 4x4 to do lorentz boost
    gsl_vector *p_ph_prime =gsl_vector_calloc(4); //create vestor to hold lorentz boosted vector

    //printf("Boost: %e, %e, %e, %e\n",gsl_blas_dnrm2(&b.vector), *(boost+0), *(boost+1), *(boost+2));
    //printf("4 Momentum to Boost: %e, %e, %e, %e\n",*(p_ph+0), *(p_ph+1), *(p_ph+2), *(p_ph+3));

    //if magnitude of fluid velocity is != 0 do lorentz boost otherwise dont need to do a boost
    if (gsl_blas_dnrm2(&b.vector) > 0)
    {
        //printf("in If\n");
        beta=gsl_blas_dnrm2(&b.vector);
        gamma=1.0/sqrt(1-pow(beta, 2.0));
        //printf("Beta: %e\tGamma: %e\n",beta,gamma );

        //initalize matrix values
        gsl_matrix_set(lambda1, 0,0, gamma);
        gsl_matrix_set(lambda1, 0,1,  -1*gsl_vector_get(&b.vector,0)*gamma);
        gsl_matrix_set(lambda1, 0,2,  -1*gsl_vector_get(&b.vector,1)*gamma);
        gsl_matrix_set(lambda1, 0,3,  -1*gsl_vector_get(&b.vector,2)*gamma);
        gsl_matrix_set(lambda1, 1,1,  1+((gamma-1)*pow(gsl_vector_get(&b.vector,0),2.0)/pow(beta,2.0) ) );
        gsl_matrix_set(lambda1, 1,2,  ((gamma-1)*(gsl_vector_get(&b.vector,0)*  gsl_vector_get(&b.vector,1)/pow(beta,2.0) ) ));
        gsl_matrix_set(lambda1, 1,3,  ((gamma-1)*(gsl_vector_get(&b.vector,0)*  gsl_vector_get(&b.vector,2)/pow(beta,2.0) ) ));
        gsl_matrix_set(lambda1, 2,2,  1+((gamma-1)*pow(gsl_vector_get(&b.vector,1),2.0)/pow(beta,2.0) ) );
        gsl_matrix_set(lambda1, 2,3,  ((gamma-1)*(gsl_vector_get(&b.vector,1)*  gsl_vector_get(&b.vector,2)/pow(beta,2.0)) ) );
        gsl_matrix_set(lambda1, 3,3,  1+((gamma-1)*pow(gsl_vector_get(&b.vector,2),2.0)/pow(beta,2.0) ) );

        gsl_matrix_set(lambda1, 1,0, gsl_matrix_get(lambda1,0,1));
        gsl_matrix_set(lambda1, 2,0, gsl_matrix_get(lambda1,0,2));
        gsl_matrix_set(lambda1, 3,0, gsl_matrix_get(lambda1,0,3));
        gsl_matrix_set(lambda1, 2,1, gsl_matrix_get(lambda1,1,2));
        gsl_matrix_set(lambda1, 3,1, gsl_matrix_get(lambda1,1,3));
        gsl_matrix_set(lambda1, 3,2, gsl_matrix_get(lambda1,2,3));

        gsl_blas_dgemv(CblasNoTrans, 1, lambda1, &p.vector, 0, p_ph_prime );

        //printf("Lorentz Boost Matrix 0: %e,%e, %e, %e\n", gsl_matrix_get(lambda1, 0,0), gsl_matrix_get(lambda1, 0,1), gsl_matrix_get(lambda1, 0,2), gsl_matrix_get(lambda1, 0,3));
        //printf("Lorentz Boost Matrix 1: %e,%e, %e, %e\n", gsl_matrix_get(lambda1, 1,0), gsl_matrix_get(lambda1, 1,1), gsl_matrix_get(lambda1, 1,2), gsl_matrix_get(lambda1, 1,3));
        //printf("Lorentz Boost Matrix 2: %e,%e, %e, %e\n", gsl_matrix_get(lambda1, 2,0), gsl_matrix_get(lambda1, 2,1), gsl_matrix_get(lambda1, 2,2), gsl_matrix_get(lambda1, 2,3));
        //printf("Lorentz Boost Matrix 3: %e,%e, %e, %e\n", gsl_matrix_get(lambda1, 3,0), gsl_matrix_get(lambda1, 3,1), gsl_matrix_get(lambda1, 3,2), gsl_matrix_get(lambda1, 3,3));

        //printf("Before Check: %e %e %e %e\n ",gsl_vector_get(p_ph_prime, 0), gsl_vector_get(p_ph_prime, 1), gsl_vector_get(p_ph_prime, 2), gsl_vector_get(p_ph_prime, 3));

        //double check vector for 0 norm condition if photon
        if (object == 'p')
        {
            //printf("In if\n");
            boosted_p=zeroNorm(gsl_vector_ptr(p_ph_prime, 0));
        }
        else
        {
            boosted_p=gsl_vector_ptr(p_ph_prime, 0);
        }

        //printf("After Check: %e %e %e %e\n ", *(boosted_p+0),*(boosted_p+1),*(boosted_p+2),*(boosted_p+3) );
    }
    else
    {
        //printf("in else");
         //double check vector for 0 norm condition
         if (object=='p')
         {
            boosted_p=zeroNorm(p_ph);
         }
         else
         {
             //if 4 momentum isnt for photon and there is no boost to be done, we dont care about normality and just want back what was passed to lorentz boost
            boosted_p=gsl_vector_ptr(&p.vector, 0);
         }
    }
    //assign values to result
    *(result+0)=*(boosted_p+0);
    *(result+1)=*(boosted_p+1);
    *(result+2)=*(boosted_p+2);
    *(result+3)=*(boosted_p+3);

    //free up memory
    //free(boosted_p);
    gsl_matrix_free (lambda1); gsl_vector_free(p_ph_prime);
}

double *zeroNorm(double *p_ph)
{
    //ensures zero norm condition of photon 4 monetum is held
    int i=0;
    double normalizing_factor=0;
    gsl_vector_view p=gsl_vector_view_array((p_ph+1), 3); //make last 3 elements of p_ph pointer into vector

    if (*(p_ph+0) != gsl_blas_dnrm2(&p.vector ) )
    {
        normalizing_factor=(gsl_blas_dnrm2(&p.vector ));
        //printf("in zero norm if\n");
        //go through and correct 4 momentum assuming the energy is correct
        /*
        for (i=1;i<4;i++)
        {
            *(p_ph+i)= ((*(p_ph+i))/(normalizing_factor))*(*(p_ph+0));
        }
         */
        *(p_ph+1)= ((*(p_ph+1))/(normalizing_factor))*(*(p_ph+0));
        *(p_ph+2)= ((*(p_ph+2))/(normalizing_factor))*(*(p_ph+0));
        *(p_ph+3)= ((*(p_ph+3))/(normalizing_factor))*(*(p_ph+0));

    }
    /*
     if (pow((*(p_ph+0)),2) != (  pow((*(p_ph+1)),2)+pow((*(p_ph+2)),2)+pow((*(p_ph+3)),2) ) )
        {
            printf("This isnt normalized in the function\nThe difference is: %e\n", pow((*(p_ph+0)),2) - (  pow((*(p_ph+1)),2)+pow((*(p_ph+2)),2)+pow((*(p_ph+3)),2) )  );
        }
    */ //normalized within a factor of 10^-53
    return p_ph;
}

int findNearestPropertiesAndMinMFP( struct photon *ph, int num_ph, int array_num, double *time_step, double *x, double  *y, double *velx,  double *vely, double *dens_lab,\
    double *temp, double *n_dens_lab, double *n_vx, double *n_vy,double *n_temp, gsl_rng * rand)
{

    int i=0, j=0, min_index=0;
    double ph_x=0, ph_y=0, ph_phi=0, dist=0, dist_min=1e12;
    double fl_v_x=0, fl_v_y=0, fl_v_z=0; //to hold the fluid velocity in MCRaT coordinates
    double ph_v_norm=0, fl_v_norm=0;
    double n_cosangle=0, n_dens_lab_tmp=0,n_vx_tmp=0, n_vy_tmp=0, n_temp_tmp=0, rnd_tracker=0 ;

        int index=0;
        double mfp=0,min_mfp=0, beta=0;

    //go through each photon and find the blocks around it and then get the distances to all of those blocks and choose the one thats the shortest distance away
    //can optimize here, exchange the for loops and change condition to compare to each of the photons is the radius of the block is .95 (or 1.05) times the min (max) photon radius
    //or just parallelize this part here
    min_mfp=1e12;
    for (i=0;i<num_ph; i++)
    {
        //printf("%e,%e\n", ((ph+i)->r0), ((ph+i)->r1));
        ph_x=pow(pow(((ph+i)->r0),2.0)+pow(((ph+i)->r1),2.0), 0.5); //convert back to FLASH x coordinate
        ph_y=((ph+i)->r2);
        //printf("ph_x:%e, ph_y:%e\n", ph_x, ph_y);
        ph_phi=atan2(((ph+i)->r1), ((ph+i)->r0));

        dist_min=1e12;//set dist to impossible value to make sure at least first distance calulated is saved
        for(j=0;j<array_num;j++)
        {
            //if the distance between them is within 3e9, to restrict number of possible calculations,  calulate the total distance between the box and photon
            if ( (fabs(ph_x- (*(x+j)))<3e9) && (fabs(ph_y- (*(y+j)))<3e9))
            {
                //printf("In if statement\n");
                dist= pow(pow(ph_x- (*(x+j)), 2.0) + pow(ph_y- (*(y+j)) , 2.0),0.5);
                //printf("Dist calculated as: %e, index: %d\n", dist, j);
                //printf("In outer if statement, OLD: %e, %d\n", dist_min, min_index);

                if((dist<dist_min))
                {
                    //printf("In innermost if statement, OLD: %e, %d\n", dist_min, min_index);
                    dist_min=dist; //save new minimum distance
                    min_index=j; //save index
                    //printf("New Min dist: %e, New min Index: %d\n", dist_min, min_index);
                }

            }
        }
        //save values
        /*
        *(n_dens_lab+i)= (*(dens_lab+min_index));
        *(n_vx+i)= (*(velx+min_index));
        *(n_vy+i)= (*(vely+min_index));
        *(n_temp+i)= (*(temp+min_index));
        */
        (n_dens_lab_tmp)= (*(dens_lab+min_index));
        (n_vx_tmp)= (*(velx+min_index));
        (n_vy_tmp)= (*(vely+min_index));
        (n_temp_tmp)= (*(temp+min_index));

        fl_v_x=(*(velx+min_index))*cos(ph_phi);
        fl_v_y=(*(velx+min_index))*sin(ph_phi);
        fl_v_z=(*(vely+min_index));

        fl_v_norm=pow(pow(fl_v_x, 2.0)+pow(fl_v_y, 2.0)+pow(fl_v_z, 2.0), 0.5);
        ph_v_norm=pow(pow(((ph+i)->p1), 2.0)+pow(((ph+i)->p2), 2.0)+pow(((ph+i)->p3), 2.0), 0.5);

        //(*(n_cosangle+i))=((fl_v_x* ((ph+i)->p1))+(fl_v_y* ((ph+i)->p2))+(fl_v_z* ((ph+i)->p3)))/(fl_v_norm*ph_v_norm ); //find cosine of the angle between the photon and the fluid velocities via a dot product
        (n_cosangle)=((fl_v_x* ((ph+i)->p1))+(fl_v_y* ((ph+i)->p2))+(fl_v_z* ((ph+i)->p3)))/(fl_v_norm*ph_v_norm ); //make one so this is defined

        beta=pow((pow((n_vx_tmp),2)+pow((n_vy_tmp),2)),0.5);
        //printf("Beta, n_vx_tmp, n_vy_tmp: %e, %e, %e\n", beta, (n_vx_tmp), (n_vy_tmp));
        //put this in to double check that random number is between 0 and 1 (exclusive)
        rnd_tracker=0;
        while ((rnd_tracker<=0) || (rnd_tracker>=1))
        {
            rnd_tracker=gsl_rng_uniform_pos(rand);
            //printf("rnd_tracker: %e\n", rnd_tracker);
        }
        mfp=(-1)*(M_P/((n_dens_lab_tmp))/THOMP_X_SECT/(1.0-beta*((n_cosangle))))*log(rnd_tracker) ; //calulate the mfp and then multiply it by the ln of a random number to simulate distribution of mean free paths
        //printf("mfp: %e\n",mfp );
        //printf("\nPhoton: %d mfp: %e  cos_angle: %e beta: %e dens_lab: %e rnd_tracker: %e\n\n", i, mfp, n_cosangle , beta,n_dens_lab_tmp, rnd_tracker );
        if (mfp<0)
        {
            //printf("\nPhoton: %d mfp: %e  cos_angle: %e beta: %e dens_lab: %e rnd_tracker: %e\n\n", i, mfp, n_cosangle , beta,n_dens_lab_tmp, rnd_tracker );
        }

        if ( mfp<min_mfp)
        {
            min_mfp=mfp;
            *(n_dens_lab)= n_dens_lab_tmp;
            *(n_vx)= n_vx_tmp;
            *(n_vy)= n_vy_tmp;
            *(n_temp)= n_temp_tmp;
            index=i;
            //printf("new min: %e for photon %d with block properties: %e, %e, %e\n", mfp, index, n_vx_tmp, n_vy_tmp, n_temp_tmp);
        }

    }

    (*time_step)=min_mfp/C_LIGHT;
    return index;

}

void updatePhotonPosition(struct photon *ph, int num_ph, double t)
{
    //move photons by speed of light

    int i=0;
    double old_position=0, new_position=0;

    for (i=0;i<num_ph;i++)
    {
            old_position= pow(  pow(ph->r0,2)+pow(ph->r1,2)+pow(ph->r2,2), 0.5 );

            ((ph+i)->r0)+=(((ph+i)->p1)/((ph+i)->p0))*C_LIGHT*t; //update x position

            ((ph+i)->r1)+=(((ph+i)->p2)/((ph+i)->p0))*C_LIGHT*t;//update y

            ((ph+i)->r2)+=(((ph+i)->p3)/((ph+i)->p0))*C_LIGHT*t;//update z

            new_position= pow(  pow(ph->r0,2)+pow(ph->r1,2)+pow(ph->r2,2), 0.5 );

            if ((new_position-old_position)/t > C_LIGHT)
            {
                printf("PHOTON NUMBER %d IS SUPERLUMINAL. ITS SPEED IS %e c.\n", i, ((new_position-old_position)/t)/C_LIGHT);
            }
            //printf("In update  function: %e, %e, %e, %e, %e, %e, %e\n",((ph+i)->r0), ((ph+i)->r1), ((ph+i)->r2), t, ((ph+i)->p1)/((ph+i)->p0), ((ph+i)->p2)/((ph+i)->p0), ((ph+i)->p3)/((ph+i)->p0) );
    }

    //printf("In update  function: %e, %e, %e, %e\n",t, ((ph)->p1)/((ph)->p0), ((ph)->p2)/((ph)->p0), ((ph)->p3)/((ph)->p0) );

}

void photonScatter(struct photon *ph, double flash_vx, double flash_vy, double fluid_temp, gsl_rng * rand)
{
    //function to perform single photon scattering
    double ph_phi=0;
    double *ph_p=malloc(4*sizeof(double)); //pointer to hold only photon 4 momentum @ start
    double *el_p_comov=malloc(4*sizeof(double));//pointer to hold the electron 4 momenta in comoving frame
    double *ph_p_comov=malloc(4*sizeof(double));//pointer to hold the comoving photon 4 momenta
    double *fluid_beta=malloc(3*sizeof(double));//pointer to hold fluid velocity vector
    double *negative_fluid_beta=malloc(3*sizeof(double));//pointer to hold negative fluid velocity vector

    /*
    printf("%p\n", ph_p);
    ph_p=calloc(4*sizeof(double),0);
    printf("%p\n", ph_p);
    el_p_comov=calloc(4*sizeof(double),0);
    printf("%p\n", el_p_comov);
    ph_p_comov=calloc(4*sizeof(double),0);
    printf("%p\n", ph_p_comov);
    fluid_beta=calloc(3*sizeof(double),0);
    printf("%p\n", fluid_beta);
    negative_fluid_beta=calloc(3*sizeof(double),0);
    printf("Done calling calloc\n");
     */
    ph_phi=atan2((ph->r1), ((ph->r0)));
    //printf("ph_phi=%e\n", ph_phi);
    //convert flash coordinated into MCRaT coordinates
    //printf("Getting fluid_beta\n");

    (*(fluid_beta+0))=flash_vx*cos(ph_phi);
    (*(fluid_beta+1))=flash_vx*sin(ph_phi);
    (*(fluid_beta+2))=flash_vy;

    //fill in photon 4 momentum
    //printf("filling in 4 momentum in photonScatter\n");
    *(ph_p+0)=(ph->p0);
    *(ph_p+1)=(ph->p1);
    *(ph_p+2)=(ph->p2);
    *(ph_p+3)=(ph->p3);
    //printf("Unscattered Photon in Lab frame: %e, %e, %e,%e\n", *(ph_p+0), *(ph_p+1), *(ph_p+2), *(ph_p+3));
    //printf("Fluid Beta: %e, %e, %e\n", *(fluid_beta+0),*(fluid_beta+1), *(fluid_beta+2));

    //first we bring the photon to the fluid's comoving frame
    lorentzBoost(fluid_beta, ph_p, ph_p_comov, 'p');
    //printf("Old: %e, %e, %e,%e\n", ph->p0, ph->p1, ph->p2, ph->p3);
    //printf("Before Scattering, In Comov_frame:\n");
    //printf("ph_comov: %e, %e, %e,%e\n", *(ph_p_comov+0), *(ph_p_comov+1), *(ph_p_comov+2), *(ph_p_comov+3));

    //second we generate a thermal electron at the correct temperature
    singleElectron(el_p_comov, fluid_temp, ph_p_comov, rand);
    //printf("el_comov: %e, %e, %e,%e\n", *(el_p_comov+0), *(el_p_comov+1), *(el_p_comov+2), *(el_p_comov+3));

    //third we perform the scattering and save scattered photon 4 monetum in ph_p_comov @ end of function
    singleComptonScatter(el_p_comov, ph_p_comov, rand);
    //printf("After Scattering, After Lorentz Boost to Comov frame: %e, %e, %e,%e\n", *(ph_p_comov+0), *(ph_p_comov+1), *(ph_p_comov+2), *(ph_p_comov+3));

    //fourth we bring the photon back to the lab frame
     *(negative_fluid_beta+0)=-1*( *(fluid_beta+0));
     *(negative_fluid_beta+1)=-1*( *(fluid_beta+1));
     *(negative_fluid_beta+2)=-1*( *(fluid_beta+2));
    lorentzBoost(negative_fluid_beta, ph_p_comov, ph_p, 'p');
    //printf("Scattered Photon in Lab frame: %e, %e, %e,%e\n", *(ph_p+0), *(ph_p+1), *(ph_p+2), *(ph_p+3));
    //printf("Old: %e, %e, %e,%e\n", ph->p0, ph->p1, ph->p2, ph->p3);
    //printf("Old: %e, %e, %e,%e\n", *(ph_p_comov+0), *(ph_p_comov+1), *(ph_p_comov+2), *(ph_p_comov+3));

    //assign the photon its new lab 4 momentum
    (ph->p0)=(*(ph_p+0));
    (ph->p1)=(*(ph_p+1));
    (ph->p2)=(*(ph_p+2));
    (ph->p3)=(*(ph_p+3));
    //printf("Done assigning values to original struct\n");

    free(el_p_comov);
    //printf("done here\n");
    free(ph_p_comov);
    //printf("done here\n");
    free(fluid_beta); // ?maybe not? getting an error checksum for freed object - object was probably modified after being freed.
    //printf("done here\n");
    free(negative_fluid_beta);
    //printf("done here\n");
    free(ph_p);
    //printf("done here\n");
    ph_p=NULL;negative_fluid_beta=NULL;ph_p_comov=NULL; el_p_comov=NULL;
}

void singleElectron(double *el_p, double temp, double *ph_p, gsl_rng * rand)
{
    //generates an electron with random energy

    double factor=0, gamma=0;
    double y_dum=0, f_x_dum=0, x_dum=0, beta_x_dum=0, beta=0, phi=0, theta=0, ph_theta=0, ph_phi=0;
    gsl_matrix *rot= gsl_matrix_calloc (3, 3); //create matrix thats 3x3 to do rotation
    gsl_vector_view el_p_prime ; //create vector to hold rotated electron 4 momentum
    gsl_vector *result=gsl_vector_alloc (3);

    //printf("Temp in singleElectron: %e\n", temp);
    if (temp>= 1e7)
    {
        //printf("In if\n");
        factor=K_B*temp/(M_EL*pow(C_LIGHT,2.0));
        y_dum=1; //initalize loop to get a random gamma from the distribution of electron velocities
        f_x_dum=0;
        while ((isnan(f_x_dum) !=0) || (y_dum>f_x_dum) )
        {

            x_dum=gsl_rng_uniform_pos(rand)*(1+100*factor);
            beta_x_dum=pow(1-(pow(x_dum, -2.0)) ,0.5);
            y_dum=gsl_rng_uniform(rand)/2.0;

            f_x_dum=pow(x_dum,2)*(beta_x_dum/gsl_sf_bessel_Kn (2, 1.0/factor))*exp(-1*x_dum/factor); //not sure if this is right is giving small values of gamma -> beta=nan
            //printf("Choosing a Gamma: xdum: %e, f_x_dum: %e, y_dum: %e\n", x_dum, f_x_dum, y_dum);
        }
        gamma=x_dum;

    }
    else
    {

        //printf("In else\n");
        factor=pow(K_B*temp/M_EL,0.5);
        //calculate a random gamma from 3 random velocities drawn from a gaussian distribution with std deviation of "factor"
        gamma=pow( 1- (pow(gsl_ran_gaussian(rand, factor)/C_LIGHT, 2)+ pow(gsl_ran_gaussian(rand, factor)/C_LIGHT, 2)+pow(gsl_ran_gaussian(rand, factor)/C_LIGHT, 2)  ) ,-0.5);
    }

    //printf("Chosen Gamma: %e\n",gamma);

    beta=pow( 1- (1/pow(  gamma,2.0 ))  ,0.5);
    //printf("Beta is: %e in singleElectron\n", beta);
    phi=gsl_rng_uniform(rand)*2*M_PI;

    y_dum=1; //initalize loop to get a random theta
    f_x_dum=0;
    while (y_dum>f_x_dum)
    {
        y_dum=gsl_rng_uniform(rand)*1.3;
        x_dum=gsl_rng_uniform(rand)*M_PI;
        f_x_dum=sin(x_dum)*(1-(beta*cos(x_dum)));
    }
    theta=x_dum;
    //printf("Beta: %e\tPhi: %e\tTheta: %e\n",beta,phi, theta);
    //fill in electron 4 momentum NOT SURE WHY THE ORDER IS AS SUCH SEEMS TO BE E/c, pz,py,px!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    *(el_p+0)=gamma*(M_EL)*(C_LIGHT);
    *(el_p+1)=gamma*(M_EL)*(C_LIGHT)*beta*cos(theta);
    *(el_p+2)=gamma*(M_EL)*(C_LIGHT)*beta*sin(theta)*sin(phi);
    *(el_p+3)=gamma*(M_EL)*(C_LIGHT)*beta*sin(theta)*cos(phi);

    //printf("Old: %e, %e, %e,%e\n", *(el_p+0), *(el_p+1), *(el_p+2), *(el_p+3));

    el_p_prime=gsl_vector_view_array((el_p+1), 3);

    //find angles of photon NOT SURE WHY WERE CHANGING REFERENCE FRAMES HERE???!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ph_phi=atan2(*(ph_p+2), *(ph_p+3)); //Double Check
    ph_theta=atan2(pow( pow(*(ph_p+2),2)+  pow(*(ph_p+3),2) , 0.5) , (*(ph_p+1)) );

    //printf("Calculated Photon phi and theta in singleElectron:%e, %e\n", ph_phi, ph_theta);

    //fill in rotation matrix to rotate around x axis to get rid of phi angle
    gsl_matrix_set(rot, 1,1,1);
    gsl_matrix_set(rot, 2,2,cos(ph_theta));
    gsl_matrix_set(rot, 0,0,cos(ph_theta));
    gsl_matrix_set(rot, 0,2,-sin(ph_theta));
    gsl_matrix_set(rot, 2,0,sin(ph_theta));
    gsl_blas_dgemv(CblasNoTrans, 1, rot, &el_p_prime.vector, 0, result);

    //printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot, 0,0), gsl_matrix_get(rot, 0,1), gsl_matrix_get(rot, 0,2));
    //printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot, 1,0), gsl_matrix_get(rot, 1,1), gsl_matrix_get(rot, 1,2));
    //printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot, 2,0), gsl_matrix_get(rot, 2,1), gsl_matrix_get(rot, 2,2));

    //printf("Middle: %e, %e, %e,%e\n", *(el_p+0), gsl_vector_get(result,0), gsl_vector_get(result,1), gsl_vector_get(result,2));

    gsl_matrix_set_all(rot,0);

    gsl_matrix_set(rot, 0,0,1);
    gsl_matrix_set(rot, 1,1,cos(-ph_phi));
    gsl_matrix_set(rot, 2,2,cos(-ph_phi));
    gsl_matrix_set(rot, 1,2,-sin(-ph_phi));
    gsl_matrix_set(rot, 2,1,sin(-ph_phi));
    gsl_blas_dgemv(CblasNoTrans, 1, rot, result, 0, &el_p_prime.vector);
    /*
    printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot, 0,0), gsl_matrix_get(rot, 0,1), gsl_matrix_get(rot, 0,2));
    printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot, 1,0), gsl_matrix_get(rot, 1,1), gsl_matrix_get(rot, 1,2));
    printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot, 2,0), gsl_matrix_get(rot, 2,1), gsl_matrix_get(rot, 2,2));
    printf("Final EL_P_vec: %e, %e, %e,%e\n", *(el_p+0), gsl_vector_get(&el_p_prime.vector,0), gsl_vector_get(&el_p_prime.vector,1), gsl_vector_get(&el_p_prime.vector,2));
    */

    //gsl_rng_free (rand);
    //printf("freeing pointers in singleElectron\n");
    gsl_matrix_free (rot);gsl_vector_free(result);
    //printf("Done freeing pointers in singleElectron\n");
}

void singleComptonScatter(double *el_comov, double *ph_comov, gsl_rng * rand)
{
    //This routine performs a Compton scattering between a photon and a moving electron.
    int i=0;
    double *el_v=malloc(3*sizeof(double));
    double *negative_el_v=malloc(3*sizeof(double));
    double *ph_p_prime=malloc(4*sizeof(double));//use this to keep track of how the ph 4 momentum changes with each rotation
    double *el_p_prime=malloc(4*sizeof(double));
    double phi0=0, phi1=0, phi=0, theta=0;
    double y_dum, f_x_dum, x_dum;
    gsl_matrix *rot0= gsl_matrix_calloc (3, 3); //create matricies thats 3x3 to do rotations
    gsl_matrix *rot1= gsl_matrix_calloc (3, 3);
    gsl_vector *result0=gsl_vector_alloc (3); //vectors to hold results of rotations
    gsl_vector *result1=gsl_vector_alloc (3);
    gsl_vector *result=gsl_vector_alloc (4);
    gsl_vector *whole_ph_p=gsl_vector_alloc (4);
    gsl_vector_view ph_p ; //create vector to hold comoving photon and electron 4 momentum
    gsl_vector_view el_p ;

    //fill in electron velocity array and photon 4 momentum
    *(el_v+0)=(*(el_comov+1))/(*(el_comov+0));
    *(el_v+1)=(*(el_comov+2))/(*(el_comov+0));
    *(el_v+2)=(*(el_comov+3))/(*(el_comov+0));
    //printf("el_v: %e, %e, %e\n", *(el_v+0), *(el_v+1), *(el_v+2));

    //lorentz boost into frame where the electron is stationary
    lorentzBoost(el_v, el_comov, el_p_prime, 'e');
    lorentzBoost(el_v, ph_comov, ph_p_prime, 'p');
    //printf("New ph_p in electron rest frame: %e, %e, %e,%e\n", *(ph_p_prime+0), *(ph_p_prime+1), *(ph_p_prime+2), *(ph_p_prime+3));

    ph_p=gsl_vector_view_array((ph_p_prime+1), 3);
    el_p=gsl_vector_view_array(el_p_prime,4);
    phi0=atan2(*(ph_p_prime+2), *(ph_p_prime+1) );
    //printf("Photon Phi: %e\n", phi0);
    //rotate the axes so that the photon incomes along the x-axis
    gsl_matrix_set(rot0, 2,2,1);
    gsl_matrix_set(rot0, 0,0,cos(-phi0));
    gsl_matrix_set(rot0, 1,1,cos(-phi0));
    gsl_matrix_set(rot0, 0,1,-sin(-phi0));
    gsl_matrix_set(rot0, 1,0,sin(-phi0));
    gsl_blas_dgemv(CblasNoTrans, 1, rot0, &ph_p.vector, 0, result0);

    /*
    printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot0, 0,0), gsl_matrix_get(rot0, 0,1), gsl_matrix_get(rot0, 0,2));
    printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot0, 1,0), gsl_matrix_get(rot0, 1,1), gsl_matrix_get(rot0, 1,2));
    printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot0, 2,0), gsl_matrix_get(rot0, 2,1), gsl_matrix_get(rot0, 2,2));
    */

    //set values of ph_p_prime equal to the result and get new phi from result
    *(ph_p_prime+1)=gsl_vector_get(result0,0);
    *(ph_p_prime+2)=0;//gsl_vector_get(result,1); //just directly setting it to 0 now?
    *(ph_p_prime+3)=gsl_vector_get(result0,2);

    phi1=atan2(gsl_vector_get(result0,2), gsl_vector_get(result0,0));

    /*
    printf("rotation 1: %e, %e, %e\n",  *(ph_p_prime+1),  *(ph_p_prime+2),  *(ph_p_prime+3));
    printf("Photon Phi: %e\n", phi1);
    printf("make sure the vector view is good: %e, %e, %e,%e\n", *(ph_p_prime+0), gsl_vector_get(&ph_p.vector,0), gsl_vector_get(&ph_p.vector,1), gsl_vector_get(&ph_p.vector,2));
    */

    //rotate around y to bring it all along x
    gsl_matrix_set(rot1, 1,1,1);
    gsl_matrix_set(rot1, 0,0,cos(-phi1));
    gsl_matrix_set(rot1, 2,2,cos(-phi1));
    gsl_matrix_set(rot1, 0,2,-sin(-phi1));
    gsl_matrix_set(rot1, 2,0,sin(-phi1));
    gsl_blas_dgemv(CblasNoTrans, 1, rot1, &ph_p.vector, 0, result1);

    /*
    printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot1, 0,0), gsl_matrix_get(rot1, 0,1), gsl_matrix_get(rot1, 0,2));
    printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot1, 1,0), gsl_matrix_get(rot1, 1,1), gsl_matrix_get(rot1, 1,2));
    printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot1, 2,0), gsl_matrix_get(rot1, 2,1), gsl_matrix_get(rot1, 2,2));
    */

    //set values of ph_p_prime equal to the result and get new phi from result
    *(ph_p_prime+1)=*(ph_p_prime+0);//why setting it to the energy?
    *(ph_p_prime+2)=gsl_vector_get(result1,1);
    *(ph_p_prime+3)=0; //just directly setting it to 0 now?

    //printf("rotation 2: %e, %e, %e, %e\n",  *(ph_p_prime+0), *(ph_p_prime+1),  *(ph_p_prime+2),  *(ph_p_prime+3));

    //generate random theta and phi angles for scattering
    phi=gsl_rng_uniform(rand)*2*M_PI;
    //printf("Phi: %e\n", phi);

    y_dum=1; //initalize loop to get a random theta
    f_x_dum=0;
    while (y_dum>f_x_dum)
    {
        y_dum=gsl_rng_uniform(rand)*1.09;
        x_dum=gsl_rng_uniform(rand)*M_PI;
        f_x_dum=sin(x_dum)*(1+pow(cos(x_dum),2));
    }
    theta=x_dum;

    //printf("Theta: %e\n", theta);

    //perform scattering and compute new 4-momenta of electron and photon
    //scattered photon 4 momentum
    gsl_vector_set(result, 0, (*(ph_p_prime+0))/(1+ (( (*(ph_p_prime+0))*(1-cos(theta)) )/(M_EL*C_LIGHT )) ) ); //DOUBLE CHECK HERE!!!!
    gsl_vector_set(result, 1, gsl_vector_get(result,0)*cos(theta) );
    gsl_vector_set(result, 2, gsl_vector_get(result,0)*sin(theta)*sin(phi) );
    gsl_vector_set(result, 3, gsl_vector_get(result,0)*sin(theta)*cos(phi) );
    //printf("%e\n", gsl_vector_get(result,0));

    //calculate electron 4 momentum OPTIMIZE HERE: DONT USE A FOR LOOP HERE!!!! Done
    //prescattered photon 4 momentum
    gsl_vector_set(whole_ph_p, 0, (*(ph_p_prime+0)));
    gsl_vector_set(whole_ph_p, 1, (*(ph_p_prime+1)));
    gsl_vector_set(whole_ph_p, 2, (*(ph_p_prime+2)));
    gsl_vector_set(whole_ph_p, 3, (*(ph_p_prime+3)));
    /*
    for (i=0;i<4;i++)
    {
        gsl_vector_set(whole_ph_p, i, (*(ph_p_prime+i)));
    }
    */
     gsl_vector_sub(whole_ph_p,result); //resut is saved into ph_p vector, unscattered-scattered 4 mometum of photon
    gsl_vector_add(&el_p.vector ,whole_ph_p);

    /*
    printf("After scattering:\n");
    printf("el_p: %e, %e, %e,%e\n", gsl_vector_get(&el_p.vector,0), gsl_vector_get(&el_p.vector,1), gsl_vector_get(&el_p.vector,2), gsl_vector_get(&el_p.vector,3));
    printf("ph_p: %e, %e, %e,%e\n", gsl_vector_get(result,0), gsl_vector_get(result,1), gsl_vector_get(result,2), gsl_vector_get(result,3));
    */

    //rotate back to comoving frame
    *(ph_p_prime+0)=gsl_vector_get(result,0);
    *(ph_p_prime+1)=gsl_vector_get(result,1); //set values of photon prime momentum from doing the scattering to use the vector view of it in dot product
    *(ph_p_prime+2)=gsl_vector_get(result,2);
    *(ph_p_prime+3)=gsl_vector_get(result,3);
    gsl_matrix_set_all(rot1,0);
    gsl_matrix_set(rot1, 1,1,1);
    gsl_matrix_set(rot1, 0,0,cos(-phi1));
    gsl_matrix_set(rot1, 2,2,cos(-phi1));
    gsl_matrix_set(rot1, 0,2,sin(-phi1));
    gsl_matrix_set(rot1, 2,0,-sin(-phi1));
    gsl_blas_dgemv(CblasNoTrans, 1, rot1, &ph_p.vector, 0, result1);

    /*
    printf("Photon Phi: %e\n", phi1);
    printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot1, 0,0), gsl_matrix_get(rot1, 0,1), gsl_matrix_get(rot1, 0,2));
    printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot1, 1,0), gsl_matrix_get(rot1, 1,1), gsl_matrix_get(rot1, 1,2));
    printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot1, 2,0), gsl_matrix_get(rot1, 2,1), gsl_matrix_get(rot1, 2,2));
    */

    //set values of ph_p_prime to result1 from undoing 2nd rotation
    *(ph_p_prime+1)=gsl_vector_get(result1,0);
    *(ph_p_prime+2)=gsl_vector_get(result1,1);
    *(ph_p_prime+3)=gsl_vector_get(result1,2);
    //printf("Undo rotation 2: %e, %e, %e, %e\n",  *(ph_p_prime+0), *(ph_p_prime+1),  *(ph_p_prime+2),  *(ph_p_prime+3));
    //ignore the electron, dont care about it, undo the first rotation
    gsl_matrix_set_all(rot0,0);
    gsl_matrix_set(rot0, 2,2,1);
    gsl_matrix_set(rot0, 0,0,cos(-phi0));
    gsl_matrix_set(rot0, 1,1,cos(-phi0));
    gsl_matrix_set(rot0, 0,1,sin(-phi0));
    gsl_matrix_set(rot0, 1,0,-sin(-phi0));
    gsl_blas_dgemv(CblasNoTrans, 1, rot0, &ph_p.vector, 0, result0);

    /*
    printf("Photon Phi: %e\n", phi0);
    printf("Rotation Matrix 0: %e,%e, %e\n", gsl_matrix_get(rot0, 0,0), gsl_matrix_get(rot0, 0,1), gsl_matrix_get(rot0, 0,2));
    printf("Rotation Matrix 1: %e,%e, %e\n", gsl_matrix_get(rot0, 1,0), gsl_matrix_get(rot0, 1,1), gsl_matrix_get(rot0, 1,2));
    printf("Rotation Matrix 2: %e,%e, %e\n", gsl_matrix_get(rot0, 2,0), gsl_matrix_get(rot0, 2,1), gsl_matrix_get(rot0, 2,2));
    */

    *(ph_p_prime+1)=gsl_vector_get(result0,0);
    *(ph_p_prime+2)=gsl_vector_get(result0,1);
    *(ph_p_prime+3)=gsl_vector_get(result0,2);
    //printf("Undo rotation 1: %e, %e, %e, %e\n",  *(ph_p_prime+0), *(ph_p_prime+1),  *(ph_p_prime+2),  *(ph_p_prime+3));
    //deboost photon to lab frame
    *(negative_el_v+0)=(-1*(*(el_v+0)));
    *(negative_el_v+1)=(-1*(*(el_v+1)));
    *(negative_el_v+2)=(-1*(*(el_v+2)));

    lorentzBoost(negative_el_v, ph_p_prime, ph_comov, 'p');
    //printf("Undo boost 1: %e, %e, %e, %e\n",  *(ph_comov+0), *(ph_comov+1),  *(ph_comov+2),  *(ph_comov+3));

    gsl_matrix_free(rot0); gsl_matrix_free(rot1);gsl_vector_free(result0);gsl_vector_free(result1);gsl_vector_free(result);
    //gsl_rng_free (rand);
    gsl_vector_free(whole_ph_p);free(ph_p_prime);free(el_p_prime);free(el_v); free(negative_el_v);
}

double averagePhotonEnergy(struct photon *ph, int num_ph)
{
    //to calculate average photon energy
    int i=0;
    double sum=0;
    for (i=0;i<num_ph;i++)
    {
        sum+=((ph+i)->p0);
    }

    return (sum*C_LIGHT)/num_ph;
}

void phScattStats(struct photon *ph, int ph_num, int *max, int *min, double *avg )
{
    int temp_max=0, temp_min=-1,  i=0;
    double sum=0;

    for (i=0;i<ph_num;i++)
    {
        sum+=((ph+i)->num_scatt);

        if (((ph+i)->num_scatt) > temp_max )
        {
            temp_max=((ph+i)->num_scatt);
            //printf("The new max is: %d\n", temp_max);
        }

        if ((i==0) || (((ph+i)->num_scatt)<temp_min))
        {
            temp_min=((ph+i)->num_scatt);
            //printf("The new min is: %d\n", temp_min);
        }
    }

    *avg=sum/ph_num;
    *max=temp_max;
    *min=temp_min;

}

void cylindricalPrep(double *gamma, double *vx, double *vy, double *dens, double *dens_lab, double *pres, double *temp, int num_array)
{
    double  gamma_infinity=100, t_comov=1*pow(10, 5), ddensity=3e-9;// the comoving temperature in Kelvin, and the comoving density in g/cm^2
    int i=0;
    double vel=pow(1-pow(gamma_infinity, -2.0) ,0.5), lab_dens=gamma_infinity*ddensity;

    for (i=0; i<num_array;i++)
    {
        *(gamma+i)=gamma_infinity;
        *(vx+i)=0;
        *(vy+i)=vel;
        *(dens+i)=ddensity;
        *(dens_lab+i)=lab_dens;
        *(pres+i)=(A_RAD*pow(t_comov, 4.0))/(3*pow(C_LIGHT, 2.0));
        *(temp+i)=pow(3*(*(pres+i))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0); //just assign t_comov
    }

}

void sphericalPrep(double *r,  double *x, double *y, double *gamma, double *vx, double *vy, double *dens, double *dens_lab, double *pres, double *temp, int num_array)
{
    double  gamma_infinity=100, lumi=1e52, r00=1e8;
    double vel=0;
    int i=0;

    for (i=0;i<num_array;i++)
    {
        if ((*(r+i)) >= (r00*gamma_infinity))
        {
            *(gamma+i)=gamma_infinity;
            *(pres+i)=(lumi*pow(r00, 2.0/3.0)*pow(*(r+i), -8.0/3.0) )/(12.0*M_PI*C_LIGHT*pow(gamma_infinity, 4.0/3.0)*pow(C_LIGHT, 2.0));
        }
        else
        {
            *(gamma+i)=(*(r+i))/r00;
            *(pres+i)=(lumi*pow(r00, 2.0))/(12.0*M_PI*C_LIGHT*pow(C_LIGHT, 2.0)*pow(*(r+i), 4.0) );
        }

        vel=pow(1-pow(*(gamma+i), -2.0) ,0.5);
        *(vx+i)=(vel*(*(x+i)))/pow(pow(*(x+i), 2)+ pow(*(y+i), 2) ,0.5);
        *(vy+i)=(vel*(*(y+i)))/pow(pow(*(x+i), 2)+ pow(*(y+i), 2) ,0.5);
        *(dens+i)=lumi/(4*M_PI*pow(*(r+i), 2.0)*pow(C_LIGHT, 3.0)*gamma_infinity*(*(gamma+i)));
        *(dens_lab+i)=*(dens+i)*(*(gamma+i));
    }

}

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 * Copyright 2018 - Matteo Ragni, Matteo Cocetti - University of Trento
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef LIBEULER_H_
#define LIBEULER_H_

#include <cblas.h>
#include "libnewton.h" /**< Custom solver for implicit step */

/**
 * @brief Callback for the ode vector field
 * The callback for the ODE model. The callback stores the output in
 * the first pointer (f). The function does not need to allocate nor free
 * the input vector, but may overflow if exceedes the dimension that is
 * declared in the euler options structure.
 * When accessing an element of the parameter, please remember that the
 * implementation keeps in mind a MATLAB-like interface. If the parameters
 * (in MATLAB) are passed like:
 * @code
 * f(t, x, u, p1, [p2, p3])
 * @endcode
 * then, in the C code, p1, p2 and p3 may be accessed in p as:
 * @code
 * rk4_float p1, p2, p3;
 * p1 = p[0][0];
 * p2 = p[1][0];
 * p3 = p[1][1];
 * @endcode
 * If you need the time step **inside** the callback, you should directly define
 * it inside or pass it as a parameter (and remeber to update it in the euler_options
 * structure).
 * @param f otput vector for the ODE
 * @param t current time for the function
 * @param x state for the ODE evaluation
 * @param u external input for the ODE evaluation
 * @param p pointer to arrays of parameters
 * @param data auxiliary data pointer to void for user data
 * @returns nothing
 * @warning The callback may produce a buffer overflow error if tries to
 *          write more than x_size elements in xdot.
 */
typedef void (*euler_ode_function)(
    double *f,
    const double t,
    const double *x,
    const double *u,
    const double **p,
    void *data);

/**
 * @brief Callback for the ode vector field Jacobian
 * The callback for the ODE Jacobian. The callback stores the output in
 * the first pointer (df). The function does not need to allocate nor free
 * the input vector, but may overflow if exceedes the dimension that is
 * declared in the euler options structure. The matrix is stored continuosly
 * inan array, with ordering as specified in options structure.
 * @param f otput vector for the ODE
 * @param t current time for the function
 * @param x state for the ODE evaluation
 * @param u external input for the ODE evaluation
 * @param p pointer to arrays of parameters
 * @param data auxiliary data pointer to void for user data
 * @returns nothing
 * @warning The callback may produce a buffer overflow error if tries to
 *          write more than x_size elements in xdot.
 */
typedef void (*euler_ode_jacobian)(
    double *f,
    const double t,
    const double *x,
    const double *u,
    const double **p,
    void *data);

/**
 * @brief Returning value for the integrator
 */
typedef enum euler_ret {
  EULER_SUCCESS = 0,  /**< Correct execution */
  EULER_EMALLOC,      /**< Memory allocation error */
  EULER_NULLPTR,      /**< Received a null pointer */
  EULER_GENERIC       /**< Generic error raised. @todo: should be divided for Newton error */
} euler_ret;

typedef struct euler_options {
  double ts;               /**< Integration step */
  double alpha;            /**< Tustin transform coefficient. \f$\alpha \in [0,1]\f$.
                                For \f$\alpha = 0\f$, no optimization is performed */
  lapack_int x_size;       /**< State dimensions. Will also be vector field size */
  lapack_int u_offset;     /**< Input offset for implicit steps. It can also be 0 */
  lapack_int ordering;     /**< Ordering: LAPACK_COL_MAJOR or LAPACK_ROW_MAJOR */
  double s_tol;            /**< Tolerance for Newton solution */
  double x_tol;            /**< Step tolerance for Newton solution */
  lapack_int max_iter;     /**< Maximum number of iterations for Newton solver */
  euler_ode_function f;    /**< Actual ODE vector field */
  euler_ode_jacobian df;   /**< Jacobian of the vector field */
  void *data;              /**< Empty space for user data */
} euler_options;

/**
 * @brief Euler step (explicit or implicit)
 *
 * The function performs the euler step.
 * If \f$ \alpha = 0\f$ the step performed is:
 * \f{
 *   x(t+h) = x(t) + f(x(t), u_{1..dim(u)}, p)
 * \f}
 * while for \f$\alpha \in (0,1]\f$ the step perfomed is the solution of
 * the following non linear problem:
 * \f{
 *   x(t+h) = x(t) + (1-\alpha) h f(x(t), u_{1..u_{off}, p} +
 *            \alpha h f(x(t+h), u_{u_{off}..dim(u)}, p)
 * \f}
 * which is an implicit integration step. The initial guess for the solution of the
 * non linear problem is the state \f$x(t)\f$.
 * @param opt pointer to scruct with options
 * @param xp next integration step
 * @param t current integration time
 * @param x current state
 * @param u control vector. May contain both current and next control input,
 *          in sequence: \f$u = [u(t), u(t+h)]\f$. The offset allows to get as
 *          input in the callback the next input.
 * @param p  pointer to arrays of parameters
 * @param data void pointer to userspace data
 * @return an exit code to check if integration step succeeded
 */
euler_ret euler(
  const euler_options *opt,
  double *xp,
  const double t,
  const double *x,
  const double *u,
  const double **p,
  void *data);

#endif /* LIBEULER_H_ */

ext: h

max_stars_repo_name: MatteoRagni/libeuler

max_issues_repo_name: MatteoRagni/libeuler

max_forks_repo_name: MatteoRagni/libeuler


text:
#ifndef FILTER_H
#define FILTER_H

#include <gsl/gsl_vector.h>

#include <utility>

void filter_fir(gsl_vector *b, gsl_vector *x, gsl_vector *res);

void filter_iir(gsl_vector *b, gsl_vector *a, gsl_vector *x, gsl_vector *res);

void filter_hpf(gsl_vector *x, const double fc);

void filter_lpf(gsl_vector *x, const double fc);

#endif // FILTER_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
#ifndef WORKING_DIRECTORY_INCLUDE
#define WORKING_DIRECTORY_INCLUDE

#include <string>

#include <gsl/string_span>

#include "config/commandLineOptions.h"
#include "config/variablesMap.h"
#include "core/task.h"

namespace execHelper {
namespace config {
class FleetingOptionsInterface;
} // namespace config
} // namespace execHelper

namespace execHelper {
namespace plugins {
using WorkingDir = std::string;
const gsl::czstring<> WORKING_DIR_KEY = "working-dir";

/**
 * \brief Extends the functionality to add the working directory configuration variables of the process to execute, using the _working-dir_ key in the configuration
 */
struct WorkingDirectory {
    /**
     * Adds the variables for this functionality to the given variables map
     */
    inline static void
    getVariables(config::VariablesMap& /*variables*/,
                 const config::FleetingOptionsInterface& /*options*/) noexcept {
        ;
    }

    /*! @copydoc JobsLong::apply(core::Task&, const config::VariablesMap&)
     */
    inline static void apply(core::Task& task,
                             const config::VariablesMap& variables) noexcept {
        auto workingDir = variables.get<WorkingDir>(WORKING_DIR_KEY);
        if(workingDir) {
            task.setWorkingDirectory(*workingDir);
        }
    }
};
} // namespace plugins
} // namespace execHelper

#endif /* WORKING_DIRECTORY_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2017 by Synge Todo <wistaria@phys.s.u-tokyo.ac.jp>
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#include <math.h>
#include <cblas.h>
#include <lapacke.h>
#include <rokko/cmatrix.h>

int imin(int x, int y) { return (x < y) ? x : y; }

int main(int argc, char** argv) {
  int n = 5;
  int i, j;
  float norm, check;
  float **a;
  float *work;

  if (argc > 1) n = atoi(argv[1]);

  // generate matrix and rhs vector
  a = alloc_smatrix(n, n);
  for (j = 0; j < n; ++j) {
    for (i = 0; i < n; ++i) {
      mat_elem(a, i, j) = imin(i+1, j+1);
    }
  }
  printf("Matrix A: ");
  fprint_smatrix(stdout, n, n, a);
  work = alloc_svector(n);

  /* calculate various norms */
  norm = LAPACKE_slange_work(LAPACK_COL_MAJOR, 'M', n, n, mat_ptr(a), n, vec_ptr(work));
  printf("element of largest absolute value = %e\n", norm);

  norm = LAPACKE_slange(LAPACK_COL_MAJOR, '1', n, n, mat_ptr(a), n);
  printf("one norm = %e\n", norm);

  norm = LAPACKE_slange(LAPACK_COL_MAJOR, 'I', n, n, mat_ptr(a), n);
  printf("infinity norm = %e\n", norm);

  norm = LAPACKE_slange(LAPACK_COL_MAJOR, 'F', n, n, mat_ptr(a), n);
  printf("Frobenius norm = %e\n", norm);

  /* check of result */
  check = 0;
  for (j = 0; j < n; ++j) {
    for (i = 0; i < n; ++i) {
      check += mat_elem(a, i, j) * mat_elem(a, i, j);
    }
  }
  if (fabsf(check - norm * norm) > 1e-6 * check) {
    fprintf(stderr, "Error: check error %e != %e\n", norm * norm, check);
    exit(255);
  }

  free_smatrix(a);
  return 0;
}

ext: c

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
/* from lapack-netlib, 3-clause BSD licensed */
/* modified (C) 2016, A. Boyle */

#include "config.h"
#include <stdio.h>
#include <math.h> /* for fabs */
#include <lapacke.h>

int main (int argc, const char * argv[])
{
    double A[] = {1, 1, 1,
                  2, 3, 4,
                  3, 5, 2,
                  4, 2, 5,
                  5, 4, 3
                 };
    double b[] = { -10, -3,
                   12, 14,
                   14, 12,
                   16, 16,
                   18, 16
                 };
    int i, j;
    /* Solve least squares problem: DGELS using row-major layout */
    int info = LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', 5, 3, 2, A, 3, b, 2);
    if (info) {
        printf("info=%d\n", info);
        exit(EXIT_FAILURE);
    }
    /* matlab: expect = A \ b */
    double const expect[] = {2, 1,
                             1, 1,
                             1, 2
                            };
    for(i = 0; i < 3; i++) {
        for(j = 0; j < 2; j++) {
            printf(" %g", b[i * 2 + j]);
        }
        printf("\n");
    }
    for(i = 0; i < 6; i++) {
        if(fabs(expect[i] - b[i]) > 4e-15) {
            printf("%d %g\n", i, fabs(expect[i] - b[i]));
            printf("FAIL\n");
            exit(EXIT_FAILURE);
        }
    }
    printf("PASS\n");
    exit(EXIT_SUCCESS);
}

max_stars_repo_name: boyle/zedhat

max_issues_repo_name: boyle/zedhat

max_forks_repo_name: boyle/zedhat


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <mpi.h>
#include "parmt_mtsearch.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "parmt_utils.h"
#include "cps_mpe.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

int parmt_obsSearch64f(const MPI_Comm globalComm,
                       const MPI_Comm obsComm, const MPI_Comm locComm,
                       const bool linObsComm, const bool linLocComm,
                       const struct parmtGeneralParms_struct parms,
                       struct parmtData_struct data,
                       struct localMT_struct mtloc,
                       double *__restrict__ phi)
{
    MPI_Win counterWin;
    const char *archiveName = NULL;
    double lagTime, *CeInv, *phiWork, *phiAll, *var, *varAll, *varOut, *wts,
           defaultWeight, wtNorm;
    int *lagsWork, *nlags, ierr, iobs, myid, nobsProcs,
        npmax, value, value0;
    bool ldefault, lwantLags;
    const double half = 0.5;
    const int master = 0;
    //------------------------------------------------------------------------//
    //
    // Have the root process determine lag information
    ierr = 0;
    MPI_Comm_rank(globalComm, &myid);
    nlags = memory_calloc32i(data.nobs);
    nobsProcs = 0;
    lagsWork = NULL;
    phiWork = NULL;
    phiAll = NULL;
    CeInv = NULL;
    var = NULL;
    varAll = NULL;
    varOut = NULL;
    wts = NULL;
    npmax = 0;
    lwantLags = false;
    if (myid == master)
    {
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            nlags[iobs] = 0;
            if (parms.lwantLags)
            {
                //lagTime = data.data[iobs].header.user0;
                //if (lagTime < zero){lagTime = parms.defaultMaxLagTime;}
                lagTime = parmt_utils_getLagTime(data.data[iobs],
                                                 parms.defaultMaxLagTime,
                                                 &ldefault);
                nlags[iobs] = (int)
                              (lagTime/data.data[iobs].header.delta + half);
            }
            lwantLags = false;
            if (nlags[iobs] > 0){lwantLags = true;}
            npmax = MAX(npmax, data.data[iobs].npts);
            if (parms.lwantLags)
            {
                printf("%s: Number of lags for waveform %d is %d\n",
                       __func__, iobs+1, nlags[iobs]);
            }
        }
    }
    MPI_Bcast(nlags, data.nobs, MPI_INT,    master, globalComm);
    MPI_Bcast(&npmax,        1, MPI_INT,    master, globalComm);
    MPI_Bcast(&lwantLags,    1, MPI_C_BOOL, master, globalComm);
    // Check for failures
    if (npmax < 1 || mtloc.nmtAll < 1 || parms.objFnType < 1 ||
        parms.objFnType > 3)
    {
        if (npmax < 1){printf("%s: No data points\n", __func__);}
        if (mtloc.nmtAll < 1){printf("%s: Empty mts\n", __func__);}
        if (parms.objFnType < 1 || parms.objFnType > 3)
        {
            printf("%s: Invalid objfn function\n", __func__);
        }
        return -1;
    }
    if (lwantLags)
    {
        if (mtloc.myid == master)
        {
            lagsWork = memory_calloc32i(data.nlocs*mtloc.nmtAll);
        }
        else
        {
            lagsWork = memory_calloc32i(1);
        }
    }
    // Set space
    var = memory_calloc64f(npmax);
    if (linObsComm){MPI_Comm_size(obsComm, &nobsProcs);}
    if (linObsComm)
    {
        if (nobsProcs > 1)
        {
            phiAll = memory_calloc64f(data.nlocs*mtloc.nmtAll);
        }
        else
        {
            if (phi == NULL)
            {
                printf("%s: phi not set on master - segfault\n", __func__);
                return -1;
            }
            phiAll = phi;
            array_zeros64f_work(data.nlocs*mtloc.nmtAll, phiAll);
        }
        phiWork = memory_calloc64f(data.nlocs*mtloc.nmtAll);
        varAll = memory_calloc64f(data.nobs*npmax);
        defaultWeight = 1.0; ///(double) data.nobs;
        wts = memory_calloc64f(data.nobs); //array_set64f(data.nobs, defaultWeight, &ierr);
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            wts[iobs] = parmt_utils_getWeight(data.data[iobs],
                                              defaultWeight,
                                              &ldefault);
        }
        // Normalize so weighted sum is unity
        wtNorm = 1.0/array_sum64f(data.nobs, wts, &ierr);
        // And rescale the weights
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            wts[iobs] = wts[iobs]*wtNorm;
        }
    }
    else
    {
        phiWork = memory_calloc64f(1);
        phiAll = memory_calloc64f(1);
    }
    // Set up the MPI RMA window
    value = 0;
    value0 = 0;
    if (linObsComm)
    {
        MPE_setKeyval(MPI_KEYVAL_INVALID);
        MPE_counterCreate(obsComm, 1, &counterWin);
    }
    MPI_Barrier(globalComm);
    // Work-steal the observations to accomodate different sizes/lags
    while (true)
    {
        // Update the counter
        if (linObsComm){MPE_counterNextval(counterWin, 0, &value);}
        // Send the counter to the locations commuincator
        if (linLocComm){MPI_Bcast(&value, 1, MPI_INT, master, locComm);}
        // Trickle the value down to the moment tensor communicator
        MPI_Bcast(&value, 1, MPI_INT, master, mtloc.comm);
        if (value > data.nobs){break;}  // We are done
        if (value == value0){continue;} // Already did this
        value0 = value;   // Update the old value
        iobs = value - 1; // Map to C based numbering to indicate observation
        if (linObsComm)
        {
            printf("%s: Processing observation %d with %d points\n", __func__,
                   iobs+1, data.data[iobs].npts);
        }
        // Figure out the lags
        lwantLags = false;
        if (nlags[iobs] > 0){lwantLags = true;}
        if (parms.objFnType == 1)
        {
            ierr = parmt_locSearchL164f(locComm,
                                        iobs, parms.blockSize,
                                        nlags[iobs], lwantLags,
                                        parms.lrescale,
                                        mtloc, &data,
                                        CeInv, phiWork, var, lagsWork);
            if (ierr != 0)
            {
                printf("%s: L1 search error for observation %d\n",
                       __func__, iobs + 1);
                return -1;
            }
            // Put the variance into the global buffer
            if (linObsComm)
            {
                array_copy64f_work(npmax, var, &varAll[iobs*npmax]);
            }
        }
        else if (parms.objFnType == 2)
        {
            ierr = parmt_locSearchXC64f(locComm, iobs, parms.blockSize,
                                        nlags[iobs], lwantLags,
                                        mtloc, &data, phiWork, lagsWork);
            if (ierr != 0)
            {
                printf("%s: XC search error for observation %d\n",
                       __func__, iobs + 1);
                return -1;
            }
        }
        else if (parms.objFnType == 3)
        {
            printf("%s: need to write this function\n", __func__);
            return -1;
        }
        // Scale by number of observations and stack
        if (linObsComm)
        {
            //printf("%d %f\n", iobs, wts[iobs]);
            cblas_daxpy(data.nlocs*mtloc.nmtAll, wts[iobs], //1.0/(double) data.nobs,
                        phiWork, 1, phiAll, 1);
        }
    }
    // Reduce the objective function
    if (linObsComm)
    {
        if (nobsProcs > 1)
        {
            MPI_Reduce(phiAll, phi, data.nlocs*mtloc.nmtAll, MPI_DOUBLE,
                       MPI_SUM, master, obsComm);
        }
        MPI_Win_free(&counterWin);
    }
    MPI_Barrier(mtloc.comm);
    // Release resources
    if (linObsComm)
    {
        if (nobsProcs > 1){memory_free64f(&phiAll);}
        // Reduce the variance
        varOut = NULL;
        if (myid == master){varOut = memory_calloc64f(npmax*data.nobs);}
        MPI_Reduce(varAll, varOut, npmax*data.nobs, MPI_DOUBLE,
                   MPI_SUM, master, obsComm);
        // Have the master write it
        if (myid == master)
        {
            if (parms.programID == PARMT_ID)
            {
                archiveName = parms.parmtArchive;
            }
            else if (parms.programID == POLARMT_ID)
            {
                archiveName = parms.polarmtArchive;
            }
            for (iobs=0; iobs<data.nobs; iobs++)
            {
                parmt_io_writeVarianceForWaveform64f(
                   //parms.resultsDir, parms.projnm, parms.resultsFileSuffix,
                   archiveName,
                   iobs,
                   data.data[iobs].npts, &varOut[iobs*npmax]);
            }
            archiveName = NULL;
        }
        memory_free64f(&varAll);
        memory_free64f(&varOut);
    }
    else
    {
        memory_free64f(&phiAll);
    }
    memory_free64f(&phiWork);
    memory_free32i(&nlags);
    memory_free32i(&lagsWork);
    memory_free64f(&CeInv);
    memory_free64f(&wts);
    phiAll = NULL;
    MPI_Barrier(globalComm);
    return 0;
}


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#ifdef MAC_INTEGRATION
#include <gtkosxapplication.h>
#endif
#include <gtk/gtk.h>
#include <gdk/gdk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

#define DEFINE_GTK_INIT(_class, _type, _name) \
_class * \
win_init_ ##_name (GtkBuilder *b, const gchar *name) \
{ \
	g_assert(gtk_builder_get_object(b, name)); \
	g_assert(GTK_IS_##_type (gtk_builder_get_object(b, name))); \
	return(GTK_ ##_type (gtk_builder_get_object(b, name))); \
}

DEFINE_GTK_INIT(GtkAdjustment, ADJUSTMENT, adjustment)
DEFINE_GTK_INIT(GtkStatusbar, STATUSBAR, status)
DEFINE_GTK_INIT(GtkDrawingArea, DRAWING_AREA, draw)
DEFINE_GTK_INIT(GtkMenuBar, MENU_BAR, menubar)
DEFINE_GTK_INIT(GtkCheckMenuItem, CHECK_MENU_ITEM, menucheck)
DEFINE_GTK_INIT(GtkMenuItem, MENU_ITEM, menuitem)
DEFINE_GTK_INIT(GtkWindow, WINDOW, window)
DEFINE_GTK_INIT(GtkLabel, LABEL, label)
DEFINE_GTK_INIT(GtkBox, BOX, box)
DEFINE_GTK_INIT(GtkButton, BUTTON, button)
DEFINE_GTK_INIT(GtkRadioButton, RADIO_BUTTON, radio)
DEFINE_GTK_INIT(GtkNotebook, NOTEBOOK, notebook)
DEFINE_GTK_INIT(GtkToggleButton, TOGGLE_BUTTON, toggle)
DEFINE_GTK_INIT(GtkEntry, ENTRY, entry)
DEFINE_GTK_INIT(GtkFileChooser, FILE_CHOOSER, filechoose)

GtkBuilder *
builder_get(const gchar *name)
{
	gchar		*file = NULL;
#ifdef MAC_INTEGRATION
	gchar		*dir;
#endif
	GtkBuilder	*builder;

	/*
	 * Look up our `glade' file as follows.
	 * If we're in MAC_INTEGRATION, then intuit whether we're in a
	 * bundle and, if so, look up within the bundle.
	 * If we're not in a bundle, look in DATADIR.
	 */
#ifdef MAC_INTEGRATION
	g_object_new(GTKOSX_TYPE_APPLICATION, NULL);
	if (NULL != (dir = gtkosx_application_get_bundle_id())) {
		g_free(dir);
		dir = gtkosx_application_get_resource_path();
		file = g_strdup_printf
			("%s" G_DIR_SEPARATOR_S
			 "share" G_DIR_SEPARATOR_S
			 "bmigrate" G_DIR_SEPARATOR_S
			 "%s", dir, name);
		g_free(dir);
	}
#endif
	if (NULL == file)
		file = g_strdup_printf(DATADIR
			G_DIR_SEPARATOR_S "%s", name);

	builder = gtk_builder_new();
	assert(NULL != builder);

	/*
	 * This should be gtk_builder_new_from_file(), but we don't
	 * support that on older versions of GTK, so do it like this.
	 */
	if ( ! gtk_builder_add_from_file(builder, file, NULL)) {
		g_free(file);
		g_object_unref(G_OBJECT(builder));
		return(NULL);
	}

	g_free(file);
	return(builder);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate