// Aripa cu profil Naca0012
// 07.08.2002
// Mihai Victor Pricop

crad=Pi/180.;
hi=75.*crad;         // Unghiul de sageata al aripii - la sfert de coarda
c=1.008930411365;    // Coarda aripii la incastrare
ce=.2*c;             // Coarda aripii la extrmitate
lambda=3.;           // Alungirea aripii
la=.5*(c+ce)*lambda; // Anvergura aripii
ud=6*crad;	     // Unghiul diedru

np=50; 	             // numarul de puncte de pe extrados
np2=3;              // numarul de puncte incepand cu cel de la BA care sunt 
                     // incluse in interpolantul bordului de atac pe intrados
                     // si in continuare  pe intrados
pi=4.*Atan(1.);


                     // Sectiunea de incastrare
y=0;                 // y-ul planului de incastrare


eps=.7;
ax=1.-2.*eps; bx=-eps;
For i In {1:np}
  xx=(np-i)/(np-1)*pi;
  x[i]=c/2*(1+Cos(xx));
  z[i]=5*0.12*(0.2969*Sqrt(x[i]) -0.126*x[i] -0.3516*x[i]^2+0.2843*x[i]^3 
-0.1015*x[i]^4);
  z[i]=z[i]/c;
EndFor

lc1=0.04;
// Extradosul profilului Naca0012
//lc[1]=.03;
//Point(1)={x[1],y,z[1],lc[1]};
For i In {1:np}
  lc[i]=.004+.1*(x[i]*(c-x[i]))^1.5;
//  lc[i]=.5*((x[i]-x[i-1])^2+(z[i]-z[i-1])^2)^0.5;
//  lc[i]=lc1*0.4;
  Point(i)={x[i],y,z[i],lc[i]};
EndFor
//For i In {np2+1:np}
//  lc[i]=.004+.1*(x[i]*(c-x[i]))^1.5;
//  lc[i]=.5*((x[i]-x[i-1])^2+(z[i]-z[i-1])^2)^0.5;
//  lc[i]=lc1;
//  Point(i)={x[i],y,z[i],lc[i]};
//EndFor


// Intradosul profilului
For i In {2:np-1}
  Point(np+i-1)={x[i],y,-z[i],lc[i]};
//  Printf("%g, %g, %g",i,x[i],-z[i]);
EndFor

Printf("Lungimea caracteristica 1 = %g",lc[1]);
Printf("Coarda aripii la incastrare %g",x[np]-x[1]);

// Sectiunea de extremitate
ye=-la/2;

// Noua origine - punctul din bordul de atac al profilului de extremitate
x0e=la/2*Sin(hi)-.25*ce;
z0e=la/2*Sin(ud); // Avem si unghi diedru
Printf("Coordonatele ba extr. : %g,%g,%g",x0e,y,z0e);

raport=ce/c;
i1=2*np-2+1; // Extradosul extremitatii
For i In {1:np}
  xe[i]=x0e+x[i]*raport;
  ze[i]=z[i]*raport;
  //Printf("%g, %g",lc[1],lc[i]*raport);
  Point(i1)={xe[i],ye,z0e+ze[i],lc[i]*raport};
  i1=newp;
EndFor

For i In {2:np-1}
  Point(i1)={xe[i],ye,z0e-ze[i],lc[i]*raport};
  i1=newp;
EndFor


baincastrare=newreg;
Spline(baincastrare)={np2:1,np+1:np+np2-1}; // Bord atac incastrare

baextremitate=newreg;
na=2*np-2; // Bord atac extremitate
Spline(baextremitate)={na+np2:na+1,na+np+1:na+np+np2-1};

// Bordul de fuga
bf=newreg;
Line(bf)={np,na+np};

// Muchie bord atac extrados
mbae=newreg;
Line(mbae)={np2,na+np2};

// Muchie bord atac intrados
mbai=newreg;
Line(mbai)={np+np2-1,na+np+np2-1};
//Attractor Line{mbai}={0.1,6,0.1};

//ba=newreg; Line(ba)={1,na+1};

// Suprafata bord atac
loopba=newreg;
Line Loop(loopba)={baincastrare,mbai,-baextremitate,-mbae};
sba=newreg;
Ruled Surface(sba)={loopba};
Physical Surface(3)={sba}; // Perete solid

//Trasnfinite Surface{100}={baincastrare,mbai,-baextremitate,-mbae};


ei=newreg; // Extrados incastrare
Spline(ei)={np2:np};

ii=newreg; // Intrados incastrare
Spline(ii)={np+np2-1:2*np-2,np};

ee=newreg; // Extrados extremitate
Spline(ee)={na+np2:na+np};

ie=newreg; // Intrados extremitate
Spline(ie)={na+np+np2-1:na+2*np-2,na+np};

// Suprafata extrados bord fuga
loopextrados=newreg;
Line Loop(loopextrados)={ei,bf,-ee,-mbae};
ebf=newreg;
Ruled Surface(ebf)={loopextrados};
Physical Surface(3)={ebf}; // Perete solid

// Suprafata intrados bord fuga
loopintrados=newreg;
Line Loop(loopintrados)={ii,bf,-ie,-mbai};
ibf=newreg;
Ruled Surface(ibf)={loopintrados};
Physical Surface(3)={ibf}; // Perete solid

