SetFactory("OpenCASCADE");
sf = 0.001;

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//  UPPER SURFACE OF THE VEHICLE
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//
// ---------------------------------------------------------------
// OGIVE (Tangent ogive formula with R = 31.75mm and L = 228.6mm)
// ---------------------------------------------------------------

R=31.75*sf;
L=228.6*sf;
rho=(R^2+L^2)/(2*R);
x1=1*sf;
x2=57.15*sf;
x3=114.3*sf;
x4=171.45*sf;
x5=228.6*sf;

y1=Sqrt(rho^2-(L-x1)^2)+R-rho;  // For later use in offset
y2=Sqrt(rho^2-(L-x2)^2)+R-rho;  // For later use in offset
y4=Sqrt(rho^2-(L-x4)^2)+R-rho;  // For later use in offset

Point(0) = {0,  0, 0};
Point(1) = {x1, 0, Sqrt(rho^2-(L-x1)^2)+R-rho};
Point(2) = {x2, 0, Sqrt(rho^2-(L-x2)^2)+R-rho};
Point(3) = {x3, 0, Sqrt(rho^2-(L-x3)^2)+R-rho};
Point(4) = {x4, 0, Sqrt(rho^2-(L-x4)^2)+R-rho};
Point(5) = {x5, 0, Sqrt(rho^2-(L-x5)^2)+R-rho};

Line(1) = {0, 1};
Spline(2) = {1, 2, 3, 4, 5};

// ---------------------------------------------------------------
// FUSELAGE with wing root
// ---------------------------------------------------------------

Point(6) = {291.069988526*sf, 0, 31.75*sf};
Point(7) = {381*sf, 0, 31.75*sf};
Point(8) = {470.930011474*sf, 0, 31.75*sf};
Point(9) = {635*sf, 0, 31.75*sf};

Line(3) = {5, 6};
Line(4) = {6, 7};
Line(5) = {7, 8};
Line(6) = {8, 9};

Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} {Curve{1}; Curve{2}; Curve{5}; Curve{3}; Curve{4}; Curve{5}; Curve{6};}

Point(16) = {381*sf, 0, -1122*sf};
Point(17) = {470.930011474*sf, 0, 0};
Point(18) = {291.069988526*sf, 0, 0};

xwr=291.069988526*sf;  ywr=31.75*sf;  // to be used later for offset

Circle(20) = {17, 16, 18};
Extrude {0, 40*sf, 0} { Curve{20}; }

BooleanDifference{ Surface{3}; Surface{5}; Delete; }{ Surface{7}; Delete; }
Recursive Delete { Surface{8}; Surface{10}; Point{16};}
Coherence;

// ---------------------------------------------------------------
// BICONVEX WING
// ---------------------------------------------------------------

rothick=(470.930011474*sf - 291.069988526*sf)*0.02;

Point(23) = {381*sf, 143.891*sf, -423*sf};
Point(24) = {414.909*sf, 143.891*sf, 0};
Point(25) = {347.091*sf, 143.891*sf, 0};
Point(26) = {381*sf, 143.891*sf, 1.3569491372*sf};

xwt=347.091*sf;  ywt=143.891*sf;  // to be used later in offset

Circle(30) = {24, 23, 26};
Circle(31) = {26, 23, 25};
Recursive Delete { Point{23}; }

Line(32) = {24, 13};
Line(33) = {26, 18};
Line(34) = {25, 14};

Curve Loop(10) = {21, -33, -30, 32};
Surface(10) = {10};
Curve Loop(12) = {33, 26, -34, -31};
Surface(11) = {12};

Point(27) = {381*sf, 143.891*sf, 0*sf};

Line(35) = {24, 27};
Line(36) = {27, 26};
Line(37) = {27, 25};

Curve Loop(14) = {30, -36, -35};
Plane Surface(12) = {14};
Curve Loop(15) = {37, -31, -36};
Plane Surface(13) = {15};

// ---------------------------------------------------------------
// REAR
// ---------------------------------------------------------------

Point(28) = {635*sf, 0*sf, 0*sf};
Line(38) = {28, 9};
Line(39) = {28, 15};

Curve Loop(16) = {38, 18, -39};
Plane Surface(14) = {16};

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//  LOWER SURFACE OF THE VEHICLE
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

// Ogive including nose, free-of-wing fuselage
Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} { Curve{9}; Curve{11}; Curve{16}; Curve{19};}
Coherence;

// Rear
Line(50) = {28, 33};
Curve Loop(21) = {39, 48, -50};
Plane Surface(19) = {21};

// Fuselage wing section
Point(34) = {381*sf, 0, -31.75*sf};

Line(51) = {32, 34};
Line(52) = {34, 31};

Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Curve{51}; Curve{52}; }

Point(38) = {381*sf, 0, 1122*sf};
Point(39) = {470.930011474*sf, 0, 0};
Point(40) = {291.069988526*sf, 0, 0};
//
Circle(58) = {39, 38, 40};

Extrude {0, 40*sf, 0} { Curve{58}; }

BooleanDifference{ Surface{20}; Surface{21}; Delete; }{ Surface{22}; Delete; }
Recursive Delete { Surface{21}; Surface{23};  Point{38}; }
Coherence;

