function fm=oz(x)
%OZ FM=OZ(X) is the nonlinear residual for the 
%   Ornstein-Zernike equations.
%
% Evaluate the nonlinearity as a sum of the substitution
% and compact parts.
%
fm=ozsub(x)+ozinteg(x);
%
% This is the substitution part. 
%
function km=ozsub(x)
global L U rho
n=length(x); n2=n/2; h=x(1:n2); c=x(n2+1:n);
nl=U.*exp(h-c);
mup=h+1- nl;
mdown=c+1-nl;
km=[mup',mdown']';
%
% This is the compact part of the nonlinearity. The upper
% component is zero, the lower is rho*L*c*h.
%
% rho and L are scalars for this problem
%
function km=ozinteg(x)
global L U rho
L3=L*L*L*rho;
n=length(x); n2=n/2; h=x(1:n2); c=x(n2+1:n);
kd=L3*convk(c,h);
kd(1)=2*kd(2)-kd(3);
kd(n2)=2*kd(n2-1)-kd(n2-2);
kup=zeros(n2,1);
km=[kup',kd']';
%
% Discrete convolution with the Hankel transform.
%
function hc=convk(h,c)
Lx=1.0;
n=length(h);
dr=Lx/n;
th=hank2(h); tc=hank2(c); fhc=th.*tc; hc=ihank2(fhc);
%
% Hankel transform using the fast sine transform.
%
function hf=hank2(f)
nf=length(f); n=nf-1; m=nf-2;
h=1/n; beta=2*(m+1)*(h^3); ff=beta;
p=1:m; p=p';
hft=ff*lsint(p.*f(2:n))./p;
hf=[0,hft',0]';
%
% Inverse Hankel transform using the fast sine transform.
%
function ihf=ihank2(f)
nf=length(f); n=nf-1; m=nf-2;
h=1/n; beta=2*(m+1)*(h^3); ff=2/(n*beta);
p=1:m; p=p'; 
hft=ff*lsint(p.*f(2:n))./p;
p1=1:n; p1=p1'; p1=p1.*p1; p1(n)=p1(n)/2;
ihf=[0,hft',0]';
% LSINT
% Fast sine transform with MATLAB's FFT.
%
function lf=lsint(f)
n=length(f);
ft=-fft([0,f']',2*n+2);
lf=imag(ft(2:n+1));