% [R, date, jd] = readjpltraj(FILENAME) reads trajectory data from a JPL file (see below) 
% and returns the (xyz) coordinates 'R', the date vector 'date' and the julian day number 'jd'
%
% The JPL file must be of the following format (example):
% 
% 2453436.949305556 = A.D. 2005-Mar-07 10:47:00.0000 (CT)
%   1.446073698267540E+08 -3.366047598472733E+07 -7.259085732325912E+01
%~
% E.i. first the julian day number, then the date and time
% next row contains the x,y and z eclitic and geocentric coordinates 
% in km. The total length of the pair of rows should be 127 bytes
% for the program to work correctly, otherwise change in the M-file!
%
% Initial version: Magnus Billvik, Swedish Institute of Space Physics,
%Uppsala, June 2005
% Revised by: Niklas Edberg, Swedish Institute of Space Physics, Uppsala,
% Oct 2005

function [R, date, jd] = readjpltraj(filename)

fp = fopen(filename,'r');    % open the earth trajectory file

disp(strcat('Reading (JPL) TRAJECTORY data from "', filename, '"...'));

i=1;
j=1;
proc=0;

% find out the number of elements in the file by dividing the bytesize with 198

%Find out the number of rows in the file, i.e. the number of elements
linedat = fgetl(fp);
linevek = fgetl(fp);
a=length(linevek);
fseek(fp,-a,'eof');
endline = fgetl(fp);
fseek(fp,0,'bof');
linedat = fgetl(fp);
linevek = fgetl(fp);
elements=1;
while strcmp(linevek,endline)==0 %Check if the line is equal to the end line, if not, count up the element
    linedat = fgetl(fp);
    linevek = fgetl(fp);
    elements = elements + 1;
end


%elements=ceil(filesize/129); %129 for matlab6 och 127 for matlab7...
disp(strcat('elements:', int2str(elements)));
fseek(fp,0,'bof');
R=zeros([3 elements]);



while (i<=elements)

    jd(i)=fscanf(fp,'%f',1);  % moon's julian date
    temp=fscanf(fp,'%*17c %16c',1);
    temp=strrep(temp,':',' ');
    [date(1,i), date(2,i), date(3,i), date(4,i)]=strread(temp,'%f %f %f %f',4);
    fgetl(fp);

    % read the vectors
    temp=fscanf(fp,'%71c',1);  % position vectors
    [R(1,i), R(2,i), R(3,i)] = strread(temp,'%f %f %f',3);
    clear(temp);

    procx=proc;
    proc=round(10*i/elements);
    if proc~=procx
        disp(strcat(int2str(proc),'0 %'));
    end
    i=i+1;
    j=j+1;

end
