main.c

#include <math.h>
#include <stdio.h>

#include <allegro.h>

#define ONE_DEG (M_PI / 180.0)

typedef struct surfdata
{
	double
		liftSlope,
		skinFriction,
		zeroLiftAoA,
		stallAngleLo, stallAngleHi,
		chord,
		flapFraction,
		aspectRatio;
} surfdata_t;

surfdata_t sd = {
	.liftSlope = 2.0 * M_PI,
	.skinFriction = 0.02,
	.zeroLiftAoA = 0.0,
	.stallAngleLo = -10.0 * ONE_DEG, .stallAngleHi = 10.0 * ONE_DEG,
	.chord = 1.0,
	.flapFraction = 0.4,
	.aspectRatio = 4.0,
};

double Lerp(double a, double b, double t)
{
	return a + (b - a) * t;
}

double Clamp01(double x)
{
	return (x > 1.0 ? 1.0 : (x < 0.0 ? 0.0 : x));
}

double Clamp(double x, double a, double b)
{
	return (x > b ? b : (x < a ? a : x));
}

double CalcFlapEffCorrection(double flapAngle)
{
	return Lerp(0.8, 0.4, (fabs(flapAngle) * (1.0/ONE_DEG) - 10.0) / 50.0);
}

double CalcLiftCoeffMaxFrac(double flapFraction)
{
	return Clamp01(1.0 - 0.5 * (flapFraction - 0.1) / 0.3);
}

double CalcTorqCoeffProp(double effAngle)
{
	return 0.25 - 0.175 * (1.0 - 2.0 * fabs(effAngle) / M_PI);
}

double FrictionAt90Deg(double flapAngle)
{
	return 1.98 - 4.26e-2 * flapAngle * flapAngle + 2.1e-1 * flapAngle;
}

typedef struct coeffparams
{
	double
		zeroLiftAoA,
		corrLiftSlope,
		stallAngleHi,
		stallAngleLo,
		flapAngle;
} coeffparams_t;

typedef struct vec3
{
	double x, y, z;
} vec3;

vec3 vec3_Lerp(const vec3 a, const vec3 b, double t)
{
	return (vec3){
		a.x + (b.x - a.x) * t,
		a.y + (b.y - a.y) * t,
		a.z + (b.z - a.z) * t,
	};
}

vec3 CalcCoeffsLowAoA(const surfdata_t *sd, const coeffparams_t p, double angleOfAttack)
{
	double liftCoeff = p.corrLiftSlope * (angleOfAttack - p.zeroLiftAoA);
	double inducedAngle = liftCoeff / (M_PI * sd->aspectRatio);
	double effAngle = angleOfAttack - p.zeroLiftAoA - inducedAngle;
	
	double tanCoeff = sd->skinFriction * cos(effAngle);
	
	double normalCoeff = (liftCoeff + sin(effAngle) * tanCoeff) / cos(effAngle);
	double dragCoeff = normalCoeff * sin(effAngle) + tanCoeff * cos(effAngle);
	double torqueCoeff = -normalCoeff * CalcTorqCoeffProp(effAngle);
	
	return (vec3){liftCoeff, dragCoeff, torqueCoeff};
}

vec3 CalcCoeffsStalled(const surfdata_t *sd, const coeffparams_t p, double angleOfAttack)
{
	double liftCoeffLowAoA;
	
	if(angleOfAttack > p.stallAngleHi)
		liftCoeffLowAoA = p.corrLiftSlope * (p.stallAngleHi - p.zeroLiftAoA);
	else
		liftCoeffLowAoA = p.corrLiftSlope * (p.stallAngleLo - p.zeroLiftAoA);
	
	double inducedAngle = liftCoeffLowAoA / (M_PI * sd->aspectRatio);
	
	double lerpParam;
	if(angleOfAttack > p.stallAngleHi)
		lerpParam = ( M_PI / 2.0 - Clamp(angleOfAttack, -M_PI / 2.0, M_PI / 2.0)) / ( M_PI / 2.0 - p.stallAngleHi);
	else
		lerpParam = (-M_PI / 2.0 - Clamp(angleOfAttack, -M_PI / 2.0, M_PI / 2.0)) / (-M_PI / 2.0 - p.stallAngleLo);
	
	inducedAngle = Lerp(0, inducedAngle, lerpParam);
	double effAngle = angleOfAttack - p.zeroLiftAoA - inducedAngle;
	
	double normalCoeff = FrictionAt90Deg(p.flapAngle) * sin(effAngle) * (
			1.0 / (0.56 + 0.44 * fabs(sin(effAngle))) - 0.41 * (1.0 - exp(-17.0 / sd->aspectRatio))
		);
		
	double tanCoeff = 0.5 * sd->skinFriction * cos(effAngle);
	
	double liftCoeff = normalCoeff * cos(effAngle) - tanCoeff * sin(effAngle);
	double dragCoeff = normalCoeff * sin(effAngle) + tanCoeff * cos(effAngle);
	double torqueCoeff = -normalCoeff * CalcTorqCoeffProp(effAngle);
	
	return (vec3){liftCoeff, dragCoeff, torqueCoeff};
}

