//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'. 


function [jd] = CL__dat_cal2jd(cal)
  // Calendar to Julian date
  // cal = Calendar date (3xN) or (6xN)
  // [year;month;day] or [year;month;day;hour;minute;second]
  // jd = Two part Julian date (2xN)
  // NB: no test of size(cal)

  // NB: hour,minute and second can be scalars
  // => see at the end: jd2 is possibly resized. 
  // %nan is returned if %nan is present in cal

// sub-function for the algorithm (MEEUS)
function [jd] = dat_cal2jd(cal) 
  hour = 0; 
  minute = 0; 
  second = 0; 

  year = cal(1,:);
  month = cal(2,:);
  day = cal(3,:);

  if (size(cal,1) == 6)
    hour = cal(4,:);
    minute = cal(5,:);
    second = cal(6,:);
  end

  // Check inputs
  // No test on day as no impact on algorithm 
  if (find(month < 1 | month > 12 | year < -4712 | ...
           month<>round(month) | year<>round(year) | day<>round(day)) <> []);
    CL__error("Invalid calendar date");
  end

  // Algorithm (see Meeus)
  // Note on int: int(1.5)=1, int(-1.5)=-1
  k = find(month <= 2); 
  year(k) = year(k) - 1; 
  month(k) = month(k)+12; 
   
  A = int(year / 100.0); 
  B = 2 - A + int(A / 4); 

  // Julian day (Julian calendar)
  jd1 = int(365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + day - 1524; 
  jd2 = ((hour*3600 + minute*60 + second)/86400 - 0.5) .* ones(jd1); 

  k = find(jd1 + jd2 >= 2299160.5); // 2299160.5 = 5 oct 0h julian calendar
  jd1(k) = jd1(k) + B(k); // gregorian calendar

  // Date is returned as a two-part Julian date
  // jd(1,:) : integer, jd(2,:): in [0,1[
  jd = [jd1 + floor(jd2) ; jd2 - floor(jd2)]; 

endfunction

// [] => []
if (cal == [])
  jd = []; 
  return; 
end

// %nan is returned if %nan is present in input
N = size(cal, 2); 
jd = %nan * ones(2,N); 

// computation for non %nan inputs
I = find(~isnan(sum(cal,"r"))); 
if (I <> [])
  jd(:,I) = dat_cal2jd(cal(:,I)); 
end

endfunction