Skip to content
Snippets Groups Projects
Commit 39a47aff authored by Guylaine Prat's avatar Guylaine Prat
Browse files

Upload New File

parent 6687140e
No related branches found
No related tags found
No related merge requests found
function DirectLocation()
clear all
%Put 'orekit-data.zip' in the current directory, if not I get an error in line 55 %org.orekit.errors.OrekitException: aucune donne d'historique UTC-TAI n'a t charge
% These seems to work if pasted to prompt.
% javaaddpath 'C:\ ... enter your path here ...\MATLAB'
% javaaddpath 'C:\.. enter your path here ...\MATLAB\orekit-7.1.jar'
% javaaddpath 'C:\.. enter your path here ...\\MATLAB\rugged-1.0.jar'
% javaaddpath 'C:\.. enter your path here... \\MATLAB\commons-math3-3.6.1.jar'
%%
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.orekit.rugged.los.*;
import org.orekit.rugged.utils.*;
%% Configure Orekit. The file orekit-data.zip must be in current dir
import org.orekit.data.ZipJarCrawler
import org.orekit.time.AbsoluteDate
import org.orekit.time.TimeScalesFactory
DM=org.orekit.data.DataProvidersManager.getInstance();
crawler=org.orekit.data.ZipJarCrawler('orekit-data.zip');
DM.clearProviders();
DM.addProvider(crawler);
%%
%The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
rawDirs = ArrayList();
for i =1:2000
%20 field of view, 2000 pixels
rawDirs.add (Vector3D(0,i*deg2rad(20)/2000,1));
end
%The instrument is oriented 10 off nadir around the X-axis, we need to rotate the viewing
%direction to obtain the line of sight in the satellite frame
losBuilder = LOSBuilder(rawDirs);
losBuilder.addTransform(FixedRotation(ParameterType.FIXED, Vector3D.PLUS_I, deg2rad(10)));
lineOfSight = losBuilder.build();
%%
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.rugged.linesensor.LinearLineDatation;
%We use Orekit for handling time and dates, and Rugged for defining the datation model:
utc = TimeScalesFactory.getUTC();
absDate = AbsoluteDate('2009-12-11T16:59:30.0', utc);
lineDatation = LinearLineDatation(absDate, 1, 20);
%%
import org.orekit.rugged.linesensor.LineSensor;
%With the LOS and the datation now defined , we can initialize a line sensor object in Rugged:
lineSensor = LineSensor('mySensor', lineDatation, Vector3D.ZERO, lineOfSight);
%%
import org.orekit.frames.*;
import org.orekit.utils.*;
%In our application, we simply need to know the name of the frames we are working with. Positions and
%velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
%inertial frame.
eme2000 = FramesFactory.getEME2000();
simpleEOP = true; % we don't want to compute tiny tidal effects at millimeter level
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
%%
satelliteQList = ArrayList();
gps=TimeScalesFactory.getGPS();
addSatelliteQ(gps, satelliteQList, '2009-12-11T16:58:42.592937', -0.340236, 0.333952, -0.844012, -0.245684);
addSatelliteQ(gps, satelliteQList, '2009-12-11T16:59:06.592937', -0.354773, 0.329336, -0.837871, -0.252281);
addSatelliteQ(gps, satelliteQList, '2009-12-11T16:59:30.592937', -0.369237, 0.324612, -0.831445, -0.258824);
addSatelliteQ(gps, satelliteQList, '2009-12-11T16:59:54.592937', -0.3836, 0.319792, -0.824743, -0.265299);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:00:18.592937', -0.397834, 0.314883, -0.817777, -0.271695);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:00:42.592937', -0.411912, 0.309895, -0.810561, -0.278001);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:01:06.592937', -0.42581, 0.304838, -0.803111, -0.284206);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:01:30.592937', -0.439505, 0.299722, -0.795442, -0.290301);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:01:54.592937', -0.452976, 0.294556, -0.787571, -0.296279);
addSatelliteQ(gps, satelliteQList, '2009-12-11T17:02:18.592937', -0.466207, 0.28935, -0.779516, -0.302131);
%%
satellitePVList = ArrayList();
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T16:58:42.592937', -726361.466, -5411878.485, 4637549.599, -2463.635, -4447.634, -5576.736);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T16:59:04.192937', -779538.267, -5506500.533, 4515934.894, -2459.848, -4312.676, -5683.906);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T16:59:25.792937', -832615.368, -5598184.195, 4392036.13, -2454.395, -4175.564, -5788.201);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T16:59:47.392937', -885556.748, -5686883.696, 4265915.971, -2447.273, -4036.368, -5889.568);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:00:08.992937', -938326.32, -5772554.875, 4137638.207, -2438.478, -3895.166, -5987.957);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:00:30.592937', -990887.942, -5855155.21, 4007267.717, -2428.011, -3752.034, -6083.317);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:00:52.192937', -1043205.448, -5934643.836, 3874870.441, -2415.868, -3607.05, -6175.6);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:01:13.792937', -1095242.669, -6010981.571, 3740513.34, -2402.051, -3460.291, -6264.76);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:01:35.392937', -1146963.457, -6084130.93, 3604264.372, -2386.561, -3311.835, -6350.751);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:01:56.992937', -1198331.706, -6154056.146, 3466192.446, -2369.401, -3161.764, -6433.531);
addSatellitePV(gps, eme2000, itrf, satellitePVList, '2009-12-11T17:02:18.592937', -1249311.381, -6220723.191, 3326367.397, -2350.574, -3010.159, -6513.056);
%%
% import org.orekit.rugged.api.AlgorithmId;
% import org.orekit.rugged.api.BodyRotatingFrameId;
% import org.orekit.rugged.api.EllipsoidId;
% import org.orekit.rugged.api.InertialFrameId;
% import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.*;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.utils.*;
% import org.orekit.utils.CartesianDerivativesFilter;
rugged=RuggedBuilder();
rugged.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
rugged.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
rugged.setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0));
rugged.setTrajectory(InertialFrameId.EME2000,...
satellitePVList, 4, CartesianDerivativesFilter.USE_P,...
satelliteQList, 4, AngularDerivativesFilter.USE_R);
rugged.addLineSensor(lineSensor);
Rugged=rugged.build();
%%
import org.orekit.bodies.GeodeticPoint;
position = lineSensor.getPosition(); % This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
firstLineDate = lineSensor.getDate(0);
los = lineSensor.getLos(firstLineDate, 0);
upLeftPoint = Rugged.directLocation(firstLineDate, position, los);
fprintf('upper left point: Lat = %8.3f , Lon = %8.3f , h = %8.3f m\n',rad2deg(upLeftPoint.getLatitude()), rad2deg(upLeftPoint.getLongitude()),upLeftPoint.getAltitude() )
end
function satelliteQList=addSatelliteQ( gps, satelliteQList ,absDate,q0, q1, q2, q3)
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.time.AbsoluteDate;
import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
attitudeDate = AbsoluteDate(absDate, gps);
rotation = Rotation(q0, q1, q2, q3, true);
rotationRate = Vector3D.ZERO;
rotationAcceleration = Vector3D.ZERO;
pair = TimeStampedAngularCoordinates(attitudeDate, rotation, rotationRate, rotationAcceleration);
satelliteQList.add(pair);
end
function satellitePVList= addSatellitePV( gps, eme2000, itrf, satellitePVList, absDate, px, py, pz, vx, vy, vz)
import org.orekit.time.AbsoluteDate;
import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.orekit.utils.TimeStampedPVCoordinates;
import org.orekit.utils.PVCoordinates;
import org.orekit.frames.Transform;
ephemerisDate = AbsoluteDate(absDate, gps);
position = Vector3D(px, py, pz);
velocity = Vector3D(vx, vy, vz);
pvITRF = PVCoordinates(position, velocity);
transform = itrf.getTransformTo(eme2000, ephemerisDate);
pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add( TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment