function [rx] = parseNSTB(filename, starthour, endhour, crccheck);
%---------------------------------------------------------
% parse NSTB formatted data to SGMP
% Ju-Yong Do,  Jan 12, 2009
%---------------------------------------------------------
% input  : NSTB filename
% output : SGMP file (.mat) and a structure containing all data
% 01/26/06 : divide and save data into smaller sub-files for speed-up
%=========================================================
%=========================================================
switch nargin()    
    case 1, starthour=0; endhour=24; crccheck=0;    %only filename is given
    case 3, crccheck=0;                             %filename and parsing time given
    case 4,                                         %all options given
    otherwise, error('Usage: parseNSTB(filename, starthour, endhour,  crccheck)');    
end
if ~( (0<=starthour) && (starthour<endhour) ), 
    error('Please, use time period that 0 <= starthour < endhour.'); 
end
parameterNSTB;
option.crc      = crccheck;                         %to save time, crc check can be skipped.
option.starttow = starthour*3600;                   %start time (sec) of day to be parsed (later converted to tow)
option.endtow   = min([endhour*3600-1 SEC_PER_DAY-1]);%end time (sec) of day to be parsed (later converted to tow)

%=========================================================
% parse NSTB data and save to SGMP format
if ~readFileInfoNSTB(filename), return; end;%general file information (sitename, rxtype, rxid, gpswk)   
fid     = fopen([filename], 'r');
h       = waitbar(0,['Searching data in ',option.fname4plot, '...']);

% parse NSTB file
while(1)    
    %find header of message (starts with sync)    
    messageId   = findMsgNSTB(fid);
    if (messageId==-1)|feof(fid), break; end; %exit if no more msg or End of file    
        
    % parse messages by calling subroutines
    switch(messageId),                
        case (1),  parseNSTBtracking(fid);  % satellite tracking data
        case (5),  parseNSTBwaas(fid);      % GEO/WAAS message    
        case (20), parseNSTBeph(fid);       % ephemeral data
        case (30), parseNSTBiono(fid);      % Klobuchar data
        case (31), parseNSTButc(fid);       % UTC data
        case (32), parseNSTBalmanac(fid);   % Almanac        
        case (40), parseNSTBtms(fid);       % TMS to TUS    
        case (41), parseNSTBtus(fid);       % TUS to TMS    
        case (50), parseNSTBsqm(fid);       % SQM data
        case (51), parseNSTBsqminfo(fid);   % SQM information
        case (100),parseNSTBwaasstatus(fid);% WAAS subsystem status data
    end    
    
    %save into .mat file (sub-files) and show parsing progress
    if saveNSTB(h, filename), break; end; %exit if parsing completed.
end
fclose(fid); close(h);

return;
%=========================================================
%=========================================================
% %recevier info
% rx.info.id          = 1;                %receiver identification number
% rx.info.sitename    = 'hawaii';         %receiver site name if applicable
% rx.info.type        = 'super';          %receiver type
% rx.info.initpos     = [0 0 0];          %initial position of receiver in XYZ (meter)
% rx.info.gpswk       = 1513;             %reference time in gps week number
% rx.info.starttow    = 0;                %start time of reception in time of week (second)
% rx.info.endtow      = 100;              %end time of reception in time of week (second)
% 
% %=========================================================
% %measurements
% %measurement set info (set of a same kind of data collected at a same time)
% rx.meas(j).type      = 'gps';            %measurement type (gps, glonass, galileo...)
% rx.meas(j).numtx     = 10;               %number of transmitters (differ from numch)
% rx.meas(j).tow       = 0;                %time of validity in time of week (second)

% %measurements per transmitter
% rx.meas(j).tx(k).prn         = 1;        %transmitter identification number. prn for GPS satellties
% rx.meas(j).tx(k).numch       = 1;        %number of channels for this transmitter
% rx.meas(j).tx(k).status      = 0;        %status flag. see NSTB document for definition
% rx.meas(j).tx(k).prL1        = 0;        %pseudo-range (meter) L1
% rx.meas(j).tx(k).crL1        = 0;        %carrier-range (meter) L1
% rx.meas(j).tx(k).dpL1        = 0;        %doppler (meter/second) L1
% rx.meas(j).tx(k).snrL1       = 0;        %signal-to-noise ratio (dBHz) L1
% rx.meas(j).tx(k).slipL1      = 0;        %cyclic slip counter L1
% rx.meas(j).tx(k).prL2        = 0;        %pseudo-range (meter) L2
% rx.meas(j).tx(k).crL2        = 0;        %carrier-range (meter) L2
% rx.meas(j).tx(k).dpL2        = 0;        %doppler (meter/second) L2
% rx.meas(j).tx(k).snrL2       = 0;        %signal-to-noise ratio (dBHz) L2
% rx.meas(j).tx(k).slipL2      = 0;        %cyclic slip counter L2
 
