diff --git a/sdpm/_src/sdpmw.i b/sdpm/_src/sdpmw.i
index 8a384a2ba996e27f0ef40b14bdd0cd21d4114ef8..b7342f7619562a564a42a888a72bbc38de687bc2 100644
--- a/sdpm/_src/sdpmw.i
+++ b/sdpm/_src/sdpmw.i
@@ -51,6 +51,7 @@ threads="1"
 %include "std_string.i"
 %include "std_vector.i"
 %include "std_map.i"
+%include "std_complex.i"
 
 // EXCEPTION
 // from: http://swig.10945.n7.nabble.com/Trapping-Swig-DirectorException-td6013.html
@@ -86,6 +87,15 @@ threads="1"
    }
 }
 
+// SDPM templates
+%template(std_vector_double) std::vector<double>;
+%template(std_vector_complex) std::vector<std::complex<double>>;
+%template(std_vector_vector_complex) std::vector<std::vector<std::complex<double>>>;
+%template(std_vector_Nodep) std::vector<sdpm::Node *>;
+%template(std_vector_Elementp) std::vector<sdpm::Element *>;
+%template(std_vector_Tagp) std::vector<sdpm::Tag *>;
+%template(std_map_string_Tagp) std::map<std::string, sdpm::Tag *>;
+
 // SDPM
 %include "sdpm.h"
 %include "sdpmObject.h"
@@ -110,12 +120,32 @@ threads="1"
 %include "sdpmCppBuf2Py.h"
 %warnfilter(+401);
 
-// SDPM template
-%template(std_vector_Nodep) std::vector<sdpm::Node *>;
-%template(std_vector_Elementp) std::vector<sdpm::Element *>;
-%template(std_vector_Tagp) std::vector<sdpm::Tag *>;
-%template(std_map_string_Tagp) std::map<std::string, sdpm::Tag *>;
+// Interface to SDPM types
+class sdpmVector3d {};
+
+%extend sdpmVector3d {
+    double getitem(int i)
+    {
+        return (*self)(i);
+    }
+    void setitem(int i, double a)
+    {
+        (*self)(i) = a;
+    }
+    std::string __str__()
+    {
+        std::ostringstream out; out << *self;  return out.str();
+    }
+
+    %pythoncode {
+def __getitem__(self, idx): 
+    return self.getitem(idx)
+def __setitem__(self, idx, v): 
+    return self.setitem(idx, v)
+    }
+}
 
