
 /*
 * This program implements the 
 * Sliding Windowed Infinite Fourier Transform (SWIFT)
 * and the slightly lower noise version, alpha-SWIFT
 * 
 * Logan L. Grado, Matthew D. Johnson, and Theoden I. Netoff
 *  Digital Object Identifier 10.1109/MSP.2017.2718039
 *  https://ieeexplore.ieee.org/document/8026592
 * 
 * A bank of digital, complex, resonaters each with the appropriate window
 * copmputes a complete fourier transform at each input sample time.
 * Sample rate in 10 KHz.

 * HARDWARE CONNECTIONS
 *  - GPIO 16 ---> VGA Hsync
 *  - GPIO 17 ---> VGA Vsync
 *  - GPIO 18 ---> 470 ohm resistor ---> VGA Green lo-bit |__ both wired to VGA Green 
 *  - GPIO 19 ---> 330 ohm resistor ---> VGA Green hi_bit |   
 *  - GPIO 20 ---> 330 ohm resistor ---> VGA Blue
 *  - GPIO 21 ---> 330 ohm resistor ---> VGA Red
 *  - RP2040 GND ---> VGA GND
 * 
 *  - GPIO 28 ADC input 
 * in this example, input cam from a laptop through a high-pass biased to Vdd/2
 * 
 * -- SPI for DAC
 * #define PIN_CS   5
 * define PIN_SCK  6
 * define PIN_MOSI 7
 * define SPI_PORT spi0
 *
 * RESOURCES USED
 *  - PIO state machines 0, 1, and 2 on PIO instance 0
 *  - DMA channels 0, 1, 2, and 3
 *  - 153.6 kBytes of RAM (for pixel color data)
 *
 * ======================================
 * Fixed decay time at 10Ksamples/sec
 * VERY roughly: Fs/(decay_samples)*0.5 = FWHM bandwidth in Hz
 * bandwidth is cnstant across frequiencies
 * decay     bandwidth  risetime
 * samples   FWHM HZ    samples
 * 2          ~3400      0  
 * 5          ~1050      1
 * 10         512       2
 * 25         205       5
 * 50         103       10
 * 100        58        10
 * 200        32        10    
 * 300        27        10
 * 400        16        10
 * 
 * ======================================
 * Propportional decay time
 * Bandwidth is proportional to frequency
 * ===
 * decay set to 2.5*Fs/F : rise time 2 samples
 * freq       500 1000 1500 2000 3000
 * bandwidth   98 165  205  400  560
 * ===
 * decay set to 5*Fs/F : rise time 2 samples
 * freq       500 1000 1500 2000 3000
 * bandwidth  58  110  160  212  313
 * ===
 * decay set to 7*Fs/F : rise time 2 samples
 * freq       500 1000 1500 2000 3000
 * bandwidth   39   78  116  154  229
 * ===
 * decay set to 10*Fs/F : rise time 2 samples
 * freq       500 1000 1500 2000 3000
 * bandwidth   32  59  81   110  160
 * ===
 * roughly speaking , the frequency/FWHM = 2*(relative_fall_time_ratio)
 * so roughly frequency/(relative_fall_time_ratio)*0.5 = FWHM in Hz
 * === r
 * roughly speaking, as the freq goes up the product of peak_amp*FWHM stays constant
 * =======================================
 * 
 */

// ==========================================
// === VGA graphics libraryphase shift
// ==========================================
#include "vga16_graphics_v3.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>

// pico specific
#include "pico/stdlib.h"
#include "hardware/pio.h"
#include "hardware/dma.h"
#include "hardware/adc.h"
#include "hardware/spi.h"

//
//#include "hardware/vreg.h"
//#include "hardware/clocks.h"
//
// ==========================================
// === protothreads globals
// ==========================================
#include "hardware/sync.h"
#include "hardware/timer.h"
#include "pico/multicore.h"
#include "string.h"
// protothreads header
#include "pt_cornell_rp2040_v1_4.h"

// unprotected global to halt display
static int run_stop = 1 ;
// restart flag
int restart_graphics = false ;


// ==========================================
// === set up timer ISR  used in this pgm
// ==========================================
// === timer alarm ========================
// !! modifiying alarm zero trashes the cpu 
//        and causes LED  4 long - 4 short
// !! DO NOT USE alarm 0
// This low-level setup is ocnsiderably faster to execute
// than the hogh-level callback

