matlab2013激活文件

clc;clear;
% INS\Gps位置速度综合（静基座）。
%初值定义区
h（1）=1000;L（1）=30*pi/180;Lmd（1）=120*pi/180;%位置误初值。依次为高度，纬度，经度
Ve（1）=300*sin（35*pi/180）;Vn（1）=300*cos（35*pi/180）;Vu（1）=0;%初始速度，依次为，东，北，天速度。
yaw（1）=600/60/60*pi/180;pitch（1）=300/60/60*pi/180;roll（1）=300/60/60*pi/180;;% 初始姿态，依次为航向，俯仰，横滚角。
Eh（1）=50;EL（1）=0.000446*pi/180;ELmd（1）=0.000446*pi/180;%位置误初值误差。依次为高度，纬度，经度误差
EVe（1）=0.6;EVn（1）=0.6;EVu（1）=0.6;%初始速度误差，依次为，东，北，天速度误差。
Eyaw（1）=0.2*pi/180;Epitch（1）=0.1*pi/180;Eroll（1）=0.1*pi/180;% 初始姿态误差，依次为航向，俯仰，横滚角误差。


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%常量定义区，
wie=7.292115147e-5;%地球自转角速度。
Re=6378254;%长半轴
e=1/298.3;%椭圆偏心率。
g0=9.7803267714;gk1=0.00193185138639;gk2=0.00669437999013;%求取地球重力加速度，与当纬度，高有关。
T=0.2;%迭代周期
g=g0*（1+gk1*sin（L（1））^2）*（1-2*h（1）/Re）/sqrt（1-gk2*sin（L（1））^2）;
Ax=g*10^-4;Ay=g*10^-4;Az=g*10^-4;%加速度计零偏
eTbx=0.1*pi/180/3600;eTby=0.1*pi/180/3600;eTbz=0.1*pi/180/3600;%陀螺常值漂移，0.1度/小时
eTrx=0.1*pi/180/3600;eTry=0.1*pi/180/3600;eTrz=0.1*pi/180/3600;%陀螺随机漂移，0.1度/小时
Trx=300;Try=300;Trz=300;Tax=300;Tay=300;Taz=300;%  陀螺相关时间，加速度相关时间。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%kalman计算X=[姿态，速度，位置，陀螺常值漂移，陀螺随机漂移，加计零偏]为状态（误差）

X（11）=Eyaw（1）;X（21）=Epitch（1）;X（31）=Eroll（1）;X（41）=EVe（1）;X（51）=EVn（1）;X（61）=EVu（1）;X（71）=EL（1）;X（81）=ELmd（1）;X（91）=Eh（1）;X（101）=eTbx;X（111）=eTby;X（121）=eTbz;X（131）=eTrx;X（141）=eTry;X（151）=eTrz;X（161）=Ax;X（171）=Ay;X（181）=Az;
P0=diag（[（1*pi/180/3600）^2（1*pi/180/3600）^2（1*pi/180/3600）^2（0.1）^2（0.1）^2（0.1）^2（0.000446）^2（0.000446）^2（50）^2（0.1*pi/180/3600）^2（0.1*pi/180/3600）^2（0.1*pi/180/3600）^2（0.1*pi/180/3600）^2（0.1*pi/180/3600）^2（0.1*pi/180/3600）^2（0.0001*g）^2（0.0001*g）^2（0.0001*g）^2]）;
Q0=diag（[（0.05*pi/180/3600）^2（0.05*pi/180/3600）^2（0.05*pi/180/3600）^2（0.1/57.2958/3600）^2/1800（0.1/57.2958/3600）^2（0.1/57.2958/3600）^2（10^-4*g）^2（10^-4*g）^2（10^-4*g）^2]）;
R0=diag（[（0.0002883*pi/180）^2（0.0002883*pi/180）^2（32）^2（0.05）^2（0.05）^2（0.05）^2]）;
H=zeros（618）;
%fu=8*randn（13601）;fn=8*randn（13601）;fe=8*randn（13601）;
for k=1:3600
    Rn=Re*（1+e*sin（L（k））^2）;%主曲率半径Rn
    Rm=Re*（1-2*e+3*e*sin（L（k））^2）;%主曲率半径Rm
    g=g0*（1+gk1*sin（L（k））^2）*（1-2*h（k）/Re）/sqrt（1-gk2*sin（L（k））^2）;
    Cnb=[cos（roll（k））*cos（yaw（k））+sin（roll（k））*sin（yaw（k））*sin（pitch（k）） sin（yaw（k））*cos（pitch（k）） sin（roll（k））*cos（yaw（k））-cos（roll（k））*sin（yaw（k））*sin（roll（k））
    -cos（roll（k））*sin（yaw（k））+sin（roll（k））*cos（yaw（k））*sin（pitch（k）） cos（yaw（k））*sin（pitch（k）） -sin（roll（k））*sin（yaw（k））-cos（roll（k））*cos（yaw（k））*sin（pitch（k））
    -sin（roll（k））*cos（pitch（k）） sin（pitch（k）） cos（roll（k））*cos（pitch（k））];
    wiet=[0;0;wie];wet=[-Vn（k）/Rm;Ve（k）/Rn;Ve（