// Wing tip shape
Point(51) = {381*sf, 143.891*sf, 423*sf};
Point(52) = {381*sf, 143.891*sf, -1.3569491372*sf};

Circle(83) = {24, 51, 52};
Circle(84) = {52, 51, 25};
Recursive Delete { Point{51}; }

Line(85) = {52, 40};
Line(86) = {52, 27};

// Wing plan
Curve Loop(22) = {63, 70, -85, -83};
Surface(23) = {22};
Curve Loop(24) = {85, 72, -64, -84};
Surface(24) = {24};

// Wing tip
Curve Loop(26) = {83, 86, -35};
Surface(25) = {26};
Curve Loop(28) = {84, -37, -86};
Surface(26) = {28};

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//  VEHICLE OFFSET
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

MAC= 152.07*sf;   // Mean Aerodynamic Chord
thick= 0.02*MAC;  // Thickness as a function of MAC
cossweep=(xwt-xwr)/Sqrt((xwt-xwr)^2+(ywt-ywr)^2);  // Cos of wing sweep angle
sinsweep=(ywt-ywr)/Sqrt((xwt-xwr)^2+(ywt-ywr)^2);  // Sin of wing sweep angle

// Dealing with the pointy nose
cosalpha=x1/Sqrt(x1^2+y1^2);
sinalpha=y1/Sqrt(x1^2+y1^2);

dx= -thick*sinalpha;  dy= 0;  dz= thick*cosalpha;
Translate {dx, dy, dz} { Duplicata { Point{0}; }}
dx= -thick*sinalpha;  dy= 0;  dz= -thick*cosalpha;
Translate {dx, dy, dz} { Duplicata { Point{0}; }}
Point(55) = {381*sf, -0.5, 0};
dx= -thick;  dy= 0;  dz= 0;
Translate {dx, dy, dz} { Duplicata { Point{0}; }}

// Middle point along the spline of the ogive
cosalpha=(x4-x2)/Sqrt((x4-x2)^2+(y4-y2)^2);
sinalpha=(y4-y2)/Sqrt((x4-x2)^2+(y4-y2)^2);
dx= -thick*sinalpha;  dy= 0;  dz= thick*cosalpha;

Translate {dx, dy, dz} { Duplicata { Point{3}; }}
Symmetry {0, 0, 1, 0} { Duplicata { Point{57}; }}

// Points along the fuselage
dx= 0;  dy= 0;  dz= thick;
Translate {dx, dy, dz} { Duplicata { Point{5}; Point{17}; }}
Translate {dx, dy, -dz} { Duplicata { Point{30}; Point{39}; }}
dx= -thick;  dy= 0;  dz= thick;
Translate {dx, dy, dz} { Duplicata { Point{41}; }}
Translate {dx, dy, -dz} { Duplicata { Point{45};}}
dx= thick;  dy= 0;  dz= thick;
Translate {dx, dy, dz} { Duplicata { Point{43}; }}
Translate {dx, dy, -dz} { Duplicata { Point{46};}}

// Points on the rear cap
dx= 0;  dy= 0;  dz= thick;
Translate {dx, dy, dz} { Duplicata { Point{9}; }}
dx= 0;  dy= 0;  dz= -thick;
Translate {dx, dy, dz} { Duplicata { Point{33}; }}

dx= thick;  dy= 0;  dz= 0;
Translate {dx, dy, dz} { Duplicata { Point{9}; Point{28}; Point{33}; }} 

Recursive Delete { Point{68}; Point{71}; }

Point(68) = {-0.02, 0, 0};
Point(71) = {-0.01, 0, 0};

// Additional points near the pointy nose
cosalpha=x1/Sqrt(x1^2+y1^2);
sinalpha=y1/Sqrt(x1^2+y1^2);

dx= -thick*sinalpha;  dy= 0;  dz= thick*cosalpha;
Translate {dx, dy, dz} { Duplicata { Point{1}; }}
dx= -thick*sinalpha;  dy= 0;  dz= -thick*cosalpha;
Translate {dx, dy, dz} { Duplicata { Point{29}; }}
Point(92) = {x1, 0, 0};

// Offset lines and surfaces by revolution and boolean
// Fuselage
Spline(87) = {72, 57, 59};
Line(88) = {59, 63};
Line(89) = {63, 60};
Line(90) = {60, 65};
Line(91) = {65, 67};

Circle(92) = {67, 9, 69};
Line(93) = {1, 72};

Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} {Curve{87}; Curve{88}; 
         Curve{89}; Curve{90}; Curve{91}; Curve{92}; Curve{93};}                  

// Creating a biconvex trace on the symmetry plane to be subsequently extruded and 
// intersected with the shell of the fuselage to create the trace of wing root on the shell
Translate {0, -31.75*sf-thick, 0} { Duplicata { Point{97}; Point{95}; }}  // Creating points for LE and TE of wing root along plane of symmetry
Point(105) = {381*sf, 0, rothick+thick};
Spline(112) = {102, 105, 101};
Extrude {0, 31.75*sf+1.1*thick, 0} {Curve{112}; }

BooleanDifference{ Surface{29}; Surface{30}; Delete; }{ Surface{34}; Delete; }
Recursive Delete { Surface{35}; Surface{37}; Point{105}; }
Coherence;

// Creating the shell for the wing tip planar cap
// Point at trailing edge and leading edge of tip
alpha=55*3.1415962/180;
dx= thick*Cos(alpha);  dy= thick;  dz= 0;
Translate {dx, dy, dz} { Duplicata { Point{24};}}
dx= -thick*Cos(alpha);  dy= thick;  dz= 0;
Translate {dx, dy, dz} { Duplicata { Point{25};}}

// Points at wing tip. Creating an spline at thick from planar tip cap along y
dx= 0;  dy= thick;  dz= 0;
Translate {dx, dy, dz} { Duplicata { Point{27};}}
dx= 0;  dy= thick;  dz= thick;
Translate {dx, dy, dz} { Duplicata { Point{26};}}
dx= 0;  dy= thick;  dz= -thick;
Translate {dx, dy, dz} { Duplicata { Point{52};}}

Line(128) = {97, 112};
Line(129) = {95, 113};
Line(130) = {107, 115};
Line(131) = {115, 114};

tiphch=(414.909*sf-381*sf);
tiphth=2*tiphch*0.015;
Point(117) = {381*sf+0.3*tiphch, 143.891*sf+thick, tiphth+thick};
Point(118) = {381*sf-0.3*tiphch, 143.891*sf+thick, tiphth+thick};

Spline(132) = {112, 117, 115};
Spline(133) = {115, 118, 113};
Curve Loop(40) = {115, 130, -132, -128};
Surface(37) = {40};
Curve Loop(42) = {111, 129, -133, -130};
Surface(38) = {42};

Point(119) = {381*sf+0.3*tiphch, 143.891*sf+thick, -tiphth-thick};
Point(120) = {381*sf-0.3*tiphch, 143.891*sf+thick, -tiphth-thick};

Spline(134) = {112, 119, 116};
Spline(135) = {116, 120, 113};
Line(136) = {112, 114};
Line(137) = {114, 116};
Line(138) = {114, 113};

Curve Loop(44) = {132, 131, -136};
Plane Surface(39) = {44};
Curve Loop(45) = {133, -138, -131};
Plane Surface(40) = {45};
Curve Loop(46) = {136, 137, -134};
Plane Surface(41) = {46};
Curve Loop(47) = {138, -135, -137};
Plane Surface(42) = {47};

// Pointy nose
Line(139) = {72, 53};
Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} { Curve{139}; }
Coherence;

Circle(147) = {123, 108, 56};
Line(150) = {123, 108};
Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} { Curve{147};Curve{150}; }
Coherence;

// Finalizing rear
Line(153) = {69, 70};
Line(154) = {70, 99};

// Delete points that were created in the lower part of the nose shell - we will not need them and erasing their creation at the beggining would require total reconstruction
Recursive Delete { Point{66};Point{62};Point{64};Point{61};Point{58};Point{73};Point{54}; }

// Extruding lines along the fuselage to create lower shell
// Nose, free fuselage
Point(126) = {381*sf, 0, -31.75*sf-thick};
Extrude {{1, 0, 0}, {0, 0, 0}, -Pi/2} { Curve{141};Curve{98};Curve{104};Curve{106};Curve{146};Curve{149};Curve{152};Curve{143};}
Coherence;
Recursive Delete { Point{71};Point{68};}

// Wing
Line(176) = {130, 126};
Line(177) = {126, 129};
Extrude {{1, 0, 0}, {0, 0, 0}, Pi/2} { Curve{176};Curve{177};}

// Creating a biconvex trace on the symmetry plane to be subsequently extruded and 
// intersected with the shell of the fuselage to create the trace of wing root on the shell
Translate {0, -31.75*sf-thick, 0} { Duplicata { Point{97}; Point{95}; }}  // Creating points for LE and TE of wing root along plane of symmetry
Point(140) = {381*sf, 0, -rothick-thick};
Spline(183) = {138, 140, 139};
Extrude {0, 31.75*sf+1.1*thick, 0} {Curve{183}; }

BooleanDifference{ Surface{54}; Surface{55}; Delete; }{ Surface{56}; Delete; }
Recursive Delete { Surface{57}; Surface{55}; Point{140};  }
Coherence;

// Lower shell wing planforms
Line(198) = {142, 116};
Curve Loop(55) = {188, 134, -198, 195};
Surface(57) = {55};
Curve Loop(57) = {198, 135, -189, 197};
Surface(58) = {57};

// Pointy nose and last lines on the rear
Line(199) = {70, 132};
Line(200) = {108, 125};

// BL lines 
Line(201) = {59, 5};
Line(202) = {11, 94};
Line(203) = {30, 128};
Line(204) = {41, 143};
Line(205) = {42, 144};
Line(206) = {45, 147};
Line(207) = {106, 17};
Line(208) = {107, 18};
Line(209) = {40, 142};
Line(210) = {39, 141};
Line(211) = {43, 145};
Line(212) = {44, 146};
Line(213) = {46, 148};
Line(214) = {33, 131};
Line(215) = {33, 132};
Line(216) = {15, 98};
Line(217) = {15, 99};
Line(218) = {28, 70};
Line(219) = {9, 67};
Line(220) = {9, 69};
Line(221) = {26, 115};
Line(222) = {52, 116};
Line(223) = {24, 112};
Line(224) = {25, 113};
Line(225) = {27, 114};

