Exception for Trowbridge-Reitz Sample, sample surface normal, cosTheta

\[\left(\left(\left(2.328306437 \cdot 10^{-10} \leq u0 \land u0 \leq 1\right) \land \left(2.328306437 \cdot 10^{-10} \leq u1 \land u1 \leq 0.5\right)\right) \land \left(0.0001 \leq alphax \land alphax \leq 1\right)\right) \land \left(0.0001 \leq alphay \land alphay \leq 1\right)\]

Math

\[\frac{1}{\sqrt{1 + \frac{\frac{1}{\frac{\cos \tan^{-1} \left(\frac{alphay}{alphax} \cdot \tan \left(\left(2 \cdot \pi\right) \cdot u1 + 0.5 \cdot \pi\right)\right) \cdot \cos \tan^{-1} \left(\frac{alphay}{alphax} \cdot \tan \left(\left(2 \cdot \pi\right) \cdot u1 + 0.5 \cdot \pi\right)\right)}{alphax \cdot alphax} + \frac{\sin \tan^{-1} \left(\frac{alphay}{alphax} \cdot \tan \left(\left(2 \cdot \pi\right) \cdot u1 + 0.5 \cdot \pi\right)\right) \cdot \sin \tan^{-1} \left(\frac{alphay}{alphax} \cdot \tan \left(\left(2 \cdot \pi\right) \cdot u1 + 0.5 \cdot \pi\right)\right)}{alphay \cdot alphay}} \cdot u0}{1 - u0}}} \]

FPCore

(FPCore (u0 u1 alphax alphay)
 :precision binary64
 (/
  1.0
  (sqrt
   (+
    1.0
    (/
     (*
      (/
       1.0
       (+
        (/
         (*
          (cos
           (atan (* (/ alphay alphax) (tan (+ (* (* 2.0 PI) u1) (* 0.5 PI))))))
          (cos
           (atan
            (* (/ alphay alphax) (tan (+ (* (* 2.0 PI) u1) (* 0.5 PI)))))))
         (* alphax alphax))
        (/
         (*
          (sin
           (atan (* (/ alphay alphax) (tan (+ (* (* 2.0 PI) u1) (* 0.5 PI))))))
          (sin
           (atan
            (* (/ alphay alphax) (tan (+ (* (* 2.0 PI) u1) (* 0.5 PI)))))))
         (* alphay alphay))))
      u0)
     (- 1.0 u0))))))

C

double code(double u0, double u1, double alphax, double alphay) {
	return 1.0 / sqrt((1.0 + (((1.0 / (((cos(atan(((alphay / alphax) * tan((((2.0 * ((double) M_PI)) * u1) + (0.5 * ((double) M_PI))))))) * cos(atan(((alphay / alphax) * tan((((2.0 * ((double) M_PI)) * u1) + (0.5 * ((double) M_PI)))))))) / (alphax * alphax)) + ((sin(atan(((alphay / alphax) * tan((((2.0 * ((double) M_PI)) * u1) + (0.5 * ((double) M_PI))))))) * sin(atan(((alphay / alphax) * tan((((2.0 * ((double) M_PI)) * u1) + (0.5 * ((double) M_PI)))))))) / (alphay * alphay)))) * u0) / (1.0 - u0))));
}

Java

public static double code(double u0, double u1, double alphax, double alphay) {
	return 1.0 / Math.sqrt((1.0 + (((1.0 / (((Math.cos(Math.atan(((alphay / alphax) * Math.tan((((2.0 * Math.PI) * u1) + (0.5 * Math.PI)))))) * Math.cos(Math.atan(((alphay / alphax) * Math.tan((((2.0 * Math.PI) * u1) + (0.5 * Math.PI))))))) / (alphax * alphax)) + ((Math.sin(Math.atan(((alphay / alphax) * Math.tan((((2.0 * Math.PI) * u1) + (0.5 * Math.PI)))))) * Math.sin(Math.atan(((alphay / alphax) * Math.tan((((2.0 * Math.PI) * u1) + (0.5 * Math.PI))))))) / (alphay * alphay)))) * u0) / (1.0 - u0))));
}

