function rho = crafrate(lat,vx,vy,height,DCMel,earthflg,vertmech)
%該函數(shù)的功能是實(shí)現(xiàn)計(jì)算運(yùn)載體的旋轉(zhuǎn)角速度向量rho，輸入量是當(dāng)?shù)氐乩砭暥萳at、載體速度值、載體高度值、方向余弦矩陣DCMel，以及函數(shù)標(biāo)志數(shù)值ea
%rthflg和向量vertmech標(biāo)志值，輸出載體旋轉(zhuǎn)角速度值。
if nargin<7, vertmech = 0; end
if nargin<6,error('insufficient number of input arguments'),end
if (earthflg ~= 0) & (earthflg ~= 1), error('EARTHFLG not specified correctly'),end

if earthflg == 0,
   rm = 6378137;
   rp = rm;
else,
   [rm,rp] = radicurv(lat);%根據(jù)earthflg符號(hào)判斷地球模型的形態(tài)：是圓形還是橢球形？如果是圓形，則模型的長短半軸都取為6378137；
   %否則是橢球形，按照radicurv函數(shù)計(jì)算地球模型的長短半軸半徑；
end

if vertmech == 0,
   Ve = vx;%東向速度
   Vn = vy;%北向速度
   rho(1,1) = -Vn/(rm + height);
   rho(2,1) = Ve/(rp + height);
   rho(3,1) = Ve*tan(lat)/(rp + height);%計(jì)算運(yùn)載體旋轉(zhuǎn)角速度
else
   h = height;
   f = 1/298.257223563;
   [b,a] = radicurv(0);
   DCMelYx = DCMel(1,3);
   DCMelYy = DCMel(2,3);
   DCMelYz = DCMel(3,3);

   rho(1,1) = -(vy/a)*(1 - h/a - f*(1-3*DCMelYy*DCMelYy-DCMelYx*DCMelYx)) ...
      - (vx/a)*(2*f*DCMelYx*DCMelYy);
   rho(2,1) = (vx/a)*(1 - h/a - f*(1-3*DCMelYx*DCMelYx-DCMelYy*DCMelYy)) ...
      + (vy/a)*(2*f*DCMelYx*DCMelYy);

   if vertmech == 1,
     rho(3,1) = 0;
   end

   if vertmech == 2,
      OMEGA = 7.292115e-5;   %Earth's rate in rad/s
      rho(3,1) = -OMEGA*DCMelYz;
   end

   if vertmech == 3,
      num = -( rho(2,1)*DCMelYy + rho(1,1)*DCMelYx );
      if lat >= 0,
         den = DCMelYz + 1;
      else
         den = DCMelYz - 1;
      end
      rho(3,1) = num/den;
   end
   
end