%=========================================================
%aiding information
% %ephemeris, clock parameters
% eph.tow     = fread(fid, 1, 'uint32')*1e-3;     %time of validity in time of week (ms->sec)
% eph.prn     = fread(fid, 1, '*uint8');          %prn number
% eph.rxtow   = fread(fid, 1, 'uint32')*1e-3;     %received time of msg (sec)
% eph.ura     = fread(fid, 1, 'uint8');           %user range accuracy (0~15)
% eph.health  = fread(fid, 1, '*uint8');          %satellite health
% eph.iodc    = fread(fid, 1, '*uint16');         %issue of data, clock
% eph.tgd     = fread(fid, 1, 'int8')*2^-31;      %group delay (sec)
% eph.toc     = fread(fid, 1, 'uint16')*2^4;      %reference time of clock parameter (second)
% eph.af2     = fread(fid, 1, 'int8')*2^-55;      %clk correction coef (s/s/s)
% eph.af1     = fread(fid, 1, 'int16')*2^-43;     %clk correction coef (s/s)
% eph.af0     = fread(fid, 1, 'int32')*2^-31;     %clk correction coef (sec)
% eph.m0      = fread(fid, 1, 'int32')*pi*2^-31;  %mean anomaly at ref time (rad)
% eph.dn      = fread(fid, 1, 'int16')*pi*2^-43;  %mean motion correction [rad/s]
% eph.e       = fread(fid, 1, 'uint32')*2^-33;    %eccentricity
% eph.sqrta   = fread(fid, 1, 'uint32')*2^-19;    %sqrt(semi-major axis) (sqrt_m)
% eph.omg0    = fread(fid, 1, 'int32')*pi*2^-31;  %longitude of ascending node (rad)
% eph.i0      = fread(fid, 1, 'int32')*pi*2^-31;  %inclination angles at ref time (rad)
% eph.w       = fread(fid, 1, 'int32')*pi*2^-31;  %argument of perigee (rad)
% eph.odot    = fread(fid, 1, 'int32')*pi*2^-43;  %rate of right ascension (rad/s)
% eph.idot    = fread(fid, 1, 'int16')*pi*2^-43;  %rate of inclination angle (rad/s)
% eph.cuc     = fread(fid, 1, 'int16')*2^-29;     %argument of latitude correciton, cosine (rad)
% eph.cus     = fread(fid, 1, 'int16')*2^-29;     %argument of latitude correciton, sine (rad)
% eph.crc     = fread(fid, 1, 'int16')*2^-5;      %radius correciton, cosine (meter)
% eph.crs     = fread(fid, 1, 'int16')*2^-5;      %radius correciton, sine (meter)
% eph.cic     = fread(fid, 1, 'int16')*2^-29;     %inclination correciton, cosine (rad)
% eph.cis     = fread(fid, 1, 'int16')*2^-29;     %inclination correciton, sine (rad)
% eph.toe     = fread(fid, 1, 'uint16')*2^4;      %reference time of ephemeris parameter (second)
% eph.iode    = fread(fid, 1, '*uint8');          %issue of data, ephemeris

% Almanac message
% almanac.tow     = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% almanac.dwn     = fread(fid, 1, 'uint8');       %difference b.t. GPS week and week of almanac (weeks)
% almanac.toa     = fread(fid, 1, 'uint8')*2^12;  %time of almanac (sec)
% almanac.numtx   = fread(fid, 1, 'uint8');       %number of satellites in the record
% almanac.prn(i)     = fread(fid, 1, '*uint8');          %prn number
% almanac.e(i)       = fread(fid, 1, 'uint16')*2^-21;    %eccentricity
% almanac.di(i)      = fread(fid, 1, 'int16')*pi*2^-19;  %delta inclination from 0.3 (rad)
% almanac.odot(i)    = fread(fid, 1, 'int16')*pi*2^-38;  %rate of right ascension (rad/s)
% almanac.sqrta(i)   = fread(fid, 1, 'uint32')*2^-11;    %sqrt(semi-major axis)
% almanac.omg0(i)    = fread(fid, 1, 'int32')*pi*2^-23;  %longitude of ascending node (rad)
% almanac.w(i)       = fread(fid, 1, 'int32')*pi*2^-23;  %argument of perigee (rad)
% almanac.m0(i)      = fread(fid, 1, 'int32')*pi*2^-23;  %mean anomaly at ref time (rad)
% almanac.af0(i)     = fread(fid, 1, 'int16')*2^-20;     %clk correction coef (sec)
% almanac.af1(i)     = fread(fid, 1, 'int16')*2^-38;     %clk correction coef (s/s)
% almanac.health(i)  = fread(fid, 1, '*uint8');          %satellite health    

