Modelling a rocket

[Glom]

I'm writing a program to model the launch of a rocket. Here's the particular function so far.

//1. Load the position and speed of the rocket on the pad
	double lslat = ( latdegbox->value() + ( latminbox->value() + latsecbox->value()/60) / 60 ) * M_PI / 180;
	if ( southbutton->set() ) {
		lslat *= -1;
	}
	double lslong = ( longdegbox->value() + ( longminbox->value() + longsecbox->value()/60) / 60 ) * M_PI / 180;
	if ( westbutton->set() ) {
		lslong *= -1;
	}
	double lselev = elevbox->value() * 0.3048;
	double t0 = ( utchourbox->value() + ( utcminbox->value() + utcsecbox->value()/60) / 60 ) * M_PI / 180;
	double x = (lselev + georad) * cos (lslat) * cos ( lslong + t0);
	double y = (lselev + georad) * cos (lslat) * sin ( lslong + t0);
	double z = (lselev + georad) * sin (lslat);
	vec r(x,y,z);
	x = - ( lselev + georad ) * cos (lslat) * georate * sin (lslong);
	y = ( lselev + georad ) * cos (lslat) * georate * cos (lslong);
	z = 0;
	vec v(x,y,z);
	lsspeedbox->value(v.mag());
	vec att = r.unit();
	double t=0;
	double tstep = 0.001;
	
	//2. Load the parameters of the launch
	double oew = oewbox->value();
	double payload = payloadbox->value();
	double zfw = oew + payload;
	double fuel = fuelbox->value();
	double mass = zfw + fuel;
	double mdot = mdotbox->value();
	double vex = vexbox->value();
	double azimuth = azimuthbox->value() * M_PI / 180;
	double pitchtime = pitchtimebox->value() * 0.3048;
	double pitchrate = pitchratebox->value();
	double pitchterm = pitchtermbox->value();
	double pitch = 0;
	vec g = r*(-mu/pow(r.mag(),3));
	vec mg = g*mass;
	towbox->value(mg.mag());
	thrustbox->value(mdot*vex);
	
	//3. Liftoff from the pad
	do {
		//mdv = cdm + mgdt
		vec g = r*(-mu/pow(r.mag(),3));
		v = v + att*vex*mdot*tstep/mass + g*tstep;
		vec gforce = att*vex*mdot/mass;
		r = r + v*tstep;
		mass -= mdot*tstep;
		t += tstep;
		wtf << t << " " << r.mag() - georad << " " << v.mag() << " " << gforce.mag()/9.81 << endl;
	}
	while ( r.mag() - georad < pitchtime && mass > zfw);

The problem bit is that upon reaching the altitude of pitchtime, the rocket is then supposed to pitch over at an inputted rate toward a final pitch at an inputted azimuth. The describe the attitude of the rocket, I created the vec att, which is a unit vector pointing in the direction of the rocket (vec is a class defined to act like a 3D physics vector). The tricky part is that I'm not sure how to alter the value of the vec att through the pitch rate and the azimuth.

[Launch window]

sorry, I can't help

[Nicolas]

you have a timestep "tstep" and a pitch rate "prate". That gives an angle "deltapitch" = prate*tstep per time interval. You turn the vector att around it's local pitch axis with this angle. To find the components of the vector, keep track of the roll, pitch and yaw angles:

pitch = pitch + detapitch

Now apply Euler angles to get the parameters of vec:

Euler angles

You can calculate the matrices B,C and D. From that, you can calculate matrix A.

vec = A.vec gives the new vector. (vec = [x y z]' in my case)

Unless I'm mistaken here, I'm writing this on the fly .

I explained things for pitch; keep track of roll and yaw in the same way and take care of order of angles and matices.

[Glom]

Thanks, Nicolas. I'll have to read through that multiple times in order to grasp it.

For the moment though, I cheated a little and eliminated azimuth control. The rocket can only be launched easterly for the moment. That was a simplification that made the maths easier. What I did was I worked out a local North unit vector. The North vector would lie in the plane of r and K, but would be perpendicular to r. The solution was obviously cross products. N = r × ( K × r ). Then I apply my trusty unit vector operator to give it unit magnitude.

The next step was to get the delta att vector. The simplest way to find it would be to imagine that delta att was perpendicular to att. That way, the magnitude would be tan(delta pitch). If I limited the direction to easterly (perpendicular to N), then the direction would be given by N × att. Since those are both unit vectors, the magnitude is unity before scaling it with the tan(delta pitch). Then I add delta att to att to get the new att, which has the right direction, but unfortunately, not unit magnitude, so I simply apply the unit operator.

Easy!

Here's the code I have for the burn of the first stage (the rocket is now three stage).

do {
		//mdv = cdm + mgdt
		vec g = r*(-mu/pow(r.mag(),3));
		vec N = (r.cross( K.cross(r) )).unit();
		string phase = "Straight ascent";
		if ( (r.mag() - georad) > pitchtime && att.arg(r) < pitchterm ) {
			att = (att + N.cross(att)*tan(pitchrate*tstep)).unit();
			phase = "Pitching";
		}
		v = v + att*vex1*mdot1*tstep/mass1 + g*tstep;
		vec gforce = att*vex1*mdot1/mass1;
		r = r + v*tstep;
		mass1 -= mdot1*tstep;
		t += tstep;
		telapsebox->value(t);
		heightbox->value(r.mag()-georad);
		speedbox->value(v.mag());
		pathbox->value(att.arg(r));
		gforcebox->value(gforce.mag());
		massbox->value(mass1);
		phasebox->value("1st stage");
		wtf << t << " " << r.mag() - georad << " " << v.mag() << " " << gforce.mag()/9.81 << " " << att.arg(r) << " " << phase << endl;
	}
	while ( mass1 > zfw1 );

As you can see, the pitching happens only when the height of the rocket is sufficient, and before the pitch reaches the terminal pitch attitude.

Other parts of difficulty, is getting the widgets to update as the program does the calculations. Also, the compiler didn't seem to like it when I tried to output the string phase in the phasebox.