# Processing to draw NACA 4 digit airfoils

I want to use processing to draw NACA 4 digit airfoils. Ive been to lots of NASA websites, read a few papers, found some python code, messed with a java applet and still got NADA. I understand basically how they are generated. For instance, a NACA 2412 airfoil means something specific. The first number is how much camber the foil has. This is a variable that goes from 0 to 9. the 4 means that the highest camber point occurs at the 40% chord(length) point. This too is a variable from 0 to 9. The 12 means that the foil is a maximum of 12% thickness, as a percentage of chord length and drawn perpendicular to the camber line. This thickness is also a variable from 1 to 40.

There are three steps to generating the outer envelope which is the actual airfoil shape. The first step is to calculate x and y points along the chord length for the camber line. I got the formula from this site:

Airfoil Tools

From that formula from the link, I have to calculate y for the camber line for every x along the chord length and calculate the first derivative(slope) at each one of those. So, hereâ€™s what I think I need to put into code:

Pick the number of points to calculate and make it an int.

Divide up the chord length from 0 to 1 by the number of points to get the x value.

So, would that look like this? Iâ€™ve only calculated the upper surface but I get this giant leg at the 40% chord and its pointed down. What did I do wrong?

``````int N = 20;  //number of points
float M = 2/10.0;
float[] slope = new float[20];
float T = 12/100.0;

float x[] = new float[20];  //divide up unit chord length by N
float y[] = new float[20];  //divide up unit chord length by N
float xU[] = new float[20];
float yU[] = new float[20];
float P = 4/10.0;
PShape s;
float a0 = 0.2969;
float a1 = -0.126;
float a2 = -0.3516;
float a3 = 0.2843;
float a4 = -0.1036;
float[] theta = new float[20];

float[] yt = new float[20];

void setup()
{
size( 650, 300, P2D); //P3D );

}

void draw() {
background(0, 0, 0);

for (int i = 0; i < N; i++){

x[i] = (1/20.0) * i;
// println(x[i]);
if (x[i] < P){

y[i] = M/(P*P) * (2*P*x[i] - x[i]*x[i]);

slope[i] = (2*M)/(P*P) * (P - x[i]);

}

if (x[i] >= P){

y[i] = M/(1-P)*(1-P) * (1 - 2*P + 2*P*x[i] - x[i]*x[i]);
slope[i] = 2*M/((1-P)*(1-P)) * (P - x[i]);

}

yt[i] = T/0.2* (sqrt(a0*x[i]) + a2*(x[i]*x[i]) + a3*(x[i]*x[i]*x[i]) + a4*(x[i]*x[i]*x[i]*x[i]));
theta[i] = atan(slope[i]);
// println(theta);
xU[i] = x[i] - yt[i] * sin(theta[i]);
yU[i] = y[i] + yt[i] * cos(theta[i]);

}  //end for loop
s = createShape();

s.beginShape();
// s.fill(125);
s.noFill();
s.stroke(0,255,0);
s.strokeWeight(3);
for (int j = 0; j<20; j++) s.vertex(200*xU[j], -200*yU[j]);
s.endShape(CLOSE);

translate(width/2,height/2);
shape(s, 0, 0);

noLoop();

}
``````

To save someone the trouble of trying the code to see the error. Hereâ€™s a shot of the problem.

If I read you right, I believe you did not take into account Processingâ€™s Y axis is inverted as zero value is located in the top edge of your sketch and it increases as you approach the bottom edge of your sketch. What you need is a cartesian coordinate system where positive increases upward. The following lines:

``````  translate(0,height);
scale(1,-1);
``````

at the beginning of `draw()` should do the trick.

Kf

2 Likes

Hmmm, that flipped it over. I added the lower surface and its really messed up now.

Hereâ€™s a comparison of what it should look like vs. how it looks if I try to generate it.

I think I might have found the problem. I println(x[I]) and it goes from 0 to 0.4875 and it should go from 0 to 1. So its not calculating anything for the whole length of the unit foil.

Edit: That didnâ€™t lead to much improvement. It does divide up the xâ€™s from 0 to 1 now, but its still screwed up.
Digging in a little further, if I plot just the camber line, the jog is still there and everything is built on top of the camber so any problem there is going to propogateâ€¦ Its like it hits the 40% point and the scale changes.

Edit2: Hereâ€™s what I tried. I went to the â€śAirfoil Toolsâ€ť link and downloaded the .dat file for the airfoil Iâ€™m trying to generate. I figured maybe the way Iâ€™m dividing up the xâ€™s is part of the problem. I now read in the xâ€™s from the file with the same number of points as my sketch. I use these xâ€™s to calculate the camber, thickness and final envelope. Thereâ€™s an improvement in that the shape is closed, but the jog is still there. I plotted the upper surface, the camber line, the lower surface and it still looks like a cow or duck, not an airfoil.

