#!/bin/sh
#this program calculates the number of calories consumed per mile travelled on
#a bicycle.  based on a program by ken roberts of columbia university.
#peter chen, 1995

E=0.249		#efficiency of human in cycling
T=0.95		#transmission efficiency of bicycle drivetrain
A1=0.0		#linear coefficient of air resistance
C_a=0.9		#Air resistance coefficient
BM_rate=1.4	#Basal Metabolism rate [Watts/kg body weight]
G=0.00		#grade of the hill (vertical_rise/odometer_distance)
HEAD=0.0	#velocity of headwind [miles/hour]
CROSS=1		#should take into account a true angelular speed
SPEED=12.0	#MPH
TF=60		#temperature in F
WM=30		#weight of bicycle (lb)
WC=165		#weight of cyclist (lb)
AC=0.4925155	#frontal area (m^2) for upright
		#0.4297982 crouch 0.3080527 racing 0.1844627 in pack
R=0.012		#rolling friction, mountain bike tires
		#0.004 racing slick 0.0051 27x1.25

while getopts a:c:g:h:m:r:s:t:x i
do
	case $i in
	a)	AC=$OPTARG;;
	c)	WC=$OPTARG;;
	g)	G=$OPTARG;;
	h)	HEAD=$OPTARG;;
	m)	WM=$OPTARG;;
	r)	R=$OPTARG;;
	s)	SPEED=$OPTARG;;
	t)	TF=$OPTARG;;
	x)	CROSS=.7;;	#Crosswinds are about 70% that of Headwinds
	?)	echo "usage: $0 [-a rontal_area] [-c cyclist_weight] [-g grade] [-h headwind] [-m machine_weight] [-r rolling_resistence] [-s speed] [-t temperature] [-x]" >&2
		exit 1;;
	esac
done
shift `expr $OPTIND - 1`

awk "
BEGIN	{
#mOs_per_kmOhr=1000/3600	#meters/second per kilometers/hour
mOs_per_miOhr=0.44704		#meters/second per miles/hour 
Nt_per_lb=4.4482
kg_per_Nt=0.102
Watts_per_CalOhr=1.163		#Watts per dietary Calories/hour
				#1 dietary Calorie == 1000 calories
#Watts_per_horsepower=745.700
air_density[1] = 1.5147		#-30C 3.6f
air_density[2] = 1.5083
air_density[3] = 1.5019
air_density[4] = 1.4955
air_density[5] = 1.4892
air_density[6] = 1.4829		#-25
air_density[7] = 1.4767
air_density[8] = 1.4706
air_density[9] = 1.4645
air_density[10] = 1.4584
air_density[11] = 1.3951	#-20
air_density[12] = 1.3896
air_density[13] = 1.3841
air_density[14] = 1.3787
air_density[15] = 1.3734
air_density[16] = 1.3680	#-15
air_density[17] = 1.3628
air_density[18] = 1.3575
air_density[19] = 1.3523
air_density[20] = 1.3472
air_density[21] = 1.3420	#-10
air_density[22] = 1.3370
air_density[23] = 1.3319
air_density[24] = 1.3269
air_density[25] = 1.3219
air_density[26] = 1.3170	#-5
air_density[27] = 1.3121
air_density[28] = 1.3072
air_density[29] = 1.3024
air_density[30] = 1.2977
air_density[31] = 1.2929	#0
air_density[32] = 1.2882
air_density[33] = 1.2835
air_density[34] = 1.2789
air_density[35] = 1.2742
air_density[36] = 1.2697	#5
air_density[37] = 1.2651
air_density[38] = 1.2606
air_density[39] = 1.2561
air_density[40] = 1.2517
air_density[41] = 1.2472	#10
air_density[42] = 1.2428
air_density[43] = 1.2385
air_density[44] = 1.2342
air_density[45] = 1.2299
air_density[46] = 1.2256	#15
air_density[47] = 1.2214
air_density[48] = 1.2171
air_density[49] = 1.2130
air_density[50] = 1.2088
air_density[51] = 1.2047	#20
air_density[52] = 1.2006
air_density[53] = 1.1965
air_density[54] = 1.1925
air_density[55] = 1.1885
air_density[56] = 1.1845	#25
air_density[57] = 1.1805
air_density[58] = 1.1766
air_density[59] = 1.1727
air_density[60] = 1.1688
air_density[61] = 1.1649	#30
air_density[62] = 1.1611
air_density[63] = 1.1573
air_density[64] = 1.1535
air_density[65] = 1.1498
air_density[66] = 1.1460	#35
air_density[67] = 1.1423
air_density[68] = 1.1387
air_density[69] = 1.1350
air_density[70] = 1.1314
air_density[71] = 1.1277	#40C 138.6f
air_density[72] = 1.1242
air_density[73] = 1.1206
air_density[74] = 1.1170
air_density[75] = 1.1135
}	END	{
T_a = (5.0/9.0)*($TF-32)	#convert temp F to C
V = mOs_per_miOhr * $SPEED	#velocity (miles/second) (speed)
Wc = Nt_per_lb * $WC		#weight of cyclist [Newtons]
Wm = Nt_per_lb * $WM		#weight of machine and clothing [Newtons]
W = Wc + Wm			#total weight
B = $BM_rate * kg_per_Nt * Wc;	#Create total coloric rate
D_a = air_density[int(T_a)+31]	#humidity and presure
A2 = ($C_a * D_a / 2) * $AC	#quadratic coefficient of air resistance
H = mOs_per_miOhr * $HEAD	#headwind
H *= $CROSS			#crosswind
F_a = A2 * (V+H) * (V+H) + $A1 * (V+H)	#force of air resistance [Newtons]
if((V+H) < 0)	F_a *= -1	#Calculate the force (+/-) of the air
F_r  =  $R * W			#Calculate the force or rolling restance
F_g  =  $G * W			#Calculate the force (+/-) of the grade
F  =  F_a + F_r + F_g		#Calculate the total force
P  =  V * F / $T		#Calculate Power in Watts
if (P > 0)			#Calculate Calories and drivetrain loss
	C  =  P / $E  +  B
else	C  =  B
#printf\"W=%d speed=$SPEED P=%f C=%f V=%f\n\",W/Nt_per_lb,P,C,V
#printf\"F_a=%f F_r=%f F_g=%f F=%f\n\",F_a,F_r,F_g,F
#printf\"R=$R A2=%f D_a=%s AC=$AC T_a=%f\n\",A2,D_a,T_a
printf\"%f\n\",(C * mOs_per_miOhr)/(V * Watts_per_CalOhr)}
" /dev/null