Examples with NslSun gaiasky.util.coord.NslSun used on opensource projects

Example 1 with NslSun

use of gaiasky.util.coord.NslSun in project gaiasky by langurmonkey.

the class AttitudeConverter method getHeliotropicAnglesRates.

/**
 * Calculate the heliotropic angles and rates for a given attitude
 *
 * @param gt
 *            Time for the attitude
 * @param att
 *            attitude
 * @return
 */
public static HeliotropicAnglesRates getHeliotropicAnglesRates(long gt, Attitude att) {
    HeliotropicAnglesRates anglesAndRates = new HeliotropicAnglesRates();
    // k is a unit vector (in ICRS) towards the north ecliptic pole:
    Vector3d k = new Vector3d(0.0, -sinObliquity, cosObliquity);
    // s is a unit vector (in ICRS) towards the nominal sun:
    NslSun sun = new NslSun();
    sun.setTime(gt);
    double cosLSun = Math.cos(sun.getSolarLongitude());
    double sinLSun = Math.sin(sun.getSolarLongitude());
    Vector3d s = new Vector3d(cosLSun, sinLSun * cosObliquity, sinLSun * sinObliquity);
    // xyz[0], xyz[1], xyz[2] are unit vectors (in ICRS) along the SRS axes:
    att.getSrsAxes(xyz);
    // m = s x z is a non-unit vector (of length sinXi) normal to the plane
    // containing s and z:
    Vector3d m = new Vector3d(s);
    m.crs(xyz[2]);
    // compute solar aspect angle xi in range [0, pi]:
    double sinXi = m.len();
    double cosXi = s.dot(xyz[2]);
    anglesAndRates.setFirstAngle(Math.atan2(sinXi, cosXi));
    // NOTE: all subsequent computations fail if sinXi = 0
    // compute revolving phase angle nu in range [-pi, pi]:
    double sinXiCosNu = k.dot(m);
    double sinXiSinNu = k.dot(xyz[2]);
    anglesAndRates.setSecondAngle(Math.atan2(sinXiSinNu, sinXiCosNu));
    // compute spin phase Omega:
    double sinXiCosOmega = -m.dot(xyz[1]);
    double sinXiSinOmega = -m.dot(xyz[0]);
    anglesAndRates.setThirdAngle(Math.atan2(sinXiSinOmega, sinXiCosOmega));
    // inertial spin rate in ICRS:
    Vector3d spin = att.getSpinVectorInIcrs();
    // subtract motion of the nominal sun to get heliotropic spin rate:
    Vector3d spinHel = new Vector3d(spin);
    spinHel.add(k.scl(-sun.getSolarLongitudeDot()));
    // scalar products with s, z, and m are used to determine the angular
    // rates:
    double sSpinHel = s.dot(spinHel);
    double zSpinHel = xyz[2].dot(spinHel);
    double mSpinHel = m.dot(spinHel);
    // d(xi)/dt:
    anglesAndRates.setFirstRate(mSpinHel / sinXi);
    // d(nu)/dt:
    anglesAndRates.setSecondRate((sSpinHel - zSpinHel * cosXi) / (sinXi * sinXi));
    // d(Omega)/dt:
    anglesAndRates.setThirdRate((zSpinHel - sSpinHel * cosXi) / (sinXi * sinXi));
    return anglesAndRates;
}

Example 2 with NslSun

use of gaiasky.util.coord.NslSun in project gaiasky by langurmonkey.

the class AttitudeConverter method getQuaternionAndRate.

/**
 * Converts heliotropic angles and rates to an attitude quaternion and its
 * derivative
 *
 * @param gt
 *            GaiaTime
 * @param h
 *            heliotropic angles and their rates in [rad] and [rad/day]
 * @return
 * @return an array of two quaternions, q (the attitude quaternion) and qDot
 *         (the time derivative of q, per day)
 */
