import math
g = -9.807
pos = [0, 1.3]
vel = [0, 0]
accel = [0, g]
force = [0, 0]
wing_area = 16*23/10000
mass = 0.01452 #kg
e = 0.01 #time increment
k = 50 #spring constant of rubber band
dx = 0.15

def launch(k, mass, dx):
    return dx*math.sqrt(k/mass)

def alpha(vel):
    return -math.atan(vel[1]/vel[0])
    
def c_lift(vel):
    angle = 180*alpha(vel)/math.pi
    if angle <= 6:
        return 0.5*angle/6
    elif (angle > 6) and (angle < 14):
        return -0.19444444444*(angle**2) + 0.3888888888888*angle - 1.244444
    else:
        return 0.3

def lift_f(vel, wing_area):
    cl = c_lift(vel) #lift coefficient
    speed = math.sqrt(vel[0]**2 + vel[1]**2)
    return 0.5*cl*1.2*(speed**2)*wing_area

def drag_f(vel, wing_area):
    dc = 0.17 #drag coeffcient
    speed = math.sqrt(vel[0]**2 + vel[1]**2)
    return 0.5*dc*1.2*(speed**2)*wing_area

def g_force(mass, g):
    return mass*g

def forces(mass, g, vel, wing_area):
    angle = alpha(vel)
    weight = g_force(mass, g)
    lift = lift_f(vel, wing_area)
    drag = drag_f(vel, wing_area)
    force[0] = -drag*math.cos(angle)
    force[1] = lift + drag*math.sin(angle) + weight
    return force

def motion(mass, accel, vel, pos, e, g, wing_area):
    force = forces(mass, g, vel, wing_area)
    accel[0] = force[0]/mass
    accel[1] = force[1]/mass
    vel[0] += e*accel[0]
    vel[1] += e*accel[1]
    pos[0] += e*vel[0]
    pos[1] += e*vel[1]
    return pos, vel

vel[0] = launch(k, mass, dx)
t = 0
while True:
    pos, vel = motion(mass, accel, vel, pos, e, g, wing_area)
    t += e
    if pos[1] <= 0:
        print(pos)
        break