% %WAAS message
% waas.tow   = 0;        %time of validity in time of week (second)
% waas.nummsg= 1;        %number of messages
% waas.prn(i)= 0;        %prn number
% waas.msg(i,:)=;        %WAAS messages

% % Iono message
% iono.tow     = fread(fid, 1, 'uint32')*1e-3;        %time of validity in time of week (ms->sec)
% iono.rxtow   = fread(fid, 1, 'uint32')*1e-3;        %received time of msg (sec)
% iono.alpha0  = fread(fid, 1, 'uint8')*2^-30;        %alpha0, amplitude coefficient (sec)
% iono.alpha1  = fread(fid, 1, 'uint8')*2^-27/pi;     %alpha1, amplitude coefficient (sec/rad)
% iono.alpha2  = fread(fid, 1, 'uint8')*2^-24/pi^2;   %alpha2, amplitude coefficient (sec/rad^2)
% iono.alpha3  = fread(fid, 1, 'uint8')*2^-24/pi^3;   %alpha3, amplitude coefficient (sec/rad^3)
% iono.beta0   = fread(fid, 1, 'uint8')*2^11;         %beta0, period coefficient (sec)
% iono.beta1   = fread(fid, 1, 'uint8')*2^14/pi;      %beta1, period coefficient (sec/rad)
% iono.beta2   = fread(fid, 1, 'uint8')*2^16/pi^2;    %beta2, period coefficient (sec/rad^2)
% iono.beta3   = fread(fid, 1, 'uint8')*2^16/pi^3;    %beta3, period coefficient (sec/rad^3)

% % utc message
% utc.tow         = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% utc.a0          = fread(fid, 1, 'uint32')*2^-30;%UTC offset coef (sec)
% utc.a1          = fread(fid, 1, 'uint32')*2^-50;%UTC offset coef (sec/sec)
% utc.dtls        = fread(fid, 1, 'uint8');       %delta time due to leap seconds (sec)
% utc.tot         = fread(fid, 1, 'uint8')*2^12;  %reference time for UTC (sec)
% utc.wnt         = fread(fid, 1, 'uint8');       %reference week number for UTC (weeks)
% utc.wnlsf       = fread(fid, 1, 'uint8');       %effectivity reference week number (weeks)
% utc.dn          = fread(fid, 1, 'uint8');       %effectivity reference day number (days)
% utc.dtlsf       = fread(fid, 1, 'uint8');       %delta time at leap seconds event (sec)

% %tms message
% tms.id          = fread(fid, 1, '*uint16');      %tms identification number
% tms.tow         = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% tms.msgcnt      = fread(fid, 1, '*uint8');      %message sequence counter (0,1, or 2)
% tms.geoprn      = fread(fid, 1, '*uint8');      %GEO prn
% tms.status      = fread(fid, 1, '*uint16');     %TMS status flags
% tms.waasmsg     = fread(fid, 32, '*uint8');     %WAAS message

% %TUS message
% tus.id      = fread(fid, 1, '*uint16');      %TUS identification number
% tus.tow     = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% tus.status  = fread(fid, 1, '*uint16');     %TUS status flags
% tus.numtms  = fread(fid, 1, '*uint8');       %number of TMSs sending message type 40s
% tus.numgps  = fread(fid, 1, '*uint8');       %number of GPS satellites tracked
% tus.starttow= fread(fid, 1, 'uint32')*1e-3; %time first bit of WAAS message broadcast to GEO in time of week (ms->sec)
% tus.endtow  = fread(fid, 1, 'uint32')*1e-3; %time last bit of WAAS message broadcast to GEO in time of week (ms->sec)
% tus.geo.prn         = fread(fid, 1, '*uint8');  %prn number
% tus.geo.waasmsg     = fread(fid, 32, '*uint8'); %FEC decoded WAAS message received from GEO
% tus.geo.pr          = fread(fid, 1, 'float64'); %pseudo-range in L1 (meter)
% tus.geo.cr          = fread(fid, 1, 'float64'); %carrier-range in L1 (meter)    
% tus.geo.dp          = fread(fid, 1, 'float32'); %doppler in L1 (meter/second)
% tus.geo.snr         = fread(fid, 1, 'float32'); %signal-to-noise ratio in L2 (dBHz)
% tus.geo.range       = fread(fid, 1, 'float64'); %computed geometric range (meter)
% tus.geo.ionodelay   = fread(fid, 1, 'float32'); %computed ionospheric delay (meter)
% tus.geo.tropodelay  = fread(fid, 1, 'float32'); %computed tropospheric delay (meter)
% tus.tms(i).id       = fread(fid, 1, '*uint16');     %tms identification number
% tus.tms(i).status   = fread(fid, 1, '*uint16');     %TMS status flags
% tus.tms(i).msgcnt   = fread(fid, 1, '*uint8');      %message sequence counter (0,1, or 2)
% tus.tms(i).rxtow    = fread(fid, 1, 'uint32')*1e-3; %time message type 40 was received from TMS (ms->sec)
% tus.gps(i).prn      = fread(fid, 1, '*uint8');  %prn number
% tus.gps(i).status   = fread(fid, 1, '*uint32'); %satellite status flags
% tus.gps(i).pr       = fread(fid, 1, 'float64'); %pseudo-range in L1 (meter)
% tus.gps(i).cr       = fread(fid, 1, 'float64'); %carrier-range in L1 (meter)    
% tus.gps(i).dp       = fread(fid, 1, 'float32'); %doppler in L1 (meter/second)
% tus.gps(i).snr      = fread(fid, 1, 'float32'); %signal-to-noise ratio in L2 (dBHz)
 