public static Quaterniond[] getQuaternionAndRate(long gt, HeliotropicAnglesRates h) {
    /**
     * SOME AXES NEED TO BE SWAPPED TO ALIGN WITH OUR REF SYS:
     * 	GLOBAL	GAIASANDBOX
     * 	Z -> Y
     * 	X -> Z
     * 	Y -> X
     */
    NslSun sun = new NslSun();
    sun.setTime(gt);
    double lSun = sun.getSolarLongitude();
    double lSunDot = sun.getSolarLongitudeDot();
    /**
     * Calculate the attitude quaternion *
     */
    Quaterniond q = new Quaterniond(Z_AXIS, OBLIQUITY_DEG);
    q.mul(new Quaterniond(Y_AXIS, Math.toDegrees(lSun)));
    q.mul(new Quaterniond(Z_AXIS, Math.toDegrees(h.getNu() - PI_HALF)));
    q.mul(new Quaterniond(X_AXIS, Math.toDegrees(PI_HALF - h.getXi())));
    q.mul(new Quaterniond(Y_AXIS, Math.toDegrees(h.getOmega())));
    /**
     * Calculate the time derivative of the attitude quaternion using (A.17)
     * in AGIS paper, based on the rates in the ICRS:
     */
    double sinLSun = Math.sin(lSun);
    double cosLSun = Math.cos(lSun);
    Vector3d zInSrs = aux1;
    zInSrs.set(Y_AXIS).mul(q);
    Vector3d sz = aux2;
    sz.set(sun.getSolarDirection(aux3)).crs(zInSrs).nor();
    double rateX = h.getNuDot() * cosLSun + h.getOmegaDot() * zInSrs.x + h.getXiDot() * sz.x;
    double rateY = -lSunDot * sinObliquity + h.getNuDot() * sinLSun * cosObliquity + h.getOmegaDot() * zInSrs.y + h.getXiDot() * sz.y;
    double rateZ = lSunDot * cosObliquity + h.getNuDot() * sinLSun * sinObliquity + h.getOmegaDot() * zInSrs.z + h.getXiDot() * sz.z;
    Quaterniond halfSpinInIcrs = new Quaterniond(0.5 * rateZ, 0.5 * rateX, 0.5 * rateY, 0.0);
    Quaterniond qDot = halfSpinInIcrs.mul(q);
    return new Quaterniond[] { q, qDot };
}

Example 3 with NslSun

use of gaiasky.util.coord.NslSun in project gaiasky by langurmonkey.

the class TransitionScanningLaw method initialize.

/**
 * Initialization mainly calculates the acceleration required for the
 * specified ramp
 */
public void initialize() {
    // constants:
    double xi = getXiRef();
    double sinXi = Math.sin(xi);
    double cosXi = Math.cos(xi);
    double Snom = NslUtil.calcSNom(xi, getTargetPrecessionRate());
    double rampDays = ramp.asDays();
    // make sure nuRef is nothing but 0.0 or PI
    setNuRef(getMode());
    // derivative of solar longitude [rad/day] at the end of the ramp
    long tEnd = getRefTime() + ramp.asNanoSecs();
    NslSun sunDir = getNominalSunVector();
    sunDir.setTime(tEnd);
    double lSunDotEnd = sunDir.getSolarLongitudeDot();
    // reference quantities for lSun, lSunDot, lSunDotDot:
    sunDir.setTime(getRefTime());
    double lSunDotRef = sunDir.getSolarLongitudeDot();
    double lSunDotDotRef = (lSunDotEnd - lSunDotRef) / rampDays;
    // calculate acceleration to reach nominal nuDot at tEnd:
    acc = 0.0;
    double sign = Math.cos(getNuRef());
    // to be sure. The feedback in the loop comes in via the acc variable.
    for (int i = 0; i < 10; i++) {
        // delta = nu (preceding) or nu - PI (following mode) at tEnd
        double delta = 0.5 * acc * rampDays * rampDays;
        double sinNu = sign * Math.sin(delta);
        acc = (Math.sqrt(Snom * Snom - 1.0 + sinNu * sinNu) + cosXi * sinNu) * lSunDotEnd / (sinXi * rampDays);
    }
    om0 = getOmegaRef();
    om1 = getTargetScanRate() * Nature.ARCSEC_TO_RAD * Nature.D_TO_S;
    om2 = -acc * cosXi / 2;
    om3 = -sign * acc * sinXi * lSunDotRef / 6;
    om4 = -sign * acc * sinXi * lSunDotDotRef / 8;
    setInitialized(true);
}

Example 4 with NslSun

use of gaiasky.util.coord.NslSun in project gaiasky by langurmonkey.

the class TransitionScanningLaw method getAttitudeNative.

/**
 * @see gaiasky.util.gaia.AnalyticalAttitudeDataServer#getAttitude(long)
 *
 * @param time
 *            - the time elapsed since the epoch of J2010 in ns (TCB)
 * @return attitude for the given time
 */
