#!/usr/bin/perl -w

# script for generating Ip to AFR table and Ip to O2% table

if (@ARGV == 1) {
    $cal = shift @ARGV;
} else {
    $cal = 170; # free air Ip current
}

$pump_pw_zero = 102;
$mult = 1024;


# pump OPA transfer function
# input: dutycycle, n = [0 .. 255]
# output: Ip [mA]
sub opa_transfer_function {
# Ip= 5V/1M * 61.9k / 61.9 Ohm * (command - 102) / 256
    my $cmd = shift;
    my $t = 5 / 1000000 * 61900 / 61.9;
    return $t * ($cmd - $pump_pw_zero) / 256;
}


# gnuplot curve fit of Bosch Ip->lambda table (section 4.4 in datasheet)
# curve fit ONLY VALID in the range lambda [0,6 ... 2,4]
# input: Ip [mA]
# output: lambda
sub curve_fit {
    my $x = shift;
    my $res = -0.0119369 * $x * $x * $x - 0.122866 * $x * $x - 0.331488 * $x + 0.329051 + 0.698059 * exp($x);
    return $res;
}

$last_idx=999;
# calculate ideal values
for ($i=0; $i<256; $i++) {
    $dc2ip[$i] = opa_transfer_function($i) * 1000; # *1000 : result in mA please!
    $lambda[$i] = curve_fit($dc2ip[$i]);

    if ($i > 144) {  #    if ($lambda[$i] >= 1.5) { # lambda 1.5 is max for the lambdacorr table
	$lambda[$i] = 0; # end of supported range
	if ($last_idx == 999) {
	    $last_idx = $i; # save the index for later
	}
    }

    $afr[$i] = $lambda[$i] * 14.7;
    $afr_mult[$i] = $afr[$i] * $mult;
}


# do a calibration calculation (simulate what the avr is doing)
for ($i=0; $i<256; $i++) {

    # first calculate the new effective Ip (measured Ip adjusted by calibration constant)
    my $ip_adjusted = ($i-$pump_pw_zero)/256 * (2*$cal + 1024) / 256;
    $ip_cal[$i] = $ip_adjusted;

    # then lookup the adjusted Ip in the Ip->lambda table
    if ($i < $last_idx) {  # table length is limited to lambda range [0,6 ... 1,5]
	$lambda_cal[$i] = curve_fit($ip_adjusted);
	$afr_cal[$i] = $lambda_cal[$i] * 14.7;
	$tblidx_cal[$i] = $ip_adjusted * 256/5 + $pump_pw_zero; # 256/5 ~= 51
    } else {
	$lambda_cal[$i] = 0;
	$afr_cal[$i] = 0;
	$tblidx_cal[$i] = 0;
    }
}

# interpolate the array (eg. do the same calculation as the AVR) and later display the deviance from ideal result
$step_size=8;
$last_element=144;

for ($i=0; $i<=$last_element; $i++) {
    my $m = $i % $step_size;
    my $previdx = $i-$m;
    my $nextidx = $i-$m+$step_size;

    if ($m != 0) {
	$afr_mult_key[$i] = 0;
	$afr_mult_interp[$i] = linterp($afr_mult[$previdx], $afr_mult[$nextidx], $m / $step_size);
    } else {
	$afr_mult_key[$i] = $afr_mult[$i];
	$afr_mult_interp[$i] = $afr_mult[$i];
    }
}

for ($i=$last_element+1; $i<256; $i++) {
    $afr_mult_key[$i] = 0;
    $afr_mult_interp[$i] = 0;
}


# now print the results
# left of the '|' is the ideal values
# right of the '|' is the calibrated values
# right of '||' is ip->afr*1024 table and table index for the calibrated ip
# right of '|||' is the interpolation as the AVR does and the error caused by interpolation

print "cal = $cal\n";
print "dc\tIp\tlambda\tafr\t  |\tIpC\tlambdaC\tafrC\t||\tafr1024\ttblidx\t|||\tinterp\tIdiff\n";
for ($i=0; $i<256; $i++) {
    if (($i % $step_size) != 0) {
	$key = ' ';
    } else {
	$key = '*';
    }

    printf("%d\t%0.3f\t%0.3f\t%0.3f\t  |\t%0.3f\t%0.3f\t%0.3f\t||\t%d\t%0.1f\t|||\t%s%0.1f\t%0.2f\n",
	   $i, $dc2ip[$i], $lambda[$i], $afr[$i], 
	   $ip_cal[$i], $lambda_cal[$i], $afr_cal[$i], int($afr_mult[$i]),
	   $tblidx_cal[$i],
	   $key, $afr_mult_interp[$i], $afr_mult[$i]-$afr_mult_interp[$i]
	   );
}

print "\nuint16_t PROGMEM ip_to_afr[] = {";
for ($i=0; $i<$last_element; $i+=$step_size) {
    my $val = int($afr_mult[$i]+0.5);
    print "$val, ";
}
my $val = int($afr_mult[$last_element]+0.5);
print "$val};\n";


sub print_gnuplot_curve_fit {
    print <<HOPLA;

# f(x) = a*x*x*x + b*x*x + c*x + d + e*exp(x)
# fit f(x) "transferfunc.dat" using 1:2 via a, b, c, d, e
# plot "transferfunc.dat" using 1:2 with points 6, f(x)

# Ip -> lambda for Bosch LSU4
-2.8	0.60 # GUESS!
-1.85	0.70
-1.08	0.80
-0.76	0.85
-0.47	0.90
0.00	1.009
0.34	1.18
0.68	1.43
0.95	1.70
1.40	2.42
#2.55	207
HOPLA
}


#------ linear interpolation
sub linterp
{
	my $a = shift;
	my $b = shift;
	my $offset = shift;	# range is 0 (a) .. 1 (b)
	return $offset * $b + (1-$offset) * $a;
}


####################################################################################
####################################################################################
####################################################################################
####################################################################################
####################################################################################

exit;

# Here is the Ip->O2% table generated

sub print_gnuplot_O2_curve_fit {
    print <<PINGPONG;
    
# f(x) = a*x*x + b*x + c
# fit f(x) "transferfunc.dat" using 1:2 via a, b, c
# plot "transferfunc.dat" using 1:2 with points 6, f(x)

# a = -0.331824
# b = 9.04668
# c = -0.0111669

# Ip -> O2% for Bosch LSU4

0.00	0.0
0.34	3.0
0.68	6.0
0.95	8.29
1.40	12.0
2.55	20.9
PINGPONG
}

# gnuplot curve fit of Bosch Ip->O2 table (section 4.4 in datasheet)
# input: Ip [mA]
# output: O2% * 100
sub curve_fit_IpO2 {
    my $x = shift;
    my $res = -0.331824 * $x * $x +  9.04668 * $x - 0.0111669;
    $res *= 100;
    return $res;
}

# input: index [0..255]
# output: O2% * 100
sub calc_O2 {
    my $ip = shift;
    my $t = opa_transfer_function($ip) * 1000;
    my $res = int(curve_fit_IpO2($t)/10 + 0.5);
    return $res;
}

# make a 6-element Ip->O2 table
print "\nuint8_t PROGMEM ip_to_o2[] = {";
$step_size = 32;
$steps = 6;

for ($i=0; $i<$steps * $step_size; $i+=$step_size) {
    my $idx = $i + $pump_pw_zero;
    my $val = calc_O2($idx);

    print "," unless $i == 0;
#    print " $idx, $val";
    print " $val";
}
print " };\n";