vec3 CalcCoeffs(const surfdata_t *sd, double angleOfAttack, double flapAngle)
{
	coeffparams_t p;
	p.flapAngle = flapAngle;
	p.corrLiftSlope = sd->liftSlope * (sd->aspectRatio / (sd->aspectRatio + 2.0 * (sd->aspectRatio + 4.0) / (sd->aspectRatio + 2.0)));
	
	double theta = acos(2.0 * sd->flapFraction - 1.0);
	double flapEff = 1.0 - (theta - sin(theta)) / M_PI;
	double deltaLift = p.corrLiftSlope * flapEff * CalcFlapEffCorrection(p.flapAngle) * p.flapAngle;
	
	double zeroLiftAoABase = sd->zeroLiftAoA;
	p.zeroLiftAoA = zeroLiftAoABase - deltaLift / p.corrLiftSlope;
	
	double stallAngleHiBase = sd->stallAngleHi;
	double stallAngleLoBase = sd->stallAngleLo;
	
	double clMaxHi = p.corrLiftSlope * (stallAngleHiBase - zeroLiftAoABase) + deltaLift * CalcLiftCoeffMaxFrac(sd->flapFraction);
	double clMaxLo = p.corrLiftSlope * (stallAngleLoBase - zeroLiftAoABase) + deltaLift * CalcLiftCoeffMaxFrac(sd->flapFraction);
	
	p.stallAngleHi = p.zeroLiftAoA + clMaxHi / p.corrLiftSlope;
	p.stallAngleLo = p.zeroLiftAoA + clMaxLo / p.corrLiftSlope;
	
	double padAngleHi = ONE_DEG * Lerp(15.0, 5.0, ( (1.0/ONE_DEG) * flapAngle + 50.0) / 100.0);
	double padAngleLo = ONE_DEG * Lerp(15.0, 5.0, (-(1.0/ONE_DEG) * flapAngle + 50.0) / 100.0);
	double padStallAngleHi = p.stallAngleHi + padAngleHi;
	double padStallAngleLo = p.stallAngleLo - padAngleLo;
	
	/*
	if(angleOfAttack < p.stallAngleHi && angleOfAttack > p.stallAngleLo){
		return CalcCoeffsLowAoA(sd, p, angleOfAttack);
	}
	
	if(angleOfAttack > padStallAngleHi || angleOfAttack < padStallAngleLo){
		return CalcCoeffsStalled(sd, p, angleOfAttack);
	}
	
	vec3 coeffsLow, coeffsStall;
	double lerpParam;
	if(angleOfAttack > p.stallAngleHi){
		coeffsLow = CalcCoeffsLowAoA(sd, p, p.stallAngleHi);
		coeffsStall = CalcCoeffsStalled(sd, p, padStallAngleHi);
		lerpParam = (angleOfAttack - p.stallAngleHi) / (padStallAngleHi - p.stallAngleHi);
	}
	else{
		coeffsLow = CalcCoeffsLowAoA(sd, p, p.stallAngleLo);
		coeffsStall = CalcCoeffsStalled(sd, p, padStallAngleLo);
		lerpParam = (angleOfAttack - p.stallAngleLo) / (padStallAngleLo - p.stallAngleLo);
	}
	
	return vec3_Lerp(coeffsLow, coeffsStall, lerpParam);
	*/
	
	vec3 coeffs;

	if(angleOfAttack < p.stallAngleHi && angleOfAttack > p.stallAngleLo){
		// Low angle of attack mode.
		coeffs = CalcCoeffsLowAoA(sd, p, angleOfAttack);
	}
	else{
		if(angleOfAttack > padStallAngleHi || angleOfAttack < padStallAngleLo){
			// Stall mode.
			coeffs = CalcCoeffsStalled(sd, p, angleOfAttack);
		}
		else{
			// Linear stitching in-between stall and low angles of attack modes.
			vec3 coeffsLow, coeffsStall;
			double lerpParam;
			
			if(angleOfAttack > p.stallAngleHi){
				coeffsLow = CalcCoeffsLowAoA(sd, p, p.stallAngleHi);
				coeffsStall = CalcCoeffsStalled(sd, p, padStallAngleHi);
				lerpParam = (angleOfAttack - p.stallAngleHi) / (padStallAngleHi - p.stallAngleHi);
			}
			else{
				coeffsLow = CalcCoeffsLowAoA(sd, p, p.stallAngleLo);
				coeffsStall = CalcCoeffsStalled(sd, p, padStallAngleLo);
				lerpParam = (angleOfAttack - p.stallAngleLo) / (padStallAngleLo - p.stallAngleLo);
			}
			coeffs = vec3_Lerp(coeffsLow, coeffsStall, lerpParam);
		}
	}
	return coeffs;
}

