File	Date	Author	Commit
inst	2010-01-05	highegg	[3030fe] improve clqinterp
src	2017-01-02	John W. Eaton	[34895a] F77_INT fixes for Octave 4.4.
.hgignore	2016-02-29	Carnë Draug	[b3bdfa] maint: add hgignore file.
COPYING	2008-06-18	highegg	[d5c96c] NLWing2: initial commit
ChangeLog	2010-01-05	highegg	[83e255] nlwing2: implement clqinterp and qcalc_extra
DESCRIPTION	2009-12-03	highegg	[6ac85e] NLWing2 1.2.0
INDEX	2009-12-03	highegg	[6ac85e] NLWing2 1.2.0
README	2008-07-18	highegg	[19fda7] update GPL2 files to GPL3
TODO	2008-06-18	highegg	[d5c96c] NLWing2: initial commit

Read Me

Global coordinate system
--------------------------------------------------------------------------------

The spatial coordinates are bound to the aircraft:
x : frontal 
y : vertical 
z : spanwise left (horizontal)
                           
                         ^ ____                                                            
                       y |/   /            
                         |   /             
                        /|  /     __       
                  _____/ | /_____/ _/     x
                 /       |      /__/ ------------->
                 \_____    _______/        
                       \   \               
                        \ \ \              
                         \ \ \             
                          \_\_\            
                             \  z           
                              \| 

x-y is the wing symmetry plane, while the wing lies along z axis. The origin is 
the center of the wing.                           

Maths & Physics
--------------------------------------------------------------------------------
The method is based on approximation the flowfield around the wing by the 
potential flow induced by an ensemble of horseshoe vortices, with bound 
segments aligned along the wing's centerline, and the free segments shedding
in the streamwise direction. Their strengths, i.e. circulations, are unknown
and need to be determined using the flow equations. With known strengths,
the flowfield velocity can be calculated at an arbitrary point using 
Biot-Savart law.

The control (unknown) spanwise quantity on i-th vortex is
  +---------------------+
  |  g_i = gamma_i/c_i  | (1)
  +---------------------+
where gamma_i is the circulation and ch_i is local chord length (this gives
approximately local cl). 

The flow equation to be satisfied locally follows from expressing
the lift created on a particular section in two ways:
from the Kutta-Joukowski law (the left-hand side) and
using the 2D section data (the right-hand side) 
  +-------------------------------------------------------+
  |       rho * v * gamma = 1/2 * rho * v^2 * cl * c      | (2)
  +-------------------------------------------------------+

giving
  +----------------------------+
  |       g = v * cl / 2       | (3)
  +----------------------------+

where
  +--------------------------+
  |       cl = cl(alfa)      |
  |    alfa = atan(vy/vx)    | (4)
  |   v = sqrt(vx^2 + vy^2)  |
  +--------------------------+

Numerics
--------------------------------------------------------------------------------
Using Biot-Savart law, it is possible to express the induced velocity at each 
collocation point caused by unit circulation on each vortex, thus obtain the
"influence" tensor. Thus, from g_i, vx_i and vy_i are obtained by applying 
appropriate influence matrices. atan(vy_i/vx_i) then gives alfa_i - see (4) 
(at this point, nonlinearity is introduced). cl_i is then obtained from alfa_i
by interpolating the provided 2D section data (a combination of spanwise and
angle-wise interpolation is used). Finally, (3) closes the cycle, arriving
at g_i again.
In this fashion, for any global angle of attack we obtain a system of nonlinear
equations. This parametric system is solved by starting at a low angle of attack
and tracking the nonlinear solution to higher angles while possible, using a
predictor-corrector strategy.
The 2D section lifts are interpolated in two dimensions: by spanwise coordinate
(combining different datasets) and by angle of attack. Instead of a naive
interpolation, a "feature" interpolation is done: at first, each lifting line is
analyzed for a zero-lift and maximum-lift angle. The lifting line is then scaled
to 0,1, and these three features (zero-lift angle, max-lift angle, and scaled
lifting line) are interpolated spanwise by linear interpolation.