ballistic.py

from math import *

g = 9.81

mass = 1
density = 1.225
drag_coeff = 0.1
area = 1
zero_velocity = 100
zero_angle = radians(45)
ls = 30
modes = 1

def menu():
    while 1:
        print(\
    """Make your choice:
    1. Start simulation
    0. Exit"""
              )
        choice = -1
        while not (int(choice) in range(4)):
            choice = input("-->")
        choice = int(choice)
        elif (choice == 1):
            options()
        elif (choice == 0):
            return
def is_int(value):
    try:
        value = int(value)
    except:
        return 0
    return 1

def is_float(value):
    try:
        value = float(value)
    except:
        return
    return 1

def print_options_menu():
    strings = ['1','2']
    str_modes = ["Without air drag", "Laminar air drag"]
    for x in range(2):
        if (modes&2**x):
            strings[x] = "*"+str_modes[x]+"*"
        else:
            strings[x] = str_modes[x]
    print(\
"""
Body:
        1: m = {mass:.3f} kg -- mass
        2: v = {zero_velocity:.3f} m/s -- zero velocity
        3: a = {zero_angle:.3f} degrees -- zero angle
        4: S = {area:.3f} m^2 -- cross-sectional area
        5: C = {drag_coeff:.3f} -- drag coeff
Environment:
        6: g = {g:.3f}m/s^2 -- acceleration of gravity
        7: p = {density:.3f} kg/m^3 -- density of medium
Simulation:
        8: Field size = {ls} symbols
Operations:
        9:  {str0}
        10: {str1}
        11: Start simulation
        0: Exit
        """.format(\
            mass=mass,\
            zero_velocity=zero_velocity,\
            zero_angle=degrees(zero_angle),\
            area=area,\
            drag_coeff=drag_coeff,\
            g=g,\
            density=density,\
            ls=ls,\
            str0 = strings[0],\
            str1 = strings[1]))

def options():
    global modes
    while(1):
        the_num = -1
        drag = 0
        print_options_menu()
        print("Choose the variable number you wish to change")
        while(1):
            the_num = input("-->")
            if (not is_int(the_num)):
                print("Not an integer is entered, try again")
                continue
            the_num = int(the_num)
            if (the_num < 0 or the_num > 12):
                print("Incorrect value entered, try again")
                continue
            break
        if (the_num == 0): break
        if (the_num == 11):
            bal_sim(modes)
            continue
        elif (the_num >= 9):
            modes = modes^2**(the_num-9)
            continue
        print("Now enter a value you wish for chosen variable to be")

        rules_negative = [0, 0, 1, 0, 0, 0, 0, 0]
        rules_upto1 = [0, 0, 0, 0, 0, 0, 0, 0]
        rules_upto90 = [0, 0, 1, 0, 0, 0, 0, 0]
        rules_int = [0, 0, 0, 0, 0, 0, 0, 1]
        rules_positive = [1, 1, 0, 1, 1, 1, 1, 1]
        rules = ["mass", "zero_velocity", "zero_angle", "area", "drag_coeff", "g", "density", "ls"]

        print("Your value must be...")
        if(rules_negative[the_num-1]):
            print(" not negative")
        if(rules_upto1[the_num-1]):
            print(" be less or equal to 1 modulo")
        if(rules_upto90[the_num-1]):
            print(" be less or equal to 90 modulo")
        if(rules_positive[the_num-1]):
            print(" positive")
        the_value = 0
        while(1):
            the_value = input("-->")
            if (not is_float(the_value)):
                print("Not a number is entered, try again")
                continue
            the_value = float(the_value)
            is_incorrect = 0
            if (rules_negative[the_num-1] and the_value < 0):  is_incorrect=1
            if ((rules_upto1[the_num-1] and the_value > 1) or (rules_negative[the_num-1] and the_value < -1)): is_incorrect=1
            if ((rules_upto90[the_num-1] and the_value > 90) or (rules_negative[the_num-1] and the_value < -90)):   is_incorrect=1
            if (rules_int[the_num-1] and not(int(the_value) == the_value)):   is_incorrect=1
            if (rules_positive[the_num-1] and the_value <= 0):   is_incorrect=1
            if (is_incorrect):
                print("Incorrect value entered, try again")
                continue
            if (rules_upto90[the_num-1]):
                the_value = radians(the_value)
            if (rules_int[the_num-1]):
                the_value = int(the_value)
            break
        globals()[rules[the_num-1]] = the_value

