constint max_particles = 1000;
struct element{
float x[max_particles]; //Transverse horizontal position
float p_x[max_particles]; //transverse horizontal momentum
float y[max_particles]; //Transverse vertical position
float p_y[max_particles]; //Transverse vertical momentum
float t[max_particles]; //Time of flight relative to ideal reference particle
float p_t[max_particles]; // (Delta)E/(p_s*c) i.e. Energy relative to reference particle divided by the product of the nominal momentum of an on-energy particle and the velocity of light
} particle; //The struct object for the elements of particles is called particle
//Define the transport matrix R
float R[6][6] = {0}; //initialise the transport matrix R to all zero values
float Beta = 0.8; //Assign a value to relativistic Beta (Obviously this will be determined by program at a later stage)
float Gamma = 1.666; //As above for relativistic Gamma
float L = 100.0; //Length of the drift space
//To start with, we need to change the value of R to match a purpose, as at the moment it is set to all zeros
//For now, we will make it equal to the drift matrix
for(int j = 0; j < 5; ++j){ //Remember that the first matrix element begins at zero so the last is the 5th element
R[j][j] = 1.0; //Sets the diagonals to 1
}
R[1][0] = R[3][2] = L;
R[5][4] = L/(Beta*Beta*Gamma*Gamma);
float X[max_particles] = {0};
float P_X[max_particles] = {0};
float Y[max_particles] = {0};
float P_Y[max_particles] = {0};
float T[max_particles] = {0};
float P_T[max_particles] = {0};
//Now we need to create variable arrays to store the results before returning them to the host
for(int p = 0; p < N; ++p){ //For all particles
X[p] = R[0][0]*x[p] + R[1][0]*p_x[p] + R[2][0]*y[p] + R[3][0]*p_y[p] + R[4][0]*t[p] + R[5][0]*p_t;
P_X[p] = R[0][1]*x[p] + R[1][1]*p_x[p] + R[2][1]*y[p] + R[3][1]*p_y[p] + R[4][1]*t[p] + R[5][1]*p_t;
Y[p] = R[0][2]*x[p] + R[1][2]*p_x[p] + R[2][2]*y[p] + R[3][2]*p_y[p] + R[4][2]*t[p] + R[5][2]*p_t;
P_Y[p] = R[0][3]*x[p] + R[1][3]*p_x[p] + R[2][3]*y[p] + R[3][3]*p_y[p] + R[4][3]*t[p] + R[5][3]*p_t;
T[p] = R[0][4]*x[p] + R[1][4]*p_x[p] + R[2][4]*y[p] + R[3][4]*p_y[p] + R[4][4]*t[p] + R[5][4]*p_t;
P_T[p] = R[0][5]*x[p] + R[1][5]*p_x[p] + R[2][5]*y[p] + R[3][5]*p_y[p] + R[4][5]*t[p] + R[5][5]*p_t;
}
the variables x, p_x, y, p_y, t, p_t have all their values read into them by a file, and I am confident that that bit works correctly.
N is the number of particles as counted when inputting the data from file.
Hope this helps, and thanks again