// Suprafata plana de extremitate
lpe=newreg;
Line Loop(lpe)={baextremitate,ie,-ee}; 
spe=newreg;
Plane Surface(spe)={lpe};
Physical Surface(3)={spe};

                     // FRONTIERA EXTERIOARA

le=.2; // lungimea caracteristica
raza=3.*c;
razai=raza;

 // Punctele ce definesc arcul de cerc din planul de incastrare - fizic
p1=newp; Point(p1)={0,0,raza,le};
p2=newp; Point(p2)={-raza,0,0,le};
p3=newp; Point(p3)={0,0,-raza,le};
pc=newp; Point(pc)={0,0,0,le}; // Centrul cercului

c1=newreg; Circle(c1)={p1,pc,p2};
c2=newreg; Circle(c2)={p2,pc,p3};


  // Punctele ce definesc arcul de cerc din planul de langa cel de extremitate
le=.3*ce/c;
raza=raza*ce/c;
dye = 0.7*ye; // distanta planului de frontiera fata de planul extremitatii (sunt paralele)

yc=ye+dye;
const=yc/ye; //Printf("const = %g",const);
xc=xe[1]*const;
zc=Sin(ud)*Fabs(ye)*const;
//Printf("xc, yc, zc = %g %g %g",xc,yc,zc);
pc=newp; Point(pc)={xc,yc,zc,le};

p4=newp; Point(p4)={xc,yc,zc+raza,le};
p5=newp; Point(p5)={xc-raza,yc,zc,le};
p6=newp; Point(p6)={xc,yc,zc-raza,le};

c3=newreg; Circle(c3)={p4,pc,p5};
c4=newreg; Circle(c4)={p5,pc,p6};

liniafrontala=newreg; Line(liniafrontala)={p2,p5};
liniasuperioara=newreg; Line(liniasuperioara)={p1,p4};
liniainferioara=newreg; Line(liniainferioara)={p3,p6};

//Suprafata frontala superioara
lsfs=newreg; Line Loop(lsfs)={c1,liniafrontala,-c3,-liniasuperioara};
sfs=newreg; Ruled Surface(sfs)={lsfs};
Physical Surface(5)={sfs}; // Intrare

//Suprafata frontala inferioara
lsfi=newreg; Line Loop(lsfi)={c2,liniainferioara,-c4,-liniafrontala};
sfi=newreg; Ruled Surface(sfi)={lsfi};
Physical Surface(5)={sfi}; // Intrare

lamonte=10.*c;
le=le*3;
p7=newp; Point(p7)={lamonte,0,razai,le};
p8=newp; Point(p8)={lamonte,0,-razai,le};
p78=newp; Point(p78)={lamonte,0,0,le};
p9=newp; Point(p9)={lamonte,yc,zc+raza,le};
p10=newp; Point(p10)={lamonte,yc,zc-raza,le};

l11=newreg; Line(l11)={p1,p7};
l12=newreg; Line(l12)={p3,p8};
l13=newreg; Line(l13)={p4,p9};
l14=newreg; Line(l14)={p6,p10};
l15=newreg; Line(l15)={p7,p9};
l16=newreg; Line(l16)={p8,p10};
l17=newreg; Line(l17)={p9,p10};
l18s=newreg; Line(l18s)={p7,p78};
l18i=newreg; Line(l18i)={p78,p8};
l19=newreg; Line(l19)={np,p78};
l20=newreg; Line(l20)={np2,p1};
l21=newreg; Line(l21)={np+np2-1,p3};

lplansuperior=newreg; Line Loop(lplansuperior)={l11,l15,-l13,-liniasuperioara};
plansuperior=newreg; Plane Surface(plansuperior)={lplansuperior};
Physical Surface(4)={plansuperior}; // Aval 

lplaninferior=newreg; Line Loop(lplaninferior)={l12,l16,-l14,-liniainferioara};
planinferior=newreg; Plane Surface(planinferior)={lplaninferior};
Physical Surface(4)={planinferior}; // Aval

lplanposterior=newreg; Line Loop(lplanposterior)={l15,l17,-l16,-l18i,-l18s};
planposterior=newreg; Plane Surface(planposterior)={lplanposterior};
Physical Surface(4)={planposterior}; // Aval

lplanlateral=newreg; Line Loop(lplanlateral)={c3,c4,l14,-l17,-l13};
planlateral=newreg; Plane Surface(planlateral)={lplanlateral};
Physical Surface(4)={planlateral};  // Aval

// Planul din planul de incastrare alcatuit din Spline-ul de bord de atac....
lplan1=newreg; Line Loop(lplan1)={c1,c2,-l21,-baincastrare,l20};
plan1=newreg; Plane Surface(plan1)={lplan1};
Physical Surface(3)={plan1}; // Planul de incastrare este perete solid

lplan2=newreg; Line Loop(lplan2)={l19,-l18s,-l11,-l20,ei};
plan2=newreg; Plane Surface(plan2)={lplan2};
Physical Surface(3)={plan2}; // Planul de incastrare este perete solid

lplan3=newreg; Line Loop(lplan3)={l19,l18i,-l12,-l21,ii};
plan3=newreg; Plane Surface(plan3)={lplan3};
Physical Surface(3)={plan3}; // Planul de incastrare este perete solid

lvolum=newreg; Surface Loop(lvolum)={sba,ebf,ibf,spe,sfs,sfi,plansuperior,planinferior,
planposterior,planlateral,plan1,plan2,plan3};
volum=newreg;
Volume(volum)={lvolum};
Physical Volume(0)={volum};