+// Interface to SDPM classes
 namespace sdpm {
 
 // So that print(obj) use << *obj
diff --git a/sdpm/models/lann.geo b/sdpm/models/lann.geo
new file mode 100644
index 0000000000000000000000000000000000000000..e61924efb70491dc11f292565e85a43a8dbfdfed
--- /dev/null
+++ b/sdpm/models/lann.geo
@@ -0,0 +1,1730 @@
+/******************************************/
+/* Gmsh geometry for                 lann */
+/* Generated by                 geoGen.py */
+/* Adrien Crovato                         */
+/* ULiege, 2018-2019                      */
+/******************************************/
+
+// --- Wing geometry ---
+// Number of spanwise stations: 8
+// Spanwise stations normalized coordinate: 0.000000 0.200000 0.325000 0.475000 0.650000 0.825000 0.950000 1.000000 
+// Chord lengths: 0.360600 0.317650 0.290710 0.258060 0.220290 0.182350 0.155340 0.144450 
+// Half-wing area: 0.252692
+// Half-wing span: 1.000000
+// Full-wing aspect ratio: 7.914766
+
+// --- Options ---
+DefineConstant[ nC = { 50, Name "Number of panel along chord" }
+                bC = { 0.05, Name "Progression along chord" }
+                rS = { 1, Name "Refinement factor on span" } ];
+
+
+// --- Wing points ---
+// -- Airfoil 0
+Point(1) = {0.360600,0.000000,-0.008200};
+Point(2) = {0.360513,0.000000,-0.008196};
+Point(3) = {0.360247,0.000000,-0.008176};
+Point(4) = {0.359802,0.000000,-0.008142};
+Point(5) = {0.359180,0.000000,-0.008095};
+Point(6) = {0.358380,0.000000,-0.008036};
+Point(7) = {0.357405,0.000000,-0.007968};
+Point(8) = {0.356257,0.000000,-0.007892};
+Point(9) = {0.354935,0.000000,-0.007809};
+Point(10) = {0.353441,0.000000,-0.007721};
+Point(11) = {0.351775,0.000000,-0.007597};
+Point(12) = {0.349941,0.000000,-0.007478};
+Point(13) = {0.347939,0.000000,-0.007369};
+Point(14) = {0.345771,0.000000,-0.007274};
+Point(15) = {0.343440,0.000000,-0.007201};
+Point(16) = {0.340948,0.000000,-0.007154};
+Point(17) = {0.338298,0.000000,-0.007140};
+Point(18) = {0.335491,0.000000,-0.007162};
+Point(19) = {0.332532,0.000000,-0.007224};
+Point(20) = {0.329423,0.000000,-0.007330};
+Point(21) = {0.326166,0.000000,-0.007484};
+Point(22) = {0.322765,0.000000,-0.007690};
+Point(23) = {0.319223,0.000000,-0.007946};
+Point(24) = {0.315545,0.000000,-0.008253};
+Point(25) = {0.311733,0.000000,-0.008608};
+Point(26) = {0.307791,0.000000,-0.009011};
+Point(27) = {0.303724,0.000000,-0.009461};
+Point(28) = {0.299534,0.000000,-0.009954};
+Point(29) = {0.295227,0.000000,-0.010489};
+Point(30) = {0.290807,0.000000,-0.011062};
+Point(31) = {0.286278,0.000000,-0.011670};
+Point(32) = {0.281644,0.000000,-0.012310};
+Point(33) = {0.276909,0.000000,-0.012979};
+Point(34) = {0.272080,0.000000,-0.013672};
+Point(35) = {0.267160,0.000000,-0.014386};
+Point(36) = {0.262154,0.000000,-0.015117};
+Point(37) = {0.257068,0.000000,-0.015861};
+Point(38) = {0.251906,0.000000,-0.016612};
+Point(39) = {0.246673,0.000000,-0.017364};
+Point(40) = {0.241374,0.000000,-0.018112};
+Point(41) = {0.236016,0.000000,-0.018848};
+Point(42) = {0.230602,0.000000,-0.019567};
+Point(43) = {0.225139,0.000000,-0.020261};
+Point(44) = {0.219631,0.000000,-0.020924};
+Point(45) = {0.214085,0.000000,-0.021548};
+Point(46) = {0.208505,0.000000,-0.022125};
+Point(47) = {0.202898,0.000000,-0.022649};
+Point(48) = {0.197268,0.000000,-0.023116};
+Point(49) = {0.191621,0.000000,-0.023526};
+Point(50) = {0.185963,0.000000,-0.023874};
+Point(51) = {0.180300,0.000000,-0.024160};
+Point(52) = {0.174637,0.000000,-0.024381};
+Point(53) = {0.168979,0.000000,-0.024540};
+Point(54) = {0.163332,0.000000,-0.024641};
+Point(55) = {0.157702,0.000000,-0.024685};
+Point(56) = {0.152095,0.000000,-0.024675};
+Point(57) = {0.146515,0.000000,-0.024615};
+Point(58) = {0.140969,0.000000,-0.024504};
+Point(59) = {0.135461,0.000000,-0.024342};
+Point(60) = {0.129998,0.000000,-0.024128};
+Point(61) = {0.124584,0.000000,-0.023861};
+Point(62) = {0.119226,0.000000,-0.023540};
+Point(63) = {0.113927,0.000000,-0.023167};
+Point(64) = {0.108694,0.000000,-0.022741};
+Point(65) = {0.103532,0.000000,-0.022264};
+Point(66) = {0.098446,0.000000,-0.021736};
+Point(67) = {0.093440,0.000000,-0.021157};
+Point(68) = {0.088520,0.000000,-0.020531};
+Point(69) = {0.083691,0.000000,-0.019860};
+Point(70) = {0.078956,0.000000,-0.019146};
+Point(71) = {0.074322,0.000000,-0.018394};
+Point(72) = {0.069793,0.000000,-0.017606};
+Point(73) = {0.065372,0.000000,-0.016787};
+Point(74) = {0.061066,0.000000,-0.015943};
+Point(75) = {0.056876,0.000000,-0.015080};
+Point(76) = {0.052809,0.000000,-0.014202};
+Point(77) = {0.048867,0.000000,-0.013314};
+Point(78) = {0.045055,0.000000,-0.012423};
+Point(79) = {0.041377,0.000000,-0.011530};
+Point(80) = {0.037835,0.000000,-0.010640};
+Point(81) = {0.034434,0.000000,-0.009757};
+Point(82) = {0.031177,0.000000,-0.008884};
+Point(83) = {0.028068,0.000000,-0.008022};
+Point(84) = {0.025109,0.000000,-0.007169};
+Point(85) = {0.022302,0.000000,-0.006325};
+Point(86) = {0.019652,0.000000,-0.005491};
+Point(87) = {0.017160,0.000000,-0.004665};
+Point(88) = {0.014829,0.000000,-0.003851};
+Point(89) = {0.012661,0.000000,-0.003050};
+Point(90) = {0.010659,0.000000,-0.002265};
+Point(91) = {0.008825,0.000000,-0.001498};
+Point(92) = {0.007159,0.000000,-0.000749};
+Point(93) = {0.005665,0.000000,-0.000001};
+Point(94) = {0.004342,0.000000,0.000764};
+Point(95) = {0.003194,0.000000,0.001565};
+Point(96) = {0.002220,0.000000,0.002420};
+Point(97) = {0.001422,0.000000,0.003344};
+Point(98) = {0.000801,0.000000,0.004327};
+Point(99) = {0.000356,0.000000,0.005353};
+Point(100) = {0.000089,0.000000,0.006406};
+Point(101) = {0.000000,0.000000,0.007472};
+Point(102) = {0.000089,0.000000,0.008535};
+Point(103) = {0.000356,0.000000,0.009582};
+Point(104) = {0.000801,0.000000,0.010599};
+Point(105) = {0.001422,0.000000,0.011576};
+Point(106) = {0.002220,0.000000,0.012497};
+Point(107) = {0.003194,0.000000,0.013354};
+Point(108) = {0.004342,0.000000,0.014150};
+Point(109) = {0.005665,0.000000,0.014892};
+Point(110) = {0.007159,0.000000,0.015585};
+Point(111) = {0.008825,0.000000,0.016236};
+Point(112) = {0.010659,0.000000,0.016849};
+Point(113) = {0.012661,0.000000,0.017424};
+Point(114) = {0.014829,0.000000,0.017956};
+Point(115) = {0.017160,0.000000,0.018442};
+Point(116) = {0.019652,0.000000,0.018878};
+Point(117) = {0.022302,0.000000,0.019263};
+Point(118) = {0.025109,0.000000,0.019597};
+Point(119) = {0.028068,0.000000,0.019887};
+Point(120) = {0.031177,0.000000,0.020135};
+Point(121) = {0.034434,0.000000,0.020345};
+Point(122) = {0.037835,0.000000,0.020521};
+Point(123) = {0.041377,0.000000,0.020667};
+Point(124) = {0.045055,0.000000,0.020785};
+Point(125) = {0.048867,0.000000,0.020877};
+Point(126) = {0.052809,0.000000,0.020948};
+Point(127) = {0.056876,0.000000,0.020999};
+Point(128) = {0.061066,0.000000,0.021032};
+Point(129) = {0.065372,0.000000,0.021045};
+Point(130) = {0.069793,0.000000,0.021038};
+Point(131) = {0.074322,0.000000,0.021013};
+Point(132) = {0.078956,0.000000,0.020968};
+Point(133) = {0.083691,0.000000,0.020903};
+Point(134) = {0.088520,0.000000,0.020816};
+Point(135) = {0.093440,0.000000,0.020708};
+Point(136) = {0.098446,0.000000,0.020578};
+Point(137) = {0.103532,0.000000,0.020425};
+Point(138) = {0.108694,0.000000,0.020249};
+Point(139) = {0.113927,0.000000,0.020053};
+Point(140) = {0.119226,0.000000,0.019836};
+Point(141) = {0.124584,0.000000,0.019600};
+Point(142) = {0.129998,0.000000,0.019347};
+Point(143) = {0.135461,0.000000,0.019074};
+Point(144) = {0.140969,0.000000,0.018781};
+Point(145) = {0.146515,0.000000,0.018466};
+Point(146) = {0.152095,0.000000,0.018129};
+Point(147) = {0.157702,0.000000,0.017769};
+Point(148) = {0.163332,0.000000,0.017384};
+Point(149) = {0.168979,0.000000,0.016974};
+Point(150) = {0.174637,0.000000,0.016537};
+Point(151) = {0.180300,0.000000,0.016074};
+Point(152) = {0.185963,0.000000,0.015582};
+Point(153) = {0.191621,0.000000,0.015065};
+Point(154) = {0.197268,0.000000,0.014524};
+Point(155) = {0.202898,0.000000,0.013962};
+Point(156) = {0.208505,0.000000,0.013380};
+Point(157) = {0.214085,0.000000,0.012782};
+Point(158) = {0.219631,0.000000,0.012168};
+Point(159) = {0.225139,0.000000,0.011539};
+Point(160) = {0.230602,0.000000,0.010896};
+Point(161) = {0.236016,0.000000,0.010240};
+Point(162) = {0.241374,0.000000,0.009571};
+Point(163) = {0.246673,0.000000,0.008893};
+Point(164) = {0.251906,0.000000,0.008207};
+Point(165) = {0.257068,0.000000,0.007517};
+Point(166) = {0.262154,0.000000,0.006824};
+Point(167) = {0.267160,0.000000,0.006131};
+Point(168) = {0.272080,0.000000,0.005440};
+Point(169) = {0.276909,0.000000,0.004752};
+Point(170) = {0.281644,0.000000,0.004068};
+Point(171) = {0.286278,0.000000,0.003389};
+Point(172) = {0.290807,0.000000,0.002717};
+Point(173) = {0.295227,0.000000,0.002055};
+Point(174) = {0.299534,0.000000,0.001404};
+Point(175) = {0.303724,0.000000,0.000767};
+Point(176) = {0.307791,0.000000,0.000145};
+Point(177) = {0.311733,0.000000,-0.000459};
+Point(178) = {0.315545,0.000000,-0.001045};
+Point(179) = {0.319223,0.000000,-0.001611};
+Point(180) = {0.322765,0.000000,-0.002156};
+Point(181) = {0.326166,0.000000,-0.002681};
+Point(182) = {0.329423,0.000000,-0.003184};
+Point(183) = {0.332532,0.000000,-0.003664};
+Point(184) = {0.335491,0.000000,-0.004121};
+Point(185) = {0.338298,0.000000,-0.004554};
+Point(186) = {0.340948,0.000000,-0.004963};
+Point(187) = {0.343440,0.000000,-0.005348};
+Point(188) = {0.345771,0.000000,-0.005706};
+Point(189) = {0.347939,0.000000,-0.006040};
+Point(190) = {0.349941,0.000000,-0.006347};
+Point(191) = {0.351775,0.000000,-0.006627};
+Point(192) = {0.353441,0.000000,-0.006881};
+Point(193) = {0.354935,0.000000,-0.007144};
+Point(194) = {0.356257,0.000000,-0.007382};
+Point(195) = {0.357405,0.000000,-0.007593};
+Point(196) = {0.358380,0.000000,-0.007776};
+Point(197) = {0.359180,0.000000,-0.007928};
+Point(198) = {0.359802,0.000000,-0.008049};
+Point(199) = {0.360247,0.000000,-0.008135};
+Point(200) = {0.360513,0.000000,-0.008186};
+// -- Airfoil 1
+Point(501) = {0.421732,0.200000,-0.005789};
+Point(502) = {0.421656,0.200000,-0.005794};
+Point(503) = {0.421422,0.200000,-0.005780};
+Point(504) = {0.421030,0.200000,-0.005751};
+Point(505) = {0.420481,0.200000,-0.005708};
+Point(506) = {0.419777,0.200000,-0.005654};
+Point(507) = {0.418918,0.200000,-0.005590};
+Point(508) = {0.417906,0.200000,-0.005519};
+Point(509) = {0.416742,0.200000,-0.005442};
+Point(510) = {0.415426,0.200000,-0.005362};
+Point(511) = {0.413959,0.200000,-0.005274};
+Point(512) = {0.412343,0.200000,-0.005196};
+Point(513) = {0.410579,0.200000,-0.005130};
+Point(514) = {0.408670,0.200000,-0.005077};
+Point(515) = {0.406616,0.200000,-0.005039};
+Point(516) = {0.404421,0.200000,-0.005016};
+Point(517) = {0.402087,0.200000,-0.005011};
+Point(518) = {0.399614,0.200000,-0.005031};
+Point(519) = {0.397007,0.200000,-0.005081};
+Point(520) = {0.394268,0.200000,-0.005168};
+Point(521) = {0.391399,0.200000,-0.005299};
+Point(522) = {0.388404,0.200000,-0.005479};
+Point(523) = {0.385284,0.200000,-0.005707};
+Point(524) = {0.382044,0.200000,-0.005980};
+Point(525) = {0.378686,0.200000,-0.006296};
+Point(526) = {0.375213,0.200000,-0.006654};
+Point(527) = {0.371630,0.200000,-0.007050};
+Point(528) = {0.367940,0.200000,-0.007484};
+Point(529) = {0.364146,0.200000,-0.007956};
+Point(530) = {0.360252,0.200000,-0.008466};
+Point(531) = {0.356262,0.200000,-0.009012};
+Point(532) = {0.352180,0.200000,-0.009595};
+Point(533) = {0.348010,0.200000,-0.010208};
+Point(534) = {0.343756,0.200000,-0.010848};
+Point(535) = {0.339422,0.200000,-0.011508};
+Point(536) = {0.335012,0.200000,-0.012183};
+Point(537) = {0.330532,0.200000,-0.012868};
+Point(538) = {0.325984,0.200000,-0.013558};
+Point(539) = {0.321375,0.200000,-0.014249};
+Point(540) = {0.316707,0.200000,-0.014936};
+Point(541) = {0.311987,0.200000,-0.015614};
+Point(542) = {0.307218,0.200000,-0.016281};
+Point(543) = {0.302406,0.200000,-0.016927};
+Point(544) = {0.297554,0.200000,-0.017545};
+Point(545) = {0.292668,0.200000,-0.018130};
+Point(546) = {0.287753,0.200000,-0.018671};
+Point(547) = {0.282814,0.200000,-0.019165};
+Point(548) = {0.277854,0.200000,-0.019609};
+Point(549) = {0.272880,0.200000,-0.020000};
+Point(550) = {0.267896,0.200000,-0.020339};
+Point(551) = {0.262907,0.200000,-0.020623};
+Point(552) = {0.257919,0.200000,-0.020852};
+Point(553) = {0.252935,0.200000,-0.021028};
+Point(554) = {0.247961,0.200000,-0.021151};
+Point(555) = {0.243001,0.200000,-0.021224};
+Point(556) = {0.238062,0.200000,-0.021248};
+Point(557) = {0.233146,0.200000,-0.021224};
+Point(558) = {0.228261,0.200000,-0.021154};
+Point(559) = {0.223409,0.200000,-0.021038};
+Point(560) = {0.218597,0.200000,-0.020878};
+Point(561) = {0.213828,0.200000,-0.020673};
+Point(562) = {0.209107,0.200000,-0.020426};
+Point(563) = {0.204440,0.200000,-0.020136};
+Point(564) = {0.199830,0.200000,-0.019802};
+Point(565) = {0.195283,0.200000,-0.019425};
+Point(566) = {0.190802,0.200000,-0.019004};
+Point(567) = {0.186393,0.200000,-0.018540};
+Point(568) = {0.182059,0.200000,-0.018035};
+Point(569) = {0.177805,0.200000,-0.017491};
+Point(570) = {0.173634,0.200000,-0.016910};
+Point(571) = {0.169552,0.200000,-0.016296};
+Point(572) = {0.165563,0.200000,-0.015651};
+Point(573) = {0.161668,0.200000,-0.014979};
+Point(574) = {0.157875,0.200000,-0.014287};
+Point(575) = {0.154184,0.200000,-0.013579};
+Point(576) = {0.150601,0.200000,-0.012859};
+Point(577) = {0.147129,0.200000,-0.012133};
+Point(578) = {0.143771,0.200000,-0.011402};
+Point(579) = {0.140531,0.200000,-0.010670};
+Point(580) = {0.137411,0.200000,-0.009938};
+Point(581) = {0.134415,0.200000,-0.009208};
+Point(582) = {0.131546,0.200000,-0.008482};
+Point(583) = {0.128807,0.200000,-0.007761};
+Point(584) = {0.126200,0.200000,-0.007045};
+Point(585) = {0.123728,0.200000,-0.006334};
+Point(586) = {0.121393,0.200000,-0.005629};
+Point(587) = {0.119198,0.200000,-0.004931};
+Point(588) = {0.117145,0.200000,-0.004240};
+Point(589) = {0.115236,0.200000,-0.003557};
+Point(590) = {0.113472,0.200000,-0.002883};
+Point(591) = {0.111856,0.200000,-0.002219};
+Point(592) = {0.110389,0.200000,-0.001564};
+Point(593) = {0.109072,0.200000,-0.000908};
+Point(594) = {0.107907,0.200000,-0.000236};
+Point(595) = {0.106896,0.200000,0.000462};
+Point(596) = {0.106038,0.200000,0.001200};
+Point(597) = {0.105335,0.200000,0.001986};
+Point(598) = {0.104788,0.200000,0.002812};
+Point(599) = {0.104396,0.200000,0.003664};
+Point(600) = {0.104161,0.200000,0.004530};
+Point(601) = {0.104082,0.200000,0.005397};
+Point(602) = {0.104161,0.200000,0.006254};
+Point(603) = {0.104396,0.200000,0.007099};
+Point(604) = {0.104788,0.200000,0.007933};
+Point(605) = {0.105335,0.200000,0.008754};
+Point(606) = {0.106038,0.200000,0.009563};
+Point(607) = {0.106896,0.200000,0.010359};
+Point(608) = {0.107907,0.200000,0.011133};
+Point(609) = {0.109072,0.200000,0.011877};
+Point(610) = {0.110389,0.200000,0.012581};
+Point(611) = {0.111856,0.200000,0.013236};
+Point(612) = {0.113472,0.200000,0.013835};
+Point(613) = {0.115236,0.200000,0.014380};
+Point(614) = {0.117145,0.200000,0.014875};
+Point(615) = {0.119198,0.200000,0.015322};
+Point(616) = {0.121393,0.200000,0.015728};
+Point(617) = {0.123728,0.200000,0.016095};
+Point(618) = {0.126200,0.200000,0.016425};
+Point(619) = {0.128807,0.200000,0.016722};
+Point(620) = {0.131546,0.200000,0.016987};
+Point(621) = {0.134415,0.200000,0.017222};
+Point(622) = {0.137411,0.200000,0.017431};
+Point(623) = {0.140531,0.200000,0.017614};
+Point(624) = {0.143771,0.200000,0.017773};
+Point(625) = {0.147129,0.200000,0.017908};
+Point(626) = {0.150601,0.200000,0.018023};
+Point(627) = {0.154184,0.200000,0.018116};
+Point(628) = {0.157875,0.200000,0.018191};
+Point(629) = {0.161668,0.200000,0.018247};
+Point(630) = {0.165563,0.200000,0.018286};
+Point(631) = {0.169552,0.200000,0.018309};
+Point(632) = {0.173634,0.200000,0.018317};
+Point(633) = {0.177805,0.200000,0.018309};
+Point(634) = {0.182059,0.200000,0.018283};
+Point(635) = {0.186393,0.200000,0.018238};
+Point(636) = {0.190802,0.200000,0.018173};
+Point(637) = {0.195283,0.200000,0.018086};
+Point(638) = {0.199830,0.200000,0.017978};
+Point(639) = {0.204440,0.200000,0.017850};
+Point(640) = {0.209107,0.200000,0.017703};
+Point(641) = {0.213828,0.200000,0.017539};
+Point(642) = {0.218597,0.200000,0.017357};
+Point(643) = {0.223409,0.200000,0.017158};
+Point(644) = {0.228261,0.200000,0.016940};
+Point(645) = {0.233146,0.200000,0.016704};
+Point(646) = {0.238062,0.200000,0.016448};
+Point(647) = {0.243001,0.200000,0.016172};
+Point(648) = {0.247961,0.200000,0.015875};
+Point(649) = {0.252935,0.200000,0.015558};
+Point(650) = {0.257919,0.200000,0.015219};
+Point(651) = {0.262907,0.200000,0.014858};
+Point(652) = {0.267896,0.200000,0.014475};
+Point(653) = {0.272880,0.200000,0.014072};
+Point(654) = {0.277854,0.200000,0.013647};
+Point(655) = {0.282814,0.200000,0.013203};
+Point(656) = {0.287753,0.200000,0.012740};
+Point(657) = {0.292668,0.200000,0.012258};
+Point(658) = {0.297554,0.200000,0.011760};
+Point(659) = {0.302406,0.200000,0.011244};
+Point(660) = {0.307218,0.200000,0.010713};
+Point(661) = {0.311987,0.200000,0.010168};
+Point(662) = {0.316707,0.200000,0.009609};
+Point(663) = {0.321375,0.200000,0.009038};
+Point(664) = {0.325984,0.200000,0.008459};
+Point(665) = {0.330532,0.200000,0.007875};
+Point(666) = {0.335012,0.200000,0.007287};
+Point(667) = {0.339422,0.200000,0.006698};
+Point(668) = {0.343756,0.200000,0.006109};
+Point(669) = {0.348010,0.200000,0.005520};
+Point(670) = {0.352180,0.200000,0.004934};
+Point(671) = {0.356262,0.200000,0.004349};
+Point(672) = {0.360252,0.200000,0.003767};
+Point(673) = {0.364146,0.200000,0.003190};
+Point(674) = {0.367940,0.200000,0.002621};
+Point(675) = {0.371630,0.200000,0.002060};
+Point(676) = {0.375213,0.200000,0.001510};
+Point(677) = {0.378686,0.200000,0.000974};
+Point(678) = {0.382044,0.200000,0.000451};
+Point(679) = {0.385284,0.200000,-0.000057};
+Point(680) = {0.388404,0.200000,-0.000549};
+Point(681) = {0.391399,0.200000,-0.001024};
+Point(682) = {0.394268,0.200000,-0.001482};
+Point(683) = {0.397007,0.200000,-0.001921};
+Point(684) = {0.399614,0.200000,-0.002337};
+Point(685) = {0.402087,0.200000,-0.002730};
+Point(686) = {0.404421,0.200000,-0.003097};
+Point(687) = {0.406616,0.200000,-0.003437};
+Point(688) = {0.408670,0.200000,-0.003750};
+Point(689) = {0.410579,0.200000,-0.004037};
+Point(690) = {0.412343,0.200000,-0.004299};
+Point(691) = {0.413959,0.200000,-0.004538};
+Point(692) = {0.415426,0.200000,-0.004753};
+Point(693) = {0.416742,0.200000,-0.004960};
+Point(694) = {0.417906,0.200000,-0.005149};
+Point(695) = {0.418918,0.200000,-0.005318};
+Point(696) = {0.419777,0.200000,-0.005465};
+Point(697) = {0.420481,0.200000,-0.005587};
+Point(698) = {0.421030,0.200000,-0.005683};
+Point(699) = {0.421422,0.200000,-0.005750};
+Point(700) = {0.421656,0.200000,-0.005786};
+// -- Airfoil 2
+Point(1001) = {0.459844,0.325000,-0.004603};
+Point(1002) = {0.459774,0.325000,-0.004604};
+Point(1003) = {0.459560,0.325000,-0.004590};
+Point(1004) = {0.459201,0.325000,-0.004562};
+Point(1005) = {0.458699,0.325000,-0.004524};
+Point(1006) = {0.458054,0.325000,-0.004477};
+Point(1007) = {0.457268,0.325000,-0.004423};
+Point(1008) = {0.456342,0.325000,-0.004363};
+Point(1009) = {0.455276,0.325000,-0.004300};
+Point(1010) = {0.454072,0.325000,-0.004237};
+Point(1011) = {0.452730,0.325000,-0.004154};
+Point(1012) = {0.451251,0.325000,-0.004080};
+Point(1013) = {0.449636,0.325000,-0.004017};
+Point(1014) = {0.447889,0.325000,-0.003968};
+Point(1015) = {0.446010,0.325000,-0.003933};
+Point(1016) = {0.444001,0.325000,-0.003914};
+Point(1017) = {0.441864,0.325000,-0.003915};
+Point(1018) = {0.439602,0.325000,-0.003938};
+Point(1019) = {0.437216,0.325000,-0.003988};
+Point(1020) = {0.434709,0.325000,-0.004070};
+Point(1021) = {0.432083,0.325000,-0.004188};
+Point(1022) = {0.429342,0.325000,-0.004346};
+Point(1023) = {0.426487,0.325000,-0.004544};
+Point(1024) = {0.423521,0.325000,-0.004782};
+Point(1025) = {0.420448,0.325000,-0.005061};
+Point(1026) = {0.417270,0.325000,-0.005381};
+Point(1027) = {0.413991,0.325000,-0.005741};
+Point(1028) = {0.410614,0.325000,-0.006141};
+Point(1029) = {0.407141,0.325000,-0.006578};
+Point(1030) = {0.403578,0.325000,-0.007052};
+Point(1031) = {0.399926,0.325000,-0.007560};
+Point(1032) = {0.396191,0.325000,-0.008102};
+Point(1033) = {0.392374,0.325000,-0.008672};
+Point(1034) = {0.388480,0.325000,-0.009267};
+Point(1035) = {0.384514,0.325000,-0.009882};
+Point(1036) = {0.380479,0.325000,-0.010511};
+Point(1037) = {0.376378,0.325000,-0.011152};
+Point(1038) = {0.372216,0.325000,-0.011798};
+Point(1039) = {0.367997,0.325000,-0.012447};
+Point(1040) = {0.363726,0.325000,-0.013094};
+Point(1041) = {0.359406,0.325000,-0.013735};
+Point(1042) = {0.355041,0.325000,-0.014364};
+Point(1043) = {0.350637,0.325000,-0.014977};
+Point(1044) = {0.346197,0.325000,-0.015565};
+Point(1045) = {0.341726,0.325000,-0.016122};
+Point(1046) = {0.337227,0.325000,-0.016641};
+Point(1047) = {0.332707,0.325000,-0.017117};
+Point(1048) = {0.328168,0.325000,-0.017547};
+Point(1049) = {0.323616,0.325000,-0.017930};
+Point(1050) = {0.319054,0.325000,-0.018265};
+Point(1051) = {0.314489,0.325000,-0.018550};
+Point(1052) = {0.309923,0.325000,-0.018784};
+Point(1053) = {0.305362,0.325000,-0.018969};
+Point(1054) = {0.300810,0.325000,-0.019105};
+Point(1055) = {0.296271,0.325000,-0.019195};
+Point(1056) = {0.291750,0.325000,-0.019239};
+Point(1057) = {0.287252,0.325000,-0.019239};
+Point(1058) = {0.282781,0.325000,-0.019196};
+Point(1059) = {0.278340,0.325000,-0.019110};
+Point(1060) = {0.273936,0.325000,-0.018984};
+Point(1061) = {0.269572,0.325000,-0.018818};
+Point(1062) = {0.265252,0.325000,-0.018612};
+Point(1063) = {0.260980,0.325000,-0.018369};
+Point(1064) = {0.256761,0.325000,-0.018086};
+Point(1065) = {0.252600,0.325000,-0.017765};
+Point(1066) = {0.248499,0.325000,-0.017406};
+Point(1067) = {0.244463,0.325000,-0.017008};
+Point(1068) = {0.240497,0.325000,-0.016575};
+Point(1069) = {0.236604,0.325000,-0.016108};
+Point(1070) = {0.232787,0.325000,-0.015609};
+Point(1071) = {0.229051,0.325000,-0.015081};
+Point(1072) = {0.225400,0.325000,-0.014526};
+Point(1073) = {0.221836,0.325000,-0.013949};
+Point(1074) = {0.218364,0.325000,-0.013352};
+Point(1075) = {0.214987,0.325000,-0.012740};
+Point(1076) = {0.211707,0.325000,-0.012117};
+Point(1077) = {0.208530,0.325000,-0.011486};
+Point(1078) = {0.205457,0.325000,-0.010850};
+Point(1079) = {0.202491,0.325000,-0.010211};
+Point(1080) = {0.199636,0.325000,-0.009572};
+Point(1081) = {0.196894,0.325000,-0.008934};
+Point(1082) = {0.194269,0.325000,-0.008300};
+Point(1083) = {0.191762,0.325000,-0.007671};
+Point(1084) = {0.189376,0.325000,-0.007048};
+Point(1085) = {0.187113,0.325000,-0.006434};
+Point(1086) = {0.184977,0.325000,-0.005828};
+Point(1087) = {0.182968,0.325000,-0.005232};
+Point(1088) = {0.181089,0.325000,-0.004640};
+Point(1089) = {0.179341,0.325000,-0.004048};
+Point(1090) = {0.177727,0.325000,-0.003449};
+Point(1091) = {0.176248,0.325000,-0.002838};
+Point(1092) = {0.174905,0.325000,-0.002212};
+Point(1093) = {0.173701,0.325000,-0.001579};
+Point(1094) = {0.172635,0.325000,-0.000948};
+Point(1095) = {0.171709,0.325000,-0.000330};
+Point(1096) = {0.170923,0.325000,0.000266};
+Point(1097) = {0.170280,0.325000,0.000836};
+Point(1098) = {0.169779,0.325000,0.001404};
+Point(1099) = {0.169421,0.325000,0.002000};
+Point(1100) = {0.169206,0.325000,0.002653};
+Point(1101) = {0.169134,0.325000,0.003395};
+Point(1102) = {0.169206,0.325000,0.004244};
+Point(1103) = {0.169421,0.325000,0.005168};
+Point(1104) = {0.169779,0.325000,0.006124};
+Point(1105) = {0.170280,0.325000,0.007069};
+Point(1106) = {0.170923,0.325000,0.007959};
+Point(1107) = {0.171709,0.325000,0.008761};
+Point(1108) = {0.172635,0.325000,0.009481};
+Point(1109) = {0.173701,0.325000,0.010135};
+Point(1110) = {0.174905,0.325000,0.010738};
+Point(1111) = {0.176248,0.325000,0.011306};
+Point(1112) = {0.177727,0.325000,0.011853};
+Point(1113) = {0.179341,0.325000,0.012374};
+Point(1114) = {0.181089,0.325000,0.012866};
+Point(1115) = {0.182968,0.325000,0.013323};
+Point(1116) = {0.184977,0.325000,0.013738};
+Point(1117) = {0.187113,0.325000,0.014110};
+Point(1118) = {0.189376,0.325000,0.014440};
+Point(1119) = {0.191762,0.325000,0.014735};
+Point(1120) = {0.194269,0.325000,0.015000};
+Point(1121) = {0.196894,0.325000,0.015239};
+Point(1122) = {0.199636,0.325000,0.015458};
+Point(1123) = {0.202491,0.325000,0.015658};
+Point(1124) = {0.205457,0.325000,0.015840};
+Point(1125) = {0.208530,0.325000,0.016003};
+Point(1126) = {0.211707,0.325000,0.016150};
+Point(1127) = {0.214987,0.325000,0.016279};
+Point(1128) = {0.218364,0.325000,0.016391};
+Point(1129) = {0.221836,0.325000,0.016487};
+Point(1130) = {0.225400,0.325000,0.016566};
+Point(1131) = {0.229051,0.325000,0.016627};
+Point(1132) = {0.232787,0.325000,0.016670};
+Point(1133) = {0.236604,0.325000,0.016695};
+Point(1134) = {0.240497,0.325000,0.016701};
+Point(1135) = {0.244463,0.325000,0.016688};
+Point(1136) = {0.248499,0.325000,0.016655};
+Point(1137) = {0.252600,0.325000,0.016601};
+Point(1138) = {0.256761,0.325000,0.016528};
+Point(1139) = {0.260980,0.325000,0.016437};
+Point(1140) = {0.265252,0.325000,0.016329};
+Point(1141) = {0.269572,0.325000,0.016207};
+Point(1142) = {0.273936,0.325000,0.016070};
+Point(1143) = {0.278340,0.325000,0.015918};
+Point(1144) = {0.282781,0.325000,0.015751};
+Point(1145) = {0.287252,0.325000,0.015567};
+Point(1146) = {0.291750,0.325000,0.015365};
+Point(1147) = {0.296271,0.325000,0.015144};
+Point(1148) = {0.300810,0.325000,0.014905};
+Point(1149) = {0.305362,0.325000,0.014645};
+Point(1150) = {0.309923,0.325000,0.014366};
+Point(1151) = {0.314489,0.325000,0.014066};
+Point(1152) = {0.319054,0.325000,0.013746};
+Point(1153) = {0.323616,0.325000,0.013405};
+Point(1154) = {0.328168,0.325000,0.013044};
+Point(1155) = {0.332707,0.325000,0.012666};
+Point(1156) = {0.337227,0.325000,0.012269};
+Point(1157) = {0.341726,0.325000,0.011856};
+Point(1158) = {0.346197,0.325000,0.011426};
+Point(1159) = {0.350637,0.325000,0.010981};
+Point(1160) = {0.355041,0.325000,0.010521};
+Point(1161) = {0.359406,0.325000,0.010047};
+Point(1162) = {0.363726,0.325000,0.009560};
+Point(1163) = {0.367997,0.325000,0.009061};
+Point(1164) = {0.372216,0.325000,0.008554};
+Point(1165) = {0.376378,0.325000,0.008038};
+Point(1166) = {0.380479,0.325000,0.007518};
+Point(1167) = {0.384514,0.325000,0.006993};
+Point(1168) = {0.388480,0.325000,0.006465};
+Point(1169) = {0.392374,0.325000,0.005936};
+Point(1170) = {0.396191,0.325000,0.005405};
+Point(1171) = {0.399926,0.325000,0.004873};
+Point(1172) = {0.403578,0.325000,0.004341};
+Point(1173) = {0.407141,0.325000,0.003812};
+Point(1174) = {0.410614,0.325000,0.003288};
+Point(1175) = {0.413991,0.325000,0.002771};
+Point(1176) = {0.417270,0.325000,0.002264};
+Point(1177) = {0.420448,0.325000,0.001769};
+Point(1178) = {0.423521,0.325000,0.001286};
+Point(1179) = {0.426487,0.325000,0.000817};
+Point(1180) = {0.429342,0.325000,0.000360};
+Point(1181) = {0.432083,0.325000,-0.000084};
+Point(1182) = {0.434709,0.325000,-0.000514};
+Point(1183) = {0.437216,0.325000,-0.000928};
+Point(1184) = {0.439602,0.325000,-0.001322};
+Point(1185) = {0.441864,0.325000,-0.001693};
+Point(1186) = {0.444001,0.325000,-0.002038};
+Point(1187) = {0.446010,0.325000,-0.002354};
+Point(1188) = {0.447889,0.325000,-0.002644};
+Point(1189) = {0.449636,0.325000,-0.002908};
+Point(1190) = {0.451251,0.325000,-0.003149};
+Point(1191) = {0.452730,0.325000,-0.003368};
+Point(1192) = {0.454072,0.325000,-0.003567};
+Point(1193) = {0.455276,0.325000,-0.003771};
+Point(1194) = {0.456342,0.325000,-0.003957};
+Point(1195) = {0.457268,0.325000,-0.004124};
+Point(1196) = {0.458054,0.325000,-0.004269};
+Point(1197) = {0.458699,0.325000,-0.004391};
+Point(1198) = {0.459201,0.325000,-0.004488};
+Point(1199) = {0.459560,0.325000,-0.004557};
+Point(1200) = {0.459774,0.325000,-0.004596};
+// -- Airfoil 3
+Point(1501) = {0.505256,0.475000,-0.003178};
+Point(1502) = {0.505194,0.475000,-0.003180};
+Point(1503) = {0.505003,0.475000,-0.003169};
+Point(1504) = {0.504685,0.475000,-0.003147};
+Point(1505) = {0.504239,0.475000,-0.003114};
+Point(1506) = {0.503667,0.475000,-0.003073};
+Point(1507) = {0.502969,0.475000,-0.003025};
+Point(1508) = {0.502147,0.475000,-0.002972};
+Point(1509) = {0.501201,0.475000,-0.002916};
+Point(1510) = {0.500132,0.475000,-0.002858};
+Point(1511) = {0.498940,0.475000,-0.002777};
+Point(1512) = {0.497628,0.475000,-0.002701};
+Point(1513) = {0.496195,0.475000,-0.002634};
+Point(1514) = {0.494643,0.475000,-0.002578};
+Point(1515) = {0.492975,0.475000,-0.002534};
+Point(1516) = {0.491192,0.475000,-0.002504};
+Point(1517) = {0.489295,0.475000,-0.002492};
+Point(1518) = {0.487287,0.475000,-0.002501};
+Point(1519) = {0.485169,0.475000,-0.002535};
+Point(1520) = {0.482944,0.475000,-0.002599};
+Point(1521) = {0.480613,0.475000,-0.002697};
+Point(1522) = {0.478179,0.475000,-0.002833};
+Point(1523) = {0.475645,0.475000,-0.003007};
+Point(1524) = {0.473012,0.475000,-0.003219};
+Point(1525) = {0.470284,0.475000,-0.003468};
+Point(1526) = {0.467463,0.475000,-0.003753};
+Point(1527) = {0.464553,0.475000,-0.004076};
+Point(1528) = {0.461554,0.475000,-0.004435};
+Point(1529) = {0.458472,0.475000,-0.004830};
+Point(1530) = {0.455309,0.475000,-0.005261};
+Point(1531) = {0.452068,0.475000,-0.005728};
+Point(1532) = {0.448751,0.475000,-0.006231};
+Point(1533) = {0.445363,0.475000,-0.006763};
+Point(1534) = {0.441907,0.475000,-0.007318};
+Point(1535) = {0.438386,0.475000,-0.007888};
+Point(1536) = {0.434804,0.475000,-0.008469};
+Point(1537) = {0.431164,0.475000,-0.009054};
+Point(1538) = {0.427470,0.475000,-0.009641};
+Point(1539) = {0.423725,0.475000,-0.010229};
+Point(1540) = {0.419933,0.475000,-0.010817};
+Point(1541) = {0.416098,0.475000,-0.011404};
+Point(1542) = {0.412224,0.475000,-0.011987};
+Point(1543) = {0.408314,0.475000,-0.012561};
+Point(1544) = {0.404372,0.475000,-0.013118};
+Point(1545) = {0.400403,0.475000,-0.013650};
+Point(1546) = {0.396410,0.475000,-0.014149};
+Point(1547) = {0.392397,0.475000,-0.014608};
+Point(1548) = {0.388368,0.475000,-0.015026};
+Point(1549) = {0.384327,0.475000,-0.015402};
+Point(1550) = {0.380278,0.475000,-0.015734};
+Point(1551) = {0.376226,0.475000,-0.016021};
+Point(1552) = {0.372173,0.475000,-0.016262};
+Point(1553) = {0.368124,0.475000,-0.016459};
+Point(1554) = {0.364083,0.475000,-0.016612};
+Point(1555) = {0.360054,0.475000,-0.016723};
+Point(1556) = {0.356041,0.475000,-0.016793};
+Point(1557) = {0.352048,0.475000,-0.016823};
+Point(1558) = {0.348079,0.475000,-0.016815};
+Point(1559) = {0.344137,0.475000,-0.016770};
+Point(1560) = {0.340228,0.475000,-0.016688};
+Point(1561) = {0.336353,0.475000,-0.016573};
+Point(1562) = {0.332518,0.475000,-0.016424};
+Point(1563) = {0.328727,0.475000,-0.016242};
+Point(1564) = {0.324982,0.475000,-0.016027};
+Point(1565) = {0.321287,0.475000,-0.015780};
+Point(1566) = {0.317647,0.475000,-0.015501};
+Point(1567) = {0.314065,0.475000,-0.015189};
+Point(1568) = {0.310544,0.475000,-0.014847};
+Point(1569) = {0.307088,0.475000,-0.014477};
+Point(1570) = {0.303700,0.475000,-0.014079};
+Point(1571) = {0.300384,0.475000,-0.013658};
+Point(1572) = {0.297142,0.475000,-0.013213};
+Point(1573) = {0.293979,0.475000,-0.012749};
+Point(1574) = {0.290897,0.475000,-0.012268};
+Point(1575) = {0.287899,0.475000,-0.011774};
+Point(1576) = {0.284988,0.475000,-0.011269};
+Point(1577) = {0.282167,0.475000,-0.010757};
+Point(1578) = {0.279439,0.475000,-0.010241};
+Point(1579) = {0.276806,0.475000,-0.009722};
+Point(1580) = {0.274272,0.475000,-0.009203};
+Point(1581) = {0.271838,0.475000,-0.008688};
+Point(1582) = {0.269507,0.475000,-0.008177};
+Point(1583) = {0.267282,0.475000,-0.007671};
+Point(1584) = {0.265164,0.475000,-0.007166};
+Point(1585) = {0.263156,0.475000,-0.006662};
+Point(1586) = {0.261259,0.475000,-0.006157};
+Point(1587) = {0.259476,0.475000,-0.005649};
+Point(1588) = {0.257808,0.475000,-0.005141};
+Point(1589) = {0.256257,0.475000,-0.004635};
+Point(1590) = {0.254824,0.475000,-0.004133};
+Point(1591) = {0.253511,0.475000,-0.003639};
+Point(1592) = {0.252319,0.475000,-0.003151};
+Point(1593) = {0.251249,0.475000,-0.002659};
+Point(1594) = {0.250303,0.475000,-0.002145};
+Point(1595) = {0.249481,0.475000,-0.001597};
+Point(1596) = {0.248784,0.475000,-0.000998};
+Point(1597) = {0.248213,0.475000,-0.000338};
+Point(1598) = {0.247768,0.475000,0.000368};
+Point(1599) = {0.247450,0.475000,0.001103};
+Point(1600) = {0.247259,0.475000,0.001847};
+Point(1601) = {0.247196,0.475000,0.002583};
+Point(1602) = {0.247259,0.475000,0.003295};
+Point(1603) = {0.247450,0.475000,0.003985};
+Point(1604) = {0.247768,0.475000,0.004656};
+Point(1605) = {0.248213,0.475000,0.005312};
+Point(1606) = {0.248784,0.475000,0.005959};
+Point(1607) = {0.249481,0.475000,0.006598};
+Point(1608) = {0.250303,0.475000,0.007225};
+Point(1609) = {0.251249,0.475000,0.007834};
+Point(1610) = {0.252319,0.475000,0.008418};
+Point(1611) = {0.253511,0.475000,0.008972};
+Point(1612) = {0.254824,0.475000,0.009490};
+Point(1613) = {0.256257,0.475000,0.009972};
+Point(1614) = {0.257808,0.475000,0.010421};
+Point(1615) = {0.259476,0.475000,0.010839};
+Point(1616) = {0.261259,0.475000,0.011226};
+Point(1617) = {0.263156,0.475000,0.011586};
+Point(1618) = {0.265164,0.475000,0.011920};
+Point(1619) = {0.267282,0.475000,0.012229};
+Point(1620) = {0.269507,0.475000,0.012515};
+Point(1621) = {0.271838,0.475000,0.012779};
+Point(1622) = {0.274272,0.475000,0.013023};
+Point(1623) = {0.276806,0.475000,0.013248};
+Point(1624) = {0.279439,0.475000,0.013455};
+Point(1625) = {0.282167,0.475000,0.013645};
+Point(1626) = {0.284988,0.475000,0.013820};
+Point(1627) = {0.287899,0.475000,0.013979};
+Point(1628) = {0.290897,0.475000,0.014122};
+Point(1629) = {0.293979,0.475000,0.014249};
+Point(1630) = {0.297142,0.475000,0.014358};
+Point(1631) = {0.300384,0.475000,0.014450};
+Point(1632) = {0.303700,0.475000,0.014523};
+Point(1633) = {0.307088,0.475000,0.014580};
+Point(1634) = {0.310544,0.475000,0.014621};
+Point(1635) = {0.314065,0.475000,0.014649};
+Point(1636) = {0.317647,0.475000,0.014665};
+Point(1637) = {0.321287,0.475000,0.014670};
+Point(1638) = {0.324982,0.475000,0.014663};
+Point(1639) = {0.328727,0.475000,0.014644};
+Point(1640) = {0.332518,0.475000,0.014609};
+Point(1641) = {0.336353,0.475000,0.014558};
+Point(1642) = {0.340228,0.475000,0.014489};
+Point(1643) = {0.344137,0.475000,0.014402};
+Point(1644) = {0.348079,0.475000,0.014297};
+Point(1645) = {0.352048,0.475000,0.014174};
+Point(1646) = {0.356041,0.475000,0.014035};
+Point(1647) = {0.360054,0.475000,0.013878};
+Point(1648) = {0.364083,0.475000,0.013704};
+Point(1649) = {0.368124,0.475000,0.013512};
+Point(1650) = {0.372173,0.475000,0.013302};
+Point(1651) = {0.376226,0.475000,0.013073};
+Point(1652) = {0.380278,0.475000,0.012826};
+Point(1653) = {0.384327,0.475000,0.012561};
+Point(1654) = {0.388368,0.475000,0.012280};
+Point(1655) = {0.392397,0.475000,0.011983};
+Point(1656) = {0.396410,0.475000,0.011673};
+Point(1657) = {0.400403,0.475000,0.011348};
+Point(1658) = {0.404372,0.475000,0.011011};
+Point(1659) = {0.408314,0.475000,0.010659};
+Point(1660) = {0.412224,0.475000,0.010293};
+Point(1661) = {0.416098,0.475000,0.009911};
+Point(1662) = {0.419933,0.475000,0.009514};
+Point(1663) = {0.423725,0.475000,0.009104};
+Point(1664) = {0.427470,0.475000,0.008682};
+Point(1665) = {0.431164,0.475000,0.008250};
+Point(1666) = {0.434804,0.475000,0.007811};
+Point(1667) = {0.438386,0.475000,0.007366};
+Point(1668) = {0.441907,0.475000,0.006916};
+Point(1669) = {0.445363,0.475000,0.006460};
+Point(1670) = {0.448751,0.475000,0.005999};
+Point(1671) = {0.452068,0.475000,0.005532};
+Point(1672) = {0.455309,0.475000,0.005060};
+Point(1673) = {0.458472,0.475000,0.004586};
+Point(1674) = {0.461554,0.475000,0.004114};
+Point(1675) = {0.464553,0.475000,0.003646};
+Point(1676) = {0.467463,0.475000,0.003187};
+Point(1677) = {0.470284,0.475000,0.002738};
+Point(1678) = {0.473012,0.475000,0.002302};
+Point(1679) = {0.475645,0.475000,0.001878};
+Point(1680) = {0.478179,0.475000,0.001466};
+Point(1681) = {0.480613,0.475000,0.001065};
+Point(1682) = {0.482944,0.475000,0.000677};
+Point(1683) = {0.485169,0.475000,0.000303};
+Point(1684) = {0.487287,0.475000,-0.000054};
+Point(1685) = {0.489295,0.475000,-0.000393};
+Point(1686) = {0.491192,0.475000,-0.000712};
+Point(1687) = {0.492975,0.475000,-0.001008};
+Point(1688) = {0.494643,0.475000,-0.001283};
+Point(1689) = {0.496195,0.475000,-0.001536};
+Point(1690) = {0.497628,0.475000,-0.001769};
+Point(1691) = {0.498940,0.475000,-0.001982};
+Point(1692) = {0.500132,0.475000,-0.002175};
+Point(1693) = {0.501201,0.475000,-0.002375};
+Point(1694) = {0.502147,0.475000,-0.002558};
+Point(1695) = {0.502969,0.475000,-0.002720};
+Point(1696) = {0.503667,0.475000,-0.002861};
+Point(1697) = {0.504239,0.475000,-0.002978};
+Point(1698) = {0.504685,0.475000,-0.003070};
+Point(1699) = {0.505003,0.475000,-0.003136};
+Point(1700) = {0.505194,0.475000,-0.003172};
+// -- Airfoil 4
+Point(2001) = {0.558558,0.650000,-0.001008};
+Point(2002) = {0.558505,0.650000,-0.001006};
+Point(2003) = {0.558342,0.650000,-0.000993};
+Point(2004) = {0.558070,0.650000,-0.000969};
+Point(2005) = {0.557690,0.650000,-0.000937};
+Point(2006) = {0.557202,0.650000,-0.000896};
+Point(2007) = {0.556606,0.650000,-0.000850};
+Point(2008) = {0.555904,0.650000,-0.000798};
+Point(2009) = {0.555097,0.650000,-0.000742};
+Point(2010) = {0.554184,0.650000,-0.000684};
+Point(2011) = {0.553167,0.650000,-0.000602};
+Point(2012) = {0.552046,0.650000,-0.000526};
+Point(2013) = {0.550823,0.650000,-0.000459};
+Point(2014) = {0.549499,0.650000,-0.000401};
+Point(2015) = {0.548075,0.650000,-0.000356};
+Point(2016) = {0.546553,0.650000,-0.000327};
+Point(2017) = {0.544933,0.650000,-0.000315};
+Point(2018) = {0.543219,0.650000,-0.000324};
+Point(2019) = {0.541411,0.650000,-0.000356};
+Point(2020) = {0.539511,0.650000,-0.000415};
+Point(2021) = {0.537522,0.650000,-0.000505};
+Point(2022) = {0.535444,0.650000,-0.000627};
+Point(2023) = {0.533281,0.650000,-0.000781};
+Point(2024) = {0.531034,0.650000,-0.000967};
+Point(2025) = {0.528705,0.650000,-0.001185};
+Point(2026) = {0.526297,0.650000,-0.001434};
+Point(2027) = {0.523812,0.650000,-0.001713};
+Point(2028) = {0.521253,0.650000,-0.002023};
+Point(2029) = {0.518622,0.650000,-0.002365};
+Point(2030) = {0.515921,0.650000,-0.002740};
+Point(2031) = {0.513154,0.650000,-0.003149};
+Point(2032) = {0.510323,0.650000,-0.003592};
+Point(2033) = {0.507431,0.650000,-0.004064};
+Point(2034) = {0.504481,0.650000,-0.004558};
+Point(2035) = {0.501475,0.650000,-0.005067};
+Point(2036) = {0.498417,0.650000,-0.005587};
+Point(2037) = {0.495310,0.650000,-0.006111};
+Point(2038) = {0.492157,0.650000,-0.006639};
+Point(2039) = {0.488960,0.650000,-0.007170};
+Point(2040) = {0.485723,0.650000,-0.007706};
+Point(2041) = {0.482449,0.650000,-0.008246};
+Point(2042) = {0.479142,0.650000,-0.008789};
+Point(2043) = {0.475805,0.650000,-0.009328};
+Point(2044) = {0.472440,0.650000,-0.009856};
+Point(2045) = {0.469052,0.650000,-0.010365};
+Point(2046) = {0.465643,0.650000,-0.010846};
+Point(2047) = {0.462218,0.650000,-0.011293};
+Point(2048) = {0.458778,0.650000,-0.011704};
+Point(2049) = {0.455329,0.650000,-0.012077};
+Point(2050) = {0.451872,0.650000,-0.012412};
+Point(2051) = {0.448413,0.650000,-0.012708};
+Point(2052) = {0.444953,0.650000,-0.012963};
+Point(2053) = {0.441497,0.650000,-0.013178};
+Point(2054) = {0.438047,0.650000,-0.013356};
+Point(2055) = {0.434608,0.650000,-0.013495};
+Point(2056) = {0.431182,0.650000,-0.013599};
+Point(2057) = {0.427773,0.650000,-0.013668};
+Point(2058) = {0.424385,0.650000,-0.013703};
+Point(2059) = {0.421021,0.650000,-0.013707};
+Point(2060) = {0.417683,0.650000,-0.013681};
+Point(2061) = {0.414376,0.650000,-0.013628};
+Point(2062) = {0.411102,0.650000,-0.013548};
+Point(2063) = {0.407866,0.650000,-0.013443};
+Point(2064) = {0.404669,0.650000,-0.013312};
+Point(2065) = {0.401515,0.650000,-0.013155};
+Point(2066) = {0.398408,0.650000,-0.012973};
+Point(2067) = {0.395350,0.650000,-0.012765};
+Point(2068) = {0.392344,0.650000,-0.012533};
+Point(2069) = {0.389394,0.650000,-0.012279};
+Point(2070) = {0.386502,0.650000,-0.012002};
+Point(2071) = {0.383671,0.650000,-0.011705};
+Point(2072) = {0.380904,0.650000,-0.011390};
+Point(2073) = {0.378204,0.650000,-0.011057};
+Point(2074) = {0.375573,0.650000,-0.010711};
+Point(2075) = {0.373013,0.650000,-0.010353};
+Point(2076) = {0.370528,0.650000,-0.009985};
+Point(2077) = {0.368120,0.650000,-0.009610};
+Point(2078) = {0.365792,0.650000,-0.009230};
+Point(2079) = {0.363545,0.650000,-0.008849};
+Point(2080) = {0.361381,0.650000,-0.008469};
+Point(2081) = {0.359304,0.650000,-0.008092};
+Point(2082) = {0.357314,0.650000,-0.007721};
+Point(2083) = {0.355414,0.650000,-0.007352};
+Point(2084) = {0.353606,0.650000,-0.006982};
+Point(2085) = {0.351892,0.650000,-0.006609};
+Point(2086) = {0.350273,0.650000,-0.006228};
+Point(2087) = {0.348751,0.650000,-0.005838};
+Point(2088) = {0.347327,0.650000,-0.005439};
+Point(2089) = {0.346002,0.650000,-0.005031};
+Point(2090) = {0.344779,0.650000,-0.004617};
+Point(2091) = {0.343659,0.650000,-0.004197};
+Point(2092) = {0.342641,0.650000,-0.003773};
+Point(2093) = {0.341728,0.650000,-0.003343};
+Point(2094) = {0.340920,0.650000,-0.002904};
+Point(2095) = {0.340219,0.650000,-0.002454};
+Point(2096) = {0.339624,0.650000,-0.001991};
+Point(2097) = {0.339136,0.650000,-0.001512};
+Point(2098) = {0.338757,0.650000,-0.001015};
+Point(2099) = {0.338485,0.650000,-0.000499};
+Point(2100) = {0.338322,0.650000,0.000039};
+Point(2101) = {0.338268,0.650000,0.000601};
+Point(2102) = {0.338322,0.650000,0.001187};
+Point(2103) = {0.338485,0.650000,0.001791};
+Point(2104) = {0.338757,0.650000,0.002406};
+Point(2105) = {0.339136,0.650000,0.003023};
+Point(2106) = {0.339624,0.650000,0.003634};
+Point(2107) = {0.340219,0.650000,0.004234};
+Point(2108) = {0.340920,0.650000,0.004817};
+Point(2109) = {0.341728,0.650000,0.005380};
+Point(2110) = {0.342641,0.650000,0.005919};
+Point(2111) = {0.343659,0.650000,0.006431};
+Point(2112) = {0.344779,0.650000,0.006913};
+Point(2113) = {0.346002,0.650000,0.007366};
+Point(2114) = {0.347327,0.650000,0.007792};
+Point(2115) = {0.348751,0.650000,0.008192};
+Point(2116) = {0.350273,0.650000,0.008569};
+Point(2117) = {0.351892,0.650000,0.008924};
+Point(2118) = {0.353606,0.650000,0.009259};
+Point(2119) = {0.355414,0.650000,0.009576};
+Point(2120) = {0.357314,0.650000,0.009876};
+Point(2121) = {0.359304,0.650000,0.010161};
+Point(2122) = {0.361381,0.650000,0.010433};
+Point(2123) = {0.363545,0.650000,0.010692};
+Point(2124) = {0.365792,0.650000,0.010936};
+Point(2125) = {0.368120,0.650000,0.011165};
+Point(2126) = {0.370528,0.650000,0.011379};
+Point(2127) = {0.373013,0.650000,0.011577};
+Point(2128) = {0.375573,0.650000,0.011760};
+Point(2129) = {0.378204,0.650000,0.011927};
+Point(2130) = {0.380904,0.650000,0.012081};
+Point(2131) = {0.383671,0.650000,0.012222};
+Point(2132) = {0.386502,0.650000,0.012350};
+Point(2133) = {0.389394,0.650000,0.012465};
+Point(2134) = {0.392344,0.650000,0.012566};
+Point(2135) = {0.395350,0.650000,0.012653};
+Point(2136) = {0.398408,0.650000,0.012725};
+Point(2137) = {0.401515,0.650000,0.012783};
+Point(2138) = {0.404669,0.650000,0.012826};
+Point(2139) = {0.407866,0.650000,0.012856};
+Point(2140) = {0.411102,0.650000,0.012873};
+Point(2141) = {0.414376,0.650000,0.012878};
+Point(2142) = {0.417683,0.650000,0.012871};
+Point(2143) = {0.421021,0.650000,0.012852};
+Point(2144) = {0.424385,0.650000,0.012821};
+Point(2145) = {0.427773,0.650000,0.012775};
+Point(2146) = {0.431182,0.650000,0.012715};
+Point(2147) = {0.434608,0.650000,0.012639};
+Point(2148) = {0.438047,0.650000,0.012547};
+Point(2149) = {0.441497,0.650000,0.012439};
+Point(2150) = {0.444953,0.650000,0.012315};
+Point(2151) = {0.448413,0.650000,0.012175};
+Point(2152) = {0.451872,0.650000,0.012018};
+Point(2153) = {0.455329,0.650000,0.011846};
+Point(2154) = {0.458778,0.650000,0.011659};
+Point(2155) = {0.462218,0.650000,0.011457};
+Point(2156) = {0.465643,0.650000,0.011243};
+Point(2157) = {0.469052,0.650000,0.011015};
+Point(2158) = {0.472440,0.650000,0.010775};
+Point(2159) = {0.475805,0.650000,0.010524};
+Point(2160) = {0.479142,0.650000,0.010260};
+Point(2161) = {0.482449,0.650000,0.009985};
+Point(2162) = {0.485723,0.650000,0.009699};
+Point(2163) = {0.488960,0.650000,0.009400};
+Point(2164) = {0.492157,0.650000,0.009087};
+Point(2165) = {0.495310,0.650000,0.008758};
+Point(2166) = {0.498417,0.650000,0.008413};
+Point(2167) = {0.501475,0.650000,0.008050};
+Point(2168) = {0.504481,0.650000,0.007673};
+Point(2169) = {0.507431,0.650000,0.007284};
+Point(2170) = {0.510323,0.650000,0.006888};
+Point(2171) = {0.513154,0.650000,0.006487};
+Point(2172) = {0.515921,0.650000,0.006084};
+Point(2173) = {0.518622,0.650000,0.005682};
+Point(2174) = {0.521253,0.650000,0.005281};
+Point(2175) = {0.523812,0.650000,0.004882};
+Point(2176) = {0.526297,0.650000,0.004487};
+Point(2177) = {0.528705,0.650000,0.004097};
+Point(2178) = {0.531034,0.650000,0.003713};
+Point(2179) = {0.533281,0.650000,0.003341};
+Point(2180) = {0.535444,0.650000,0.002981};
+Point(2181) = {0.537522,0.650000,0.002637};
+Point(2182) = {0.539511,0.650000,0.002311};
+Point(2183) = {0.541411,0.650000,0.002003};
+Point(2184) = {0.543219,0.650000,0.001710};
+Point(2185) = {0.544933,0.650000,0.001432};
+Point(2186) = {0.546553,0.650000,0.001168};
+Point(2187) = {0.548075,0.650000,0.000916};
+Point(2188) = {0.549499,0.650000,0.000678};
+Point(2189) = {0.550823,0.650000,0.000455};
+Point(2190) = {0.552046,0.650000,0.000248};
+Point(2191) = {0.553167,0.650000,0.000059};
+Point(2192) = {0.554184,0.650000,-0.000111};
+Point(2193) = {0.555097,0.650000,-0.000289};
+Point(2194) = {0.555904,0.650000,-0.000450};
+Point(2195) = {0.556606,0.650000,-0.000594};
+Point(2196) = {0.557202,0.650000,-0.000719};
+Point(2197) = {0.557690,0.650000,-0.000823};
+Point(2198) = {0.558070,0.650000,-0.000905};
+Point(2199) = {0.558342,0.650000,-0.000965};
+Point(2200) = {0.558505,0.650000,-0.000999};
+// -- Airfoil 5
+Point(2501) = {0.611690,0.825000,0.001248};
+Point(2502) = {0.611646,0.825000,0.001246};
+Point(2503) = {0.611511,0.825000,0.001251};
+Point(2504) = {0.611286,0.825000,0.001262};
+Point(2505) = {0.610971,0.825000,0.001279};
+Point(2506) = {0.610567,0.825000,0.001300};
+Point(2507) = {0.610074,0.825000,0.001324};
+Point(2508) = {0.609493,0.825000,0.001350};
+Point(2509) = {0.608825,0.825000,0.001376};
+Point(2510) = {0.608069,0.825000,0.001402};
+Point(2511) = {0.607227,0.825000,0.001459};
+Point(2512) = {0.606300,0.825000,0.001512};
+Point(2513) = {0.605287,0.825000,0.001558};
+Point(2514) = {0.604191,0.825000,0.001596};
+Point(2515) = {0.603012,0.825000,0.001625};
+Point(2516) = {0.601752,0.825000,0.001643};
+Point(2517) = {0.600412,0.825000,0.001648};
+Point(2518) = {0.598993,0.825000,0.001638};
+Point(2519) = {0.597496,0.825000,0.001611};
+Point(2520) = {0.595924,0.825000,0.001563};
+Point(2521) = {0.594277,0.825000,0.001491};
+Point(2522) = {0.592557,0.825000,0.001395};
+Point(2523) = {0.590766,0.825000,0.001272};
+Point(2524) = {0.588906,0.825000,0.001122};
+Point(2525) = {0.586978,0.825000,0.000942};
+Point(2526) = {0.584985,0.825000,0.000732};
+Point(2527) = {0.582928,0.825000,0.000492};
+Point(2528) = {0.580810,0.825000,0.000221};
+Point(2529) = {0.578632,0.825000,-0.000078};
+Point(2530) = {0.576396,0.825000,-0.000405};
+Point(2531) = {0.574106,0.825000,-0.000757};
+Point(2532) = {0.571763,0.825000,-0.001133};
+Point(2533) = {0.569369,0.825000,-0.001531};
+Point(2534) = {0.566926,0.825000,-0.001950};
+Point(2535) = {0.564439,0.825000,-0.002386};
+Point(2536) = {0.561907,0.825000,-0.002837};
+Point(2537) = {0.559335,0.825000,-0.003302};
+Point(2538) = {0.556725,0.825000,-0.003778};
+Point(2539) = {0.554078,0.825000,-0.004263};
+Point(2540) = {0.551399,0.825000,-0.004755};
+Point(2541) = {0.548689,0.825000,-0.005253};
+Point(2542) = {0.545952,0.825000,-0.005753};
+Point(2543) = {0.543189,0.825000,-0.006250};
+Point(2544) = {0.540404,0.825000,-0.006738};
+Point(2545) = {0.537599,0.825000,-0.007211};
+Point(2546) = {0.534778,0.825000,-0.007661};
+Point(2547) = {0.531942,0.825000,-0.008084};
+Point(2548) = {0.529095,0.825000,-0.008478};
+Point(2549) = {0.526240,0.825000,-0.008841};
+Point(2550) = {0.523379,0.825000,-0.009172};
+Point(2551) = {0.520515,0.825000,-0.009469};
+Point(2552) = {0.517651,0.825000,-0.009732};
+Point(2553) = {0.514790,0.825000,-0.009960};
+Point(2554) = {0.511934,0.825000,-0.010156};
+Point(2555) = {0.509087,0.825000,-0.010319};
+Point(2556) = {0.506252,0.825000,-0.010451};
+Point(2557) = {0.503430,0.825000,-0.010552};
+Point(2558) = {0.500626,0.825000,-0.010626};
+Point(2559) = {0.497840,0.825000,-0.010673};
+Point(2560) = {0.495078,0.825000,-0.010697};
+Point(2561) = {0.492340,0.825000,-0.010701};
+Point(2562) = {0.489630,0.825000,-0.010686};
+Point(2563) = {0.486951,0.825000,-0.010652};
+Point(2564) = {0.484305,0.825000,-0.010600};
+Point(2565) = {0.481694,0.825000,-0.010530};
+Point(2566) = {0.479122,0.825000,-0.010441};
+Point(2567) = {0.476591,0.825000,-0.010335};
+Point(2568) = {0.474103,0.825000,-0.010211};
+Point(2569) = {0.471661,0.825000,-0.010071};
+Point(2570) = {0.469267,0.825000,-0.009916};
+Point(2571) = {0.466923,0.825000,-0.009746};
+Point(2572) = {0.464633,0.825000,-0.009563};
+Point(2573) = {0.462398,0.825000,-0.009368};
+Point(2574) = {0.460220,0.825000,-0.009161};
+Point(2575) = {0.458101,0.825000,-0.008946};
+Point(2576) = {0.456044,0.825000,-0.008723};
+Point(2577) = {0.454051,0.825000,-0.008494};
+Point(2578) = {0.452123,0.825000,-0.008259};
+Point(2579) = {0.450263,0.825000,-0.008021};
+Point(2580) = {0.448472,0.825000,-0.007779};
+Point(2581) = {0.446753,0.825000,-0.007536};
+Point(2582) = {0.445106,0.825000,-0.007291};
+Point(2583) = {0.443533,0.825000,-0.007046};
+Point(2584) = {0.442037,0.825000,-0.006799};
+Point(2585) = {0.440618,0.825000,-0.006551};
+Point(2586) = {0.439277,0.825000,-0.006301};
+Point(2587) = {0.438017,0.825000,-0.006048};
+Point(2588) = {0.436838,0.825000,-0.005790};
+Point(2589) = {0.435742,0.825000,-0.005524};
+Point(2590) = {0.434730,0.825000,-0.005245};
+Point(2591) = {0.433802,0.825000,-0.004952};
+Point(2592) = {0.432960,0.825000,-0.004641};
+Point(2593) = {0.432204,0.825000,-0.004312};
+Point(2594) = {0.431536,0.825000,-0.003967};
+Point(2595) = {0.430955,0.825000,-0.003605};
+Point(2596) = {0.430462,0.825000,-0.003228};
+Point(2597) = {0.430059,0.825000,-0.002836};
+Point(2598) = {0.429745,0.825000,-0.002426};
+Point(2599) = {0.429520,0.825000,-0.001997};
+Point(2600) = {0.429385,0.825000,-0.001544};
+Point(2601) = {0.429340,0.825000,-0.001067};
+Point(2602) = {0.429385,0.825000,-0.000563};
+Point(2603) = {0.429520,0.825000,-0.000040};
+Point(2604) = {0.429745,0.825000,0.000495};
+Point(2605) = {0.430059,0.825000,0.001032};
+Point(2606) = {0.430462,0.825000,0.001563};
+Point(2607) = {0.430955,0.825000,0.002082};
+Point(2608) = {0.431536,0.825000,0.002585};
+Point(2609) = {0.432204,0.825000,0.003072};
+Point(2610) = {0.432960,0.825000,0.003543};
+Point(2611) = {0.433802,0.825000,0.003995};
+Point(2612) = {0.434730,0.825000,0.004429};
+Point(2613) = {0.435742,0.825000,0.004844};
+Point(2614) = {0.436838,0.825000,0.005241};
+Point(2615) = {0.438017,0.825000,0.005621};
+Point(2616) = {0.439277,0.825000,0.005982};
+Point(2617) = {0.440618,0.825000,0.006326};
+Point(2618) = {0.442037,0.825000,0.006655};
+Point(2619) = {0.443533,0.825000,0.006970};
+Point(2620) = {0.445106,0.825000,0.007273};
+Point(2621) = {0.446753,0.825000,0.007566};
+Point(2622) = {0.448472,0.825000,0.007850};
+Point(2623) = {0.450263,0.825000,0.008125};
+Point(2624) = {0.452123,0.825000,0.008391};
+Point(2625) = {0.454051,0.825000,0.008645};
+Point(2626) = {0.456044,0.825000,0.008889};
+Point(2627) = {0.458101,0.825000,0.009121};
+Point(2628) = {0.460220,0.825000,0.009341};
+Point(2629) = {0.462398,0.825000,0.009550};
+Point(2630) = {0.464633,0.825000,0.009747};
+Point(2631) = {0.466923,0.825000,0.009932};
+Point(2632) = {0.469267,0.825000,0.010106};
+Point(2633) = {0.471661,0.825000,0.010268};
+Point(2634) = {0.474103,0.825000,0.010419};
+Point(2635) = {0.476591,0.825000,0.010558};
+Point(2636) = {0.479122,0.825000,0.010685};
+Point(2637) = {0.481694,0.825000,0.010801};
+Point(2638) = {0.484305,0.825000,0.010905};
+Point(2639) = {0.486951,0.825000,0.010998};
+Point(2640) = {0.489630,0.825000,0.011078};
+Point(2641) = {0.492340,0.825000,0.011147};
+Point(2642) = {0.495078,0.825000,0.011204};
+Point(2643) = {0.497840,0.825000,0.011251};
+Point(2644) = {0.500626,0.825000,0.011286};
+Point(2645) = {0.503430,0.825000,0.011311};
+Point(2646) = {0.506252,0.825000,0.011326};
+Point(2647) = {0.509087,0.825000,0.011331};
+Point(2648) = {0.511934,0.825000,0.011325};
+Point(2649) = {0.514790,0.825000,0.011308};
+Point(2650) = {0.517651,0.825000,0.011277};
+Point(2651) = {0.520515,0.825000,0.011232};
+Point(2652) = {0.523379,0.825000,0.011173};
+Point(2653) = {0.526240,0.825000,0.011099};
+Point(2654) = {0.529095,0.825000,0.011011};
+Point(2655) = {0.531942,0.825000,0.010911};
+Point(2656) = {0.534778,0.825000,0.010800};
+Point(2657) = {0.537599,0.825000,0.010679};
+Point(2658) = {0.540404,0.825000,0.010547};
+Point(2659) = {0.543189,0.825000,0.010405};
+Point(2660) = {0.545952,0.825000,0.010251};
+Point(2661) = {0.548689,0.825000,0.010086};
+Point(2662) = {0.551399,0.825000,0.009909};
+Point(2663) = {0.554078,0.825000,0.009720};
+Point(2664) = {0.556725,0.825000,0.009517};
+Point(2665) = {0.559335,0.825000,0.009299};
+Point(2666) = {0.561907,0.825000,0.009065};
+Point(2667) = {0.564439,0.825000,0.008814};
+Point(2668) = {0.566926,0.825000,0.008548};
+Point(2669) = {0.569369,0.825000,0.008267};
+Point(2670) = {0.571763,0.825000,0.007973};
+Point(2671) = {0.574106,0.825000,0.007667};
+Point(2672) = {0.576396,0.825000,0.007350};
+Point(2673) = {0.578632,0.825000,0.007026};
+Point(2674) = {0.580810,0.825000,0.006698};
+Point(2675) = {0.582928,0.825000,0.006369};
+Point(2676) = {0.584985,0.825000,0.006043};
+Point(2677) = {0.586978,0.825000,0.005722};
+Point(2678) = {0.588906,0.825000,0.005407};
+Point(2679) = {0.590766,0.825000,0.005100};
+Point(2680) = {0.592557,0.825000,0.004801};
+Point(2681) = {0.594277,0.825000,0.004510};
+Point(2682) = {0.595924,0.825000,0.004228};
+Point(2683) = {0.597496,0.825000,0.003955};
+Point(2684) = {0.598993,0.825000,0.003694};
+Point(2685) = {0.600412,0.825000,0.003444};
+Point(2686) = {0.601752,0.825000,0.003207};
+Point(2687) = {0.603012,0.825000,0.002983};
+Point(2688) = {0.604191,0.825000,0.002774};
+Point(2689) = {0.605287,0.825000,0.002578};
+Point(2690) = {0.606300,0.825000,0.002398};
+Point(2691) = {0.607227,0.825000,0.002233};
+Point(2692) = {0.608069,0.825000,0.002083};
+Point(2693) = {0.608825,0.825000,0.001915};
+Point(2694) = {0.609493,0.825000,0.001763};
+Point(2695) = {0.610074,0.825000,0.001628};
+Point(2696) = {0.610567,0.825000,0.001511};
+Point(2697) = {0.610971,0.825000,0.001414};
+Point(2698) = {0.611286,0.825000,0.001338};
+Point(2699) = {0.611511,0.825000,0.001284};
+Point(2700) = {0.611646,0.825000,0.001254};
+// -- Airfoil 6
+Point(3001) = {0.649731,0.950000,0.002532};
+Point(3002) = {0.649694,0.950000,0.002533};
+Point(3003) = {0.649579,0.950000,0.002537};
+Point(3004) = {0.649388,0.950000,0.002544};
+Point(3005) = {0.649119,0.950000,0.002553};
+Point(3006) = {0.648775,0.950000,0.002563};
+Point(3007) = {0.648355,0.950000,0.002574};
+Point(3008) = {0.647860,0.950000,0.002585};
+Point(3009) = {0.647291,0.950000,0.002594};
+Point(3010) = {0.646647,0.950000,0.002601};
+Point(3011) = {0.645930,0.950000,0.002649};
+Point(3012) = {0.645139,0.950000,0.002694};
+Point(3013) = {0.644277,0.950000,0.002734};
+Point(3014) = {0.643343,0.950000,0.002768};
+Point(3015) = {0.642339,0.950000,0.002795};
+Point(3016) = {0.641266,0.950000,0.002813};
+Point(3017) = {0.640124,0.950000,0.002821};
+Point(3018) = {0.638915,0.950000,0.002815};
+Point(3019) = {0.637640,0.950000,0.002795};
+Point(3020) = {0.636300,0.950000,0.002756};
+Point(3021) = {0.634897,0.950000,0.002698};
+Point(3022) = {0.633432,0.950000,0.002618};
+Point(3023) = {0.631907,0.950000,0.002515};
+Point(3024) = {0.630322,0.950000,0.002390};
+Point(3025) = {0.628680,0.950000,0.002242};
+Point(3026) = {0.626982,0.950000,0.002071};
+Point(3027) = {0.625230,0.950000,0.001877};
+Point(3028) = {0.623425,0.950000,0.001657};
+Point(3029) = {0.621570,0.950000,0.001409};
+Point(3030) = {0.619666,0.950000,0.001128};
+Point(3031) = {0.617715,0.950000,0.000812};
+Point(3032) = {0.615718,0.950000,0.000460};
+Point(3033) = {0.613679,0.950000,0.000078};
+Point(3034) = {0.611598,0.950000,-0.000328};
+Point(3035) = {0.609479,0.950000,-0.000751};
+Point(3036) = {0.607323,0.950000,-0.001183};
+Point(3037) = {0.605132,0.950000,-0.001618};
+Point(3038) = {0.602908,0.950000,-0.002057};
+Point(3039) = {0.600653,0.950000,-0.002499};
+Point(3040) = {0.598371,0.950000,-0.002946};
+Point(3041) = {0.596062,0.950000,-0.003399};
+Point(3042) = {0.593730,0.950000,-0.003857};
+Point(3043) = {0.591377,0.950000,-0.004316};
+Point(3044) = {0.589004,0.950000,-0.004772};
+Point(3045) = {0.586615,0.950000,-0.005218};
+Point(3046) = {0.584211,0.950000,-0.005650};
+Point(3047) = {0.581796,0.950000,-0.006062};
+Point(3048) = {0.579371,0.950000,-0.006450};
+Point(3049) = {0.576938,0.950000,-0.006813};
+Point(3050) = {0.574501,0.950000,-0.007146};
+Point(3051) = {0.572061,0.950000,-0.007447};
+Point(3052) = {0.569622,0.950000,-0.007714};
+Point(3053) = {0.567184,0.950000,-0.007947};
+Point(3054) = {0.564752,0.950000,-0.008151};
+Point(3055) = {0.562326,0.950000,-0.008327};
+Point(3056) = {0.559911,0.950000,-0.008478};
+Point(3057) = {0.557507,0.950000,-0.008607};
+Point(3058) = {0.555118,0.950000,-0.008715};
+Point(3059) = {0.552745,0.950000,-0.008802};
+Point(3060) = {0.550392,0.950000,-0.008870};
+Point(3061) = {0.548060,0.950000,-0.008919};
+Point(3062) = {0.545751,0.950000,-0.008950};
+Point(3063) = {0.543469,0.950000,-0.008965};
+Point(3064) = {0.541215,0.950000,-0.008964};
+Point(3065) = {0.538991,0.950000,-0.008950};
+Point(3066) = {0.536800,0.950000,-0.008924};
+Point(3067) = {0.534643,0.950000,-0.008887};
+Point(3068) = {0.532524,0.950000,-0.008840};
+Point(3069) = {0.530444,0.950000,-0.008782};
+Point(3070) = {0.528404,0.950000,-0.008712};
+Point(3071) = {0.526408,0.950000,-0.008631};
+Point(3072) = {0.524457,0.950000,-0.008539};
+Point(3073) = {0.522552,0.950000,-0.008435};
+Point(3074) = {0.520697,0.950000,-0.008323};
+Point(3075) = {0.518893,0.950000,-0.008203};
+Point(3076) = {0.517140,0.950000,-0.008077};
+Point(3077) = {0.515442,0.950000,-0.007945};
+Point(3078) = {0.513800,0.950000,-0.007810};
+Point(3079) = {0.512216,0.950000,-0.007673};
+Point(3080) = {0.510690,0.950000,-0.007536};
+Point(3081) = {0.509225,0.950000,-0.007400};
+Point(3082) = {0.507822,0.950000,-0.007265};
+Point(3083) = {0.506482,0.950000,-0.007128};
+Point(3084) = {0.505208,0.950000,-0.006984};
+Point(3085) = {0.503999,0.950000,-0.006826};
+Point(3086) = {0.502857,0.950000,-0.006651};
+Point(3087) = {0.501783,0.950000,-0.006455};
+Point(3088) = {0.500779,0.950000,-0.006238};
+Point(3089) = {0.499845,0.950000,-0.006006};
+Point(3090) = {0.498983,0.950000,-0.005761};
+Point(3091) = {0.498193,0.950000,-0.005506};
+Point(3092) = {0.497475,0.950000,-0.005246};
+Point(3093) = {0.496831,0.950000,-0.004976};
+Point(3094) = {0.496262,0.950000,-0.004697};
+Point(3095) = {0.495767,0.950000,-0.004406};
+Point(3096) = {0.495347,0.950000,-0.004101};
+Point(3097) = {0.495004,0.950000,-0.003780};
+Point(3098) = {0.494736,0.950000,-0.003444};
+Point(3099) = {0.494545,0.950000,-0.003094};
+Point(3100) = {0.494430,0.950000,-0.002730};
+Point(3101) = {0.494391,0.950000,-0.002353};
+Point(3102) = {0.494430,0.950000,-0.001965};
+Point(3103) = {0.494545,0.950000,-0.001565};
+Point(3104) = {0.494736,0.950000,-0.001153};
+Point(3105) = {0.495004,0.950000,-0.000729};
+Point(3106) = {0.495347,0.950000,-0.000293};
+Point(3107) = {0.495767,0.950000,0.000155};
+Point(3108) = {0.496262,0.950000,0.000609};
+Point(3109) = {0.496831,0.950000,0.001062};
+Point(3110) = {0.497475,0.950000,0.001509};
+Point(3111) = {0.498193,0.950000,0.001942};
+Point(3112) = {0.498983,0.950000,0.002357};
+Point(3113) = {0.499845,0.950000,0.002754};
+Point(3114) = {0.500779,0.950000,0.003135};
+Point(3115) = {0.501783,0.950000,0.003503};
+Point(3116) = {0.502857,0.950000,0.003859};
+Point(3117) = {0.503999,0.950000,0.004205};
+Point(3118) = {0.505208,0.950000,0.004542};
+Point(3119) = {0.506482,0.950000,0.004869};
+Point(3120) = {0.507822,0.950000,0.005187};
+Point(3121) = {0.509225,0.950000,0.005496};
+Point(3122) = {0.510690,0.950000,0.005796};
+Point(3123) = {0.512216,0.950000,0.006087};
+Point(3124) = {0.513800,0.950000,0.006369};
+Point(3125) = {0.515442,0.950000,0.006642};
+Point(3126) = {0.517140,0.950000,0.006905};
+Point(3127) = {0.518893,0.950000,0.007159};
+Point(3128) = {0.520697,0.950000,0.007404};
+Point(3129) = {0.522552,0.950000,0.007639};
+Point(3130) = {0.524457,0.950000,0.007865};
+Point(3131) = {0.526408,0.950000,0.008082};
+Point(3132) = {0.528404,0.950000,0.008289};
+Point(3133) = {0.530444,0.950000,0.008486};
+Point(3134) = {0.532524,0.950000,0.008673};
+Point(3135) = {0.534643,0.950000,0.008850};
+Point(3136) = {0.536800,0.950000,0.009017};
+Point(3137) = {0.538991,0.950000,0.009173};
+Point(3138) = {0.541215,0.950000,0.009318};
+Point(3139) = {0.543469,0.950000,0.009453};
+Point(3140) = {0.545751,0.950000,0.009579};
+Point(3141) = {0.548060,0.950000,0.009695};
+Point(3142) = {0.550392,0.950000,0.009803};
+Point(3143) = {0.552745,0.950000,0.009901};
+Point(3144) = {0.555118,0.950000,0.009990};
+Point(3145) = {0.557507,0.950000,0.010069};
+Point(3146) = {0.559911,0.950000,0.010137};
+Point(3147) = {0.562326,0.950000,0.010195};
+Point(3148) = {0.564752,0.950000,0.010242};
+Point(3149) = {0.567184,0.950000,0.010280};
+Point(3150) = {0.569622,0.950000,0.010307};
+Point(3151) = {0.572061,0.950000,0.010325};
+Point(3152) = {0.574501,0.950000,0.010335};
+Point(3153) = {0.576938,0.950000,0.010334};
+Point(3154) = {0.579371,0.950000,0.010324};
+Point(3155) = {0.581796,0.950000,0.010303};
+Point(3156) = {0.584211,0.950000,0.010271};
+Point(3157) = {0.586615,0.950000,0.010227};
+Point(3158) = {0.589004,0.950000,0.010171};
+Point(3159) = {0.591377,0.950000,0.010103};
+Point(3160) = {0.593730,0.950000,0.010023};
+Point(3161) = {0.596062,0.950000,0.009930};
+Point(3162) = {0.598371,0.950000,0.009825};
+Point(3163) = {0.600653,0.950000,0.009706};
+Point(3164) = {0.602908,0.950000,0.009573};
+Point(3165) = {0.605132,0.950000,0.009425};
+Point(3166) = {0.607323,0.950000,0.009262};
+Point(3167) = {0.609479,0.950000,0.009083};
+Point(3168) = {0.611598,0.950000,0.008887};
+Point(3169) = {0.613679,0.950000,0.008674};
+Point(3170) = {0.615718,0.950000,0.008445};
+Point(3171) = {0.617715,0.950000,0.008199};
+Point(3172) = {0.619666,0.950000,0.007937};
+Point(3173) = {0.621570,0.950000,0.007662};
+Point(3174) = {0.623425,0.950000,0.007380};
+Point(3175) = {0.625230,0.950000,0.007094};
+Point(3176) = {0.626982,0.950000,0.006809};
+Point(3177) = {0.628680,0.950000,0.006529};
+Point(3178) = {0.630322,0.950000,0.006254};
+Point(3179) = {0.631907,0.950000,0.005985};
+Point(3180) = {0.633432,0.950000,0.005723};
+Point(3181) = {0.634897,0.950000,0.005468};
+Point(3182) = {0.636300,0.950000,0.005221};
+Point(3183) = {0.637640,0.950000,0.004981};
+Point(3184) = {0.638915,0.950000,0.004751};
+Point(3185) = {0.640124,0.950000,0.004531};
+Point(3186) = {0.641266,0.950000,0.004322};
+Point(3187) = {0.642339,0.950000,0.004126};
+Point(3188) = {0.643343,0.950000,0.003941};
+Point(3189) = {0.644277,0.950000,0.003769};
+Point(3190) = {0.645139,0.950000,0.003610};
+Point(3191) = {0.645930,0.950000,0.003464};
+Point(3192) = {0.646647,0.950000,0.003332};
+Point(3193) = {0.647291,0.950000,0.003172};
+Point(3194) = {0.647860,0.950000,0.003028};
+Point(3195) = {0.648355,0.950000,0.002900};
+Point(3196) = {0.648775,0.950000,0.002790};
+Point(3197) = {0.649119,0.950000,0.002698};
+Point(3198) = {0.649388,0.950000,0.002625};
+Point(3199) = {0.649579,0.950000,0.002573};
+Point(3200) = {0.649694,0.950000,0.002542};
+// -- Airfoil 7
+Point(3501) = {0.664862,1.000000,0.002901};
+Point(3502) = {0.664827,1.000000,0.002902};
+Point(3503) = {0.664721,1.000000,0.002906};
+Point(3504) = {0.664542,1.000000,0.002910};
+Point(3505) = {0.664293,1.000000,0.002915};
+Point(3506) = {0.663973,1.000000,0.002920};
+Point(3507) = {0.663582,1.000000,0.002925};
+Point(3508) = {0.663122,1.000000,0.002930};
+Point(3509) = {0.662592,1.000000,0.002933};
+Point(3510) = {0.661994,1.000000,0.002936};
+Point(3511) = {0.661327,1.000000,0.002986};
+Point(3512) = {0.660592,1.000000,0.003039};
+Point(3513) = {0.659790,1.000000,0.003091};
+Point(3514) = {0.658922,1.000000,0.003138};
+Point(3515) = {0.657988,1.000000,0.003179};
+Point(3516) = {0.656990,1.000000,0.003209};
+Point(3517) = {0.655928,1.000000,0.003226};
+Point(3518) = {0.654804,1.000000,0.003227};
+Point(3519) = {0.653618,1.000000,0.003212};
+Point(3520) = {0.652373,1.000000,0.003178};
+Point(3521) = {0.651068,1.000000,0.003125};
+Point(3522) = {0.649706,1.000000,0.003050};
+Point(3523) = {0.648287,1.000000,0.002954};
+Point(3524) = {0.646813,1.000000,0.002833};
+Point(3525) = {0.645286,1.000000,0.002689};
+Point(3526) = {0.643708,1.000000,0.002518};
+Point(3527) = {0.642078,1.000000,0.002321};
+Point(3528) = {0.640400,1.000000,0.002100};
+Point(3529) = {0.638675,1.000000,0.001854};
+Point(3530) = {0.636904,1.000000,0.001586};
+Point(3531) = {0.635090,1.000000,0.001298};
+Point(3532) = {0.633233,1.000000,0.000990};
+Point(3533) = {0.631337,1.000000,0.000662};
+Point(3534) = {0.629402,1.000000,0.000317};
+Point(3535) = {0.627431,1.000000,-0.000047};
+Point(3536) = {0.625426,1.000000,-0.000429};
+Point(3537) = {0.623389,1.000000,-0.000828};
+Point(3538) = {0.621321,1.000000,-0.001241};
+Point(3539) = {0.619225,1.000000,-0.001667};
+Point(3540) = {0.617102,1.000000,-0.002103};
+Point(3541) = {0.614955,1.000000,-0.002546};
+Point(3542) = {0.612787,1.000000,-0.002994};
+Point(3543) = {0.610598,1.000000,-0.003442};
+Point(3544) = {0.608392,1.000000,-0.003885};
+Point(3545) = {0.606170,1.000000,-0.004320};
+Point(3546) = {0.603935,1.000000,-0.004740};
+Point(3547) = {0.601689,1.000000,-0.005142};
+Point(3548) = {0.599434,1.000000,-0.005523};
+Point(3549) = {0.597172,1.000000,-0.005881};
+Point(3550) = {0.594905,1.000000,-0.006213};
+Point(3551) = {0.592637,1.000000,-0.006517};
+Point(3552) = {0.590368,1.000000,-0.006792};
+Point(3553) = {0.588102,1.000000,-0.007038};
+Point(3554) = {0.585840,1.000000,-0.007257};
+Point(3555) = {0.583585,1.000000,-0.007451};
+Point(3556) = {0.581338,1.000000,-0.007622};
+Point(3557) = {0.579103,1.000000,-0.007771};
+Point(3558) = {0.576881,1.000000,-0.007899};
+Point(3559) = {0.574675,1.000000,-0.008007};
+Point(3560) = {0.572487,1.000000,-0.008096};
+Point(3561) = {0.570318,1.000000,-0.008167};
+Point(3562) = {0.568171,1.000000,-0.008221};
+Point(3563) = {0.566049,1.000000,-0.008259};
+Point(3564) = {0.563953,1.000000,-0.008283};
+Point(3565) = {0.561885,1.000000,-0.008296};
+Point(3566) = {0.559847,1.000000,-0.008298};
+Point(3567) = {0.557842,1.000000,-0.008292};
+Point(3568) = {0.555871,1.000000,-0.008276};
+Point(3569) = {0.553937,1.000000,-0.008251};
+Point(3570) = {0.552040,1.000000,-0.008217};
+Point(3571) = {0.550184,1.000000,-0.008172};
+Point(3572) = {0.548370,1.000000,-0.008116};
+Point(3573) = {0.546599,1.000000,-0.008051};
+Point(3574) = {0.544874,1.000000,-0.007977};
+Point(3575) = {0.543195,1.000000,-0.007895};
+Point(3576) = {0.541566,1.000000,-0.007807};
+Point(3577) = {0.539987,1.000000,-0.007714};
+Point(3578) = {0.538460,1.000000,-0.007616};
+Point(3579) = {0.536987,1.000000,-0.007514};
+Point(3580) = {0.535568,1.000000,-0.007408};
+Point(3581) = {0.534206,1.000000,-0.007298};
+Point(3582) = {0.532901,1.000000,-0.007185};
+Point(3583) = {0.531655,1.000000,-0.007068};
+Point(3584) = {0.530470,1.000000,-0.006947};
+Point(3585) = {0.529346,1.000000,-0.006819};
+Point(3586) = {0.528284,1.000000,-0.006686};
+Point(3587) = {0.527286,1.000000,-0.006545};
+Point(3588) = {0.526352,1.000000,-0.006396};
+Point(3589) = {0.525484,1.000000,-0.006235};
+Point(3590) = {0.524682,1.000000,-0.006061};
+Point(3591) = {0.523947,1.000000,-0.005872};
+Point(3592) = {0.523280,1.000000,-0.005667};
+Point(3593) = {0.522681,1.000000,-0.005445};
+Point(3594) = {0.522151,1.000000,-0.005209};
+Point(3595) = {0.521691,1.000000,-0.004958};
+Point(3596) = {0.521301,1.000000,-0.004693};
+Point(3597) = {0.520981,1.000000,-0.004416};
+Point(3598) = {0.520732,1.000000,-0.004124};
+Point(3599) = {0.520554,1.000000,-0.003813};
+Point(3600) = {0.520447,1.000000,-0.003481};
+Point(3601) = {0.520412,1.000000,-0.003124};
+Point(3602) = {0.520447,1.000000,-0.002742};
+Point(3603) = {0.520554,1.000000,-0.002338};
+Point(3604) = {0.520732,1.000000,-0.001916};
+Point(3605) = {0.520981,1.000000,-0.001482};
+Point(3606) = {0.521301,1.000000,-0.001041};
+Point(3607) = {0.521691,1.000000,-0.000599};
+Point(3608) = {0.522151,1.000000,-0.000159};
+Point(3609) = {0.522681,1.000000,0.000273};
+Point(3610) = {0.523280,1.000000,0.000691};
+Point(3611) = {0.523947,1.000000,0.001092};
+Point(3612) = {0.524682,1.000000,0.001471};
+Point(3613) = {0.525484,1.000000,0.001833};
+Point(3614) = {0.526352,1.000000,0.002185};
+Point(3615) = {0.527286,1.000000,0.002533};
+Point(3616) = {0.528284,1.000000,0.002884};
+Point(3617) = {0.529346,1.000000,0.003241};
+Point(3618) = {0.530470,1.000000,0.003601};
+Point(3619) = {0.531655,1.000000,0.003956};
+Point(3620) = {0.532901,1.000000,0.004301};
+Point(3621) = {0.534206,1.000000,0.004628};
+Point(3622) = {0.535568,1.000000,0.004934};
+Point(3623) = {0.536987,1.000000,0.005220};
+Point(3624) = {0.538460,1.000000,0.005493};
+Point(3625) = {0.539987,1.000000,0.005756};
+Point(3626) = {0.541566,1.000000,0.006015};
+Point(3627) = {0.543195,1.000000,0.006273};
+Point(3628) = {0.544874,1.000000,0.006530};
+Point(3629) = {0.546599,1.000000,0.006782};
+Point(3630) = {0.548370,1.000000,0.007027};
+Point(3631) = {0.550184,1.000000,0.007263};
+Point(3632) = {0.552040,1.000000,0.007487};
+Point(3633) = {0.553937,1.000000,0.007700};
+Point(3634) = {0.555871,1.000000,0.007903};
+Point(3635) = {0.557842,1.000000,0.008094};
+Point(3636) = {0.559847,1.000000,0.008275};
+Point(3637) = {0.561885,1.000000,0.008446};
+Point(3638) = {0.563953,1.000000,0.008606};
+Point(3639) = {0.566049,1.000000,0.008758};
+Point(3640) = {0.568171,1.000000,0.008899};
+Point(3641) = {0.570318,1.000000,0.009032};
+Point(3642) = {0.572487,1.000000,0.009156};
+Point(3643) = {0.574675,1.000000,0.009271};
+Point(3644) = {0.576881,1.000000,0.009376};
+Point(3645) = {0.579103,1.000000,0.009472};
+Point(3646) = {0.581338,1.000000,0.009559};
+Point(3647) = {0.583585,1.000000,0.009636};
+Point(3648) = {0.585840,1.000000,0.009704};
+Point(3649) = {0.588102,1.000000,0.009763};
+Point(3650) = {0.590368,1.000000,0.009812};
+Point(3651) = {0.592637,1.000000,0.009851};
+Point(3652) = {0.594905,1.000000,0.009881};
+Point(3653) = {0.597172,1.000000,0.009902};
+Point(3654) = {0.599434,1.000000,0.009914};
+Point(3655) = {0.601689,1.000000,0.009916};
+Point(3656) = {0.603935,1.000000,0.009910};
+Point(3657) = {0.606170,1.000000,0.009894};
+Point(3658) = {0.608392,1.000000,0.009868};
+Point(3659) = {0.610598,1.000000,0.009831};
+Point(3660) = {0.612787,1.000000,0.009783};
+Point(3661) = {0.614955,1.000000,0.009721};
+Point(3662) = {0.617102,1.000000,0.009646};
+Point(3663) = {0.619225,1.000000,0.009556};
+Point(3664) = {0.621321,1.000000,0.009451};
+Point(3665) = {0.623389,1.000000,0.009331};
+Point(3666) = {0.625426,1.000000,0.009194};
+Point(3667) = {0.627431,1.000000,0.009041};
+Point(3668) = {0.629402,1.000000,0.008871};
+Point(3669) = {0.631337,1.000000,0.008684};
+Point(3670) = {0.633233,1.000000,0.008481};
+Point(3671) = {0.635090,1.000000,0.008262};
+Point(3672) = {0.636904,1.000000,0.008027};
+Point(3673) = {0.638675,1.000000,0.007780};
+Point(3674) = {0.640400,1.000000,0.007523};
+Point(3675) = {0.642078,1.000000,0.007259};
+Point(3676) = {0.643708,1.000000,0.006993};
+Point(3677) = {0.645286,1.000000,0.006727};
+Point(3678) = {0.646813,1.000000,0.006463};
+Point(3679) = {0.648287,1.000000,0.006203};
+Point(3680) = {0.649706,1.000000,0.005948};
+Point(3681) = {0.651068,1.000000,0.005700};
+Point(3682) = {0.652373,1.000000,0.005461};
+Point(3683) = {0.653618,1.000000,0.005231};
+Point(3684) = {0.654804,1.000000,0.005011};
+Point(3685) = {0.655928,1.000000,0.004803};
+Point(3686) = {0.656990,1.000000,0.004606};
+Point(3687) = {0.657988,1.000000,0.004422};
+Point(3688) = {0.658922,1.000000,0.004250};
+Point(3689) = {0.659790,1.000000,0.004090};
+Point(3690) = {0.660592,1.000000,0.003942};
+Point(3691) = {0.661327,1.000000,0.003805};
+Point(3692) = {0.661994,1.000000,0.003680};
+Point(3693) = {0.662592,1.000000,0.003522};
+Point(3694) = {0.663122,1.000000,0.003381};
+Point(3695) = {0.663582,1.000000,0.003257};
+Point(3696) = {0.663973,1.000000,0.003151};
+Point(3697) = {0.664293,1.000000,0.003062};
+Point(3698) = {0.664542,1.000000,0.002993};
+Point(3699) = {0.664721,1.000000,0.002942};
+Point(3700) = {0.664827,1.000000,0.002912};
+
+// --- Wing lines ---
+// -- Airfoil 0
+Spline(1) = {1:101};
+Spline(2) = {101:200, 1};
+// -- Airfoil 1
+Spline(11) = {501:601};
+Spline(12) = {601:700, 501};
+// -- Airfoil 2
+Spline(21) = {1001:1101};
+Spline(22) = {1101:1200, 1001};
+// -- Airfoil 3
+Spline(31) = {1501:1601};
+Spline(32) = {1601:1700, 1501};
+// -- Airfoil 4
+Spline(41) = {2001:2101};
+Spline(42) = {2101:2200, 2001};
+// -- Airfoil 5
+Spline(51) = {2501:2601};
+Spline(52) = {2601:2700, 2501};
+// -- Airfoil 6
+Spline(61) = {3001:3101};
+Spline(62) = {3101:3200, 3001};
+// -- Airfoil 7
+Spline(71) = {3501:3601};
+Spline(72) = {3601:3700, 3501};
+// -- Planform 0
+Line(101) = {1,501};
+Line(102) = {101,601};
+// -- Planform 1
+Line(111) = {501,1001};
+Line(112) = {601,1101};
+// -- Planform 2
+Line(121) = {1001,1501};
+Line(122) = {1101,1601};
+// -- Planform 3
+Line(131) = {1501,2001};
+Line(132) = {1601,2101};
+// -- Planform 4
+Line(141) = {2001,2501};
+Line(142) = {2101,2601};
+// -- Planform 5
+Line(151) = {2501,3001};
+Line(152) = {2601,3101};
+// -- Planform 6
+Line(161) = {3001,3501};
+Line(162) = {3101,3601};
+
+// --- Wing line loops and surfaces ---
+// -- Planform 0
+Line Loop(1) = {1,102,-11,-101};
+Line Loop(2) = {2,101,-12,-102};
+Surface(1) = {1};
+Surface(2) = {2};
+// -- Planform 1
+Line Loop(11) = {11,112,-21,-111};
+Line Loop(12) = {12,111,-22,-112};
+Surface(11) = {11};
+Surface(12) = {12};
+// -- Planform 2
+Line Loop(21) = {21,122,-31,-121};
+Line Loop(22) = {22,121,-32,-122};
+Surface(21) = {21};
+Surface(22) = {22};
+// -- Planform 3
+Line Loop(31) = {31,132,-41,-131};
+Line Loop(32) = {32,131,-42,-132};
+Surface(31) = {31};
+Surface(32) = {32};
+// -- Planform 4
+Line Loop(41) = {41,142,-51,-141};
+Line Loop(42) = {42,141,-52,-142};
+Surface(41) = {41};
+Surface(42) = {42};
+// -- Planform 5
+Line Loop(51) = {51,152,-61,-151};
+Line Loop(52) = {52,151,-62,-152};
+Surface(51) = {51};
+Surface(52) = {52};
+// -- Planform 6
+Line Loop(61) = {61,162,-71,-161};
+Line Loop(62) = {62,161,-72,-162};
+Surface(61) = {61};
+Surface(62) = {62};
+
+// --- Mesh ---
+Transfinite Line{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62, 71, 72} = nC + 1 Using Bump bC;
+Transfinite Line{101, 102} = 5 * rS + 1;
+Transfinite Line{111, 112} = 3 * rS + 1;
+Transfinite Line{121, 122} = 4 * rS + 1;
+Transfinite Line{131, 132} = 4 * rS + 1;
+Transfinite Line{141, 142} = 4 * rS + 1;
+Transfinite Line{151, 152} = 3 * rS + 1;
+Transfinite Line{161, 162} = 1 * rS + 1;
+Transfinite Surface{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62};
+Recombine Surface{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62};
+
+// --- Wing physical groups ---
+Physical Line("te") = {101, 111, 121, 131, 141, 151, 161};
+Physical Surface("wing") = {2, 12, 22, 32, 42, 52, 62, 1, 11, 21, 31, 41, 51, 61};
diff --git a/sdpm/models/oneraM6.geo b/sdpm/models/oneraM6.geo
index 418d940091afa3989699463a178221d2a97ba642..3821d3951f25d60e3d10da68ef77daa78df76552 100644
--- a/sdpm/models/oneraM6.geo
+++ b/sdpm/models/oneraM6.geo
@@ -2,8 +2,8 @@
 
 // Mesh parameters
 // domain and mesh
-DefineConstant[ wNc = { 50, Name "Number of panel along chord" }
-                wNs = { 10, Name "Number of panel along span" }
+DefineConstant[ nC = { 50, Name "Number of panel along chord" }
+                nS = { 10, Name "Number of panel along span" }
                 bC = { 0.05, Name "Progression along chord" }
                 pS = { 1.0, Name "Progression along span" } ];
 
@@ -315,13 +315,13 @@ Surface(11) = {-11};
 Surface(12) = {-12};
 
 //// MESH
-Transfinite Line{1,2,3,4} = wNc + 1 Using Bump bC;
-Transfinite Line{61,62} = wNs + 1 Using Progression pS;
+Transfinite Line{1,2,3,4} = nC + 1 Using Bump bC;
+Transfinite Line{61,62} = nS + 1 Using Progression pS;
 Transfinite Surface{11,12};
 Recombine Surface{11,12};
 
 //// PHYSICAL GROUPS
 // Trailing edge
-Physical Line("wTe") = {61};
+Physical Line("te") = {61};
 // Wing:
 Physical Surface("wing") = {11,12};
diff --git a/sdpm/src/sdpm.h b/sdpm/src/sdpm.h
index 4fe47a5edf1206f069ebf5f34591b9c07a7435d7..c580d78cdd219228b30bf29c8a507342ead3eb87 100644
--- a/sdpm/src/sdpm.h
+++ b/sdpm/src/sdpm.h
@@ -32,12 +32,17 @@
 #endif
 
 #include <Eigen/Dense>
+#include <complex>
 
 using sdpmDouble = double;
+using sdpmComplex = std::complex<double>;
 using sdpmMatrixXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, Eigen::Dynamic>;
 using sdpmVectorXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, 1>;
 using sdpmMatrix3d = Eigen::Matrix<sdpmDouble, 3, 3>;
 using sdpmVector3d = Eigen::Matrix<sdpmDouble, 3, 1>;
+using sdpmMatrixXcd = Eigen::Matrix<sdpmComplex, Eigen::Dynamic, Eigen::Dynamic>;
+using sdpmVectorXcd = Eigen::Matrix<sdpmComplex, Eigen::Dynamic, 1>;
+using sdpmVector3cd = Eigen::Matrix<sdpmComplex, 3, 1>;
 
 namespace sdpm
 {
@@ -60,6 +65,7 @@ class GmshExport;
 class Problem;
 class Body;
 class Wake;
+class Motion;
 class Edge;
 class Builder;
 class Solver;
diff --git a/sdpm/src/sdpmBody.cpp b/sdpm/src/sdpmBody.cpp
index 9f9e496c59ed099efa90316eddf1c581e4178880..1cae6e3cbc799edfb6ca45bf7c8987abb7acbf85 100644
--- a/sdpm/src/sdpmBody.cpp
+++ b/sdpm/src/sdpmBody.cpp
@@ -21,10 +21,19 @@
 #include "sdpmQuad4.h"
 #include "sdpmTag.h"
 #include "sdpmGmshExport.h"
+#include <unordered_set>
 #include <iostream>
 using namespace sdpm;
 
-Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double xF) : Group(msh, name), cl(0), cd(0), cs(0), cm(0)
+Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
+           double shedLgth, size_t nShedDiv,
+           size_t nM, size_t nF,
+           double ux) : Group(msh, name),
+                        _cl(0), _cd(0), _cs(0), _cm(0),
+                        _cl1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cd1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cs1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cm1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.)))
 {
     // Create wake
     std::cout << "Creating wake... " << std::flush;
@@ -43,30 +52,34 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
     auto it = tags.find(wkName);
     if (it == tags.end())
     {
+        // sanity check
+        if (shedLgth <= 0.)
+        {
+            std::stringstream err;
+            err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
+            throw std::runtime_error(err.str());
+        }
         // create tag and add it to the mesh
         Tag *tagp = new Tag(tags.size() + 1, wkName, 2);
         _msh.addTag(wkName, tagp);
-        // translate TE nodes (along x-coordinate)
-        std::vector<Node *> wkNodes;
-        for (auto n : teNodes)
+        // create wake nodes and elements
+        std::vector<Node *> wkNodes = teNodes;
+        for (size_t i = 0; i < 10 * nShedDiv; ++i)
         {
-            sdpmVector3d pos = n->getCoords();
-            // sanity check
-            if (xF <= pos(0))
+            // translate TE nodes along x-coordinate
+            sdpmVector3d delta((i + 1) * shedLgth / nShedDiv, 0., 0.);
+            for (auto n : teNodes)
             {
-                std::stringstream err;
-                err << "sdpm::Body: zero or negative length wake requested for " << *_tag << "!\n";
-                throw std::runtime_error(err.str());
+                Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords() + delta);
+                wkNodes.push_back(nN);
+                _msh.addNode(nN);
+            }
+            // create elements
+            for (size_t j = 0; j < teNodes.size() - 1; ++j)
+            {
+                std::vector<Node *> qnodes = {wkNodes[i * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j + 1], wkNodes[i * teNodes.size() + j + 1]};
+                _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
             }
-            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, sdpmVector3d(xF, pos(1), pos(2)));
-            _msh.addNode(nN);
-            wkNodes.push_back(nN);
-        }
-        // create wake elements
-        for (size_t i = 0; i < wkNodes.size() - 1; ++i)
-        {
-            std::vector<Node *> qnodes = {teNodes[i], wkNodes[i], wkNodes[i + 1], teNodes[i + 1]};
-            _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
         }
         // duplicate TE nodes
         std::map<Node *, Node *> teMap;