// BL surfaces
// Nose
Curve Loop(59) = {145, 151, 118, 142};
Plane Surface(59) = {59};
Curve Loop(61) = {151, 200, -148};
Plane Surface(60) = {61};
Curve Loop(63) = {119, 143, 146, 152};
Plane Surface(61) = {63};
Curve Loop(64) = {152, 200, -149};
Plane Surface(62) = {64};
Curve Loop(65) = {173, 175, 166, 170};
Plane Surface(63) = {65};
Curve Loop(66) = {170, 200, -168};
Plane Surface(64) = {66};

// Ogive
Curve Loop(67) = {142, 140, 201, -120};
Plane Surface(65) = {67};
Curve Loop(69) = {202, -141, -143, 121};
Plane Surface(66) = {69};
Curve Loop(70) = {174, 203, -157, -175};
Plane Surface(67) = {70};

// Fuselage before wing
Curve Loop(71) = {203, -156, -202, 43};
Plane Surface(68) = {71};
Curve Loop(72) = {202, -95, 201, 10};
Plane Surface(69) = {72};
Curve Loop(73) = {206, -191, -203, 66};
Plane Surface(70) = {73};
Curve Loop(74) = {55, 205, -180, -202};
Plane Surface(71) = {74};
Curve Loop(75) = {204, -179, 201, 54};
Plane Surface(72) = {75};

//Fuselage with wing
Curve Loop(76) = {204, 178, -205, -53};
Surface(73) = {76};
Curve Loop(78) = {205, 190, -206, -65};
Surface(74) = {78};
Curve Loop(80) = {204, 184, 207, -61};
Plane Surface(75) = {80};
Curve Loop(81) = {207, 22, -208, -110};
Plane Surface(76) = {81};
Curve Loop(82) = {207, 59, 211, -186};
Plane Surface(77) = {82};
Curve Loop(83) = {211, 181, -212, -56};
Surface(78) = {83};
Curve Loop(85) = {192, -213, -67, 212};
Surface(79) = {85};
Curve Loop(87) = {209, -177, -210, 52};
Surface(80) = {87};
Curve Loop(105) = {213, 194, -210, -69};
Plane Surface(94) = {105};
Curve Loop(106) = {210, 196, -206, -71};
Plane Surface(95) = {106};


// Fuselage after wing
Curve Loop(89) = {211, 182, -219, -57};
Plane Surface(81) = {89};
Curve Loop(90) = {58, 216, -183, -212};
Plane Surface(82) = {90};
Curve Loop(91) = {68, 214, -193, -213};
Plane Surface(83) = {91};

// Rear
Curve Loop(92) = {219, 103, -216, -18};
Plane Surface(84) = {92};
Curve Loop(93) = {92, -220, 219};
Plane Surface(85) = {93};
Curve Loop(94) = {161, -214, -48, 216};
Plane Surface(86) = {94};
Curve Loop(96) = {215, -164, -214};
Plane Surface(87) = {96};
Curve Loop(97) = {215, -199, -218, 50};
Plane Surface(88) = {97};
Curve Loop(98) = {218, -153, -220, -38};
Plane Surface(89) = {98};
Curve Loop(99) = {218, 154, -217, -39};
Plane Surface(90) = {99};
Curve Loop(100) = {106, -217, 216};
Plane Surface(91) = {100};
Curve Loop(101) = {105, -217, -18, 220};
Surface(92) = {101};
Curve Loop(103) = {217, 163, -215, -48};
Surface(93) = {103};

// Wing
Curve Loop(110) = {187, 208, -60, 212};
Surface(96) = {110};
Curve Loop(112) = {212, -195, -209, -70};
Surface(97) = {112};
Curve Loop(114) = {209, -197, -205, -72};
Surface(98) = {114};
Curve Loop(116) = {205, -185, 208, 62};
Surface(99) = {116};
Curve Loop(118) = {208, -33, 221, -130};
Surface(100) = {118};
Curve Loop(120) = {221, 131, -225, 36};
Plane Surface(101) = {120};
Curve Loop(121) = {225, 137, -222, 86};
Plane Surface(102) = {121};
Curve Loop(122) = {222, -198, -209, -85};
Plane Surface(103) = {122};
Curve Loop(123) = {63, 212, 188, -223};
Plane Surface(104) = {123};
Curve Loop(124) = {64, 205, 189, -224};
Plane Surface(105) = {124};
Curve Loop(126) = {223, 132, -221, -30};
Surface(106) = {126};
Curve Loop(128) = {221, 133, -224, -31};
Surface(107) = {128};
Curve Loop(130) = {224, -135, -222, 84};
Surface(108) = {130};
Curve Loop(132) = {222, -134, -223, 83};
Surface(109) = {132};

