evalcp.m

航天工程工具箱

function cp=evalcp(comp)
%EVALCP  Evaluate center of pressure CP.
%   CP = EVALCP(COMP)  evaluates the center of pressure CP for
%   a rocket by applying the Barrowman method (NARAM-8 1966).
%   The following criteria have to be fulfilled in order
%   to get reliable results when using this method:
%
%      1. The angle-of-attack of the rocket is near zero
%         (less than 10 deg).
%      2. The speed of the rocket is much less than the
%         speed of sound.
%      3. The air flow over the rocket is smooth and does
%         not change rapidly.
%      4. The rocket is thin compared to its length.
%      5. The nose of the rocket comes smoothly to a point,
%         (however it can be approximated).
%      6. The rocket is an axially symmetric rigid body.
%      7. The fins are thin flat plates.
%
%   The distance CP is counted from the tip of the nose.
%   COMP is a (Nx11) matrix where each row in COMP can be
%   any of the following rows in the table below:
%
%   TYPE SUBTYPE POSITION LENGTH ARG5 ARG6 ARG7 ARG8 ARG9 REL-CG MASS
%     1     1     delta    len    0    d    -    -    -    rpos   w
%     1     2     delta    len    0    d    -    -    -    rpos   w
%     1     3     delta    len    0    d    -    -    -    rpos   w
%     2     -      pos     len    0    d    -    -    -    rpos   w
%     3     1      pos     len   d1   d2    -    -    -    rpos   w
%     3     2      pos     len   d1   d2    -    -    -    rpos   w
%     4     -      pos      a     b    R    S    m    n    rpos   w
%     5     -      pos     len    -    d    -    -    n    rpos   w
%
%   where the last two columns REL-CG and MASS is ignored by this
%   function (they are used by EVALCG). Dashes are representing NaNs.
%   The meaning of TYPE and SUBTYPE is given in the following table:
%
%      TYPE SUBTYPE      COMPONENT
%        1     1      :   Nose, conical
%        1     2      :   Nose, ogive
%        1     3      :   Nose, parabolic
%        2     -      :   Body tube
%        3     1      :   Conical shoulder
%        3     2      :   Conical boattail
%        4     -      :   Fins
%        5     -      :   Other (internal)
%
%   where TYPE=2,5 is ignored by this function.
%   The meaning of the abbreviations in the first table above is the
%   following:
%
%      delta : For approximations of the nose cone
%              this is the distance from the real tip and the
%              approximated tip
%      pos   : Position of component counted from real tip of nose
%      rpos  : Position of relative CG counted from the top of the
%              component
%      len   : Length of component (or length of approx nose cone)
%      d     : Diameter at the nose base or diameter of tube or
%              "other"
%      d1    : Upper diameter
%      d2    : Lower diameter
%      n     : Number of fins (3,4 or 6), or number of parts
%      w     : Weight/mass of component
%
%   and for the fins we have that:
%
%             R    S
%       |   :<->|<-->:
% ......|...:...|\...:....
%  ^    |   :   |#\  :   ^
%  |    |   :   |##\ :   | m
%  |    |   :   |###\:   |
%  | a  |   :   |####\...v
%  |    |   :   |####|   ^
%  |    |   :   |####|   | b
% .v....|_______|####|   |
%                 "##|   |
%                   "|...v
%
%             
%   See also EVALCG, RCKPLOT, EVALDCM, SMARG.

% References : Jim Barrowman, "Calculating the Center of Pressure"
%                             Centuri Technical Information Report 33.

% Copyright (c) 2003-06-22, B. Rasmus Anthin.
% Revision 2003-06-23, 2003-06-24, 2003-06-26,
%          2003-07-23.

if sum(comp(:,1)==1)~=1
   error('Must have exactly one nose cone.');
end
d=comp(comp(:,1)==1,6);

for i=1:size(comp,1)
   pos=comp(i,3);
   switch(comp(i,1))
   case 1
      dL=pos;
      L=comp(i,4);
      CN(i)=2;
      switch(comp(i,2))
      case 1
         x(i)=2/3*L;
      case 2
         x(i)=0.446*L;
      case 3
         x(i)=0.5*L;
      end
      x(i)=x(i)-dL;
   case 3
      L=comp(i,4);
      d1=comp(i,5);
      d2=comp(i,6);
      CN(i)=2*[(d2/d)^2-(d1/d)^2];
      x(i)=pos+L/3*[1+(1-d1/d2)/(1-(d1/d2)^2)];
   case 4
      a=comp(i,4);
      b=comp(i,5);
      R=comp(i,6);
      S=comp(i,7);
      m=comp(i,8);
      n=comp(i,9);
      l=sqrt(S^2+(m+.5*b-.5*a)^2);
      p=R/S;
      if any(n==[3:4])
         Kfb=1+p/(1+p);
      elseif n==6
         Kfb=1+p/(1+p)/2;
      else
         error('Must have 3, 4 or 6 fins.')
      end
      CNf=4*n*(S/d)^2/(1+sqrt(1+(2*l/(a+b))^2));
      CN(i)=CNf*Kfb;
      x(i)=pos+m*(a+2*b)/(a+b)/3+(a+b-a*b/(a+b))/6;
   end
end

cp=evaldcm(x,CN);