Calculating Satellite Positions(Error msg: ??? Subscripted assignment dimension mismatch.)

Hi guys,

can anyone explain why i'm getting a subscripted assignment dimension mismatch when i try to write X,Y,Z sat positions (10 x 1) matrices into an array for each sat number? My codes are provided below. I cant seem to output my position values to the array created(satpos) .......

    GPSsecs=rem(dayssince/7,1)*7*24*60*60

     %EXTRACT INFORMATION FROM RINEX FILE - 

    rinexe('Rinex835.N', 'rinex_n.dat');
    eph = get_eph('rinex_n.dat');
    %Satellite numbers i am interested in from the rinex nav file
    svsprn = [1 3 4 13 16 19 20 23 24 27];

    for t=1:length(svsprn)
        icol(t) = find_eph(eph,svsprn(t),GPSsecs);

    end 

    for t=1:length(icol)
        ephsat(:,t)=eph(:,icol(t));
    end

    %Satellite positions are computed here

    %Satellite psudopositions for G1 G3 G4 G13 G16 G19 G20 G23 G24 G27 %(satellites)  extracted from rinex.835.O file

    PRsat = [24504451.001 24203490.140 24845984.603 22118003.935 22035927.410 25215698.785 20835960.556 20899808.315 24270950.683 22983577.273];

    %tgd values extracted manually from Rinex835.N file
    tgd= [-.372529029846E-08 -.419095158577E-08 -.605359673500E-08 -.116415321827E-07 -.977888703346E-08 -.139698386192E-07 -.698491930962E-08 -.214204192162E-07 -.139698386192E-08 -.465661287308E-08];

    %Constants  
    c=3e8;
    GM = 3.986005e14;  % earth's universal gravitational
                       % parameter m^3/s^2
    Omega_dot = 7.2921151467e-5; % earth rotation rate, rad/s
    F= -4.442807633e-10; %constant, sec/m^(1/2)

    % Loop starts here
     for t=1:length(PRsat)

    %Set up array for sat positions to be written into later on

    satpos=[1 0 0 0;3 0 0 0; 4 0 0 0;13 0 0 0;16 0 0 0;19 0 0 0;20 0 0 0;23 0 0 0;24 0 0 0;27 0 0 0];

%SATPOS Calculation of X,Y,Z coordinates at time t
  %       for given ephemeris eph

%extract all required variables 
      svprn=ephsat(1,t);
      toe=ephsat(18,t);
      A=(ephsat(4,t));
      deltaN=ephsat(5,t);
      M0=ephsat(3,t);
      ecc=ephsat(6,t);
      af0=ephsat(19,t);
      af1=ephsat(20,t);
      af2=ephsat(2,t);
      toc=ephsat(21,t);
      omega=ephsat(7,t);
      cuc=ephsat(8,t);
      cus=ephsat(9,t);
      crc=ephsat(10,t);
      crs=ephsat(11,t);
      i0=ephsat(12,t);
      idot=ephsat(13,t);
      cic=ephsat(14,t);
      cis=ephsat(15,t);
      Omega0=ephsat(16,t);
      Omegadot=ephsat(17,t);

% Finding satellite's position
  time= GPSsecs-(PRsat(t))/c; %GPS time  
  %extract toe
  tk = time-toe;
%clock correction
  if (tk> 302400)
  tk=tk-604800;
  elseif tk<-302400
  tk=tk+604800;
  end
%Extract semi-major axis and initial mean motion
  A = (sqrt(A))^2;
  n0 = sqrt(GM/A^3);
%find  mean motion n, mean anomaly Mk and Ek
  n = n0+deltaN;
  M_k = M0+n*tk;
%reduce mean anomaly to between 0 and 360deg
  M = rem(M_k+2*pi,2*pi);
  wrapTo2pi(M);
  E = M;
%iterate to compute eccentric anomaly
  for i = 1:10
     E_old = E;
     E = M+ecc*sin(E);
     dE = rem(E-E_old,2*pi);
     if abs(dE) < 1.e-12
        break;
     end
  end
