Modeling Thrust Vectors in Orbiter

[N_Molson]

But I can look at the problem again the next days. I'm also interested in finding a physical/general solution to this problem.

Oops missed that one, let's say instead :

SetAnimation(anim_yaw_eng1, (1.0+1.0) / 4 + 0.5);

So that

1+1 = 2
2/4 = 0.5
0.5+0.5 = 1

right ?

Mmh no because if yaw = 1 the output is 0.75 only... mmhh... how to deal with that ?

@Face: any idea ?

What I need is a function that caps one value to 0.5 when the other is equal to 0.5. What am I missing ?

 

[N_Molson]

So I explained my case on a french forum specialized in maths, where you have bored high school math teachers hunting for problems to solve ? ... and they found the answer without sweating ! Again something in my math toolbox was missing : here the "extremum" or "maxima and minima" concept (https://en.wikipedia.org/wiki/Maxima_and_minima).

And we have that in C++ of course (maybe I should include math.h, it works on my end but maybe not on systems that don't have VS installed ?). So the animation state can be defined by :

anim_state = (1 + max(Ylvl, Rlvl)) / 2

where Ylvl is my yaw channel input and Rlvl is my roll channel input.

Edit : no, still not there yet.

 

[Face]

@Face: any idea ?

Sorry, I'm already out of the loop a bit. I don't really understand what you tried to achieve there. The operation there means that when you gradually change a positive yaw command to a positive roll command, your animation state will stay fully deflected most of the time.
I'd suggest to derive the animation states from the resulting thruster direction vector instead, this way you only have to do the math once, independent of the actual force vector calculation.

And we have that in C++ of course (maybe I should include math.h, it works on my end but maybe not on systems that don't have VS installed ?). So the animation state can be defined by :

anim_state = (1 + max(Ylvl, Rlvl)) / 2

The standard max() and min() functions in VS are macro definitions in WinDef.h . Should work everywhere.
However, there is the STL with its std::max() and std::min() functions, which are a different kind of beast. You'd need <algorithm> for that, and I am not sure if this means that the redistributable for the compiled platform toolset is necessary then.
Try to get to the declaration of your used function, if it gets you into WinDef.h with the simple macro, it should be as good as it gets.

 

[N_Molson]

You'd need <algorithm> for that, and I am not sure if this means that the redistributable for the compiled platform toolset is necessary then.

Yes, they mentionned that on the maths forum, say that in C++ I could write something like :

#include <algorithm>

using namespace std;

double anim(double x, double y)
{
    return clamp((1.+x+y)/2, 0. , 1.)
}

 

[N_Molson]

So I've got the code working in Orbiter ; the nozzles move consistently with the Pitch, Yaw, Roll inputs and there's no "priority" involved :

// Engines gimbal animations

    Plvl = GetControlSurfaceLevel(AIRCTRL_ELEVATOR);
    Ylvl = GetControlSurfaceLevel(AIRCTRL_RUDDER);
    Rlvl = GetControlSurfaceLevel(AIRCTRL_AILERON);

    Y_mix = (1 + min(1, max(-1, Ylvl + Rlvl))) / 2;
    R_mix = fabs(1 + min(1, max(-1, Ylvl + -Rlvl))) / 2;
  
    SetAnimation(anim_pitch_eng1, Plvl / 2 + 0.5);
    SetAnimation(anim_pitch_eng2, Plvl / 2 + 0.5);
  
    SetAnimation(anim_pitch_glow1, Plvl / 2 + 0.5);
    SetAnimation(anim_pitch_glow2, Plvl / 2 + 0.5);

    SetAnimation(anim_yaw_eng1, 1-Y_mix);
    SetAnimation(anim_yaw_eng2, 1-R_mix);

    SetAnimation(anim_yaw_glow1, 1 - R_mix); // engines and glow effects anims are "inverted", would be a good idea to fix that
    SetAnimation(anim_yaw_glow2, 1 - Y_mix);

I kept track of the 1-Y_mix ; 1-R_mix values with a DebugString, they both stay included between 0 and 1, so I'd say job done ! The thrust vectoring code posted above was working well on its own, I haven't touched it.

Tell me if you see ways of improving the code, or if I made some mistake. The "1-value" thing is messy, it could certainly simplified but it works and given my math skills I prefer to leave it that way. But if you know how to do it, please tell me.

I think we're good with both the thrust vectors and the anims. Give me the green light and I'll mark the thread as "solved for 2 engines".

Special thanks to the french guys of "l'île aux maths" forum (translation : "the maths island") : https://www.ilemaths.net/ especially "lematheuxmatou" (translation : the good-at-math tomcat) and "LittleFox" (no translation required ?). They nailed it with ease.

Edit : ah, not quite yet. I have to apply the same technique to the engine rotation matrixes, right now we have Y+R/2 and that's not good, it means that full yaw or roll is only 0.5 ! Should be doable.

 

[N_Molson]

So fixing the engines matrixes was easy. Control is of course now better and the vectors are still perfectly normal. It will be very easy to transpose it to the second stage. So here's the code for a pair of math-powered engines :

// we need separate matrixes to implement the roll channel (mixed with yaw)

    // engine 1 matrixes

    Ym_mix = min(1, max(-1, Ylvl + Rlvl));
    Rm_mix = min(1, max(-1, Ylvl + -Rlvl));

    MATRIX3 matrix_x_eng1 = {
        1, 0, 0,
        0, cos(GIMBAL_ANGLE*-Plvl),-sin(GIMBAL_ANGLE*-Plvl),
        0, sin(GIMBAL_ANGLE*-Plvl), cos(GIMBAL_ANGLE*-Plvl)
    };
    MATRIX3 matrix_y_eng1 = {
        cos(GIMBAL_ANGLE*Rm_mix),0,sin(GIMBAL_ANGLE*Rm_mix),
        0,1,0,
        -sin(GIMBAL_ANGLE*Rm_mix),0,cos(GIMBAL_ANGLE*Rm_mix),
    };

    MATRIX3 matrix_xy_eng1 = mul(matrix_x_eng1, matrix_y_eng1);

    VECTOR3 output_vec_eng1 = tmul(matrix_xy_eng1, YF77_dir);

    SetThrusterDir(th_YF77[0], output_vec_eng1);

    // engine 2 matrixes

    MATRIX3 matrix_x_eng2 = {
    1, 0, 0,
    0, cos(GIMBAL_ANGLE*-Plvl),-sin(GIMBAL_ANGLE*-Plvl),
    0, sin(GIMBAL_ANGLE*-Plvl), cos(GIMBAL_ANGLE*-Plvl)
    };
    MATRIX3 matrix_y_eng2 = {
        cos(GIMBAL_ANGLE*Ym_mix),0,sin(GIMBAL_ANGLE*Ym_mix),
        0,1,0,
        -sin(GIMBAL_ANGLE*Ym_mix),0,cos(GIMBAL_ANGLE*Ym_mix),
    };

    MATRIX3 matrix_xy_eng2 = mul(matrix_x_eng2, matrix_y_eng2);

    VECTOR3 output_vec_eng2 = tmul(matrix_xy_eng2, YF77_dir);

    SetThrusterDir(th_YF77[1], output_vec_eng2);

    sprintf(oapiDebugString(), "Y: %.3f ||R: %.3f", Ym_mix, Rm_mix);

Next step : apply the same idea to the "cross" 4 radial boosters configuration.

 

[N_Molson]

Applying the "extremum" concept to boosters matrixes was equally easy. Vectors are good. How powerful is that thing ! I begin to feel the true power of maths ?‍

The stack control is now very good. I use a full gimbal range of only +-5° ; giving we have 10 engines on liftoff it is more than enough. So I can start writing an ascent program

 

[N_Molson]

OK I'm having a slightly different problem with the Proton rocket.

The first stage features 6 engines that form an "hexagon" pattern like this :



Firstly, how can I define the "rotation" in 60 DEG steps for the 6 engines ? I tried with matrixes but without success. Because the "reference vector" is always _V(0,0,1) as all engines point straight down. So their direction is the same but their rotation is different.

Secondly, those engines gimbal only on 1 axis. What do you think is the most efficient configuration to achieve pitch, yaw and roll control ?

I tried that, but obviously it doesn't work as expected :

where vec_RD253dir is a constant _V(0,0,1)

// Engine 2 - 60 DEG rotation

    MATRIX3 matrix_rot_eng2 = {
    cos(1.0472),0,sin(1.0472),
    -sin(1.0472),0,cos(1.0472),
    0, 0, 1
    };

    MATRIX3 matrix_eng2 = {
    cos(GIMBAL_ANGLE*Plvl),0,sin(GIMBAL_ANGLE*Plvl),
    0,1,0,
    -sin(GIMBAL_ANGLE*Plvl),0,cos(GIMBAL_ANGLE*Plvl),
    };

    MATRIX3 matrix_total2 = mul(matrix_rot_eng2, matrix_eng2);
    VECTOR3 output_vec_eng2 = tmul(matrix_total2, vec_RD253dir);

    SetThrusterDir(th_RD253[1], output_vec_eng2);

 

[Urwumpe]

Didn't the Proton use a combination of differential thrust and TVC in the tangential plane? At least for roll control, the only real option you have is making the rotation axis radial (outward or inward from the core).

The rotation matrix is quite easy (in degrees):

sin(60°*n), cos(60°*n), 0,
cos(60° *n), -sin(60°*n), 0,
0, 0, 1

with n being the number of the engine.

Of course, you can also simplify this to just produce the rotation axis of a radial gimbal to _V(sin(n*DEG60), cos(n*DEG60), 0), with const double DEG60 = PI/3.0;
again, n is just the number of the engine, 0 should produce the gimbal for the top engine.

 

[N_Molson]

Didn't the Proton use a combination of differential thrust and TVC in the tangential plane? At least for roll control, the only real option you have is making the rotation axis radial (outward or inward from the core).

Now that's very interesting information, could you find a source about this ? That would be fantastic. Differential thrust is definitively a way to control pitch & yaw.

 

[Urwumpe]

Now that's very interesting information, could you find a source about this ? That would be fantastic. Differential thrust is definitively a way to control pitch & yaw.

Not that it would be needed... also it looks like it wasn't used:

"The engines in this stage can swivel tangentially up to 7.0° from the neutral position, providing full thrust vector control."

 

[N_Molson]

Yes we looked into Russian documents with MaxBuzz the exact number is "7 degrees 30 minutes" or 7.5°

Yes the tangential direction makes sense, else no roll control is possible at all.

Still have issues getting the right angles, modified my first matrix from the code above by yours, but I can't get that 60° rotation working.

 

[N_Molson]

I think I have the rotation working, please check my maths/code :

VECTOR3 vec_RD253dir = _V(0, 0, 1);

   // *1 = 60°, *2 = 120°....


    MATRIX3 matrix_rot_eng2 = {
    cos(1.0472*1), -sin(1.0472*1), 0,
    sin(1.0472*1), cos(1.0472*1), 0,
    0, 0, 1
    };

    MATRIX3 matrix_eng2y = {
    cos(GIMBAL_ANGLE*Plvl),0,-sin(GIMBAL_ANGLE*Plvl),
    0, 1, 0,
    sin(GIMBAL_ANGLE*Plvl),0,cos(GIMBAL_ANGLE*Plvl),
    };

    MATRIX3 matrix_total = mul(matrix_eng2y, matrix_rot_eng2);
    VECTOR3 output_vec_eng2 = tmul(matrix_total, vec_RD253dir);

    SetThrusterDir(th_RD253[1], output_vec_eng2);

 

[Urwumpe]

Why so complicated? Just let the compiler calculate those constants for you, thats more accurate, easier to read and avoids many typos.

 

[N_Molson]

Like this ?

    MATRIX3 matrix_rot_eng2 = {
    cos(PI / 3 * 1), -sin(PI / 3 * 1), 0,
    sin(PI / 3 * 1),  cos(PI / 3 * 1), 0,
    0, 0, 1,
    };

 

[Urwumpe]

For example, but since you need this six times for six engines, I would just define PI/3 as a const double in your program and name it in a way that YOU know its a 60° angle. Your compiler would turn this version now into a single constant anyway - while the later version would likely be produced as many different intermediate constants that have to be loaded, which is also a bit slower as just loading the constant into a floating point register once.

 

[N_Molson]

Yes, no need to ask the computer to calculate 6 times the very same value ?

 

[Urwumpe]

Yes, no need to ask the computer to calculate 6 times the very same value ?

More important: Tell the compiler what you really want to do, instead of relying that he guesses correctly about your intentions.

And remembering that you are human: I wouldn't know that 1.047198 is PI/3 after a few minutes. Not sure how your memory is there, but I can image you have better ideas what to think about than trying to remember why these digits are in that place and what you wanted to do with them...

 

[ChrisRowland]

I set up a ring of engines like this:

/// <summary>
/// set up a ring of raptors in the raptors[] array, assume the number is a multiple of 2 and generates them in pairs so
/// that any even number of thrusters will not change the mean centre of pressure.
/// </summary>
/// <param name="num">number of engines in the ring</param>
/// <param name="radius">radius of the ring m</param>
/// <param name="offset">angle by which the ring is offset rad</param>
/// <param name="initidx">the index of the first one in the raptors array</param>
void Superheavy::RaptorRing(int num, double radius, double offset, int initidx)
{
    double d = PI / num;
    for (int i = initidx; i < num + initidx; i++)
    {
        VECTOR3 pos = { radius * sin(i * d + offset), radius * cos(i * d + offset), RAP_Z };
        raptors[i] = CreateThruster(pos, _V(0,0,1), RAP_THRUST, hProp, RAP_ISP_SL);
        AddExhaust(raptors[i], 40, 1, pos, _V(0, 0, -1));
        i++;
        pos.x = -pos.x; pos.y = -pos.y;            // second raptor diagonally opposite
        raptors[i] = CreateThruster(pos, _V(0,0,1), RAP_THRUST, hProp, RAP_ISP_SL);
        AddExhaust(raptors[i], 40, 1, pos, _V(0, 0, -1));
    }
}

I'm offsetting the position, not rotating the engine. I've currently got 3 rings, an inner ring of 4, a middle ring of 8 and an outer ring of 16.

 

[Urwumpe]

Shouldn't d be defined as

double d = 2.0 * PI / num;

PI alone is just 180°, your angular steps would be only half as large as they should be for a full circle.