function [data, wavenumbers, header] = JCSD_Matlab_Read3(filename, pulses);
% Loads data from *.dat (which is created by running JCSD_Data3.m) into
% Matlab.
%
% [data, wavenumbers, header] = JCSD_Matlab_Read(filename, pulses);
% inputs
% filename - datafile to be opened.  Should be the output of JCSD_Data3.m.
% pulses - a vector that includes the pulse number you wish to take from
%          the datafile and add to a matrix in matlab.
%  
% outputs
% data - Data values form an [Ax1024] matrix where A is the number of
%        pulses selected.
% wavenumbers - [1x1024] matrix of the wavenumbers
% header includes the following:
%   energy - Recorded energy for the selected pulse (in mJ).
%   time - Format: HHMMSSTTT where H=hour (can be 1 or 2 characters in
%        length), M=minutes, S=seconds, and T=milliseconds
%   pulseNum - Pulse number recorded by JCSD
%   lat and long - latitude and longitude from PLGR recorded by JCSD
%   DetAlarm
%   range - Value of range from rangefinder
%   RangerOK
%
% Version 1.0
% Written By: Phillip G. Wilcox
% Origination Date: 04 October 2006
% This program is designed to read in files from JCSD_Data v3.0
% 
% 21 June 2007 - Version 2.0
% Updated to read in files from JCSD_Data v3.0.  The following is added:
%       -Pulse Number
%       -Latitude and Longitude
%       -DetAlarm
%       -SensorRange and RangerOK
%
% 2 July 2007 - Version 3.0
% Increased speed by preallocating outputs and writing in columns.
% Updates so output is a structure called header to clean up output args.
% Also changed the order of input arguments so pulses does not have to be
%       entered.
tic

if nargin < 1
    filename = 'output.dat';
elseif nargin == 1
    pulses = 1;
elseif nargin > 2
    fprintf('Error.  Too many arguements.  Run help for more information.\n')
    return
end


% Open file
fid = fopen(filename, 'r','l');
if fid == -1
    fprintf('File Not Found\n')
    return
end

% Size of one line of data (data+energy+time)
dataSize = 1032;
% Calculate the total number of laser shots.
fseek(fid,0,'eof');
fileSize = ftell(fid);
fileSize = fileSize-(dataSize*4);
% 2072 = 1024 uint16 (Data) + 3 float32 (energy,lat,long) + 2 uint32
% (pulsenum,time) + 1 uint16 (range) + 2 uint8 (detAlarm,RangerOK)
scanTotal = fileSize/2072;

% Check to see if pulses is used as a vector to list the pulses desired, 
% or if it is a single number used as a spacing between pulses.
TotalPulse = size(pulses);

% If the number of pulses are not specified then the program will read in
% all pulses.
if nargin == 1
    pulses = [1:1:scanTotal];
elseif TotalPulse(1,2) == 1
    pulses = [1:pulses:scanTotal];
end

% Recheck the size to count the total number of pulses to be taken
TotalPulse = size(pulses);
% Initialize counter
pulsecount = 0;

% Read in wavenumbers
fseek(fid, 0, 'bof');
wavenumbers = [fread(fid, 1024, 'float32')]';

% Move pointer to beginning of data and start collecting the selected
% pulses.
fseek(fid,32,'cof');

% Preallocate the Data
data = zeros(1024,TotalPulse(2));
header.energy = zeros(TotalPulse(2),1);
header.time = zeros(TotalPulse(2),1);
header.pulseNum = zeros(TotalPulse(2),1);
header.lat = zeros(TotalPulse(2),1);
header.long = zeros(TotalPulse(2),1);
header.DetAlarm = zeros(TotalPulse(2),1);
header.range = zeros(TotalPulse(2),1);
header.RangerOK = zeros(TotalPulse(2),1);

for i=1:TotalPulse(2)
    %Make sure that the user has not entered a laser shot that does not
    %exist and dispay an error message if they do.
    if pulses(i) > scanTotal
        display(sprintf('Error: Data does not contain scan #%d.  Highest scan is #%d.',pulses(i),scanTotal));
        fclose(fid);
        return;
    else
        fseek(fid,(pulses(i)-pulsecount-1)*2072,'cof');
        data(:,i) = fread(fid,1024,'uint16');
        header(i,1).energy = fread(fid,1,'float32');
        header(i,1).time = fread(fid,1,'uint32');
        header(i,1).pulseNum = fread(fid,1,'uint32');
        header(i,1).lat = fread(fid,1,'float32');
        header(i,1).long = fread(fid,1,'float32');
        header(i,1).DetAlarm = fread(fid,1,'uint8');
        header(i,1).range = fread(fid,1,'uint16');
        header(i,1).RangerOK = fread(fid,1,'uint8');
        pulsecount = pulses(i);

    end
end

% Transpose Data
data = data';
% header.energy = header.energy';
% header.time = header.time';
% header.pulseNum = header.pulseNum';
% header.lat = header.lat';
% header.long = header.long';
% header.DetAlarm = header.DetAlarm';
% header.range = header.range';
% header.RangerOK = header.RangerOK';

toc
fclose(fid);