[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Orekit Users] orekit and Matlab

Interface to Matlab was an interesting idea so I tried that too, there seems to be some bugs in some matlab versions where the javaclasspath can not be set in a script, and needs to be stated directly at the prompt. This was very confusing, but when cutn' paste the javaclasspath commands to the prompt, integration between orekit and matlab seems to work out fine.

Another matlab example, of slave propagator, basically from tutorials:

% 
% /* Copyright 2002-2010 CS Communication & Systèmes
%  * Licensed to CS Communication & Systèmes (CS) under one or more
%  * contributor license agreements.  See the NOTICE file distributed with
%  * this work for additional information regarding copyright ownership.
%  * CS licenses this file to You under the Apache License, Version 2.0
%  * (the "License"); you may not use this file except in compliance with
%  * the License.  You may obtain a copy of the License at
%  *
%  *   http://www.apache.org/licenses/LICENSE-2.0
%  *
%  * Unless required by applicable law or agreed to in writing, software
%  * distributed under the License is distributed on an "AS IS" BASIS,
%  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
%  * See the License for the specific language governing permissions and
%  * limitations under the License.
%  */
% 
% Translated from SlaveMode.java to Matlab by Petrus Hyvönen 2011 as an
% example of how to access orekit from matlab
% the jars orekit-5.0.jar and commons-math-2.2.jar and orekit-data.zip is
% in current matlab dir.


______
% These seems to work if pasted to prompt.
javaaddpath 'C:\ ... enter your path here ...\MATLAB'
javaaddpath 'C:\.. enter your path here ...\MATLAB\orekit-5.0.jar'
javaaddpath 'C:\.. enter your path here ...\\MATLAB\commons-math-2.2.jar'
%javaaddpath 'C:\.. enter your path here ...\\MATLAB\orekit-data.zip' %not sure if this is needed


% Check your java pwd
foo=java.io.File('.');
foo.getCanonicalPath()

%% do the imports 
import org.orekit.errors.OrekitException
import org.orekit.frames.Frame
import org.orekit.frames.FramesFactory
import org.orekit.orbits.KeplerianOrbit
import org.orekit.orbits.Orbit
import org.orekit.propagation.SpacecraftState
import org.orekit.propagation.analytical.KeplerianPropagator
import org.orekit.data.DataProvidersManager
import org.orekit.data.ZipJarCrawler
import org.orekit.time.AbsoluteDate
import org.orekit.time.TimeScalesFactory

%% Configure Orekit. The file orekit-data.zip must be in current dir
DM=org.orekit.data.DataProvidersManager.getInstance()
crawler=org.orekit.data.ZipJarCrawler('orekit-data.zip')
DM.clearProviders()
DM.addProvider(crawler)

%% Initial orbit parameters
a = 24396159;    % semi major axis in meters
e = 0.72831215;  % eccentricity
i = (7.0)/180*pi; % inclination
omega = (180)/180*pi; % perigee argument
raan = (261)/180*pi; %right ascension of ascending node
lM = 0.0; % mean anomaly

%% Set inertial frame
inertialFrame = FramesFactory.getEME2000()

%% Initial date in UTC time scale
utc = TimeScalesFactory.getUTC();
initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc)

%% Setup orbit propagator
%gravitation coefficient
mu =  3.986004415e+14

%Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM, KeplerianOrbit.MEAN_ANOMALY, inertialFrame, initialDate, mu)

%Simple extrapolation with Keplerian motion
kepler = KeplerianPropagator(initialOrbit);

%Set the propagator to slave mode (could be omitted as it is the default mode)
kepler.setSlaveMode()

%% Setup propagation time
%Overall duration in seconds for extrapolation
duration = 180*60.0
            
%Stop date
finalDate =  AbsoluteDate(initialDate, duration, utc)
            
%Step duration in seconds
stepT = 30.0