Edit3: I tried one more thing. I set the camber variable and camber high point variable to zero. That means the slope of the camber is also zero. It looks more like an airfoil except it is truncated at about 20%. I donâ€™t understand why it doesnâ€™t draw between 0 and 1 when the xâ€™s go from 0 to 1. I donâ€™t know. Heres what it looks like.

I got something that works but as soon as I move the horizontal slider, it is back to drawing the upper and lower at an angle other than level.

``````//https://discourse.processing.org/t/processing-to-draw-naca-4-digit-airfoils/5739

int N = 81;  //number of points
float M = -0.10000f; //2/100.0;
float[] dydx = new float[81];
float T = 12/100.0;
float dT;

float x[] = new float[81];  //divide up unit chord length by N
float y[] = new float[81];  //divide up unit chord length by N
float xU[] = new float[81];
float yU[] = new float[81];
float xL[] = new float[81];
float yL[] = new float[81];
float P = 0.4; //40/100.0;
PShape s,s2,s3,s4,s5;
float a0 = 0.2969;
float a1 = -0.126;
float a2 = -0.3516;
float a3 = 0.2843;
float a4 = -0.1036;
float[] theta = new float[81];

float[] yt = new float[81];
String[] coordinates = new String[81];

float[] beta = new float[81];

Table table;
TableRow row;

String nameOfFile = "Airfoil";
float Dat[]= new float[81];

Controls controls;
HorizontalControl controlX;
int showControls;
boolean draggingZoomSlider = false;
boolean released = true;
float zoom = -200f;
float tzoom = 120.37037f;
float velocityX = 0;
float velocityY = 0;

void setup()
{
//  fullScreen(P3D);
size( 650, 300, P3D);

controls = new Controls();
controlX = new HorizontalControl();
showControls = 1;
table = new Table();
makeCSVFileRows("",Dat);
WriteDataToCSVFile(Dat, "");
WriteDataToCSVFile(Dat, "");
saveTable(table, "data/"+nameOfFile+".csv");
getOrdinates();
airfoiltest(xU,yU);
airfoiltest(xL,yL);
println ("done.");
}

void draw() {

if (mousePressed) {
if( (showControls == 1) && (controls.isZoomSliderEvent(mouseX, mouseY)) || ( showControls == 1 && controlX.isZoomSliderEvent(mouseX,mouseY))) {
draggingZoomSlider = true;

zoom = controls.getZoomValue(mouseY);

tzoom = controlX.getZoomValue(mouseX,mouseY);

M = map(M,-0.10,0.10,-264.00,264.00); //,-0.10,0.10);
T = map(T,-12,12,-tzoom,tzoom);
}
else if (!draggingZoomSlider) {
if (released == true){
velocityX += (mouseY-pmouseY) * 0.01;
velocityY -= (mouseX-pmouseX) * 0.01;

}
}
}

background(0, 0, 0);

if (showControls == 1) {
controls.render();
controlX.render();
}

controls.updateZoomSlider(zoom);
controlX.updateZoomSlider(tzoom);

translate(width/2-450/2, height/2);

scale(1.5,1);
//makeCSVFile("Airfoil");

makeCSVFileRows(nameOfFile,Dat);

dT = sqrt((x[80]-x[0])*(x[80]-x[0]) + (y[80]-y[0])*(y[80]-y[0]));
dT = asin(y[80]/dT);
//  println(dT*1000);
println(zoom);

shape(airfoiltest(xU,yU), 0, 0);

shape(airfoiltest(xL,yL), 0, 0);
//noLoop();
}

PShape airfoiltest(float[] xvalue, float[] yvalue) {

s = createShape();
s.beginShape();
// s.fill(125);
s.noFill();
s.stroke(0,255,0);
s.strokeWeight(3);

for (int i = 0; i<coordinates.length-1; i++){
s.vertex(300*xvalue[i],-300* (yvalue[i]));
}

s.endShape(CLOSE);
return s;
}

//println("Loaded " + coordinates.length + " coordinates:");
//   println( ": ", String.format("%.5f %.5f", "", ""));
for (int i = 0; i < coordinates.length; i++) {
String[] ordinate = splitTokens(coordinates[i], ",");
if(i >81 && i <82){
x[i] = float(ordinate[i]);
y[i] = float(ordinate[i]);
println(i, ": ", String.format("%.5f %.5f", x[i], y[i]));
}
//  println(i, ": ", String.format("%.5f %.5f", x[i], y[i]));
}
}