// BL Volumes
Surface Loop(1) = {1, 59, 43, 33, 45, 61};
Volume(1) = {1};
Surface Loop(2) = {45, 60, 44, 62};
Volume(2) = {2};
Surface Loop(3) = {53, 50, 63, 15, 52, 61};
Volume(3) = {3};
Surface Loop(4) = {52, 51, 64, 62};
Volume(4) = {4};
Surface Loop(5) = {33, 27, 65, 2, 69, 66};
Volume(5) = {5};
Surface Loop(6) = {66, 67, 16, 46, 68, 53};
Volume(6) = {6};
Surface Loop(7) = {71, 73, 72, 28, 4, 69};
Volume(7) = {7};
Surface Loop(8) = {74, 71, 47, 70, 17, 68};
Volume(8) = {8};
Surface Loop(9) = {75, 34, 9, 73, 99, 76};
Volume(9) = {9};
Surface Loop(10) = {76, 77, 7, 36, 96, 78};
Volume(10) = {10};
Surface Loop(11) = {74, 56, 95, 22, 80, 98};
Volume(11) = {11};
Surface Loop(12) = {94, 54, 20, 80, 97, 79};
Volume(12) = {12};
Surface Loop(13) = {79, 83, 18, 48, 82, 86};
Volume(13) = {13};
Surface Loop(14) = {84, 81, 31, 6, 82, 78};
Volume(14) = {14};

// Missing surfaces at rear
Curve Loop(160) = {153, 154, -105};
Plane Surface(136) = {160};
Curve Loop(161) = {199, -163, -154};
Plane Surface(137) = {161};

// BL Volumes
Surface Loop(15) = {89, 136, 14, 92, 90};
Volume(15) = {15};
Surface Loop(16) = {88, 137, 19, 90, 93};
Volume(16) = {16};
Surface Loop(17) = {87, 49, 93, 91, 86};
Volume(17) = {17};
Surface Loop(18) = {92, 91, 32, 85, 84};
Volume(18) = {18};
Surface Loop(19) = {96, 100, 104, 37, 10, 106};
Volume(19) = {19};
Surface Loop(20) = {23, 57, 104, 103, 109, 97};
Volume(20) = {20};
Surface Loop(21) = {24, 58, 105, 103, 98, 108};
Volume(21) = {21};
Surface Loop(22) = {38, 11, 105, 100, 99, 107};
Volume(22) = {22};

// Missing surfaces at wing tip
Curve Loop(164) = {224, -138, -225, 37};
Plane Surface(140) = {164};
Curve Loop(165) = {225, -136, -223, 35};
Plane Surface(141) = {165};

// BL Volumes
Surface Loop(23) = {101, 140, 107, 40, 13};
Volume(23) = {23};
Surface Loop(24) = {140, 108, 102, 42, 26};
Volume(24) = {24};
Surface Loop(25) = {101, 141, 106, 39, 12};
Volume(25) = {25};
Surface Loop(26) = {141, 25, 41, 109, 102};
Volume(26) = {26};

Coherence;

Delete { Surface{33}; Surface{52}; Surface{50}; Surface{53}; Surface{15}; Surface{1}; Surface{2}; Surface{27};}
Delete { Surface{44}; Surface{51}; Surface{43}; Surface{45}; Surface{47}; Surface{17}; Surface{4}; Surface{28};}
Delete { Surface{46}; Surface{16}; Surface{31}; Surface{32}; Surface{49}; Surface{48}; Surface{18}; Surface{6};}

// ---------------------------------------------------------------
// Near-Field
// ---------------------------------------------------------------

//x0 = 381*sf;  z0 = 0*sf;
//maxA = 635*sf;  minA = 400*sf;
//
alphaD = -20;
alpha = alphaD*3.1415962/180;
//
//x1 = maxA;    z1 = 0;
//x2 = 0;       z2 = + minA;
//x3 = - maxA;  z3 = 0;
//x4 = 0;       z4 = - minA;
//x5 = x0;      y5 = + minA;
//
//Point(149) = {x1+x0,  0, z1+z0, 10*sf};
//Point(150) = {x2+x0,  0, z2+z0, 10*sf};
//Point(151) = {x3+x0,  0, z3+z0, 10*sf};
//Point(152) = {x0, 0, 0};
//
//Ellipse(226) = {151, 152, 151, 150};
//Ellipse(227) = {149, 152, 149, 150};
//
//Extrude {{1, 0, 0}, {0, 0, 0}, -Pi} {Curve{226}; Curve{227};}
//
//Rotate {{0, 1, 0}, {381*sf, 0, 0}, alpha} {
//  Surface{153}; Surface{152}; Curve{227}; Curve{226}; Curve{232}; Curve{230}; 
//}
//
//Curve Loop(179) = {233, -236, 237, -234};
//Curve Loop(180) = {148, -168, -166, 157, 191, -196, -194, 193, 164, -199, -153, -92, -182, -186, -184, -179, -140, 145};
//Plane Surface(154) = {179, 180};
//Surface Loop(27) = {154, 153, 152, 134, 139, 136, 137, 138, 133, 36, 54, 34, 56, 130, 128, 146, 151, 111, 143, 117, 144, 38, 58, 41, 42, 40, 39, 57, 37};
//
//Volume(27) = {27};
//
//Coherence;
//
//Delete { Surface{153}; Surface{152}; }
//
// -----------------------------------------------------------//----
// FAR-FIELD
// -----------------------------------------------------------//----
x6 = 2500*sf; z6 = 0*sf;
maxB = 3000*sf;  minB = 1500*sf;
//
x7 = - maxB;  z7 = 0;
x8 = 0;       z8 = minB;
x9 = 0;       z9 = - minB;
x10 = x6;     y10 = + minB;

