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
Best Answer