int main(int argc,char **argv)
{
	allegro_init();
	set_color_depth(32);
	set_gfx_mode(GFX_AUTODETECT_WINDOWED, 800, 600, 0, 0);
	install_keyboard();
	install_mouse();
	
	BITMAP *db = create_bitmap(SCREEN_W, SCREEN_H), *db2 = create_bitmap(SCREEN_W, SCREEN_H);
	clear_to_color(db, 0xFFFFFF);
	
	BITMAP *chart[3] = {
		load_bitmap("cl.pcx", 0),
		load_bitmap("cd.pcx", 0),
		load_bitmap("cm.pcx", 0),
	};
	
	for(int i = 0; i < 3; i++){
		BITMAP *tmp = chart[i];
		chart[i] = create_bitmap(db->w, db->h / 3);
		stretch_blit(tmp, chart[i], 0, 0, tmp->w, tmp->h, 0, 0, chart[i]->w, chart[i]->h);
		destroy_bitmap(tmp);
	}
	
	for(int i = 0; i < 8; i++)
		vline(db, i * db->w / 8, 0, db->h - 1, 0xAAAAAA);
	hline(db, 0, 100, db->w - 1, 0xAAAAAA);
	hline(db, 0, 200, db->w - 1, 0xAAAAAA);
	hline(db, 0, 200+200/3, db->w - 1, 0xAAAAAA);
	hline(db, 0, 200+2*200/3, db->w - 1, 0xAAAAAA);
	hline(db, 0, 400, db->w - 1, 0xAAAAAA);
	hline(db, 0, 500, db->w - 1, 0xAAAAAA);
	
	for(int j = 0; j < 3; j++){
		//FILE *f = fopen(fn[j], "w");
		for(double angle = -180.0 * ONE_DEG; angle < 180.0 * ONE_DEG; angle += ONE_DEG){
			double flapAngle[] = {
				0.0,
				25.0 * ONE_DEG,
				50.0 * ONE_DEG,
				-25.0 * ONE_DEG,
				-50.0 * ONE_DEG,
			};
			//fprintf(f, "%f", angle * (1.0/ONE_DEG));
			for(int i = 0; i < 5; i++){
				union {
					double comp[3];
					vec3 coeffs;
				} u;
				u.coeffs = CalcCoeffs(&sd, angle, flapAngle[i]);
				
				double limits[3][2] = {
					{-1.6, 1.6},
					{0, 1.5},
					{-0.6, 0.6},
				};
				double v = (u.comp[j] - limits[j][0]) / (limits[j][1] - limits[j][0]);
				int col[5] = {
					0x000000,
					0xFF00,
					0xFF0000,
					0xFFFF,
					0xFF
				};
				int h3rd = db->h / 3;
				set_clip_rect(db, 0, j * h3rd, db->w, (j + 1) * h3rd - 1);
				circle(db, (angle * (1.0/ONE_DEG) + 180.0)*(db->w / 360.0), (j + 1) * h3rd - v * h3rd, 2, col[i]);
				
				//fprintf(f, ";%f", u.comp[j]);
			}
			//fprintf(f,"\n");
		}
		//fclose(f);
	}
	while(!key[KEY_ESC]){
		set_trans_blender(255,255,255,(256 * mouse_x)/SCREEN_W);
		draw_sprite(db2, db, 0, 0);
		for(int i = 0; i < 3; i++){
			draw_trans_sprite(db2, chart[i], 0, i * db->h/3);
		}
		draw_sprite(screen, db2, 0, 0);
	}
	return 0;
}
END_OF_MAIN();