#define ALARM_NUM 1
#define ALARM_IRQ TIMER0_IRQ_1
// ISR interval will be 100 uSec
// 10 KHz
volatile int alarm_period = 100 ;
float Fs = 10000 ;
//
// the actual ISR
void compute_sample(void);
//
static void alarm_irq(void) {
    // mark ISR entry
    gpio_put(15,1);
    // Clear the alarm irq
    hw_clear_bits(&timer_hw->intr, 1u << ALARM_NUM);
    // arm the next interrupt
    // Write the lower 32 bits of the target time to the alarm to arm it
    timer_hw->alarm[ALARM_NUM] = timer_hw->timerawl + alarm_period ;
    //
    // The actual computation
    compute_sample();

    // mark ISR exit
    gpio_put(15,0);
}

// set up the timer alarm ISR
static void alarm_in_us(uint32_t delay_us) {
    // Enable the interrupt for our alarm (the timer outputs 4 alarm irqs)
    hw_set_bits(&timer_hw->inte, 1u << ALARM_NUM);
    // Set irq handler for alarm irq
    irq_set_exclusive_handler(ALARM_IRQ, alarm_irq);
    // Enable the alarm irq
    irq_set_enabled(ALARM_IRQ, true);
    // Enable interrupt in block and at processor
    // Alarm is only 32 bits 
    uint64_t target = timer_hw->timerawl + delay_us;
    // Write the lower 32 bits of the target time to the alarm which
    // will arm it
    timer_hw->alarm[ALARM_NUM] = (uint32_t) target;   
}

// ===========================================
// ADC setup 
// ===========================================

float adc_sample ;
float adc_mean ;
void ADC_setup(void){
  adc_init();
  adc_gpio_init(28);
  // p26 is ADC input 0 - gpio 28 is 2
  adc_select_input(2);
  // free run
  adc_run(1);
  // result is in adc_hw->result
}

// ========================================
// === spi setup 
// =======================================
//SPI configurations
#define PIN_CS   5
#define PIN_SCK  6
#define PIN_MOSI 7
#define SPI_PORT spi0

// constant to tell SPI DAC what to do
// prepend to each 12-bit sample
#define DAC_config_chan_A 0b0011000000000000
// B-channel, 1x, active
#define DAC_config_chan_B 0b1011000000000000
uint16_t DAC_data ; 
void spi_setup(void){
    // Initialize SPI channel (channel, baud rate set to 20MHz)
    // connected to spi DAC
    spi_init(SPI_PORT, 20000000) ;
    // Format (channel, data bits per transfer, polarity, phase, order)
    spi_set_format(SPI_PORT, 16, 0, 0, 0);
    // Map SPI signals to GPIO ports
    //gpio_set_function(PIN_MISO, GPIO_FUNC_SPI);
    gpio_set_function(PIN_SCK, GPIO_FUNC_SPI);
    gpio_set_function(PIN_MOSI, GPIO_FUNC_SPI);
    gpio_set_function(PIN_CS, GPIO_FUNC_SPI) ;
}
//

// ===========================================
// dSP definitions 
// ===========================================
// The Sliding Windowed Infinite Fourier Transform -- with alpha-SWIFT
// Logan L. Grado, Matthew D. Johnson, and Theoden I. Netoff
// Digital Object Identifier 10.1109/MSP.2017.2718039
//file:///home/bruce/Downloads/10078055.pdf

// === SWIFT setup
// can be any number, limited by the speed of the cpu
// set to number desired frequencies in the filter bank
#define N_FREQ  24 //201   

// A list of the frequencys to transform
// need not be uniform
float F[N_FREQ] ;

// the complex spectral estimate
// frequencies are normailsed and run 0 to pi
// at 10khz sample rate that is 0 to 5kz
float Xr[N_FREQ], Xi[N_FREQ] ;
float gain[N_FREQ] ;
float dk_rate[N_FREQ] ;
float rise_rate[N_FREQ] ;

// additive components of the alpha-SWIFT
float slow_Xr[N_FREQ], slow_Xi[N_FREQ] ;
float fast_Xr[N_FREQ], fast_Xi[N_FREQ] ;

// and the mag of X
float Fmag[N_FREQ], lp_Fmag[N_FREQ];

