Oddly specific question for help regarding a satellite in a shadow

[spaceiscoolright]

TL;DR

Vallado calculation to find shadow entry and exit points in an elliptical orbit.

Hello fellow Space People,

Currently trying to create a function in MATLAB to input an R and V vector and have it spit out all the orbital elements as well as shadow entry/exit locations. I have the Vallado "Fundamentals of Astrodynamics" book and I've used it religiously. I searched on here and only found help with circular orbits. I have those totally figured out, however, with the code that is free online, for some reason, it will not work with elliptical orbits based on my modeling software.

There are a couple of bugs in it out the gate, I'm pretty sure I fixed most of them. I model the orbits in GMAT and they don't match what the MATLAB code is giving me. I calculate circular orbits with different code (Because it's easier) with the same sun position vector and it works. I calculated the sun position vector the way Vallado instructs. So I'm not sure where the disconnect is. Does anyone have experience with this kind of thing? Most of the literature I find all comes back to Vallado so I'm not sure if there's a better way.

Thanks guys.

 

[BiGyElLoWhAt]

If the issue is with the equation, try using keplars law for orbits, area swept out over a time interval is constant.

 

[BiGyElLoWhAt]

So something like R*V*dt that should get you your elipse, fit it, and then calculate your axiis if that's what you're looking for.

 

[Filip Larsen]

You have implemented in Matlab equations found in Vallado for elliptical orbits, but then find that the trajectory itself doesn't match a same-parameter orbit modeled with GMAT, is that correct? Or is it only the sun transits that are wrong?

If you keep getting stuck on this perhaps you can list up which equations you are using from Vallado and the equivalent part of your Matlab code. As you no doubt know, without specific details its difficult for anyone here to provide you with other than general pointers that may or may not already be obvious to you.

 

[spaceiscoolright]

So something like R*V*dt that should get you your elipse, fit it, and then calculate your axiis if that's what you're looking for.

The issue is not knowing the orbit (if I'm understanding your explanation correctly). And I'm using Kepler to calculate the delta t but the problem is the input values for E are incorrect. I'm adding more specifics based on the question below you if a more detailed explanation is warranted.

 

[spaceiscoolright]

You have implemented in Matlab equations found in Vallado for elliptical orbits, but then find that the trajectory itself doesn't match a same-parameter orbit modeled with GMAT, is that correct? Or is it only the sun transits that are wrong?

If you keep getting stuck on this perhaps you can list up which equations you are using from Vallado and the equivalent part of your Matlab code. As you no doubt know, without specific details its difficult for anyone here to provide you with other than general pointers that may or may not already be obvious to you.

The sun transit is incorrect. The orbital elements match perfectly with what I'm calculating but the function is supposed to spit out shadow entry and exit points. The function itself is referenced in the book and downloadable from the internet, but it has some bugs in it out the gate. For example its output is the true anomaly, but it labels it as the Eccentric anomaly. And the code I'm using works like 75% of the time, just in some specific cases it is giving me more outputs than I can use computer logic to sift through.

IDK if the below code is too long for me to post on here for help but this is the part that is having issues.

P & Q are orientation vectors based off inclination, RAAN and Argument of Periapsis

RSun is the eath-sun distance in X,Y,Z and that has been validated

A0-A4 are coefficients in a quadratic equation that shows if the satellite is in shadow or not based on if it's positive or negative, so we find the zeros with "quartic" which obviously as a 4th order equation, it has 4 zeros.

In the quadratic equation "x" is equal to the cosine of the true anomaly, so by taking the inverse cosine of the zero points, you get 4 potential crossing points.

After that, the "for" loop first checks to see if beta*cos(true anomaly + zeta*sin(true anomaly) < 0 then goes back and recalculates the original A0-A4 quadratic equation for immediately before and immediately after the true anomaly that has passed the test thus far ("before" and "after" in the code). If the sign switches, in theory it's a transition point from shadow to sun or vice versa. Somehow, this will still output 4 angles even with the check in place. I've even put in two more logic checks and it still doesn't narrow it down. Based on GMAT results it does output two of the correct angles but I can't logically tell the computer which ones I want it to pick. I'm trying to make it work for any input, not just one particular solution I'm solving for.

Sorry if it's too long/complicated

    beta = dot(P_,RSun)/norm(RSun);
    zeta = dot(Q_,RSun)/norm(RSun);
    A0 = ((r_e/p)^4)*(ecc^4) - 2*((r_e/p)^2)*(zeta^2 - beta^2)*ecc^2 + (beta^2 + zeta^2)^2;
    A1 = 4*((r_e/p)^4)*ecc^3 - 4*((r_e/p)^2)*(zeta^2 - beta^2)*ecc;
    A2 = 6*((r_e/p)^4)*ecc^2 - 2*((r_e/p)^2)*(zeta^2 - beta^2) - 2*((r_e/p)^2)*(1-zeta^2)*ecc^2 + 2*(zeta^2 - beta^2)*(1 - zeta^2) - 4*(beta^2)*(zeta^2);
    A3 = 4*((r_e/p)^4)*ecc - 4*((r_e/p)^2)*(1-zeta^2)*ecc;
    A4 = (r_e/p)^4 - 2*((r_e/p)^2)*(1 - zeta^2) + (1 - zeta^2)^2;

    [r1r,~,r2r,~,r3r,~,r4r,~] = quartic( A0,A1,A2,A3,A4,'R' );
    nu(1) = acosd(r1r);
    nu(2) = acosd(r2r);
    nu(3) = acosd(r3r);
    nu(4) = acosd(r4r);
    k = 1;
    z=1;
    ii=1;
    for x = 1:length(nu)
        check(x) = beta*cosd(nu(x)) + zeta*sind(nu(x));
        if check(x) < 0
            nugood(k) = nu(x);
            %k = k + 1;
            disp('----------------------------------------------------------------------')
            disp(strcat(['Valid Eccentric Anomaly (beta*cos(nu) + zeta*sin(nu) < 0): ',num2str(nu(x)),' deg']))
            before = A0*cosd(nugood(k)-0.01)^4 + A1*cosd(nugood(k)-0.01)^3 + A2*cosd(nugood(k)-0.01)^2 + A3*cosd(nugood(k)-0.01) + A4;
            after = A0*cosd(nugood(k)+0.01)^4 + A1*cosd(nugood(k)+0.01)^3 + A2*cosd(nugood(k)+0.01)^2 + A3*cosd(nugood(k)+0.01) + A4;
            k = k + 1;
            if before > 0 && after < 0
                disp('Entering Shadow for This Eccentric Anomaly')
                Een(z) = nu(x);
                z=z+1;
            elseif before < 0 && after > 0
                disp('Exiting Shadow for This Eccentric Anomaly')
                Eex(ii) = nu(x);
                ii=ii+1;
            else
                disp('Error - whether entering or exiting is not clear')
            end
        end

[Code tags added by the Mentors]

 

[Filip Larsen]

From a quick look, your coefficients do appear to match those of equation 5-5. I was not able to find any information on the quartic function and and with no access to Matlab I'm not able to run it. I assume that the error you mention is that you for some problem values never get into the first or second branch for any of the ν-solutions. If that's the case I would recommend focus on validating that quartic really do what you think it does, if you haven't already.

Vallado also have a section "Numerical Shadow Analysis" where he mentions some methods useful for finding solutions that perhaps is another option for you if everything else fails.

Edit: just noticed that you have "check(x) < 0" for your "is zero" test. I would expect it to be something like "abs(check(x)) < eps" with eps being a suficiently small number, say 1E-6, depending on the accuracy of the solution coming out of quartic.

 

[spaceiscoolright]

Thanks for the recommendations! We were trying not to use the numerical shadow analysis ones because the slow speed of the iterative calculations and it's easier to incorporate J2 into the function. The quartic function is supposed to just find the zeros of an equation if you input the coefficients. I'll try and see if I can calculate it a different way and get a new result.

I like the new check parameters. Giving that a try as well!

 

[Filip Larsen]

You may also want to consider detecting the entry/exit type from the sign of the differential of the polynomial since you have all the coefficients anyway and it can found as a single evaluation of a 3rd order polynomial instead of two evaluations of a 4th order polynomial. Differentials at (or close to) zero probably means a (nearly) point-wise entry and exit.