//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Mon Jun 24 20:05:15 PDT 2002
// Last Modified: Wed Jul 10 13:50:11 PDT 2002
// Filename:      ...sig/examples/all/chordplot.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/humdrum/chordplot.cpp
// Syntax:        C++; museinfo
//
// Description:   Determine the root of a pitch set according to
//                different root-interval class sizes.
// 

#define MYPI 3.141592653589793238462643383279

#include "humdrum.h"
#include <string.h>
#include <math.h>

// function declarations
void   checkOptions             (Options& opts, int argc, char* argv[]);
void   example                  (void);
void   usage                    (const char* command);
void   printScores              (Array<double>& rootscores);
int    analyzeChromatic         (Array<double>& rootscores, HumdrumFile& infile,
                                   double theta1, double theta2, int startline, 
                                   int stopline, int debugQ = 0);
double base40ToBase7            (int pitch);
char   base7pcToName            (int value);
void   identifyErrorsChromatic  (HumdrumFile& infile, int type, int debugQ);
void   getDistancesChromatic1   (Array<double>& distances, double theta);
void   getDistancesChromatic2   (Array<double>& distances, double theta1,
                                   double theta2);

// global variables
Options      options;            // database for command-line arguments
int          debugQ     = 0;     // used with --debug option
int          appendQ    = 0;     // used with -a option
int          rootQ      = 0;     // used with -R option
double       theta1     = 60.0;  // used with -t option
double       theta2     = 60.0;  // used with -u option
int          invertQ    = 0;     // used with -i option
int          distanceQ  = 0;     // used with -d option
int          notesQ     = 0;     // used with -n option

///////////////////////////////////////////////////////////////////////////

int main(int argc, char* argv[]) {
   HumdrumFile infile;

   // process the command-line options
   checkOptions(options, argc, argv);

   infile.clear();

   // if no command-line arguments read data file from standard input
   if (options.getArgCount() < 1) {
      infile.read(cin);
   } else {
      infile.read(options.getArg(1));
   }

   if (errorsQ) {
      if (options.getBoolean("theta2")) {
         identifyErrorsChromatic(infile, 2, debugQ);
      } else {
         identifyErrorsChromatic(infile, 1, debugQ);
      }
   } else {

      //if (options.getBoolean("theta2")) {
      // analyzeChromatic(rootscores, infile, theta1, theta2, 0, 0, -1, debugQ);
      // printScores(rootscores);
      //} else {
      // analyzeChromatic(rootscores, infile, theta1, theta2, 0, 0, -1, debugQ);
      // printScores(rootscores);
      //}

   } 

   return 0;
}


///////////////////////////////////////////////////////////////////////////


//////////////////////////////
//
// findErrorsX -- compare root scores value to root in data.
//

void findErrorsX(HumdrumFile& infile, int debugQ) {
   int oldline = 0;
   int init = 0;
   int best;
   int correct = 0;
   int i;
   int j;
   Array<double> rootscores(7);
   rootscores.setAll(0);
   for (i=0; i<infile.getNumLines(); i++) {
      if (!infile[i].isData()) {
         continue;
      }
      for (j=0; j<infile[i].getFieldCount(); j++) {
         if (strcmp(infile[i].getExInterp(j), "**root") == 0) {
            if (strcmp(infile[i][j], ".") == 0) {
               break;
            } 
            if (strchr(infile[i][j], '_') != NULL) {  // cont. ties
               break;
            }
            if (strchr(infile[i][j], ']') != NULL) {  // tie ends
               break;
            }
            if (init == 0) {
               init = 1;
               oldline = i;
               /////////////////////////
               // for diatonic data:
               // correct = (int)base40ToBase7(
               //       Convert::kernToBase40(infile[i][j])) % 7;
               // correct = (correct + 70) % 7;
               /////////////////////////
               correct = Convert::kernToBase40(infile[i][j]) % 40;
               break;
            }
            if (options.getBoolean("theta2")) {
               best=analyzeChromatic(rootscores, infile, theta1, 
                     theta2, oldline, i-1, correct, debugQ);
            } else {
               best=analyzeChromatic(rootscores, infile, theta1, theta2, 
                     oldline, i-1, debugQ);
            }
            // best = ((int)base40ToBase7(best) + 70) % 7;
            best = (best + 40) % 40;
            if (best != correct) {
               cout << "Error on line " << oldline + 1
                    << "\tcorrect=" << base7pcToName(correct)
                    << "\tvalue=" << base7pcToName(best) 
                    << endl;
            }
            oldline = i;
            // correct = (int)base40ToBase7(Convert::kernToBase40(infile[i][j]))%7;
            correct = Convert::kernToBase40(infile[i][j]) % 40;
            // correct = (correct + 70) % 7;
            break;
         }
      }
      
   }
   best=analyzeChromatic(rootscores,infile,theta,oldline,i-1,debugQ);
   if (best != correct) {
      char zbuffer[128] = {0};
      cout << "Error on line " << oldline + 1
           << "\tcorrect=" << Convert::base40ToKern(zbuffer, correct+4*40) 
           << endl;
      cout << "\tvalue=" << Convert::base40ToKern(zbuffer, best+4*40) << endl;

//         << "\tcorrect=" << base7pcToName(correct)
//         << "\tvalue=" << base7pcToName(best) << endl;
   }

}