@@ -77,17 +90,19 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
             teMap[n] = nN;
         }
         // create wake
-        _wake = new Wake(_msh, wkName, _tag);
+        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, ux);
         // swap nodes
-        for (auto p : _wake->getMap())
+        std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
+        for (auto p : _wake->getElMap())
+            lwEl.insert(p.second.second);
+        for (auto e : lwEl)
         {
-            Element *lwEl = p.second.second;
-            std::vector<Node *> lwNods = lwEl->getNodes();
+            std::vector<Node *> lwNods = e->getNodes();
             for (size_t i = 0; i < lwNods.size(); ++i)
             {
                 try
                 {
-                    lwEl->setNode(i, teMap.at(lwNods[i]));
+                    e->setNode(i, teMap.at(lwNods[i]));
                 }
                 catch (const std::out_of_range &)
                 {
@@ -111,10 +126,21 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
         for (auto ten : teNodes)
             for (auto wn : _nodes)
                 if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
+                {
                     iTeMap[wn] = ten;
+                    break;
+                }
         // create wake and print
-        _wake = new Wake(_msh, wkName, _tag, iTeMap);
+        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, ux, iTeMap);
         std::cout << *_wake << "already existing, nothing done." << std::endl;
+        // check if parameters match existing wake geometry
+        std::vector<Element *> wEls = _wake->getElements();
+        if (std::abs(wEls[0]->getNodes()[1]->getCoords()(0) - wEls[0]->getNodes()[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || wEls.size() != 10 * nShedDiv * (teNodes.size() - 1))
+        {
+            std::stringstream err;
+            err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
+            throw std::runtime_error(err.str());
+        }
     }
 
     // Associate each node to its adjacent elements
@@ -122,16 +148,44 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
         for (auto n : e->getNodes())
             _neMap[n].push_back(e);
 
-    // Size load vectors
-    nLoads.resize(_nodes.size());
+    // Size vectors
+    _nLoads.resize(_nodes.size());
+    _nLoads1.resize(nM, std::vector<std::vector<sdpmVector3cd>>(nF, std::vector<sdpmVector3cd>(_nodes.size(), sdpmVector3cd::Zero())));
 }
 
 Body::~Body()
 {
     delete _wake;
+    for (auto m : _motions)
+        delete m;
+    _motions.clear();
     std::cout << "~sdpm::Body()\n";
 }
 
+/**
+ * @brief Add motion
+ */
+void Body::addMotion(double ux, double amp, double xr)
+{
+    _motions.push_back(new Motion(ux, amp, xr));
+}
+
+/**
+ * @brief Compute the motion-induced velocity on the body surface
+ */
+void Body::computeVelocity(size_t m)
+{
+    std::vector<Element *> elems = getElements();
+    std::vector<sdpmVector3cd> ucS(elems.size()), ucH(elems.size());
+    for (size_t i = 0; i < elems.size(); ++i)
+    {
+        ucS[i] = _motions[m]->computeSteady(*elems[i]);
+        ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
+    }
+    _iVelocityS.push_back(ucS);
+    _iVelocityH.push_back(ucH);
+}
+
 /**
  * @brief Return a vector of unique nodes belonging to the tag
  */
diff --git a/sdpm/src/sdpmBody.h b/sdpm/src/sdpmBody.h
index c3fa2c8076c18e81d750fe2f7e2862ebffe736d2..ae7e9244b4483f90948822ecea76adeda2a21fb3 100644
--- a/sdpm/src/sdpmBody.h
+++ b/sdpm/src/sdpmBody.h
@@ -20,6 +20,7 @@
 #include "sdpm.h"
 #include "sdpmGroup.h"
 #include "sdpmWake.h"
+#include "sdpmMotion.h"
 #include <map>
 
 namespace sdpm
@@ -29,20 +30,31 @@ namespace sdpm
  * @brief Manage a lifting body immersed in the fluid
  * @authors Adrien Crovato
  */
-class SDPM_API Body : public sdpm::Group
+class SDPM_API Body : public Group
 {
+    // Geometry
     std::vector<Node *> _nodes;                      ///< nodes of the surface
     std::map<Node *, std::vector<Element *>> _neMap; ///< map between nodes and adjacent elements
     Wake *_wake;                                     ///< wake attached to the lifting body
+    // Motion
+    std::vector<Motion *> _motions;                      ///< body motions
+    std::vector<std::vector<sdpmVector3cd>> _iVelocityS; ///< steady motion induced velocity
+    std::vector<std::vector<sdpmVector3cd>> _iVelocityH; ///< harmonic motion induced velocity
+    // Nodal aerodynamic load (normalized by dynamic pressure)
+    std::vector<sdpmVector3d> _nLoads;                             ///< steady nodal aerodynamic load
+    std::vector<std::vector<std::vector<sdpmVector3cd>>> _nLoads1; ///< unsteady first harmonic nodal aerodynamic load
+    //  Coefficients
+    sdpmDouble _cl;                             ///< steady lift coefficient
+    sdpmDouble _cd;                             ///< steady drag coefficient
+    sdpmDouble _cs;                             ///< steady side force coefficient
+    sdpmDouble _cm;                             ///< steady pitch moment coefficient (positive nose-up)
+    std::vector<std::vector<sdpmComplex>> _cl1; ///< unsteady first harmonic lift coefficient
+    std::vector<std::vector<sdpmComplex>> _cd1; ///< unsteady first harmonic drag coefficient
+    std::vector<std::vector<sdpmComplex>> _cs1; ///< unsteady first harmonic side force coefficient
+    std::vector<std::vector<sdpmComplex>> _cm1; ///< unsteady first harmonic pitch moment coefficient (positive nose-up)
 
 public:
-    sdpmDouble cl;                    ///< lift coefficient
-    sdpmDouble cd;                    ///< drag coefficient
-    sdpmDouble cs;                    ///< sideforce coefficient
-    sdpmDouble cm;                    ///< pitch moment coefficient (positive nose-up)
-    std::vector<sdpmVector3d> nLoads; ///< nodal aerodynamic load (normalized by dynamic pressure)
-
-    Body(Mesh &msh, std::string const &name, std::string const &teName, double xF);
+    Body(Mesh &msh, std::string const &name, std::string const &teName, double shedLgth, size_t nShedDiv = 1, size_t nM = 0, size_t nF = 0, double ux = 1.);
     ~Body();
 
 private:
@@ -50,12 +62,38 @@ private:
 
 public:
     std::vector<Node *> const &getNodes() const { return _nodes; }
+    std::vector<sdpmVector3d> const &getLoads() const { return _nLoads; }
+    std::vector<sdpmVector3cd> const &getUnsteadyLoads(size_t imd, size_t ifq) const { return _nLoads1[imd][ifq]; }
+    sdpmDouble getLiftCoeff() const { return _cl; }
+    sdpmDouble getDragCoeff() const { return _cd; }
+    sdpmDouble getSideforceCoeff() const { return _cs; }
+    sdpmDouble getMomentCoeff() const { return _cm; }
+    sdpmComplex getUnsteadyLiftCoeff(size_t imd, size_t ifq) const { return _cl1[imd][ifq]; }
+    sdpmComplex getUnsteadyDragCoeff(size_t imd, size_t ifq) const { return _cd1[imd][ifq]; }
+    sdpmComplex getUnsteadySideforceCoeff(size_t imd, size_t ifq) const { return _cs1[imd][ifq]; }
+    sdpmComplex getUnsteadyMomentCoeff(size_t imd, size_t ifq) const { return _cm1[imd][ifq]; }
+    void addMotion(double ux, double amp, double xr);
 #ifndef SWIG
     std::map<Node *, std::vector<Element *>> const &getMap() const
     {
         return _neMap;
     }
     Wake const &getWake() const { return *_wake; }
+    size_t getNumMotions() const { return _motions.size(); }
+    sdpmDouble getAmplitude(size_t imd) const { return _motions[imd]->getAmplitude(); }
+    std::vector<sdpmVector3cd> const &getSteadyVelocity(size_t imd) const { return _iVelocityS[imd]; };
+    std::vector<sdpmVector3cd> const &getHarmonicVelocity(size_t imd) const { return _iVelocityH[imd]; };
+    void setLoads(std::vector<sdpmVector3d> const &nl) { _nLoads = nl; }
+    void setUnsteadyLoads(size_t imd, size_t ifq, std::vector<sdpmVector3cd> const &nl) { _nLoads1[imd][ifq] = nl; }
+    void setLiftCoeff(sdpmDouble cl) { _cl = cl; }
+    void setDragCoeff(sdpmDouble cd) { _cd = cd; }
+    void setSideforceCoeff(sdpmDouble cs) { _cs = cs; }
+    void setMomentCoeff(sdpmDouble cm) { _cm = cm; }
+    void setUnsteadyLiftCoeff(size_t imd, size_t ifq, sdpmComplex cl) { _cl1[imd][ifq] = cl; }
+    void setUnsteadyDragCoeff(size_t imd, size_t ifq, sdpmComplex cd) { _cd1[imd][ifq] = cd; }
+    void setUnsteadySideforceCoeff(size_t imd, size_t ifq, sdpmComplex cs) { _cs1[imd][ifq] = cs; }
+    void setUnsteadyMomentCoeff(size_t imd, size_t ifq, sdpmComplex cm) { _cm1[imd][ifq] = cm; }
+    void computeVelocity(size_t imd);
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/sdpm/src/sdpmBuilder.cpp b/sdpm/src/sdpmBuilder.cpp
index 367ee52547ea25be3548c072cf6edde996910eef..d720fcca8177c3d1f2a4e9991b587e9c86cae5fe 100644
--- a/sdpm/src/sdpmBuilder.cpp
+++ b/sdpm/src/sdpmBuilder.cpp
@@ -26,10 +26,10 @@ using namespace sdpm;
 sdpmMatrix3d Builder::computeTransfoMat(bool symy, Element const *ej)
 {
     // Compute tangents and normal
-    sdpmVector3d u = (ej->getNodes()[0]->getCoords() - ej->getCg()).normalized();
-    sdpmVector3d v_ = ej->getNodes()[1]->getCoords() - ej->getCg();
+    sdpmVector3d u = (ej->getCompressibleCoords()[0] - ej->getCompressibleCg()).normalized();
+    sdpmVector3d v_ = ej->getCompressibleCoords()[1] - ej->getCompressibleCg();
     sdpmVector3d v = (v_ - (u.dot(v_)) * u).normalized();
-    sdpmVector3d n = ej->getNormal();
+    sdpmVector3d n = ej->getCompressibleNormal();
 
     // Assemble into transformation matrix
     sdpmMatrix3d g2l;
@@ -49,32 +49,43 @@ sdpmMatrix3d Builder::computeTransfoMat(bool symy, Element const *ej)
  */
 sdpmMatrixXd Builder::gatherCoords(bool symy, Element const *ej)
 {
-    // Get nodes
-    std::vector<Node *> nods = ej->getNodes();
-    size_t nV = nods.size();
+    // Get compressible coordinates
+    std::vector<sdpmVector3d> ccoords = ej->getCompressibleCoords();
+    size_t nV = ccoords.size();
 
     // Get coordinates
     Eigen::MatrixXd coords(3, nV);
     for (size_t j = 0; j < nV; ++j)
-        coords.col(j) = nods[nV - 1 - j]->getCoords();
+        coords.col(j) = ccoords[nV - 1 - j];
     if (symy)
         coords.row(1) = -coords.row(1);
     return coords;
 }
 
+/**
+ * @brief Compute vector joining body panel ej to body panel ei
+ */
+sdpmVector3d Builder::computeDistVec(bool symy, Element const *ei, Element const *ej)
+{
+    sdpmVector3d cgSrc = ej->getCompressibleCg();
+    sdpmVector3d cgTgt = ei->getCompressibleCg();
+    sdpmVector3d dist = cgTgt - cgSrc;
+    if (symy)
+        dist(1) = cgTgt(1) + cgSrc(1);
+    return dist;
+}
+
 /**
  * @brief Compute doublet AIC and source AIC from body panel ej to body panel ei
+ *
+ * Katz and Plotkin, Low-speed aerodynamics, 2001
+ * Hess and Smith, Calculation of potential flow about arbitrary bodies, 1967
  */
 std::vector<sdpmDouble> Builder::computeAIC(bool symy, Element const *ei, Element const *ej, Eigen::MatrixXd const &sCoords, sdpmMatrix3d const &g2l)
 {
     // Compute distance between panels and target panel coordinates in local axes
-    sdpmVector3d cgSrc = ej->getCg();
-    sdpmVector3d cgTgt = ei->getCg();
-    sdpmVector3d distance = cgTgt - cgSrc;
-    if (symy)
-        distance(1) = cgTgt(1) + cgSrc(1);
-    distance = g2l * distance;
-    sdpmVector3d tCoords = g2l * ei->getCg();
+    sdpmVector3d distance = g2l * computeDistVec(symy, ei, ej);
+    sdpmVector3d tCoords = g2l * ei->getCompressibleCg();
 
     // Compute coefficients
     int nV = static_cast<int>(ej->getNodes().size());
diff --git a/sdpm/src/sdpmBuilder.h b/sdpm/src/sdpmBuilder.h
index d56934355477a4e04d3f9f3881fdb5ef45b65e14..a816fd08cc026e5f8323d4969d3cce17b0bf0bc4 100644
--- a/sdpm/src/sdpmBuilder.h
+++ b/sdpm/src/sdpmBuilder.h
@@ -33,6 +33,7 @@ class SDPM_API Builder
 public:
     static sdpmMatrix3d computeTransfoMat(bool symy, Element const *ej);
     static sdpmMatrixXd gatherCoords(bool symy, Element const *ej);
+    static sdpmVector3d computeDistVec(bool symy, Element const *ei, Element const *ej);
     static std::vector<sdpmDouble> computeAIC(bool symy, Element const *ei, Element const *ej, sdpmMatrixXd const &trsfCoordSrc, sdpmMatrix3d const &g2l);
 };
 
diff --git a/sdpm/src/sdpmElement.cpp b/sdpm/src/sdpmElement.cpp
index 53b092f9f54ae5224bdc6b4b1061659299f98df3..06bd015992177b48cb3dffbb7f601870b0fadca5 100644
--- a/sdpm/src/sdpmElement.cpp
+++ b/sdpm/src/sdpmElement.cpp
@@ -38,20 +38,24 @@ SDPM_API std::ostream &operator<<(std::ostream &out, ElType const &obj)
 }
 } // namespace sdpm
 
-Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes)
+Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes), _toDate(false)
 {
+    _coords.resize(_nodes.size());
     _tag->addElements(this);
 }
 
 /**
  * @brief Update precomputed values
  */
