// cet exemple est dans la doc de Freefem++ 
// 
// it was in teh very first version of MacGFem
//
// Computation of the potential flow around a NACA0012 airfoil.
// The method of decomposition is used to apply the Joukowski condition
// The solution is seeked in the form psi0 + beta psi1 and beta is
// adjusted so that the pressure is continuous at the trailing edge
// janvier 2003
// avril 2008 version ++
// mars 2011  version thanks to Geordie MacBain
//
// ./FreeFem++-CoCoa3 nacaPYL.edp 
// fichier genere  cff.data P.gp

verbosity = 0;
wait=0;
real Wing=99,beta,alpha;
int  np=64*2;      
real pis180=pi/180.; // angle conversion
real eps=0.001;      // for the derivative 
real R0=50;         // size of teh domain 


// def NACA0012
func real extra(real x){
    return 0.17735*sqrt(x)-0.075597*x 
       - 0.212836*(x^2)+0.17363*(x^3)-0.06254*(x^4);}
       
func real intra(real x){
	return   -(0.17735*sqrt(x)-0.075597*x 
        -0.212836*(x^2)+0.17363*(x^3)-0.06254*(x^4));}
       
// external an zoomed domains
border C(t=0,2*pi) { x=R0*cos(t)+0.5; y=R0*sin(t);} 
border c(t=0,2*pi) { x=1* cos(t)+0.5; y= 1*sin(t); }       

border Splus(t=0,1){ x = t; y = extra(t); label=Wing;} 
border Sminus(t=1,0){ x =t; y = intra(t); label=Wing;} 
   
// mesh generation   
mesh Th= buildmesh(C(np/2)+Splus(np)+Sminus(np)); 
mesh Zoom = buildmesh(c(np)+Splus(np)+Sminus(np));
plot (Th);

// FEM spaces
fespace Vh(Th,P2); 
Vh psi,psi0,psi1,w,phi0; 
Vh u,v,cp;
fespace ZVh(Zoom,P2);

// zoom 
plot(Zoom,wait=0);


alpha=20;
solve pot(phi0,w)=   
   int2d(Th)(dx(phi0)*dx(w)+dy(phi0)*dy(w)) 
    + on(C,phi0= x*cos(alpha*pis180) +y*sin(alpha*pis180));
    
 ZVh Zphi0=phi0; 
  
 plot(Zphi0,cmm=" sans aucune circulation  a="+alpha+" degres,"); 


{
// values are strored in cff.data
	ofstream ff("cff.data");

// loop in angle of attack
alpha=20;
for(alpha=5;alpha<120;alpha=alpha+5)
{	
// compute the potential without Joukowski
solve potential(psi0,w)=   
   int2d(Th)(dx(psi0)*dx(w)+dy(psi0)*dy(w)) 
    + on(C,psi0= y*cos(alpha*pis180) -sin(alpha*pis180)*x)    // alpha lift
    + on(Wing,psi0=0); 
    
 ZVh Zpsi0=psi0; 
 real circZ0;
  circZ0 = int1d (Zoom, c) (((dx(Zpsi0)*N.x + dy(Zpsi0)*N.y)));

 plot(Zpsi0,cmm=" sans circulation  a="+alpha+" degres,"+ " "+ circZ0); 

// prepare Joukowski correction
solve potentialG(psi1,w)=   
   int2d(Th)(dx(psi1)*dx(w)+dy(psi1)*dy(w)) 
    + on(C,psi1=0)    // alpha lift
    + on(Wing,psi1=1); 
    
 ZVh Zpsi1=psi1;
 real circZ1;
  circZ1 = int1d (Zoom, c) (((dx(Zpsi1)*N.x + dy(Zpsi1)*N.y)));
   
 plot(Zpsi1,cmm=" que circulation  a="+alpha+" degres"+ circZ1);

// continuity of pressure/ velocity at trailing edge
beta = psi0(1-eps,eps)+psi0(1-eps,-eps);		
beta = -beta / (psi1(1-eps,eps)+ psi1(1-eps,-eps)-2);

cout << "beta=" <<beta<<" " <<  psi0(1-eps,eps) << " " << psi1(1-eps,eps) << "\n" ;
 
// reconstructing the total solution 
psi = beta*psi1+psi0;  
ZVh  Zpsi=psi; 

// computing the circulation by integration of the tangential velocity
  real circ,circZ,circW;
  circZ =  int1d (Zoom, c) (((dx(Zpsi)*N.x + dy(Zpsi)*N.y)));
  circ  =  int1d (Th  , C) (((dx(psi )*N.x + dy(psi )*N.y)));
  circW =  int1d (Th  , C) (((dx(psi )*N.x + dy(psi )*N.y)));
 
plot(Zpsi,bw=true,nbiso=50,
    cmm=" Solution totale avec circulation \n "+ "beta=" + beta + "\n c= "
     + circ + ", czoom= " + circZ + " c wing =" + circW); 
 
// vitesse
ZVh Zu=u,Zv=v,Zcp;
Zu=dy(Zpsi);
Zv=-dx(Zpsi);
// coeff de pression
cp = -1.0*(dx(psi)^2 + dy(psi)^2); 
Zcp=cp;
plot([Zu,Zv],Zpsi,Zcp,cmm=" psi u et v ",coef=.1);
//plot(cp,cmm=" Cp  ",wait=wwait); 


//macro ut (psi) (-(dx (psi) * N.x + dy (psi) * N.y)) // tangential velocity
//macro circulation (psi, patch) int1d (Th, patch) (-ut (psi))//EOM

 Vh psii;
 psii = psi-(y*cos(alpha*pis180) -sin(alpha*pis180)*x); 
//plot(psii,bw=true,nbiso=50,
//    cmm=" Solution totale avec circulation \n "+ "beta=" + beta + ", c= " + circ + ", czoom= " + circZ); 
 
 


// save values on the wing 
int n=100;
{ ofstream gnu("P.gp");  
for (int i=0;i<=n;i++) 
{ x = (i*1.)/n; 
 real ywe=extra(x); 
 real ywi=intra(x);  
gnu<< x <<" "<<  ywe << " " << ywi << " " << cp(x,ywe)-cp(1,0)  << " "<< cp(x,ywi)-cp(1,0)  <<" "  << endl; }
 }

n=200; 
{ ofstream gnu("logP.gp");  
for (int i=0;i<=n;i++) 
{ x = (i*R0)/n; 
  y = 0;  
gnu<< x <<" "<<  psii(x,y)  << endl; }
 }
 

 cout <<"alpha   -beta -----------------------\n" ;
 cout <<  alpha <<" " <<  -beta<< " " << circ << " " << circZ << "\n" ;
// save in the file
 ff <<  alpha <<" " <<  -beta<< " " << circ << " " << circZ << "\n" ;	
}}

cout <<"Fin -----------------------\n" ;