FGI-GSRx icon indicating copy to clipboard operation
FGI-GSRx copied to clipboard

B1I GEO satellite velocity ?

Open xixi88888888 opened this issue 1 year ago • 3 comments

Hello, author. Thank you very much for sharing your code; it's really well-written. However, I still have a few questions. Could you please clarify why there was no coordinate axis rotation correction applied when calculating the velocity of the BeiDou B1I GEO satellite? Is this a mistake? I look forward to your response. Thank you!

xixi88888888 avatar Dec 14 '24 07:12 xixi88888888

Dear Xixi,

Thank you for your query. We do not usually use BeiDou B1I GEO satellites in PVT computation, and hence this part is kept open to be addressed. Please feel free to offer your own modifications to include this correction. Can you share your update/fix that can be incorporated in a newer version of FGI-GSRx?

Thank you.

Zahid

zahidulbhuiyan avatar Jan 07 '25 13:01 zahidulbhuiyan

First, I believe the original velocity calculation formula for BeiDou MEO and IGSO satellites is incorrect, and I have made the following modifications. Please refer to see if it is correct:

% satVelocity(1) = -Ldot*(dY) - (dtemp3 + dtemp2)sL+ XpdotcL ; %The velocity calculation is incorrect.The correction is as follows: satVelocity(1) = -Ldot*(dY) - (dtemp3 - dtemp2)sL+ XpdotcL ; satVelocity(2) = Ldot*(dX) + (dtemp3 - dtemp2)cL+ XpdotsL ; satVelocity(3) = dtempIdot + YpdotsI; satVelocity(4) = dDeltaFreq;

Next, here is the code I wrote for calculating the velocity of BeiDou GEO satellites. Please review and provide feedback!

satVelocity_copy(1:4)=satVelocity(1:4); if(geoSV) % For GEO satellite position computation: minus5degreeInRadian = ((gpsPi*(-5))/180); R_z = [ cos(WGS84oetk) sin(WGS84oetk) 0; -sin(WGS84oetk) cos(WGS84oetk) 0; 0 0 1]; R_x = [1 0 0; 0 cos(minus5degreeInRadian) sin(minus5degreeInRadian); 0 -sin(minus5degreeInRadian) cos(minus5degreeInRadian)]; satPositions = (R_z*R_x)*satPositions';

% For GEO Velocity  computation:
minus5degreeInRadian = ((gpsPi*(-5))/180);
R_z_x=(R_z*R_x);
R_z_x_dot=[-WGS84oe*sin(WGS84oe*tk)  WGS84oe*cos(WGS84oe*tk)*cos(minus5degreeInRadian)  WGS84oe*cos(WGS84oe*tk)*sin(minus5degreeInRadian);
           -WGS84oe*cos(WGS84oe*tk) -WGS84oe*sin(WGS84oe*tk)*cos(minus5degreeInRadian) -WGS84oe*sin(WGS84oe*tk)*sin(minus5degreeInRadian);
                0                               0                                           0
          ];
satVelocity(1:3) = R_z_x_dot*satPositions_copy'+ R_z_x*satVelocity_copy(1:3)';

end

xixi88888888 avatar Jan 17 '25 09:01 xixi88888888

Hi, Xixi!

I checked the modifications but came to a conclusion that the MEO satellite velocity is correct as it is in the current implementation of the receiver.

The addition of GEO velocity appears to be correct except for a change to the longitude of ascending node rate Ldot. I would make the following changes to the mulSatPos codes:

  1. The longitude of ascending node rate Ldot should be in an inertial coordinate system that aligns with the BDC ECEF-system at ephemeris reference time. It would then match how the GEO satellite position is computed.
if(geoSV)
    L  = dOmega0 + tk * dOmegaDot;
    L = L - WGS84oe * dToe;
    Ldot = dOmegaDot;                    % Earth's angular velocity is removed for GEO
else
    L  = dOmega0 + tk * (dOmegaDot - WGS84oe);
    L = L - WGS84oe * dToe;
    Ldot = dOmegaDot - WGS84oe;
end
  1. Going from the inertial system to BDCS as you have suggested:
if(geoSV)
    minus5degreeInRadian = -5*gpsPi/180;
    R_x =     [1    0                           0;
               0    cos(minus5degreeInRadian)   sin(minus5degreeInRadian);
               0    -sin(minus5degreeInRadian)  cos(minus5degreeInRadian)];
    R_z =     [cos(WGS84oe*tk)  sin(WGS84oe*tk) 0;
               -sin(WGS84oe*tk) cos(WGS84oe*tk) 0;
               0                0               1];
    R_z_dot = [-WGS84oe*sin(WGS84oe*tk)    WGS84oe*cos(WGS84oe*tk)  0;
               -WGS84oe*cos(WGS84oe*tk)    -WGS84oe*sin(WGS84oe*tk) 0;
               0                           0                        0];

    % Product matrices
    R_z_x = R_z*R_x;
    R_z_x_dot = R_z_dot*R_x;

    % Compute satellite position and velocity
    satVelocity(1:3) = R_z_x_dot*satPositions'+ R_z_x*satVelocity(1:3)';
    satPositions = R_z_x*satPositions';
end

Best regards, Into

paakkoi avatar Jan 27 '25 09:01 paakkoi