-void Element::update()
+void Element::update(sdpmDouble beta)
 {
-    // Center of gravity
-    _cg = Eigen::Vector3d::Zero();
-    for (auto n : _nodes)
-        _cg += n->getCoords();
+    // Scale nodes and compute center of gravity
+    _cg = sdpmVector3d::Zero();
+    for (size_t i = 0; i < _nodes.size(); ++i)
+    {
+        _coords[i] = _nodes[i]->getCoords().cwiseQuotient(sdpmVector3d(beta, 1., 1.));
+        _cg += _coords[i];
+    }
     _cg /= static_cast<double>(_nodes.size());
 }
 
diff --git a/sdpm/src/sdpmElement.h b/sdpm/src/sdpmElement.h
index 938aad731823b74008115e8a0581b36114d3dc6e..a42c4105e8a462664721386da86dcb4585495b3f 100644
--- a/sdpm/src/sdpmElement.h
+++ b/sdpm/src/sdpmElement.h
@@ -51,10 +51,15 @@ protected:
     size_t _id;                 ///< element number
     Tag *_tag;                  ///< physical tag to which element belongs
     std::vector<Node *> _nodes; ///< nodes of the element
+    bool _toDate;               ///< whether the element values are up-to-date
     // Geometry
-    sdpmVector3d _cg;  ///< element center of gravity
-    sdpmDouble _vol;   ///< element volume
-    sdpmVector3d _nrm; ///< unit normal vector
+    std::vector<sdpmVector3d> _coords; ///< scaled (compressible) coordinates
+    sdpmVector3d _cg0;                 ///< element center of gravity
+    sdpmDouble _srf0;                  ///< element surface area
+    sdpmVector3d _nrm0;                ///< unit normal vector
+    sdpmVector3d _cg;                  ///< compressible element center of gravity
+    sdpmDouble _srf;                   ///< compressible element surface area
+    sdpmVector3d _nrm;                 ///< compressible unit normal vector
 
 public:
 #ifndef SWIG
