```/*
DEMO4.I
Airfoil demo

\$Id\$
*/
/*    Copyright (c) 1995.  The Regents of the University of California.

func demo4(mono=)
/* DOCUMENT demo4
or demo4, mono=1
solves for the flow past a 2D airfoil using Kutta-Jakowski theory.
The colors represent static pressure (darker is lower pressure),
red lines are streamlines.
Solutions for various angles of attack are shown by animation.
With the mono=1 keyword, only the streamlines are shown.  (On a
monochrome terminal, the pressure makes the streamlines invisible.)
*/
{
require, "movie.i";
local nx, ny;
imax= 50;
movie, display;
display, imax;
}

func display(i)
{
attack= 28.*double(imax-i)/(imax-1);
solve, attack, 2.0, 0.2;
if (i==1) {
extern cn, cx;
cn= 0.8*min(pressure);
cx= max(pressure);
}
limits, -2.8, 2.4, -2.4, 2.8;
plmesh, z.im,z.re;
if (!mono) plf, zncen(pressure), cmin=cn,cmax=cx;
plc, potential.im, marks=0,color="red", levs=span(-2.5,3.5,33);
plg,z.im(,1),z.re(,1), marks=0;
return i<imax;
}

func solve (attack, chord, thick)
{
extern w, z, potential, velocity, pressure; /* outputs */
attack*= pi/180.;
a= 0.5*chord;
w= get_mesh(a, attack);
/* the log(w) term adds just enough circulation to move the rear
stagnation point to the trailing edge -- the Kutta condition */
potential= jakowski(w, a) + 2i*a*sin(attack)*log(w);
dpdw= 1.-(a/w)^2 + 2i*a*sin(attack)/w;
emith= exp(-1i*attack);
a= a*emith;
b= thick*emith;  /* could add camber here too someday? */
w-= b;
a-= b;
z= jakowski(w, a);
dzdw= 1.-(a/w)^2;
velocity= conj(dpdw/dzdw);
pressure= 0.5*(1.0-abs(velocity)^2);
}

func get_mesh (a, attack)
{
/* get a mesh in the w-plane
-- the plane in which the airfoil is a circle
the mesh splits at the point which will become the trailing edge */
if (is_void(nx)) nx= 120;
if (is_void(ny)) ny= 30;
a= abs(a);
theta= span(1.e-9, 2*pi-1.e-9, nx) - attack;
r= a/span(1.0+1.e-9, 0.25, ny)(-,);
return r*exp(1i*theta);
}

func jakowski (z, a)
{
/* Jakowski transform - circle of radius a into slot of length 4a */
return z + a*a/z;
}

func ijakowski (z, a)
{
/* Inverse Jakowski - slot back into circle (unused here) */
z*= 0.5;
a= complex(a);
sgn= sign((z*conj(a)).re);
return z + sgn(z);
}
```