void WriteDataToCSVFile(float array[], String nameOfData) {

for (int i = 0; i< array.length; i++) {
table.setFloat(i, nameOfData, array[i]);
}//for
}//func

void makeCSVFileRows(String name, float array[]) {
for (int i = 0; i< array.length; i++) {
}
}

void getOrdinates(){

for (int i = 0; i < N; i++){

M = zoom/1000f;
T = tzoom/1000f;

//      String[] ordinate = splitTokens(coordinates[i], " ");
if(i >= 0){
x[i] = (1 - cos(beta[i]))/2;
}
if(i < 40){
if (x[i] >= 0){
y[i] = (M/(P*P) * (2*P*x[i] - x[i]*x[i]));
dydx[i] = (2*M)/(P*P) * (P - x[i]);
}
}
if(i >=40){
if (x[i] <= 1.00000f){
y[i] = (M/(1-P)*(1-P) * (1 - 2*P + 2*P*x[i] - x[i]*x[i]))* 2.75; //*2.686;
dydx[i] = (2*M/((1-P)*(1-P)) * (P - x[i]));
}
}
}  //end for loop

for (int i = 0; i < 81; i++){
if (i >= 0 && i < 81) {
yt[i] = (T/0.2* (sqrt(a0*x[i])+ a1*x[i] + a2*(x[i]*x[i]) + a3*(x[i]*x[i]*x[i]) + a4*(x[i]*x[i]*x[i]*x[i])));
theta[i] = (atan((dydx[i])));
xU[i] = x[i] - yt[i] * -(sin(radians(theta[i])));
//  WriteDataToCSVFile(xU,"");
yU[i] = ((y[i] + yt[i]  * (cos(radians(theta[i]))))); //*atan((y[i]-yt[i])/x[i]);
//  WriteDataToCSVFile(yU,"");
xL[i] = x[i] + yt[i] * sin(radians(theta[i]));
//  WriteDataToCSVFile(xL,"");
yL[i] = (y[i] - yt[i] * cos(radians(theta[i]))); //*asin((yt[i]*slope[i]));
//  WriteDataToCSVFile(yL,"");
}
}

}

void mouseReleased() {

if (released == true  && draggingZoomSlider == true){

}

draggingZoomSlider = false;
}
``````

Put this in a new tab:

``````/*

Kepler Visualization - Controls

GUI controls added by Lon Riesberg, Laboratory for Atmospheric and Space Physics
lon@ieee.org

April, 2012

Current release consists of a vertical slider for zoom control.  The slider can be toggled
on/off by pressing the 'c' key.

Slide out controls that map to the other key bindings is currently being implemented and
will be released soon.

*/

class Controls {

int barWidth;
int barX;                          // x-coordinate of zoom control
int minY, maxY;                    // y-coordinate range of zoom control
float minZoomValue, maxZoomValue;  // values that map onto zoom control
float valuePerY;                   // zoom value of each y-pixel
int sliderY;                       // y-coordinate of current slider position
float sliderValue;                 // value that corresponds to y-coordinate of slider
int sliderWidth, sliderHeight;
int sliderX;                       // x-coordinate of left-side slider edge

Controls () {

barX = 40;
barWidth = 15;

minY = 40;
maxY = minY + height/3 - sliderHeight/2;

minZoomValue = height - height;
maxZoomValue = height;   // 300 percent
valuePerY = (maxZoomValue - minZoomValue) / (maxY - minY);

sliderWidth = 25;
sliderHeight = 10;
sliderX = (barX + (barWidth/2)) - (sliderWidth/2);
sliderValue = minZoomValue;
sliderY = minY;
}

void render() {

// strokeWeight(1.5);
strokeWeight(1);
//  stroke(105, 105, 105);  // fill(0xff33ff99);
//   stroke(0xff33ff99);  // fill(0xff33ff99);  0xffff0000
stroke(0xffff0000);

// zoom control bar
fill(0, 0, 0, 0);

rect(barX, minY, barWidth, maxY-minY);

// slider
// fill(105, 105, 105); //0x3300FF00
fill(0xffff0000); // 0xff33ff99//0x3300FF00
rect(sliderX, sliderY, sliderWidth, sliderHeight);
}

float getZoomValue(int y) {
if ((y >= minY) && (y <= (maxY - sliderHeight/2))) {
sliderY = (int) (y - (sliderHeight/2));
if (sliderY < minY) {
sliderY = minY;
}
sliderValue = (y - minY) * valuePerY + minZoomValue;
}
return sliderValue;
}

void updateZoomSlider(float value) {
int tempY = (int) (value / valuePerY) + minY;
if ((tempY >= minY) && (tempY <= (maxY-sliderHeight))) {
sliderValue = value;
sliderY = tempY;
}
}

boolean isZoomSliderEvent(int x, int y) {
int slop = 50;  // number of pixels above or below slider that's acceptable.  provided for ease of use.
int sliderTop = (int) (sliderY - (sliderHeight/2)) - slop;
int sliderBottom = sliderY + sliderHeight + slop;
return ((x >= sliderX) && (x <= (sliderX    + sliderWidth)) && (y >= sliderTop)  && (y <= sliderBottom) || draggingZoomSlider );
}
}
``````

