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

Upload New File

parent 39a47aff
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/site/markdown/tutorials/matlab-examples/InverseLocation.m 0 → 100644
+ 183
− 0
View file @ fbac55fa
function InverseLocation()
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 org.orekit.rugged.errors.RuggedException;
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
% Initialize Orekit, assuming an orekit-data folder is in user home directory
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);
%% Building of the LOS
%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();
%% Configure the Datation
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);
%% Initialize the Line Sensor using LOS and Datation
import org.orekit.rugged.linesensor.LineSensor;
%With the LOS and the datation now defined , we can initialize a line sensor object in Rugged:
sensorName = 'mySensor';
lineSensor = LineSensor(sensorName, lineDatation, Vector3D.ZERO, lineOfSight);
%% Configure Frames (EME2000 and ITRF)
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);
%% Reading Quaternions
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);
%% Reading Potitions/Velocities
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);
%% Configure and Build Rugged
import org.orekit.rugged.api.*;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.utils.*;
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();
%% Geodetic Point creation and Inverse Location
%Point defined by its latitude, longitude and altitude
import org.orekit.bodies.GeodeticPoint;
latitude = deg2rad(34.27452);
longitude = deg2rad(-99.56804);
altitude = 0.0;
gp=GeodeticPoint(latitude, longitude, altitude);
% Search the sensor pixel seeing point
% the obteined result must be (1000,1000).
% the above geodetic point coordinates are obteined using DirecLocation.m to
% project the image point (1000,1000) on the ground.
minLine = 950;
maxLine = 1050;
sensorPixel = Rugged.inverseLocation(sensorName, gp, minLine, maxLine);
%% Test of the result and Print outputs
%we need to test if the sensor pixel is found in the prescribed lines
%otherwise the sensor pixel is empty
if (isempty(sensorPixel))
fprintf('Sensor Pixel is Empty: point cannot be seen between the prescribed line numbers\n')
else
fprintf('Sensor Pixel found : line = %5.3f, pixel = %5.3f \n', sensorPixel.getLineNumber(), sensorPixel.getPixelNumber())
%Find the date at which the sensor sees the ground point
dateLine = Rugged.dateLocation(sensorName, gp, minLine, maxLine);
fprintf(['Date at which the sensor sees the ground point :', char(dateLine.toString()),'\n'])
end
end
%% read Quaternions Function
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
%% read Position/Velocity Function
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% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment