/**
 * A short piece of code to read the tycho.dat file which is the combbinded data 
 * from all the other files
 * Reads each line in the file and prints and image
 *
 * (c) Nathan Bergey 01-22-2010
 * Availible under Creative Commons Attricbution 3.0 license.
 */

BufferedReader reader;
float an;
float ngpRa;
float ngpDec;
float cos_ngpDec;
float sin_ngpDec;
static float SMALL = 1e-20;
PFont font;

/* Just a setup and no "draw" loop because we are only making one image
 */
void setup() 
{
  size(4096,2048);
  background(0);
  smooth();
  ellipseMode(CENTER);
  
  /* Change this to a font that exists on your machine to print
   * numbers on the grid. See processing's documentaion on createFont() 
   *
  font = createFont("Futura LT Bold", 40);
  textFont(font);
  */

  /* Setup for converting to Galactic Coord. */
  an          = radians(32.93192);   // Galactic long of asc node on equator
  ngpRa       = radians(192.85948);  // RA of North Galactic Pole
  ngpDec      = radians(27.12825);   // Dec of North Galactic Pole
  cos_ngpDec  = cos(ngpDec);
  sin_ngpDec  = sin(ngpDec);

  /* Open the data file. The file must be in your data folder in 
   * your sketch. To add a file from inside processing go to the menu
   * [Sketch > Add File...] 
   */
  reader = createReader("tycho.dat");
  
  // Do the work
  printStars();
  
  /* Uncomment this to draw a grid with labels on top of the image */
  //printGrid();
  
  // Save an image. It shows up in your sketch folder
  save("tycho.png");
}

/**
 * This draws a pixel or small circle to the screen.
 */
void drawStar(float ra, float dec, float vmag, float bmag)
{
  float x, y;
  float index = bmag - vmag;
  color c = indexToRgb(index, vmag);  // Get a color for the star
  stroke(c);
  fill(c);
  
  // map to screen coordinates
  x = map(ra, 0, 360, 0, width);
  y = map(dec, -90, 90, height, 0);
  
  
  /* Uncomment this to draw in Galactic Coordinates.
   *
  PVector galacticCoord = Eq2Gal(radians(ra), radians(dec));
  x = galacticCoord.x;
  
  // Put the center of the Galaxy at the center of the image
  if(x > 180)
    x -= 360;
  
  // map to screen coordinates
  x = map(x, -180, 180, 0, width);
  y = map(galacticCoord.y, -90, 90, height, 0);
  */
  
  // For bright stars draw a cicle, otherwise just a pixel
  if (vmag < 3)
  {
    float s = map(vmag, 3, -1, 2, 3);
    ellipse(x, y, s, s);
  }
  else
  {
    point(x, y);
  }
}

/**
 * Convert Equtorial Coordinates to Galactic Coordinates
 * Based on code from libastro. You are not expected to understand this.
 */
PVector Eq2Gal(float ra, float dec)
{
  float sin_dec, cos_dec, a, cos_a, sin_a, b, square, c, d; 
  float lat_gal;
  float lon_gal;
  
  cos_dec = cos(dec);
  sin_dec = sin(dec);
  a = ra - ngpRa;
  cos_a = cos(a);
  sin_a = sin(a);
  b = cos_a;
  square = (cos_dec*cos_ngpDec*b) + (sin_dec*sin_ngpDec);

  // Galactic Latitude
  lat_gal = asin(square);

  c = sin_dec - (square*sin_ngpDec);
  d = cos_dec*sin_a*cos_ngpDec;
  if (abs(d) < SMALL)
    d = SMALL;

  // Galactic Longitude
  lon_gal = atan(c/d) + an;
    
  if (d < 0) lon_gal += PI;
  if (lon_gal < 0) lon_gal += TWO_PI;
  if (lon_gal > TWO_PI) lon_gal -= TWO_PI;
  
  PVector galCoord = new PVector(degrees(lon_gal), degrees(lat_gal));
  return galCoord;
}

/**
 * Draws a grid and numbers (coordinates) on the screen
 */
void printGrid()
{
  float x, y;
  stroke(100);
  strokeWeight(2);
  fill(255);
  
  // x for Equtorial
  for(int i = 0; i <= 24; i+= 2)
  { 
    x = map(i, 0, 24, 0, width);
    line(x, 0, x, height);
    text(i, x - 60, 50);
  }

  /* Uncomment for galactic coordinates
   *
  // x for Galactic
  for(int i = -180; i <= 180; i+= 30)
  { 
    x = map(i, -180, 180, 0, width);
    line(x, 0, x, height);
    text(i, x - 60, 50);
  }
  */
 
  // y
  for(int i = -90; i < 90; i += 30)
  {
    y = map(i, -90, 90, height, 0);
    line(0, y, width, y);
    text(i, 10, y + 50);
  }
}

/**
 * Parse the dat file and draw stars
 */
void printStars()
{
  int fileLength = 2539913;  // On *nix try: $ wc -l tycho.dat
  
  // Loop through file
  for (int i = 0; i < fileLength; i++)
  {
    String line;
    try 
    {
      line = reader.readLine();
    } 
    catch (IOException e) 
    {
      e.printStackTrace();
      line = null;
    }
    
    if (line == null) 
    {
      // Stop reading because of an error or file is empty
    } 
    else 
    {
      String[] pieces = split(line, "|");  // Pipe delimited
      
      float Ra = float(pieces[2]); 
      float Dec = float(pieces[3]);
      float vmag = float(pieces[19]);
      float bmag = float(pieces[17]);
     
      drawStar(Ra, Dec, vmag, bmag);
    }
  }
}

/**
 * Method for mapping color index to a Color object
 * Based on http://www.vendian.org/mncharity/dir3/starcolor/
 * numbers are more of less made up but match
 * vendian.org's results pretty well.
 */
color indexToRgb(float index, float vmag)
{
  float bright = map(vmag, 18, -1, 0, 255);
  int r = indexToRed(index);
  int g = indexToGreen(index);
  int b = indexToBlue(index);
  
  color c = color(r, g, b, bright);
  return c;
}

int indexToRed(float index)
{  
  float a = 132.206;
  float b = 206.618;
  
  float r = (a*index) + b;
  r = constrain(r, 0, 255);
  
  return round(r);
}

int indexToGreen(float index)
{  
  if (index < 0.4)
  {
    float a = 92.9412;
    float b = 216.647;
    float g = (a*index) + b;
    g = constrain(g, 0, 255);
    return round(g);
  }
  else if (index <= 1.5)
  {
    float a = -33.7549;
    float b = 262.957;
    float g = (a*index) + b;
    g = constrain(g, 0, 255);
    return round(g);
  }
  else
  {
    float a = -564.678;
    float b = 210.636905;
    float c = 1.55;
    float g = a*((index-c)*(index-c)) + b;
    g = constrain(g, 0, 255);
    return round(g);
  }
}

int indexToBlue(float index)
{  
  if (index < 1.65)
  {
    float a = -84.1162;
    float b = 284.527;
    float bl = (a*index) + b;
    bl = constrain(bl, 0, 255);
    return round(bl);
  }
  else
  {
    float a = -530.0;
    float b = 1020.0;
    float bl = (a*index) + b;
    bl= constrain(bl, 0, 255);
    return round(bl);
  }
}