%% Perform propagation
%Extrapolation loop
cpt = 1;
extrapDate = initialDate
while extrapDate.compareTo(finalDate) <= 0
    currentState = kepler.propagate(extrapDate);
    fprintf('step %d: time %s %s\n', cpt, char(currentState.getDate()), char(currentState.getOrbit()))
    coord=currentState.getPVCoordinates.getPosition;
    P(:,cpt)=[coord.getX coord.getY coord.getZ]';
    extrapDate = AbsoluteDate(extrapDate, stepT, utc);
    cpt=cpt+1;
end
%%
figure;
plot3(P(1,:),P(2,:),P(3,:));

____

On Thu, May 19, 2011 at 1:33 PM, Brian H <brianh4321@hotmail.com> wrote:
Concerning the matlab interface(?)

While its somewhat crude, here is my early test script of the first tutorial. If I remember correctly you have to set the orekit data path using the System call below. There was also a problem with matlab's imported version of the commons-math library being somewhat out of date. This was fixed by downloading the current version and adding that directory at the top of your classpath.txt located in the MATLAB directory. Anyway, here is a copy of the original mfile. It was my first attempt at the matlab-java interface, pardon it's crudeness it was meant as an exercise.

import java.lang.*
import java.lang.Enum
import java.text.*
import java.util.*

import org.apache.commons.*
import org.apache.commons.math.*
import org.apache.commons.math.geometry.*
import org.apache.commons.math.util.*
import org.apache.commons.math.util.MathUtils.*
import org.orekit.bodies.*
import org.orekit.bodies.BodyShape.*
import org.orekit.bodies.GeodeticPoint.*
import org.orekit.bodies.OneAxisEllipsoid.*
import org.orekit.errors.OrekitException.*
import org.orekit.frames.*
import org.orekit.frames.Frame.*
import org.orekit.frames.FramesFactory.*
import org.orekit.frames.TopocentricFrame.*
import org.orekit.frames.LocalOrbitalFrame.*
import org.orekit.orbits.*
import org.orekit.orbits.CartesianOrbit.*
import org.orekit.orbits.Orbit.*
import org.orekit.propagation.*
import org.orekit.propagation.Propagator.*
import org.orekit.propagation.analytical.*
import org.orekit.propagation.analytical.KeplerianPropagator.*
import org.orekit.time.*
import org.orekit.time.AbsoluteDate.*
import org.orekit.time.TimeScale.*
import org.orekit.time.TimeScalesFactory.*
import org.orekit.utils.*
import org.orekit.utils.PVCoordinates.*
import org.orekit.data.*
import org.orekit.data.DataProvidersManager.*
import LOFTypeWrapper
System.setProperty('orekit.data.path','orekit-data')


symbols = DecimalFormatSymbols(Locale.US);
d3 = DecimalFormat('0.000', symbols);
utc = TimeScalesFactory.getUTC();
initialDate = AbsoluteDate(2008, 10, 01, 0, 0, 00.000, utc);
mu =  3.986004415e+14; % gravitation coefficient
inertialFrame = FramesFactory.getEME2000();
posisat = Vector3D(-6142438.668, 3492467.560, -25767.25680);
velosat = Vector3D(505.8479685, 942.7809215, 7435.922231);
pvsat = PVCoordinates(posisat, velosat);
initialOrbit = CartesianOrbit(pvsat, inertialFrame, initialDate, mu);
kepler = KeplerianPropagator(initialOrbit);
lof = LocalOrbitalFrame(inertialFrame, LOFTypeWrapper.getType(1), kepler, 'LOF');

ae =  6378137.0; % equatorial radius in meter
f  =  1.0 / 298.257223563; % flattening
ITRF2005 = FramesFactory.getITRF2005(true); % terrestrial frame at an arbitrary date

earth = OneAxisEllipsoid(ae, f, ITRF2005);

