function [outputArg1,outputArg2] = ctSent_setpram(N,window)
%CTSENT_SBAS Summary of this function goes here
%   Detailed explanation goes here
savename=['ps',num2str(psver)];
incsavename=['inc',num2str(psver)];
hgtsavename=['hgt',num2str(psver)];
phsavename=['ph',num2str(psver)];
n_ds = length(ij);
matlab_version = version('-release');           % [DB] getting the matlab version
matlab_version = str2num(matlab_version(1:4));  % [DB] the year

day_yyyymmdd = date;
year=floor(day_yyyymmdd/10000);
month=floor((day_yyyymmdd-year*10000)/100);
monthday=day_yyyymmdd-year*10000-month*100;
slave_day=datenum(year,month,monthday);
[slave_day,day_ix]=sort(slave_day);
day_yyyymmdd=day_yyyymmdd(day_ix);

master_day_yyyymmdd = str2double(info_mas.date);
year=floor(master_day_yyyymmdd/10000);
month=floor((master_day_yyyymmdd-year*10000)/100);
monthday=master_day_yyyymmdd-year*10000-month*100;
master_day=datenum(year,month,monthday);
master_ix=sum(master_day>slave_day)+1;
day = slave_day; % insert master day 
n_image=size(day,1);
n_ifg = length(order);

