function  ctSent_extractSLC
%Extract SLC
%Zhangfeng Ma, Hohai University, jspcmazhangfeng@hhu.edu.cn
%2019/07/03
[HeadFile,OrbitFile] = ctSent_findSLC;
fprintf('---------Apply Orit File & Update SLC--------- \n');
fprintf('Zhangfeng Ma, Hohai University 2019/07/10 \n');
fprintf('jspcmazhangfeng@hhu.edu.cn  \n');
fprintf('jspcmazhangfeng@gmail.com  \n');
fprintf('--------Import Initialized Parameters------- \n');
getparm;
pause(1);
parm = load('parms.mat');
date = zeros(length(HeadFile),1);
for i=1:length(HeadFile)
fprintf('----------Apply Orbit File---------- \n');
fileDir  = [parm.workDir,'\',HeadFile(i).acquisitionDate];
Measurement = cell(3,1);
if ~exist(fileDir,'dir')
    mkdir(fileDir);
end
if  strcmp(parm.polarization,'VV')
    Measurement{1} = HeadFile(i).filepart.IW1VV.SLC;
    Measurement{2} = HeadFile(i).filepart.IW2VV.SLC;
    Measurement{3} = HeadFile(i).filepart.IW3VV.SLC;
    info(i).IW(1)  = orbit_info( HeadFile(i).filepart.IW1VV.ANNO,OrbitFile{i},3);
    info(i).IW(2)  = orbit_info( HeadFile(i).filepart.IW2VV.ANNO,OrbitFile{i},3);
    info(i).IW(3)  = orbit_info( HeadFile(i).filepart.IW3VV.ANNO,OrbitFile{i},3);
    info(i).IW(1).Swathname = 'IW1';
    info(i).IW(2).Swathname = 'IW2';
    info(i).IW(3).Swathname = 'IW3';
    info(i).date = HeadFile(i).acquisitionDate;
    date(i) = str2double(info(i).date); 
    [info(i).IW(1).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(1),3);
    [info(i).IW(2).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(2),3);
    [info(i).IW(3).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(3),3);
    warning off;
    info(i) = defineBurst_updateBurst(info(i),Measurement,parm.rectangle);
%     [savepath,infopath]  = writeBurstData(HeadFile(i).filepart.IW1VV.SLC,info);
%     slcpath(i).savepath = savepath;
%     slcpath(i).infopath = infopath;
%     slcpath(i).date = HeadFile(i).acquisitionDate;
elseif strcmp(parm.polarization,'VH')
    Measurement{1} = HeadFile(i).filepart.IW1VH.SLC;
    Measurement{2} = HeadFile(i).filepart.IW2VH.SLC;
    Measurement{3} = HeadFile(i).filepart.IW3VH.SLC;
    info(i).IW(1)  = orbit_info( HeadFile(i).filepart.IW1VH.ANNO,OrbitFile{i},3);
    info(i).IW(2)  = orbit_info( HeadFile(i).filepart.IW2VH.ANNO,OrbitFile{i},3);
    info(i).IW(3)  = orbit_info( HeadFile(i).filepart.IW3VH.ANNO,OrbitFile{i},3);
    info(i).IW(1).Swathname = 'IW1';
    info(i).IW(2).Swathname = 'IW2';
    info(i).IW(3).Swathname = 'IW3';
    info(i).date = HeadFile(i).acquisitionDate;
    date(i) = str2double(info(i).date); 
    [info(i).IW(1).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(1),3);
    [info(i).IW(2).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(2),3);
    [info(i).IW(3).orbit_para , ~ ] = ApplyOrbit(OrbitFile{i},info(i).IW(3),3);
    info(i) = defineBurst_updateBurst(info(i),Measurement,parm.rectangle);
%     [savepath,infopath]  = writeBurstData(HeadFile(i).filepart.IW3VH.SLC,info);
%     slcpath(i).savepath = savepath;
%     slcpath(i).infopath = infopath;
%     slcpath(i).date = HeadFile(i).acquisitionDate;
end
fprintf(['        ',HeadFile(i).acquisitionDate ' Finished \n']);
end
m = 1;
for i = 1:length(info)
    nsubswath(i) = length(info(i).IW);
    if nsubswath(i)
        orbitfile{m} = OrbitFile{i};
        m = m+1;
    end
end
info(nsubswath==0) = [];
date(nsubswath==0) = [];
[dateUnique,~,number] = unique(date);
for i = 1:length(dateUnique)
    idx = find(number==i);
    info_out(i) = stitchInfo(info,idx);
end
setparm('masterdate',num2str(dateUnique(ceil(length(dateUnique)/2))));
parm = load('parms.mat');
for i = 1:length(dateUnique)
    datestr = info_out(i).date;
    if strcmp(datestr,parm.masterdate)
        masterIndex = i;
    end
end
for i = 1:length(dateUnique)
    idx = find(number==i);
    for j = 1:length(info_out(i).IW)
    tic;
    [info_out(i).IW(j).orbit_para , ~ ] = ApplyOrbit(orbitfile{idx(1)},info_out(i).IW(j),3);
    info_out(masterIndex).IW(j).BurstUse(i,:) = defineOvl(info_out(masterIndex).IW(j),info_out(i).IW(j));
    t = toc;
    fprintf(['      Applying oribt file of ',info_out(i).IW(j).Swathname,' in ',info_out(i).date,' finished in %d s \n'],floor(t));
    end
%     save([parm.workDir,'\',num2str(dateUnique(i)),'\info'],'-struct','info_out');
end
slcpath = writeBurstData(info_out);
save 'slcpath.mat' slcpath;
end

function [info] = defineBurst_updateBurst(info,Measurement,rectangle)
%info: Head File
info = ctSent_locate(info,rectangle);
useNot = false(3,1);
for i = 1:3
[startBurstIndex, endBurstIndex] = BurstArea(info.IW(i),rectangle);
if startBurstIndex==0
    useNot(i) = 1; 
    continue;
end
fprintf('       Updating Burst from %d - %d to 1 - %d \n',startBurstIndex,endBurstIndex,endBurstIndex-startBurstIndex+1);
for j = startBurstIndex:endBurstIndex
   info.IW(i).slcpath{j - startBurstIndex + 1} = Measurement{i};
   info.IW(i).ROI(j - startBurstIndex + 1).firstSample = info.IW(i).sensor_para.current_window.left;
   info.IW(i).ROI(j - startBurstIndex + 1).numSample = info.IW(i).sensor_para.samplesPerBurst;
   info.IW(i).ROI(j - startBurstIndex + 1).firstLine = info.IW(i).sensor_para.current_window.top+(j-1)*info.IW(i).sensor_para.linesPerBurst;
   info.IW(i).ROI(j - startBurstIndex + 1).numLine = info.IW(i).sensor_para.linesPerBurst;
end
info.IW(i).sensor_para.first_line_time = info.IW(i).sensor_para.burstFirstLineTime(startBurstIndex);
info.IW(i).sensor_para.last_line_time = info.IW(i).sensor_para.burstLastLineTime(endBurstIndex);
info.IW(i).first_vector_time = info.IW(i).subSwath.BurstFirstLineTime(startBurstIndex);
info.IW(i).last_vector_time = info.IW(i).subSwath.BurstLastLineTime(endBurstIndex); 
%Update Info
t_Az0 = info.IW(i).sensor_para.burstFirstLineTime(startBurstIndex:endBurstIndex);
t_Az1 = info.IW(i).sensor_para.burstLastLineTime(startBurstIndex:endBurstIndex);
info.IW(i).sensor_para.burstFirstLineTime = t_Az0;
info.IW(i).sensor_para.burstLastLineTime = t_Az1;
info.IW(i).subSwath.burstFirstLineTime = t_Az0;
info.IW(i).subSwath.burstLastLineTime = t_Az1;
info.IW(i).sensor_para.t_Azi1 = t_Az0(1);
t_Az0MJD = info.IW(i).subSwath.BurstFirstLineTime(startBurstIndex:endBurstIndex);
t_Az1MJD = info.IW(i).subSwath.BurstLastLineTime(startBurstIndex:endBurstIndex);
info.IW(i).subSwath.BurstFirstLineTime = t_Az0MJD;
info.IW(i).subSwath.BurstLastLineTime = t_Az1MJD;
info.IW(i).subSwath.numOfLines = info.IW(i).subSwath.linesPerBurst*(endBurstIndex-startBurstIndex+1);
ValidSample0 = info.IW(i).subSwath.firstValidSample(startBurstIndex:endBurstIndex,:);
info.IW(i).subSwath.firstValidSample = ValidSample0;
ValidSample1 = info.IW(i).subSwath.lastValidSample(startBurstIndex:endBurstIndex,:);
info.IW(i).subSwath.lastValidSample = ValidSample1;
ValidLine0 = info.IW(i).subSwath.firstValidLine(startBurstIndex:endBurstIndex);
info.IW(i).subSwath.firstValidLine = ValidLine0;
ValidLine1 = info.IW(i).subSwath.lastValidLine(startBurstIndex:endBurstIndex);
info.IW(i).subSwath.lastValidLine = ValidLine1;
info.IW(i).subSwath.firstLineTime = info.IW(i).subSwath.BurstFirstLineTime(1);
info.IW(i).subSwath.lastLineTime = info.IW(i).subSwath.BurstLastLineTime(end);
info.IW(i).subSwath.numOfBursts = endBurstIndex-startBurstIndex+1;
info.IW(i).sensor_para.current_window.top = (startBurstIndex-1)*info.IW(i).sensor_para.linesPerBurst+1;
info.IW(i).sensor_para.current_window.bottom = endBurstIndex*info.IW(i).sensor_para.linesPerBurst;
info.IW(i).subSwath.latitude = info.IW(i).burst_lat(:,startBurstIndex:endBurstIndex);
info.IW(i).subSwath.longitude = info.IW(i).burst_lon(:,startBurstIndex:endBurstIndex);
info.IW(i).subSwath.incidence = info.IW(i).subSwath.incidenceAngle(startBurstIndex:endBurstIndex,:);
info.IW(i).subSwath.azFmRateList = info.IW(i).subSwath.azFmRateList(startBurstIndex:endBurstIndex);
info.IW(i).subSwath.DCEstimate = info.IW(i).subSwath.DCEstimate(startBurstIndex:endBurstIndex);
end
info.IW(useNot) = [];

end
function info_out = stitchInfo(info,idx)
info_out.date = info(idx(1)).date;
for i = 1:length(info(idx(1)).IW)
    firstlinetime = zeros(length(idx),1);
    lastlinetime = zeros(length(idx),1);
    info_out.IW(i).sensor_para.burstFirstLineTime = [];
    info_out.IW(i).sensor_para.burstLastLineTime = [];
    info_out.IW(i).subSwath.latitude = [];
    info_out.IW(i).subSwath.longitude = [];
    info_out.IW(i).subSwath.incidence = [];
    info_out.IW(i).subSwath.burstFirstLineTime = [];
    info_out.IW(i).subSwath.burstLastLineTime = [];
    info_out.IW(i).subSwath.BurstFirstLineTime = [];
    info_out.IW(i).subSwath.BurstLastLineTime = [];
    info_out.IW(i).subSwath.numOfLines = 0;
    info_out.IW(i).subSwath.numOfBursts = 0;
    info_out.IW(i).subSwath.firstValidSample = [];
    info_out.IW(i).subSwath.lastValidSample = [];
    info_out.IW(i).subSwath.firstValidLine = [];
    info_out.IW(i).subSwath.lastValidLine = [];
    info_out.IW(i).subSwath.azFmRateList = [];
    info_out.IW(i).subSwath.DCEstimate = [];
    info_out.IW(i).subSwath.firstValidPixel=[];
    info_out.IW(i).slcpath = [];
    info_out.IW(i).ROI = [];
    info_out.IW(i).first_vector_time = [];
    info_out.IW(i).last_vector_time = [];
    for j = 1:length(idx)
        firstlinetime(j) = info(idx(j)).IW(i).sensor_para.first_line_time;
        lastlinetime(j) = info(idx(j)).IW(i).sensor_para.first_line_time;
    end
%     for k = 1:length(idx)-1
%          if lastlinetime(idx(k))<firstlinetime(idx(k))
%              fprintf('       Images cannot be stitched!!! \n');
%              break;
%          end
%     end
    [~,tid]=sort(firstlinetime);
    for j = 1:length(idx)
        firstlinetime(j) = info(idx(tid(j))).IW(i).sensor_para.first_line_time;
        lastlinetime(j) = info(idx(tid(j))).IW(i).sensor_para.first_line_time;
        info_out.IW(i).ROI = horzcat(info_out.IW(i).ROI,info(idx(tid(j))).IW(i).ROI);
        info_out.IW(i).subSwath.firstValidPixel = info(idx(tid(j))).IW(i).subSwath.firstValidPixel;
        info_out.IW(i).sensor_para.burstFirstLineTime = horzcat(info_out.IW(i).sensor_para.burstFirstLineTime,info(idx(tid(j))).IW(i).sensor_para.burstFirstLineTime);
        info_out.IW(i).sensor_para.burstLastLineTime = horzcat(info_out.IW(i).sensor_para.burstLastLineTime,info(idx(tid(j))).IW(i).sensor_para.burstLastLineTime);
        info_out.IW(i).subSwath.latitude =  vertcat(info_out.IW(i).subSwath.latitude,info(idx(tid(j))).IW(i).subSwath.latitude);
        info_out.IW(i).subSwath.longitude =  vertcat(info_out.IW(i).subSwath.longitude,info(idx(tid(j))).IW(i).subSwath.longitude);
        info_out.IW(i).subSwath.incidence =  vertcat(info_out.IW(i).subSwath.incidence,info(idx(tid(j))).IW(i).subSwath.incidence);
        info_out.IW(i).subSwath.burstFirstLineTime = horzcat(info_out.IW(i).subSwath.burstFirstLineTime,info(idx(tid(j))).IW(i).subSwath.burstFirstLineTime);
        info_out.IW(i).subSwath.burstLastLineTime = horzcat(info_out.IW(i).subSwath.burstLastLineTime,info(idx(tid(j))).IW(i).subSwath.burstLastLineTime);
        info_out.IW(i).subSwath.BurstFirstLineTime = horzcat(info_out.IW(i).subSwath.BurstFirstLineTime,info(idx(tid(j))).IW(i).subSwath.BurstFirstLineTime);
        info_out.IW(i).subSwath.BurstLastLineTime = horzcat(info_out.IW(i).subSwath.BurstLastLineTime,info(idx(tid(j))).IW(i).subSwath.BurstLastLineTime);
        info_out.IW(i).subSwath.numOfLines  = info_out.IW(i).subSwath.numOfLines + info(idx(tid(j))).IW(i).subSwath.numOfLines;
        info_out.IW(i).subSwath.firstValidSample = vertcat(info_out.IW(i).subSwath.firstValidSample,info(idx(tid(j))).IW(i).subSwath.firstValidSample);
        info_out.IW(i).subSwath.lastValidSample = vertcat(info_out.IW(i).subSwath.lastValidSample,info(idx(tid(j))).IW(i).subSwath.lastValidSample);
        info_out.IW(i).subSwath.firstValidLine = horzcat(info_out.IW(i).subSwath.firstValidLine,info(idx(tid(j))).IW(i).subSwath.firstValidLine);
        info_out.IW(i).subSwath.lastValidLine = horzcat(info_out.IW(i).subSwath.lastValidLine,info(idx(tid(j))).IW(i).subSwath.lastValidLine);
        info_out.IW(i).subSwath.azFmRateList = horzcat(info_out.IW(i).subSwath.azFmRateList,info(idx(tid(j))).IW(i).subSwath.azFmRateList);
        info_out.IW(i).subSwath.DCEstimate = horzcat(info_out.IW(i).subSwath.DCEstimate,info(idx(tid(j))).IW(i).subSwath.DCEstimate);
        info_out.IW(i).subSwath.numOfBursts  = info_out.IW(i).subSwath.numOfBursts + info(idx(tid(j))).IW(i).subSwath.numOfBursts;
        info_out.IW(i).slcpath = horzcat(info_out.IW(i).slcpath ,info(idx(tid(j))).IW(i).slcpath);
%         info_out.IW(i).ROI = horzcat(info_out.IW(i).ROI,info(idx(tid(j))).IW(i).ROI);
    end
info_out.IW(i).sensor_para.Wavelength = info(idx(tid(1))).IW(i).sensor_para.Wavelength; 
info_out.IW(i).sensor_para.range_band_width = info(idx(tid(1))).IW(i).sensor_para.range_band_width; 
info_out.IW(i).sensor_para.azimuth_band_width = info(idx(tid(1))).IW(i).sensor_para.azimuth_band_width;
info_out.IW(i).sensor_para.current_window = info(idx(1)).IW(i).sensor_para.current_window; 
info_out.IW(i).sensor_para.t_Range1 = info(idx(1)).IW(i).sensor_para.t_Range1; 
info_out.IW(i).sensor_para.RSR2x = info(idx(1)).IW(i).sensor_para.RSR2x; 
info_out.IW(i).sensor_para.t_Azi1 = info(idx(1)).IW(i).sensor_para.t_Azi1; 
info_out.IW(i).sensor_para.PRF = info(idx(1)).IW(i).sensor_para.PRF;
info_out.IW(i).sensor_para.ApproxCenterOriginal = info(idx(1)).IW(i).sensor_para.ApproxCenterOriginal;
info_out.IW(i).sensor_para.line_time_interval = info(idx(1)).IW(i).sensor_para.line_time_interval; 
info_out.IW(i).sensor_para.rangePixelSpacing = info(idx(1)).IW(i).sensor_para.rangePixelSpacing; 
info_out.IW(i).sensor_para.azimuthPixelSpacing = info(idx(1)).IW(i).sensor_para.azimuthPixelSpacing; 
info_out.IW(i).sensor_para.linesPerBurst = info(idx(1)).IW(i).sensor_para.linesPerBurst; 
info_out.IW(i).sensor_para.first_line_time = info(idx(1)).IW(i).sensor_para.first_line_time; 
info_out.IW(i).sensor_para.last_line_time = info(idx(end)).IW(i).sensor_para.last_line_time; 
info_out.IW(i).first_vector_time = info(idx(1)).IW(i).first_vector_time; 
info_out.IW(i).last_vector_time = info(idx(end)).IW(i).last_vector_time; 
info_out.IW(i).sensor_para.samplesPerBurst = info(idx(1)).IW(i).sensor_para.samplesPerBurst;
info_out.IW(i).subSwath.numOfSamples = info(idx(1)).IW(i).subSwath.numOfSamples;
info_out.IW(i).subSwath.azimuthTimeInterval = info(idx(1)).IW(i).subSwath.azimuthTimeInterval;
info_out.IW(i).subSwath.slrTimeToFirstPixel = info(idx(1)).IW(i).subSwath.slrTimeToFirstPixel;
info_out.IW(i).subSwath.slrTimeToLastPixel = info(idx(1)).IW(i).subSwath.slrTimeToLastPixel;
info_out.IW(i).subSwath.linesPerBurst = info(idx(1)).IW(i).subSwath.linesPerBurst;
info_out.IW(i).subSwath.samplesPerBurst = info(idx(1)).IW(i).subSwath.samplesPerBurst ;
info_out.IW(i).subSwath.azimuthSteeringRate = info(idx(1)).IW(i).subSwath.azimuthSteeringRate;
info_out.IW(i).subSwath.firstLineTime = info(idx(1)).IW(i).subSwath.firstLineTime;
info_out.IW(i).subSwath.lastLineTime = info(idx(end)).IW(i).subSwath.lastLineTime;
info_out.IW(i).subSwath.rangePixelSpacing = info(idx(end)).IW(i).subSwath.rangePixelSpacing;
info_out.IW(i).Swathname = info(idx(end)).IW(i).Swathname;
end   
end
function [slcpath] = writeBurstData(info)
parm = load('parms');
for i = 1:length(info)
    datestr = info(i).date;
    if strcmp(datestr,parm.masterdate)
        masterIndex = i;
    end
end
for k = 1:length(info)
    fprintf(['----------Start Writing ',info(k).date,'----------','\n']);
for i = 1:length(info(masterIndex).IW)
    fileDir  = [parm.workDir,'\',info(k).date,'\',info(k).IW(i).Swathname];
    if ~exist(fileDir,'dir')
       mkdir(fileDir);
    end
    nBurst = info(masterIndex).IW(i).subSwath.numOfBursts;
for j = 1:nBurst
    if info(masterIndex).IW(i).BurstUse(k,j) == 0
        continue;
    end  
    filePath = info(k).IW(i).slcpath{info(masterIndex).IW(i).BurstUse(k,j)};
    roi = info(k).IW(i).ROI(info(masterIndex).IW(i).BurstUse(k,j));
    fprintf([' Reading :',info(k).IW(i).Swathname, ' Burst %d of SLC\n'],j);
    data = readSent1Data(filePath,roi);
    fprintf([' Writing :',info(k).IW(i).Swathname, ' Burst %d of SLC\n'],j);
    slcpath(k).savepath{i,j} = [parm.workDir,'\',info(k).date,'\',info(k).IW(i).Swathname,'\Burst',num2str(j),'.slc'];
    slcpath(k).date = info(k).date;
    writecpx(slcpath(k).savepath{i,j},conj(data'));
end
end
infosave = info(k);
slcpath(k).infopath = [parm.workDir,'\',info(k).date,'\info'];
save(slcpath(k).infopath,'infosave');
end
end
function BurstUse = defineOvl(IWm,IWs)
mlatitude = IWm.subSwath.latitude;
mlongitude = IWm.subSwath.longitude;
slatitude = IWs.subSwath.latitude;
slongitude = IWs.subSwath.longitude;
mburst = IWm.subSwath.numOfBursts;
sburst = IWs.subSwath.numOfBursts;
BurstUse = zeros(mburst,1);
for iovl = 1:mburst
    latBurst = [mlatitude(:,iovl);mlatitude(1,iovl)];
    lonBurst = [mlongitude(:,iovl);mlongitude(1,iovl)];
    [latBurst, lonBurst] = poly2cw(latBurst,lonBurst);
i = 1;
 while i <= sburst
      lat = [slatitude(:,i);slatitude(1,i)];
      lon = [slongitude(:,i);slongitude(1,i)];
      [lat, lon] = poly2cw(lat,lon);
     poly1 = polyshape(lat,lon);
     poly2 = polyshape(latBurst,lonBurst);
%      [xb, yb] = polybool('intersection', lat, lon, latBurst, lonBurst);
     polyi = intersect(poly1,poly2);
%      [xb, yb] = poly2cw(xb, yb);
     pera = abs((area(polyi) - area(poly2)))/area(poly2);
%      perb = abs((polyarea(xb,yb) - polyarea(latBurst,lonBurst)))/polyarea(latBurst,lonBurst);
     if pera<0.2&BurstUse(iovl)==0
         BurstUse(iovl) = i;
     end
      i = i+1;
 end
end
end
function [startBurstIndex, endBurstIndex]= BurstArea(info,rectangle)
% [Lat,Lon] = initial_radar_code(info_m,mburstIndex);
latitude = info.burst_lat;
longitude = info.burst_lon;
polyrect = polyshape([rectangle(2),rectangle(2), rectangle(1), rectangle(1)],...
    [rectangle(3), rectangle(4),rectangle(4), rectangle(3)]);
nburst = info.subSwath.numOfBursts;
i = 1;
startBurstIndex = 0;
endBurstIndex = 0;
 while i <= nburst
      lat = latitude(:,i);
      lon = longitude(:,i);
      polyBurst = polyshape(lat,lon);
      polyintersect = intersect(polyBurst,polyrect);
      if startBurstIndex == 0 & polyintersect.NumRegions ~= 0
          startBurstIndex = i;endBurstIndex = i;
      elseif  endBurstIndex ~= 0 & polyintersect.NumRegions ~= 0
          endBurstIndex = i;
      end
      i = i+1;
 end
end
function [ info_m ] = orbit_info( patha,orbit_path,polyDegree)
tic
try
    xdoc = xmlread(patha);
catch
    error('Error in Sentinel-1 Annotation file(): failed to read XML file ''%s''!', a);
end
first_line_time = xdoc.getElementsByTagName('startTime').item(0).getFirstChild.getData;
last_line_time = xdoc.getElementsByTagName('stopTime').item(0).getFirstChild.getData;
lineinterval = xdoc.getElementsByTagName('azimuthTimeInterval').item(0).getFirstChild.getData;
t_Range1 = xdoc.getElementsByTagName('slantRangeTime').item(0).getFirstChild.getData;
item = xdoc.getElementsByTagName('rangeProcessing').item(0);
range_band_width = getDataByTagName(item, 'totalBandwidth', 1);
item = xdoc.getElementsByTagName('azimuthProcessing').item(0);
azimuth_band_width = getDataByTagName(item, 'totalBandwidth', 1);
item = xdoc.getElementsByTagName('productInformation').item(0);
radar_frequency = getDataByTagName(item, 'radarFrequency', 1);
item = xdoc.getElementsByTagName('downlinkInformation').item(0);
pulse_repeat_frequency = getDataByTagName(item, 'prf', 1);
item = xdoc.getElementsByTagName('productInformation').item(0);
range_sampling_rate = getDataByTagName(item, 'rangeSamplingRate', 1);
rangePixelSpacing = xdoc.getElementsByTagName('rangePixelSpacing').item(0).getFirstChild.getData;
azimuthPixelSpacing = xdoc.getElementsByTagName('azimuthPixelSpacing').item(0).getFirstChild.getData;
lines = xdoc.getElementsByTagName('numberOfLines').item(0).getFirstChild.getData;
Samples = xdoc.getElementsByTagName('numberOfSamples').item(0).getFirstChild.getData;
Heading = xdoc.getElementsByTagName('platformHeading').item(0).getFirstChild.getData;
%--------------------------------------------------%
info.PRF= str2double(char(pulse_repeat_frequency));
info.first_line_time=getSeconds(char(first_line_time));
info.last_line_time=getSeconds(char(last_line_time));
info_m.first_vector_time=UTC2MJD(char(first_line_time));
info_m.last_vector_time=UTC2MJD(char(last_line_time));
info.Centriod_time=(info.first_line_time+info.last_line_time)/2;
info.RSR2x = str2double(char(range_sampling_rate))*2;
info.first_range=str2double(char(t_Range1))/2;
info.line_time_interval=str2double(char(lineinterval));
info.radar_frequency=str2double(char(radar_frequency));
info.rangePixelSpacing=str2double(char(rangePixelSpacing));
%info.azimuth_height= str2double(char(lines));
%info.range_width= str2double(char(Samples));
info.azimuth_height=str2double(char(lines));
info.range_width= str2double(char(Samples));
linesPerBurst = str2double(char((xdoc.getElementsByTagName('linesPerBurst').item(0).getFirstChild.getData)));
samplesPerBurst=str2double(char((xdoc.getElementsByTagName('samplesPerBurst').item(0).getFirstChild.getData)));
Nodes= xdoc.getElementsByTagName('burst');
Num_burst=Nodes.getLength();
for i = 1:Num_burst
  item = xdoc.getElementsByTagName('burst').item(i-1);
  burstLineTime(i)  = getDataByTagName(item, 'azimuthTime', 1);
  burstFirstLineTime(i) = getSeconds(char(burstLineTime(i)));
  burstLastLineTime(i) = burstFirstLineTime(i) +(linesPerBurst - 1) *info.line_time_interval;
  BurstFirstLineTime(i) = MJD(char(burstLineTime(i)));
  BurstLastLineTime(i) = BurstFirstLineTime(i) +(linesPerBurst - 1) *info.line_time_interval/86400;
end

info_m.subSwath=getlocgrid( xdoc );
info.coord=BLH2XYZ(mean(info_m.subSwath.latitude(:)),mean(info_m.subSwath.longitude(:)),0);
info_m.subSwath.burstFirstLineTime=burstFirstLineTime;
info_m.subSwath.burstLastLineTime=burstLastLineTime;
info_m.subSwath.BurstFirstLineTime=BurstFirstLineTime;
info_m.subSwath.BurstLastLineTime=BurstLastLineTime;
info_m.subSwath.numOfSamples=str2double(char((xdoc.getElementsByTagName('numberOfSamples').item(0).getFirstChild.getData)));
info_m.subSwath.numOfLines = str2double(char((xdoc.getElementsByTagName('numberOfLines').item(0).getFirstChild.getData)));
info_m.subSwath.azimuthTimeInterval=info.line_time_interval;
info_m.subSwath.slrTimeToFirstPixel = str2double(char((xdoc.getElementsByTagName('slantRangeTime').item(0).getFirstChild.getData)))/2;
info_m.subSwath.slrTimeToLastPixel = info_m.subSwath.slrTimeToFirstPixel +(info_m.subSwath.numOfSamples- 1) * info.rangePixelSpacing / 299792458.0;
info_m.subSwath.linesPerBurst=linesPerBurst;
info_m.subSwath.samplesPerBurst=samplesPerBurst;
info_m.subSwath.numOfBursts=Num_burst; 
info_m.subSwath.azimuthSteeringRate = str2double(char((xdoc.getElementsByTagName('azimuthSteeringRate').item(0).getFirstChild.getData)));
[ firstValidsample,lastValidsample,firstValidLine,lastValidLine,firstValidPixel,lastValidPixel] = getValidSample( xdoc,info_m.subSwath );
info_m.subSwath.firstValidSample=firstValidsample;
info_m.subSwath.lastValidSample=lastValidsample;
info_m.subSwath.firstValidLine=firstValidLine;
info_m.subSwath.lastValidLine=lastValidLine;
info_m.subSwath.firstValidPixel=firstValidPixel;
info_m.subSwath.lastValidPixel=lastValidPixel;
info_m.subSwath.azFmRateList = getAzimuthFmRateList(xdoc);
info_m.subSwath.DCEstimate = getDCEstimate( xdoc);
% [orbit_para,OSV] = ApplyOrbit( orbit_path,info ,polyDegree);
% info_m.orbit_para = orbit_para;
% info_m.subSwath.OSV = OSV;%get_orbit(orbit_path);

info_m.subSwath.dopplerCentroid = info_m.subSwath.DCEstimate(1).c0;
info_m.subSwath.firstLineTime=info_m.first_vector_time;
info_m.subSwath.lastLineTime=info_m.last_vector_time;
info_m.sensor_para.Wavelength=299792458/info.radar_frequency;
info_m.sensor_para.range_band_width=str2double(range_band_width);
info_m.sensor_para.azimuth_band_width=str2double(azimuth_band_width);
current_window.top=1;% also used during polyval
current_window.bottom=info.azimuth_height;% also used during polyval
current_window.left=1;% also used during polyval
current_window.right=info.range_width;% also used during polyval
info_m.sensor_para.current_window=current_window;
info_m.sensor_para.t_Range1=str2double(char(t_Range1))/2;
info_m.sensor_para.RSR2x=info.RSR2x;
info_m.sensor_para.t_Azi1=info.first_line_time;
info_m.sensor_para.PRF=info.PRF;
info_m.sensor_para.ApproxCenterOriginal.x=info.coord.X;
info_m.sensor_para.ApproxCenterOriginal.y=info.coord.Y;
info_m.sensor_para.ApproxCenterOriginal.z=info.coord.Z;
info_m.sensor_para.line_time_interval=info.line_time_interval;
info_m.sensor_para.rangePixelSpacing=info.rangePixelSpacing;
info_m.sensor_para.azimuthPixelSpacing=str2double(char(azimuthPixelSpacing));
info_m.sensor_para.Centriod_time=info.Centriod_time;
info_m.sensor_para.linesPerBurst=linesPerBurst;
info_m.sensor_para.first_line_time=info.first_line_time;
info_m.sensor_para.last_line_time=info.last_line_time;
info_m.sensor_para.samplesPerBurst=samplesPerBurst;
info_m.sensor_para.burstFirstLineTime=burstFirstLineTime;
info_m.sensor_para.burstLastLineTime=burstLastLineTime;
info_m.sensor_para.Heading = Heading;
info_m.subSwath.rangePixelSpacing=info_m.sensor_para.rangePixelSpacing;
t=toc;
fprintf('       Read Annotation File of Sentinel-1 Finished in %d s \n',floor(t));
end
function azFmRateList = getAzimuthFmRateList(xdoc)
item = xdoc.getElementsByTagName('azimuthFmRate');
num=getLength(item);
for i=1 : num
    item = xdoc.getElementsByTagName('azimuthFmRate').item(i-1);
    azi  = getDataByTagName(item, 'azimuthTime', 1);
    t0  = getDataByTagName(item, 't0', 1);
    azimuthFmRatePolynomial  = getDataByTagName(item, 'azimuthFmRatePolynomial', 1);
    coeff= regexp(char(azimuthFmRatePolynomial), '\s+', 'split');
    coeff=str2double(coeff);
    t0=str2double(char(t0));
    azimuthTime = getSeconds(char(azi));
    azFmRateList(i).time=azimuthTime;
    azFmRateList(i).t0=t0;
    azFmRateList(i).c0=coeff(1);
    azFmRateList(i).c1=coeff(2);
    azFmRateList(i).c2=coeff(3);
end
end
function [ azFmRateList ] = getDCEstimate( xdoc)
item = xdoc.getElementsByTagName('dcEstimate');
num=getLength(item);
for i=1 : num
    item = xdoc.getElementsByTagName('dcEstimate').item(i-1);
    azi  = getDataByTagName(item, 'azimuthTime', 1);
    t0  = getDataByTagName(item, 't0', 1);
    azimuthFmRatePolynomial  = getDataByTagName(item, 'dataDcPolynomial', 1);
    coeff= regexp(char(azimuthFmRatePolynomial), '\s+', 'split');
    coeff=str2double(coeff);
    t0=str2double(char(t0));
    azimuthTime = getSeconds(char(azi));
    azFmRateList(i).time=azimuthTime;
    azFmRateList(i).t0=t0;
    azFmRateList(i).c0=coeff(1);
    azFmRateList(i).c1=coeff(2);
    azFmRateList(i).c2=coeff(3);
end
end
function [ orbit_para , OSV ] = ApplyOrbit( orbit_path,info ,polyDegree)
xDoc=xmlread(orbit_path);
xRoot=xDoc.getDocumentElement();
Nodes= xRoot.getElementsByTagName('UTC');
Nodes_X= xRoot.getElementsByTagName('X');
Nodes_Y= xRoot.getElementsByTagName('Y');
Nodes_Z= xRoot.getElementsByTagName('Z');
Nodes_VX= xRoot.getElementsByTagName('VX');
Nodes_VY= xRoot.getElementsByTagName('VY');
Nodes_VZ= xRoot.getElementsByTagName('VZ');
Num_Nodes=Nodes.getLength();
Node1=Nodes.item(0);
Node1_X=Nodes_X.item(0);
Node1_Y=Nodes_Y.item(0);
Node1_Z=Nodes_Z.item(0);
Node1_VX=Nodes_VX.item(0);
Node1_VY=Nodes_VY.item(0);
Node1_VZ=Nodes_VZ.item(0);
TIME = cell(Num_Nodes,1);
P_X = TIME;P_Y = TIME;P_Z = TIME;
V_X = TIME;V_Y = TIME;V_Z = TIME;
for i=1:Num_Nodes
    time_UTC=char(Node1.getFirstChild.getData);
    X_unit=char(Node1_X.getTextContent());
    Y_unit=char(Node1_Y.getTextContent());
    Z_unit=char(Node1_Z.getTextContent());
    VX_unit=char(Node1_VX.getTextContent());
    VY_unit=char(Node1_VY.getTextContent());
    VZ_unit=char(Node1_VZ.getTextContent());
    TIME{i,1}=time_UTC;
    P_X{i,1}=X_unit;
    P_Y{i,1}=Y_unit;
    P_Z{i,1}=Z_unit;
    V_X{i,1}=VX_unit;
    V_Y{i,1}=VY_unit;
    V_Z{i,1}=VZ_unit;
Node1=Nodes.item(i);
Node1_X=Nodes_X.item(i);
Node1_Y=Nodes_Y.item(i);
Node1_Z=Nodes_Z.item(i);
Node1_VX=Nodes_VX.item(i);
Node1_VY=Nodes_VY.item(i);
Node1_VZ=Nodes_VZ.item(i);
end
timeUTC   = getSeconds_UTC(char(TIME)); % <s>
inputX  = str2double(cellstr(P_X)); % <m>
inputY   = str2double(cellstr(P_Y)); % <m>
inputZ   = str2double(cellstr(P_Z)); % <m>
inputVX  = str2double(cellstr(V_X)); % <m>
inputVY   = str2double(cellstr(V_Y)); % <m>
inputVZ   = str2double(cellstr(V_Z)); % <m>
OSV.TIME = timeUTC;
OSV.P_X = inputX;
OSV.P_Y = inputY;
OSV.P_Z = inputZ;
inputMJDtime   = getMJD(char(TIME));
delta=1/86400;
numExtraVectors=10;
firstVectorTime = info.first_vector_time-numExtraVectors*delta;
lastVectorTime = info.last_vector_time+numExtraVectors*delta;
numVectors = floor((lastVectorTime - firstVectorTime)/delta);
Time = zeros(numVectors,1);X = Time; Y = Time; Z = Time;
VX = Time; VY = Time; VZ = Time;
for  i = 0: (numVectors-1)
     time = firstVectorTime + i/86400;
     Time(i+1)=info.sensor_para.first_line_time-10+i;
     [X(i+1) ,Y(i+1), Z(i+1),VX(i+1) ,VY(i+1), VZ(i+1)]=getOrbitData(inputX,inputY,inputZ,...
         inputVX,inputVY,inputVZ,inputMJDtime,time,polyDegree);
end
orbit_para.orbit_points(:, 1)=Time;
orbit_para.orbit_points(:, 2)=X;
orbit_para.orbit_points(:, 3)=Y;
orbit_para.orbit_points(:, 4)=Z;
orbit_para.orbit_points(:, 5)=VX;
orbit_para.orbit_points(:, 6)=VY;
orbit_para.orbit_points(:, 7)=VZ;
end
function [ sec ] = getSeconds( time )

sec = str2double(cellstr(time(:,12:13)))*3600 ...
      + str2double(cellstr(time(:,15:16)))*60 ...
      + str2double(cellstr(time(:,18:26)));

% end function getSeconds()
end
function [ data ] = getDataByTagName( item, tagname, i )

thisItem = item.getElementsByTagName(tagname).item(i-1);
data = thisItem.getFirstChild.getData;

% end function getDataByTagName()
end
function [ orbit ] = get_orbit( orbit_path)
xDoc=xmlread(orbit_path);
xRoot=xDoc.getDocumentElement();
Nodes= xRoot.getElementsByTagName('UTC');
Nodes_X= xRoot.getElementsByTagName('X');
Nodes_Y= xRoot.getElementsByTagName('Y');
Nodes_Z= xRoot.getElementsByTagName('Z');
Num_Nodes=Nodes.getLength();
Node1=Nodes.item(0);
Node1_X=Nodes_X.item(0);
Node1_Y=Nodes_Y.item(0);
Node1_Z=Nodes_Z.item(0);
TIME = zeros(Num_Nodes,1);
P_X = TIME;
P_Y = TIME;
P_Z = TIME;
for i=1:Num_Nodes
    time_UTC=char(Node1.getFirstChild.getData);
    X_unit=char(Node1_X.getTextContent());
    Y_unit=char(Node1_Y.getTextContent());
    Z_unit=char(Node1_Z.getTextContent());
    TIME(i)=getSeconds_UTC(time_UTC);
    P_X(i)=str2double(X_unit);
    P_Y(i)=str2double(Y_unit);
    P_Z(i)=str2double(Z_unit);
Node1=Nodes.item(i);
Node1_X=Nodes_X.item(i);
Node1_Y=Nodes_Y.item(i);
Node1_Z=Nodes_Z.item(i);
end
orbit.TIME=TIME;
orbit.P_X=P_X;
orbit.P_Y=P_Y;
orbit.P_Z=P_Z;
end
function [mjd ] = UTC2MJD( time )
year = str2double(cellstr(time(:,1:4)));
month= str2double(cellstr(time(:,6:7)));
day= str2double(cellstr(time(:,9:10)));
hour= str2double(cellstr(time(:,12:13)));
minute= str2double(cellstr(time(:,15:16)));
seconds= str2double(cellstr(time(:,18:26)));
mjd=mjuliandate(year,month,day,hour,minute,seconds);
end
function [ mjd ] = MJD( time )
year = str2double(cellstr(time(:,1:4)));
month= str2double(cellstr(time(:,6:7)));
day= str2double(cellstr(time(:,9:10)));
hour= str2double(cellstr(time(:,12:13)));
minute= str2double(cellstr(time(:,15:16)));
seconds= str2double(cellstr(time(:,18:26)));
mjd=mjuliandate(year,month,day,hour,minute,seconds);
end
function [ subSwath ] = getlocgrid( xdoc )
item = xdoc.getElementsByTagName('geolocationGridPoint');
numOfGeoLocationGridPoints=getLength(item);
numOfGeoPointsPerLine = 0;
line = 0;
for i=1 : numOfGeoLocationGridPoints
    item = xdoc.getElementsByTagName('geolocationGridPoint').item(i-1);
    if (numOfGeoPointsPerLine == 0) 
        line = str2double(getDataByTagName(item, 'line', 1));
        numOfGeoPointsPerLine=numOfGeoPointsPerLine+1;
    elseif (line == str2double(getDataByTagName(item, 'line', 1)))
           numOfGeoPointsPerLine=numOfGeoPointsPerLine+1;
    else 
         break;
    end
end
numOfGeoLines = numOfGeoLocationGridPoints / numOfGeoPointsPerLine;
%         if (numOfGeoLines <= 9) 
%             numOfGeoLines = 9 + 1;%subSwath.numOfBursts + 1;
%         end
        subSwath.numOfGeoLines = numOfGeoLines;
        subSwath.numOfGeoPointsPerLine = numOfGeoPointsPerLine;
        subSwath.azimuthTime = zeros(numOfGeoLines,numOfGeoPointsPerLine);
        subSwath.slantRangeTime = zeros(numOfGeoLines,numOfGeoPointsPerLine);
        subSwath.latitude = zeros(numOfGeoLines,numOfGeoPointsPerLine);
        subSwath.longitude = zeros(numOfGeoLines,numOfGeoPointsPerLine);
        subSwath.incidenceAngle = zeros(numOfGeoLines,numOfGeoPointsPerLine);
        k = 0;
        for order=1 : numOfGeoLocationGridPoints
            item = xdoc.getElementsByTagName('geolocationGridPoint').item(order-1);
            i = floor(k / numOfGeoPointsPerLine);
            j = k - i * numOfGeoPointsPerLine;
            subSwath.azimuthTime(i+1,j+1) = getSeconds(char(getDataByTagName(item, 'azimuthTime', 1)));%getMJD()*Constants.secondsInDay;
            subSwath.slantRangeTime(i+1,j+1) = str2double(getDataByTagName(item, 'slantRangeTime', 1))/ 2.0;
            subSwath.latitude(i+1,j+1) = str2double(getDataByTagName(item, 'latitude', 1));
            subSwath.longitude(i+1,j+1) = str2double(getDataByTagName(item, 'longitude', 1));
            subSwath.incidenceAngle(i+1,j+1) =str2double(getDataByTagName(item, 'incidenceAngle', 1));
            k=k+1;
        end
end
function [X, Y ,Z,VX,VY,VZ] = getOrbitData(X,Y,Z,VX,VY,VZ,MJDtime, utc,polyDegree )
n=length(MJDtime);
%firstVectorTime = info.first_vector_time-numExtraVectors*delta;
%lastVectorTime = info.last_vector_time+numExtraVectors*delta;
%numVectors = floor((lastVectorTime - firstVectorTime) / delta);
t0=MJDtime(1);tN=MJDtime(n);numVecPolyFit=polyDegree+1;
%t0=TIME(1);tN=TIME(n);numVecPolyFit=polyDegree+1;
halfNumVecPolyFit = numVecPolyFit / 2;
vecIdx = floor(((utc - t0) / (tN - t0) * (n-1)));
vectorIndices = zeros(numVecPolyFit,1);
if (vecIdx <= halfNumVecPolyFit - 1)
    
            for i=1:numVecPolyFit
                vectorIndices(i) = i;
            end
         elseif (vecIdx >= n - halfNumVecPolyFit) 
            for i = 1:numVecPolyFit 
                vectorIndices(i) = n - numVecPolyFit + i;
            end
         else 
            for i = 0: numVecPolyFit-1
                vectorIndices(i+1) = vecIdx - halfNumVecPolyFit  + i+2;
            end
end
timeArray = zeros(numVecPolyFit,1);
xPosArray =timeArray;
yPosArray = timeArray;
zPosArray = timeArray;
VxPosArray = timeArray;
VyPosArray = timeArray;
VzPosArray = timeArray;
for  i=1: numVecPolyFit
            timeArray(i) = MJDtime(vectorIndices(i)) - t0;%MJD
            xPosArray(i) = X(vectorIndices(i));
            yPosArray(i) = Y(vectorIndices(i));
            zPosArray(i) = Z(vectorIndices(i));
            VxPosArray(i) = VX(vectorIndices(i));
            VyPosArray(i) = VY(vectorIndices(i));
            VzPosArray(i) = VZ(vectorIndices(i));
end
A=createVandermondeMatrix(timeArray,polyDegree);
[L,U] = lu(A);
rhsX=U\(L\xPosArray);
rhsY=U\(L\yPosArray);
rhsZ=U\(L\zPosArray);
rhsVX=U\(L\VxPosArray);
rhsVY=U\(L\VyPosArray);
rhsVZ=U\(L\VzPosArray);
normalizedTime = utc - t0;
val_X= poly_Val( normalizedTime,rhsX );
val_Y= poly_Val( normalizedTime,rhsY );
val_Z= poly_Val( normalizedTime,rhsZ );
val_VX= poly_Val( normalizedTime,rhsVX );
val_VY= poly_Val( normalizedTime,rhsVY );
val_VZ= poly_Val( normalizedTime,rhsVZ );
X=val_X;
Y=val_Y;
Z=val_Z;
VX=val_VX;
VY=val_VY;
VZ=val_VZ;
end
function [ sec ] = getSeconds_UTC( time )

sec = str2double(cellstr(time(:,16:17)))*3600 ...
      + str2double(cellstr(time(:,19:20)))*60 ...
      + str2double(cellstr(time(:,22:30)));
% end function getSeconds()

end
function [ mjd ] = getMJD( time )
year = str2double(cellstr(time(:,5:8)));
month= str2double(cellstr(time(:,10:11)));
day= str2double(cellstr(time(:,13:14)));
hour= str2double(cellstr(time(:,16:17)));
minute= str2double(cellstr(time(:,19:20)));
seconds= str2double(cellstr(time(:,22:30)));
mjd=mjuliandate(year,month,day,hour,minute,seconds);
end
function [ array ] = createVandermondeMatrix( time,polydegree )
n=length(time);
array=zeros(n,polydegree+1);
for i=1:n
    for j=1:(polydegree+1)
        array(i,j)=power(time(i),j-1);
    end
end

end
function [ val ] = poly_Val( t,coeff )
val=0.0;
i = length(coeff);
        while ( i > 0)
            val = val .* t + coeff(i);
            i=i-1;
        end
end
function [ firstValidSample,lastValidSample,firstValidLine,lastValidLine,firstValidPixel,lastValidPixel] = getValidSample( xdoc,subSwath )
item = xdoc.getElementsByTagName('firstValidSample');
%item1 = xdoc.getElementsByTagName('lastValidSample');
numOfGeoLocationGridPoints=getLength(item);
for i=1 : numOfGeoLocationGridPoints
    firstsample = xdoc.getElementsByTagName('firstValidSample').item(i-1).getFirstChild.getData;
    lastsample = xdoc.getElementsByTagName('lastValidSample').item(i-1).getFirstChild.getData;
    S = regexp(char(firstsample), '\s+', 'split');
    L = regexp(char(lastsample), '\s+', 'split');
    n=length(S);
    firstValidSample(i,1:n)=str2double(S);
    lastValidSample(i,1:n)=str2double(L);
end
firstValidPixel = 0;
lastValidPixel = subSwath.numOfSamples;
[~,n]=size(firstValidSample);
for k=0:subSwath.numOfBursts-1
firstValidLineIdx = -1;
lastValidLineIdx = -1;
                for lineIdx = 0:n-1
                    if (firstValidSample(k+1,lineIdx+1) ~= -1) 
                        if (firstValidSample(k+1,lineIdx+1) > firstValidPixel) 
                            firstValidPixel = firstValidSample(k+1,lineIdx+1);
                        end

                        if (firstValidLineIdx == -1) 
                            firstValidLineIdx = lineIdx;
                            lastValidLineIdx = lineIdx;
                         else 
                            lastValidLineIdx=lastValidLineIdx+1;
                        end
                    end
                end
                for lineIdx = 0:n-1
                    if (lastValidSample(k+1,lineIdx+1) ~= -1 &&...
                            lastValidSample(k+1,lineIdx+1) < lastValidPixel) 
                        lastValidPixel = lastValidSample(k+1,lineIdx+1);
                    end
                end
                firstValidLine(k+1)= firstValidLineIdx;
                lastValidLine(k+1) = lastValidLineIdx;

end
end
function data = readSent1Data( dataFile, varargin )
%Extract SLC From Measument/TIFF File
%Zhangfeng Ma, Hohai University, jspcmazhangfeng@hhu.edu.cn
%2019/07/02
narginchk( 1, 3 );
[fileInfo, fid] = getTIFFinfo( dataFile );
roi = struct(); 
verbose = false; 
for ni = 1:nargin-1
    if isstruct( varargin{ni} )
        roi = varargin{ni};
    elseif islogical( varargin{ni} )
        verbose = varargin{ni};
    else
        error( ['Optional input #%d has an invalid type. Please refer to the ' ...
            'help text of this function for a description of inputs.'], ni );
    end
end
Nrg = fileInfo.numPixels;
Naz = fileInfo.numLines;
roiFields = {'firstLine', 1;
             'numLine', Naz;
             'firstSample', 1;
             'numSample', Nrg};
fieldFlag = isfield( roi, roiFields(:,1) ); 
missingFields = roiFields(~fieldFlag,:);
for cnt = 1:size( missingFields, 1 )
    roi.(missingFields{cnt,1}) = missingFields{cnt,2};
end
clear fieldFlag missingFields roiFields;
if roi.firstLine <= 0 || roi.firstLine > Naz
    warning( ['The provided index for the first azimuth line in the ROI is ' ...
        'invalid: %d. Using default (first sample).'], roi.firstLine );
    roi.firstLine = 1;
end
if roi.numLine < 1 || roi.firstLine - 1 + roi.numLine > Naz
    NazRoi = Naz - roi.firstLine + 1; 
    warning( ['The provided number of azimuth lines in the ROI is invalid: ' ...
        '%d. Using %d lines: %d (total number of azimuth lines) - %d (first ' ...
        'line index) + 1.'], roi.numLine, NazRoi, Naz, roi.firstLine );
    roi.numLine = NazRoi; clear NazRoi;
end
if roi.firstSample <= 0 || roi.firstSample > Nrg
    warning( ['The provided index for the first range sample in the ROI is ' ...
        'invalid: %d. Using default (first sample).'], roi.firstSample );
    roi.firstSample = 1;
end
if roi.numSample < 1 || roi.firstSample - 1 + roi.numSample > Nrg
    NrgRoi = Nrg - roi.firstSample + 1; 
    warning( ['The provided number of range samples in the ROI is invalid: ' ...
        '%d. Using %d samples: %d (total number of range samples) - %d ' ...
        '(first sample index) + 1.'], roi.numSample, NrgRoi, Nrg, ...
        roi.firstSample );
    roi.numSample = NrgRoi; clear NrgRoi;
end
clear Nrg Naz;
if verbose
    fprintf( 'Reading data from the file' );
end

roiLines = roi.firstLine + (0:roi.numLine-1);
roiBytesOffset = fileInfo.linesBytesOffset(roiLines) + ...
    (roi.firstSample - 1) * fileInfo.bytesPerSample * (1 + ...
    fileInfo.complexFlag);
clear roiLines;

dataDim = [roi.numSample * (1 + fileInfo.complexFlag), roi.numLine];
NbytesLine = roi.numSample * (1 + fileInfo.complexFlag) * ...
    fileInfo.bytesPerSample; % number of bytes per line
NbytesSkip = roiBytesOffset(2:end) - (roiBytesOffset(1:end-1) + NbytesLine);
clear NbytesLine;

if numel( unique( NbytesSkip ) ) == 1
    NbytesSkip = unique( NbytesSkip );
    if verbose
        fprintf( ' (all data read at once)... ' );
    end
    try
        fseek( fid, roiBytesOffset(1), -1 );
    catch
        error( ['Error setting the pointer at the right position in data ' ...
            'file: %s.'], dataFile );
    end
    if NbytesSkip > 0
        dtype = sprintf( '%d*%s', roi.numSample * (1 + ...
            fileInfo.complexFlag), fileInfo.dataType );
    else
        dtype = fileInfo.dataType;
    end
    data = single( fread( fid, dataDim, dtype, NbytesSkip ) );
    clear dtype;
else
    if verbose
        fprintf( ':\n' );
    end
    data = zeros( dataDim, fileInfo.dataType );
    for nl = 1:roi.numLine
        if verbose
            perc = round( (nl - 1) / roi.numLine * 100 );
            fprintf( '  %d%% complete\r', perc );
            clear perc;
        end
        try
            fseek( fid, roiBytesOffset(nl), -1 );
        catch
            error( ['Error setting the pointer at the right position in data ' ...
                'file: %s.'], dataFile );
        end
        data(:,nl) = fread( fid, roi.numSample * (1 + fileInfo.complexFlag), ...
            fileInfo.dataType );
    end
    data = single( data );
    if verbose
        fprintf( '\n' );
    end
end
clear NbytesSkip roiBytesOffset;
fclose( fid );
if fileInfo.complexFlag
    data = data(1:2:end,:) + 1j * data(2:2:end,:);
end
clear fileInfo;
if verbose
    fprintf( 'Done!\n' );
end
end
function [tiffInfo, varargout] = getTIFFinfo( inpFile )

%Extract Tiff Info From Measument/TIFF File
%Zhangfeng Ma, Hohai University, jspcmazhangfeng@hhu.edu.cn
%2019/07/01
fid = fopen( inpFile, 'r', 'ieee-le' );
assert( fid >= 0, 'Failed to open the input file for reading: %s', inpFile );

order = fread( fid, 2, 'uint8=>char' ); 
if strcmp( order', 'MM' )

    fclose( fid ); 
    fid = fopen( inpFile, 'r', 'ieee-be' ); 
    byteOrder = 'be'; 
elseif strcmp( order', 'II' )
    byteOrder = 'le'; 
else
    fclose( fid );
    error( 'Invalid byte ordering in the input file.' );
end
[~, ~, endian] = computer;
boFlag = ( (endian == 'L' && order(1) == 'M') || ...
    (endian == 'B' && order(1) == 'I') );
clear order endian;

fileType = fread( fid, 1, 'int16' ); 
switch fileType
    case 42
        offsetType = 'uint32';
        ifdOffset = fread( fid, 1, 'uint32' ); 
    case 43
        offsetBytes = fread( fid, 1, 'int16' ); 
        assert( offsetBytes == 8, ['The bytes size of offsets is %d, but 8 ' ...
            'is expected for BigTIFF files.'], offsetBytes );
        clear offsetBytes;
        val = fread( fid, 1, 'int16' );
        assert( val == 0, ['The field after the bytes size of offsets should ' ...
            'be 0 in BigTIFF files, but is %d instead.'], val );
        clear val;
        offsetType = 'uint64';
        ifdOffset = fread( fid, 1, 'uint64' ); 
    otherwise
        fclose( fid );
        error( 'The input file is not in TIFF format.' );
end
val = fseek( fid, ifdOffset, -1 );
assert( val >= 0, ['Error positioning the pointer at the image file ' ...
    'directory (IFD) in the TIFF file.'] );
clear val;
if fileType == 42
    numFields = fread( fid, 1, 'uint16' ); 
elseif fileType == 43
    numFields = fread( fid, 1, 'uint64' );
end
ifdfield(numFields).tag = [];
ifdfield(numFields).type = [];
ifdfield(numFields).count = [];
ifdfield(numFields).offset = [];
for cnt = 1:numFields
    ifdField(cnt).tag = fread( fid, 1, 'uint16' );
    ifdField(cnt).type = fread( fid, 1, 'uint16' );
    ifdField(cnt).count = fread( fid, 1, offsetType );
    ifdField(cnt).offset = fread( fid, 1, offsetType );
end
clear numFields;
numPixels = getTIFFfield( fid, ifdField, 256, offsetType, boFlag );
numLines = getTIFFfield( fid, ifdField, 257, offsetType, boFlag );
numBitsSamp = getTIFFfield( fid, ifdField, 258, offsetType, boFlag );
offsetStrips = getTIFFfield( fid, ifdField, 273, offsetType, boFlag );
numSampPixel = getTIFFfield( fid, ifdField, 277, offsetType, boFlag );
numRowsStrips = getTIFFfield( fid, ifdField, 278, offsetType, boFlag );
numBytesSamp = numBitsSamp / 8;
switch numBytesSamp
    case 1
        dataType = 'uint8';
        cplxFlag = 0;
    case 2
        if numSampPixel == 1
            dataType = 'uint16';
            cplxFlag = 0;
        elseif numSampPixel == 2
            dataType = 'int16';
            cplxFlag = 1;
        else
            error( 'Unsupported number of samples per pixel, %d.', ...
                numSampPixel );
        end
    case 4
        numBytesSamp = 2;
        dataType = 'int16';
        cplxFlag = 1;
    case 8
        numBytesSamp = 4;
        dataType = 'single';
        cplxFlag = 1;
    otherwise
        error( 'Unsupported number of bytes per sample, %d.', numBytesSamp );
end
clear numSampPixel numBitsSamp;
offsetRows = (0:numRowsStrips-1)' * numBytesSamp * (1 + cplxFlag) * numPixels;
linesBytesOffset = bsxfun( @plus, offsetRows, offsetStrips(:)' );
linesBytesOffset = linesBytesOffset(:)'; 
clear offsetRows offsetStrips;


% create the info structure
tiffInfo.endian = byteOrder;
tiffInfo.numPixels = single( numPixels );
tiffInfo.numLines = single( numLines );
tiffInfo.dataType = dataType;
tiffInfo.complexFlag = cplxFlag;
tiffInfo.bytesPerSample = numBytesSamp;
tiffInfo.linesBytesOffset = linesBytesOffset;
clear byteOrder numPixels numLines cplxFlag dataType numButesSamp ...
    linesBytesOffset;


% return pointer to open TIFF file or close it
if nargout == 2
    % pointer requested, return it to caller
    varargout{1} = fid;
else
    % pointer not requested, close the file
    fclose( fid );
end

%%% end of 'getTIFFinfo'
end


%%% beginning of internal function 'getTIFFfield'
function value = getTIFFfield( fid, ifdData, tag, offsetType, boFlag )
% Get the value or array of values for the TIFF field with the specified
% tag, using the information in the image file directory (IFD).
%
% Inputs:
%   - fid: pointer to the open TIFF file.
%   - ifdData: structure containing information about the IFD.
%   - tag: desired field tag (decimal).
%   - offsetType: data type for offsets.
%   - boFlag: flag to change the data byte ordering.
%
% Outputs:
%   - value: extracted value(s) for the TIFF field.
%
% Required functions (not part of MATLAB): none
%
% Author: Louis-Philippe Rousseau (Université Laval)
% Created: May 2014
% Updated: June 2014, October 2017, Novermber 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% find the requested tag in the IFD data structure
idx = find( [ifdData.tag] == tag ); % index of the tag
assert( ~isempty( idx ), 'Tag %s is not present in the IFD.', num2str( tag ) );

% check the number of values (count) for the field
if ifdData(idx).count == 1
    % only one value for the field, the offset is the value
    value = ifdData(idx).offset;

    % set the byte size of value (4 or 8)
    if strcmp( offsetType, 'uint32' )
        % 4 bytes offsets (TIFF)
        bsize = 4;
    elseif strcmp( offsetType, 'uint64' )
        % 8 bytes offsets (BigTIFF)
        bsize = 8;
    else
        % invalid offset type, return error
        error( 'Invalid bytes size for offsets, %s.', offsetType );
    end

    % shift bits, if the flag for byte order change is set
    if boFlag
        switch ifdData(idx).type
            case {1, 2}
                % data type is uint8 or ascii, shift 3 or 7 bytes to the left
                value = bitshift( value, -(bsize-1)*8 );
            case 3
                % data type is uint16, shift 2 or 6 bytes to the left
                value = bitshift( value, -(bsize-2)*8 );
            case 4
                % data type is uint32, shift 0 or 4 bytes to the left
                value = bitshift( value, -(bsize-4)*8 );
            otherwise
                % data type is uint64 (BigTIFF only), do not shift
        end
    end
    clear bsize;
else
    % array of values for the field

    % use offset to move to the right position in the file
    val = fseek( fid, ifdData(idx).offset, -1 );
    assert( val >= 0, ['Error setting the pointer at the right position in ' ...
        'the TIFF file.'] );
    clear val;

    % read the values according to count and data type
    switch ifdData(idx).type
        case 1
            % data type is uint8
            value = fread( fid, ifdData(idx).count, 'uint8' );
        case 2
            % data type is ascii
            value = fread( fid, ifdData(idx).count, 'int8' );
        case 3
            % data type is uint16
            value = fread( fid, ifdData(idx).count, 'uint16' );
        case 4
            % data type is uint32
            value = fread( fid, ifdData(idx).count, 'uint32' );
        case 16
            % data type is uint64 (BigTIFF only)
            value = fread( fid, ifdData(idx).count, 'uint64' );
        otherwise
            % data type is not supported, return error
            error( 'Unsupported data type %d', ifdData(idx).type );
    end
end
clear idx;
end
function [HeadFile,OrbitFile] = ctSent_findSLC

%Find SLC
%Zhangfeng Ma, Hohai University, jspcmazhangfeng@hhu.edu.cn
%2019/07/03
fprintf('--------------- CTSenT  v1.0 ---------------- \n');
fprintf('Zhangfeng Ma, Hohai University 2019/07/10 \n');
fprintf('jspcmazhangfeng@hhu.edu.cn  \n');
fprintf('jspcmazhangfeng@gmail.com  \n');
parm = load('parms');
SLCpath = parm.SLCPATH;
Orbitpath = parm.ORBITPATH;
[dirOutput] = dir(fullfile(SLCpath));
[date,satellite,orbitpath] = findOrbit(Orbitpath);
for i=3:length(dirOutput)
fprintf('----------Reading SLC Directories---------- \n');
HeadFile(i-2) = parseSafeManifest( fullfile(SLCpath,dirOutput(i).name) );
OrbitFile{i-2} = orbitpath{find(strcmp(date,HeadFile(i-2).acquisitionDate)&strcmp(satellite,HeadFile(i-2).Satellite))};
fprintf([HeadFile(i-2).acquisitionDate,' import successfully ! \n']);
fprintf(['Satellite: ',HeadFile(i-2).Satellite, '\n' ]);
fprintf(['Mode: ',HeadFile(i-2).product, '\n' ]);
fprintf(['Track: ',HeadFile(i-2).Ascend_Descend, '\n' ]);
fprintf('Finding Orbit File... \n');
fprintf(['Orbit File for ', HeadFile(i-2).acquisitionDate, ' Found !', '\n']);
pause(1);
end
end 
function [date,satellite,orbitpath]= findOrbit(Orbitpath)
[dirOutput] = dir(fullfile(Orbitpath));
for i=3:length(dirOutput)
    date{i-2} = datestr((datenum(dirOutput(i).name(43:50),'yyyymmdd')+...
        datenum(dirOutput(i).name(59:66),'yyyymmdd'))/2,'yyyymmdd');
    orbitpath{i-2} = fullfile(dirOutput(i).folder,dirOutput(i).name);
    if strcmp(dirOutput(i).name(3),'A')
        satellite{i-2} = 'Sentinel-1A';
    else
        satellite{i-2} = 'Sentinel-1B';
    end
end
end
function maniHead = parseSafeManifest( dsetPath )
% Extract information from the SAFE manifest file of a Sentinel-1 dataset.
%Zhangfeng Ma, Hohai University, jspcmazhangfeng@hhu.edu.cn
%2019/07/02
maniFile = fullfile( dsetPath, 'manifest.safe' ); % manifest file
try
    xdoc = xmlread( maniFile );
catch err
    rethrow( err );
end
xroot = xdoc.getDocumentElement;
if ~strcmp( xroot.getNodeName, 'xfdu:XFDU' )
  
    error( 'The input file is not a SAFE manifest file.' );
end
meta = xroot.getElementsByTagName('metadataSection').item(0);
platform = lower( getXMLnode( 'safe:familyName', meta, 'str' ) );
platform = regexprep( platform, '(\<\w)', '${upper($1)}' );
assert( strcmpi( platform, 'sentinel-1' ), ['The platform in the manifest ' ...
    'file is not supported: %s'], platform );
platform = strcat( platform, getXMLnode( 'safe:number', meta, 'str' ) );
pass = lower( getXMLnode( 's1:pass', meta, 'str' ) );
prodType = getXMLnode( 's1sarl1:productType', meta, 'str' );
prodTest = any( strcmpi( prodType, {'SLC', 'GRD'} ) ); 
assert( prodTest, ['The product type in the manifest file is not ' ...
    'supported: %s'], prodType );
clear prodTest;
acqMode = getXMLnode( 's1sarl1:mode', meta, 'str' );
node = meta.getElementsByTagName('s1sarl1:swath');
numSwath = node.getLength; 
swathID = cell( 1, numSwath ); 
for cnt = 1:numSwath
    swathID{cnt} = char( node.item(cnt-1).getTextContent );
end
node = meta.getElementsByTagName('s1sarl1:transmitterReceiverPolarisation');
numPol = node.getLength; 
pol = cell( 1, numPol );

for cnt = 1:numPol
    pol{cnt} = char( node.item(cnt-1).getTextContent );
end
numChan = length( pol ) * length( swathID );
dateStrFmt = 'yyyy-mm-ddTHH:MM:SS.FFFFFF';
acqT0 = getXMLnode( 'safe:startTime', meta, 'dateStr', dateStrFmt );
acqTE = getXMLnode( 'safe:stopTime', meta, 'dateStr', dateStrFmt );
data = xroot.getElementsByTagName('dataObjectSection').item(0);
dataFile = cell( 1, numChan ); cntD = 1;
prodAnn = cell( 1, numChan ); cntP = 1;
calAnn = cell( 1, numChan ); cntC = 1;
noiseAnn = cell( 1, numChan ); cntN = 1;
for cnt = 1:data.getLength
    if data.item(cnt-1).hasAttributes
        repID = char( data.item(cnt-1).getAttribute('repID') );
        loc = data.item(cnt-1).getElementsByTagName('fileLocation').item(0);
        if strcmp( repID, 's1Level1MeasurementSchema' )
            str = char( loc.getAttribute('href') );
            dataFile{cntD} = str(3:end);
            cntD = cntD + 1;
        elseif strcmp( repID, 's1Level1ProductSchema' )
            str = char( loc.getAttribute('href') );
            prodAnn{cntP} = str(3:end);
            cntP = cntP + 1;
        elseif strcmp( repID, 's1Level1CalibrationSchema' )
            str = char( loc.getAttribute('href') );
            calAnn{cntC} = str(3:end);
            cntC = cntC + 1;
        elseif strcmp( repID, 's1Level1NoiseSchema' )
            str = char( loc.getAttribute('href') );
            noiseAnn{cntN} = str(3:end);
            cntN = cntN + 1;
        end
    end
end
maniHead.fileDir = fileparts( maniFile );
maniHead.Satellite = platform; 
maniHead.Ascend_Descend = pass; 
maniHead.product = prodType; 
maniHead.acquisition = acqMode; 
maniHead.swath = swathID; 
maniHead.polarization = pol; 
maniHead.StartTimeUTC = acqT0;
maniHead.StopTimeUTC = acqTE;
maniHead.slcFile = dataFile; 
maniHead.headFile = prodAnn; 
maniHead.calFile = calAnn; 
maniHead.noiseFile = noiseAnn; 
maniHead.filepart = applyFilePart(maniHead);
maniHead.acquisitionDate = datestr(maniHead.StartTimeUTC,'yyyymmdd');
end
function filepart = applyFilePart(maniHead)
slcFile = maniHead.slcFile;
headFile = maniHead.headFile;
for i = 1:length(slcFile)
    polari = slcFile{i}(25:26);
    iw = str2double(slcFile{i}(19));
     if strcmp(polari,'vv') & iw == 1
         filepart.IW1VV.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     elseif strcmp(polari,'vv') & iw == 2
         filepart.IW2VV.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     elseif strcmp(polari,'vv') & iw == 3
         filepart.IW3VV.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     elseif strcmp(polari,'vh') & iw == 1
         filepart.IW1VH.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     elseif strcmp(polari,'vh') & iw == 2
         filepart.IW2VH.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     elseif strcmp(polari,'vh') & iw == 3
         filepart.IW3VH.SLC = fullfile(maniHead.fileDir,maniHead.slcFile{i});
     end
end
for i = 1:length(headFile)
    polari = headFile{i}(24:25);
    iw = str2double(headFile{i}(18));
     if strcmp(polari,'vv') & iw == 1
         filepart.IW1VV.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     elseif strcmp(polari,'vv') & iw == 2
         filepart.IW2VV.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     elseif strcmp(polari,'vv') & iw == 3
         filepart.IW3VV.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     elseif strcmp(polari,'vh') & iw == 1
         filepart.IW1VH.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     elseif strcmp(polari,'vh') & iw == 2
         filepart.IW2VH.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     elseif strcmp(polari,'vh') & iw == 3
         filepart.IW3VH.ANNO = fullfile(maniHead.fileDir,maniHead.headFile{i});
     end
end
end
function data = getXMLnode( tag, pnode, varargin )
narginchk( 2, 5 );
cnt = 0; % extract first tag occurence
type = 'str'; % extract the tag data as a string
dateFmt = ''; % empty date format string


dtypeList = {'str', 'dbl', 'dblArr', 'int', 'intArr', 'uint', 'uintArr', ...
    'dateStr', 'dateVec'};


%dateStrList = {'yyyy', 'mm', 'dd', 'HH', 'MM', 'SS', 'FFF', 'FFFFFF'};


dateVecList = {'year', 'mon', 'day', 'hour', 'min', 'sec', 'msec', 'usec'};


if nargin > 2
    % optional inputs provided, assign each value to the correct variable
    for idx = 1:nargin-2
        if isnumeric( varargin{idx} )
            % current input is the occurence number for the tag
            cnt = varargin{idx};
        elseif any( strcmpi( varargin{idx}, dtypeList ) )
            % current input is the data type for the tag
                type = varargin{idx};
        else
            % current input is the date format string, or is invalid
            dateFmt = varargin{idx};
        end
    end
end

try
    % success
    ndata = pnode.getElementsByTagName(tag).item(cnt).getTextContent;
catch
    % fail, return error
    error( 'Non-existent node with name "%s"', tag );
end


if strcmpi( type, 'str' )
    % string
    data = char( ndata );
elseif strncmpi( type, 'date', 4 )
    % date vector or string

    % convert tag data according to the provided type
    if strcmpi( type, 'dateStr' )
        % date string type, convert to string
        ndata = char( ndata );
    else
        % date vector type, convert data vector to doubles
        ndata = str2num( ndata );
    end

    if isempty( dateFmt )
        % date format string not provided, return converted data directly
        data = ndata;
    else
        % date format string provided, use it to correctly convert data
        if strcmpi( type, 'dateStr' )
            % date string type, process string

            % find sub-second portion of the date string, if any
            subSecIdx = strfind( dateFmt, 'F' ); % sub-second part of the string

            % convert string to date vector
            if isempty( subSecIdx )
                % no sub-second part in date string, use it directly
                data = datevec( ndata, dateFmt );
            else
                % sub-second part present in date string, handle it separately
                data = datevec( ndata(1:subSecIdx(1)-1), ...
                    dateFmt(1:subSecIdx(1)-1) );
                data(6) = data(6) + str2double( ndata(subSecIdx) ) / ...
                    10^(length(subSecIdx));
            end
        else
            % date vector type, process data vector

            % split the format string at the commas
            dFields = strsplit( dateFmt, ',' );

            % construct the date vector by looping over the elements
            data = zeros( 1, 6 ); % initialize the date vector
            for nel = 1:6
                % search the current element in the format string
                elIdx = strcmp( dFields, dateVecList{nel} );

                % if it exists, copy the current element in the date vector
                if ~isempty( elIdx )
                    if nel < 6
                        % for the first five elements, round the data
                        data(nel) = round( ndata(elIdx) );
                    else
                        % for the last element, also check for sub-seconds data
                        msecIdx = strcmp( dFields, 'msec' ); % milliseconds
                        usecIdx = strcmp( dFields, 'usec' ); % microseconds
                        if ~isempty( usecIdx )
                            % microsecond element present in date vector
                            data(6) = round( ndata(elIdx) ) + ...
                                ndata(usecIdx) / 1e6;
                        elseif ~isempty( msecIdx )
                            % millisecond element present in date vector
                            data(6) = round( ndata(elIdx) ) + ...
                                ndata(msecIdx) / 1e3;
                        else
                            % no sub-second component in date vector
                            data(6) = ndata(elIdx);
                        end
                        clear msecIdx usecIdx;
                    end
                end
                clear elIdx;
            end
        end
    end
else
    % numeric value or array

    % check if an array type is specified
    if ~strcmpi( type(end-2:end), 'Arr' )
        % single value, convert to double for now
        tmpDbl = str2double( ndata );
    else
        % array, convert to an array of doubles for now
        tmpDbl = str2num( ndata );
    end

    % convert double value(s) to another numeric type, if necessary
    if strncmpi( type, 'uint', 4 )
        % unsigned integer(s)
        data = uint64( tmpDbl );
    elseif strncmpi( type, 'int', 3 )
        % signed integer(s)
        data = int64( tmpDbl );
    else
        % double(s)
        data = tmpDbl;
    end
    clear tmpDbl;
end
clear ndata;
end
function [Lat,Lon] = initial_radar_code(info_m,mburstIndex,dem_read)
%CTSENT_RADARCODE Summary of this function goes here
%   Detailed explanation goes here
parm = load('parms');
orbit_Degree = parm.orbit_degree;
morbit = orbit(info_m,orbit_Degree);
[DEM_ex,Ellipse] = computeDEMbox(info_m,mburstIndex,dem_read);
[time_m] = xyz2t(Ellipse,info_m,morbit);
[mline,mpixel] = ta2lp( info_m,time_m,mburstIndex);
ratio = info_m.sensor_para.rangePixelSpacing/info_m.sensor_para.azimuthPixelSpacing;
[mline_grid,mpixel_grid] = gridlinepixel(info_m);
Lat = trinterp(mline,mpixel,DEM_ex.lat,mline_grid,mpixel_grid,ratio);
Lon = trinterp(mline,mpixel,DEM_ex.lon,mline_grid,mpixel_grid,ratio);
end
function Inc = Theta(M, P)
[m,n] = size(M.X);
Master=[M.X(:),M.Y(:),M.Z(:)];
Point=[P.X(:),P.Y(:),P.Z(:)];
r1 = Master - Point;
Inc = vec_angle(Point, r1);
Inc = reshape(Inc,[m,n]);
end
function normvec = normalvec(PointVec)
normvec = sqrt(PointVec(:,1).^2+PointVec(:,2).^2+PointVec(:,3).^2);
end
function angle_of_vec = vec_angle(vec1, vec2)
   in_vec = sum(vec1 .* vec2,2);
   angle_of_vec = acos(in_vec ./ (normalvec(vec1).*normalvec(vec2)));
end
function [vq] = trinterp(x,y,v,xq,yq,ratio)
[m,n] = size(xq);
X = x(:);Y = y(:)*ratio;Z = v(:);
xq = xq(:);yq = yq(:);
TRI = delaunayTriangulation(X,Y);
ID = pointLocation(TRI,xq,yq.*ratio);
Vx0 = X(TRI(ID,1));Vx1 = X(TRI(ID,2));Vx2 = X(TRI(ID,3));
Vy0 = Y(TRI(ID,1))./ratio;Vy1 = Y(TRI(ID,2))./ratio;Vy2 = Y(TRI(ID,3))./ratio;
Vz0 = Z(TRI(ID,1));Vz1 = Z(TRI(ID,2));Vz2 = Z(TRI(ID,3));
xkj = Vx1 - Vx0;
ykj = Vy1 - Vy0;
xlj = Vx2 - Vx0;
ylj = Vy2 - Vy0;
f = 1 ./ (xkj.*ylj - ykj.*xlj);
zj = Vz0;zk = Vz1;zl = Vz2;
zkj = zk - zj;
zlj = zl - zj;
a = -f.*(ykj.*zlj - zkj.*ylj);
b = -f.*(zkj.*xlj - xkj.*zlj);
c = -a.*Vx1 - b.*Vy1 + zk;
zout = a.* xq + b.* yq + c;
vq = reshape(zout,[m,n]);
% [ti,bc] = pointLocation(TRI,xq(:),yq(:));
% bc(isnan(ti),:)=[];ti(isnan(ti))=[];
% triVals = v(TRI(ti,:));
% vq = dot(bc',triVals')';
% vq = reshape(vq,size(xq));
end
function [ line,pixel ] = ta2lp( info_m,t,burstIndex)
%PRF=1.717128973878037e+03;		
line=1+(t.timeAzimuth-info_m.subSwath.burstFirstLineTime(burstIndex))/info_m.sensor_para.line_time_interval;
pixel=1+(t.timeRange-info_m.sensor_para.t_Range1)*299792458/info_m.subSwath.rangePixelSpacing;
% pixel=1+(t.timeRange-info_m.sensor_para.t_Range1)*info_m.sensor_para.RSR2x;
end
function [ dopplertime,rangetime] = getZeroDopplerTime( info ,earthPoint)
  [numOrbitVec,~] = size(info.orbit_para.orbit_points);
  lineTimeInterval=info.sensor_para.line_time_interval;
  wavelength=info.sensor_para.Wavelength;
  firstVecTime = 0.0;
  secondVecTime = 0.0;
  firstVecFreq = 0.0;
  secondVecFreq = 0.0;
        for i = 0:numOrbitVec-1
            sensorPosition.x=info.orbit_para.orbit_points(i+1,2);
            sensorPosition.y=info.orbit_para.orbit_points(i+1,3);
            sensorPosition.z=info.orbit_para.orbit_points(i+1,4);
            sensorVelocity.x=info.orbit_para.orbit_points(i+1,5);
            sensorVelocity.y=info.orbit_para.orbit_points(i+1,6);
            sensorVelocity.z=info.orbit_para.orbit_points(i+1,7);
            currentFreq = getDopplerFrequency(earthPoint, sensorPosition, sensorVelocity, wavelength);
            if (i == 0 || firstVecFreq * currentFreq > 0) 
                firstVecTime = info.orbit_para.orbit_points(i+1,1);
                firstVecFreq = currentFreq;
             else 
                secondVecTime = info.orbit_para.orbit_points(i+1,1);
                secondVecFreq = currentFreq;
                break;
            end
        end

        if (firstVecFreq * secondVecFreq >= 0.0) 
            disp( 'NonValidZeroDopplerTime');
        end
        lowerBoundTime = firstVecTime;
        upperBoundTime = secondVecTime;
        lowerBoundFreq = firstVecFreq;
        upperBoundFreq = secondVecFreq;
        diffTime = abs(upperBoundTime - lowerBoundTime);
        absLineTimeInterval = abs(lineTimeInterval);%abs(lineTimeInterval)/86400;
        totalIterations = floor(diffTime/ absLineTimeInterval) + 1;
        numIterations = 0;
        while (diffTime > absLineTimeInterval && numIterations <= totalIterations) 

            midTime = (upperBoundTime + lowerBoundTime) / 2.0;
            [sensorPosition, sensorVelocity] = getPositionVelocity(midTime,info);
            midFreq = getDopplerFrequency(earthPoint, sensorPosition, sensorVelocity, wavelength);

            if (midFreq * lowerBoundFreq > 0.0) 
                lowerBoundTime = midTime;
                lowerBoundFreq = midFreq;
            elseif (midFreq * upperBoundFreq > 0.0) 
                upperBoundTime = midTime;
                upperBoundFreq = midFreq;
            elseif (midFreq == 0) 
                dopplertime= midTime;
                break;
            end

            diffTime = abs(upperBoundTime - lowerBoundTime);
            numIterations=numIterations+1;
        end


        dopplertime= lowerBoundTime - lowerBoundFreq * (upperBoundTime - lowerBoundTime) / (upperBoundFreq - lowerBoundFreq);
        rangetime = computeSlantRange(earthPoint,dopplertime,info);
end
function Dop=getDopplerFrequency(earthPoint, sensorPosition, sensorVelocity, wavelength)
xDiff = earthPoint.X - sensorPosition.x;
yDiff = earthPoint.Y - sensorPosition.y;
zDiff = earthPoint.Z - sensorPosition.z;
distance = sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
Dop=2.0 * (sensorVelocity.x * xDiff + sensorVelocity.y * yDiff + sensorVelocity.z * zDiff) / (distance * wavelength);
end
function rangetime = computeSlantRange(earthPoint,dopplerTime,info)
[sensorPos, ~]=getPositionVelocity(dopplerTime,info);
xDiff = sensorPos.x - earthPoint.X;
yDiff = sensorPos.y - earthPoint.Y;
zDiff = sensorPos.z - earthPoint.Z;
slantRange=sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
rangetime = slantRange/299792458;
end
function [line,pixel] = determinZeroDoppler(info,Ellipse,burstIndex)
[m,n] = size(Ellipse.X);
line = zeros(m,n);
pixel = line;
for i = 1:m
    for j = 1:n
        earthPoint.X = Ellipse.X(i,j);earthPoint.Y = Ellipse.Y(i,j);earthPoint.Z = Ellipse.Z(i,j);
        [ dopplertime,rangetime ] = getZeroDopplerTime( info ,earthPoint);
        line(i,j) = (dopplertime - info.subSwath.burstFirstLineTime(burstIndex)) / info.subSwath.azimuthTimeInterval+1;
        pixel(i,j) = 1+(rangetime-info.sensor_para.t_Range1)*299792458/info.subSwath.rangePixelSpacing;
    end
end
end
function [line,pixel] = gridlinepixel(info)
azimuthline = (1:info.sensor_para.linesPerBurst)';
rangepixel = 1:info.sensor_para.samplesPerBurst;
line = repmat(azimuthline,1,length(rangepixel));
pixel = repmat(rangepixel,length(azimuthline),1);
end
function [DEM_ex,Ellipse] = computeDEMbox(info,burstIndex,dem_read)
longitude = info.subSwath.longitude(burstIndex:burstIndex+1,:);
lonMin = min(longitude(:));lonMax = max(longitude(:));
latitude = info.subSwath.latitude(burstIndex:burstIndex+1,:);
latMin = min(latitude(:));latMax = max(latitude(:));
LatitudeLimitLo = min(dem_read.lat(:));LatitudeLimitHi = max(dem_read.lat(:));
LongitudeLimitLo = min(dem_read.lon(:));LongitudeLimitHi = max(dem_read.lon(:));
deltaDegree = abs(dem_read.lat(1,1)-dem_read.lat(2,1));
latMin = latMin - 200*deltaDegree;latMax = latMax + 200*deltaDegree;
lonMin = lonMin - 200*deltaDegree;lonMax = lonMax + 200*deltaDegree;
[RasterLine,RasterWidth] = size(dem_read.lat);
upperLeft.y = (latMin-LatitudeLimitLo)*RasterLine/(LatitudeLimitHi-LatitudeLimitLo)+1;
upperLeft.x = (lonMin-LongitudeLimitLo)*RasterWidth/(LongitudeLimitHi-LongitudeLimitLo)+1;
lowerRight.y = (latMax-LatitudeLimitLo)*RasterLine/(LatitudeLimitHi-LatitudeLimitLo)+1;
lowerRight.x = (lonMax-LongitudeLimitLo)*RasterWidth/(LongitudeLimitHi-LongitudeLimitLo)+1;
latMinIdx = floor(upperLeft.y);
latMaxIdx = floor(lowerRight.y);
lonMinIdx = floor(upperLeft.x);
lonMaxIdx = floor(lowerRight.x);
if latMinIdx <= 0
    latMinIdx = 1;
    fprintf('Warning: DEM cannot cover the whole data region !!! \n');
end
if lonMinIdx <= 0
    lonMinIdx = 1;
    fprintf('Warning: DEM cannot cover the whole data region !!! \n');
end
if latMaxIdx > RasterLine
    latMaxIdx = RasterLine;
    fprintf('Warning: DEM cannot cover the whole data region !!! \n');
end
if lonMaxIdx > RasterWidth
    lonMaxIdx = RasterWidth;
    fprintf('Warning: DEM cannot cover the whole data region !!! \n');
end
GeoBoxlat = floor(latMin):0.1:ceil(latMax);
GeoBoxlon = floor(lonMin):0.1:ceil(lonMax);
Lat_n = length(GeoBoxlat);
Lon_n = length(GeoBoxlon);
lat =repmat(GeoBoxlat',1,Lon_n);
lon =repmat(GeoBoxlon,Lat_n,1);
lon(lon<0) = abs(lon(lon<0))+180;
heightCorrection = geoidheight(lat,lon);
heightCorrection = reshape(heightCorrection,[Lat_n Lon_n]);
lon =repmat(GeoBoxlon,Lat_n,1);
DEM_ex.lat = dem_read.lat(latMinIdx:latMaxIdx,lonMinIdx:lonMaxIdx);
DEM_ex.lon = dem_read.lon(latMinIdx:latMaxIdx,lonMinIdx:lonMaxIdx);
DEM_ex.z = dem_read.z(latMinIdx:latMaxIdx,lonMinIdx:lonMaxIdx);
EGM96correction = griddata(lat,lon,heightCorrection,DEM_ex.lat,DEM_ex.lon);
DEM_ex.z = double(DEM_ex.z)+EGM96correction;
% [m,n] = size(DEM_ex.z);
% [pixel,line] = meshgrid(1:0.1:n,1:0.1:m);
% z = ba_interp2(DEM_ex.z,pixel,line,'linear');
% lat = ba_interp2(DEM_ex.lat,pixel,line,'linear');
% lon = ba_interp2(DEM_ex.lon,pixel,line,'linear');
% DEM_ex.lat = lat;DEM_ex.lon = lon;DEM_ex.z = z;
Ellipse = BLH2XYZ(DEM_ex.lat,DEM_ex.lon,DEM_ex.z);
end
function phase = computeDerampDemodPhase(info,burstIndex)
[ dopplerRate] = computeDopplerRate( info,burstIndex );
RangeDependDopplerRate=computeRangeDependentDopplerRate(info.subSwath,burstIndex);
dopplerCentroid  = computeDopplerCentroid( info.subSwath,burstIndex);
referenceTime= computeReferenceTime( dopplerCentroid,RangeDependDopplerRate,info.subSwath );
temp =  ((0:(info.subSwath.linesPerBurst-1)).*info.subSwath.azimuthTimeInterval)';
ta = repmat(temp,1,info.sensor_para.samplesPerBurst);
deramp = -pi .* dopplerRate .* (ta - referenceTime).^2;
demod = -2 .* pi .* dopplerCentroid .* ta;
phase = deramp + demod;
end
function derampDemodSLC = performDerampDemod(slcData, derampDemodPhase)
phaseCplx = complex(cos(derampDemodPhase),sin(derampDemodPhase));
derampDemodSLC = slcData.*conj(phaseCplx);
end
function [ dopplerRate] = computeDopplerRate( info,burstIndex )
%morbit=orbit(info,3);
RangeDependentDopplerRate=computeRangeDependentDopplerRate(info.subSwath,burstIndex);
azTime = (info.sensor_para.first_line_time + info.sensor_para.last_line_time)/2.0;
[~,velocity] = getPositionVelocity(azTime,info);
v = getVelocity(velocity); 
steeringRate = info.subSwath.azimuthSteeringRate * (pi / 180.0);
krot = 2*v*steeringRate/info.sensor_para.Wavelength;
dopplerRate= RangeDependentDopplerRate * krot...
                            ./ (RangeDependentDopplerRate - krot);
end
function RangeDependentDopplerRate = computeRangeDependentDopplerRate(subSwath,burstIndex)
azFmRateList = subSwath.azFmRateList;
x = 1:subSwath.samplesPerBurst;
slrt = getSlantRangeTime(x,subSwath)*2;
dt = slrt - azFmRateList(burstIndex).t0;
RangeDependentDopplerRate =...
   azFmRateList(burstIndex).c0 + azFmRateList(burstIndex).c1.*dt + azFmRateList(burstIndex).c2.*dt.^2;
end
function [ dopplerCentroid ] = computeDopplerCentroid( subSwath,burstIndex )
dcBurstList  = subSwath.DCEstimate;
x = 1:subSwath.samplesPerBurst;
slrt = getSlantRangeTime(x,subSwath)*2;
dt = slrt - dcBurstList(burstIndex).t0;
dopplerCentroid = dcBurstList(burstIndex).c0 + dcBurstList(burstIndex).c1.*dt + dcBurstList(burstIndex).c2*dt.^2;
end
function [ referenceTime ] = computeReferenceTime( dopplerCentroid,rangeDependDopplerRate,subSwath )

tmp1 = subSwath.linesPerBurst * subSwath.azimuthTimeInterval / 2.0;
tmp2 = tmp1 + dopplerCentroid(subSwath.firstValidPixel+1) /...
       rangeDependDopplerRate(subSwath.firstValidPixel+1);
referenceTime = tmp2 - dopplerCentroid ./ rangeDependDopplerRate;

end
function slrt=getSlantRangeTime(x,subSwath)
 slrt=subSwath.slrTimeToFirstPixel +(x-1) * subSwath.rangePixelSpacing / 299792458;
end
function [Velocity]=getVelocity(POS)
Velocity=sqrt(POS.x*POS.x+POS.y*POS.y+POS.z*POS.z);
end
function [position,velocity] = getPositionVelocity(time,info)
nv = 8;
orbitStateVectors=info.orbit_para.orbit_points;
[m,~]=size(orbitStateVectors);
dt = (orbitStateVectors(m,1)-orbitStateVectors(1,1)) / (m-1);
            if (m <= nv) 
                i0 = 0;
                iN = m - 1;
            else 
                i0 = max(floor(((time - orbitStateVectors(1,1) )/ dt)) - nv/2 + 1, 0);
                iN = min(i0 + nv - 1, m - 1);
                if iN < m - 1
                i0 = i0;
                else
                i0=iN - nv + 1;
                end
            end
                position.x = 0;
                position.y = 0;
                position.z = 0;
                velocity.x = 0;
                velocity.y = 0;
                velocity.z = 0;
            for i = i0:iN
                orbI.time_mjd = orbitStateVectors(i+1,1);
                orbI.x_pos = orbitStateVectors(i+1,2);
                orbI.y_pos = orbitStateVectors(i+1,3);
                orbI.z_pos = orbitStateVectors(i+1,4);
                orbI.x_vel = orbitStateVectors(i+1,5);
                orbI.y_vel = orbitStateVectors(i+1,6);
                orbI.z_vel = orbitStateVectors(i+1,7);
                weight = 1;
                for j = i0:iN
                    if (j ~= i) 
                        time2 = orbitStateVectors(j+1,1);
                        weight = weight * (time - time2) / (orbI.time_mjd - time2);
                    end
                end
                    position.x = position.x + weight * orbI.x_pos;
                    position.y = position.y + weight * orbI.y_pos;
                    position.z = position.z + weight * orbI.z_pos;
                    velocity.x = velocity.x + weight * orbI.x_vel;
                    velocity.y = velocity.y + weight * orbI.y_vel;
                    velocity.z = velocity.z + weight * orbI.z_vel;
            end
end
function index=computeCornerBasedIndex(x,y,sourceRectangle)  
x=x + 0.5; y=y + 0.5;
index.x = x;
index.y = y;
index.x0=sourceRectangle(1);
index.y0=sourceRectangle(2);
index.width = sourceRectangle(3);
index.height = sourceRectangle(4);
i0 = floor(x);
j0 = floor(y);
di = x - (i0 + 0.5);
dj = y - (j0 + 0.5);
index.i0 = i0;
index.j0 = j0;
iMax = index.width;
        if (di >= 0) 
            i1 = i0 + 1;
            if (i0 < 0)
            index.i(1) =  0 ;
            elseif(i0 > iMax) 
            index.i(1) =iMax ;
            else index.i(1) =i0;
            end
            if (i1 < 0)
            index.i(2) =  0 ;
            elseif(i1 > iMax) 
            index.i(2) =iMax ;
            else index.i(2) =i1;
            end
            index.ki(1) = di;
         else 
            i1 = i0 - 1;
            if (i1 < 0)
            index.i(1) =  0 ;
            elseif(i1 > iMax) 
            index.i(1) =iMax ;
            else index.i(1) =i1;
            end
            if (i0 < 0)
            index.i(2) =  0 ;
            elseif(i0 > iMax) 
            index.i(2) =iMax ;
            else index.i(2) =i0;
            end
            index.ki(1) = di + 1;
        end

        jMax = index.height;
        if (dj >= 0) 
            j1 = j0 + 1;
            if (j0 < 0)
            index.j(1) =  0 ;
            elseif(j0 > jMax) 
            index.j(1) =jMax ;
            else index.j(1) =j0;
            end
            if (j1 < 0)
            index.j(2) =  0 ;
            elseif(j1 > jMax) 
            index.j(2) =jMax; 
            else index.j(2) =j1;
            end
            index.kj(1) = dj;
         else 
            j1 = j0 - 1;
            if (j0 < 0)
            index.j(1) =  0 ;
            elseif(j1 > jMax) 
            index.j(1) =jMax ;
            else index.j(1) =j1;
            end
            if (j0 < 0)
            index.j(2) =  0 ;
            elseif(j0 > jMax) 
            index.j(2) =jMax ;
            else index.j(2) =j0;
            end
            index.kj(1) = dj + 1;
        end
end
function value=resample(raster,index)
x = [ floor(index.i(1)), floor( index.i(2))];
y = [ floor(index.j(1)), floor( index.j(2))];
samples(1,1) = raster(y(1),x(1));
samples(1,2) = raster(y(1),x(2));
samples(2,1) = raster(y(2),x(1));
samples(2,2) = raster(y(2),x(2));
if (samples(1,1)~=0&&samples(1,2)~=0&&samples(2,1)~=0&&samples(2,2)~=0)
ki = index.ki(1);
kj = index.kj(1);
value= samples(1,1)  * (1 - ki) * (1 - kj) +...
       samples(1,2) * ki * (1 - kj) +...
       samples(2,1)  * (1 - ki) * kj +...
       samples(2,2) * ki * kj;
else
    value=samples(1,1);
end
end
function [time] = xyz2t(pointOnEllips,info,coeff)
if nargin < 5
    CRITERTIM = 10e-11;
end
SOL=299792458;
if nargin < 4
    MAXITER = 50;
end
[m,n]=size(pointOnEllips.X);
timeAzimuth =ones(m,n).*line2ta(info.sensor_para.current_window.bottom/2,info.sensor_para.t_Azi1,info.sensor_para.line_time_interval); 
%solution = zeros(m,n,MAXITER);
for iter = 1:MAXITER
            satellitePosition = getXYZ(timeAzimuth,coeff);
            satelliteVelocity = getXYZDot(timeAzimuth,coeff);
            satelliteAcceleration = getXYZDotDot(timeAzimuth,coeff);
            delta.X = pointOnEllips.X-satellitePosition.X;
            delta.Y = pointOnEllips.Y-satellitePosition.Y;
            delta.Z = pointOnEllips.Z-satellitePosition.Z;
            solution = -eq1_Doppler(satelliteVelocity, delta) ./ eq1_Doppler_dt(delta, satelliteVelocity, satelliteAcceleration);
            timeAzimuth =timeAzimuth + solution ;
            if max(max(abs(solution))) < CRITERTIM
                break;
            end
            if (iter >= MAXITER) 
            disp('Iteration Numer exceeds 50...');
                break;
            end
end
        
satellitePosition = getXYZ(timeAzimuth,coeff);
delta.X = pointOnEllips.X-satellitePosition.X;
delta.Y = pointOnEllips.Y-satellitePosition.Y;
delta.Z = pointOnEllips.Z-satellitePosition.Z;
timeRange = sqrt(delta.X.^2 +delta.Y.^2+ delta.Z.^2)./ SOL;
time.timeAzimuth=timeAzimuth;
time.timeRange=timeRange;
end
function [POS]=getXYZ(azTime,coeff)
azTimeNormal = (azTime - coeff.time(floor(length(coeff.time) / 2)+1))./10;
% norm=(t-(t(floor(length(t)/2)+1)))./10;
%azTimeNormal = azTime;
POS.X=polyVal1D(azTimeNormal, coeff.X);
POS.Y=polyVal1D(azTimeNormal, coeff.Y);
POS.Z=polyVal1D(azTimeNormal, coeff.Z);
end
function [sum]=polyVal1D(aziTime, coeffs)
sum = 0.0;
d=length(coeffs);
        while d > 0
            sum =sum.* aziTime;
            sum =sum+ coeffs(d);
            d=d-1;
        end
end
function [POS]=getXYZDot(azTime,coeff)
azTimeNormal = (azTime - coeff.time(floor(length(coeff.time) / 2)+1))./10 ;
%azTimeNormal = azTime;
DEGREE = length(coeff.X)-1;
x = coeff.X(2);
y = coeff.Y(2);
z = coeff.Z(2);
        for  i = 2:DEGREE
            powT = i .*power(azTimeNormal, i - 1);
            x = x+coeff.X(i+1) .* powT;
            y = y+coeff.Y(i+1) .* powT;
            z = z+coeff.Z(i+1) .* powT;
        end
POS.X=x/10;
POS.Y=y/10; 
POS.Z=z/10;
end
function  [POS]=getXYZDotDot(azTime,coeff)
poly_degree = length(coeff.X)-1;
azTimeNormal = (azTime - coeff.time(floor(length(coeff.time) / 2))+1)./10 ;
%azTimeNormal = azTime;
x=0; y=0; z=0;
        for  i = 2: poly_degree
            powT = ((i - 1) * i) .* power(azTimeNormal, i - 2);
            x =x+ coeff.X(i+1) .* powT;
            y =y+ coeff.Y(i+1) .* powT;
            z =z+ coeff.Z(i+1) .* powT;
        end
POS.X=x/100;
POS.Y=y/100;
POS.Z=z/100;
end
function [POS]=eq1_Doppler(satVelocity,pointOnEllips)
POS=satVelocity.X.*pointOnEllips.X+satVelocity.Y.*pointOnEllips.Y+satVelocity.Z.*pointOnEllips.Z;
end
function [POS]=eq1_Doppler_dt(pointEllipsSat,satVelocity,satAcceleration)
     POS=  satAcceleration.X.*pointEllipsSat.X+satAcceleration.Y.*pointEllipsSat.Y+satAcceleration.Z.*pointEllipsSat.Z ...
       - satVelocity.X.*satVelocity.X - satVelocity.Y.*satVelocity.Y - satVelocity.Z.*satVelocity.Z;
end
function [Coeff] = orbit(info,degree)
orbitStateVectors = info.orbit_para.orbit_points;
numStateVectors = length(orbitStateVectors);
data.time= orbitStateVectors(:,1);
%data.Time=orbitStateVectors(11:37,1);
data.X = orbitStateVectors(:,2);
data.Y = orbitStateVectors(:,3);
data.Z= orbitStateVectors(:,4);
poly_degree = degree;
Coeff=computeCoefficients(data,poly_degree);
end
function norm=normalize( t) 
%t.sub(t.get(t.length / 2)).div(10.0);
norm=(t-(t(floor(length(t)/2)+1)))./10;
end
function coeff=computeCoefficients(data,poly_degree)
coeff.X = polyFitNormalized(data.time, data.X, poly_degree);
coeff.Y = polyFitNormalized(data.time, data.Y, poly_degree);
coeff.Z = polyFitNormalized(data.time, data.Z, poly_degree);
coeff.time=data.time;
end
function coeff=polyFitNormalized(time, data, poly_degree)
time=normalize(time);
numOfPoints=length(time);
numOfUnkowns=poly_degree+1;
if(numOfPoints<numOfUnkowns)
    disp('para too cheap');
end
%Amat=[time.^3, time.^2, time ,ones(numOfPoints,1)];
Amat=[ones(numOfPoints,1),  time, time.^2 ,time.^3];
N=Amat'*Amat;
rhs=Amat'*data;
coeff=SymmetricBandedCholesky(N,rhs); 
%coeff=N\rhs;
end
function [x]=SymmetricBandedCholesky(A,q)

[m,n]=size(A);
if (m~=n)
    error('Coefficient Matrix A Must be Square')
end
C=compactstorage(A);
B= compactchol(C);
L=unpack(B);L=L';
y=bandforward(L,q);%solves Ly=b
x=bandback(L',y);%solves L'x=y
end
%q=[5,3.55,2.81428571428571,2.34642857142857,2.01746031746032];
function B = compactstorage(A)

dim=size(A);
if ~(dim(1)==dim(2))
    error('A must be square')
end
if (all((all(A)~=all(A'))))
    error('A must be symmetric')
end
if ~(all(eig(A))> 0)
    error('Matrix is at least not positive definite')
end
c=find(A(1,1:dim(1))~=0);
B=zeros(dim(1),c(end));
n=dim(1);p=c(end)-1;
for i=1:n
if i<=n-p
for j=i:p+i
B(i,j-i+1)=A(i,j);
end
else 
for j=i:n
B(i,j-i+1)=A(i,j);
end
end
end
end
function C = compactchol(B)
[m,p] = size(B);  p = p - 1;
C = B;
for k = 1:m
    last = min(k+p,m) - k + 1;
    for j = 2:last
        i = k + j - 1;
        C(i,1:last-j+1) = C(i,1:last-j+1) - ...
            ((C(k,j))/C(k,1))*C(k,j:last); 
    end
    C(k,:) = C(k,:)/sqrt(C(k,1));
end
C(end-(p-1):end,end) = 0;
end
function L=unpack(B)
[m,n]=size(B);
L=zeros(m,n);
for i=1:m
for j=1:n
ind=i+j-1;
if ind<=m
L(i,ind)=B(i,j);
end
end
end
end
function y= bandforward(R,b)
dim=size(R);c=find(R(1:dim(1),1)~=0);d=c(end)-1;
n= length(b); y=b(:);
y(1)=y(1)/R(1,1);
for k=2: n
km= max (1,k-d);
y(k )=(y(k)-dot(R(k,km:k -1),y(km:k -1)))/ R(k,k);
end
end
function x= bandback(R,y)
dim=size(R);c=find(R(1,1:dim(1))~=0);d=c(end)-1;
n=length(y); x=y (:);
x(n)=x(n)/R(n,n);
for k=n-1: -1:1
kp= min(n,k+d);
x(k )=(x(k)-dot(R(k,k +1: kp ),x(k +1: kp )))/ R(k,k);
end
end
function ta = line2ta(line, ta1, deltat)
ta = ta1+(line - 1) *deltat;
end