//////////////////////////////
//
// identifyErrorsChromatic -- compare root scores value to root in
//     humdrum file and print out errors.
//

void identifyErrorsChromatic(HumdrumFile& infile, int type, int debugQ) {
   int oldline = 0;
   int init = 0;
   int best;          // best root according to algorithm
   int correct = 0;   // correct base40 root
   int i;
   int j;

   // get interval distances:
   Array<double> distances;
   switch (type) {
      case 1:
         getDistancesChromatic1(distances, theta1);
         break;
      case 2:
         getDistancesChromatic2(distances, theta1, theta2);
         break;
      default:
         cout << "Error: unknown type: " << type << endl;
   }

 
   // go though the file and compare all chords to **root info:
   Array<double> rootscores(40);
   rootscores.setAll(0);
   for (i=0; i<infile.getNumLines(); i++) {
      if (!infile[i].isData()) {
         continue;
      }
      for (j=0; j<infile[i].getFieldCount(); j++) {
         if (strcmp(infile[i].getExInterp(j), "**root") == 0) {
            if (strcmp(infile[i][j], ".") == 0) {
               break;
            } 
            if (strchr(infile[i][j], '_') != NULL) {  // cont. ties
               break;
            }
            if (strchr(infile[i][j], ']') != NULL) {  // tie ends
               break;
            }
            if (init == 0) {
               init = 1;
               oldline = i;
               correct = Convert::kernToBase40(infile[i][j]) % 40;
               break;
            }
            best=analyzeChromatic(rootscores, infile, theta1, 
                  theta2, type, oldline, i-1, correct, debugQ);
            if (best != correct) {
               char zbuffer[128] = {0};
               cout << "Error on line " << oldline + 1
                    << "\tcorrect=" 
                    << Convert::base40ToKern(zbuffer, correct+4*40) 
                    << endl;
               cout << "\tvalue=" 
                    << Convert::base40ToKern(zbuffer, best+4*40) << endl;
            }
            oldline = i;
            correct = Convert::kernToBase40(infile[i][j]) % 40;
            break;
         }
      }
      
   }
   best=analyzeChromatic(rootscores,infile,theta,oldline,i-1,debugQ);
   if (best != correct) {
      char zbuffer[128] = {0};
      cout << "Error on line " << oldline + 1
           << "\tcorrect=" << Convert::base40ToKern(zbuffer, correct+4*40) 
           << endl;
      cout << "\tvalue=" << Convert::base40ToKern(zbuffer, best+4*40) << endl;

   }

}



//////////////////////////////
//
// analyzeRootProbability2 --
//

