b = 20.6
c = 0.4125
T0=278
moho=40
z0=0
T(z) = T0 + b*(z-z0) - 0.5*c*(z-z0)**2
effic(z) = (T(z)-T(z0))/T(z)
work(z) = flux(z) * (T(z)-T(z0))/T(z)
## constant times...
flux(z) = b - c*(z-z0)
set xrange [0:moho]
plot T(x)-273
pause -1 'ready'
plot work(x)
pr "## optimum depthb is 15 km."
pr flux(15)/flux(0)
pr "there, flux is 0.7 of surface flux. Kelvin temp is 541"
pr T(15)
pr "Celsius: 268"
pr T(15)-273
pr "the flux in question is 50mW/m^2 at the surface"
pr "efficiency of ideal heat engine is 0.486 there"
pr effic(15)
pr "Here is the final answer in work done as mW/m^2: 17"
pr effic(15)* flux(15)/flux(0) * 50
pr "could conceivably have the cold sink at cooler temp than T(z0)=5?"