matlab实现导航卫星系统中计算多普勒频移

用matlab实现导航系统中计算多普勒频移，输入接收机所在处的经纬高，（f0,d0,h0）以及间隔时间(t0)，和信号对应波长频率(f),k为分析的毫秒数（此处设为1即可），利用拉格朗日差值和单位观测向量的计算实现对多普勒频移的求解。

function fd = exchange（f0d0h0t0fk）
%输入接收机所在处的纬、径、高（f0d0h0）以及间隔时间（t0），和信号对应波长频率（f）k为分析的毫秒数（此处设为1即可）
%坐标转换，打开存放坐标的文件，再将之从‘大地坐标系’转换成‘地心地固坐标系’。B是地心地固坐标系。
%注意：此程序每次只能在两个相邻的数据点中间插入所需量的数据
%t0= 0.001;%s
%f0 = 114.3316;
%d0 = 30.5772;
%h0 = 3963; %mile   地球椭球体长半径（赤道半径）     1mile=1.609km
%f = 1023000;
 %   format long g;
    A = xlsread（‘MDatas_lla.xlsx‘）;
  %k是分析的毫秒数
  %  k = 1;
%===========================坐标转换开始====================================
    M = size（A1）;
    X = zeros（1M）;
    Y = zeros（1M）;
    Z = zeros（1M）;
    a = 3963;   % 长半径（赤道半径）   6378.2km
    b = 3950;   % 短半径（极轴半径）   6356.8km
    e2 = （a^2 - b^2）/a^2;
    B = zeros（M3）;
    for numb = 1:M
    %N是卯酉圈曲率半径
        N = a/sqrt（1 - e2*（sin（A（numb1）*pi/180））^2）;
        X（1numb） = （N + A（numb3）） * cos（A（numb1）*pi/180） * cos（A（numb2）*pi/180）;
        Y（1numb） = （N + A（numb3）） * cos（A（numb1）*pi/180） * sin（A（numb2）*pi/180）;
        Z（1numb） = （N*（1 - e2） + A（numb3）） * sin（A（numb1）*pi/180）;
        B（numb1） = X（1numb）;
        B（numb2） = Y（1numb）;
        B（numb3） = Z（1numb）;
    end
%本地坐标转换成地心地固坐标，C是本地坐标的地心地固坐标。
%n0也是卯酉圈曲率半径
    n0 = a/sqrt（1 - e2*（sin（（f0）*pi/180））^2）;
    x0 = （n0 + h0）*cos（f0*pi/180）*cos（d0*pi/180）;
    y0 = （n0 + h0）*cos（f0*pi/180）*sin（d0*pi/180）;
    z0 = （n0*（1-e2） + h0）*sin（f0*pi/180）;
    C =