// the per-sample complex phase shift at each frequecy
// omega = 2*pi*F/Fs so OHM = cos(omega) + j*sin(omega)
float slow_OHMr[N_FREQ], slow_OHMi[N_FREQ] ;
float fast_OHMr[N_FREQ], fast_OHMi[N_FREQ] ;

// plot parameters
float freq_spacing = 10.0f ;
// current frequency
float current_freq ;
short current_x ;
short draw_start_flag = true;
// and the gain for the plot
float  plot_gain = 0.001f ;
// color list
short color[N_FREQ], last_plot_point[N_FREQ] ;

short avail_colors[8] = { 3, 5, 7, 8, 10, 11, 12, 15} ;

//  dds to output filter bank content and graphics
#define twoPi 6.283185f
unsigned int dds_accum[N_FREQ], accum1;
unsigned int dds_inc[N_FREQ], inc1;
float dds_phase[N_FREQ] ;
float dds_phase_inc[N_FREQ] ;
float dds_last_phase[N_FREQ] ;
float sine_table[1024] ;
float current_freq ;
int phase_fix_counter, phase_fix_interval ;

// audio input
float input_sample ;

// dac ouptut signal
float test_out ;

//====================================
// === magnitude approx good to about +/-2%
// see https://en.wikipedia.org/wiki/Alpha_max_plus_beta_min_algorithm
#define min(a,b) (((a)<(b))?(a):(b))
#define max(a,b) (((a)>(b))?(a):(b))
//====================================
float magnitude(float real, float img){
  float mmax, mmin ;
  float c1 = 0.89820f ;
  float c2 = 0.48597f ;
  
  mmin = min(fabsf(real), fabsf(img)) ; 
  mmax = max(fabsf(real), fabsf(img)) ;
  // 
  return max(mmax, ((mmax*c1) + (mmin*c2) )); 
}