@@ -69,7 +74,7 @@ public:
 
 #ifndef SWIG
     // Management
-    virtual void update();
+    virtual void update(sdpmDouble beta);
     // Setters
     void setNode(size_t i, Node *n) { _nodes[i] = n; }
     // Getters
@@ -77,15 +82,23 @@ public:
     size_t const &getIdRef() const { return _id; }
     Tag const *getTag() const { return _tag; }
     std::vector<Node *> const &getNodes() const { return _nodes; }
-    sdpmVector3d const &getCg() const { return _cg; }
-    sdpmDouble getVolume() const { return _vol; }
-    sdpmVector3d const &getNormal() const { return _nrm; }
+    sdpmVector3d const &getCg() const { return _cg0; }
+    sdpmDouble getSurfaceArea() const { return _srf0; }
+    sdpmVector3d const &getNormal() const { return _nrm0; }
+    inline std::vector<sdpmVector3d> const &getCompressibleCoords() const;
+    inline sdpmVector3d const &getCompressibleCg() const;
+    inline sdpmDouble getCompressibleSurfaceArea() const;
+    inline sdpmVector3d const &getCompressibleNormal() const;
     // Utilities
     virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const;
     virtual void write(std::ostream &out) const override;
 #endif
 };
 
+#ifndef SWIG
+#include "sdpmElement.inl"
+#endif
+
 } // namespace sdpm
 
 #endif // SDPMELEMENT_H