int analyzeChromatic(Array<double>& rootscores, HumdrumFile& infile, 
      double theta1, double theta2, int startline, int stopline, int debugQ) {

   // get interval distances
   Array<double> distances;
   if (chromatic2Q) {
      getDistancesChromatic1(distances, theta1, theta2);
   } else {
      getDistancesChromatic1(distances, theta1);
   }

   // extract note data
   Array<double> absbeat;
   Array<int>    pitches;
   Array<double> durations;
   Array<double> levels;
   infile.getNoteArray(absbeat, pitches, durations, levels, startline,
      stopline);

   char buffer[64] = {0};
   if (debugQ) {
      cout << "====================" << endl;
      for (i=startline; i<=stopline; i++) {
         cout << infile[i] << "\n";
      }
      cout << "--------------------" << endl;
      cout << "Start line=" << startline+1 
           << "\tEnd line=" << stopline+1 << endl;
      for (i=0; i<pitches.getSize(); i++) {
         cout << "pitch = " 
              << Convert::base40ToKern(buffer, pitches[i]) << endl;
      }
   }

   rootscores.setSize(40);
   rootscores.setAll(0.0);
   int j;
   for (i=0; i<rootscores.getSize(); i++) {
      for (j=0; j<pitches.getSize(); j++) {
         rootscores[i] += distance[(pitches[j]-i+400)%40];
      }
   }

   int min = 0;
   for (i=1; i<rootscores.getSize(); i++) {
      if (rootscores[i] < rootscores[min]) {
         min = i;
      }
   }

   if (debugQ) {
      cout << "BEST ROOT = " << min << endl;
   }

   return min;
}



//////////////////////////////
//
// printScores --
//

void printScores(Array<double>& rootscores) {
   int i;
   int findex;
   char buffer[64] = {0};
   for (i=0; i<rootscores.getSize(); i++) {
      findex = (17 + i * 23) % 40;
      Convert::base40ToKern(buffer, findex+4*40);
      if (buffer[0] == '\0') {
         continue;
      }
      cout << buffer << "\t" << i << "\t" << findex
           << "\t";
      if (invertQ) {
         cout << 1.0/rootscores[(findex-2+40)%40] << endl;
      } else {
         cout << rootscores[(findex-2+40)%40] << endl;
      }
   }
}



//////////////////////////////
//
// checkOptions -- validate and process command-line options.
//

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("a|append=b",     "append analysis to data in output");   
   opts.define("i|invert=b",     "invert scores");   
   opts.define("d|distance=b",   "display interval distances");   
   opts.define("t|theta|theta1=d:60.0", "chromatic plane angle");   
   opts.define("u|theta2=d:60.0", "chromatic plane angle 2");   
   opts.define("e|errors=b",     "identify root errors in the input file");   
   opts.define("n|notes=b",      "display note set for each chord");   

   opts.define("debug=b",        "trace input parsing");   
   opts.define("author=b",       "author of the program");   
   opts.define("version=b",      "compilation information"); 
   opts.define("example=b",      "example usage"); 
   opts.define("h|help=b",       "short description"); 
   opts.process(argc, argv);
   
   // handle basic options:
   if (opts.getBoolean("author")) {
      cout << "Written by Craig Stuart Sapp, "
           << "craig@ccrma.stanford.edu, June 2002" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 18 June 2002" << endl;
      cout << "compiled: " << __DATE__ << endl;
      cout << MUSEINFO_VERSION << endl;
      exit(0);
   } else if (opts.getBoolean("help")) {
      usage(opts.getCommand());
      exit(0);
   } else if (opts.getBoolean("example")) {
      example();
      exit(0);
   }

   debugQ    =  opts.getBoolean("debug");
   errorsQ   =  opts.getBoolean("errors");
   distanceQ =  opts.getBoolean("distance");
   notesQ    =  opts.getBoolean("notes");
   appendQ   =  opts.getBoolean("append");
   theta1    =  opts.getDouble("theta");
   theta2    =  opts.getDouble("theta2");
   invertQ   =  opts.getBoolean("invert");
}



//////////////////////////////
//
// example -- example usage of the maxent program
//

void example(void) {
   cout <<
   "                                                                        \n"
   << endl;
}



//////////////////////////////
//
// usage -- gives the usage statement for the quality program
//

void usage(const char* command) {
   cout <<
   "                                                                        \n"
   << endl;
}



//////////////////////////////
//
// base40ToBase7 -- Convert pitch to diatonic pitch
//   .0 = natural
//   .1 = sharp
//   .2 = double sharp
//   .9 = flat
//   .8 = double flat
//

double base40ToBase7(int pitch) {
   int octave = pitch / 40;
   if (octave < 0) {
      // rest
      return -1;
   }

   switch (pitch % 40) {
      case  0:  return octave * 7 + 0.9;     // C--
      case  1:  return octave * 7 + 0.8;     // C-
      case  2:  return octave * 7 + 0.0;     // C
      case  3:  return octave * 7 + 0.1;     // C#
      case  4:  return octave * 7 + 0.2;     // C##
      case  5:  return -1;                   // X
      case  6:  return octave * 7 + 1.9;     // D--
      case  7:  return octave * 7 + 1.8;     // D-
      case  8:  return octave * 7 + 1.0;     // D
      case  9:  return octave * 7 + 1.1;     // D#
      case 10:  return octave * 7 + 1.2;     // D##
      case 11:  return -1;                   // X
      case 12:  return octave * 7 + 2.9;     // E--
      case 13:  return octave * 7 + 2.8;     // E-
      case 14:  return octave * 7 + 2.0;     // E
      case 15:  return octave * 7 + 2.1;     // E#
      case 16:  return octave * 7 + 2.2;     // E##
      case 17:  return octave * 7 + 3.9;     // F--
      case 18:  return octave * 7 + 3.8;     // F-
      case 19:  return octave * 7 + 3.0;     // F
      case 20:  return octave * 7 + 3.1;     // F#
      case 21:  return octave * 7 + 3.2;     // F##
      case 22:  return -1;                   // X
      case 23:  return octave * 7 + 4.9;     // G--
      case 24:  return octave * 7 + 4.8;     // G-
      case 25:  return octave * 7 + 4.0;     // G
      case 26:  return octave * 7 + 4.1;     // G#
      case 27:  return octave * 7 + 4.2;     // G##
      case 28:  return -1;                   // X
      case 29:  return octave * 7 + 5.9;     // A--
      case 30:  return octave * 7 + 5.8;     // A-
      case 31:  return octave * 7 + 5.0;     // A
      case 32:  return octave * 7 + 5.1;     // A#
      case 33:  return octave * 7 + 5.2;     // A##
      case 34:  return -1;                   // X
      case 35:  return octave * 7 + 6.9;     // B--
      case 36:  return octave * 7 + 6.8;     // B-
      case 37:  return octave * 7 + 6.0;     // B
      case 38:  return octave * 7 + 6.1;     // B#
      case 39:  return octave * 7 + 6.2;     // B##
   }

   return -1;

}



//////////////////////////////
//
// base7pcToName -- base 7 pitch class name: C D E F G A B
//

char base7pcToName(int value) {
   value = (value + 70) % 7;
   switch (value) {
      case 0:  return 'C';
      case 1:  return 'D';
      case 2:  return 'E';
      case 3:  return 'F';
      case 4:  return 'G';
      case 5:  return 'A';
      case 6:  return 'B';
      default: return 'X';
   }
   return 'X';
}


//////////////////////////////
//
// analyzeRootProbability3 -- Chromatic intervals with two parameters:
//     theta1 and theta2.  
//

int analyzeRootProbability3(Array<double>& rootscores, HumdrumFile& infile, 
      double t1, double t2, int startline, int stopline, int root, int debugQ) {

// p1 = perfect 5th distance
static int p1[40] = {
     6, /* C-- */  5, /* C-  */  0, /* C   */  4, /* C#  */  5, /* C## */
 10000, /* X   */
     4, /* D-- */  3, /* D-  */  2, /* D   */  3, /* D#  */  4, /* D## */
 10000, /* X   */
     6, /* E-- */  1, /* E-  */  1, /* E   */  5, /* E#  */  6, /* E## */
     5, /* F-- */  4, /* F-  */  3, /* F   */  3, /* F#  */  4, /* F## */
 10000, /* X   */
     7, /* G-- */  2, /* G-  */  1, /* G   */  2, /* G#  */  6, /* G## */
 10000, /* X   */
     5, /* A-- */  4, /* A-  */  3, /* A   */  4, /* A#  */  5, /* A## */
 10000, /* X   */
     3, /* B-- */  2, /* B-  */  2, /* B   */  3, /* B#  */  7  /* B## */
};
   

// p2 = major or minor third distance; negative are for majors
static int p2[40] = {
    -5, /* C-- */ -3, /* C-  */  0, /* C   */  1, /* C#  */  3, /* C## */
 10000, /* X   */
    -4, /* D-- */ -2, /* D-  */  0, /* D   */  2, /* D#  */  4, /* D## */
 10000, /* X   */
    -4, /* E-- */ -1, /* E-  */  1, /* E   */  2, /* E#  */  4, /* E## */
    -5, /* F-- */ -3, /* F-  */ -1, /* F   */  1, /* F#  */  3, /* F## */
 10000, /* X   */
    -5, /* G-- */ -2, /* G-  */  0, /* G   */  2, /* G#  */  3, /* G## */
 10000, /* X   */
    -4, /* A-- */ -2, /* A-  */  0, /* A   */  2, /* A#  */  4, /* A## */
 10000, /* X   */
    -3, /* B-- */ -1, /* B-  */  1, /* B   */  3, /* B#  */  4  /* B## */
};

// p3 = theta 1 or theta 2
static double p3[40] = {
     1, /* C-- */  1, /* C-  */  0, /* C   */  2, /* C#  */  2, /* C## */
 10000, /* X   */
     1, /* D-- */  1, /* D-  */  0, /* D   */  2, /* D#  */  2, /* D## */
 10000, /* X   */
     1, /* E-- */  1, /* E-  */  2, /* E   */  2, /* E#  */  2, /* E## */
     1, /* F-- */  1, /* F-  */  1, /* F   */  2, /* F#  */  2, /* F## */
 10000, /* X   */
     1, /* G-- */  1, /* G-  */  0, /* G   */  2, /* G#  */  2, /* G## */
 10000, /* X   */
     1, /* A-- */  1, /* A-  */  0, /* A   */  2, /* A#  */  2, /* A## */
 10000, /* X   */
     1, /* B-- */  1, /* B-  */  2, /* B   */  2, /* B#  */  2  /* B## */
};


   // convert thetas to radians:
   theta1 = t1 * MYPI / 180.0;
   theta2 = t2 * MYPI / 180.0;

   int i;
   double theta;
   double avalue;
   double bvalue = 0.0;
   double distance[40];
   for (i=0; i<40; i++) {
      avalue = p1[i] * sin(MYPI - theta1 - theta2)/sin(theta2);
      if (p2[i] < 0) {
         bvalue = -p2[i];
      } else {
         bvalue = p2[i] * sin(theta1)/sin(theta2);
      }
      if (p3[i] == 1) {
         theta = theta1;
      } else if (p3[i] == 2) {
         theta = theta2;
      } else if (p3[i] == 0) {
         theta = 0;
      } else {
         theta = MYPI;
      }
      distance[i] = sqrt(avalue * avalue + bvalue * bvalue -
            2.0 * avalue * bvalue * cos(theta));
   }
   if (distanceQ) {
      char buffer[64] = {0};
      cout << "# theta1 = " << theta1 * 180.0 / MYPI  << "\n";
      cout << "# theta2 = " << theta2 * 180.0 / MYPI  << "\n";
      for (i=0; i<40; i++) {
         cout << Convert::base40ToKern(buffer, i+4*40);
         cout << "\t" << distance[i] << "\n";
      }
      exit(0);
   }

   // extract note data
   Array<double> absbeat;
   Array<int>    pitches;
   Array<double> durations;
   Array<double> levels;
   infile.getNoteArray(absbeat, pitches, durations, levels, startline,
      stopline);

   char buffer[64] = {0};
   if (debugQ) {
      cout << "====================" << endl;
      for (i=startline; i<=stopline; i++) {
         cout << infile[i] << "\n";
      }
      cout << "--------------------" << endl;
      cout << "Start line=" << startline+1 
           << "\tEnd line=" << stopline+1 << endl;
      for (i=0; i<pitches.getSize(); i++) {
         cout << "pitch = " 
              << Convert::base40ToKern(buffer, pitches[i]) << endl;
      }
   }

   if (notesQ) {
      for (i=0; i<pitches.getSize(); i++) {
         cout << Convert::base40ToKern(buffer, pitches[i]);
         if (i<pitches.getSize()-1) {
            cout << " ";
         }
      }
      cout << "\n";
   }

   rootscores.setSize(40);
   rootscores.setAll(0.0);
   int j;
   for (i=0; i<rootscores.getSize(); i++) {
      for (j=0; j<pitches.getSize(); j++) {
         rootscores[i] += distance[(pitches[j]-i+400)%40];
      }
   }

   int min = 0;
   for (i=1; i<rootscores.getSize(); i++) {
      if (rootscores[i] < rootscores[min]) {
         min = i;
      }
   }

   if (debugQ) {
      cout << "BEST ROOT = " << min << endl;
   }

   return min;
}




