main36.m

仿真了地球杂波的情况。考虑了地球自转带来的多普勒旋转等影响。仿真结果和理论相一致。可作为星载雷达地球杂波仿真的参考

Downside_Calculate_Result_Set1 = zeros(1, length(All_Intersect_Data_Phi_Radian));
Downside_Calculate_Result_Set2 = zeros(1, length(All_Intersect_Data_Phi_Radian));
% FSolve的选项
FSolve_Options = optimset('Display','off');  % Turn off Display
for w1 = 1 : 1 : length(All_Intersect_Data_Phi_Radian)
    % 按照第一种方式求解
    % 赋初值给Phi，这个是精确值
    Phi_Radian_Tmp = All_Intersect_Data_Phi_Radian(w1);
    % 赋初值给Theta
    Theta_Radian_Tmp = Downside_Intersect_Data_Theta_Radian(w1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun1, Theta_Radian_Tmp, FSolve_Options);
    Downside_Theta_Search_Result_Set1(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun1(Theta_Radian_Tmp);
    Downside_Calculate_Result_Set1(1, w1) = Tmp;

    % 按照第二种方式求解
    % 赋初值给Phi，这个是精确值
    Phi_Radian_Tmp = All_Intersect_Data_Phi_Radian(w1);
    % 赋初值给Theta
    Theta_Radian_Tmp = Downside_Intersect_Data_Theta_Radian(w1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun2, Theta_Radian_Tmp, FSolve_Options);
    Downside_Theta_Search_Result_Set2(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun2(Theta_Radian_Tmp);
    Downside_Calculate_Result_Set2(1, w1) = Tmp;
end;
figure;
subplot(2,2,1);
plot(All_Intersect_Data_Phi_Radian, Downside_Calculate_Result_Set1);
title('逆向验证结果(1)');
subplot(2,2,2);
plot(All_Intersect_Data_Phi_Radian, Downside_Theta_Search_Result_Set1);
xlabel('Phi');
ylabel('Theta');
title('对应求解得到的Theta(1)');
subplot(2,2,3);
plot(All_Intersect_Data_Phi_Radian, Downside_Calculate_Result_Set2);
title('逆向验证结果(2)');
subplot(2,2,4);
plot(All_Intersect_Data_Phi_Radian, Downside_Theta_Search_Result_Set2);
xlabel('Phi()');
ylabel('Theta()');

% 选择第二种方式的精确求解结果
Accurate_Intersect_Data_Phi_Radian = All_Intersect_Data_Phi_Radian;
% 获得的角度值
Accurate_Intersect_Data_Upside_Theta_Radian = Upside_Theta_Search_Result_Set2;
Accurate_Intersect_Data_Downside_Theta_Radian = Downside_Theta_Search_Result_Set2;
% 逆向验证的计算结果
Upside_Calculate_Result_Set = Upside_Calculate_Result_Set2;
Downside_Calculate_Result_Set = Downside_Calculate_Result_Set2;


% 根据逆向验证的结果判断是否需要对求解结果进行前后扩展，使得相交区域完整
% 对前部和后部进行足够的扩展
Extend_Width_Tmp = 20;
Accurate_Intersect_Data_Phi_Radian_Tmp = zeros(1, Extend_Width_Tmp*2+...
    length(Accurate_Intersect_Data_Phi_Radian));
Accurate_Intersect_Data_Upside_Theta_Radian_Tmp = zeros(size(Accurate_Intersect_Data_Phi_Radian_Tmp));
Accurate_Intersect_Data_Downside_Theta_Radian_Tmp = zeros(size(Accurate_Intersect_Data_Phi_Radian_Tmp));
Upside_Calculate_Result_Set_Tmp = zeros(size(Accurate_Intersect_Data_Phi_Radian_Tmp));
Downside_Calculate_Result_Set_Tmp = zeros(size(Accurate_Intersect_Data_Phi_Radian_Tmp));

Accurate_Intersect_Data_Phi_Radian_Tmp(Extend_Width_Tmp+1 : end-Extend_Width_Tmp) = ...
    Accurate_Intersect_Data_Phi_Radian(1:end); 
Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(Extend_Width_Tmp+1 : end-Extend_Width_Tmp) = ...
    Accurate_Intersect_Data_Upside_Theta_Radian(1:end);  
Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(Extend_Width_Tmp+1 : end-Extend_Width_Tmp) = ...
    Accurate_Intersect_Data_Downside_Theta_Radian(1:end);  
Upside_Calculate_Result_Set_Tmp(Extend_Width_Tmp+1 : end-Extend_Width_Tmp) = ...
    Upside_Calculate_Result_Set(1:end);  
Downside_Calculate_Result_Set_Tmp(Extend_Width_Tmp+1 : end-Extend_Width_Tmp) = ...
    Downside_Calculate_Result_Set(1:end);  

% 对Phi坐标范围进行前后的扩展
Phi_Radian_Step_Tmp = Accurate_Intersect_Data_Phi_Radian(2) - ...
    Accurate_Intersect_Data_Phi_Radian(1);
Tmp_Front = Accurate_Intersect_Data_Phi_Radian(1) - (1 : 1 : Extend_Width_Tmp)*...
    Phi_Radian_Step_Tmp;
Tmp_Rear = Accurate_Intersect_Data_Phi_Radian(end) + (1 : 1 : Extend_Width_Tmp)*...
    Phi_Radian_Step_Tmp;
Accurate_Intersect_Data_Phi_Radian_Tmp(1 : Extend_Width_Tmp) = Tmp_Front;
Accurate_Intersect_Data_Phi_Radian_Tmp(end-Extend_Width_Tmp+1 : end) = Tmp_Rear;
% 对Theta坐标的前部进行扩展
for w1 =  Extend_Width_Tmp : -1 : 1
    % 赋初值给Phi，这个是精确值
    Phi_Radian_Tmp = Accurate_Intersect_Data_Phi_Radian_Tmp(w1);
    % 对上部进行扩展
    % 赋初值给Theta
    Theta_Radian_Tmp = Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(w1+1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun2, Theta_Radian_Tmp, FSolve_Options);
    Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun2(Theta_Radian_Tmp);
    Upside_Calculate_Result_Set_Tmp(1, w1) = Tmp;

    % 对下部进行扩展
    % 赋初值给Theta
    Theta_Radian_Tmp = Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(w1+1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun2, Theta_Radian_Tmp, FSolve_Options);
    Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun2(Theta_Radian_Tmp);
    Downside_Calculate_Result_Set_Tmp(1, w1) = Tmp;
end;

% 对Theta坐标的后部进行扩展
for w1 =  (length(Accurate_Intersect_Data_Phi_Radian_Tmp)-Extend_Width_Tmp+1) : ...
        -1 : length(Accurate_Intersect_Data_Phi_Radian_Tmp)
    % 赋初值给Phi，这个是精确值
    Phi_Radian_Tmp = Accurate_Intersect_Data_Phi_Radian_Tmp(w1);
    % 对上部进行扩展
    % 赋初值给Theta
    Theta_Radian_Tmp = Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(w1-1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun2, Theta_Radian_Tmp, FSolve_Options);
    Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun2(Theta_Radian_Tmp);
    Upside_Calculate_Result_Set_Tmp(1, w1) = Tmp;

    % 对下部进行扩展
    % 赋初值给Theta
    Theta_Radian_Tmp = Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(w1-1);
    [Theta_Radian_Tmp,Fval,exitflag] = fsolve(@CalThetaFun2, Theta_Radian_Tmp, FSolve_Options);
    Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(1, w1) = Theta_Radian_Tmp;
    % 对求解结果进行验证
    Tmp = CalThetaFun2(Theta_Radian_Tmp);
    Downside_Calculate_Result_Set_Tmp(1, w1) = Tmp;
end;

Tmp = Accurate_Intersect_Data_Upside_Theta_Radian_Tmp - Accurate_Intersect_Data_Downside_Theta_Radian_Tmp;
[Index_Tmp] = find(Tmp > 0.001);
Index_Tmp = [Index_Tmp(1)-1 Index_Tmp Index_Tmp(end)+1];
% Index_Tmp = (1 : 1 : length(Accurate_Intersect_Data_Upside_Theta_Radian_Tmp));

Accurate_Intersect_Data_Phi_Radian = Accurate_Intersect_Data_Phi_Radian_Tmp(1, Index_Tmp);
Accurate_Intersect_Data_Upside_Theta_Radian = Accurate_Intersect_Data_Upside_Theta_Radian_Tmp(1, Index_Tmp);
Accurate_Intersect_Data_Downside_Theta_Radian = Accurate_Intersect_Data_Downside_Theta_Radian_Tmp(1, Index_Tmp);
Upside_Calculate_Result_Set = Upside_Calculate_Result_Set_Tmp(1, Index_Tmp);
Downside_Calculate_Result_Set = Downside_Calculate_Result_Set_Tmp(1, Index_Tmp);
Accurate_Intersect_Data_Upside_Theta_Radian(1) = ...
    Accurate_Intersect_Data_Downside_Theta_Radian(1);
Accurate_Intersect_Data_Upside_Theta_Radian(end) = ...
    Accurate_Intersect_Data_Downside_Theta_Radian(end);

clear Upside_Theta_Search_Result_Set1 Downside_Theta_Search_Result_Set1 ...
    Upside_Calculate_Result_Set1 Downside_Calculate_Result_Set1;
clear Upside_Theta_Search_Result_Set2 Downside_Theta_Search_Result_Set2 ...
    Upside_Calculate_Result_Set2 Downside_Calculate_Result_Set2;
clear Accurate_Intersect_Data_Phi_Radian_Tmp Accurate_Intersect_Data_Upside_Theta_Radian_Tmp ...
    Accurate_Intersect_Data_Downside_Theta_Radian_Tmp Upside_Calculate_Result_Set_Tmp ...
    Downside_Calculate_Result_Set_Tmp;
    

% 处理便于显示，没有数据处理意义上的作用
Accurate_Intersect_Data_Phi_Degree = Accurate_Intersect_Data_Phi_Radian * 180 / pi;
Accurate_Intersect_Data_Upside_Theta_Degree = Accurate_Intersect_Data_Upside_Theta_Radian * 180 / pi;
Accurate_Intersect_Data_Downside_Theta_Degree = Accurate_Intersect_Data_Downside_Theta_Radian * 180 / pi;

figure;
subplot(2,1,1);
plot(Accurate_Intersect_Data_Phi_Degree, Accurate_Intersect_Data_Upside_Theta_Degree, 'r');
hold on;
plot(Accurate_Intersect_Data_Phi_Degree, Accurate_Intersect_Data_Downside_Theta_Degree, 'b');
grid on;
xlabel('Phi');
ylabel('Theta');
title('获得的精确经纬度');
% max(Accurate_Intersect_Data_Phi_Degree) - min(Accurate_Intersect_Data_Phi_Degree)
% max(Accurate_Intersect_Data_Upside_Theta_Degree) - min(Accurate_Intersect_Data_Downside_Theta_Degree)
clear Accurate_Intersect_Data_Phi_Degree Accurate_Intersect_Data_Phi_Degree ...
    Accurate_Intersect_Data_Upside_Theta_Degree Accurate_Intersect_Data_Downside_Theta_Degree;

subplot(2,1,2);
plot(Accurate_Intersect_Data_Phi_Radian, Upside_Calculate_Result_Set, 'r');
hold on;
plot(Accurate_Intersect_Data_Phi_Radian, Downside_Calculate_Result_Set, 'b');
grid on;
title('逆向检验误差');
clear Upside_Calculate_Result_Set Downside_Calculate_Result_Set;

% 根据得到的经度和纬度打印图形
X_Tmp = Const_Earth_Re * cos(Accurate_Intersect_Data_Upside_Theta_Radian) .*...
    cos(Accurate_Intersect_Data_Phi_Radian);
Y_Tmp = Const_Earth_Re * cos(Accurate_Intersect_Data_Upside_Theta_Radian) .*...
    sin(Accurate_Intersect_Data_Phi_Radian);
Z_Tmp = Const_Earth_Re * sin(Accurate_Intersect_Data_Upside_Theta_Radian);
Axis_Range = [-Const_Earth_Re Const_Earth_Re -Const_Earth_Re Const_Earth_Re ...
    -Const_Earth_Re Const_Earth_Re];
figure;
plot3(X_Tmp, Y_Tmp, Z_Tmp);
axis(Axis_Range);
title('将所得到的经度纬度用图形表示出来');
hold on;
X_Tmp = Const_Earth_Re * cos(Accurate_Intersect_Data_Downside_Theta_Radian) .*...
    cos(Accurate_Intersect_Data_Phi_Radian);
Y_Tmp = Const_Earth_Re * cos(Accurate_Intersect_Data_Downside_Theta_Radian) .*...
    sin(Accurate_Intersect_Data_Phi_Radian);
Z_Tmp = Const_Earth_Re * sin(Accurate_Intersect_Data_Downside_Theta_Radian);
Axis_Range = [-Const_Earth_Re Const_Earth_Re -Const_Earth_Re Const_Earth_Re ...
    -Const_Earth_Re Const_Earth_Re];
plot3(X_Tmp, Y_Tmp, Z_Tmp);
axis(Axis_Range);
grid on;
clear X_Tmp Y_Tmp Z_Tmp;

%% +++++++++++++++++++++++++++++++++++++
% 根据求取的经度和纬度方位，获得经度和纬度的分割间距
Num_Tmp = 100;
% 距离对Theta的求导
Delta_Rst_Theta_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
% 距离对Phi的求导
Delta_Rst_Phi_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
% 速度对Theta的求导
Delta_Vst_Theta_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
% 速度对Phi的求导
Delta_Vst_Phi_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
% 相对径向距离
Rst_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
% 相对径向速度
Vst_Result = zeros(length(Accurate_Intersect_Data_Phi_Radian), Num_Tmp);
for w1 = 1 : 1 : length(Accurate_Intersect_Data_Phi_Radian)
    Target_Phi = Accurate_Intersect_Data_Phi_Radian(w1);
    for w2 = 1 : 1 : Num_Tmp
        Target_Theta = Accurate_Intersect_Data_Downside_Theta_Radian(w1) + ...
            ((Accurate_Intersect_Data_Upside_Theta_Radian(w1) - ...
            Accurate_Intersect_Data_Downside_Theta_Radian(w1)) / Num_Tmp) * (w2-1);

        [Delta_Rst_Theta, Delta_Rst_Phi, Delta_Vst_Theta, Delta_Vst_Phi, Rst, Vst] = ...
            Calculate_R_V_DeltaR_DeltaV(Satellite_r, Const_Earth_Re, Satellite_Theta_s_Radian, ...
            Satellite_Phi_s_Radian, Satellite_Theta_i_Radian, Earth_We, Satellite_Wp, ...
            Target_Theta, Target_Phi);

        % 保存相关计算结果
        Delta_Rst_Theta_Result(w1, w2) = Delta_Rst_Theta;
        Delta_Rst_Phi_Result(w1, w2) = Delta_Rst_Phi;
        Delta_Vst_Theta_Result(w1, w2) = Delta_Vst_Theta;
        Delta_Vst_Phi_Result(w1, w2) = Delta_Vst_Phi;
        Rst_Result(w1, w2) = Rst;
        Vst_Result(w1, w2) = Vst;
    end;
end;

% 将相对速度变换导数转换为相对多普勒频移
Delta_Fd_Theta_Result = -2*Delta_Vst_Theta_Result / Const_Wave_Length;
Delta_Fd_Phi_Result = -2*Delta_Vst_Phi_Result / Const_Wave_Length;
clear Delta_Vst_Theta_Result Delta_Vst_Theta_Result;

figure;
subplot(3,2,1);
mesh(Rst_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Rst');
subplot(3,2,4);
mesh(Vst_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Vst');

subplot(3,2,2);
mesh(Delta_Rst_Theta_Result);
xlabel('Phi');
ylabel('Theta');
zlabel('Delt_R/Delt_Theta');
subplot(3,2,5);