Sous Vide PID controller

#include <science.h>
#include <avr/io.h>
#include <avr/interrupt.h>
 
Temperature t( PIN3 );
enum { heater = 9, zerocrossdetector = PIN2 };
 
Stats stats;
 
enum
{
  SAMPLE_TIME_IN_MS = 300,
  MILLISECONDS_PER_SECOND = 1000,
  OUTMIN = 0,
  OUTMAX = 255,
  PULSE = 4
};
 
class PID
{
  double kp;                  // Proportional coefficient
  double ki;                  // Integral coefficient
  double kd;                  // Derivative coefficient
  double previous_input;
  unsigned long lastTime;
  long output;
 
public:
 
  double Integral;
  double Target;
 
  PID( double Kp, double Ki, double Kd, double target )
  {
    Integral = 0;
    previous_input = target;
    output = 0;
 
    Target = target;
 
    double seconds_per_sample = double( SAMPLE_TIME_IN_MS )
      / MILLISECONDS_PER_SECOND;
 
    kp = Kp;
    ki = Ki * seconds_per_sample;
    kd = Kd / seconds_per_sample;
 
    lastTime = millis() - SAMPLE_TIME_IN_MS;        
  }
 
  int
  out( double input )
  {
    unsigned long now = millis();
    unsigned long timeChange = now - lastTime;
 
    if( timeChange >= SAMPLE_TIME_IN_MS )
    {
      double error = Target - input;
 
      Integral += ki * error;
 
      double Proportional = kp * error;
      double Derivative = kd * (previous_input - input);
 
      previous_input = input;
 
      if( Integral > OUTMAX )
        Integral = OUTMAX;
      else if( Integral < OUTMIN )
        Integral = OUTMIN;
 
      output = long( Proportional + Integral + Derivative );
 
      if( output > OUTMAX )
        output = OUTMAX;
      else if( output < OUTMIN )
        output = OUTMIN;
 
      lastTime = now;
    }
    return output;
  }
};
 
unsigned long int crossings = 0;
unsigned long int compares = 0;
unsigned long int overflows = 0;
 
void
zeroCrossingInterrupt()
  TCCR1B = 0x04;               // start timer with divide by 256 input
  TCNT1 = 0;                   // reset timer - count from zero
  crossings++;
}
 
ISR(TIMER1_COMPA_vect)
{
  digitalWrite( heater, HIGH );  // set triac gate to high
  TCNT1 = 65536 - PULSE;       // trigger pulse width
  compares++;
}
 
ISR(TIMER1_OVF_vect)
{
  digitalWrite( heater, LOW ); // turn off triac gate
  TCCR1B = 0x00;             // disable timer stop unintended triggers
  overflows++;
}
 
PID pid( 8, .5, 40, 133 );
 
void
setup()
{
  Serial.begin( 115200 );
 
  pinMode( zerocrossdetector, INPUT );        // zero cross detect
  digitalWrite( zerocrossdetector, HIGH );    // enable pull-up resistor
  pinMode( heater, OUTPUT );                  // triac gate control
  digitalWrite( heater, LOW );
 
  //
  // set up Timer1 
  // (see ATMEGA 328 data sheet pg 134 for more details)
  //
  OCR1A = 100;      // initialize the comparator
  TIMSK1 = 0x03;    // enable comparator A and overflow interrupts
  TCCR1A = 0x00;    // timer control registers set for
  TCCR1B = 0x00;    // normal operation, timer disabled
 
  attachInterrupt( 0, zeroCrossingInterrupt, RISING );    
  // IRQ0 is pin 2. Call zeroCrossingInterrupt on rising signal
}
 
void
leading_spaces( double num )
{
  if( num < 10 )
    Serial.print( "  " );
  else if( num < 100 )
    Serial.print( " " );
  Serial.print( num );
}
 
void
loop()
{
  int num = t.how_many_sensors();
 
  if( num <= 0 )
  {
    Serial.println( "No sensors" );
    delay( 2000 );
    return;
  }
  
  t.read( 0 );
 
  double f = t.fahrenheit( 0 );
 
  if( f < -100 )      // Sometimes we get spurious readings of -196 degrees
    return;
    
  stats.sample( f );
 
  int out = pid.out( f );
 
  OCR1A = out;                        // Set heater temperature
  
  Serial.print( f, 4 );
  Serial.print( ", " );
  Serial.print( stats.get_mean(), 6 );
  Serial.print( ", " );
  Serial.print( stats.get_sd(), 6 );
 
  Serial.print( ", " );
  leading_spaces( out );
  
  Serial.print( ", " );
  leading_spaces( pid.Integral );
  
  Serial.print( ", " );
  Serial.print( pid.Target );
 
  Serial.println();
 
  delay( 300 );
}