gav = gravity_local(6400, 6400, 9.8, 1, 388400); %starting altitude from center of the earth, altitude of the gravity acceleration
%from center of the earth, gravity acceleration, altitude increment in kilometers, final altitude 
aw = flight_drag( 63.62, .045, 0, 1); %frontal srea m^2 of craft, drag coefiecient, starting altitude, altitude increments km)
vel = (2*aw/104400 + 2*gav/104400)^(1/2); %velocity m/s, craft mass equals 104400
acc = .5*(vel^2)/3100; %acceleration in m/s^2 based on velocity, 3100 is acceleration distance
gf = acc/9.8; % 1 g equals 9.8 m/s^2 this estimates g expirnienced by passengers
eng = energy_consumption(104400, acc); %104400=craft mass
kwh = eng / 3600000 * 1000; %watt=volt*amp, amp=coloum/seconds, one hour is 3600 seconds, 1000 volt is just an example
plcost = .093 * kwh; %energy costs 9.3 cents per kWh, estimated in dollars
display(gav);
display(aw);
display(vel);
display(acc);
display(eng);


function gav = gravity_local(current_altitude, starting_altitude, g1, inc, target_altitude) %determines work against gravity
    gav = 0;
   while current_altitude < target_altitude
    g2 = g1 / (current_altitude/starting_altitude)^2; %determines gravitational pull based on elevation, as they are proportional
    gav = 104400 * g2 * inc + gav; %determines the work done every kilometer then adds it to the total
    current_altitude = current_altitude + inc;
   end

end

function aw = flight_drag(frontal_area, drag_coefficient, bottom_elevation, altinc) %determines the work done agains drag
    aw = 0;
    gav = gravity_local((6400+bottom_elevation), 6400, 9.8, 1, 388400);
    while bottom_elevation < 16 %determines work done in troposphere
        vel = (2*aw/104400 + 2*(gav)/104400)^(1/2);
        time = altinc / vel;
        aw = 1/2 * drag_coefficient * .362 * vel * frontal_area * time + aw;
        bottom_elevation = bottom_elevation + altinc;
    end
    gav = gravity_local((6400+bottom_elevation), 6400, 9.8, 1, 388400);
    while (16 <= bottom_elevation) && ( bottom_elevation < 50)%determines work done in stratosphere
        vel = (2*aw/104400 + 2*(gav)/104400)^(1/2);
        time = altinc / vel;
        aw = 1/2 * drag_coefficient * .00108 * vel * frontal_area * time + aw;
        bottom_elevation = bottom_elevation + altinc;
    end
end

function eng = energy_consumption(mass, acceleration) %energy consumption based on coulombs formula
    eng = mass * acceleration * 12.566 * 9.0e+09 * .001^2;
end