// ==================================================
// === graphics demo -- RUNNING on core 0
// ==================================================
static PT_THREAD (protothread_graphics(struct pt *pt)) {
    PT_BEGIN(pt);

    // Draw some filled rectangles
    fillRect(20, 1, 400, 50, LIGHT_BLUE); // blue box
    drawTextVGA437(100, 20, "FFT spectrogram and SWIFT", WHITE, LIGHT_BLUE) ;

    fillRect(435, 1, 200, 50, YELLOW); // green box
    // Write some text
     // Draw some filled rectangles
    fillRect(0, 0, 640, 480, BLACK); // erase
    fillRect(64, 0, 176, 50, BLUE); // blue box
    fillRect(250, 0, 176, 50, YELLOW); // red box
    fillRect(435, 0, 176, 50, GREEN); // green box

    // Write some text
    drawTextTiny8(65, 0, "Raspberry Pi Pico 2", WHITE, BLUE) ;
    drawTextGLCD(65, 10, "SWIFT spectrum", WHITE, BLUE) ;
    drawTextGLCD(65, 20, "Bruce Land brl4@cornell.edu", WHITE, BLUE) ;
    drawTextGLCD(65, 30, "Hunter Adams vha3@cornell.edu", WHITE, BLUE) ;
    //
    drawTextGLCD(438, 10, "Protothreads rp2040/2350 1.4", BLACK, GREEN) ;
    drawTextGLCD(438, 20, "One CPU; clock 150 Mhz", BLACK, GREEN) ;
    drawTextGLCD(438, 30, "vga16_v3 driver", BLACK, GREEN) ;
    drawTextGLCD(438, 40, "Double Buffer 60Hz", BLACK, GREEN) ;
    //
    drawTextVGA437(260, 5, "SWIFT Transform", BLACK, YELLOW) ;
     drawTextGLCD(260, 20, "Filter Bank", BLACK, YELLOW) ;
    //

    //copy_buffer_to_other() ;

    // =====================================
    // =====================================
    // This is where you set up the desired frequency
    // analysis that you want
    // You must set up each seqarate resonator
    // --set up F, gain, decay and rise time  arrays
    
    for(int i=0; i<N_FREQ; i++){
      // choose to make a log-spaced spectrum
      //F[i] = 1000.0f * log10((float) (10*i+2)) ;
      if (i==0) F[0] = 100.0f ;
      else  F[i] = F[i-1] * 1.175 ; // for 12 bands:1.41f ; // 
      
      // set the gain of each frequency
      // this seems to equalize fairly well
      gain[i] = F[i] / 100.0f;

      // set the decay time constant for each freq
      // sert to five times the period of F[i]
      #define dk_ratio 3.0f
      dk_rate[i] = exp(-1.0f/(dk_ratio * Fs /F[i])) ; // samples

      // set the rise rate consrtant for each freq
      // bigger lowers noise, but attenuates higher freq
      // 0.1 is very small, 10.0 is big
      #define rise_time 1.0f
      rise_rate[i] = exp(-1.0f/rise_time) ; // samples
      // =====================================
      // =====================================
      
      // compute complex frequenc for phasor
      // omega = 2*pi*F/Fs so OHM = cos(omega) + j*sin(omega)
      // and premultiply by rate constants
      float omega = twoPi * F[i]/Fs ;
      // we alos want the phase/sample in dds units
     //>> dds_phase_inc[i] = omega * 256.0f / twoPi;
      // omega is in ratians/timestep 
      // for resunthesis want the units to be (fraction of 256)/timestep
      //dds_phase_inc[i]= (int)(omega * 256.0f /  6.28318f) ;

      slow_OHMr[i] = dk_rate[i] * cos(omega) ;
      slow_OHMi[i] = dk_rate[i] * sin(omega); 
      fast_OHMr[i] = rise_rate[i] * cos(omega) ;
      fast_OHMi[i] = rise_rate[i] * sin(omega);      
    }

    // ======================================
    // set up output DDS frequencies
     for(int i=0; i<N_FREQ; i++){
        dds_inc[i] = F[i] * pow(2,32)/ Fs ;
     }

     inc1 = 500.0 * pow(2,32)/ Fs ;

     // ======================================

    ////
    while(true) {

      // restaert logic 
        PT_YIELD_UNTIL(pt, restart_graphics || draw_start_flag); 
      if (restart_graphics){
            // reset the flag
            restart_graphics = false ;    
            draw_start_flag = true ;   
            // restarts the current thread at PT_BEGIN(pt);
            // as if it has not executed before
            PT_RESTART(pt) ;
        }

        clearLowFrame(60, BLACK) ;

      // wait rfor start
      //PT_YIELD_UNTIL(pt, draw_start_flag); 
      draw_start_flag = false ;
      current_x = 50 ;
      current_freq = 50.0f ;

      // the main loop sets a frequency, 
      // waits for steady state response
      // plots the magnitude of each resonator at that freq
      // increments the frequency ubtil reach max scan freq

      while (current_freq < 5000.0f){
        //>>inc1 = current_freq * pow(2,32)/ Fs ;
        // now yield to let filters reach steady state
        // at the new frequency
        PT_YIELD_usec(1000 * (int)(Fs/current_freq)) ;

        //  plot the power spectra
        for(int i=0; i<N_FREQ; i++){
          // the scale factors will change depending on recording conditions
          //Fmag[i] = magnitude(Xr[i], Xi[i]) * plot_gain ;
          Fmag[i] = lp_Fmag[i] * gain[i] * plot_gain ;
          short plot_point = -(short)Fmag[i]+270 ;
          fillRect(current_x, plot_point, 2, 2, avail_colors[i & 0x07]) ;  
          // sconnect the dots, execpt for the very first point
          if(current_x > 50){
             drawLine(current_x-1, last_plot_point[i], current_x, plot_point, avail_colors[i & 0x07] ) ;
          }
          // magnify low freq spectrum horizontally
          if ((current_freq <1000) && (current_x > 50)){
            fillRect((current_x-50)*5+50, plot_point+200, 2, 2, avail_colors[i & 0x07]) ; 
            drawLine((current_x-50)*5+50-5, last_plot_point[i]+200, (current_x-50)*5+50, plot_point+200, avail_colors[i & 0x07] ) ;
          }
          // need last point to plot connecting lines
          last_plot_point[i] = plot_point ;        }
        //

        current_freq += freq_spacing ;
        current_x += 1 ;

      }
      #define vert_x 275
      drawTextTiny8(5000/(int)freq_spacing+40, vert_x, "|", RED, BLACK);
      drawTextAscii(5000/(int)freq_spacing+46, vert_x, "5000", WHITE, BLACK);
      drawTextTiny8(1000/(int)freq_spacing+40, vert_x, "|", RED, BLACK);
      drawTextAscii(1000/(int)freq_spacing+46, vert_x, "1000", WHITE, BLACK);
      drawTextTiny8(2000/(int)freq_spacing+40, vert_x, "|", RED, BLACK);
      drawTextAscii(2000/(int)freq_spacing+46, vert_x, "2000", WHITE, BLACK);
      drawTextTiny8(3000/(int)freq_spacing+40, vert_x, "|", RED, BLACK);
      drawTextAscii(3000/(int)freq_spacing+46, vert_x, "3000", WHITE, BLACK);
      drawTextTiny8(4000/(int)freq_spacing+40, vert_x, "|", RED, BLACK);
      drawTextAscii(4000/(int)freq_spacing+46, vert_x, "4000", WHITE, BLACK);
      //
      drawTextAscii(1000/(int)freq_spacing*5+20, 470, "1000", WHITE, BLACK);
      // test freame timing
      // gpio_put(14,1);

       char str[80] ;
      // // plot gain
      sprintf(str, "Plot gain: %5.5f  ", plot_gain);
      drawTextTiny8(150, 60, str, WHITE, BLACK);
      // //
   }
   PT_END(pt);
} // graphics thread


