// Computation of the potential flow around circle. 
// The solution is seeked in the form psi0 + gamma psi1 and gamma  
//janvier 2002
wait=0;

border a(t=0,2*pi) { x=5*cos(t);  y=5*sin(t) };	// approximates infinity
 
border b(t=0,2) { 				// the circle airfoil
	  x = cos(t*pi); 
	  y = sin(-t*pi); } 
   
mesh th = buildmesh(a(30)+b(70));

solve (psi0) with A(0){							// compute the potential with Joukowski
	pde(psi0) -laplace(psi0) = 0;
	on(a)psi0=y-0.*x; 							// 0 percent lift
	on(b)psi0= 0;
}; 			
//plot(psi0);	

solve(psi1) with A(1){							// prepare Joukowski correction
	pde(psi1) -laplace(psi1)= 0; 
	on(a) psi1 = 0; 
	on(b) psi1 = 1;
};
// plot(psi1);

gamma := -8;		//  circulation
 

femp1 psi = gamma*psi1+psi0; 
plot(psi);
cp = -1.0*(dx(psi)^2 + dy(psi)^2); 
u=dy(psi);
v=-dx(psi);
save("uu.txt",u);
save("vv.txt",v);
save("pp.txt",cp);
savemesh("toto.mesh");
plot(cp);