//  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'.// Internal function: Julian to calendar date 

function [cal] = CL__dat_jd2cal(jd)
  // Julian date to calendar date
  // jd = Two-part Julian date. jd(1,:) preferably an integer and jd(2,:) preferably in [0,1[ (2xN)
  // cal = Calendar date (6xN) [year;month;day;hour;minute;second]

  
// sub-function for the algorithm (MEEUS)
function [cal] = dat_jd2cal(jd)

  if (find(jd < 0) <> [])
    CL__error("Julian date cannot be negative");
  end

  jd1 = jd(1,:);
  jd2 = jd(2,:);
  
  // Algorithm (see Meeus) 
  // Note on int: int(1.5)=1, int(-1.5)=-1
  jd = jd1 + jd2 + 0.5; 
  Z = int(jd); 
  F = jd - Z;

  A = Z; 
  ALPHA = int((Z - 1867216.25)/36524.25);
  k = find(Z >= 2299161); 
  A(k) = Z(k) + 1 + ALPHA(k) - int(ALPHA(k)/4); 

  B = A + 1524;
  C = int((B - 122.1)/365.25);
  D = int(365.25 * C);
  E = int((B - D)/30.6001);

  RD = B - D - int(30.6001*E) + F; // decimal day

  day = floor(RD);

  month = E - 1; 
  k = find(E == 14 | E == 15); 
  month(k) = E(k) - 13; 

  year = C - 4716; 
  k = find(month == 1 | month == 2); 
  year(k) = C(k) - 4715; 
  
  // Retrieve decimal part of the day
  rh = (jd1-floor(jd1)) + (jd2-floor(jd2)) + 0.5;  
  rh = (rh - floor(rh))*24; // <=> F (see above = fractional part from 00:00) 
  hours = floor(rh);
  minutes = floor((rh-hours)*60); 
  seconds = ((rh-hours)*60 - minutes)*60; 

  // output
  cal = [year; month; day; hours; minutes; seconds];

endfunction


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

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

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

endfunction