@Override
public Attitude getAttitudeNative(long time) {
    if (!isInitialized()) {
        initialize();
    }
    // t = time in [days] from tRef
    double t = (time - getRefTime()) * 1e-9 / Nature.D_TO_S;
    double tNorm = t / ramp.asDays();
    // tNorm must by in [-eps, |ramp|+eps] for a valid t
    if (tNorm < -eps || tNorm > 1.0 + eps) {
        throw new RuntimeException("TSL requested for time outside of ramp: t = " + t + " days, ramp = " + ramp.asDays() + " days");
    }
    // nu(t) is calculated for constant acceleration
    double nu = getNuRef() + 0.5 * acc * t * t;
    double nuDot = acc * t;
    // omega(t) is calculated to fourth order
    double omega = om0 + t * (om1 + t * (om2 + t * (om3 + t * om4)));
    double omegaDot = om1 + t * (2 * om2 + t * (3 * om3 + t * 4 * om4));
    NslSun sunDir = getNominalSunVector();
    // the sunDir uses the same reference epoch as this class, so we can pass the ns directly
    sunDir.setTime(time);
    // convert heliotropic angles and rates to quaternion and rate
    Quaterniond[] qr = AttitudeConverter.heliotropicToQuaternions(sunDir.getSolarLongitude(), getXiRef(), nu, omega, sunDir.getSolarLongitudeDot(), nuDot, omegaDot);
    return new ConcreteAttitude(time, qr[0], qr[1], false);
}

Example 5 with NslSun

use of gaiasky.util.coord.NslSun in project gaiasky by langurmonkey.

the class ConcreteAttitude method getHeliotropicAnglesRates.

/**
 */
public HeliotropicAnglesRates getHeliotropicAnglesRates() {
    HeliotropicAnglesRates anglesAndRates = new HeliotropicAnglesRates();
    // k is a unit vector (in ICRS) towards the north ecliptic pole:
    double obliquity = Coordinates.OBLIQUITY_RAD_J2000;
    double cosObliquity = Math.cos(obliquity);
    double sinObliquity = Math.sin(obliquity);
    Vector3d k = new Vector3d(0.0, -sinObliquity, cosObliquity);
    // s is a unit vector (in ICRS) towards the nominal sun:
    NslSun sun = new NslSun();
    double cosLSun = Math.cos(sun.getSolarLongitude());
    double sinLSun = Math.sin(sun.getSolarLongitude());
    Vector3d s = new Vector3d(cosLSun, sinLSun * cosObliquity, sinLSun * sinObliquity);
    // xyz[0], xyz[1], xyz[2] are unit vectors (in ICRS) along the SRS axes:
    getSrsAxes(xyz);
    // m = s x z is a non-unit vector (of length sinXi) normal to the plane
    // containing s and z:
    Vector3d m = aux;
    m.set(s);
    m.crs(xyz[2]);
    // compute solar aspect angle xi in range [0, pi]:
    double sinXi = m.len();
    double cosXi = s.dot(xyz[2]);
    anglesAndRates.setFirstAngle(Math.atan2(sinXi, cosXi));
    // NOTE: all subsequent computations fail if sinXi = 0
    // compute revolving phase angle nu in range [-pi, pi]:
    double sinXiCosNu = k.dot(m);
    double sinXiSinNu = k.dot(xyz[2]);
    anglesAndRates.setSecondAngle(Math.atan2(sinXiSinNu, sinXiCosNu));
    // compute spin phase Omega:
    double sinXiCosOmega = -m.dot(xyz[1]);
    double sinXiSinOmega = -m.dot(xyz[0]);
    anglesAndRates.setThirdAngle(Math.atan2(sinXiSinOmega, sinXiCosOmega));
    // inertial spin rate in ICRS:
    Vector3d spin = getSpinVectorInIcrs();
    // subtract motion of the nominal sun to get heliotropic spin rate:
    Vector3d spinHel = new Vector3d(spin);
    spinHel.scaleAdd(-sun.getSolarLongitudeDot(), k);
    // scalar products with s, z, and m are used to determine the angular
    // rates:
    double sSpinHel = s.dot(spinHel);
    double zSpinHel = xyz[2].dot(spinHel);
    double mSpinHel = m.dot(spinHel);
    // d(xi)/dt:
    anglesAndRates.setFirstRate(mSpinHel / sinXi);
    // d(nu)/dt:
    anglesAndRates.setSecondRate((sSpinHel - zSpinHel * cosXi) / (sinXi * sinXi));
    // d(Omega)/dt:
    anglesAndRates.setThirdRate((zSpinHel - sSpinHel * cosXi) / (sinXi * sinXi));
    return anglesAndRates;
}