%find true anomaly 
   E = rem(E+2*pi,2*pi);
   wrapTo2pi(E);

% %reducing ecc anomaly to between 0 and 360 deg

% E= acos((ecc+cos(v_k)/(1+(ecc*cos(v_k)))));
%relativistic correction term
  dtr= F* ecc*sqrt(A)*sin(E);
  dts= af0+(af1*(time-toc))+(af2*(time-toc)^2)+dtr;

%correction for single frequency
  dts=dts-tgd;

%correct gpstime
  time=time-dts;
%revise tk
  tk = time-toe;
  if (tk> 302400)
  tk=tk-604800;
  elseif tk<-302400
  tk=tk+604800;
  end
  M_k = M0+n*tk;
%reduce mean anomaly to between 0 and 360deg
  M = rem(M_k+2*pi,2*pi);
  wrapTo2pi(M);
  E = M;
%iterate to compute eccentric anomaly
  for i = 1:10
     E_old = E;
     E = M+ecc*sin(E);
     dE = rem(E-E_old,2*pi);
     wrapTo2pi(dE);
     if abs(dE) < 1.e-12
        break;
     end
  end
%compute true anomaly and arguement of latitude phi
  v_k = atan2(sqrt(1-ecc^2)*sin(E), cos(E)-ecc);
  wrapTo2pi(v_k);
  phi = v_k+omega;
  phi = rem(phi,2*pi);
  wrapTo2pi(phi);

%find 2nd harmonic pertubations
  deltauk=cuc*cos(2*phi)+cus*sin(2*phi);
  deltark=crc*cos(2*phi)+crs*sin(2*phi);
  deltaik=cic*cos(2*phi)+cis*sin(2*phi);

%correct arguement of latitude
  u = phi+ deltauk;
%correct of radius
  r = A*(1-ecc*cos(E)) + deltark ;

%correct inclination
  i = i0+ (idot*tk)  + deltaik;
  %angle between ascending node and greenwich meridian
  Omega = Omega0+(Omegadot-Omega_dot)*tk-Omega_dot*toe;
%reduced to between 0 and 360deg
  Omega = rem(Omega+2*pi,2*pi);
  wrapTo2pi(Omega);
%position in orbital space
   x_1 = ((cos(u'))*r);
   y_1 = (sin(u')*r);

%compute satellite  X,Y and Z coordinates

    x1 = x_1.*cos(Omega(t))-((y_1).*cos(i(t))*sin(Omega(t)));
    y1 =  x_1*sin(Omega(t))+(y_1*cos(i(t))*cos(Omega(t)));
    z1 = y_1*sin(i(t));

    satpos(2,1) =(x1); <-----stuck here!
    satpos(2,1) = y1;
    satpos(3,1) = z1;
    satpos(4,1) = dts ;

     end

[bym]

please format your code using {}

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 1

[Sean de Wolski]

dbstop if error

run your program, it'll take you the line of the error. When it gets there, evaluate the right-hand-side to see how big it is, and then the indices of the left hand side to see how big a block you're trying to stuff it into.

[Dorairaj]

thanks for the answer..i did that and eventually i got the satpos writing values for X Y and Z and dts when i did the following bit

% Loop starts here

for t=1:length(PRsat)

%Set up array for sat positions to be written into later on

satpos=[1 0 0 0;3 0 0 0; 4 0 0 0;13 0 0 0;16 0 0 0;19 0 0 0;20 0 0 0;23 0 0

.....

.......%compute satellite X,Y and Z coordinates

x1 = x_1.*cos(Omega(t))-((y_1).*cos(i(t))*sin(Omega(t)));

y1 = x_1*sin(Omega(t))+(y_1*cos(i(t))*cos(Omega(t)));

z1 = y_1*sin(i(t));

satpos(t,2) =(x1);

satpos(t,3) = y1;

satpos(t,4) = z1;

satpos(t,5) = dts ;

end