#include<iostream>
#include<math.h>
#include "GeneralFunctions.h"

using namespace std;

void cal2jdn(int a[7], double &jdn)
{
  double secs;
  int i, j, k, x, z;
  
  x = (a[1] - 14) / 12;
  z = 1461 * (a[0] + 4800 + x) / 4;
  i = 367 * (a[1] - 2 - x * 12) / 12;
  j = (a[0] + 4900 + x) / 100;
  k = 3 * j / 4;
  jdn = double(a[2]) - 32076.0 + double(z) + double(i)
    - double(k) + 1.0;
  
  secs = 3600.0 * double(a[3]) + 60.0 * double(a[4]) + double(a[5])
    + double(a[6]) / 1000.0;
  if(secs > 43200.0)
  {
    secs = secs - 43200.0;
    secs = secs / 86400.0;
    jdn += secs;
  }
  else
  {
    secs = secs / 86400.0;
    jdn += secs - 0.5;
  }
  return;
}

void jdn2cal(double jdn, int cal[7])
{
  const int dp400y = 146097;
  const int dp4000y = 1460970 + 31;  // Fudged Days
  
  int x, y, z, m, d, i;
  int jdnint = int(jdn+0.5);
  double usec;
  double fracjdn = jdn - jdnint;
  double secsofday = 86400.0 * fracjdn;
  int houraddon = 12;
  
  if(secsofday >= 43200.0)
  {
    jdnint++;
    fracjdn -= 0.5;
    secsofday -= 43200.0;
    houraddon = 0;
  }
  
  x = jdnint + 68569;
  z = 4 * x / dp400y;
  i = (dp400y * z + 3) / 4;
  x -= i;
  y = 4000 * (x + 1) / dp4000y;
  i = 1461 * y / 4;
  x += 31 - i;
  m = 80 * x / 2447;
  i = 2447 * m / 80;
  d = x - i;
  x = m / 11;
  m += 2 - 12 * x;
  y = 100 * (z - 49) + y + x;
  cal[0] = y;
  cal[1] = m;
  cal[2] = d;
/*	cal[3] = int(secsofday) / 3600;
  cal[4] = (int(secsofday) - cal[3] * 3600) / 60;
  cal[5] = cal[3] * 3600 + cal[4] * 60;
//	cal[5] = int(secsofday) - cal[3] * 3600 - cal[4] * 60;
cal[5] = int(secsofday) - cal[5];
usec = (secsofday - cal[3] * 3600 - cal[4] * 60 - cal[5]) * 1000;*/
  
  
  cal[3] = int(fracjdn * 24.0 + 0.5);
  fracjdn = fracjdn - double(cal[3]) / 24.0;
  cal[4] = int(fracjdn * 1440.0 + 0.5);
  fracjdn -= double(cal[4]) / 1440.0;
  cal[5] = int(fracjdn * 86400.0 + 0.5);
  fracjdn -= double(cal[5]) / 86400;
  usec = fracjdn * 86400000000.0;
  
  cal[6] = int(usec + 0.5);
  cal[3] += houraddon;
  
  return;
}

int cal2doy(int *calDate){
  double cjd;
  double yjd;
  int doy;
  int ycal[7] = { calDate[0], 1, 1, 0, 0, 0, 0 };
  //Calculate elapsed day of year (Jan 1 = day 001)
  cal2jdn(calDate, cjd);
  cal2jdn(ycal, yjd);
  doy = int (cjd - yjd + 1);
  
  return calDate[0]*1000 + doy;
//if ($ejd < 100) { $ejd = '0' . "$ejd"}
//elsif ($ejd < 10) { $ejd = '00' . "$ejd"};
  
//return $ejd;
}

void doy2cal(int yyyydoy, int *calDate)
{
  int ML[13] = {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  int doy = yyyydoy % 1000;
  int year = yyyydoy / 1000;
  int day, T, mon;
  mon = 0;
  day = T = doy;
  if(!(year % 4))  ML[2]++;
  while (T > 0) {
    mon++;
    day = T;
    T -= ML[mon];
  }
  calDate[0] = year;
  calDate[1] = mon;
  calDate[2] = day;
  calDate[3] = 0;
  calDate[4] = 0;
  calDate[5] = 0;
  calDate[6] = 0;
  
  return;
}




int search(EventVar <double> A, double nmcTP){
  int type;
  //cout << "ASIZE:  " << A.size() << endl;
  if(A[2]>A[1]) type=1;
  else type=0;
  //cout << "Type: " << type << endl;
  //if Event is a Sunset, make the TP Altitudes Ascending
  int j;
  EventVar <double> B(A);
  //cout << "BSIZE:  " << B.size() << endl;
  if(type == 0)
    for(int i=0; i<B.size(); i++)
      B[i] = A[A.size()-i-1];
  j = search(B, nmcTP, 0, B.size()-1);
  if(type == 0 && j != -1)
    j = A.size()-j-1;
  return j;
}

// int search(EventVar <double> A, double key, int l, int r){
//         int i;
//         //cout << "A0: " << A[0] << endl;
//         //cout << "An: " << A[r] << endl;
//         //cout << "Key: " << key << endl;
//         while(r >= l){
//                 i = (l+r)/2;
//                 //cout << "I:  " << i << "  Ai:  " << A[i] << "  L:  " << l << "  Al:  " << A[l] << "  R:  " << r << "  Ar:  " << A[r] << endl;
//                 //if(A[i] >= key-TPEPSILON && A[i]<= key+TPEPSILON) return i;
//                 if(DOUBLE_EQ(A[i], key)) return i;
// 		else if(i == 0) return -1; 
//                 else if(A[i] > key)
//                         r = i;
//                 else    l = i;
//         }
//         return (r+l)/2;
// }

int search(EventVar <double> A, double key, int l, int r)
{
  int i;
  unsigned long counter = 0;
  //cout << "A0: " << A[0] << endl;
  //cout << "An: " << A[r] << endl;
  //cout << "Key: " << key << endl;
  while(r >= l){
    i = (l+r)/2;
    //cout << "I:  " << i << "  Ai:  " << A[i] << "  L:  " << l << "  Al:  " << A[l] << "  R:  " << r << "  Ar:  " << A[r] << endl;
    //if(A[i] >= key-TPEPSILON && A[i]<= key+TPEPSILON) return i;
    if(DOUBLE_EQ(A[i], key)) 
      return i;
    else if(i == 0) 
      return -1; 
    else if(A[i] > key)
      r = i;
    else    
      l = i;
    counter++;
    if(counter >= A.size())
      return -1;                
  }
  
  return (r+l)/2;
}