2-D Kinematics with accurate air drag, Javascript

Language: Layout:

//Written by Max Gao

//all units are SI, and assumes flat earth
//INPUT
//some are later updated
var v_tot=27.73;//starting total velocity
var v_x=0;//starting x comp
var v_y=0;//starting y comp
var s_x=0;//start x displacement
var s_y=0;//start y displacement
var ts=0.1;//time step per calculation
var m=0.13//mass
var theta=45//launch angle from horizontal

var A=0.00423205605;//cross-sectional area of projectile
var Cd=0.9;//drag coeffecient of projectile
var time=300;//how long to run the simulation

//Calculation variables
var a_x=0;//net x acceleration
var a_y=0;//net y acceleration
var t=0;//time elapsed
var p=0;//atmospheric pressure
var T=0;//absolute temp
var L=0;//temperature lapse rate
var h=0;//minimum height of band
var b=0;//current band
var pt=1;//pressure trigger if falling down (positive if falling down)
var g=0;//gravitational acceleration
var pa_tot=0;//total drag accel
var pa_x=0;//atmospheric drag x acceleration
var pa_y=0;//atmospheric drag y acceleration

//components calculated once
v_x=Math.cos(theta*3.1415926/180)*v_tot;
v_y=Math.sin(theta*3.1415926/180)*v_tot;

print("time(s) ","x-accel(m/s2) ","y-accel(m/s2) ","x-vel(m/s) ","y-vel(m/s)","x-disp(m) ","y-disp(m) ","pressure(pa) ","x-drag-accel(m/s2) ","y-drag-accel(m/s2) ","g-accel(m/s2) ");

for(t=0; t<time; t=t+ts){
    if(0<=s_y){
    
    
    if(s_y<11000){//band 0
        b=0;
        T=288.15;
        h=0;
        L=-0.0065;
        p=101325*Math.pow((1-(L*s_y)/(T+L*s_y)), (0.0342/L));
    }
    if(11000<=s_y&&s_y<20000){//band 1
        b=1;
        T=216.65;
        h=11000;
        L=0;
        p=22632.1*Math.exp(-(0.0342*(s_y-h))/T);
    }
    if(20000<=s_y&&s_y<32000){//band 2
        b=2;
        T=216.65;
        h=20000;
        L=0.001;
        p=5474.89*Math.pow((1-(L*(s_y-h))/(T+L*(s_y-h))), (0.0342/L));
    }
    if(32000<=s_y&&s_y<47000){//band 3
        b=3;
        T=228.65;
        h=32000;
        L=0.0028;
        p=868.02*Math.pow((1-(L*(s_y-h))/(T+L*(s_y-h))), (0.0342/L));
    }
    if(47000<=s_y&&s_y<51000){//band 4
        b=4;
        T=270.65;
        h=47000;
        L=0;
        p=110.91*Math.exp(-(0.0342*(s_y-h))/T);
    }
    if(51000<=s_y&&s_y<71000){//band 5
        b=5;
        T=270.65;
        h=51000;
        L=-0.0028;
        p=66.94*Math.pow((1-(L*(s_y-h))/(T+L*(s_y-h))), (0.0342/L));
    }
    if(71000<=s_y&&s_y<86000){//band 6
        b=6;
        T=214.65;
        h=71000;
        L=-0.002;
        p=3.96*Math.pow((1-(L*(s_y-h))/(T+L*(s_y-h))), (0.0342/L));
    }
    if(86000<=s_y){//band 6
        b=7;
        T=386;
        p=0;
        h=s_y;
    }
    
    /*
    if(v_y<0){
        pt=-1;
    }
    if(v_y>=0){
        pt=1;
    }*/
        
        v_tot=Math.sqrt(v_x*v_x+v_y*v_y);//calculate total velocity 
    
        //check if component forces are conserved when dealing with atmospheric pressure.
        //https://en.wikipedia.org/wiki/Drag_equation
        //https://en.wikipedia.org/wiki/Barometric_formula#Density_equations
    g=-3.98E14/(s_y*s_y+1.28E8*s_y+4.096E13); //gravity. despite a flat earth, gravity peters out
    pa_tot=0//-((p*v_tot*v_tot*A*Cd*0.00174)/(T+L*(s_y-h)))/(m); //atmospheric pressure accel total (force/mass) with air density subbed over (density = pressure / (gas constant * temp))
    pa_x=pa_tot*(v_x/v_tot);
    pa_y=pa_tot*(v_y/v_tot);//component will be positive if y is negative
        //end of dicey zone
    
    a_x=pa_x;
    a_y=pa_y+g;
        
    print(t, a_x, a_y, v_x, v_y, s_x, s_y, p, pa_x, pa_y, g);//display numbers
	s_x=s_x+0.5*ts*ts*a_x+ts*v_x;//s=s0+ut+0.5at^2
    s_y=s_y+0.5*ts*ts*a_y+ts*v_y;
	v_x=ts*a_x+v_x;//a=dv/dt
    v_y=ts*a_y+v_y;
    }
    
    else{
        print("Landed!");
    }
}

[ + ] Show input