diff --git a/sdpm/src/sdpmElement.inl b/sdpm/src/sdpmElement.inl
new file mode 100644
index 0000000000000000000000000000000000000000..33ef73bd7f1b157185ff860ff52e2d097cc6171d
--- /dev/null
+++ b/sdpm/src/sdpmElement.inl
@@ -0,0 +1,59 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * @brief Return the compressible coordinates
+ */
+inline std::vector<sdpmVector3d> const &Element::getCompressibleCoords() const
+{
+    if (_toDate)
+        return _coords;
+    else
+        throw std::runtime_error("Element::getCompressibleCoords is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the center of gravity (compressible coordinates)
+ */
+inline sdpmVector3d const &Element::getCompressibleCg() const
+{
+    if (_toDate)
+        return _cg;
+    else
+        throw std::runtime_error("Element::getCompressibleCg is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the element volume (compressible coordinates)
+ */
+inline sdpmDouble Element::getCompressibleSurfaceArea() const
+{
+    if (_toDate)
+        return _srf;
+    else
+        throw std::runtime_error("Element::getCompressibleSurfaceArea is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the unit normal vector (compressible coordinates)
+ */
+inline sdpmVector3d const &Element::getCompressibleNormal() const
+{
+    if (_toDate)
+        return _nrm;
+    else
+        throw std::runtime_error("Element::getCompressibleNormal is invalid because Element is not up-to-date!\n");
+}
\ No newline at end of file
diff --git a/sdpm/src/sdpmLine2.cpp b/sdpm/src/sdpmLine2.cpp
index 21f7b6bda853a7739d0b2481f508779c725d24eb..efa473a6c749461286bf2efa9bd6034c266dcc11 100644
--- a/sdpm/src/sdpmLine2.cpp
+++ b/sdpm/src/sdpmLine2.cpp
@@ -20,15 +20,21 @@ using namespace sdpm;
 
 Line2::Line2(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
 {
-    update();
+    // Compute unscaled values
+    update(1.);
+    _cg0 = _cg;
+    _srf0 = _srf;
+    _nrm0 = _nrm;
+    _toDate = false;
 }
 
 /**
  * @brief Update precomputed values and gradients
  */
-void Line2::update()
+void Line2::update(double beta)
 {
-    Element::update();                                               // cg
-    _vol = (_nodes[1]->getCoords() - _nodes[0]->getCoords()).norm(); // volume
-    _nrm = sdpmVector3d::Zero();                                     // unit normal
+    Element::update(beta);                   // cg
+    _srf = (_coords[1] - _coords[0]).norm(); // volume
+    _nrm = sdpmVector3d::Zero();             // unit normal
+    _toDate = true;
 }
diff --git a/sdpm/src/sdpmLine2.h b/sdpm/src/sdpmLine2.h
index 8eea05ec80a8b98099f460f7cd4e52057314ee9e..1022666d7a3b531f096fa9dede491220749035b2 100644
--- a/sdpm/src/sdpmLine2.h
+++ b/sdpm/src/sdpmLine2.h
@@ -35,7 +35,7 @@ public:
 
     virtual ElType getType() const override { return ElType::LINE2; }
 #ifndef SWIG
-    virtual void update() override;
+    virtual void update(double beta) override;
 #endif
 };
 
diff --git a/sdpm/src/sdpmMotion.cpp b/sdpm/src/sdpmMotion.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1339b924f0eb0a06fe3fbf22e507c844d7563fcb
--- /dev/null
+++ b/sdpm/src/sdpmMotion.cpp
@@ -0,0 +1,43 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "sdpm.h"
+#include "sdpmMotion.h"
+#include "sdpmElement.h"
+
+using namespace sdpm;
+
+Motion::Motion(sdpmDouble xvel, sdpmDouble amp, sdpmDouble xref) : sdpmObject(),
+                                                                   _xVel(xvel), _amp(amp), _xRef(xref)
+{
+}
+
+/**
+ * @brief Compute the steady part of the motion induced velocity on an element of the body
+ */
+sdpmVector3cd Motion::computeSteady(Element const &e)
+{
+    return sdpmVector3cd(0., 0., _xVel * _amp);
+}
+
+/**
+ * @brief Compute the harmonic part of the motion induced velocity on an element of the body
+ */
+sdpmVector3cd Motion::computeHarmonic(Element const &e)
+{
+    sdpmVector3d cg = e.getCompressibleCg();
+    return sdpmVector3cd(-_amp * cg(2), 0., _amp * (cg(0) - _xRef));
+}
diff --git a/sdpm/src/sdpmMotion.h b/sdpm/src/sdpmMotion.h
new file mode 100644
index 0000000000000000000000000000000000000000..3973d4390d48d62680ebc1c35d46fd3bff1502ba
--- /dev/null
+++ b/sdpm/src/sdpmMotion.h
@@ -0,0 +1,48 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef SDPMMOTION_H
+#define SDPMMOTION_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+
+namespace sdpm
+{
+
+/**
+ * @brief Hold the parameters defning the motion of a body and compute the induced velocity
+ * @authors Adrien Crovato
+ * @todo support several types of motion and combinations
+ */
+class SDPM_API Motion : public sdpmObject
+{
+    sdpmDouble _xVel;
+    sdpmDouble _amp;
+    sdpmDouble _xRef;
+
+public:
+    Motion(sdpmDouble xvel, sdpmDouble amp, sdpmDouble xref);
+    ~Motion(){};
+
+    sdpmDouble getAmplitude() const { return _amp; }
+    sdpmVector3cd computeSteady(Element const &e);
+    sdpmVector3cd computeHarmonic(Element const &e);
+};
+
+} // namespace sdpm
+
+#endif // SDPMMOTION_H
diff --git a/sdpm/src/sdpmProblem.cpp b/sdpm/src/sdpmProblem.cpp
index ad141a3006e2da1070c1fe4a9f0f6b5ea831b6da..50fcb8f1c58fc77187f1f78069429dbce10b03fe 100644
--- a/sdpm/src/sdpmProblem.cpp
+++ b/sdpm/src/sdpmProblem.cpp
@@ -21,14 +21,13 @@
 #include <iostream>
 using namespace sdpm;
 
-Problem::Problem(Mesh &msh, double aoa, double aos, double mch,
-                 double sref, double cref, double xref, double yref, double zref,
-                 bool ysym) : _msh(msh), _mach(mch), _sRef(sref), _cRef(cref), _ySym(ysym)
+Problem::Problem(Mesh &msh, size_t nmods, std::vector<double> const &fs) : _msh(msh),
+                                                                           _aoa(0.), _aos(0.), _mach(0.), _beta(1.0), _sound(0.), _uinf(sdpmVector3d::Zero()),
+                                                                           _dDir(sdpmVector3d::Zero()), _sDir(sdpmVector3d::Zero()), _lDir(sdpmVector3d::Zero()),
+                                                                           _sRef(0.), _cRef(0.), _xRef(sdpmVector3d::Zero()), _ySym(false),
+                                                                           _nModes(nmods), _omega(fs),
+                                                                           _hasAero(false), _hasGeo(false)
 {
-    updateFreestream(aoa, aos);
-    _xRef(0) = xref;
-    _xRef(1) = yref;
-    _xRef(2) = zref;
 }
 
 Problem::~Problem()
@@ -45,18 +44,38 @@ void Problem::addBody(Body &b)
 }
 
 /**
- * @brief Update the freestream values
+ * @brief Update the freestream aerodynamic values
  */
-void Problem::updateFreestream(double aoa, double aos)
+void Problem::updateFreestream(double aoa, double aos, double mch, double ui)
 {
-    // Angles and velocity
+    // Angles
     _aoa = aoa;
     _aos = aos;
-    _uinf = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
-    // Directions
+    // Compressiblity
+    _mach = mch;
+    _beta = std::sqrt(1 - _mach * _mach);
+    _sound = ui / _mach;
+    // Velocity and directions
+    _uinf = ui * sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
     _dDir = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
     _sDir = sdpmVector3d(-cos(_aoa) * sin(_aos), cos(_aos), -sin(_aoa) * sin(_aos));
-    _lDir = sdpmVector3d(-sin(_aoa), 0, cos(_aoa));
+    _lDir = sdpmVector3d(-sin(_aoa), 0., cos(_aoa));
+    // Flag
+    _hasAero = true;
+}
+
+/**
+ * @brief Set the reference geometric values
+ */
+void Problem::setGeometry(double sref, double cref, double xref, double yref, double zref, bool ysym)
+{
+    // Reference values
+    _sRef = sref;
+    _cRef = cref;
+    _xRef = sdpmVector3d(xref, yref, zref);
+    _ySym = ysym;
+    // Flag
+    _hasGeo = true;
 }
 
 /**
@@ -64,13 +83,12 @@ void Problem::updateFreestream(double aoa, double aos)
  */
 void Problem::updateElements()
 {
-    // Update surface Jacobian
     for (auto b : _bodies)
     {
         for (auto e : b->getElements())
-            e->update();
+            e->update(_beta);
         for (auto e : b->getWake().getElements())
-            e->update();
+            e->update(_beta);
     }
 }
 
@@ -80,10 +98,15 @@ void Problem::updateElements()
 void Problem::check() const
 {
     // Sanity checks
+    if (!_hasAero)
+        throw std::runtime_error("sdpm::Problem: aerodynamic parameters have not been provided!\n");
+    if (!_hasGeo)
+        throw std::runtime_error("sdpm::Problem: geometric parameters have not been provided!\n");
     if (_bodies.empty())
         throw std::runtime_error("sdpm::Problem: no lifting bodies provided!\n");
-    if (_mach != 0)
-        throw std::runtime_error("sdpm::Problem: compressible computation feature not implemented yet!\n");
+    for (auto b : _bodies)
+        if (b->getNumMotions() != _nModes)
+            throw std::runtime_error("sdpm::Problem: Body and Problem do not have the same number of modes!\n");
 
     // Surface elements
     for (auto b : _bodies)
@@ -102,6 +125,7 @@ void Problem::write(std::ostream &out) const
     out << "sdpm::Problem definition"
         << "\n\tAngle of attack: " << _aoa
         << "\n\tAngle of sideslip: " << _aos
+        << "\n\tMach number: " << _mach
         << "\n\tReference area: " << _sRef
         << "\n\tReference length: " << _cRef
         << "\n\tReference point: [" << _xRef.transpose() << "]"
diff --git a/sdpm/src/sdpmProblem.h b/sdpm/src/sdpmProblem.h
index 4a2dac8afd76724a7468e8348c09d4297dc1fd05..bee2d6f505ceb6871b8061fca84f9bcc6856f847 100644
--- a/sdpm/src/sdpmProblem.h
+++ b/sdpm/src/sdpmProblem.h
@@ -35,22 +35,29 @@ class SDPM_API Problem : public sdpmObject
     Mesh &_msh;                  ///< mesh data structure
     std::vector<Body *> _bodies; ///< bodies in fluid
     // Aerodynamic
-    sdpmDouble _aoa; ///< Angle of attack
-    sdpmDouble _aos; ///< Angle of sideslip
-    double _mach;    ///< Mach number
+    sdpmDouble _aoa;    ///< Angle of attack
+    sdpmDouble _aos;    ///< Angle of sideslip
+    sdpmDouble _mach;   ///< Mach number
+    sdpmDouble _beta;   ///< Compressibility factor
+    sdpmDouble _sound;  ///< Speed of sound
+    sdpmVector3d _uinf; ///< Freestream velocity
+    sdpmVector3d _dDir; ///< Direction of drag force
+    sdpmVector3d _sDir; ///< Direction of side force
+    sdpmVector3d _lDir; ///< Direction of lift force
     // Geometry
-    double _sRef;          ///< Reference surface area
-    double _cRef;          ///< Reference chord length
-    Eigen::Vector3d _xRef; ///< Reference center point for moment computation
-    bool _ySym;            ///< True if mesh is symmetric with respect to y-plane
-    // Freestream
-    sdpmVector3d _uinf; ///< freestream velocity
-    sdpmVector3d _dDir; ///< direction of drag force
-    sdpmVector3d _sDir; ///< direction of side force
-    sdpmVector3d _lDir; ///< direction of lift force
+    sdpmDouble _sRef;   ///< Reference surface area
+    sdpmDouble _cRef;   ///< Reference chord length
+    sdpmVector3d _xRef; ///< Reference center point for moment computation
+    bool _ySym;         ///< True if mesh is symmetric with respect to y-plane
+    // Unsteady
+    size_t _nModes;             ///< Number of modal motions
+    std::vector<double> _omega; ///< Circular frequencies
+    // Check flags
+    bool _hasAero; ///< True if aerodynamic parameters have been updated
+    bool _hasGeo;  ///< True if geometric parameters have been updated
 
 public:
-    Problem(Mesh &msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref, bool ysym);
+    Problem(Mesh &msh, size_t nmods = 0, std::vector<double> const &fs = {});
     ~Problem();
 
 #ifndef SWIG
@@ -62,17 +69,23 @@ public:
     std::vector<Element *> const &getElements() const { return _msh.getElements(); }
     sdpmDouble getAoa() const { return _aoa; }
     sdpmDouble getAos() const { return _aos; }
+    sdpmDouble getMach() const { return _mach; }
+    sdpmDouble getCompressibility() const { return _beta; }
+    sdpmDouble getSpeedOfSound() const { return _sound; }
     sdpmDouble getRefArea() const { return _sRef; }
     sdpmDouble getRefLgth() const { return _cRef; }
     sdpmVector3d const &getRefCtr() const { return _xRef; }
     sdpmDouble getSym() const { return _ySym; }
+    std::vector<double> const &getFreqs() const { return _omega; }
+    size_t getNumModes() const { return _nModes; }
     sdpmVector3d const &getFrstrmVel() const { return _uinf; }
     sdpmVector3d const &getDragDir() const { return _dDir; }
     sdpmVector3d const &getSideDir() const { return _sDir; }
     sdpmVector3d const &getLiftDir() const { return _lDir; }
 #endif
     void addBody(Body &b);
-    void updateFreestream(double aoa, double aos);
+    void updateFreestream(double aoa, double aos, double mch, double ui);
+    void setGeometry(double sref, double cref, double xref, double yref, double zref, bool ysym);
 
 #ifndef SWIG
     void check() const;
diff --git a/sdpm/src/sdpmQuad4.cpp b/sdpm/src/sdpmQuad4.cpp
index 8281830ea14be254ce467522ea074efaa64dd96f..3425915fbbb280a79a1dd9cbaf0e252eee6494bb 100644
--- a/sdpm/src/sdpmQuad4.cpp
+++ b/sdpm/src/sdpmQuad4.cpp
@@ -20,19 +20,25 @@ using namespace sdpm;
 
 Quad4::Quad4(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
 {
-    update();
+    // Compute unscaled values
+    update(1.);
+    _cg0 = _cg;
+    _srf0 = _srf;
+    _nrm0 = _nrm;
+    _toDate = false;
 }
 
 /**
  * @brief Update precomputed values
  */
-void Quad4::update()
+void Quad4::update(double beta)
 {
-    Element::update();                                       // cg
-    _tgt0 = _nodes[2]->getCoords() - _nodes[0]->getCoords(); // unit tangent (n0->n2)
-    _tgt1 = _nodes[3]->getCoords() - _nodes[1]->getCoords(); // unit tangent (n1->n4)
-    _vol = 0.5 * _tgt0.cross(_tgt1).norm();                  // volume
-    _nrm = _tgt0.cross(_tgt1).normalized();                  // unit normal
+    Element::update(beta);                  // cg
+    _tgt0 = _coords[2] - _coords[0];        // unit tangent (n1->n3)
+    _tgt1 = _coords[3] - _coords[1];        // unit tangent (n2->n4)
+    _srf = 0.5 * _tgt0.cross(_tgt1).norm(); // volume
+    _nrm = _tgt0.cross(_tgt1).normalized(); // unit normal
+    _toDate = true;
 }
 
 /**
@@ -57,14 +63,14 @@ sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
     sdpmMatrix3d jac = sdpmMatrix3d::Zero(), ijac = sdpmMatrix3d::Zero();
     // 1) compute first 2 rows: dXi, dEta (same as 2D, but with 3D coordinates)
     for (size_t i = 0; i < _nodes.size(); ++i)
-        jac.block(0, 0, 2, 3) += dsf.col(i) * _nodes[i]->getCoords().transpose();
+        jac.block(0, 0, 2, 3) += dsf.col(i) * _coords[i].transpose();
     // 2) compute thrid rows: dZeta = dXi X dEta
     jac.row(2) = jac.row(0).cross(jac.row(1));
     // 3) inverse Jacobian
     ijac = jac.inverse();
 
     // Get solution at element nodes
-    Eigen::VectorXd ue(_nodes.size());
+    sdpmVectorXd ue(_nodes.size());
     for (size_t i = 0; i < _nodes.size(); ++i)
         ue(i) = u[_nodes[i]->getId() - 1];
 
diff --git a/sdpm/src/sdpmQuad4.h b/sdpm/src/sdpmQuad4.h
index 4d9a948a80b7dfcf75ba1f346141c4b751a3041f..98ee4361c2185513b9d66564b98955a755d18e54 100644
--- a/sdpm/src/sdpmQuad4.h
+++ b/sdpm/src/sdpmQuad4.h
@@ -38,7 +38,7 @@ public:
 
     virtual ElType getType() const override { return ElType::QUAD4; }
 #ifndef SWIG
-    virtual void update() override;
+    virtual void update(double beta) override;
     virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const override;
 #endif
 };
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
index ca8e0bd57f635ed206c793abd7a4fe9623629ae3..e830ff15f927a04362a7f1b62761276fef4da3b4 100644
--- a/sdpm/src/sdpmSolver.cpp
+++ b/sdpm/src/sdpmSolver.cpp
@@ -21,6 +21,7 @@
 #include "sdpmMesh.h"
 #include "sdpmGmshExport.h"
 #include "sdpmResults.h"
+#include "sdpm_extras.h"
 #include <iostream>
 #include <iomanip>
 #include <fstream>
@@ -37,27 +38,41 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
     std::cout << "*******************************************************************************" << std::endl;
 
     // Check the problem
-    pbl.check();
+    _pbl.check();
 
-    // Count panel, set up global->local id map and resize AIC matrices
+    // Count panel and set up global->local id map
     _nP = 0;
     int i = 0;
-    for (auto body : pbl.getBodies())
+    for (auto body : _pbl.getBodies())
     {
         _nP += body->getElements().size();
         for (auto e : body->getElements())
             _rows[e] = i++;
     }
-    _A = sdpmMatrixXd::Zero(_nP, _nP);
-    _B = sdpmMatrixXd::Zero(_nP, _nP);
 
-    // Set up variables
+    // Get sizes
+    size_t nm = _pbl.getNumModes();
+    size_t nf = _pbl.getFreqs().size();
     size_t ne = _pbl.getElements().size();
+
+    // Resize AIC matrics
+    _A0 = sdpmMatrixXd(_nP, _nP);
+    _B0 = sdpmMatrixXd(_nP, _nP);
+    _A.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+    _B.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+
+    // Set up variables
     _phi.resize(ne, 0.);
     _tau.resize(ne, 0.);
     _mu.resize(ne, 0.);
     _u.resize(ne, sdpmVector3d::Zero());
     _cp.resize(ne, 0.);
+    _u1.resize(nm, std::vector<std::vector<sdpmVector3cd>>(nf, std::vector<sdpmVector3cd>(ne, sdpmVector3cd::Zero())));
+    _cp1.resize(nm, std::vector<std::vector<sdpmComplex>>(nf, std::vector<sdpmComplex>(ne, sdpmComplex(0., 0.))));
+    _cl1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cd1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cs1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+    _cm1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
 
     // Display configuration
     std::cout << "--- Physical groups ---\n";
@@ -67,12 +82,16 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
     std::cout << "--- Freestream conditions ---\n"
               << std::setprecision(5)
               << "Angle of attack: " << _pbl.getAoa() * 180 / 3.14159 << " deg\n"
-              << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n";
+              << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n"
+              << "Mach number: " << _pbl.getMach() << "\n"
+              << "Velocity: " << _pbl.getFrstrmVel().norm() << "\n";
     std::cout << "--- Reference data ---\n"
-              << "Reference area: " << _pbl.getRefArea() << "\n"
-              << "Reference length: " << _pbl.getRefLgth() << "\n"
-              << "Reference origin: (" << _pbl.getRefCtr().transpose() << ")"
-              << "\n"
+              << "Surface area: " << _pbl.getRefArea() << "\n"
+              << "Chord length: " << _pbl.getRefLgth() << "\n"
+              << "Origin: (" << _pbl.getRefCtr().transpose() << ")\n";
+    std::cout << "--- Unsteady parameters ---\n"
+              << "Number of frequencies: " << _pbl.getFreqs().size() << "\n"
+              << "Number of motions: " << _pbl.getNumModes() << "\n"
               << std::endl;
 }
 
@@ -86,34 +105,51 @@ Solver::~Solver()
  */
 void Solver::update()
 {
-    // Update elements
+    // Update elements and get values
     _pbl.updateElements();
+    std::vector<double> omega = _pbl.getFreqs();
+    sdpmDouble mach = _pbl.getMach();
+    sdpmDouble sound = _pbl.getSpeedOfSound();
+    sdpmDouble beta = _pbl.getCompressibility();
+    sdpmComplex img(0, 1);
 
     // Compute AIC
     std::cout << "Computing AIC matrices... " << std::flush;
     std::vector<Body *> bodies = _pbl.getBodies();
-    _A.setZero();
-    _B.setZero();
+    _A0.setZero();
+    _B0.setZero();
+    std::fill(_A.begin(), _A.end(), sdpmMatrixXcd::Zero(_nP, _nP));
+    std::fill(_B.begin(), _B.end(), sdpmMatrixXcd::Zero(_nP, _nP));
     for (auto jbody : bodies)
     {
         // body
         for (auto ej : jbody->getElements())
         {
+            sdpmDouble nx = ej->getCompressibleNormal()(0);                /// x-component of normal vector
             sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ej);      // global to local transformation matrix
             sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ej); // coordinates of source panel in local axes
             for (auto ibody : bodies)
             {
                 for (auto ei : ibody->getElements())
                 {
-                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l);
-                    _A(_rows.at(ei), _rows.at(ej)) = mutau[0];
-                    _B(_rows.at(ei), _rows.at(ej)) = mutau[1];
+                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l); // steady AIC
+                    _A0(_rows.at(ei), _rows.at(ej)) = mutau[0];
+                    _B0(_rows.at(ei), _rows.at(ej)) = mutau[1];
+                    for (size_t fk = 0; fk < omega.size(); ++fk)
+                    {
+                        sdpmDouble om = omega[fk] / sound / beta;
+                        sdpmVector3d r = Builder::computeDistVec(false, ei, ej);
+                        sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+                        _A[fk](_rows.at(ei), _rows.at(ej)) = (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+                        _B[fk](_rows.at(ei), _rows.at(ej)) = mutau[1] * exp;
+                    }
                 }
             }
         }
         // wake
         Wake const &wake = jbody->getWake();
-        std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+        std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+        std::map<Element *, double> lagMap = wake.getLagMap();
         for (auto ew : wake.getElements())
         {
             sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ew);      // global to local transformation matrix
@@ -122,9 +158,18 @@ void Solver::update()
             {
                 for (auto ei : ibody->getElements())
                 {
-                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l);
-                    _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
-                    _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+                    std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l); // steady AIC
+                    _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
+                    _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+                    for (size_t fk = 0; fk < omega.size(); ++fk)
+                    {
+                        sdpmDouble om = omega[fk] / sound / beta;
+                        sdpmVector3d r = Builder::computeDistVec(false, ei, ew);
+                        sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+                        sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+                        _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) += waktrm;
+                        _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= waktrm;
+                    }
                 }
             }
         }
@@ -137,21 +182,31 @@ void Solver::update()
             // body
             for (auto ej : jbody->getElements())
             {
+                sdpmDouble nx = ej->getCompressibleNormal()(0);               /// x-component of normal vector
                 sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ej);      // global to local transformation matrix
                 sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ej); // coordinates of source panel in local axes
                 for (auto ibody : bodies)
                 {
                     for (auto ei : ibody->getElements())
                     {
-                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l);
-                        _A(_rows.at(ei), _rows.at(ej)) -= mutau[0];
-                        _B(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l); // steady AIC
+                        _A0(_rows.at(ei), _rows.at(ej)) -= mutau[0];
+                        _B0(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+                        for (size_t fk = 0; fk < omega.size(); ++fk)
+                        {
+                            sdpmDouble om = omega[fk] / sound / beta;
+                            sdpmVector3d r = Builder::computeDistVec(true, ei, ej);
+                            sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+                            _A[fk](_rows.at(ei), _rows.at(ej)) -= (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+                            _B[fk](_rows.at(ei), _rows.at(ej)) -= mutau[1] * exp;
+                        }
                     }
                 }
             }
             // wake
             Wake const &wake = jbody->getWake();
-            std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+            std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+            std::map<Element *, double> lagMap = wake.getLagMap();
             for (auto ew : wake.getElements())
             {
                 sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ew);      // global to local transformation matrix
@@ -160,9 +215,18 @@ void Solver::update()
                 {
                     for (auto ei : ibody->getElements())
                     {
-                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l);
-                        _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
-                        _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+                        std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l); // steady AIC
+                        _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
+                        _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+                        for (size_t fk = 0; fk < omega.size(); ++fk)
+                        {
+                            sdpmDouble om = omega[fk] / sound / beta;
+                            sdpmVector3d r = Builder::computeDistVec(true, ei, ew);
+                            sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+                            sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+                            _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= waktrm;
+                            _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) += waktrm;
+                        }
                     }
                 }
             }
