Great submission. It would be nice to have cdf and inversion cdf for the vmpdf functions. Here's what I wrote for my needs
function p = circ_vmcdf(alpha, thetahat, kappa)
%integrates the pdf from an angle of -pi to an angle alpha
F = @(x)circ_vmpdf(x, thetahat, kappa);
p = quad(F,-pi(),alpha);
function theta = circ_vminv(p, thetahat, kappa)
%computes the inverse of the abovecirc_vmcdf.
fun =@(alpha)(circ_vmcdf(alpha, thetahat, kappa)-p);
theta = fzero(fun,[-pi pi]);
There appears to be a mistake in this code, sorry. Its quite clear if the following data is entered at the prompt:
[0 0.1 1 1.9]
[0.2 0.2 0.8 0.8]
I suggest its related to the max/min lines. I suggest that you evaluate which t are valid for entering the view port and then eveluate which is the max entering t and min entering t from each set of valid t for each. Also parallel lines should be checked too or risk div0.