  function int = risint (x0, int0)
% function int = risint (x0, int0)
% risint computes integral needed to find modes of graded index waveguide,
% Equation (8.125) in the book by N. N. Rao.
% Input:
%  x0   position for the end of the integration interval
%  int0 target value of integral (an integer + 1/2)
% Output: 
%  int  integral minus target
%
% To find the end x0 (between 0 and 1) of the integration interval 
% that yields  int = 0 for int0 = 1/2, write
% [x0, int]  =  fzero ('risint', [0. 1.], [], 1/2)

% A gaussian distribution with an rms halfwidth of 1/2 wavelength,
% approaching 1 far from the waveguide, 
% and equal to 2 at the center of the waveguide.
% risqd =  1. + normpdf(x, 0, 0.5)./normpdf(0, 0, 0.5);
% Then compute integral, from x=0 to x=x0, and subtract target value
% Do not integrate all the way to x=x0, because of singularity there.
int   =  4.*quad ('risqrt', 0, 0.9999.*x0, [], [], x0) - int0;
return
