diff --git a/library/lib_exec_node.ks b/library/lib_exec_node.ks
new file mode 100644
index 0000000..4615904
--- /dev/null
+++ b/library/lib_exec_node.ks
@@ -0,0 +1,50 @@
+@lazyglobal off.
+
+// This file is currently just a placeholder for a more robust version,
+// uploaded solely because it is a very basic function for KSP.
+// I would be glad if anyone helps to develop it (add staging,
+// perhaps improve burntime calculations, add special cases such
+// as not enough fuel, etc.)
+
+function MaxShipThrust{ 
+	// this function is needed as a workaround
+	// to kOS bug regarding ship thrust.
+	local engs is list().
+	list engines in engs.
+	local mth is 0.
+	for eng in engs{
+		if eng:VISP > 0 {
+			set mth to mth + eng:MAXTHRUST.
+		}
+	}
+	return mth.
+}
+
+function exec_node{
+	parameter nd.
+	
+	local mth is MaxShipThrust().
+	lock burntime to nd:deltav:mag/(mth/mass).
+
+	lock steering to nd.
+	wait until nd:eta<burntime/2.
+	local timeout is time:seconds+max(burntime*2,100). // just a guess
+
+	lock throttle to cos(vang(nd:deltav,ship:facing:vector)).
+	// Thanks to cosine, we burn only when we face the correct direction.
+
+	wait until burntime<5 or time:seconds>timeout. 
+	// last 5 seconds of burn should be slower
+	local scale is nd:deltav:mag.
+
+	lock throttle to nd:deltav:mag/scale*
+		cos(vang(nd:deltav,ship:facing:vector)).
+	// We burn slower and slower.
+	
+	local timeout is time:seconds+100.
+	wait until nd:deltav:mag<scale*0.01 or time:seconds>timeout.
+	lock throttle to 0.
+
+	unlock throttle.
+	unlock steering.
+}
diff --git a/library/lib_orbital_maneuvers.ks b/library/lib_orbital_maneuvers.ks
new file mode 100644
index 0000000..2b7959e
--- /dev/null
+++ b/library/lib_orbital_maneuvers.ks
@@ -0,0 +1,259 @@
+@lazyglobal off.
+
+// This file contains a few functions that concentrate on creating
+// maneuver nodes completing most common KSP tasks.
+
+// This function creates a node (but doesn't place it) at time t, with
+// specified delta-v vector.
+function make_node_t_deltav{
+	parameter
+		t,
+		dv.
+	
+	local pro_unit is velocityat(ship,t):orbit:normalized.
+	local rad_unit is -VXCL(pro_unit,body:position-positionat(ship,t))
+		:normalized.
+	local nor_unit is VCRS(pro_unit,rad_unit).
+	return node(t,rad_unit*dv,nor_unit*dv,pro_unit*dv).
+}
+
+// This function gets time as a parameter, and returns a maneuver
+// node representing velocity change leading to circular orbit
+// with that burn happening at specified time. E.g. if you want to 
+// circularize at apoapsis, use: 
+// circularize_at_time(time:seconds+eta:apoapsis).
+function circularize_at_time{
+	parameter t.
+
+	local r is (body:position-positionat(ship,t)).
+	local actual_vel is velocityat(ship,t):orbit.
+	local expected_vel is VXCL(r,actual_vel):normalized*sqrt(body:mu/r:mag).
+	local dv is expected_vel-actual_vel.
+	return make_node_t_deltav(t,dv).
+}
+
+// This function gets two arguments: orbit and true anomaly in degrees
+// and then returns time remaining until orbiting object will pass 
+// through that true anomaly
+function time_to_true_anomaly{
+	parameter
+		orb,
+		a2.
+	
+	local e is orb:eccentricity.
+	local a1 is orb:trueanomaly. // the current one
+	local f is sqrt((1-e)/(1+e)).
+
+	local e1 is 2*arctan(f*tan(a1/2)).
+	local e2 is 2*arctan(f*tan(a2/2)).
+	// e1 and e2 are eccentric anomalies at these times in degrees
+	local m1 is constant():pi/180*e1-e*sin(e1).
+	local m2 is constant():pi/180*e2-e*sin(e2).
+	// m1 and m2 are mean anomalies at these times in radians
+	local n is 2*constant():pi/orb:period.
+	// n is mean angular velocity
+	local t1 is m1/n.
+	local t2 is m2/n.
+	// t1 and t2 are times with regard to some, non-disclosed epoch,
+	// at which orbitable will be at true anomaly a1 or a2 respectively
+	local diff is t2-t1.
+	if diff<0{ // ETA must be positive, so switch to next orbit
+		set diff to orb:period+diff.
+	}
+	else if diff>orb:period{ // we can do one full orbit less
+		set diff to diff-orb:period.
+	}
+	return diff.
+}
+
+// This function takes two orbits CONTAINED IN THE SAME PLANE and returns
+// a list of at most three elements:
+// - first one is a "status" - if it's equal to -1, then obt1 is fully
+// contained within obt2, vice versa for +1. Otherwise, it is equal to
+// 0 and the list contains two more elements
+// - two other elements, existing only if orbits cross at some points,
+// represent the time remaining until ship orbiting obt1 will pass common
+// orbit point. This should be the point at which you would do the burn
+// to randez-vous
+function get_orbit_intersections{
+	parameter
+		obt1,
+		obt2.
+
+	local e1 is obt1:eccentricity.
+	local e2 is obt2:eccentricity.
+	local a1 is obt1:semimajoraxis.
+	local a2 is obt2:semimajoraxis.
+	local B is a1/a2*(1-e1*e1)/(1-e2*e2).
+	local delta is obt2:argumentofperiapsis-obt1:argumentofperiapsis.
+	local gamma is arctan2(e1-B*e2*cos(delta), -B*e2*sin(delta)).
+	local S is (1-B)*cos(gamma)/B/e2/sin(delta).
+	if S > 1-0.00001{
+		return list(1).
+	}
+	else if S < -1+0.00001{
+		return list(-1).
+	}
+	local beta1 is arcsin(S)-gamma.
+	local beta2 is 180-arcsin(S)-gamma.
+	return list(
+		0, 
+		time_to_true_anomaly(obt1,beta1),
+		time_to_true_anomaly(obt1,beta2)
+	).
+}
+
+// This function will put a randez-vous node at orbit intersection 
+// between your orbit and tgt's orbit, if it exists. It assumes orbits
+// share common plane.
+function match_orbits{
+	parameter tgt.
+	
+	local l is get_orbit_intersections(obt,tgt:obt).
+	if l[0]<>0{
+		return 0.
+	}
+	local t is l[1].
+	if l[2]<t{
+		set t to l[2].
+	}
+
+	local actual_vel is velocityat(ship,time:seconds+t):orbit.
+	local intersection_r is body:position-positionat(ship,time:seconds+t).
+	
+	local tgt_r is body:position-tgt:position.
+	local tgt_normal is VCRS(tgt:velocity:orbit,tgt_r).
+	local side is VCRS(tgt_normal,tgt_r).
+	local true_anomaly is VANG(intersection_r,tgt_r).
+	if VDOT(side,intersection_r)<0{
+		set true_anomaly to -true_anomaly.
+	}
+	set true_anomaly to true_anomaly+tgt:obt:trueanomaly.
+	local t2 is time_to_true_anomaly(tgt:obt,true_anomaly).
+
+	local expected_vel is velocityat(tgt,time:seconds+t2):orbit.
+	return make_node_t_deltav(time:seconds+t,expected_vel-actual_vel).
+}
+
+// This function takes an orbit as an argument. It should be IDENTICAL to
+// your ship's orbit, other than phase (the other ship or body can be in
+// different position in the orbit). The function will then return a node
+// that in N orbits (where N is integer argument) will match your position.
+function match_orbital_phase{
+	parameter
+		orb,
+		n.
+	
+	local t1 is time_to_true_anomaly(obt,0).
+	local t2 is time_to_true_anomaly(orb,0).
+	
+	local period is obt:period.
+	local dt is mod(t2-t1+period*1.5,period)-period/2.
+	// now dt is time difference [-period/2,period/2] between target's
+	// and my periapsis
+	local h is 2*obt:semimajoraxis*(1+dt/n/period)^(2/3)-
+		obt:periapsis-2*obt:body:radius.
+	return change_opposite_height(time:seconds+t1,h).
+}
+
+// This function returns a node that will at specified time raise or lower
+// opposite orbit's side to specified height. If called at Pe, it will
+// simply change Ap to specified height and vice versa.
+function change_opposite_height{
+	parameter
+		t,
+		height.
+	
+	local r is (body:position-positionat(ship,t)).
+	local actual_vel is velocityat(ship,t):orbit.
+	local expected_vel is VXCL(r,actual_vel):normalized*sqrt(
+		2*body:mu*(height+body:radius)/r:mag/(r:mag+height+body:radius)
+	).
+	return make_node_t_deltav(t,expected_vel-actual_vel).
+}
+
+// This function rotates ship's orbit by burning normally at specified time
+// BY given amomunt (not TO angle, but BY angle).
+function change_inclination{
+	parameter
+		t,
+		ang.
+	
+	local r is (body:position-positionat(ship,t)).
+	local actual_vel is velocityat(ship,t):orbit.
+	local normal_vel is VCRS(actual_vel,r):normalized*actual_vel:mag.
+	return make_node_t_deltav(t,actual_vel*(cos(ang)-1)+normal_vel*sin(ang)).
+}
+
+// This function returns a maneuver node aligning orbital planes to
+// orbitable passed as an argument.
+function align_orbital_plane{
+	parameter other.
+
+	local my_r is (body:position-ship:position).
+	local my_normal is VCRS(obt:velocity:orbit,my_r).
+	
+	local other_r is (body:position-other:position).
+	local other_normal is VCRS(other:obt:velocity:orbit,other_r).
+
+	local target_place is VCRS(my_normal,other_normal):normalized.
+	local pe_place is positionat(ship,time:seconds+eta:periapsis)-body:position.
+	local my_place is ship:position-body:position.
+
+	local vel_at_pe is velocityat(ship,time:seconds+eta:periapsis):orbit.
+
+	local true_anomaly is VANG(target_place,pe_place).
+	if VDOT(target_place,vel_at_pe)<0{
+		set true_anomaly to -true_anomaly.
+	}
+	local true_anomaly_of_me is VANG(my_place,pe_place).
+	if VDOT(my_place,vel_at_pe)<0{
+		set true_anomaly_of_me to -true_anomaly_of_me.
+	}
+	local diff is mod(true_anomaly-true_anomaly_of_me+360,360).
+	if diff>180{
+		set true_anomaly to true_anomaly-180.
+		// switch descending to ascending node or vice versa to closer one
+	}
+	local t is time_to_true_anomaly(obt,true_anomaly).
+
+	local inc is VANG(my_normal,other_normal).
+	local vel_at_t is velocityat(ship,time:seconds+t):orbit.
+	if VDOT(other_normal,vel_at_t)>0{
+		set inc to -inc.
+	}
+
+	return change_inclination(time:seconds+t,inc).
+}
+
+// This function takes as argument an orbit and a height above ground.
+// It returns a list of at most three elements:
+// - first one is a status - if -1, then the orbit's apoapsis is too low,
+// if +1, orbit's periapsis is too high, otherwise it's zero
+// - the next two elements contain time remaining until orbiting object
+// passes through specified height
+function time_to_pass_height{
+	parameter
+		orb,
+		h.
+	
+	set h to h+orb:body:radius.
+	local e is orb:eccentricity.
+	local p is orb:semimajoraxis*(1-e*e).
+	if e=0{
+		set e to 0.0000001.
+	}
+	local c is (h-p)/e/h.
+	if c>1{
+		return list(-1).
+	}
+	if c<-1{
+		return list(1).
+	}
+	local theta is 180-arccos(c).
+	return list(
+		0,
+		time_to_true_anomaly(orb,theta),
+		time_to_true_anomaly(orb,-theta)
+	).
+}
diff --git a/unit_tests/test_lib_orbital_maneuvers.ks b/unit_tests/test_lib_orbital_maneuvers.ks
new file mode 100644
index 0000000..6fc56df
--- /dev/null
+++ b/unit_tests/test_lib_orbital_maneuvers.ks
@@ -0,0 +1,173 @@
+@lazyglobal off.
+
+run lib_orbital_maneuvers.
+run lib_exec.
+
+// Please run this test while being in an orbit around Kerbin, which has
+// Pe>30km (which is any stable orbit), Ap<30Mm (inside Mun's orbit is 
+// enough). Preferably, to test weird scenarios, make your orbit
+// inclined and eccentric.
+
+function om_unit_test1{
+	// This test tries to circularize at many different points of the orbit.
+	local i is 0.
+	until i>obt:period{
+		local nd is circularize_at_time(time:seconds+i).
+		add nd.
+		local deviation is round(nd:orbit:apoapsis-nd:orbit:periapsis,1).
+		if deviation>100{ // a reasonable margin
+			return "Orbit was not circularized - (Ap-Pe)="+deviation+"m>100m".
+		}
+		wait 0.01.
+		remove nd.
+		set i to i+obt:period/100.
+	}
+	return "".
+}
+
+function om_unit_test2{
+	// This test will try to predict ETA:periapsis and ETA:apoapsis via
+	// time_to_true_anomaly
+	local deviation is ETA:periapsis-time_to_true_anomaly(obt,0).
+	if abs(deviation)>5{
+		return "Difference between ETA:periapsis and time predicted by"+
+			"library is "+deviation+"s>5s".
+	}
+	local deviation is ETA:apoapsis-time_to_true_anomaly(obt,180).
+	if abs(deviation)>5{
+		return "Difference between ETA:apoapsis and time predicted by"+
+			"library is "+deviation+"s>5s".
+	}
+	// Note that this test might fail even though it works correctly, if
+	// you are very close to Pe or Ap, so that ETA fluctuates between 0
+	// and obt:period.
+	return "".
+}
+
+function om_unit_test3{
+	// This test will try to set its periapsis at exactly 30km at different
+	// positions in orbit.
+	local i is 0.
+	until i>obt:period{
+		local nd is change_opposite_height(time:seconds+i,30000).
+		add nd.
+		local deviation is round(nd:orbit:periapsis-30000,1).
+		if abs(deviation)>5{ // a reasonable margin
+			return "The script tried to set its periapsis at 30km, but"+
+				" instead set it at "+nd:orbit:periapsis/1000+"km".
+		}
+		wait 0.01.
+		remove nd.
+		set i to i+obt:period/100.
+	}
+	local i is 0.
+	until i>obt:period{
+		local nd is change_opposite_height(time:seconds+i,30000000).
+		add nd.
+		local deviation is round(nd:orbit:apoapsis-30000000,1).
+		if abs(deviation)>1000{ // a reasonable margin
+			return "The script tried to set its periapsis at 30Mm, but"+
+				" instead set it at "+nd:orbit:apoapsis/1000000+"Mm".
+		}
+		wait 0.01.
+		remove nd.
+		set i to i+obt:period/100.
+	}
+	return "".
+}
+
+function om_unit_test4{
+	local result is true.
+	set target to body("Minmus").
+	local nd is align_orbital_plane(target).
+	add nd.
+	local minmus_vel is target:velocity:orbit.
+	local minmus_pos is target:position-body:position.
+	local minmus_normal is VCRS(minmus_vel,minmus_pos).
+
+	local expected_vel is nd:orbit:velocity:orbit.
+	local expected_pos is nd:orbit:position-body:position.
+	local expected_normal is VCRS(expected_vel,expected_pos).
+
+	local vel_before_burn is velocityat(ship,nd:eta):orbit.
+	
+	local deviation is abs(obt:period-nd:orbit:period).
+	if deviation>10{ // 10s of margin
+		return "Plane switch maneuver changed orbital period (deviation="+
+			deviation+"s)".
+	}
+	local ang is VANG(expected_normal,minmus_normal).
+	if ang>2{ // two degrees of margin
+		return "Plane switch maneuver did not place at the same plane - "+
+			"deviation of planes is "+ang+" degrees, which is more than 2.".
+	}
+	if nd:eta>obt:period/2+20{
+		return "Plane switch maneuver was set at wrong inclination node.".
+	}
+	wait 1.
+	remove nd.
+	return "".
+}
+
+function om_unit_test5{
+	local result is true.
+	local ap is obt:apoapsis.
+	local pe is obt:periapsis.
+	local res is time_to_pass_height(obt,pe-1000).
+	if res[0]<>1{
+		return "Found a moment when height=periapsis-1000 (impossible)".
+	}
+	local res is time_to_pass_height(obt,ap+1000).
+	if res[0]<>-1{
+		return "Found a moment when height=apoapsis+1000 (impossible)".
+	}
+
+	local res is time_to_pass_height(obt,(ap+pe)/2).
+	if res[0]<>0{
+		return "Haven't found a moment when height=SMA".
+	}
+	local res is time_to_pass_height(obt,pe+1).
+	if res[0]<>0{
+		return "Haven't found a moment when height=Pe+1m".
+	}
+	else{
+		// This test may give false positive depending on eccentricity
+		// of your orbit. Try to test using eccentricity > 0.1
+		local deviation is abs(res[1]-res[2]).
+		if deviation>obt:period/100{
+			return "Two moments when height=Pe+1m were too separated"+
+				" in time ("+deviation+"s)".
+		}
+	}
+	return "".
+}
+
+function om_unit_tests{
+	local fail_reason is "".
+	local res is "".
+	local i is 1.
+	until i>5{
+		set res to evaluate("om_unit_test"+i+"()").
+		if res<>""{
+			set fail_reason to res.
+			break.
+		}
+		wait 1.
+		set i to i+1.
+	}
+
+	print " ".
+	print "Overall result:".
+	if fail_reason=""{
+		print "All tests passed SUCCESSFULLY".
+	}
+	else{
+		print "Some tests FAILED".
+		print "Reasone of the first FAIL:".
+		print "".
+		print fail_reason.
+	}
+	print " ".
+}
+
+om_unit_tests().