def approximate_input(to_reach, error, constants, variable, formula, goes_down, default_value=0, RPV=0, var_max=0, fun_RPV=0, fun_min=0):
    var = default_value
    dvar = 1
    dm = 0
    current_vars = constants
    current_vars[variable] = var
    fx = eval(formula, globals(), current_vars)
    while 1:
        delta_var = var
        if (abs(fx-to_reach) <= error):
            return var
        if (fx < to_reach and not goes_down):
            if (dm==2):
                dvar = dvar / 2
            dm = 1
            if (RPV and var+dvar > var_max):
                var = var - dvar
            else:
                var = var + dvar
            current_vars[variable] = var
        elif (fx > to_reach and goes_down):
            if (dm==2):
                dvar = dvar / 2
            dm = 1
            if (RPV and var+dvar > var_max):
                var = var - dvar
            else:
                var = var + dvar
            current_vars[variable] = var
        elif (fx > to_reach and not goes_down):
            if (dm==1):
                dvar = dvar / 2
            dm = 2
            if (RPV and var+dvar > var_max):
                var = var + dvar
            else:
                var = var - dvar
            current_vars[variable] = var
        elif (fx < to_reach and goes_down):
            if (dm==1):
                dvar = dvar / 2
            dm = 2
            if (RPV and var+dvar > var_max):
                var = var + dvar
            else:
                var = var - dvar
            current_vars[variable] = var
        fx = eval(formula, globals(), current_vars)
        delta_var = delta_var - var
        delta_var = var - delta_var
        if (fun_RPV):
            if (fx < fun_min):
                var = var + delta_var


def lar1hfyt(t, k, vy): #Laminar Air drag First Half Y(t)
    kt = k*t
    return k**-2 * (g * (-kt - e**-kt + 1) + vy * (k - k*e**-kt))

def lar2hfyt(t, k, my):
    return (my*k**2 + g*k*t - g*e**(k*t)+g)/k**2

def larfxt (t, k, vx): #Laminar Air drag x(t)
    return (vx-vx*e**(-k*t))/k

def larftx (x, k, vx): #Laminar Air drag t(x)
    return (log(vx / (vx-x*k)) / k)

def tarfvt (t, k, vy):
    anti_gk = sqrt(k/g)
    gk = sqrt(g*k)
    antik2gk = sqrt(g/k)

    return -antik2gk * tan(t*gk - atan(vy*anti_gk))