Point(158) = {x6, 0, z6, 50*sf};
Point(159) = {x7+x6,  0, z7+z6, 300*sf};
Point(160) = {x8+x6,  0, z8+z6, 50*sf};

Ellipse(238) = {159, 158, 159, 160};

Extrude {{1, 0, 0}, {0, 0, 0}, -Pi} {Curve{238}; }

Line(242) = {160, 158};
Line(243) = {158, 161};

//Rotate {{0, 1, 0}, {381*sf, 0, 0}, alpha} {
//   Surface{157}; Curve{238}; Curve{241}; Curve{243}; Curve{240}; Curve{242}; 
//}

Rotate {{0, 1, 0}, {381*sf, 0, 0}, alpha} {
   Surface{152}; Curve{238}; Curve{241}; Curve{243}; Curve{240}; Curve{242}; 
}

Curve Loop(177) = {242, 243, -245};
//+
Plane Surface(153) = {177};

Curve Loop(178) = {246, 242, 243, -247};
//+
Curve Loop(179) = {193, 164, -199, -153, -92, -182, -186, -184, -179, -140, 145, 148, -168, -166, 157, 191, -196, -194};
//+
Plane Surface(154) = {178, 179};
//+
Surface Loop(27) = {152, 154, 153, 134, 139, 136, 137, 138, 133, 54, 36, 34, 56, 128, 130, 146, 151, 111, 143, 117, 144, 37, 57, 58, 38, 41, 39, 42, 40};
//+
Volume(27) = {27};


//Curve Loop(186) = {245, -243, -242};
//Plane Surface(158) = {186};
//Curve Loop(188) = {242, 243, -247, 246};
//Curve Loop(189) = {233, -236, 237, -234};
//Plane Surface(159) = {188, 189};
//Surface Loop(28) = {158, 157, 159, 155, 156};
//
//Volume(28) = {28};
//
Delete { Surface{152}; }


//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//  MESH DEFINITIONS
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
coarse= -0.3;
medium= +0.0;
fine= +0.3;
ref= coarse;

//ref
//ogivep=90+Ceil(ref*90);         // Ogive longitudinally
//chordwise=15+Ceil(ref*15);      // Half of the chord
//wingspan=30+Ceil(ref*30);       // Wing leading edge, trailing edge and mid-chord
//fusbefor=10+Ceil(ref*10);       // Fuselage longitudinally before wing
//fusafter=25+Ceil(ref*25);       // Fuselage longitudinally after wing
//circumf=35+Ceil(ref*35);        // Half of circumferential points
//wingtip=7+Ceil(ref*7);          // Wing tip cap
//rearcap=10+Ceil(ref*10);        // Rear cap
//blpoints=25+Ceil(ref*25);       // Points in the boundary layer
//nosetip=20+Ceil(ref*20);
//

ogivep=150+Ceil(ref*150);         // Ogive longitudinally
chordwise=150+Ceil(ref*200);      // Half of the chord
wingspan=100+Ceil(ref*150);       // Wing leading edge, trailing edge and mid-chord
fusbefor=50+Ceil(ref*50);       // Fuselage longitudinally before wing
fusafter=100+Ceil(ref*100);       // Fuselage longitudinally after wing
circumf=70+Ceil(ref*70);        // Half of circumferential points
wingtip=15+Ceil(ref*15);          // Wing tip cap
rearcap=20+Ceil(ref*20);        // Rear cap
blpoints=25+Ceil(ref*25);       // Points in the boundary layer
nosetip=40+Ceil(ref*40);
//
dprog=0.15;
blprogrF=1+dprog;   // Progression in BL in forward direction
blprogrB=1-dprog;   // Progression in BL in backward direction

dprogcirc1 = 0; //0.09;
dprogcirc1F=1+dprogcirc1;   // Progression in circumferential in forward direction
dprogcirc1B=1-dprogcirc1;   // Progression in circumferential in backward direction

dprogcirc2 = 0; //0.15;
dprogcirc2F=1+dprogcirc2;   // Progression in circumferential in forward direction
dprogcirc2B=1-dprogcirc2;   // Progression in circumferential in backward direction

dproglong = 0; //0;
dproglongF=1+dproglong;   // Progression in longitudinal fuselage-wing in forward direction
dproglongB=1-dproglong;   // Progression in longitudinal fuselage-wing in backward direction

dprogch = 0.05;
dprogchF=1+dprogch;   // Progression in chord in forward direction
dprogchB=1-dprogch;   // Progression in chord in backward direction

dprogws = 0;
dprogwsF=1+dprogws;   // Progression in WS in forward direction
dprogwsB=1-dprogws;   // Progression in WS in backward direction

// Ogive vehicle and shell
Transfinite Curve {140, 120, 121, 141, 174, 157} = ogivep Using Progression 1.05;
Transfinite Curve {95, 10, 156, 43, 139, 117, 155, 139, 172, 144, 165} = circumf Using Progression 1;
Transfinite Surface {146};
Transfinite Surface {148};
Transfinite Surface {149};
Transfinite Surface {151};
Recombine Surface {146, 148, 149, 151};