And this in another new tab:

``````
/*
I modified this so the slider is horizontal.  That gives me a vertical for
tweaking altitude and horizontal for right ascension/longitude
*/

/*

Kepler Visualization - Controls

GUI controls added by Lon Riesberg, Laboratory for Atmospheric and Space Physics
lon@ieee.org

April, 2012

Current release consists of a vertical slider for zoom control.  The slider can be toggled
on/off by pressing the 'c' key.

Slide out controls that map to the other key bindings is currently being implemented and
will be released soon.

*/

class HorizontalControl {

int barHeight;
int barY;                          // y-coordinate of zoom control
int minX, maxX;                    // x-coordinate range of zoom control
float minZoomValue, maxZoomValue;  // values that map onto zoom control
float valuePerX;                   // zoom value of each y-pixel
int sliderY;                       // y-coordinate of current slider position
float sliderValue;                 // value that corresponds to y-coordinate of slider
int sliderWidth, sliderHeight;
int sliderX;                       // x-coordinate of left-side slider edge

HorizontalControl () {

barY = 15; //40;
barHeight = 40; //15;

minX = 40;
maxX = minX + width/3 - sliderWidth/2;

minZoomValue = width - width;
maxZoomValue = width;   // 300 percent
valuePerX = (maxZoomValue - minZoomValue) / (maxX - minX);

sliderWidth = 10; //25;
sliderHeight = 25; //10;
// sliderY = (barY + (barHeight/2)) - (sliderHeight/2);
sliderY = (barY - (sliderHeight/2)) + (barHeight/2);
sliderValue = minZoomValue;
sliderX = minX;
}

void render() {
pushMatrix();

// strokeWeight(1.5);
strokeWeight(1);
//  stroke(105, 105, 105);  // fill(0xff33ff99);
//   stroke(0xff33ff99);  // fill(0xff33ff99);  0xffff0000
stroke(0xffff0000);

// zoom control bar
fill(0, 0, 0, 0);

rect(minX,barHeight + height - height/4,maxX-minX, barY );
// rect(maxX-minX, barHeight/2,minX,barY + height - height/4 );

// slider
// fill(105, 105, 105); //0x3300FF00
fill(0xffff0000); // 0xff33ff99//0x3300FF00

rect(sliderX, sliderY + height - height/4 + sliderHeight/2 , sliderWidth, sliderHeight);

popMatrix();

}

float getZoomValue(int x, int y) {
if ((x >= minX) && (x <= (maxX - sliderWidth/2)) && (y > (height - height/3))) {
sliderX = (int) (x - (sliderWidth/2));
if (sliderX < minX) {
sliderX = minX;
}
sliderValue = (x - minX) * valuePerX + minZoomValue;
}
return sliderValue;
}

void updateZoomSlider(float value) {
int tempX = (int) (value / valuePerX) + minX;

if ( (tempX >= minX) && (tempX <= (maxX+sliderWidth))  ) {
sliderValue = value;
sliderX = tempX;
}
}

/*   boolean isZoomSliderEvent(int x, int y) {
int slop = 50;  // number of pixels above or below slider that's acceptable.  provided for ease of use.
int sliderTop = (int) (sliderY - (sliderHeight/2)) - slop;
int sliderBottom = sliderY + sliderHeight + slop;
return ((x >= sliderX) && (x <= (sliderX    + sliderWidth)) && (y >= sliderTop)  && (y <= sliderBottom) || draggingZoomSlider );
} */

boolean isZoomSliderEvent(int x, int y) {
int slop = 50;  // number of pixels above or below slider that's acceptable.  provided for ease of use.
int sliderLeft = (int) (sliderX - (sliderWidth/2)) - slop;
int sliderRight = sliderX + sliderWidth + slop;
//  return ((y >= sliderY + height - height/4) && (y <= (sliderY + height - height/4    + sliderHeight)) && (x >= sliderLeft)  && (x <= sliderRight) || draggingZoomSlider );
return ((y >= sliderY + height - height/4 - sliderHeight/2) && (y <= (sliderY + height - height/4 + sliderHeight*2 )) && (x >= sliderLeft )  && (x <= sliderRight ) || draggingZoomSlider );
}
}
``````

The best I can do right now:

https://github.com/HarryDrones/NACA4Series