//////////////////////////////
//
// getDistancesChromatic1 -- 
//

void getDistancesChromatic1(Array<double>& distances, theta) {
   static int xpos[40] = {
        7, /* C-- */  7, /* C-  */  0, /* C   */  7, /* C#  */  7, /* C## */
    10000, /* X   */
        4, /* D-- */  4, /* D-  */  4, /* D   */  4, /* D#  */  4, /* D## */
    10000, /* X   */
        8, /* E-- */  1, /* E-  */  1, /* E   */  8, /* E#  */  8, /* E## */
        5, /* F-- */  5, /* F-  */  5, /* F   */  5, /* F#  */  5, /* F## */
    10000, /* X   */
        9, /* G-- */  2, /* G-  */  2, /* G   */  2, /* G#  */  9, /* G## */
    10000, /* X   */
        6, /* A-- */  6, /* A-  */  6, /* A   */  6, /* A#  */  6, /* A## */
    10000, /* X   */
        3, /* B-- */  3, /* B-  */  3, /* B   */  3, /* B#  */ 10  /* B## */
   };

   static int ypos[40] = {
       -5, /* C-- */ -3, /* C-  */  0, /* C   */  1, /* C#  */  3, /* C## */
    10000, /* X   */
       -4, /* D-- */ -2, /* D-  */  0, /* D   */  2, /* D#  */  4, /* D## */
    10000, /* X   */
       -4, /* E-- */ -1, /* E-  */  1, /* E   */  2, /* E#  */  4, /* E## */
       -5, /* F-- */ -3, /* F-  */ -1, /* F   */  1, /* F#  */  3, /* F## */
    10000, /* X   */
       -5, /* G-- */ -2, /* G-  */  0, /* G   */  2, /* G#  */  3, /* G## */
    10000, /* X   */
       -4, /* A-- */ -2, /* A-  */  0, /* A   */  2, /* A#  */  4, /* A## */
    10000, /* X   */
       -3, /* B-- */ -1, /* B-  */  1, /* B   */  3, /* B#  */  4  /* B## */
   };


   // convert theta to radians:
   theta = theta * MYPI / 180.0;

   distance.setSize(40);
   double cost2 = cos(theta) * cos(theta);
   double sint2 = sin(theta) * sin(theta);
   int i;
   for (i=0; i<40; i++) {
      distance[i] = sqrt(cost2*xpos[i]*xpos[i] + sint2*ypos[i]*ypos[i]);
   }

   if (distanceQ) {
      char buffer[64] = {0};
      cout << "# theta = " << theta * 180.0 / MYPI  << "\n";
      for (i=0; i<40; i++) {
         cout << Convert::base40ToKern(buffer, i+4*40);
         cout << "\t" << distance[i] << "\n";
      }
      exit(0);
   }

}




//////////////////////////////
//
// getDistancesChromatic2 --
//