// ==================================================
// === user's serial input thread on core 1
// ==================================================
// serial_read an serial_write do not block any thread
// except this one
static PT_THREAD (protothread_serial(struct pt *pt))
{
    PT_BEGIN(pt);
      
      static char *cmd, *arg1, *arg2, *arg3, *arg4;

      while(1) {
        // prompt the user for input
        sprintf(pt_serial_out_buffer, "cmd> ") ;
        // spawn a thread to do the non-blocking serial write
        // to print the pt_serial_out_buffer string
        // this output routine YIELDs between each character printed.
        serial_write ; 
        // wait for user
        serial_read ;
        // tokenize the input string
        //  --- NOTE that strtok is not re-entrant!
        //  ---   Use the reentrant, thread-safe, version if you need to use it
        //   ---  in moe than one thread!
        cmd = strtok(pt_serial_in_buffer, "  ");
        arg1 = strtok(NULL, "  ");
        arg2 = strtok(NULL, "  ");
        arg3 = strtok(NULL, "  ");
        arg4 = strtok(NULL, "  ");

        //wait
        sleep_ms(10);
        // commmands:
        // restart
        // gain value -- sets the vertical size of spectrum
        // freq value -- sets the frequency bin spacing in Hz
        // 
        // dk r value -- sets the longer time constant (decay time) to a multiple of the dds period
        // dk f value -- sets the longer time constant (decay time) to a fixed number of samples
        // 

        // restart graphics thread
        if(strcmp(cmd,"restart")==0){
          restart_graphics = true ;  
        }
        //
        // plot gain
        else if(strcmp(cmd,"gain")==0){
          sscanf(arg1,"%f", &plot_gain);  
          draw_start_flag = true ;  
        }
        
        // freq delta
        else if(strcmp(cmd,"freq")==0){
          sscanf(arg1,"%f", &freq_spacing);  
          restart_graphics = true ;  
        }
        //
        //
        else{
          printf("Huh?\n\r") ;
        }
        //
        // NEVER exit while
      } // END WHILE(1)
  PT_END(pt);
} // serial thread

