function [xp,yp,psi]=corner(n,mu)
%	calculate psi=r^m1 f_1(theta)+mu*r^m2 f_2(theta )
%	according to Weinbaum's solution
	m1=1.5445058;
	m2=1.9085318;
	a=3*m1*pi/4;
	ct=cos(a)/sin(a);
	a=3*m2*pi/4;
	tc=sin(a)/cos(a);
	h=0.5/n;
	m=2*n+1;
	for i=1:m,
		for j=1:m,
			x=(i-1-n)*h;
			y=(j-1-n)*h;
			r=(x*x+y*y)^0.5;
			t=atan2(y,x)-pi/4;
			f1=cos(m1*t)+ct*cos((m1-2)*t);
			f2=sin(m2*t)-tc*sin((m2-2)*t);
			psi(j,i)=f1*r^m1+mu*f2*r^m2;
			xp(j,i)=x;
			yp(j,i)=y;
			if x<0.
				if y<0. 
					psi(j,i)=0.;
				end
			end
		end
	end
levels=[-0.01 -0.005 -0.0000001 0.0000001 0.025 0.05 0.075 0.1 0.2 0.3];
contour(xp,yp,psi,levels);