The climb rate depends on the excess power which is available after drag has been subtracted from net thrust. If the airplane stays at the same polar point while climbing, it needs to accelerate in order to compensate for the decrease in air density. Therefore, besides drag also this acceleration work needs to be subtracted before the remaining thrust can be used for climbing.
First let's clarify terms:
x$_g$, y$_g$, z$_g$ : Earth-fixed coordinate system
x$_f$, y$_f$, z$_f$ : Airplane-fixed coordinate system
x$_k$, y$_k$, z$_k$ : Kinetic coordinate system where x is the direction of movement
L$\;\;$ : Lift
D$\;\;$ : Drag
T$\;\;$ : Thrust
m$\;\:$ : mass
$\alpha\;\;$ : Angle of attack (between the x-axes of the airplane-fixed and kinetic coordinate systems)
$\gamma\;\;$ : Flight path angle (between the x-axes of the earth-fixed and kinetic coordinate systems)
$\sigma\;\:$ : Thrust angle relative to the airplane-fixed coordinate system
$v_{\infty}$ : Airspeed
The polar point should be the one for optimum climb speed. There is also one for optimum climb angle, but this simplification is justified. It also helps to make the math easier, since propeller aircraft climb best at the polar point where minimum power is required to maintain flight. This is at
$$c_L = \sqrt{3\cdot c_{D0}\cdot AR\cdot\pi\cdot\epsilon}$$
with
$c_L\;\;$: Lift coefficient
$c_{D0}$ : Zero-lift drag coefficient
$AR$ : Wing aspect ratio
$\epsilon\;\;$ : Wing efficiency factor
The zero-lift drag coefficient of propeller aircraft is around 0.025 to 0.04, with the high value for fixed-gear aircraft and the lower for those with retractable gear. It increases slightly with altitude due to the decrease of the Reynolds number from the drop in temperature. Here you need to pick a value which is appropriate for each specific aircraft.
Staying at the same polar point also means that weight will influence only the speed at which the aircraft climbs best, not the lift coefficient. The speed $v$ will change with the square root of the weight difference, because
$$v = \sqrt{\frac{m\cdot g}{\frac{\rho}{2}\cdot S_{ref}\cdot c_L}}$$
with $S_{ref}$ being the reference area of the aircraft and $\rho$ the air density.
Next to the correction term $C$ for acceleration. It depends on the local speed of sound, the gas constant for humid air $R_h$ and the temperature gradient (lapse rate $\Gamma$) of the atmosphere. This answer explains in detail how it is calculated and I repeat here only the result for standard atmospheric conditions:
$$C = 1 - 0.13335\cdot Ma^2 + \frac{(1+0.2\cdot Ma^2)^{3.5}-1}{(1+0.2\cdot Ma^2)^{2.5}}$$
with $Ma$ being the ratio between flight speed and local speed of sound.
Now your climb speed $v_z$ becomes $$v_z = \frac{v}{C}\cdot sin\gamma = \frac{v}{C}\cdot\frac{T\cdot cos(\sigma)-D}{m\cdot g} = \frac{P\cdot\eta_{prop}\cdot cos(\sigma) - D\cdot v}{C\cdot m\cdot g}$$
with $\eta_{Prop}$ the propeller efficiency and $P$ the engine brake power at the given altitude and throttle setting.
This leaves a bunch of unknown variables in order to correctly calculate the climb rate:
- engine power
- aircraft zero-lift drag coefficient
- propeller efficiency
Therefore, it will be best to look up the possible climb speeds at several altitudes and power settings from each POH and to interpolate between those values. Or you settle for an approximation and use rule-of-thumb values for the unknown parameters.
- for $\epsilon$ assume 0.8
- for $\sigma$ assume zero
- for $c_{D0}$ assume 0.026 at low and 0.03 at high altitude for retracted gear and 0.035 at low and 0.04 at high altitude for fixed gear.
- for $D$ use $\left(c_{D0} + \frac{c_L^2}{AR\cdot\pi\cdot\epsilon} \right) \cdot\frac{\rho\cdot v^2\cdot S_{ref}}{2}$
- for $\eta_{Prop}$ use 0.75 for a fixed-pitch and 0.8 for a constant speed prop.
- for normally aspirated engines reduce power proportionally with density. For turbocharged engines assume constant power up to their critical height and reduce power in proportion to density above that. Let the users of your program set the throttle setting themselves.
Where you have performance charts available, compare your results with published figures and tweak the variables such that you get a good fit. For example, look at the published optimum climb speed and adjust $c_{D0}$ until your result, taken from the optimum lift coefficient, agrees. And so on. This should give you very useable results.
...and find some known values for climb-rate for given pressures/temps/weights and use curve-fitting to come up with simple formulae which are still within a few percent of what I would get with these charts.