// ==================================================
// === dsp ISR -- 
// ==================================================
// 
void compute_sample(void){
    // DDS
    // make the sine wave to send to filter
    accum1 += inc1 ;
    // update main waveform
    //sine_table[accum1>>24] ;
    input_sample = sine_table[accum1>>22] ;

    // convert ADC to float
    adc_sample = (float) (adc_hw->result) ;
    
    // low pass to get mean
    // then subtract the mean (highpass)
    // high pass at around 0.1 second, so one Hz
    adc_mean = 0.999f * adc_mean + 0.001 * adc_sample ;
    adc_sample -= adc_mean ;

    input_sample = adc_sample ;

    // 
   // zero the output sum
    test_out = 0.0f ;

    for(int i=0; i<N_FREQ; i++){   //N_FREQ
      float temp ;
      // the slow arm of the swift
      temp = (slow_Xr[i] * slow_OHMr[i] - slow_Xi[i] * slow_OHMi[i]) + input_sample ; //-
      slow_Xi[i] = (slow_Xr[i] * slow_OHMi[i] + slow_Xi[i] * slow_OHMr[i]) ;//-
      slow_Xr[i] = temp ;
      //
      // the fast arm of the swift
      temp = (fast_Xr[i] * fast_OHMr[i] - fast_Xi[i] * fast_OHMi[i]) + input_sample ; //-
      fast_Xi[i] = (fast_Xr[i] * fast_OHMi[i] + fast_Xi[i] * fast_OHMr[i]) ;//-
      fast_Xr[i] = temp ;
      //
      // combine the two to make X
      Xr[i] = slow_Xr[i] - fast_Xr[i] ;
      Xi[i] = slow_Xi[i] - fast_Xi[i] ;
      //
      // low pass over about 50 samplwes (10 mSec)
      float temp_mag =  magnitude(Xr[i], Xi[i]) ;
      lp_Fmag[i] = 0.95f * lp_Fmag[i] + 0.05f * temp_mag;

      // now test output -- simple sum output for testing
      /// +++ this works ====
       //test_out += Xr[i] ; // 
      // ====================

      // actuial phase of the resonator
      float temp_phase = fmodf( (atan2f(Xi[i], Xr[i]) + twoPi), twoPi) * 1023.0f / twoPi ;
      // actual phase diff of resonator
      dds_phase_inc[i] =  (temp_phase >= dds_last_phase[i])? temp_phase - dds_last_phase[i] :
                                                            (1023.0f - dds_last_phase[i]) + temp_phase ;
      // copy for next pass
      dds_last_phase[i] = temp_phase ;
      // now we can get dds phase 
      #define shift_ratio 0.85f
      dds_phase[i]= fmodf(shift_ratio * dds_phase_inc[i] + dds_phase[i], 1023.0f) ;
      // sum up sine waves
      test_out += sine_table[(int)dds_phase[i]] *  temp_mag ;
      
    } // end the per-filter loop

    //   up date the phase adjustor
     //phase_fix_counter = ++phase_fix_counter % N_FREQ ;
    // phase_fix_interval++ ;

    // scale may need to change depending on source amp
    // offset converts signed to offset binary
    DAC_data = DAC_config_chan_A | (((int)(test_out*0.0001f) + 2048) & 0xfff)   ; 
    spi0_hw->dr = DAC_data ;
    DAC_data = DAC_config_chan_B | (((int)(test_out*0.0001f) + 2048) & 0xfff)   ; 
    // NOTE === nonblocking SPI write ===
    spi0_hw->dr = DAC_data ;
    //blocking out: spi_write16_blocking(SPI_PORT, &DAC_data, 1) ;
    
}

// ========================================
// === core 1 main -- started in main below
// ========================================
void core1_main(){ 
  //
  //  === add threads  ====================
  // for core 1
  pt_add_thread(protothread_serial) ;
  pt_add_thread(protothread_graphics);
  //
  // === initalize the scheduler ==========
  pt_schedule_start ;
  // NEVER exits
  // ======================================
}

// ========================================
// === core 0 main
// ========================================
int main(){
  // set the clock
  //set_sys_clock_khz(250000, true); // 171us
  // start the serial i/o
  stdio_init_all() ;
  // announce the threader version on system reset
  printf("\n\rProtothreads RP2040/2350 v1.4 two-core\n\r");

  // Initialize the VGA screen
  initVGA() ;

  // the spi for DAC
  spi_setup() ;

  // === congure the ADC =========
  ADC_setup() ;

  //
  // fire off interrupt
  alarm_in_us(alarm_period);

  gpio_init(15) ;	
  gpio_set_dir(15, GPIO_OUT) ;
  gpio_init(14) ;	
  gpio_set_dir(14, GPIO_OUT) ;

  // start core 1 threads
  multicore_reset_core1();
  multicore_launch_core1(&core1_main);

  // dds table for testing
  for (int i=0; i<1024; i++) {
    // sine table is in natural +1/-1 range
    sine_table[i] = 1000.0f*(sin(2*3.14159*(float)i/1024)) ;
  }

  // === config threads ========================
  // for core 0
  //
  //pt_add_thread(protothread_graphics);
  //pt_add_thread(protothread_toggle25);
  //
  // === initalize the scheduler ===============
  pt_schedule_start ;
  // NEVER exits
  // ===========================================
} // end main