Transfinite Curve {192, 67, 65, 190} = circumf Using Progression 1; //1.1;
Transfinite Curve {181, 56, 53, 178} = circumf Using Progression 1; //0.9;




// Fuselage before wing and shell
Transfinite Curve { 180, 55} = fusbefor Using Progression 0.9;
Transfinite Curve {179, 54, 66, 191} = fusbefor Using Progression 0.95;
Transfinite Surface {128};
Transfinite Surface {129};
Transfinite Surface {131};
Transfinite Surface {130};
Recombine Surface {128, 129, 131, 130};



// Fuselage at wing and shell
Transfinite Curve {177, 52, 110, 22} = circumf Using Progression dprogcirc1B;
Transfinite Curve {197} = chordwise Using Progression dprogchF;
Transfinite Curve {185, 62, 72} = chordwise Using Progression dprogchB;
Transfinite Curve {187, 60, 70} = chordwise Using Progression dprogchF;
Transfinite Curve {195} = chordwise Using Progression dprogchB;

Transfinite Curve { 69, 194} = chordwise Using Progression dproglongF;
Transfinite Curve {186, 59} = chordwise Using Progression dproglongB;
Transfinite Curve {184, 61} = chordwise Using Progression dproglongF;
Transfinite Curve {71, 196} = chordwise Using Progression dproglongB;


Transfinite Surface {34};
Transfinite Surface {9};
Transfinite Surface {36};
Transfinite Surface {7};
Transfinite Surface {20};
Transfinite Surface {54};
Transfinite Surface {22};
Transfinite Surface {56};
Recombine Surface {36, 7, 9, 34, 56, 22, 20, 54};

// Fuselage after wing, rear cap and shells
Transfinite Curve {18, 103, 105, 48, 161, 163} = circumf Using Progression 1;

Transfinite Curve { 183, 58} = fusafter Using Bump 0.02;
Transfinite Curve { 182, 57, 68, 193} = fusafter Using Bump 0.07;

Transfinite Curve {38, 154, 39, 199, 50} = rearcap Using Progression 0.9;
Transfinite Curve {153} = rearcap Using Progression 1.2;
Transfinite Curve {92, 106, 164} = 10 Using Progression 1.2;
Transfinite Surface {134};
Transfinite Surface {135};
Transfinite Surface {132};
Transfinite Surface {133};
Transfinite Surface {139};
Transfinite Surface {138};
Recombine Surface {134, 135, 132, 133, 138, 139};
Transfinite Surface {14} = {28, 15, 9};
Transfinite Surface {136} = {70, 69, 99};
Recombine Surface {14, 136};
Transfinite Surface {137} = {70, 132, 99};
Transfinite Surface {19} = {28, 33, 15};
Recombine Surface {137, 19};


// Wing and shells
Transfinite Curve {131, 137, 36, 86} = wingtip Using Progression 1;
Transfinite Curve {63, 33, 64} = wingspan Using Bump 0.005; 
Transfinite Curve {188, 189} = wingspan Using Bump 0.005;
Transfinite Curve {33, 130, 85, 198} = wingspan Using Bump 0.1;

Transfinite Curve { 133, 135, 138, 31, 84, 37} = chordwise Using Progression dprogchB;
Transfinite Curve {132, 134, 30, 83, 136, 35} = chordwise Using Progression dprogchF;
Transfinite Surface {37};
Transfinite Surface {10};
Transfinite Surface {23};
Transfinite Surface {57};
Transfinite Surface {38};
Transfinite Surface {11};
Transfinite Surface {24};
Transfinite Surface {58};
Recombine Surface {37, 10, 23, 57, 58, 24, 11, 38};
Transfinite Surface {40} = {113, 115, 114};
Transfinite Surface {39} = {112, 115, 114};
Transfinite Surface {41} = {112, 116, 114};
Transfinite Surface {42} = {113, 116, 114};
Recombine Surface {39, 41, 40, 42};
Transfinite Surface {25} = {24, 52, 27};
Transfinite Surface {26} = {25, 52, 27};
Transfinite Surface {13} = {25, 26, 27};
Transfinite Surface {12} = {24, 26, 27};
Recombine Surface {25, 12, 13, 26};

// Pointy nose
Transfinite Curve {145, 146, 166} = 5 Using Progression 1;
Transfinite Curve {118, 119, 173} = nosetip Using Bump 0.2;
Transfinite Curve {148, 149, 168} = 20 Using Progression 1;
Transfinite Curve {200} = blpoints Using Progression 1;//blprogrF;

// All lines and surfaces across the BL
Transfinite Curve { 175, 142, 143, 203, 202,206, 205, 204,  209, 210, 213, 211, 212, 219, 220, 218, 217, 216, 215, 214, 224, 221, 223, 222, 225} = blpoints Using Progression blprogrF;
Transfinite Curve {152, 201, 207, 208} = blpoints Using Progression blprogrB;
Transfinite Curve {200} = blpoints Using Progression 1.15;
Transfinite Curve {170, 151} = blpoints Using Progression 0.85;
Transfinite Surface {76};
Transfinite Surface {80};
Transfinite Surface {112};
Transfinite Surface {116};
Transfinite Surface {66};
Transfinite Surface {147};
Transfinite Surface {150};
Transfinite Surface {68};
Transfinite Surface {69};
Recombine Surface {68, 69, 150, 147, 66, 116, 112};

