function [f,omega,a,x,p1]=riley(cp,n)

% implementation of Riley's method to compute
% Sadoviskii vortex
%
%   [f,omega,a,x,p1]=riley(cp,n)
% 
%   cp = pressure coefficient
%   n  = number of points for solution
%
%   f     = boundary shape (so plot(x,f) will draw shape)
%   omega = value of vorticity
%   a     = vector of coefficents A_n in (7.65)
%   x     = x locations of points in f
%   p1    = vortex sheet strength
%
	n2=n;
	for i=1:n,
		a(1,1)=0.;
		da(1)=0.;
		b(i,1)=0.;
		e(i,1)=0.;
		x(i)=-cos((i-1)*pi/n2);
		for j=2:n,
			da(i)=0.;
			b(i,j)=cos(j*(i-1)*pi/(n2));
			u=sin((j-1)*(i-1)*pi/(n2))/(j-1);
			v=sin((j+1)*(i-1)*pi/(n2))/(j+1);
			e1(i,j)=(u-v)/2;
		end;
	end;
	x(1)=-1.;
	for j=2:2:n,
		a(j,1)=4*cp/j;
	end;
	for loop=1:100,
		s=e1*a;
		t=b*a;
		for i=1:n,
		c(i,1)=s(i)*s(i)/8;
			for j=2:n,
				c(i,j)=b(i,j)+a(1,1)*s(i)*e1(i,j)/4;
				r(i,1)=cp-t(i)-a(1,1)*s(i)*s(i)/8;
			end;
		end;
		da=inv(c)*r;
		a=a+da;
		a(1,1)
	end;
	f=e1*a;
	p1=b*a;
	omega=sqrt(a(1,1))