% %sqm message
% sqm.tow     = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% sqm.numprn  = fread(fid, 1, '*uint8');      %number of PRNs
% sqm.prn(i).prn      = fread(fid, 1, '*uint8');      %prn number
% sqm.prn(i).numch    = fread(fid, 1, '*uint8');      %number of channels for this prn
% sqm.prn(i).ch(j).accum1     = fread(fid, 1, 'uint32');     %accumulation 1   
% sqm.prn(i).ch(j).accum2     = fread(fid, 1, 'uint32');     %accumulation 2   
% sqm.prn(i).ch(j).accum3     = fread(fid, 1, 'uint32');     %accumulation 3   
% sqm.prn(i).ch(j).accum4     = fread(fid, 1, 'uint32');     %accumulation 4   
% sqm.prn(i).ch(j).accum5     = fread(fid, 1, 'uint32');     %accumulation 5   
% sqm.prn(i).ch(j).sync       = fread(fid, 1, '*uint8');     %sync flag (0=in sync, 1=not sync)   

% %sqminfo message
% sqminfo.tow     = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% sqminfo.rxconfig= fread(fid, 1, '*uint8');      %receiver configuration 0=narrow, 1=ultranarrow
% sqminfo.numch   = fread(fid, 1, '*uint8');      %number of channels
% sqminfo.ch(i).id        = fread(fid, 1, '*uint8');      %channel number
% sqminfo.ch(i).accumtime = fread(fid, 1, 'uint32')*1e-3; %accumulation time (ms->sec)
% sqminfo.ch(i).numintch  = fread(fid, 1, '*uint8');      %number of internal channels for this channel (<=2) 
% sqminfo.ch(i).intch(j).a1xval   = fread(fid, 1, 'int8');    %correlator spacing in C/A chips from punctual for accumulation 1
% sqminfo.ch(i).intch(j).a2xval   = fread(fid, 1, 'int8');    %correlator spacing in C/A chips from punctual for accumulation 2
% sqminfo.ch(i).intch(j).a3xval   = fread(fid, 1, 'int8');    %correlator spacing in C/A chips from punctual for accumulation 3
% sqminfo.ch(i).intch(j).a4xval   = fread(fid, 1, 'int8');    %correlator spacing in C/A chips from punctual for accumulation 4
% sqminfo.ch(i).intch(j).a5xval   = fread(fid, 1, 'int8');    %correlator spacing in C/A chips from punctual for accumulation 5                
% sqminfo.ch(i).intch(j).statusa  = fread(fid, 1, '*uint8');  %status A
% sqminfo.ch(i).intch(j).statusb  = fread(fid, 1, '*uint8');  %status B

% %waas subsystem status message
% waasstatus.weiid    = fread(fid, 1, '*uint16');     %WEI id
% waasstatus.tow      = fread(fid, 1, 'uint32')*1e-3; %time of validity in time of week (ms->sec)
% waasstatus.nummsg   = fread(fid, 1, '*uint8');      %number of messages
% waasstatus.subid(i)     = fread(fid, 1, '*uint16');     %waas subsystem id
% waasstatus.net(i)       = fread(fid, 1, '*uint8');      %network status
% waasstatus.master(i)    = fread(fid, 1, '*uint8');      %master status
% waasstatus.composite(i) = fread(fid, 1, '*uint8');      %composite fault code status    
% waasstatus.p43(i)       = fread(fid, 1, '*uint8');      %43p status    
% waasstatus.rx(i)        = fread(fid, 1, '*uint8');      %receiver status    
% waasstatus.cs(i)        = fread(fid, 1, '*uint8');      %cesium status    
% waasstatus.pcu(i)       = fread(fid, 1, '*uint8');      %power control unit status    