ParticleMath.cs

static Matrix4F identityF = Matrix4.Identity.ToFloat();

// This is a combination of a scale, rotate and translate matrix. Rather than calculating three different matrices and multiplying
// them together, we can save float multiplications and additions by using this function
public static Matrix4F ParticleMatrix(float s, float rX, float rZ, Vector3 d)
{
    var sX = GetArrayedSinCheap(rX);
    var sZ = GetArrayedSinCheap(rZ);
    var cX = GetArrayedCosCheap(rX);
    var cZ = GetArrayedCosCheap(rZ);

    Matrix4F m = identityF;

    m.M11 = s * cZ;
    m.M12 = s * sZ;

    m.M21 = -s * cX * sZ;
    m.M22 = s * cX * cZ;
    m.M23 = s * sX;

    m.M31 = s * sX * sZ;
    m.M32 = -s * sX * cZ;
    m.M33 = s * cX;

    m.M41 = (float)d.X;
    m.M42 = (float)d.Y;
    m.M43 = (float)d.Z;

    return m;
}

private const float epsilon = 0.001f;
private static float[] arrayedSinCheap;
private const float halfPI = (float)Math.PI / 2;

// Faster than the Math.Sin() function, uses a direct lookup in a float array of pre-calculated Math.Sin() values
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static float GetArrayedSinCheap(float t)
{
    if (t < 0)
        return -arrayedSinCheap[(int)(-t % (Math.PI * 2) / epsilon)];

    return arrayedSinCheap[(int)(t % (Math.PI * 2) / epsilon)];
}

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static float GetArrayedCosCheap(float t)
{
    t += halfPI;

    if (t < 0)
        return -arrayedSinCheap[(int)(-t % (Math.PI * 2) / epsilon)];

    return arrayedSinCheap[(int)(t % (Math.PI * 2) / epsilon)];
}

// Initialise the arrayedSinCheap array - only happens once
public static void InitialiseTrigonometry()
{
    int elements = (int)(Math.PI * 2 / epsilon) + 1;
    arrayedSinCheap = new float[elements];

    int j = 0;
    for (double angleRadians = 0; angleRadians <= Math.PI * 2;
        angleRadians += epsilon)
    {
        arrayedSinCheap[j] = (float)Math.Sin(angleRadians);
        j++;
    }
}