I should estimate the effect of B fields on particle deflection using \$F=V \times B \times Force\$.
I have a 2D mesh with \$n_x\$ and \$n_y\$ which show the number of points in each direction. \$B\$ is the magnetic field which is a 2D mesh arrays (with dimension \$n_x \times n_y= 1600 \times 900 = 1440000\$) and \$V\$ is the velocity of particles and is a 1D array (with dimension 431605462). I also have two 1D arrays (with the same dimension as the \$V\$ array) which show the x and y coordinates of each particles. I called them
Since \$V\$ is a 1D array and \$B\$ is a 2D array I can not directly multiply \$V\$ with \$B\$. To solve this issue I created a 2D zero array for the \$V\$ and using 2 for-loops I find particles which stand in each cell of the mesh. I consider the average velocity of all particles as a velocity component of each cell. Because of the large dimensionality of the mesh, if I keep the original dimension, the code needs several days to complete the \$V\$ array since it needs to be repeated 1440000 times. As a solution, I compressed \$V\$ with an optional variable called
divider. If I use
divider = 10 I guess the code needs 2 days to be finished and if I use
divider = 100 the code will be finished in half an hour but now the precision of the results reduces significantly.
You can see the code in the following.
nx,ny= 1600, 900 xmin = -30.0e-6 xmax = 50.0e-6 ymin = -45.0e-6 ymax = 45.0e-6 #### create a 2D velocity array divider = 100 m_p= 1.6726219e-27 # (kg unit) proton mass nx_new = nx/divider ny_new = ny/divider Vx = np.zeros((ny_new, nx_new)) Vy = np.zeros((ny_new, nx_new)) res_x = (xmax - xmin)/nx_new # resolution in x direction res_y = (ymax - ymin)/ny_new # resolution in y direction # Find corresponding indices for i in range( 0,ny_new ): for j in range( 0,nx_new): xx = i*res_x yy = j*res_y ##### all particles with different Px confined in a specific cell Px_pro_lim = Px_pro[(grid_pro_x >= xx ) & (grid_pro_x < xx+res_x ) & (grid_pro_y >= yy ) & (grid_pro_y < yy+res_y ) ] Vx[i,j] = np.mean(Px_pro_lim)/m_p #print 'Vx[i,j]= ' ,Vx[i,j]