// the following variables must be set on entry:
// NACA4_th     thickness in percent of chord
// NACA4_ch     aerofoil chord
// NACA4_le_x,y,z   leading edge coordinates

// constants from NASA TM4741
NACA4_a0 =  0.2969 ;
NACA4_a1 = -0.1260 ;
NACA4_a2 = -0.3537 ; // 0.3516 on http://www.aerospaceweb.org/question/airfoils/q0100.shtml but then the sum of the corfficients is not 1.0 !
NACA4_a3 =  0.2843 ;
NACA4_a4 = -0.1015 ;

NACA4_npts = 500;
NACA4_ch = 1.0;
NACA4_le_x = 0.0 ;
NACA4_le_y = 0.0 ;
NACA4_le_z = 0.0 ;

NACA4_Points[] = {} ;
NACA4_Splines[] = {} ;
xu[] = {};
yu[] = {};

// Definition du Naca####
camber    = 6/100;    // First Digit : maximum camber in percentage of the chord (m)
position  = 4/10; // Second Digit : position of the maximum camber in tenths of the chord (p)
Thickness = 12/100;   // Two last Digits : thickness of the airfoil in percentage of the chord (t)

// We here suppose that NACA4_npts is even...

// First half of the NACA starting from the trailing edge
For i In {0:NACA4_npts/2}
    x  = 0.5*(1+Cos(i/NACA4_npts*Pi));
    x  = x*((2-4*position)*x+4*position-1); // x[i=NACA4_npts/2] is now situated at the maximum camber
    yc = camber/(1-position)^2*(1-2*position+2*position*x-x^2);
    yt =5.0*Thickness*(NACA4_a0*Sqrt(x) + 
       	  x*(NACA4_a1 + x*(NACA4_a2 + x*(NACA4_a3 + x*NACA4_a4)))) ;
    theta = Atan(2*camber/(1-position)^2*(position-x));
    xu[i] = x-yt*Sin(theta);
    yu[i] = yc+yt*Cos(theta);
    pt_oppose = 2*NACA4_npts-i;
    xu[pt_oppose] = x+yt*Sin(theta);
    yu[pt_oppose] = yc-yt*Cos(theta);
EndFor

// Second Half
For i In {NACA4_npts/2+1:NACA4_npts}
    x  = 0.5*(1+Cos(i/NACA4_npts*Pi));
    x  = x*((2-4*position)*x+4*position-1);
    inter = x/position;
    yc = camber*inter*(2-inter);
    yt = 5.0*Thickness*(NACA4_a0*Sqrt(x) + 
       	  x*(NACA4_a1 + x*(NACA4_a2 + x*(NACA4_a3 + x*NACA4_a4)))) ;
    theta = Atan(2*camber/position*(1-inter));
    xu[i] = NACA4_ch*(x-yt*Sin(theta))+NACA4_le_x;
    yu[i] = NACA4_ch*(yc+yt*Cos(theta))+NACA4_le_y;
    pt_oppose = 2*NACA4_npts-i;
    xu[pt_oppose] = NACA4_ch*(x+yt*Sin(theta))+NACA4_le_x;
    yu[pt_oppose] = NACA4_ch*(yc-yt*Cos(theta))+NACA4_le_y;
EndFor

// Computation of the Characteristic Length such that the defined 
// points and only them are on the boundary...
distance[0]= Sqrt((xu[1]-xu[0])^2 +(yu[1]-yu[0])^2);
For i In {1:NACA4_npts*2-1}
    distance[i]= Sqrt((xu[i+1]-xu[i])^2 +(yu[i+1]-yu[i])^2);
    length[i] = 0.5*(distance[i]+distance[i-1]);
EndFor
length[0] = 0.5*(distance[0]+distance[NACA4_npts*2-1]);

// Creation des Points
For i In {0:NACA4_npts*2-1}
    Point(i+1) = {xu[i],yu[i],NACA4_le_z,length[i]};
EndFor

// Creation des lignes
j=NACA4_npts*2+1;
lNaca[0]=j;
For i In {0:NACA4_npts*2-2}
    Line(j) = {i+1,i+2};
    j++;
    lNaca[i+1]=j;	
EndFor
Line(j) = {NACA4_npts*2,1};


// Line Loop
ll = newll;
Line Loop(ll)={lNaca[]};


// FarField
i=0; ext = 1.0;
pCercle[i] = newp;Point(pCercle[i]) = {11,0,0,ext};
i++;pCercle[i] = newp;Point(pCercle[i]) = {1,0,0,ext};
i++;pCercle[i] = newp;Point(pCercle[i]) = {1,10,0,ext};
quartCercle = newl;
Circle(quartCercle) = {pCercle[]};
c1[] = Rotate{{0,0,1},{1,0,0},Pi/2} {Duplicata{Line{quartCercle};}};
c2[] = Rotate{{0,0,1},{1,0,0},Pi} {Duplicata{Line{quartCercle,c1[0]};}};
ligneExt[] = {quartCercle,c1[0],c2[0],c2[1]};
cercleExt = newll; 
Line Loop(cercleExt) = {ligneExt[]};
Physical Line(2) = {ligneExt[]};


// Surface
surf = news; 
/*
Printf("%g",surf);
Printf("%g %g",cercleExt,ll);
For i In {0:# ligneExt[]-1}
    Printf("%g %g",i,ligneExt[i]);
EndFor	
For i In {0:# lNaca[]-1}
    Printf("%g %g %g",i,lNaca[i],length[i]);
EndFor	
*/
Plane Surface(surf) = {cercleExt,ll};
Physical Surface(0) = surf;


