Actually, ignore the inverse_cdf function I have provided. It should generate a vlaue for kappa and it needs adjusting for values of thetahat other than zero.

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);
end

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]);
end

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.