def bal_sim(mode=1):
    vx = zero_velocity * cos(zero_angle)
    vy = zero_velocity * sin(zero_angle)
    the_trajectory_map = [[' ' for y in range(ls+1)] for x in range(ls+1)]
    info = dict()
    sd = 0
    if (mode&2**0):
        if (not g):
            # Special case of no gravity
            print("Not implemented for zero-gravity yet")
            return
        fly_time = 2*vy / g #Fly time
        mx = vx*fly_time
        my = (vy**2) / (2*g)
        true_my = (zero_velocity**2) / (2*g)
        sx = mx / ls
        sy = true_my / ls
        sd = max(sx, sy)
        for yy in range(ls+1):
            current_y_rounded = yy*sd

            D = vy**2 - 2*g*current_y_rounded
            if (D < 0): continue
            D = D**0.5
            t1= (vy - D) / g
            t2= (vy + D) / g

            corresponding_fx1 = (vx*t1)
            corresponding_fx1 = round(corresponding_fx1 / sd)
            corresponding_fx2 = (vx*t2)
            corresponding_fx2 = round(corresponding_fx2 / sd)
            if (0 <= corresponding_fx1 < ls+1):
                the_trajectory_map[yy][corresponding_fx1] = 'X'
            if (0 <= corresponding_fx2 < ls+1):
                the_trajectory_map[yy][corresponding_fx2] = 'X'
        for xx in range(ls+1):
            current_x_rounded = xx*sd
            t1 = current_x_rounded / vx
            corresponding_fy1 = (vy*t1 - (0.5*g*t1**2))
            corresponding_fy1 = round(corresponding_fy1 / sd)
            if (0 <= corresponding_fy1 < ls+1):
                the_trajectory_map[corresponding_fy1][xx] = 'X'
        info['arlmx'] = mx
        info['arlmy'] = my
        info['arltif'] = fly_time
    if (mode&2**1): #Laminar drag
        k = 0.5*area*drag_coeff*density/mass
        #y(t) = (g (-k t-e^(-k t)+1)+v_0 (k-k e^(-k t)))/k^2 -- first half, C is Vy0
        #y(t) = (C k^2+g k t-g e^(k t)+g)/k^2 -- second half, C is My
        #x(t) = (C-C e^(-k t))/k, C is Vx0
        myt = log(1 + k*vy / g) / k #Maximum Y Time
        my  = lar1hfyt(myt, k, vy) #Maximum Y
        warmy = 0.5 * vy**2 / g #Without Air Drag Maximum Y
        sy = warmy / ls
        sht = approximate_input(sy / 2, sy / 2, locals(), "t", "lar2hfyt(t, k, my)", 1)
        tif = sht + myt #Time in flight
        mx = (vx-vx*e**(-k*tif))/k
        warmx = vx*tif
        sx = warmx / ls
        if (not sd):
            sd = max(sx, sy)

        for xx in range(ls+1):
            cx = xx * sd
            if (cx > mx): continue
            t = larftx(cx, k, vx)
            y0 = 0
            if (t<=myt):
                y0 = lar1hfyt(t, k, vy)
            else:
                y0 = lar2hfyt(t-myt, k, my)
            y0 = round(y0 / sd)
            the_trajectory_map[y0][xx] = 'X'
        for yy in range(ls+1):
            cy = yy * sd
            if (cy > my):  continue
            t0 = approximate_input(cy, sd / 2, locals(), "t", "lar1hfyt(t, k, vy)", 0, 0, 1, myt)
            t1 = approximate_input(cy, sd / 2, locals(), "t", "lar2hfyt(t, k, my)", 1, 0, 1, tif)
            x0 = larfxt(t0, k, vx)
            x1 = larfxt(myt+t1, k, vx)
            x0 = round(x0 / sd)
            x1 = round(x1 / sd)
            the_trajectory_map[yy][x0] = 'X'
            the_trajectory_map[yy][x1] = 'X'
        true_mx = warmx
        true_my = warmy
        fly_time = tif
        info['larmx'] = mx
        info['larmy'] = my
        info['lartif'] = tif

    for yy in range(ls, -1, -1):
        for xx in range(ls+1):
            print(the_trajectory_map[yy][xx], end='')
        print(": {0:<} m".format(int(yy*sd)))
    if (mode&2**0):
        print("Dragless mode info:")
        print(" Maximum length: ", info['arlmx'], ' meters')
        print(" Maximum height: ", info['arlmy'], ' meters')
        print(" In flight for: ", info['arltif'], ' seconds')
    if (mode&2**1):
        print("Laminar drag mode info:")
        print(" Maximum length: ", info['larmx'], ' meters')
        print(" Maximum height: ", info['larmy'], ' meters')
        print(" In flight for: ", info['lartif'], ' seconds')

menu()