heading=str2double(info_mas.sensor_para.Heading);
setparm('heading',heading,1);
lambda = info_mas.sensor_para.Wavelength;
setparm('lambda',lambda,1);
hgt = freadCt(geo.hgt,lenth,width);
hgt = hgt(window(1):window(2),window(3):window(4));
hgt1=hgt(sub2ind(size(hgt),ij(:,2),ij(:,3)));
stamps_save(hgtsavename,hgt1);
clear hgt
lon = freadCt(geo.lon,lenth,width);
lon = lon(window(1):window(2),window(3):window(4));
lon1=lon(sub2ind(size(lon),ij(:,2),ij(:,3)));
lat = freadCt(geo.lat,lenth,width);
lat = lat(window(1):window(2),window(3):window(4));
lat1=lat(sub2ind(size(lat),ij(:,2),ij(:,3)));
lonlat = [lon1,lat1];
clear lon;clear lat; clear lon1;clear lat1;
ll0 = (max(lonlat)+ min(lonlat))/2;
xy = llh2local(lonlat',ll0)*1000;
xy = xy'; 
theta=(180-heading)*pi/180;
if theta>pi
    theta=theta-2*pi;
end
rotm=[cos(theta),sin(theta); -sin(theta),cos(theta)];
xy=xy';
xynew=rotm*xy; % rotate so that scene axes approx align with x=0 and y=0
if max(xynew(1,:))-min(xynew(1,:))<max(xy(1,:))-min(xy(1,:)) &...
   max(xynew(2,:))-min(xynew(2,:))<max(xy(2,:))-min(xy(2,:))
    xy=xynew; % check that rotation is an improvement
    disp(['Rotating by ',num2str(theta*180/pi),' degrees']);
end      
xy=single(xy');
[dummy,sort_ix]=sortrows(xy,[2,1]); % sort in ascending y order
xy=xy(sort_ix,:);
xy=[[1:n_ds]',xy];
xy(:,2:3)=round(xy(:,2:3)*1000)/1000; % round to mm
lonlat=lonlat(sort_ix,:);
temp = ij(:,2);
ij(:,2)=temp(sort_ix);
temp = ij(:,3);
ij(:,3)=temp(sort_ix);
clear temp;
for i = 1:n_image
    temp = ph(:,i);
    ph(:,i) = temp(sort_ix);
end
clear temp;
inc = freadCt(geo.inc,lenth,width);
inc = inc(window(1):window(2),window(3):window(4));
data_inc = inc(sub2ind(size(inc),ij(:,2),ij(:,3)));
clear inc;
inc = data_inc*pi./180;
stamps_save(incsavename,inc)
ps_parms_initial;
bpsavename=['bp',num2str(psver)];
baseDiff = load('baseline');
baseDiffpath = baseDiff.baseline;
bperp_mat = zeros(n_ds,n_image,'single');
for i = 2:length(baseDiffpath)
    bperpmat  =  real(readcpx(baseDiffpath(i).savepath,width,lenth));
    bperpmat= bperpmat(window(1):window(2),window(3):window(4));
    bperp_mat(:,i) = bperpmat(sub2ind(size(bperpmat),ij(:,2),ij(:,3)));
    fprintf('Extracting prep baseline... %d\n',i-1);
end
bperp = squeeze(mean(bperp_mat));
bperp = bperp';
bperp_mat(:,master_ix)=[];
stamps_save(bpsavename,bperp_mat);
setparm('heading',heading,1);
setparm('lambda',lambda,1);
setparm('unwrap_gold_n_win',10);
setparm('unwrap_gold_alpha',1);
ph_rc = ph;
stamps_save(phsavename,ph);
stamps_save('rc1',ph_rc);
save(savename,'ij','lonlat','xy','day','bperp','master_day','master_ix','n_image','n_ifg','n_ps','ll0','master_ix');
save psver psver
end

function  [order,date,intfdate] = sequential_netwrok(stacksize,N,slc)
%N: degree
%stack size
m=1;
 for i = 1:stacksize
     date(i) = str2double(slc(i).savepath(end-11:end-4));
 end
for i=1:stacksize-1
    for j=i+1:i+N
        if j>stacksize
            continue;
        end
   order(m,1)=i;order(m,2)=j; 
   m=m+1;
    end
end
intfdate = date(order);
end
function [slcstack] = ctSent_readDeb(window)
slcpath = load('debSlc.mat');
debSlc = slcpath.slcsavepath;
if nargin<1
    window =[1 slcpath.Lines 1 slcpath.Pixels];
    % row_start,row_end,column_start,column_end
end
slcstack = zeros(window(2)-window(1)+1,window(4)-window(3)+1,length(debSlc),'single');
for i=1:length(debSlc)
    [data, ~] = freadbk(debSlc(i).savepath,slcpath.Lines,'mph',window(1),window(2),window(3),window(4));
    %data = data(window(1):window(2),window(3):window(4));
    slcstack(:,:,i) = single(data);
    fprintf('Reading SLC %d / %d \n',i,length(debSlc));
end
end
function [SHP]=SHP_SelPoint(mlistack,CalWin,Alpha)
if nargin < 3
    Alpha = 0.05;
end
if nargin < 2
    CalWin = [11 11];
end
if nargin < 1
    help SHP_SelPoint
    return;
end
tic;
if length(size(mlistack))~=3
    error('Please input 3D matrix...');
end
[nlines,nwidths,npages] = size(mlistack);
mlistack=single(mlistack);
%Parameter prepare:
SHP.CalWin = CalWin;
RadiusRow=(CalWin(1)-1)/2;
RadiusCol=(CalWin(2)-1)/2;
InitRow=(CalWin(1)+1)/2; % InitRow is CenterRow
InitCol=(CalWin(2)+1)/2; % InitCol is CenterCol
LRT_nl = 3;
LRT_nw = 3;
if RadiusRow<LRT_nl
    LRT_nl=1;
end
if RadiusCol<LRT_nw
    LRT_nw=1;
end
% Select pixels randomly
% LRT_index = randsample(CalWin(1)*CalWin(1),(2*LRT_nl+1)*(2*LRT_nw+1)); 
%Statistical threshold:
CR_lo = finv(Alpha/2,2*npages,2*npages);
CR_up = finv(1-Alpha/2,2*npages,2*npages);
Galpha_L = gaminv(Alpha/2,npages,1);
Galpha_U = gaminv(1-Alpha/2,npages,1);
%Edeg mirror-image
% mlistack = padarray(mlistack,[RadiusRow RadiusCol],'symmetric');
meanmli = mean(mlistack,3);
PixelInd=false(CalWin(1)*CalWin(2),nlines*nwidths);
PixelSub=zeros(CalWin(1)*CalWin(2),nlines*nwidths,'single');
%estimate SHPs
num=1;
p=1;
all = nlines*nwidths;
all_step = floor(all/100);
    for kk=1:nwidths
        for ll=1:nlines      
            %Initial estimation (Likelihood-ratio test)
            [flag] = InBox(nlines,nwidths,CalWin(1),CalWin(2),ll,kk);
            if ~flag
                num = num +1;
            else
            temp = meanmli(ll-LRT_nl:ll+LRT_nl,kk-LRT_nw:kk+LRT_nw);
            T = meanmli(ll,kk)./temp;
            T = T>CR_lo&T<CR_up;
            SeedPoint = mean(temp(T));
            %iteration (Gamma Confidence interval)
            MeanMatrix = meanmli(ll-RadiusRow:ll+RadiusRow,kk-RadiusCol:kk+RadiusCol);
            [pixel,line]=meshgrid(kk-RadiusCol:kk+RadiusCol,ll-RadiusRow:ll+RadiusRow);
            ind = sub2ind(size(meanmli),line(:),pixel(:));
            SeedPoint = MeanMatrix>Galpha_L*SeedPoint/npages&MeanMatrix<Galpha_U*SeedPoint/npages; %check membership
            SeedPoint(InitRow,InitCol)=true;
            %connection
            LL = bwlabel(SeedPoint); %double
            PixelFlag = LL(:)==LL(InitRow,InitCol);
            PixelInd(:,num) = PixelFlag.*ind;
            PixelSub(:,num) = single(PixelFlag).*ind;
            num=num+1;
            end
            if num == all_step * p;
                disp(['progress: ', num2str(1*p),'%']);
                p = p+1;
            end
        end
    end
%SHPs map 
SHP.PixelInd=PixelInd;
BroNum = sum(PixelInd,1);
BroNum = reshape(BroNum(:),[nlines,nwidths]);
SHP.BroNum = single((BroNum-1));  
SHP.PixelSub = PixelSub;
t=toc;
disp(['SHP_SelPoint operation completed in ',int2str(t),' second(s).']);
end
function A = freadCt(savepath,lines,width)
fileID = fopen(savepath);
A = fread(fileID,[lines width],['float','=>single']);
fclose(fileID);
end
function [Bperpmat] = ctSent_computeBnormal(refinfo,info_m,info_s,dem_read)
parm = load('parms');
orbit_Degree = parm.orbit_degree;
morbit = orbit(info_m,orbit_Degree);
sorbit = orbit(info_s,orbit_Degree);
[DEM_ex,Ellipse] = computeDEMbox(dem_read);
[time_m] = xyz2t(Ellipse,info_m,morbit);
Msat = getXYZ(time_m.timeAzimuth,morbit);
[time_s] = xyz2t(Ellipse,info_s,sorbit);
Ssat = getXYZ(time_s.timeAzimuth,sorbit);
[Bperp] = BBparBperpTheta(Msat,Ellipse,Ssat);
timeref = xyz2t(Ellipse,refinfo);
[mline,mpixel] = ta2lp(info_m,time_m,firstlinetime,firstrangetime);
[mline_grid,mpixel_grid] = gridlinepixel(lines,pixels);
[Bperpmat] = griddata(mline,mpixel,Bperp,mline_grid,mpixel_grid);
end
function [Bperp] = BBparBperpTheta(M, P, S)
[m,n] = size(M.X);
Master=[M.X(:),M.Y(:),M.Z(:)];
Point=[P.X(:),P.Y(:),P.Z(:)];
Slave=[S.X(:),S.Y(:),S.Z(:)];
B = normalvec(Master - Slave); 
range1 = normalvec(Master - Point);
range2 = normalvec(Slave - Point);
Bpar = range1 - range2;
r1 = Master - Point;
r2 = Slave - Point;
Inc = vec_angle(Point, r1);
Bperp = B.^2 - Bpar.^2;
Bperp(Bperp < 0.0) = 0;
Bperp(Inc > vec_angle(Point, r2)) = sqrt(Bperp(Inc > vec_angle(Point, r2)));
Bperp(Inc <= vec_angle(Point, r2)) = 0 - sqrt(Bperp(Inc <= vec_angle(Point, r2)));
Bperp = reshape(Bperp,[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 [ line,pixel ] = ta2lp( info_m,t,burstFirstLineTime,t_Range1)
%PRF=1.717128973878037e+03;		
line=1+(t.timeAzimuth-burstFirstLineTime)/info_m.sensor_para.line_time_interval;
pixel=1+(t.timeRange-t_Range1)*299792458/info_m.subSwath.rangePixelSpacing;
% pixel=1+(t.timeRange-info_m.sensor_para.t_Range1)*info_m.sensor_para.RSR2x;
end
function [line,pixel] = gridlinepixel(line1,pixel1)

azimuthline = (1:line1)';
rangepixel = 1:pixel1;
line = repmat(azimuthline,1,length(rangepixel));
pixel = repmat(rangepixel,length(azimuthline),1);
end
function [DEM_ex,Ellipse] = computeDEMbox(dem_read)
latMin = min(dem_read.lat(:));latMax = max(dem_read.lat(:));
lonMin = min(dem_read.lon(:));lonMax = max(dem_read.lon(:));
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;
DEM_ex.lon = dem_read.lon;
DEM_ex.z = dem_read.z;
EGM96correction = griddata(lat,lon,heightCorrection,DEM_ex.lat,DEM_ex.lon);
DEM_ex.z = double(DEM_ex.z)+EGM96correction;
Ellipse = BLH2XYZ(DEM_ex.lat,DEM_ex.lon,DEM_ex.z);
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