Transfinite Surface {72};
Transfinite Surface {71};
Transfinite Surface {70};
Transfinite Surface {74};
Transfinite Surface {73};
Transfinite Surface {75};
Transfinite Surface {95};
Transfinite Surface {99};
Transfinite Surface {98};
Transfinite Surface {94};
Transfinite Surface {97};
Transfinite Surface {77};
Transfinite Surface {96};
Transfinite Surface {78};
Transfinite Surface {79};

Recombine Surface {72, 71, 70, 74, 73, 75, 99, 95, 98, 97, 94, 80, 76, 77, 96, 78, 79};
Transfinite Surface {83};
Transfinite Surface {82};
Transfinite Surface {81};
Transfinite Surface {84};
Transfinite Surface {86};
Recombine Surface {83, 82, 81, 84, 86};

Transfinite Surface {88};
Transfinite Surface {87} = {33, 131, 132};
Transfinite Surface {89};
Transfinite Surface {85} = {9, 69, 67};
Transfinite Surface {90};
Transfinite Surface {91} = {15, 98, 99};
Transfinite Surface {92};
Transfinite Surface {93};
Recombine Surface {88, 89, 85, 92, 93, 91, 90, 87};

Transfinite Surface {104};
Transfinite Surface {105};
Transfinite Surface {100};
Transfinite Surface {103};
Transfinite Surface {107};
Transfinite Surface {140};
Transfinite Surface {108};
Transfinite Surface {109};
Transfinite Surface {141};
Transfinite Surface {106};
Transfinite Surface {101};
Transfinite Surface {102};
Recombine Surface {104, 100, 103, 105, 108, 140, 107, 109, 141, 106, 101, 102};

Transfinite Volume{5};
Transfinite Volume{6};
Transfinite Volume{8};
Transfinite Volume{7};
Transfinite Volume{11};
Transfinite Volume{12};
Transfinite Volume{9};
Transfinite Volume{10};
Transfinite Volume{14};
Transfinite Volume{13};
Transfinite Volume{16};
Transfinite Volume{15};
Transfinite Volume{18};
Transfinite Volume{17};
Transfinite Volume{19};
Transfinite Volume{20};
Transfinite Volume{22};
Transfinite Volume{21};
Transfinite Volume{26};
Transfinite Volume{25};
Transfinite Volume{23};
Transfinite Volume{24};

// Near-field
//Transfinite Curve {232} = 180+Ceil(ref*180) Using Progression 1;
//Transfinite Curve {236, 237} = 120+Ceil(ref*120) Using Progression 0.98;
//Transfinite Curve {233, 234}     = 180+Ceil(ref*180) Using Progression 1;
//
//// Far-field
//Transfinite Curve {246,247} = 100+Ceil(ref*100) Using Progression 0.98;
//Transfinite Curve {245}     = 60+Ceil(ref*60) Using Progression 1;
//Transfinite Curve {242, 243}         = 90+Ceil(ref*90) Using Progression 1;

Physical Surface("Vehicle") = {14, 19, 135, 132, 7, 20, 22, 9, 131, 129, 149, 148, 118, 110, 10, 23, 24, 11, 13, 26, 12, 25};
Physical Surface("Symmetry") = {159, 154, 81, 77, 75, 72, 147, 59, 142, 63, 145, 150, 70, 95, 94, 83, 87, 88, 89, 85};
Physical Surface("Exit") = {158};
Physical Surface("Farfield") = {160};
Physical Volume("Fluid") = {28, 27, 25, 26, 23, 24, 19, 20, 22, 21, 18, 15, 14, 17, 16, 10, 13, 9, 12, 7, 11, 8, 5, 6, 1, 2, 3, 4};
//+

/*
Field[1] = Attractor;
Field[1].FacesList = {136,137,139,138,133,134,36,54,56,34,128,130,151,146,111,117,143,144,37,57,38,58,40,42,39,41}; 
////
Field[2] = Threshold;
Field[2].IField = 1;
Field[2].LcMin = 0.005;
Field[2].LcMax = 1;
Field[2].DistMin = 0.03;// > lcmin
Field[2].DistMax = 1;//100.0;
Field[2].Sigmoid = 1;

Background Field = 2;

// Don't extend the elements sizes from the boundary inside the domain
Mesh.CharacteristicLengthExtendFromBoundary = 0;
*/
Field[3] = Attractor;
Field[3].FacesList = {136,137,139,138,133,134,36,54,56,34,128,130,151,146,111,117,143,144,37,57,38,58,40,42,39,41};
//
Field[4] = Threshold;
Field[4].IField = 3;
Field[4].LcMin = 0.0025;
Field[4].LcMax = 0.1;
Field[4].DistMin = 0;
Field[4].DistMax = 1.5;
Field[4].Sigmoid = 1;

//Field[5] = Min;
//Field[5].FieldsList = {2,4};

Background Field = 4;

// Don't extend the elements sizes from the boundary inside the domain
Mesh.CharacteristicLengthExtendFromBoundary = 0;//+
//+

//+

//+