Python

def code(u0, u1, alphax, alphay):
	return 1.0 / math.sqrt((1.0 + (((1.0 / (((math.cos(math.atan(((alphay / alphax) * math.tan((((2.0 * math.pi) * u1) + (0.5 * math.pi)))))) * math.cos(math.atan(((alphay / alphax) * math.tan((((2.0 * math.pi) * u1) + (0.5 * math.pi))))))) / (alphax * alphax)) + ((math.sin(math.atan(((alphay / alphax) * math.tan((((2.0 * math.pi) * u1) + (0.5 * math.pi)))))) * math.sin(math.atan(((alphay / alphax) * math.tan((((2.0 * math.pi) * u1) + (0.5 * math.pi))))))) / (alphay * alphay)))) * u0) / (1.0 - u0))))

Julia

function code(u0, u1, alphax, alphay)
	return Float64(1.0 / sqrt(Float64(1.0 + Float64(Float64(Float64(1.0 / Float64(Float64(Float64(cos(atan(Float64(Float64(alphay / alphax) * tan(Float64(Float64(Float64(2.0 * pi) * u1) + Float64(0.5 * pi)))))) * cos(atan(Float64(Float64(alphay / alphax) * tan(Float64(Float64(Float64(2.0 * pi) * u1) + Float64(0.5 * pi))))))) / Float64(alphax * alphax)) + Float64(Float64(sin(atan(Float64(Float64(alphay / alphax) * tan(Float64(Float64(Float64(2.0 * pi) * u1) + Float64(0.5 * pi)))))) * sin(atan(Float64(Float64(alphay / alphax) * tan(Float64(Float64(Float64(2.0 * pi) * u1) + Float64(0.5 * pi))))))) / Float64(alphay * alphay)))) * u0) / Float64(1.0 - u0)))))
end

MATLAB

function tmp = code(u0, u1, alphax, alphay)
	tmp = 1.0 / sqrt((1.0 + (((1.0 / (((cos(atan(((alphay / alphax) * tan((((2.0 * pi) * u1) + (0.5 * pi)))))) * cos(atan(((alphay / alphax) * tan((((2.0 * pi) * u1) + (0.5 * pi))))))) / (alphax * alphax)) + ((sin(atan(((alphay / alphax) * tan((((2.0 * pi) * u1) + (0.5 * pi)))))) * sin(atan(((alphay / alphax) * tan((((2.0 * pi) * u1) + (0.5 * pi))))))) / (alphay * alphay)))) * u0) / (1.0 - u0))));
end

Wolfram

code[u0_, u1_, alphax_, alphay_] := N[(1.0 / N[Sqrt[N[(1.0 + N[(N[(N[(1.0 / N[(N[(N[(N[Cos[N[ArcTan[N[(N[(alphay / alphax), $MachinePrecision] * N[Tan[N[(N[(N[(2.0 * Pi), $MachinePrecision] * u1), $MachinePrecision] + N[(0.5 * Pi), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision] * N[Cos[N[ArcTan[N[(N[(alphay / alphax), $MachinePrecision] * N[Tan[N[(N[(N[(2.0 * Pi), $MachinePrecision] * u1), $MachinePrecision] + N[(0.5 * Pi), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(alphax * alphax), $MachinePrecision]), $MachinePrecision] + N[(N[(N[Sin[N[ArcTan[N[(N[(alphay / alphax), $MachinePrecision] * N[Tan[N[(N[(N[(2.0 * Pi), $MachinePrecision] * u1), $MachinePrecision] + N[(0.5 * Pi), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision] * N[Sin[N[ArcTan[N[(N[(alphay / alphax), $MachinePrecision] * N[Tan[N[(N[(N[(2.0 * Pi), $MachinePrecision] * u1), $MachinePrecision] + N[(0.5 * Pi), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(alphay * alphay), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * u0), $MachinePrecision] / N[(1.0 - u0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]