void getDistancesChromatic2(Array<double>& distances, double theta1, 
   double theta2) {

   // p1 = perfect 5th distance
   static int p1[40] = {
        6, /* C-- */  5, /* C-  */  0, /* C   */  4, /* C#  */  5, /* C## */
    10000, /* X   */
        4, /* D-- */  3, /* D-  */  2, /* D   */  3, /* D#  */  4, /* D## */
    10000, /* X   */
        6, /* E-- */  1, /* E-  */  1, /* E   */  5, /* E#  */  6, /* E## */
        5, /* F-- */  4, /* F-  */  3, /* F   */  3, /* F#  */  4, /* F## */
    10000, /* X   */
        7, /* G-- */  2, /* G-  */  1, /* G   */  2, /* G#  */  6, /* G## */
    10000, /* X   */
        5, /* A-- */  4, /* A-  */  3, /* A   */  4, /* A#  */  5, /* A## */
    10000, /* X   */
        3, /* B-- */  2, /* B-  */  2, /* B   */  3, /* B#  */  7  /* B## */
   };

   // p2 = major or minor third distance; negative are for majors
   static int p2[40] = {
       -5, /* C-- */ -3, /* C-  */  0, /* C   */  1, /* C#  */  3, /* C## */
    10000, /* X   */
       -4, /* D-- */ -2, /* D-  */  0, /* D   */  2, /* D#  */  4, /* D## */
    10000, /* X   */
       -4, /* E-- */ -1, /* E-  */  1, /* E   */  2, /* E#  */  4, /* E## */
       -5, /* F-- */ -3, /* F-  */ -1, /* F   */  1, /* F#  */  3, /* F## */
    10000, /* X   */
       -5, /* G-- */ -2, /* G-  */  0, /* G   */  2, /* G#  */  3, /* G## */
    10000, /* X   */
       -4, /* A-- */ -2, /* A-  */  0, /* A   */  2, /* A#  */  4, /* A## */
    10000, /* X   */
       -3, /* B-- */ -1, /* B-  */  1, /* B   */  3, /* B#  */  4  /* B## */
   };

   // p3 = theta 1 or theta 2
   static double p3[40] = {
        1, /* C-- */  1, /* C-  */  0, /* C   */  2, /* C#  */  2, /* C## */
    10000, /* X   */
        1, /* D-- */  1, /* D-  */  0, /* D   */  2, /* D#  */  2, /* D## */
    10000, /* X   */
        1, /* E-- */  1, /* E-  */  2, /* E   */  2, /* E#  */  2, /* E## */
        1, /* F-- */  1, /* F-  */  1, /* F   */  2, /* F#  */  2, /* F## */
    10000, /* X   */
        1, /* G-- */  1, /* G-  */  0, /* G   */  2, /* G#  */  2, /* G## */
    10000, /* X   */
        1, /* A-- */  1, /* A-  */  0, /* A   */  2, /* A#  */  2, /* A## */
    10000, /* X   */
        1, /* B-- */  1, /* B-  */  2, /* B   */  2, /* B#  */  2  /* B## */
   };

   // convert thetas to radians:
   theta1 = theta1 * MYPI / 180.0;
   theta2 = theta2 * MYPI / 180.0;

   int i;
   double theta;
   double avalue;
   double bvalue = 0.0;
   distance.setSize(40);
   for (i=0; i<40; i++) {
      avalue = p1[i] * sin(MYPI - theta1 - theta2)/sin(theta2);
      if (p2[i] < 0) {
         bvalue = -p2[i];
      } else {
         bvalue = p2[i] * sin(theta1)/sin(theta2);
      }
      if (p3[i] == 1) {
         theta = theta1;
      } else if (p3[i] == 2) {
         theta = theta2;
      } else if (p3[i] == 0) {
         theta = 0;
      } else {
         theta = MYPI;
      }
      distance[i] = sqrt(avalue * avalue + bvalue * bvalue -
            2.0 * avalue * bvalue * cos(theta));
   }

   if (distanceQ) {
      char buffer[64] = {0};
      cout << "# theta1 = " << theta1 * 180.0 / MYPI  << "\n";
      cout << "# theta2 = " << theta2 * 180.0 / MYPI  << "\n";
      for (i=0; i<40; i++) {
         cout << Convert::base40ToKern(buffer, i+4*40);
         cout << "\t" << distance[i] << "\n";
      }
      exit(0);
   }
}