using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Numerics;

namespace XCraft
{
	struct ProfileCoefficients
	{
		public float Cl;
		public float Cd;
		public float Cm25;
		public float stallDevelop; // flow separation infromation (0 to 1)
	};

	class ProfileDescription
	{
		public float liftCurveSlope;
		public float zeroLiftAoA;
		public float maxCl;
		public float minCl;
		public float Cm25_0; //quarter chord pitching moment coeff. at zero lift angle
		public float Cd_0; //drag coeff at zero lift angle of atack
		public float stallTransitionAngleDelta;//angle delta needed for full stall develop when critical angle of atack is exceeded
		public float maxStallDragCoeff;//max pressure drag coefficient (experimental: 2.0 for NACA0012 Re 1800000)
		public float maxStallLiftCoeff; //max lift coefficent when stalling (experimental: 1.1 for NACA0012 Re 1800000)
		public float reverseEfficiency;//reverse lift coeff efficiency (experimental: 0.667 for NACA0012 Re 1800000)

		public ProfileDescription()
		{
			liftCurveSlope = (float)Math.PI * 2.0f;
			zeroLiftAoA = 0.0f;//symetric profile
			maxCl = 1.3f;
			minCl = -1.3f;
			Cm25_0 = 0.0f;
			Cd_0 = 0.007f;
			stallTransitionAngleDelta = 0.1745f;
			maxStallDragCoeff = 2.0f;
			maxStallLiftCoeff = 1.1f;
			reverseEfficiency = 0.667f;
		}

		public static ProfileDescription MakeDefault()
		{
			return new ProfileDescription();
		}
	}

	class Profile
	{
        private struct DirDescription
        {
            public float stallAnglePos;
            public float stallAngleNeg;
            public float stallTransitionAngleDeltaInv;
            public float Clmax;
            public float Clmin;
        };

        // direction data
        private DirDescription fwd;
        private DirDescription rev;
        
		private float liftCurveSlope;
		private float zeroLiftAngle;
		private float Cd_0;
		private float Cm25_0;
		private float maxStallLift;
		private float maxStallDrag;
        

        private ProfileCoefficients GetCoefficientsInternal(float angleRad, ref DirDescription dir)
		{
			float effAngle = angleRad - zeroLiftAngle;
			
			float stallDevelop = Math.Max(
							M.Clamp((angleRad - dir.stallAnglePos) * dir.stallTransitionAngleDeltaInv, 0.0f, 1.0f),
							M.Clamp((dir.stallAngleNeg - angleRad) * dir.stallTransitionAngleDeltaInv, 0.0f, 1.0f)
						);

			float Cl_stall = (float) Math.Sin(effAngle*2) * maxStallLift;
			stallDevelop = (1.0f - (float)Math.Cos(stallDevelop*Math.PI)) * 0.5f;
			float Cl_norm = M.Clamp((effAngle) * liftCurveSlope, dir.Clmin, dir.Clmax);

			ProfileCoefficients ret = new ProfileCoefficients();

			ret.Cl = M.Lerp(Cl_norm, Cl_stall, stallDevelop);			
			ret.Cd = Cd_0 + stallDevelop * maxStallDrag * (float) ((1.0 - Math.Cos(effAngle*2.0)) * 0.5f);
			ret.stallDevelop = stallDevelop;

			return ret;
		}
		
		/// Angle of atack in range <-PI .. PI> 
		/// axialBack is backward unit vector of wing section, normal is section up vector
		/// Vfree is velocity of fluid free stream (it's not necessary to be normalized)
		public static float AngleOfAtack(Vector2 Vfree, Vector2 normal)
		{		
			Vector2 v = Vector2.Normalize(Vfree);
			Vector2 back = new Vector2(-normal.Y, normal.X); 
			float d = M.Clamp(Vector2.Dot(v, back), -1.0f, 1.0f); //because of float inaccuracy (acos(1.000001) is NaN)
			float a = (float) Math.Acos(d);
			if (Vector2.Dot(normal, v) < 0.0f)
				a = -a;
			return a;
		}

		public ProfileCoefficients GetCoefficients(float angleRad)
		{
			float effAngle = angleRad - zeroLiftAngle;
			bool forward = Math.Abs(effAngle) < (float) (Math.PI * 0.5);

			ProfileCoefficients coeffs;

			if (forward)
				coeffs = GetCoefficientsInternal(angleRad, ref fwd);
			else
			{
				float angle2 = angleRad - (float) Math.PI;
				if (angle2 < -Math.PI)
					angle2 = angle2 + (float) Math.PI * 2.0f;

				coeffs = GetCoefficientsInternal(angle2, ref rev);
			}

			// Cm25 calculation from Cd and Cl
			//
			// forces point of action approximation is measured as distance (r) from quarter chord point
			float cosEffAngle = (float) Math.Cos(effAngle);
			float r = M.Lerp(
				(1.0f - ((cosEffAngle + 1.0f) * 0.5f)) * 0.5f,
				forward ? coeffs.stallDevelop * 0.25f : 0.5f - coeffs.stallDevelop * 0.25f,
				0.5f
				);

			float rD = (float) Math.Sin(effAngle) * r;
			float rL = cosEffAngle * r;
			float moment = rD * coeffs.Cd + rL * coeffs.Cl;

			float baseMoment = Cm25_0 * (1.0f - coeffs.stallDevelop);
			if (!forward)
				baseMoment = -baseMoment;
			coeffs.Cm25 = -moment + baseMoment;

			return coeffs;
		}
			
		public Profile(ProfileDescription pd)
		{
			fwd.stallAnglePos = pd.maxCl / pd.liftCurveSlope;
			fwd.stallAngleNeg = pd.minCl / pd.liftCurveSlope;
			fwd.stallTransitionAngleDeltaInv = 1.0f / pd.stallTransitionAngleDelta;
			fwd.Clmax = pd.maxCl;
			fwd.Clmin = pd.minCl;

			liftCurveSlope = pd.liftCurveSlope;
			zeroLiftAngle = pd.zeroLiftAoA;
			Cd_0 = pd.Cd_0;
			Cm25_0 = pd.Cm25_0;
			maxStallLift = pd.maxStallLiftCoeff;
			maxStallDrag = pd.maxStallDragCoeff;

			rev.Clmax = Math.Sign(pd.maxCl) * Math.Abs(pd.minCl) * pd.reverseEfficiency;
			rev.Clmin = Math.Sign(pd.minCl) * Math.Abs(pd.maxCl) * pd.reverseEfficiency;
			rev.stallTransitionAngleDeltaInv = (1.0f) / (pd.stallTransitionAngleDelta * pd.reverseEfficiency);
			rev.stallAnglePos = rev.Clmax / pd.liftCurveSlope;
			rev.stallAngleNeg = rev.Clmin / pd.liftCurveSlope;
		}
	}
}