Code for simulating the elytra item in Minecraft

ElytraModel.java

/**
* An accurate simulation of the elytra item as of Minecraft 15w41b
*/
public class ElytraModel
{	
	public int glideTime;
	public int damageTaken;
	
	public double posX;
	public double posY;
	public double posZ;
	
	public double velX;
	public double velY;
	public double velZ;
	
	/**
	* Simulates a Minecraft tick (20 per second).
	* The pitch and yaw are the look direction of the player.
	*/
	public void tick(boolean isCreative, double yaw, double pitch)
	{
		if (!isCreative && (this.glideTime + 1) % 20 == 0)
		{
			this.damageTaken++;
		}
		
		//I did some simplifing of the folowing to reduce the number of negatives and trig functions
		double yawcos = Math.cos(-yaw - Math.PI);
		double yawsin = Math.sin(-yaw - Math.PI);
		double pitchcos = Math.cos(pitch);
		double pitchsin = Math.sin(pitch);
		
		double lookX = yawsin * -pitchcos;
		double lookY = -pitchsin;
		double lookZ = yawcos * -pitchcos;
		
		double hvel = Math.sqrt(velX * velX + velZ * velZ);
		double hlook = pitchcos; //Math.sqrt(lookX * lookX + lookZ * lookZ)
		double sqrpitchcos = pitchcos * pitchcos; //In MC this is multiplied by Math.min(1.0, Math.sqrt(lookX * lookX + lookY * lookY + lookZ * lookZ) / 0.4), don't ask me why, it should always =1
		
		//From here on, the code is identical to the code found in net.minecraft.entity.EntityLivingBase.moveEntityWithHeading(float, float) or rq.g(float, float) in obfuscated 15w41b
		this.velY += -0.08 + sqrpitchcos * 0.06;
		
		if (this.velY < 0 && hlook > 0)
		{
			double yacc = this.velY * -0.1 * sqrpitchcos;
			this.velY += yacc;
			this.velX += lookX * yacc / hlook;
			this.velZ += lookZ * yacc / hlook;
		}
		if (pitch < 0)
		{
			double yacc = hvel * -pitchsin * 0.04;
			this.velY += yacc * 3.5;
			this.velX -= lookX * yacc / hlook;
			this.velZ -= lookZ * yacc / hlook;
		}
		if (hlook > 0)
		{
			this.velX += (lookX / hlook * hvel - this.velX) * 0.1;
			this.velZ += (lookZ / hlook * hvel - this.velZ) * 0.1;
		}
		
		this.velX *= 0.99;
		this.velY *= 0.98;
		this.velZ *= 0.99;
		
		this.posX += this.velX;
		this.posY += this.velY;
		this.posZ += this.velZ;
		
		this.glideTime++;
	}
	
	/** 
	* Checks if the player is currently in a gliding state.
	* As you can see, if the player is in creative, they will remain gliding even if on the ground. They will stop gliding once they move (but that functionality is not shown here).
	*/
	public boolean isGliding(boolean isCreative, boolean isOnGround, float fallDistance)
	{
		if (isCreative)
		{	
			return glideTime > 0;
		}
		else
		{
			return !isOnGround && fallDistance >= 1.0f;
		}
	}
}

Mathematica Elytra.txt

BuildElytraData[pos_, vel_] := <|"time" -> 0, "pos" -> pos, 
  "vel" -> vel|>
BuildElytraData[] := BuildElytraData[{0, 0, 0}, {0, 0, 0}]
ElytraTick[dat_, \[Theta]_, \[Phi]_] := 
  Module[{out = dat, vel, pitchcos, pitchsin, sqrpitchcos, look, hvel,
     hlook, yacc}, vel = dat[["vel"]]/20; pitchcos = Cos[\[Phi]]; 
   pitchsin = Sin[\[Phi]]; 
   look = {Sin[-\[Theta] - \[Pi]]*-pitchcos, -pitchsin, 
     Cos[-\[Theta] - \[Pi]]*-pitchcos};
   hvel = Sqrt[vel[[1]]*vel[[1]] + vel[[3]]*vel[[3]]];
   hlook = pitchcos;
   sqrpitchcos = pitchcos*pitchcos;
   vel[[2]] += -0.08 + sqrpitchcos*0.06;
   If[vel[[2]] < 0 && hlook > 0, yacc = vel[[2]]*-0.1*sqrpitchcos; 
    vel[[2]] += yacc; vel[[1]] += look[[1]]*yacc/hlook; 
    vel[[3]] += look[[3]]*yacc/hlook]; 
   If[\[Phi] < 0 && hlook > 0, yacc = hvel*-pitchsin*0.04; 
    vel[[2]] += yacc*3.5; vel[[1]] -= look[[1]]*yacc/hlook; 
    vel[[3]] -= look[[3]]*yacc/hlook]; 
   If[hlook > 0, vel[[1]] += (look[[1]]/hlook*hvel - vel[[1]])*0.1; 
    vel[[3]] += (look[[3]]/hlook*hvel - vel[[3]])*0.1];
   vel *= {0.99, 0.98, 0.99}; out["pos"] += vel; out["vel"] = vel*20; 
   out["time"] += 1/20; out];
ElytraSim[initdat_, \[Theta]f_, \[Phi]f_, time_] := 
 Module[{out = {initdat}, t, dat = initdat}, 
  While[dat[["time"]] <= time, t = dat[["time"]]; 
   dat = ElytraTick[dat, \[Theta]f[t], \[Phi]f[t]]; 
   AppendTo[out, dat]]; out]
PlotElytraFlight[datlist_, options_] := 
 Module[{points, x, y, z, rangeinfo, range, ranges}, 
  points = datlist[[All, "pos"]][[All, {1, 3, 2}]]; 
  rangeinfo = (range = MinMax[#]; 
      Join[range, {Mean[range], range[[2]] - range[[1]]}]) & /@ 
    Transpose[points]; 
  ListPointPlot3D @@ 
   Join[{points, BoxRatios -> (Max[#, 1] & /@ rangeinfo[[All, 4]]), 
     AxesLabel -> {"x", "z", "y"}}, options]]
FindFlightEq[initdat_, \[Theta]_, \[Phi]_, err_] := 
 Module[{lastdat, dat, con = True}, dat = initdat; 
  While[con, lastdat = dat; dat = ElytraTick[dat, \[Theta], \[Phi]]; 
   con = Max[Abs[dat[["vel"]] - lastdat[["vel"]]]] > err]; dat]
-----------------------------------------------------------------------
                               Examples
-----------------------------------------------------------------------

PlotElytraFlight[ElytraSim[BuildElytraData[{0, 0, 0}, {0, 0, 0}], 25 # Degree &, 70 Cos[#] Degree &, 20], {Filling -> Bottom, ImageSize -> Large}]

ListLinePlot[Transpose[Table[Join[{{1}, {1}}*180.*\[Phi]/\[Pi], Partition[FindFlightEq[BuildElytraData[], 0, \[Phi], 0.001][["vel"]][[{3, 2}]]*{1, -1}, 1], 2], {\[Phi], -\[Pi]/2, \[Pi]/2, \[Pi]/64}]], AxesLabel -> {"Pitch", None}, PlotLabel -> "Glide Velocity", ImageSize -> Large, PlotLegends -> {"x", "y"}]