@@ -172,7 +236,7 @@ void Solver::update()
 }
 
 /**
- * @brief Solve the potential equation
+ * @brief Solve the steady potential equation
  */
 void Solver::run()
 {
@@ -182,21 +246,61 @@ void Solver::run()
     sdpmVectorXd etau(_nP);
     for (auto body : _pbl.getBodies())
         for (auto e : body->getElements())
-            etau(_rows[e]) = ui.dot(e->getNormal());
+            etau(_rows[e]) = ui.dot(e->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
     std::cout << "done." << std::endl;
 
     std::cout << "solving for doublets... " << std::flush;
     sdpmVectorXd emu(_nP);
-    emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A).solve(_B * etau);
+    emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A0).solve(_B0 * etau);
     std::cout << "done." << std::endl;
 
     std::cout << "computing flow and loads on bodies... " << std::flush;
-    computeFlow(etau, emu);
-    computeLoad();
+    post(etau, emu);
     std::cout << "done." << std::endl
               << std::endl;
 }
 
+/**
+ * @brief Solve the unsteady potential equation
+ * @todo remove output
+ */
+void Solver::runUnsteady()
+{
+    std::cout << "Running unsteady SDPM solver... " << std::endl;
+    for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
+    {
+        // Compute motion-induced velocity
+        for (auto body : _pbl.getBodies())
+            body->computeVelocity(imd);
+        // Solve
+        for (size_t ifq = 0; ifq < _pbl.getFreqs().size(); ++ifq)
+        {
+            std::cout << "Mode #" << imd << ", frequency #" << ifq << std::endl;
+            std::cout << "computing sources... " << std::flush;
+            sdpmVectorXcd etau(_nP);
+            for (auto body : _pbl.getBodies())
+            {
+                std::vector<sdpmVector3cd> uc0 = body->getSteadyVelocity(imd);
+                std::vector<sdpmVector3cd> uc1 = body->getHarmonicVelocity(imd);
+                std::vector<Element *> elems = body->getElements();
+                for (size_t i = 0; i < elems.size(); ++i)
+                    etau(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFreqs()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+            }
+            std::cout << "done." << std::endl;
+
+            std::cout << "solving for doublets... " << std::flush;
+            sdpmVectorXcd emu(_nP);
+            emu = Eigen::PartialPivLU<sdpmMatrixXcd>(_A[ifq]).solve(_B[ifq] * etau);
+            std::cout << "done." << std::endl;
+
+            std::cout << "computing flow and loads on bodies... " << std::flush;
+            postUnsteady(imd, ifq, etau, emu);
+            std::cout << "done." << std::endl;
+        }
+    }
+    std::cout << std::endl;
+}
+
 /**
  * @brief Write the results
  */
@@ -221,10 +325,10 @@ void Solver::save(GmshExport const &writer, std::string const &suffix)
     {
         outfile << "# Body - " << b->getName() << " (" << b->getElements().size() << " elements)" << std::endl;
         outfile << std::fixed << std::setprecision(5)
-                << "Cl = " << b->cl << "\n"
-                << "Cd = " << b->cd << "\n"
-                << "Cs = " << b->cs << "\n"
-                << "Cm = " << b->cm << std::endl
+                << "Cl = " << b->getLiftCoeff() << "\n"
+                << "Cd = " << b->getDragCoeff() << "\n"
+                << "Cs = " << b->getSideforceCoeff() << "\n"
+                << "Cm = " << b->getMomentCoeff() << std::endl
                 << std::endl;
     }
     outfile.close();
@@ -232,12 +336,12 @@ void Solver::save(GmshExport const &writer, std::string const &suffix)
 }
 
 /**
- * @brief Compute flow solution on bodies
+ * @brief Compute steady flow functionals and loads on bodies
  */
-void Solver::computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
+void Solver::post(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
 {
     // Store source and doublet strength and compute perturbation potential on elements
-    sdpmVectorXd ephi = -_A * emu - _B * etau;
+    sdpmVectorXd ephi = -_A0 * emu - _B0 * etau;
     std::vector<Body *> bodies = _pbl.getBodies();
     for (auto body : bodies)
     {
@@ -255,35 +359,27 @@ void Solver::computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
     std::vector<sdpmDouble> muNodes(_pbl.getNodes().size(), 0);
     for (auto b : bodies)
     {
-        // associate the value to its element
-        std::map<Element *, sdpmDouble> mapMu;
-        for (auto e : b->getElements())
-            mapMu[e] = _mu[e->getId() - 1];
-        // average at nodes
         std::map<Node *, std::vector<Element *>> neMap = b->getMap();
         for (auto n : b->getNodes())
             for (auto e : neMap.at(n))
-                muNodes[n->getId() - 1] += mapMu.at(e) / neMap.at(n).size();
+                muNodes[n->getId() - 1] += _mu[e->getId() - 1] / neMap.at(n).size();
     }
 
     // Compute surface flow functionals
     sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
+    sdpmDouble mi = _pbl.getMach();        // freestream speed of sound
     for (auto body : bodies)
     {
         for (auto e : body->getElements())
         {
-            sdpmVector3d u = ui - _tau[e->getId() - 1] * e->getNormal() - e->computeGradient(muNodes);
-            _u[e->getId() - 1] = u;             // velocity
-            _cp[e->getId() - 1] = 1 - u.dot(u); // pressure coefficient
+            sdpmVector3d u = ui - (_tau[e->getId() - 1] * e->getCompressibleNormal() + e->computeGradient(muNodes)).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+            sdpmDouble cpl = 1 - u.dot(u) / ui.dot(ui);
+            _u[e->getId() - 1] = u;                                      // velocity
+            _cp[e->getId() - 1] = cpl + 0.5 * 0.5 * cpl * cpl * mi * mi; // nonlinear pressure coefficient
         }
     }
-}
 
-/**
- * @brief Compute aerodynamic load coefficients on bodies
- */
-void Solver::computeLoad()
-{
+    // Compute aerodynamic loads
     sdpmDouble sRef = _pbl.getRefArea();      // reference surface area
     sdpmDouble cRef = _pbl.getRefLgth();      // reference length
     sdpmVector3d xRef = _pbl.getRefCtr();     // reference center
@@ -296,36 +392,140 @@ void Solver::computeLoad()
     _cm = 0;
     for (auto b : _pbl.getBodies())
     {
-        // Reset
-        std::fill(b->nLoads.begin(), b->nLoads.end(), sdpmVector3d::Zero());
-
         // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element and average at node
         std::vector<Node *> nods = b->getNodes();
         std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        std::vector<sdpmVector3d> nLoads(nods.size(), sdpmVector3d::Zero());
         for (size_t i = 0; i < nods.size(); ++i)
             for (auto e : neMap.at(nods[i]))
-                b->nLoads[i] += -_cp[e->getId() - 1] * e->getVolume() * e->getNormal() / e->getNodes().size();
+                nLoads[i] += -_cp[e->getId() - 1] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+        b->setLoads(nLoads);
 
         // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
-        // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
-        b->cm = 0;
+        // force coefficients along x, y and z directions and pitching moment around y axis
+        sdpmDouble bcm = 0;
         sdpmVector3d cf(0, 0, 0);
         for (size_t i = 0; i < nods.size(); ++i)
         {
             sdpmVector3d l = nods[i]->getCoords() - xRef; // lever arm
-            sdpmVector3d cfi = b->nLoads[i];              // normalized force
+            sdpmVector3d cfi = nLoads[i];                 // normalized force
             cf += cfi;
-            b->cm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+            bcm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
         }
+        b->setMomentCoeff(bcm);
         // rotate to flow direction
-        b->cd = cf.dot(dragDir) / sRef;
-        b->cs = cf.dot(sideDir) / sRef;
-        b->cl = cf.dot(liftDir) / sRef;
+        b->setDragCoeff(cf.dot(dragDir) / sRef);
+        b->setSideforceCoeff(cf.dot(sideDir) / sRef);
+        b->setLiftCoeff(cf.dot(liftDir) / sRef);
+        // compute total
+        _cl += b->getLiftCoeff();
+        _cd += b->getDragCoeff();
+        _cs += b->getSideforceCoeff();
+        _cm += bcm;
+    }
+}
+
+/**
+ * @brief Compute unsteady flow functionals and loads on bodies
+ */
+void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu)
+{
+    // Transfer doublet at nodes
+    std::vector<sdpmDouble> muNodesRe(_pbl.getNodes().size(), 0), muNodesIm(_pbl.getNodes().size(), 0);
+    for (auto b : _pbl.getBodies())
+    {
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        for (auto n : b->getNodes())
+        {
+            for (auto e : neMap.at(n))
+            {
+                muNodesRe[n->getId() - 1] += emu(_rows.at(e)).real() / neMap.at(n).size();
+                muNodesIm[n->getId() - 1] += emu(_rows.at(e)).imag() / neMap.at(n).size();
+            }
+        }
+    }
+
+    // Compute flow and loads
+    sdpmVector3d ui = _pbl.getFrstrmVel();
+    double omega = _pbl.getFreqs()[ifq];
+    sdpmDouble sRef = _pbl.getRefArea();
+    sdpmDouble cRef = _pbl.getRefLgth();
+    sdpmVector3d xRef = _pbl.getRefCtr();
+    double ca = cos(_pbl.getAos());
+    double sa = sin(_pbl.getAos());
+    _cl1[imd][ifq] = sdpmComplex(0., 0.);
+    _cd1[imd][ifq] = sdpmComplex(0., 0.);
+    _cs1[imd][ifq] = sdpmComplex(0., 0.);
+    _cm1[imd][ifq] = sdpmComplex(0., 0.);
+    for (auto b : _pbl.getBodies())
+    {
+        // compute velocity
+        std::vector<sdpmVector3cd> uc0 = b->getSteadyVelocity(imd);
+        std::vector<sdpmVector3cd> uc1 = b->getHarmonicVelocity(imd);
+        std::vector<Element *> elems = b->getElements();
+        for (size_t i = 0; i < elems.size(); ++i)
+            _u1[imd][ifq][elems[i]->getId() - 1] = uc0[i] + uc1[i] * sdpmComplex(0., 1.) * omega - (etau[_rows[elems[i]]] * elems[i]->getCompressibleNormal() + (elems[i]->computeGradient(muNodesRe) + sdpmComplex(0, 1) * elems[i]->computeGradient(muNodesIm))).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+        // compute pressure
+        std::vector<std::vector<sdpmComplex>> cp012(elems.size(), std::vector<sdpmComplex>(3));
+        for (size_t i = 0; i < elems.size(); ++i)
+        {
+            // si = phi_t + 1/2 conv(u_k,u_k)_k - 1/2 u_inf^2
+            sdpmVectorXcd si(5);
+            si << 0., std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), -0.5 * ui.dot(ui), -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]), 0.;
+            for (size_t j = 0; j < 3; ++j)
+            {
+                sdpmVector3cd utrm(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(j), _u[elems[i]->getId() - 1](j), _u1[imd][ifq][elems[i]->getId() - 1](j));
+                si += 0.5 * convolve(utrm, utrm);
+            }
+            // ecp = -2 / u_inf^2 * si + 1/u_inf^2/a_inf^2 * conv(si, si)
+            sdpmVectorXcd ecp(9);
+            ecp << 0., 0., -2 * si / ui.dot(ui), 0., 0.;
+            ecp += convolve(si, si) / (_pbl.getSpeedOfSound() * _pbl.getSpeedOfSound() * ui.dot(ui));
+            for (size_t u = 0; u < 3; ++u)
+                cp012[i][u] = ecp(4 + u);
+            _cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
+        }
+
+        // compute loads on nodes
+        std::vector<Node *> nods = b->getNodes();
+        std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+        std::vector<std::vector<sdpmVector3cd>> nLoads012(nods.size(), std::vector<sdpmVector3cd>(3, sdpmVector3cd::Zero()));
+        for (size_t i = 0; i < nods.size(); ++i)
+            for (auto e : neMap.at(nods[i]))
+                for (size_t u = 0; u < 3; ++u)
+                    nLoads012[i][u] += -cp012[_rows.at(e)][u] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+        std::vector<sdpmVector3cd> nLoads1(nods.size());
+        for (size_t i = 0; i < nods.size(); ++i)
+            nLoads1[i] = nLoads012[i][1];
+        b->setUnsteadyLoads(imd, ifq, nLoads1);
+
+        // compute moment coefficient and total load on body
+        sdpmComplex bcm1(0., 0.);
+        std::vector<sdpmVector3cd> cf(3, sdpmVector3cd::Zero());
+        for (size_t i = 0; i < nods.size(); ++i)
+        {
+            sdpmVector3d l = nods[i]->getCoords() - xRef;
+            for (size_t u = 0; u < 3; ++u)
+                cf[u] += nLoads012[i][u];
+            bcm1 += (l(2) * nLoads012[i][1](0) - l(0) * nLoads012[i][1](2)) / (sRef * cRef);
+        }
+        b->setUnsteadyMomentCoeff(imd, ifq, bcm1);
+
+        // compute load coefficients, only valid for small AoAs
+        sdpmVector3cd alpha(b->getAmplitude(imd), _pbl.getAoa(), b->getAmplitude(imd));
+        sdpmVectorXcd cfx(5);
+        cfx << cf[2].conjugate()(0), cf[1].conjugate()(0), cf[0](0), cf[1](0), cf[2](0);
+        sdpmVectorXcd cfz(5);
+        cfz << cf[2].conjugate()(2), cf[1].conjugate()(2), cf[0](2), cf[1](2), cf[2](2);
+        b->setUnsteadyDragCoeff(imd, ifq, (cf[1](0) * ca + cf[1](1) * sa + convolve(cfz, alpha)(4) * ca) / sRef);
+        b->setUnsteadySideforceCoeff(imd, ifq, (-cf[1](0) * sa + cf[1](1) * ca - convolve(cfz, alpha)(4) * sa) / sRef);
+        b->setUnsteadyLiftCoeff(imd, ifq, (-convolve(cfx, alpha)(4) + cf[1](2)) / sRef);
+
         // compute total
