function a = solvebrod(alpha,k)
%
% a = solvebrd(k) calculates an iterated solution for the
% first k coefficients A(1)->A(2*k-1) of Brodetsky's A(t) function
% when separation occurs at an angle alpha (degrees)
%
% A fixed number of iterations (20) is preset
%
	iterations=20;
	s=[0:pi/(2*(k-1)):pi/2];
	for i=1:k,
		a(i)=0;
	end
	a(1)=-1;
	alpha1=alpha*pi/180.;
	for j=1:k-1,
		for i=1:k,
			n=2*i-1;
			g(j,i)=cos(n*s(j))/cos(s(j));
		end
	end
	for i=1:k,
		n=2*i-1;
		g(k,i)=-(-1)^i/n;
	end
	for k1=1:iterations,
		for j=1:k-1,
			f(j)=1.;
			for i=1:k,
				n=2*i-1;
				f(j)=f(j)*exp(a(i)*cos(n*s(j))/n);
			end
			if j==1 
				f(j)=2*f(j);
			else
				f(j)=sin(s(j))*sqrt((1+cos(s(j)))/(1-cos(s(j))))*f(j);
			end
			for i=1:k,
				mat(j,i)=g(j,i)+(2*i-1)*f(j)*(-1)^i;
				
			end
			b(j)=0.;
		end
		for i=1:k,
			mat(k,i)=g(k,i);
		end
		b(k)=-alpha1;
		a=inv(mat)*b';
	end