function pixel_ll = RS_backward_compute_rpc(f_in_name, lsh, value_bijin)
%author  lizijia
%2022-4-28
% 读取遥感RPB文件，输入点的[line,sample]，获取像素点对应的影像经纬度值。
%当然会需要SRTM，此处不作为输入，在主程序中进行读取，取值

%f_in_name： 输入RPB,RPC文件名称
%lsh： [line,sample,height]
%value_bijin: 用于迭代更新用的。初值此工具默认使用[lat_offset,lon_offset]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用法1： pixel_ll=RS_backward_compute_rpc(f_in_name,lsh)%默认使用lat_offset和lon_offset为value_bijin
%用法2：
%pixel_ll=RS_backward_compute_rpc(f_in_name,lsh,value_bijin)%value_bijin为迭代出来更接近真值的该点经纬度值

line = lsh(1);
sample = lsh(2);
height = lsh(3);

% 1 读取RPB文件内容
fid1=fopen(f_in_name,'r');

read_one_line = fgetl(fid1);
read_one_line = fgetl(fid1);
read_one_line = fgetl(fid1);
read_one_line = fgetl(fid1);
read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',});
errBias=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',});
errRand=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',});
lineOffset =char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',});
sampOffset=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',});
latOffset =char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
longOffset=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
heightOffset=char(temp(2));


read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
linescale=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
samplescale=char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
latScale  =char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
longScale =char(temp(2));

read_one_line = fgetl(fid1);
temp=strsplit(read_one_line,{';','=',')',','});
heightScale =char(temp(2));

read_one_line = fgetl(fid1);
for i=1:20
    read_one_line = fgetl(fid1);
    temp=strsplit(read_one_line,{';','=',')',','});
    lineNumCoef{i}=temp(1);
end
read_one_line = fgetl(fid1);
for i=1:20
    read_one_line = fgetl(fid1);
    temp=strsplit(read_one_line,{';','=',')',','});
    lineDenCoef{i}=temp(1);
end
read_one_line = fgetl(fid1);
for i=1:20
    read_one_line = fgetl(fid1);
    temp=strsplit(read_one_line,{';','=',')',','});
    sampNumCoef{i}=temp(1);
end
read_one_line = fgetl(fid1);
for i=1:20
    read_one_line = fgetl(fid1);
    temp=strsplit(read_one_line,{';','=',')',','});
    sampDenCoef{i}=temp(1);
end
fclose(fid1);

num_l_coef=[];
for i=1:20
    num_l_coef=[num_l_coef,str2num(char(lineNumCoef{i}))]  ;
end
den_l_coef=[];
for i=1:20
    den_l_coef=[den_l_coef,str2num(char(lineDenCoef{i}))]  ;
end

num_s_coef=[];
for i=1:20
    num_s_coef=[num_s_coef,str2num(char(sampNumCoef{i}))]  ;
end
den_s_coef=[];
for i=1:20
    den_s_coef=[den_s_coef,str2num(char(sampDenCoef{i}))] ;
end


% nargin 表示参数个数
if nargin==3
    value_bijin=value_bijin;
elseif nargin==2
    value_bijin=[str2double(latOffset) str2double(longOffset)];
else
    error('输入有错啊！');
end

%value_bijin=[str2double(latOffset) str2double(longOffset)];
% X=(sample-str2double(sampOffset))/str2double(samplescale) ;
% Y=(line-str2double(lineOffset))/str2double(linescale);
H=(height-str2double(heightOffset))/str2double(heightScale);
syms lat lon
%li=[lat, H]
%  P==((lat-str2double(latOffset))/str2double(latScale));
%  L==((lon-str2double(longOffset))/str2double(longScale));
vars=[lat lon];
% xishu=[1,((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale)),H,((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),H*H,((lat-str2double(latOffset))/str2double(latScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*H,H*H*H];
equns=[line==((num_l_coef*[1,((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale)),H,((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),H*H,((lat-str2double(latOffset))/str2double(latScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*H,H*H*H]')/(den_l_coef*[1,((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale)),H,((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),H*H,((lat-str2double(latOffset))/str2double(latScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*H,H*H*H]'))*str2double(char(linescale))+str2double(char(lineOffset)),sample==(num_s_coef*[1,((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale)),H,((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),H*H,((lat-str2double(latOffset))/str2double(latScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*H,H*H*H]')/(den_s_coef*[1,((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale)),H,((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),H*H,((lat-str2double(latOffset))/str2double(latScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale)),((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lon-str2double(longOffset))/str2double(longScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale)),((lat-str2double(latOffset))/str2double(latScale))*H*H,((lon-str2double(longOffset))/str2double(longScale))*((lon-str2double(longOffset))/str2double(longScale))*H,((lat-str2double(latOffset))/str2double(latScale))*((lat-str2double(latOffset))/str2double(latScale))*H,H*H*H]')*str2double(char(samplescale))+str2double(char(sampOffset))];
[lat,lon]=vpasolve(equns,vars,value_bijin);

%[lat,lon]=solve(equns,vars,'ReturnConditions', false,'Real',true);

% P=(line-str2num(latOffset))/str2num(latScale);
% L=(sample-str2num(longOffset))/str2num(longScale);
% H=(height-str2num(heightOffset))/str2num(heightScale);
% X,Y的归一化公式
%X=(sample-str2num(sampOffset))/str2num(samplescale);
%Y=(line-str2num(lineOffset))/str2num(linescale);


%line=((num_l_coef*xishu')/(den_l_coef*xishu'))*str2num(char(linescale))+str2num(char(lineOffset));
%sample=(num_s_coef*xishu')/(den_s_coef*xishu')*str2num(char(samplescale))+str2num(char(sampOffset));

pixel_ll=[double(lat),double(lon)];
end


