reflection matrix for any angle illustration
Let's say very small bullet travel with vel_x and vel_y. It hits a big ball surface and bounce back with a new_vel_x and new_vel_y; 
We only know the where the bullet hit the surface (bullet_x,bullet_y), the center of the ball (ball_x, ball_y), and the velocity of the bullet vel_x and vel_y. how can I get the reflected new velocity: new_vel_x and new_vel_y.
PS, I know the reflection matrix for reflection about y-axis, x-axis, and y=x, but I don't know the reflection matrix at any angle. 
 A: Reflection across a plane is just twice subtracting the projection to the surface normal. Let's say $\vec{n}$ is the normalised (unit) normal to the surface, which you can compute from the bullet contact and the ball center. What you need:
$$\vec{v}\mapsto \vec{v}-2\vec{n}(\vec{n}\cdot \vec{v})=\underbrace{(I-2\vec{n}\otimes\vec{n})}_{A}\vec{v}$$
Here, the notation simply means
$$A=\begin{bmatrix} 1-2n_x n_x & -2n_x n_y\\ -2n_x n_y & 1-2 n_y n_y\end{bmatrix}$$
So just get the unit normal vector and you have the matrix. That is, if you need the matrix at all - it's actually easier to just use the first equation I wrote down. You just need the dot product and then subtract the double projection, basically, v-=2*n*dot(v,n) in pseudo-code.
A: If reflection line is y=mx, then the reflection matrix is:
| a11  a12 |
| a21  a22 |

where
a11 = (1-m*m)/(m*m+1)
a12 = (2*m)/(m*m+1)
a21 = (2*m)/(m*m+1)
a22 = (-1+m*m)/(m*m+1)

the bullet speed is v (vel_x, vel_y), then new velocity (new_x, new_y) is:
new_x = a11*v_x + a12*v_y
new_y = a21*v_x + a22*v_y

This is a test on python:
# assume that the bullet speed is v (vel_x, vel_y)
vel_x = 1
vel_y = 4

# assume that where the bullet hits the surface is P (x,y)
x = 15
y = 15

# assume that the ball center is located at B (ball_x, ball_y)
ball_x = 30
ball_y = 15


if(ball_y-y != 0):

    # if the line for reflection is y=mx, the value m is:
    m = -1*(ball_x-x)/(ball_y-y)

    # this is the reflection matrix based on the value m
    a11 = (1-m*m)/(m*m+1)
    a12 = (2*m)/(m*m+1)
    a21 = (2*m)/(m*m+1)
    a22 = (-1+m*m)/(m*m+1)

    new_x = a11*v_x + a12*v_y
    new_y = a21*v_x + a22*v_y

else:
    new_x = v_x
    new_y = -v_y