% Station
longitude = Math.toRadians(45.);
latitude  = Math.toRadians(25.);
altitude  = 0.;
station = GeodeticPoint(latitude, longitude, altitude);
staF = TopocentricFrame(earth, station, 'station1');
System.out.println('          time           doppler (m/s)');
finalDate = AbsoluteDate(initialDate, 6000, utc);
extrapDate = initialDate;
while (extrapDate.compareTo(finalDate) <= 0) 
   
    % We can simply get the position and velocity of station in LOF frame at any time
    pv = staF.getTransformTo(lof, extrapDate).transformPVCoordinates(PVCoordinates.ZERO);
   
    % And then calculate the doppler signal
    doppler = Vector3D.dotProduct(pv.getPosition(), pv.getVelocity()) / pv.getPosition().getNorm();
    abc= d3.format(doppler);
    def= String('     ');
    ghi= extrapDate.toString();
    jkl= ghi.concat(def.concat(abc));
    System.out.println(jkl);
   
    extrapDate = AbsoluteDate(extrapDate, 600, utc);
   
end
 


> Date: Thu, 19 May 2011 09:46:30 +0200
> From: luc.maisonobe@c-s.fr
> To: orekit-users@orekit.org
> Subject: Re: [Orekit Users] orekit and Jython

>
> beowulf zhang <beowulf.zhang@gmail.com> a écrit :
>
> > Very good! Thank you...
> >
> > I konw the m language of Matlab support java programming now, but I
> > was failure to test the orekit examples. I want to know whether anyone
> > are interesting testing the orekit examples in Matlab or not.
>
> If you want to test these examples, I would be interested to know if
> they work.
> We could set up some wiki pages (yes, we know, the wiki is quite empty
> for now) showing how to interface with several scripting languages.
>
> Luc
>
> >
> >
> > 2011/5/17, Petrus Hyvönen <petrus.hyvonen@gmail.com>:
> >> Hi all,
> >>
> >> I was playing around with the idea of using Python for mission design, and
> >> found Jython and the orekit library. As new to both Python and Java, I had
> >> some initial troubles of getting things to work, primary the access to
> >> orekit-data, what was current dir and javapaths etc. I am now using Eclipse
> >> with pydev extension, and can add the orekit jars and orekit-data to the
> >> eclipse project.
> >>
> >> As an exercise I took the liberty to translate a few of the orekit examples
> >> from java to jython, and copy them below as it could be of interest for
> >> others. I have been exploring scipy before (based on regular python) and it
> >> contains very good plotting and array routines, which is not available in
> >> jython. So next step would be to find a suitable java plotting and arrays
> >> library to be able to represent the results.
> >>
> >> Regards
> >> /Petrus
> >>
> >> First, a jython translation of the SlaveMode.java:
> >> __________________________________________
> >>
> >> # -*- coding: utf-8 -*-
> >> '''
> >> /* Copyright 2002-2010 CS Communication & Systèmes
> >> * Licensed to CS Communication & Systèmes (CS) under one or more
> >> * contributor license agreements. See the NOTICE file distributed with
> >> * this work for additional information regarding copyright ownership.
> >> * CS licenses this file to You under the Apache License, Version 2.0
> >> * (the "License"); you may not use this file except in compliance with
> >> * the License. You may obtain a copy of the License at
> >> *
> >> * http://www.apache.org/licenses/LICENSE-2.0
> >> *
> >> * Unless required by applicable law or agreed to in writing, software
> >> * distributed under the License is distributed on an "AS IS" BASIS,
> >> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
> >> * See the License for the specific language governing permissions and
> >> * limitations under the License.
> >> */
> >>
> >> Translated from SlaveMode.java to jython by Petrus Hyvönen 2011-05-04
> >> '''
> >>
> >> # orekit.jar,orekit-data.zip and commons-maths in CLASSPATH through eclipse
> >> project
> >> import java, os
> >>
> >> from org.orekit.errors import OrekitException
> >> from org.orekit.frames import Frame
> >> from org.orekit.frames import FramesFactory
> >> from org.orekit.orbits import KeplerianOrbit
> >> from org.orekit.orbits import Orbit
> >> from org.orekit.propagation import SpacecraftState
> >> from org.orekit.propagation.analytical import KeplerianPropagator
> >> from org.orekit.data import DataProvidersManager
> >> from org.orekit.data import ZipJarCrawler
> >> from org.orekit.time import AbsoluteDate
> >> from org.orekit.time import TimeScalesFactory
> >>
> >> from math import radians
> >>
> >> # Configure Orekit. The file orekit-data.zip must be in current dir
> >> DM = DataProvidersManager.getInstance()
> >> crawler=ZipJarCrawler("orekit-data.zip")
> >> DM.clearProviders()
> >> DM.addProvider(crawler)
> >>
> >> #Initial orbit parameters
> >> a = 24396159 # semi major axis in meters
> >> e = 0.72831215 # eccentricity
> >> i = radians(7.0)# inclination
> >> omega = radians(180) # perigee argument
> >> raan = radians(261) #right ascension of ascending node
> >> lM = 0.0 # mean anomaly
> >>
> >> #Inertial frame
> >> inertialFrame = FramesFactory.getEME2000()
> >>
> >> #Initial date in UTC time scale
> >> utc = TimeScalesFactory.getUTC();
> >> initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc)
> >>
> >> #gravitation coefficient
> >> mu = 3.986004415e+14
> >>
> >> #Orbit construction as Keplerian
> >> initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM,
> >> KeplerianOrbit.MEAN_ANOMALY,
> >> inertialFrame, initialDate, mu)
> >>
> >> #Simple extrapolation with Keplerian motion
> >> kepler = KeplerianPropagator(initialOrbit)
> >>
> >> #Set the propagator to slave mode (could be omitted as it is the default
> >> mode)
> >> kepler.setSlaveMode()
> >>
> >> #Overall duration in seconds for extrapolation
> >> duration = 90*60.0
> >>
> >> #Stop date
> >> finalDate = AbsoluteDate(initialDate, duration, utc)
> >>
> >> #Step duration in seconds
> >> stepT = 30.0
> >>
> >> #Extrapolation loop
> >> cpt = 1
> >> extrapDate = initialDate
> >> while extrapDate.compareTo(finalDate) <= 0:
> >> currentState = kepler.propagate(extrapDate)
> >> print "step %d: time %s %s" % (cpt, currentState.getDate(),
> >> currentState.getOrbit())
> >> extrapDate = AbsoluteDate(extrapDate, stepT, utc)
> >> cpt=cpt+1
> >>
> >>
> >> ____________________________
> >>
> >> An interesting example of jython inheritance of java class:
> >> _____________________________
> >>
> >>
> >> # -*- coding: utf-8 -*-
> >> '''
> >> /* Copyright 2002-2010 CS Communication & Syst?mes
> >> * Licensed to CS Communication & Syst?mes (CS) under one or more
> >> * contributor license agreements. See the NOTICE file distributed with
> >> * this work for additional information regarding copyright ownership.
> >> * CS licenses this file to You under the Apache License, Version 2.0
> >> * (the "License"); you may not use this file except in compliance with
> >> * the License. You may obtain a copy of the License at
> >> *
> >> * http://www.apache.org/licenses/LICENSE-2.0
> >> *
> >> * Unless required by applicable law or agreed to in writing, software
> >> * distributed under the License is distributed on an "AS IS" BASIS,
> >> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
> >> * See the License for the specific language governing permissions and
> >> * limitations under the License.
> >> */
> >>
> >> Translated from SlaveMode.java to jython by Petrus Hyv?nen 2011-05-04
> >> '''
> >> # orekit and common maths in CLASSPATH through eclipse project
> >> import java, os
> >>
> >> from org.orekit.data import DataProvidersManager
> >> from org.orekit.data import ZipJarCrawler
> >> from org.apache.commons.math.geometry import Vector3D
> >> from org.orekit.bodies import BodyShape
> >> from org.orekit.bodies import GeodeticPoint
> >> from org.orekit.bodies import OneAxisEllipsoid
> >> from org.orekit.errors import OrekitException;
> >> from org.orekit.frames import Frame
> >> from org.orekit.frames import FramesFactory
> >> from org.orekit.frames import TopocentricFrame
> >> from org.orekit.orbits import KeplerianOrbit
> >> from org.orekit.orbits import Orbit
> >> from org.orekit.propagation import Propagator
> >> from org.orekit.propagation import SpacecraftState
> >> from org.orekit.propagation.analytical import KeplerianPropagator
> >> from org.orekit.propagation.events import ElevationDetector
> >> from org.orekit.propagation.events import EventDetector
> >> from org.orekit.time import AbsoluteDate
> >> from org.orekit.time import TimeScalesFactory
> >> from org.orekit.utils import PVCoordinates
> >>
> >> from math import degrees, radians, pi
> >>
> >> # Configure Orekit
> >> DM = DataProvidersManager.getInstance()
> >> crawler=ZipJarCrawler("orekit-data.zip")
> >> DM.clearProviders()
> >> DM.addProvider(crawler)
> >>
> >> # Initial state definition: date, orbit
> >> initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000,
> >> TimeScalesFactory.getUTC())
> >> mu = 3.986004415e+14
> >> inertialFrame = FramesFactory.getEME2000() # inertial frame for orbit
> >> definition
> >> position = Vector3D(-6142438.668, 3492467.560, -25767.25680)
> >> velocity = Vector3D(505.8479685, 942.7809215, 7435.922231)
> >> pvCoordinates = PVCoordinates(position, velocity)
> >> initialOrbit = KeplerianOrbit(pvCoordinates, inertialFrame, initialDate, mu)
> >>
> >> # Propagator : consider a simple keplerian motion (could be more elaborate)
> >> kepler = KeplerianPropagator(initialOrbit)
> >>
> >> #Earth and frame
> >> ae = 6378137.0 # // equatorial radius in meter
> >> f = 1.0 / 298.257223563 #; // flattening
> >> ITRF2005 = FramesFactory.getITRF2005() #; // terrestrial frame at an
> >> arbitrary date
> >> earth = OneAxisEllipsoid(ae, f, ITRF2005)
> >>
> >> # Station
> >> longitude = radians(45.0)
> >> latitude = radians(25.0)
> >> altitude = 0.0
> >> station1 = GeodeticPoint(latitude, longitude, altitude)
> >> sta1Frame = TopocentricFrame(earth, station1, "station1")
> >>
> >> # Event definition
> >> maxcheck = 1.0
> >> elevation = radians(5.0)
> >>
> >>
> >> class VisibilityDetector(ElevationDetector):
> >> # Class for handling the eventOccured java. Example of subclassing
> >> # a java class in jython
> >> def __init__(self, maxCheck, elevation, topo):
> >> ElevationDetector.__init__(self,maxCheck, elevation, topo)
> >>
> >> def eventOccurred(self, s, increasing):
> >> if (increasing):
> >> print "Visibility on", self.topocentricFrame.getName(),"begins
> >> at" , s.getDate()
> >> else:
> >> print "Visibility on", self.topocentricFrame.getName(), "ends
> >> at" , s.getDate()
> >> return self.CONTINUE
> >>
> >> sta1Visi = VisibilityDetector(maxcheck, elevation, sta1Frame)
> >>
> >> #Add event to be detected
> >> kepler.addEventDetector(sta1Visi)
> >>
> >> #Propagate from the initial date to the first raising or for the fixed
> >> duration
> >> finalState = kepler.propagate(initialDate.shiftedBy(1500.0))
> >>
> >> print "Final state : " , finalState.getDate().durationFrom(initialDate)
> >>
> >> _________________________________
> >>
> >>
> >>
> >
> >
>
>
>
> ----------------------------------------------------------------
> This message was sent using IMP, the Internet Messaging Program.
>
>



--
_____________________________________________
Petrus Hyvönen, Uppsala, Sweden
Mobile Phone/SMS:+46 73 803 19 00