Commit cbfa28e3 authored by Thomas Paulet's avatar Thomas Paulet
Browse files

Remove Abstract class for rediffused radiation pressure force model and

complete regular version of the Knocke force model. Field version is to
be implemented.
parent 62eb331f
/* Copyright 2002-2020 CS GROUP
* Licensed to CS GROUP (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.
*/
package org.orekit.forces.radiation;
import java.util.stream.Stream;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.orekit.forces.AbstractForceModel;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldEventDetector;
/** Base class for rediffused force models.
* This kind of force model takes in account either Earth albedo, Earth IR emission or both.
* @see KnockeRediffusedForceModel
* @author Thomas Paulet
*/
public abstract class AbstractRediffusedForceModel extends AbstractForceModel {
/** Constructor. */
protected AbstractRediffusedForceModel() {
// does nothing...
}
/** {@inheritDoc} */
@Override
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public Stream<EventDetector> getEventsDetectors() {
return Stream.of();
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
return Stream.of();
}
}
......@@ -16,15 +16,26 @@
*/
package org.orekit.forces.radiation;
import java.util.stream.Stream;
import org.hipparchus.Field;
import org.hipparchus.RealFieldElement;
import org.hipparchus.analysis.polynomials.PolynomialFunction;
import org.hipparchus.analysis.polynomials.PolynomialsUtils;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.forces.AbstractForceModel;
import org.orekit.frames.Frame;
import org.orekit.frames.Transform;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.FieldEventDetector;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
......@@ -36,9 +47,9 @@ import org.orekit.utils.ParameterDriver;
* This model is based on "EARTH RADIATION PRESSURE EFFECTS ON SATELLITES", 1988, by P. C. Knocke, J. C. Ries, and B. D. Tapley.
* </p> <p>
* This model represents the effects fo radiation pressure coming from the Earth.
* It considers Solar radiation which have been reflected by Earth (albedo) and Earth infrared emissions.
* The planet is considered as a sphere and is divided into segments.
* Each segment is considered as a plane which emits radiation according to Lambert's law.
* It considers Solar radiation which has been reflected by Earth (albedo) and Earth infrared emissions.
* The planet is considered as a sphere and is divided into elementary areas.
* Each elementary area is considered as a plane and emits radiation according to Lambert's law.
* The flux the satellite receives is then equal to the sum of the elementary fluxes coming from Earth.
* </p> <p>
* The radiative model of the satellite, and its ability to diffuse, reflect or absorb radiation is handled
......@@ -48,7 +59,10 @@ import org.orekit.utils.ParameterDriver;
*/
public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
public class KnockeRediffusedForceModel extends AbstractForceModel {
/** Earth rotation around Sun pulsation in rad/sec. */
private static final double EARTH_AROUND_SUN_PULSATION = 2 * FastMath.PI / Constants.JULIAN_YEAR;
/** Coefficient for solar irradiance computation. */
private static final double ES_COEFF = 4.5606E-6;
......@@ -89,51 +103,54 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
/** Sun model. */
private final ExtendedPVCoordinatesProvider sun;
/** Central body model. */
private final double equatorialRadius;
/** Spacecraft. */
private final RadiationSensitive spacecraft;
/** Angular equatorial latitude separation for surface element computation in rad. */
private final double dLatitude;
/** Angular resolution for emissivity and albedo computation in rad. */
private final double angularResolution;
/** Angular longitude separation for surface element computation in rad. */
private final double dLongitude;
/** Counter for debug. */
// TODO
private int nb_visibility;
/** Earth equatorial radius in m. */
private double equatorialRadius;
/** Constructor.
* @param sun Sun model
* @param equatorialRadius Central Body model
* @param spacecraft Spacecraft
* @param numberOfLatitudeSlices number of latitude slices for elementary area cut
* @param numberOfLongitudeSlices number of latitude slices for elementary area cut
* @param equatorialRadius the Earth equatorial radius in m
* @param angularResolution angular resolution in rad
*/
public KnockeRediffusedForceModel (final ExtendedPVCoordinatesProvider sun,
final double equatorialRadius,
final RadiationSensitive spacecraft,
final int numberOfLatitudeSlices,
final int numberOfLongitudeSlices) {
this.sun = sun;
this.equatorialRadius = equatorialRadius;
this.spacecraft = spacecraft;
this.dLatitude = FastMath.PI / numberOfLatitudeSlices;
this.dLongitude = 2 * FastMath.PI / numberOfLongitudeSlices;
this.nb_visibility = 0;
final double equatorialRadius,
final double angularResolution) {
this.sun = sun;
this.spacecraft = spacecraft;
this.equatorialRadius = equatorialRadius;
this.angularResolution = angularResolution;
}
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s,
final double[] parameters) {
public boolean dependsOnPositionOnly() {
return false;
}
/** {@inheritDoc} */
@Override
public Stream<EventDetector> getEventsDetectors() {
return Stream.of();
}
this.nb_visibility = 0;
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(final Field<T> field) {
return Stream.of();
}
// Initialize equatorial latitude it goes within [-PI / 2; PI / 2]
double currentLatitude = -FastMath.PI / 2 - dLatitude / 2;
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s,
final double[] parameters) {
// Get date
final AbsoluteDate date = s.getDate();
......@@ -144,65 +161,59 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
// Get satellite position
final Vector3D satellitePosition = s.getPVCoordinates().getPosition();
// Initialize rediffused flux
Vector3D rediffusedFlux = new Vector3D(0.0, 0.0, 0.0);
// Slice spherical Earth into elementary areas
while (currentLatitude + dLatitude <= FastMath.PI / 2) {
// Next latitude portion
currentLatitude += dLatitude;
// Initialize longitude, it goes within [0; 2 * PI]
double currentLongitude = -dLongitude / 2;
// Get spherical Earth model
final OneAxisEllipsoid earth = new OneAxisEllipsoid(equatorialRadius, 0.0, frame);
while (currentLongitude + dLongitude <= 2 * FastMath.PI) {
// Project satellite on Earth as geodetic point and vector
final GeodeticPoint satelliteAsGeodeticPoint = earth.transform(satellitePosition, frame, date);
final Vector3D projectedToGround = earth.projectToGround(satellitePosition, date, frame);
// Next longitude portion
currentLongitude += dLongitude;
// Get elementary vector east for Earth browsing using rotations
final Vector3D east = satelliteAsGeodeticPoint.getEast();
// Compute Earth element projected area vector
final Vector3D projectedAreaVector = computeProjectedAreaVector(satellitePosition,
currentLatitude, currentLongitude,
dLatitude, dLongitude);
// Initialize rediffused flux with elementary flux coming from the circular area around the projected satellite
final double centerArea = 2 * FastMath.PI * equatorialRadius * equatorialRadius *
(1 - FastMath.cos(angularResolution));
Vector3D rediffusedFlux = computeElementaryFlux(s, projectedToGround, earth, centerArea);
// Check if satellite sees the elementary area
if (projectedAreaVector.getNorm() > 0.0) {
// Sectorize the part of Earth which is seen by the satellite into crown sectors with constant angular resolution
for (double eastAxisOffset = angularResolution * 3 / 2;
eastAxisOffset < FastMath.asin(equatorialRadius / satellitePosition.getNorm());
eastAxisOffset = eastAxisOffset + angularResolution) {
// Compute Earth emissivity and albedo
final double e = computeEmissivity(date, currentLatitude);
// Build rotation transformations to get first crown elementary sector center
final Transform eastRotation = new Transform(date,
new Rotation(east, eastAxisOffset, RotationConvention.VECTOR_OPERATOR));
double a = 0.0;
// Get first elementary crown sector center
final Vector3D firstCrownSectorCenter = eastRotation.transformPosition(projectedToGround);
// Check if elementary area is in day light
final double sunElevation = FastMath.abs(Vector3D.angle(
computeSphericalElementaryVector(currentLatitude, currentLongitude),
sun.getPVCoordinates(date, frame).getPosition()));
// Browse the entire crown
for (double radialAxisOffset = angularResolution / 2;
radialAxisOffset < FastMath.PI * 2;
radialAxisOffset = radialAxisOffset + angularResolution) {
if (sunElevation < FastMath.PI / 2 ) {
// Build rotation transformations to get elementary area center
final Transform radialRotation = new Transform(date,
new Rotation(projectedToGround, radialAxisOffset, RotationConvention.VECTOR_OPERATOR));
// Elementary area is in day light
//a = computeAlbedo(date, currentLatitude);
a = 0.3;
}
// Get current elementary crown sector center
final Vector3D currentCenter = radialRotation.transformPosition(firstCrownSectorCenter);
// Compute solar flux
final double solarFlux = computeSolarFlux(date, frame);
// Compute current elementary crown sector area, it results of the integration of an elementary crown sector
// over the angular resolution
final double sectorArea = equatorialRadius * equatorialRadius *
2 * angularResolution * FastMath.sin(0.5 * angularResolution) * FastMath.sin(eastAxisOffset);
// Compute elementary area contribution to rediffused flux
final double albedoAndIR = (a * solarFlux * FastMath.cos(sunElevation) +
e * solarFlux / 4) / Constants.SPEED_OF_LIGHT;
final Vector3D elementaryFlux = projectedAreaVector.scalarMultiply(albedoAndIR);
// Add elementary contribution to total rediffused flux
rediffusedFlux = rediffusedFlux.add(elementaryFlux);
}
// Add current sector contribution to total rediffused flux
rediffusedFlux = rediffusedFlux.add(computeElementaryFlux(s, currentCenter, earth, sectorArea));
}
}
return spacecraft.radiationPressureAcceleration(date, frame, satellitePosition, s.getAttitude().getRotation(),
s.getMass(), rediffusedFlux, parameters);
}
/** {@inheritDoc} */
@Override
public <T extends RealFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
......@@ -211,6 +222,7 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
return null;
}
/** {@inheritDoc} */
@Override
public ParameterDriver[] getParametersDrivers() {
......@@ -231,8 +243,8 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
// Compute 1rst Legendre polynomial coeficient
final double A1 = C0 +
C1 * FastMath.cos(Constants.IERS96_EARTH_ANGULAR_VELOCITY * deltaT) +
C2 * FastMath.sin(Constants.IERS96_EARTH_ANGULAR_VELOCITY * deltaT);
C1 * FastMath.cos(EARTH_AROUND_SUN_PULSATION * deltaT) +
C2 * FastMath.sin(EARTH_AROUND_SUN_PULSATION * deltaT);
// Get 1rst and 2nd order Legendre polynomials
final PolynomialFunction firstLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(1);
......@@ -262,8 +274,8 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
// Compute 1rst Legendre polynomial coeficient
final double E1 = K0 +
K1 * FastMath.cos(Constants.IERS96_EARTH_ANGULAR_VELOCITY * deltaT) +
K2 * FastMath.sin(Constants.IERS96_EARTH_ANGULAR_VELOCITY * deltaT);
K1 * FastMath.cos(EARTH_AROUND_SUN_PULSATION * deltaT) +
K2 * FastMath.sin(EARTH_AROUND_SUN_PULSATION * deltaT);
// Get 1rst and 2nd order Legendre polynomials
final PolynomialFunction firstLegendrePolynomial = PolynomialsUtils.createLegendrePolynomial(1);
......@@ -293,75 +305,75 @@ public class KnockeRediffusedForceModel extends AbstractRediffusedForceModel {
return ES_COEFF * Constants.SPEED_OF_LIGHT / (earthSunDistance * earthSunDistance);
}
/** Compute Earth element projected area vector.
* @param satellitePosition satellite position with respect to Earth center
* @param phi the surface element center equatorial latitude in rad
* @param psi the surface element center longitude in rad
* @param dPhi the equatorial latitude angular separation in rad
* @param dPsi the longitude angular separation in rad
* @return the Earth element projected area vector
/** Compute elementary rediffused flux on satellite.
* @param state the current spacecraft state
* @param elementCenter the position of the considered area center
* @param earth the Earth model
* @param elementArea the area of the current element
* @return the rediffused flux from considered element on the spacecraft
*/
private Vector3D computeProjectedAreaVector(final Vector3D satellitePosition,
final double phi,
final double psi,
final double dPhi,
final double dPsi) {
private Vector3D computeElementaryFlux(final SpacecraftState state,
final Vector3D elementCenter,
final OneAxisEllipsoid earth,
final double elementArea) {
// Get Earth suface element center
final Vector3D elementaryAreaCenterPosition = computeSphericalElementaryVector(phi, psi).
scalarMultiply(equatorialRadius);
// Get satellite position
final Vector3D satellitePosition = state.getPVCoordinates().getPosition();
// Get current date
final AbsoluteDate date = state.getDate();
// Get frame
final Frame frame = state.getFrame();
// Get satellite viewing angle as seen from current elementary area
final double alpha = Vector3D.angle(elementCenter, satellitePosition);
// Compute Earth surface element center - satellite vector
final Vector3D r = satellitePosition.subtract(elementaryAreaCenterPosition);
// Check that satellite sees the current area
if (FastMath.abs(alpha) < FastMath.PI / 2) {
// Get satellite elevation angle
final double alpha = FastMath.abs(Vector3D.angle(elementaryAreaCenterPosition, r));
// Get current elementary area center latitude
final double currentLatitude = earth.transform(elementCenter, frame, date).getLatitude();
// Check if satellite sees the surface element
if (alpha < FastMath.PI / 2) {
// Compute Earth emissivity value
final double e = computeEmissivity(date, currentLatitude);
addVisibility();
// Initialize albedo
double a = 0.0;
// Compute surface element area
final double dA = equatorialRadius * equatorialRadius *
FastMath.cos(phi) * dPhi * dPsi;
// Check if elementary area is in day light
final double sunAngle = Vector3D.angle(elementCenter, sun.getPVCoordinates(date, frame).getPosition());
// Get Earth surface element center - satellite distance
if (FastMath.abs(sunAngle) < FastMath.PI / 2 ) {
// Elementary area is in day light, compute albedo value
a = computeAlbedo(date, currentLatitude);
}
// Compute solar flux
final double solarFlux = computeSolarFlux(date, frame);
// Compute elementary area contribution to rediffused flux
final double albedoAndIR = a * solarFlux * FastMath.cos(sunAngle) +
e * solarFlux / 4;
// Compute elementary area - satellite vector and distance
final Vector3D r = satellitePosition.subtract(elementCenter);
final double rNorm = r.getNorm();
// Compute Earth element projected area vector
return r.scalarMultiply(dA * FastMath.cos(alpha) / (FastMath.PI * rNorm * rNorm * rNorm));
// Compute attenuated projected elemetary area vector
final Vector3D projectedAreaVector = r.scalarMultiply(elementArea * FastMath.cos(alpha) /
(FastMath.PI * rNorm * rNorm * rNorm));
// Compute elementary radiation flux from current elementary area
return projectedAreaVector.scalarMultiply(albedoAndIR / Constants.SPEED_OF_LIGHT);
} else {
// Spacecraft does not see the elementary area
return new Vector3D(0.0, 0.0, 0.0);
}
}
/** Compute spherical elementary area vector.
* @param phi the surface element center equatorial latitude in rad
* @param psi the surface element center longitude in rad
* @return the spherical elementary area vector
*/
private Vector3D computeSphericalElementaryVector(final double phi,
final double psi) {
// Get angle sinuses and cosinuses
final double sinPhi = FastMath.sin(phi);
final double cosPhi = FastMath.cos(phi);
final double sinPsi = FastMath.sin(psi);
final double cosPsi = FastMath.cos(psi);
// Get Earth suface element center
final double x = cosPhi * cosPsi;
final double y = cosPhi * sinPsi;
final double z = sinPhi;
return new Vector3D(x, y, z);
}
//TODO
public void addVisibility() {
final int current = this.nb_visibility;
this.nb_visibility = current + 1;
}
}
Supports Markdown
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