C# code for 4 and 5 digit NACA airfoils

Update: The naca5 code is wrong, updated (python) code can be found on github
Just in case its useful for somebody. C# code to generate 4 or 5 digit NACA airfoils. It requires a cubic spline interpolation method which you can easily write yourself or find on google code (please respect the licenses though)

–Dirk

//4 digit NACA

	// Conversions
	double cl = chordLength.ConvertToDouble("m");
	double naca1 = int.Parse(number.Substring(0,1));
	double naca2 = int.Parse(number.Substring(1,1));
	double naca34 = int.Parse(number.Substring(2,2));

	double dl;
	double ml;
	double pl;

	ml = naca1/100;
	pl = naca2/10;
	dl = naca34/100;

	// Number of points
	int n = 100;
	int m =n/2+1;

	// Airfoil-Drop
	double a0 = 1.4845;
	double a1 = -0.63;
	double a2 = -1.7580;
	double a3 = 1.4215;
	double a4 = -0.5075;

	// Create the point objects
	Point2D pt = new Point2D();
	Point2D[] ptList = new Point2D[n];

	// Calculate the points ...
	double x;
	double y;
	double yt;
	double dy;
	double theta;
	double xu;
	double xl;
	double yu;
	double yl;

	// TrailingEdge
	pt.X=1*cl;
	pt.Y=0;
	ptList[0] = new Point2D(pt);

	// Between
	for (int i=1;i<(m-1);i++) {
		x=((double) (1.0-0.5*(1.0-Math.Cos(Math.PI*(1.0*i)/(1.0*m-1.0)))));
		yt=dl*(a0*Math.Sqrt(x)+x*(a1+x*(a2+x*(a3+x*a4))));

        if (x<=pl){
	             y = (ml/Math.Pow(pl,2))*(2.0*pl-x)*x;
	             dy = (ml/Math.Pow(pl,2))*(2.0*pl-2.0*x);
	    } else {
	             y = (ml/Math.Pow((1.0-pl),2))*((1.0-2.0*pl)+2.0*pl*x-Math.Pow(x,2));
	             dy = (ml/Math.Pow((1.0-pl),2))*(2.0*pl-2.0*x);
	    }
	    theta = Math.Atan(dy);
	    xu = x-yt*Math.Sin(theta);
	    yu = y+yt*Math.Cos(theta);
	    xl = x+yt*Math.Sin(theta);
	    yl = y-yt*Math.Cos(theta);

		//upper point
		pt.X = xu*cl;
		pt.Y = yu*cl;
		ptList[i] = new Point2D(pt);

		//lower point
		pt.X = xl*cl;
		pt.Y = yl*cl;
		ptList[n-i] = new Point2D(pt);
	}

	// LeadingEdge point
	pt.X=0.0;
	pt.Y=0.0;
	ptList[m-1] = new Point2D(pt);

	return ptList;

//5 digit NACA

	// Conversions
	double cl = chordLength.ConvertToDouble("m");
	double naca1 = int.Parse(number.Substring(0,1));
	double naca23 = int.Parse(number.Substring(1,2));
	double naca45 = int.Parse(number.Substring(3,2));

	double dl;
	double cld;
	double pl;
	double ml;
	double k1;

	cld = naca1*0.75/10;
	pl = 0.5*naca23/100;
	dl = naca45/100;

	// Number of points
	int n = 100;
	int m =n/2+1;

	// Airfoil-Drop
	double a0 = 1.4845;
	double a1 = -0.63;
	double a2 = -1.7580;
	double a3 = 1.4215;
	double a4 = -0.5075;

	// Create the point objects
	Point2D pt = new Point2D();
	Point2D[] ptList = new Point2D[n];

	// Calculate the points ...
	double x;
	double y;
	double yt;
	double dy;
	double theta;
	double xu;
	double xl;
	double yu;
	double yl;

	// TrailingEdge
	pt.X = 1*cl;
	pt.Y = 0.0;

	ptList[0] = new Point2D(pt);

	double[] P = new double[]{0.05,0.1,0.15,0.2,0.25};
	double[] M = new double[]{0.0580,0.1260,0.2025,0.2900,0.3910};
	double[] K = new double[]{361.4,51.64,15.957,6.643,3.230};

	double[] mll = Interpolate(P,M,new double[]{pl});
	ml = mll[0];

	double[] k1l = Interpolate(M,K,new double[]{ml});
	k1 = k1l[0];

	for (int i=1;i<(m-1);i++) {
		x=((double) (1.0-0.5*(1.0-Math.Cos(Math.PI*(1.0*i)/(1.0*m-1.0)))));
		yt=dl*(a0*Math.Sqrt(x)+x*(a1+x*(a2+x*(a3+x*a4))));

		if (x<=pl)
	        {
	             y = cld/0.3 * (1.0/6.0)*k1*( Math.Pow(x,3)-3.0*ml*Math.Pow(x,2)+ml*ml*(3.0-ml)*x );
	             dy = (1.0/6.0)*k1*( 3.0*Math.Pow(x,2)-6.0*ml*x+ml*ml*(3.0-ml) );
	        }
	        else
	        {
	             y = cld/0.3 * (1.0/6.0)*k1*Math.Pow(ml,3)*(1.0-x);
	             dy = -(1.0/6.0)*k1*Math.Pow(ml,3);
	        }

	        theta = Math.Atan(dy);
	        xu = x-yt*Math.Sin(theta);
	        yu = y+yt*Math.Cos(theta);
	        xl = x+yt*Math.Sin(theta);
	        yl = y-yt*Math.Cos(theta);

		//upper point
		pt.X = xu*cl;
		pt.Y = yu*cl;
		ptList[i] = new Point2D(pt);

		//lower point
		pt.X = xl*cl;
		pt.Y = yl*cl;
		ptList[n-i] = new Point2D(pt);
	}

	// LeadingEdge point
	pt.X=0.0;
	pt.Y=0.0;
	ptList[m-1] = new Point2D(pt);

	return ptList;