list engines in englist.
set g0 to 9.80665.

function getthrust{
	parameter englist.
	set sumthrust to 0.
	for eng in englist
		set sumthrust to sumthrust + eng:thrust.

	return sumthrust.
}

//not tested with multiple engines.
function getweightedisp{
	parameter englist.
	local liquid_density is .005.   //tonne/unit
	local solid_density is .0075.
	local mono_density is .004.

	set summf to 0.
	set ispmf to 0.

	for eng in englist{
		set massflow to 0.
		if eng:allowshutdown{
			if eng:name = "omsengine"
				set massflow to eng:fuelflow*mono_density.
			else if eng:name = "ionengine"
				. //ignoring ion engines right now. electricity is part of flow
			else set massflow to eng:fuelflow*liquid_density.
		}
		else{//srb
			set massflow to eng:fuelflow*solid_density.
		}
		set ispmf to ispmf + eng:isp * massflow.
		set summf to summf + massflow.
	}
	if summf
		return ispmf/summf.
	else
		return 0.
}

function getdensity{
	parameter alti.
	return 1.
}

function shipvec{
	parameter vector, veccolor, text, vscale is .5, data is vector:mag, precision is 2.

	VECDRAWARGS(ship:POSITION, vector*vscale, veccolor, text+round(vector:mag,precision), 1, TRUE).
}

function forcetest{
	parameter freq.
	//get initial data
	set t1 to time:seconds.
	set m1 to ship:mass.
	set v1 to ship:velocity:orbit.
	set thrust1 to ship:facing:forevector * getthrust(englist).

	wait freq.

	set t2 to time:seconds.
	set m2 to ship:mass.
	set v2 to ship:velocity:orbit.
	set thrust2 to ship:facing:forevector * getthrust(englist).

	//perform all calculations
	set wisp to getweightedisp(englist).
	set ev to -ship:facing:forevector * (wisp*g0).
	set accel to (v2-v1)/(t2-t1).
	set massflow to (m2-m1)/(t2-t1).
	set totalforce to m1*accel + ev*massflow - thrust2.
	set gravity to ship:BODY:POSITION:normalized *(body:mu/(altitude+body:radius)^2).
	set weight to m2*gravity.
	set drag to totalforce - thrust2 - weight.//set drag to the total force that is not weight or thrust. (assuming non changing thrust, inc gimbal...)
	set cdA to 0.
	if altitude < body:atm:height
		set CdA to (2*drag:mag)/(getdensity(altitude)*v2:mag^2).

	//print and draw stuff
	clearvecdraws().
	clearscreen.
	print "Isp: " + round(wisp,2).
	print "Drag: "+round(drag:mag,2).
	print " CdA: "+round(CdA,5).
	print "dynp: "+round(ship:dynamicpressure,3).
	shipvec(v2,red,"v2: ",.02).
	shipvec(accel,yellow, "Accel: ",3).
	shipvec(totalforce, green, "#Force: ").
	shipvec(weight, red, "Weight: ").
	shipvec(drag,yellow,"Drag: ").
	if abs(thrust2:mag) > .1
		shipvec(thrust2, blue, "Thrust: ").
	if abs(ev:mag) > .1
		shipvec(ev,white,"ExhVel: ",.01).
}

function main{
	until false
		forcetest(.25).
}
main().