Processing to draw NACA 4 digit airfoils

[HackinHarry]

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();

}

[HackinHarry]

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

[kfrajer]

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

[HackinHarry]

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.

[HackinHarry]

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.

[HackinHarry]

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);      
  
             load_file();
        M = map(M,-0.10,0.10,-264.00,264.00); //,-0.10,0.10);
        T = map(T,-12,12,-tzoom,tzoom);
     }    
           // MousePress - Rotation Adjustment
  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);
  thread( "getOrdinates");

  rotate(radians(8));

        // rotate(radians(dT));
  shape(airfoiltest(xU,yU), 0, 0);
  //rotate(radians(-dT));
  rotate(radians(-17));

  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;
}

void load_file() {
  
coordinates = loadStrings("Airfoil.csv");  // also could use loadtable
//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) {
  
 table.addColumn(nameOfData, Table.FLOAT);
 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++) {
   row = table.addRow();
 }
}

void getOrdinates(){
   
for (int i = 0; i < N; i++){
  
  M = zoom/1000f;
  T = tzoom/1000f;
  
//  coordinates = loadStrings(csvfile);  // also could use loadtable 
//      String[] ordinate = splitTokens(coordinates[i], " ");
  if(i >= 0){
       beta[i] = (radians(180)/81.0) * i;
       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