-        _cl += b->cl;
-        _cd += b->cd;
-        _cs += b->cs;
-        _cm += b->cm;
+        _cd1[imd][ifq] += b->getUnsteadyDragCoeff(imd, ifq);
+        _cs1[imd][ifq] += b->getUnsteadySideforceCoeff(imd, ifq);
+        _cl1[imd][ifq] += b->getUnsteadyLiftCoeff(imd, ifq);
+        _cm1[imd][ifq] += b->getUnsteadyMomentCoeff(imd, ifq);
     }
 }
 
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
index 16e0ebc1f83a95161aac741bf71c6a7819a83bbc..a4120442c1c712f37c9fbc6de4482083aa4980ce 100644
--- a/sdpm/src/sdpmSolver.h
+++ b/sdpm/src/sdpmSolver.h
@@ -36,19 +36,27 @@ class SDPM_API Solver : public sdpmObject
     // AICs
     size_t _nP;                     ///< number of body panels
     std::map<Element *, int> _rows; ///< map linking an element to a matrix cell
-    sdpmMatrixXd _A;                ///< body-body doublet matrix
-    sdpmMatrixXd _B;                ///< body-body source matrix
+    sdpmMatrixXd _A0;               ///< zero-frequency body-body doublet matrix
+    sdpmMatrixXd _B0;               ///< zero-frequency body-body source matrix
+    std::vector<sdpmMatrixXcd> _A;  ///< unsteady body-body doublet matrix
+    std::vector<sdpmMatrixXcd> _B;  ///< unsteady body-body source matrix
     // Solution
-    std::vector<sdpmDouble> _phi; ///< perturbation potential
-    std::vector<sdpmDouble> _tau; ///< source strength
-    std::vector<sdpmDouble> _mu;  ///< doublet strength
-    std::vector<sdpmVector3d> _u; ///< velocity
-    std::vector<sdpmDouble> _cp;  ///< pressure coefficient
+    std::vector<sdpmDouble> _phi;                             ///< perturbation potential
+    std::vector<sdpmDouble> _tau;                             ///< source strength
+    std::vector<sdpmDouble> _mu;                              ///< doublet strength
+    std::vector<sdpmVector3d> _u;                             ///< steady velocity
+    std::vector<sdpmDouble> _cp;                              ///< steady pressure coefficient
+    std::vector<std::vector<std::vector<sdpmVector3cd>>> _u1; ///< unsteady first harmonic velocity
+    std::vector<std::vector<std::vector<sdpmComplex>>> _cp1;  ///< unsteady first harmonic pressure coefficient
     // Aerodynamic loads
-    sdpmDouble _cl; ///< lift coefficient
-    sdpmDouble _cd; ///< drag coefficient
-    sdpmDouble _cs; ///< sideforce coefficient
-    sdpmDouble _cm; ///< pitch moment coefficient (positive nose-up)
+    sdpmDouble _cl;                             ///< steady lift coefficient
+    sdpmDouble _cd;                             ///< steady drag coefficient
+    sdpmDouble _cs;                             ///< steady sideforce coefficient
+    sdpmDouble _cm;                             ///< steady pitch moment coefficient (positive nose-up)
+    std::vector<std::vector<sdpmComplex>> _cl1; ///< unsteady first harmonic lift coefficient
+    std::vector<std::vector<sdpmComplex>> _cd1; ///< unsteady first harmonic drag coefficient
+    std::vector<std::vector<sdpmComplex>> _cs1; ///< unsteady first harmonic side force coefficient
+    std::vector<std::vector<sdpmComplex>> _cm1; ///< unsteady first harmonic pitch moment coefficient (positive nose-up)
 
 public:
     Solver(Problem &pbl);
@@ -58,9 +66,14 @@ public:
     sdpmDouble getDragCoeff() const { return _cd; }
     sdpmDouble getSideforceCoeff() const { return _cs; }
     sdpmDouble getMomentCoeff() const { return _cm; }
+    sdpmComplex getUnsteadyLiftCoeff(size_t m, size_t k) const { return _cl1[m][k]; }
+    sdpmComplex getUnsteadyDragCoeff(size_t m, size_t k) const { return _cd1[m][k]; }
+    sdpmComplex getUnsteadySideforceCoeff(size_t m, size_t k) const { return _cs1[m][k]; }
+    sdpmComplex getUnsteadyMomentCoeff(size_t m, size_t k) const { return _cm1[m][k]; }
 
     void update();
     void run();
+    void runUnsteady();
     void save(GmshExport const &writer, std::string const &suffix = "");
 
 #ifndef SWIG
@@ -68,8 +81,8 @@ public:
 #endif
 
 private:
-    void computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
-    void computeLoad();
+    void post(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
+    void postUnsteady(size_t m, size_t k, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu);
 };
 
 } // namespace sdpm
diff --git a/sdpm/src/sdpmWake.cpp b/sdpm/src/sdpmWake.cpp
index 27a78aa39ec0bbf19e42c5920733a82a93ddbc95..cfc935823385e743ac9e0497b84a390e0ebbf0b1 100644
--- a/sdpm/src/sdpmWake.cpp
+++ b/sdpm/src/sdpmWake.cpp
@@ -22,122 +22,97 @@
 #include <unordered_set>
 using namespace sdpm;
 
-#include <iostream>
-
-/**
- * @brief Create the wake in the case where the elements on the pressure side still have their original (suction side) trailing edge nodes
- */
-Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag) : Group(msh, name)
+Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe, double ux, std::map<Node *, Node *> const &iTeMap) : Group(msh, name)
 {
     // Create TE edges
-    createEdges();
+    std::vector<Element *> elems = _tag->getElements();
+    std::unordered_set<Edge *, EquEdge, EquEdge> teEdges;
+    for (size_t i = 0; i < nTe; ++i)
+    {
+        std::vector<Node *> nods = elems[i]->getNodes();
+        Edge *ed = new Edge({nods[0], nods[3]});
+        ed->el0 = elems[i];
+        teEdges.insert(ed);
+    }
 
-    // Find wing elements having an edge on the TE
-    for (auto e : wingTag->getElements())
+    // Find wing elements having an edge on the TE...
+    // ... the elements on the pressure side still have their original (suction side) trailing edge nodes
+    if (iTeMap.empty())
     {
-        std::vector<Node *> nods = e->getNodes();
-        size_t idx = 0;
-        for (size_t i = 0; i < nods.size(); ++i)
+        for (auto e : wingTag->getElements())
         {
-            // build wing edge and check if wing edge is a TE edge
-            std::vector<Node *> eN(2);
-            for (size_t j = 0; j < 2; ++j)
-                eN[j] = nods[(idx + j) % nods.size()];
-            Edge ed(eN);
-            auto it = _teEdges.find(&ed);
-            if (it != _teEdges.end())
+            std::vector<Node *> nods = e->getNodes();
+            size_t idx = 0;
+            for (size_t i = 0; i < nods.size(); ++i)
             {
-                // check normal direction
-                if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
-                    (*it)->el1 = e; // suction side
-                else
-                    (*it)->el2 = e; // pressure side
-                break;
+                // build wing edge and check if wing edge is a TE edge
+                std::vector<Node *> eN(2);
+                for (size_t j = 0; j < 2; ++j)
+                    eN[j] = nods[(idx + j) % nods.size()];
+                Edge ed(eN);
+                auto it = teEdges.find(&ed);
+                if (it != teEdges.end())
+                {
+                    // check normal direction
+                    if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
+                        (*it)->el1 = e; // suction side
+                    else
+                        (*it)->el2 = e; // pressure side
+                    break;
+                }
+                idx++;
             }
-            idx++;
         }
     }
-
-    // Map wake to wing elements
-    createMap(wingTag);
-}
-
-/**
- * @brief Create the wake in the case where the elements on the pressure side already have their correct (pressure side) trailing edge nodes
- */
-Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap) : Group(msh, name)
-{
-    // Create TE edges
-    createEdges();
-
-    // Find wing elements having an edge on the TE
-    for (auto e : wingTag->getElements())
+    // ... the elements on the pressure side already have their correct (pressure side) trailing edge nodes
+    else
     {
-        std::vector<Node *> nods = e->getNodes();
-        size_t idx = 0;
-        for (size_t i = 0; i < nods.size(); ++i, ++idx)
+        for (auto e : wingTag->getElements())
         {
-            // build wing edge and check if wing edge is a TE edge
-            std::vector<Node *> eN(2);
-            for (size_t j = 0; j < 2; ++j)
-                eN[j] = nods[(idx + j) % nods.size()];
-            Edge ed(eN);
-            auto it = _teEdges.find(&ed);
-            // try with actual nodes
-            if (it != _teEdges.end())
-            {
-                (*it)->el1 = e; // suction side
-                break;
-            }
-            // otherwise try again, but by replacing the nodes by their corresponding upper nodes
-            else
+            std::vector<Node *> nods = e->getNodes();
+            size_t idx = 0;
+            for (size_t i = 0; i < nods.size(); ++i, ++idx)
             {
                 // build wing edge and check if wing edge is a TE edge
                 std::vector<Node *> eN(2);
-                try
-                {
-                    for (size_t j = 0; j < 2; ++j)
-                        eN[j] = iTeMap.at(nods[(idx + j) % nods.size()]);
-                }
-                catch (const std::out_of_range &)
-                {
-                    continue;
-                }
+                for (size_t j = 0; j < 2; ++j)
+                    eN[j] = nods[(idx + j) % nods.size()];
                 Edge ed(eN);
-                auto it = _teEdges.find(&ed);
-                if (it != _teEdges.end())
+                auto it = teEdges.find(&ed);
+                // try with actual nodes
+                if (it != teEdges.end())
                 {
-                    (*it)->el2 = e; // pressure side
+                    (*it)->el1 = e; // suction side
                     break;
                 }
+                // otherwise try again, but by replacing the nodes by their corresponding upper nodes
+                else
+                {
+                    // build wing edge and check if wing edge is a TE edge
+                    std::vector<Node *> eN(2);
+                    try
+                    {
+                        for (size_t j = 0; j < 2; ++j)
+                            eN[j] = iTeMap.at(nods[(idx + j) % nods.size()]);
+                    }
+                    catch (const std::out_of_range &)
+                    {
+                        continue;
+                    }
+                    Edge ed(eN);
+                    auto it = teEdges.find(&ed);
+                    if (it != teEdges.end())
+                    {
+                        (*it)->el2 = e; // pressure side
+                        break;
+                    }
+                }
             }
         }
     }
 
-    // Map wake to wing elements
-    createMap(wingTag);
-}
-
-/**
- * @brief Create set of TE edges
- */
-void Wake::createEdges()
-{
-    for (auto e : _tag->getElements())
-    {
-        std::vector<Node *> nods = e->getNodes();
-        Edge *ed = new Edge({nods[0], nods[3]});
-        ed->el0 = e;
-        _teEdges.insert(ed);
-    }
-}
-
-/**
- * @brief Create a map between the wake and wing elements
- */
-void Wake::createMap(Tag const *wingTag)
-{
-    for (auto ed : _teEdges)
+    // Map wake elements on the TE to wing elements
+    for (auto ed : teEdges)
     {
         if (ed->el1 == NULL || ed->el2 == NULL)
         {
@@ -148,6 +123,19 @@ void Wake::createMap(Tag const *wingTag)
         _wwMap[ed->el0] = std::pair<Element *, Element *>(ed->el1, ed->el2);
         delete ed;
     }
+    // Map the rest of the wake and fill time lag map
+    sdpmDouble dt = (elems[0]->getNodes()[1]->getCoords()(0) - elems[0]->getNodes()[0]->getCoords()(0)) / ux;
+    size_t nw = elems.size() / nTe; // number of shed wake rows
+    for (size_t j = 0; j < nTe; ++j)
+    {
+        _lagMap[elems[j]] = dt;
+        std::pair<Element *, Element *> p = _wwMap.at(elems[j]);
+        for (size_t i = 1; i < nw; ++i)
+        {
+            _wwMap[elems[i * nTe + j]] = p;
+            _lagMap[elems[i * nTe + j]] = (i + 1) * dt;
+        }
+    }
 }
 
 void Wake::write(std::ostream &out) const
diff --git a/sdpm/src/sdpmWake.h b/sdpm/src/sdpmWake.h
index 338892bdf8c082fa7cf3e685981d6685b8642368..61966bc7ab6c54b74e08524a46c366c113ff8f3c 100644
--- a/sdpm/src/sdpmWake.h
+++ b/sdpm/src/sdpmWake.h
@@ -20,8 +20,6 @@
 #include "sdpm.h"
 #include "sdpmGroup.h"
 #include "sdpmTag.h"
-#include "sdpmEdge.h"
-#include <unordered_set>
 #include <map>
 
 namespace sdpm
@@ -33,24 +31,19 @@ namespace sdpm
  */
 class SDPM_API Wake : public Group
 {
-    std::unordered_set<Edge *, EquEdge, EquEdge> _teEdges;       ///< trailing edge edges
     std::map<Element *, std::pair<Element *, Element *>> _wwMap; ///< wing-wake map
+    std::map<Element *, double> _lagMap;                         ///< time lag map
 
 public:
-    Wake(Mesh &msh, std::string const &name, Tag *const &wingTag);
-    Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap);
+    Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe, double ux, std::map<Node *, Node *> const &iTeMap = {});
     ~Wake() {}
 
-private:
-    void createEdges();
-    void createMap(Tag const *tag);
-
-public:
 #ifndef SWIG
-    std::map<Element *, std::pair<Element *, Element *>> const &getMap() const
+    std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const
     {
         return _wwMap;
     }
+    std::map<Element *, double> const &getLagMap() const { return _lagMap; }
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/sdpm/src/sdpm_extras.h b/sdpm/src/sdpm_extras.h
index 981800ab51678eabbf80313255baf2fed749c833..ef7b770e225d485c4dba669769115325ea3ce23d 100644
--- a/sdpm/src/sdpm_extras.h
+++ b/sdpm/src/sdpm_extras.h
@@ -18,12 +18,8 @@
 #define SDPM_EXTRAS_H
 
 #include "sdpm.h"
-
-#include <iostream>
-#include <fstream>
+#include <ostream>
 #include <vector>
-#include <string>
-#include <Eigen/Sparse>
 
 namespace std
 {
@@ -48,19 +44,21 @@ ostream &operator<<(ostream &out, vector<T> const &v)
 namespace sdpm
 {
 /**
- * @brief save CSR matrix to a file
- *
- * matlab usage: load K.txt; K=spconvert(K)
+ * @brief convolve two vectors
  */
-SDPM_API void tomatlab(std::string const &fname, Eigen::SparseMatrix<double> const &A)
+SDPM_API sdpmVectorXcd convolve(sdpmVectorXcd const &u, sdpmVectorXcd const &v)
 {
-    std::ofstream f(fname.c_str());
-
-    for (int k = 0; k < A.outerSize(); ++k)
-        for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
-            f << it.row() + 1 << ' ' << it.col() + 1 << ' ' << it.value() << '\n';
-
-    f.close();
+    int m = static_cast<int>(u.size());
+    int n = static_cast<int>(v.size());
+    sdpmVectorXcd w = sdpmVectorXcd::Zero(n + m - 1);
+    for (int i = 0; i < w.size(); ++i)
+    {
+        int jmin = (i - n + 1 >= 0) ? i - n + 1 : 0;
+        int jmax = (i - m + 1 < 0) ? i : m - 1;
+        for (int j = jmin; j <= jmax; ++j)
+            w(i) += u(j) * v(i - j);
+    }
+    return w;
 }
 } // namespace sdpm
 
diff --git a/sdpm/tests/lann.py b/sdpm/tests/lann.py
new file mode 100644
index 0000000000000000000000000000000000000000..031b1833e259d070a98fdbca43ceb27d999674d2
--- /dev/null
+++ b/sdpm/tests/lann.py
@@ -0,0 +1,86 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2023 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# LANN wing
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import math
+import cmath
+
+def main():
+    # geometry parameters
+    wnc = 20 # number of chordwise panels
+    pars = {'nC': wnc, 'bC': 0.25, 'rS': 1}
+    # flow and time parameters
+    minf = 0.621 # freestream Mach number
+    aoa = 0.59 * math.pi / 180
+    omega = 2 * math.pi * 24 # circular frequency
+    amp = 0.5 * 0.25 * math.pi / 180 # semi pitching amplitude
+    uinf = 0.5 * 0.361 * omega / 0.133; # freestream velocity
+
+    # start timer
+    tms = sdpm.Timers()
+    tms['total'].start()
+    # mesh
+    tms['msh'].start()
+    msh = MeshLoader(build_fpath('models/lann.geo', __file__, 1)).run(pars)
+    tms['msh'].stop()
+    # problem
+    tms['pbl'].start()
+    pbl = sdpm.Problem(msh, 1, [omega])
+    pbl.setGeometry(0.253, 0.361, 0.224, 0., 0., True)
+    pbl.updateFreestream(aoa, 0., minf, uinf)
+    wing = sdpm.Body(msh, 'wing', 'te', 0.361, wnc, 1, 1, uinf*math.cos(aoa))
+    wing.addMotion(uinf*math.cos(aoa), amp, 0.224/math.sqrt(1-minf*minf))
+    pbl.addBody(wing)
+    tms['pbl'].stop()
+    ## solver
+    tms['sol'].start()
+    sol = sdpm.Solver(pbl)
+    sol.update()
+    sol.run()
+    sol.runUnsteady()
+    sol.save(sdpm.GmshExport(msh))
+    tms['sol'].stop()
+    # end timer
+    tms['total'].stop()
+    print(tms)
+
+    # test
+    tests = CTests()
+    tests.add(CTest('CL0', sol.getLiftCoeff(), 0.258, 5e-2)) # MATLAB: 0.265
+    tests.add(CTest('CD0', sol.getDragCoeff(),  -0.0029, 5e-2)) # MATLAB: -0.0033
+    tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.008, 5e-2))
+    tests.add(CTest('CM0', sol.getMomentCoeff(), -0.079, 5e-2))
+    tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.54, 5e-2)) # MATLAB: 1.62
+    tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.66, 5e-2)) # MATLAB: 2.6
+    tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.038, 5e-2)) # MATLAB: 0.041
+    tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 14.6, 5e-2)) # MATLAB: 16.8
+    tests.add(CTest('Abs(CS1)', abs(sol.getUnsteadySideforceCoeff(0, 0)) / amp / math.pi, 0.037, 5e-2))
+    tests.add(CTest('Ang(CS1)', cmath.phase(sol.getUnsteadySideforceCoeff(0, 0)) * 180 / math.pi,  -10.0, 5e-2))
+    tests.add(CTest('Abs(CM1)', abs(sol.getUnsteadyMomentCoeff(0, 0)) / amp / math.pi, 0.33, 5e-2))
+    tests.add(CTest('Ang(CM1)', cmath.phase(sol.getUnsteadyMomentCoeff(0, 0)) * 180 / math.pi, -170.5, 5e-2))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()
diff --git a/sdpm/tests/m6.py b/sdpm/tests/m6.py
index eff04896bb4159e048b434ff765e7edc573aa31a..528d4b88023489dbac242f33c4cda4a38b71ec2e 100644
--- a/sdpm/tests/m6.py
+++ b/sdpm/tests/m6.py
@@ -25,9 +25,10 @@ import math
 
 def main():
     # geometry parameters
-    pars = {'wNc': 100, 'wNs': 10, 'bC': 0.25, 'pS': 1.0}
+    pars = {'nC': 100, 'nS': 10, 'bC': 0.25, 'pS': 1.0}
     # flow parameters
     aoa = 3.06 * math.pi / 180
+    minf = 0.
 
     # start timer
     tms = sdpm.Timers()
@@ -38,8 +39,10 @@ def main():
     tms['msh'].stop()
     # problem
     tms['pbl'].start()
-    pbl = sdpm.Problem(msh, aoa, 0, 0, 0.7528, 0.64, 0, 0, 0, True)
-    wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(0.7528, 0.64, 0., 0., 0., True)
+    pbl.updateFreestream(aoa, 0., minf, 1.)
+    wing = sdpm.Body(msh, 'wing', 'te', 0.805)
     pbl.addBody(wing)
     tms['pbl'].stop()
     ## solver
diff --git a/sdpm/tests/n12.py b/sdpm/tests/n12.py
index d6fc9c446b518cc9061fb60effb8a594ce39bb66..25bb02de214f06abc93d571f946947f5e56f86a5 100644
--- a/sdpm/tests/n12.py
+++ b/sdpm/tests/n12.py
@@ -26,7 +26,7 @@ import math
 def main():
     # geometry parameters
     sym = 0
-    span = 5 # true span if sym=0, half span if sym=1
+    span = 5 if sym==0 else 2.5 # span
     pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10}
     # flow parameters
     aoa = 3 * math.pi / 180
@@ -40,8 +40,10 @@ def main():
     tms['msh'].stop()
     # problem
     tms['pbl'].start()
-    pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
-    wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(span, 1., 0., 0., 0., bool(sym))
+    pbl.updateFreestream(aoa, 0., 0., 1.)
+    wing = sdpm.Body(msh, 'wing', 'wTe', 1.)
     pbl.addBody(wing)
     tms['pbl'].stop()
     ## solver
diff --git a/sdpm/tests/n12tail.py b/sdpm/tests/n12tail.py
index aee915250f7dd7287f075e0251344d1c401c60fb..d76c7ed4445494e9c1c43e03d380690e5596fd71 100644
--- a/sdpm/tests/n12tail.py
+++ b/sdpm/tests/n12tail.py
@@ -26,9 +26,9 @@ import math
 def main():
     # geometry parameters
     sym = 0
-    span = 5 # wing (true span if sym=0, half span if sym=1)
-    tspan = 3 # tail (true span if sym=0, half span if sym=1)
-    pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10, 'tSpn': tspan, 'tNc': 100, 'tNs': 6}
+    wspan = 5 if sym==0 else 2.5 # wing span
+    tspan = 3 if sym==0 else 1.5 # tail span
+    pars = {'symY': sym, 'wSpn': wspan, 'wNc': 100, 'wNs': 10, 'tSpn': tspan, 'tNc': 100, 'tNs': 6}
     # flow parameters
     aoa = 3 * math.pi / 180
 
@@ -41,9 +41,11 @@ def main():
     tms['msh'].stop()
     # problem
     tms['pbl'].start()
-    pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
-    wing = sdpm.Body(msh, 'wing', 'wTe', 10)
-    tail = sdpm.Body(msh, 'tail', 'tTe', 15)
+    pbl = sdpm.Problem(msh)
+    pbl.setGeometry(wspan, 1., 0., 0., 0., bool(sym))
+    pbl.updateFreestream(aoa, 0., 0., 1.)
+    wing = sdpm.Body(msh, 'wing', 'wTe', 1.0)
+    tail = sdpm.Body(msh, 'tail', 'tTe', 0.5)
     pbl.addBody(wing)
     pbl.addBody(tail)
     tms['pbl'].stop()
@@ -60,14 +62,14 @@ def main():
 
     # test
     tests = CTests()
-    tests.add(CTest('Wing CL', wing.cl, 0.230, 5e-2))
-    tests.add(CTest('Wing CD', wing.cd, 0.003, 5e-2))
-    tests.add(CTest('Wing CS', wing.cs, 0.000, 5e-2))
-    tests.add(CTest('Wing CM', wing.cm, -0.056, 5e-2))
-    tests.add(CTest('Tail CL', tail.cl, 0.045, 5e-2))
-    tests.add(CTest('Tail CD', tail.cd, 0.001, 5e-2))
-    tests.add(CTest('Tail CS', tail.cs, 0.000, 5e-2))
-    tests.add(CTest('Tail CM', tail.cm, -0.229, 5e-2))
+    tests.add(CTest('Wing CL', wing.getLiftCoeff(), 0.230, 5e-2))
+    tests.add(CTest('Wing CD', wing.getDragCoeff(), 0.003, 5e-2))
+    tests.add(CTest('Wing CS', wing.getSideforceCoeff(), 0.000, 5e-2))
+    tests.add(CTest('Wing CM', wing.getMomentCoeff(), -0.056, 5e-2))
+    tests.add(CTest('Tail CL', tail.getLiftCoeff(), 0.045, 5e-2))
+    tests.add(CTest('Tail CD', tail.getDragCoeff(), 0.001, 5e-2))
+    tests.add(CTest('Tail CS', tail.getSideforceCoeff(), 0.000, 5e-2))
+    tests.add(CTest('Tail CM', tail.getMomentCoeff(), -0.229, 5e-2))
     tests.run()
 
     # eof