How to calculate force of air resistance in 2D projectile motion - Computational Science in C#

Question

I am writing a program that will output a table of various stats on a 100th of a second interval of throwing a ball at a specific angle with a specific initial speed. I am having trouble calculating the air resistance. I have my Code below. Level 1 is solved, witch is modeling the 2D projectile motion without air resistance. I am stuck on level 2 witch is the same thing, but I need to factor in air resistance.

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

namespace Homework2ProjectileMotion
{
class Program
{
    private const double interval = 0.001;
    private const int mass = 5;
    private const double initialSpeed = 10;
    private const double angle = Math.PI / 4;
    private const double startTime = 0;

    private void CalcAll(StreamWriter sw, double drag)
    {
        double speed = initialSpeed;
        double xPos = 0;
        double pXPos = 0;
        double yPos = 0;
        double pYPos = 0;
        double zPos = 0;
        double pZPos = 0;
        double distance = 0;
        double xVel = 10 * Math.Cos(angle);
        double yVel = 0;
        double zVel = 10 * Math.Sin(angle);
        double xAcc = 0;
        double yAcc = 0;
        double zAcc = -9.8;
        double mAcc = 9.8;
        double xAir = 0;
        double yAir = 0;
        double zAir = 0;
        double air = 0;
        double slope = 0;
        double vector = 0;
        sw.WriteLine("Time\tx\ty\tz\tDistance\tvx\tvy\tvz\tSpeed\tax\tay\taz\tm_Acc\n" + startTime + "\t" + xPos + "\t" + yPos + "\t" + zPos + "\t" + distance + "\t" + xVel + "\t" + yVel + "\t" + zVel + "\t" + speed + "\t" + xAcc + "\t" + yAcc + "\t" + zAcc + "\t" + mAcc);
        for (double i=startTime+interval;zPos>=0;i+=interval)
        {
            air = drag * speed * speed;
            if (xPos == 0)
            {
                slope = 1;
            }
            else
            {
                slope = (zPos - pZPos) / (xPos - pXPos);

            }
            vector = Math.Atan(slope);
            xAir = air * Math.Cos(vector);
            zAir = air * Math.Sin(vector);
            xAcc = ((xAcc * mass) - xAir) / mass;
            yAcc = ((yAcc * mass) - yAir) / mass;
            zAcc = ((zAcc * mass) - zAir) / mass;
            mAcc = Math.Sqrt(xAcc * xAcc + zAcc * zAcc);
            xVel = xVel + (xAcc * interval);
            yVel = yVel + (yAcc * interval);
            zVel = zVel + (zAcc * interval);
            speed = Math.Sqrt(xVel * xVel + zVel * zVel);
            pXPos = xPos;
            pYPos = yPos;
            pZPos = zPos;
            xPos = xPos + xVel * interval;
            yPos = yPos + yVel * interval;
            zPos = zPos + zVel * interval;
            distance = Math.Sqrt(xPos * xPos + zPos * zPos);
            sw.WriteLine(i + "\t" + xPos + "\t" + yPos + "\t" + zPos + "\t" + distance + "\t" + xVel + "\t" + yVel + "\t" + zVel + "\t" + speed + "\t" + xAcc + "\t" + yAcc + "\t" + zAcc + "\t" + mAcc);
        }
        sw.Close();
    }

    private void Level1()
    {
        StreamWriter writer = File.CreateText("Level1.txt");
        CalcAll(writer, 0);
        writer.Close();
    }
    private void Level2(double dragCoe)
    {
        StreamWriter writer = File.CreateText("Level2.txt");
        CalcAll(writer, dragCoe);
        writer.Close();
    }
    static void Main(string[] args)
    {
        Program run = new Program();
        //run.Level1();
        run.Level2(0.3);
    }
}

}

Here is an image of the first few outputs and a graph of the horizontal vs the vertical:

Answer 1

The formula I use to calculate drag is using vectors:

Fd = -.5*rho*|V|.*V*Cd*A;

or:

Vt = sqrt((2*m*g)/(Cd*A*rho)); 
y = Vt^2/(2*g)*ln((V0^2+Vt^2)/((V^2*Vt^2));
x = (Vt^2/g)*ln((Vt^2+g*U0*t)/Vt^2);

Where,

Vt is the terminal velocity

g is gravity

m is the mass

U0 is the initial x velocity (initialSpeed*cos(launchAngle))

V0 is the initial vertical velocity

rho is the air density

V is the velocity vector

|V| is the absolute value of the velocity vector

Cd is the drag coefficient (0.47 for spheres)

A is the cross sectional area (radius of sphere^2*pi)

Keep in mind, Fd is a force:

acceleration = (Force of gravity + force of drag )/ mass

because:

F = ma
a = F/m
Fnet = Fd+Fg;
a = Fnet/m

Answer 2

Your acceleration computation is wrong. The differential equation is

x'' = -k*v*x'
z'' = -k*v*y' - g

where v=speed is the length of the velocity vector.

Thus you should compute

speed = Math.Sqrt(xVel * xVel + zVel * zVel);
xAcc = (xGrav*mass-drag*speed*xVel)/mass;
zAcc = (zGrav*mass-drag*speed*zVel)/mass;

with constants xGrav=0 and zGrav=-9.81 (y as x) and then as you did with the Euler step.