Calculating air resistance

[icarusuk 0]

Ok Ive run into some problems with my arty script.  I have it all mathmatically setup for perfect conditions, however when I give a projectile and ideal velocity, it falls short, because of drag!

How can I find out the effects of drag, so i know what velocity I really need to use to get a projectile to land exactly on a specified target?

[Guest BratZ]

If I knew someone with lets say a piece of artillery.

Then after multiple tests,we could have these figures.

I don't recall any lists offhand with bulletdrop (maybe)

but artillery is fairly new

[Guest]

Two options:

1) Write a scirpt that solves the differential equation

</span><table border="0" align="center" width="95%" cellpadding="3" cellspacing="1"><tr><td>Code Sample </td></tr><tr><td id="CODE">

ds^2/dt^2 = -g + k(ds/dt)^2/m

<span id='postcolor'>

Get the coefficient of drag by empirical testing.

2) Use a neural network. I'm going to release the CoC Neural Network pack soon, but I need more time to write some examples and documentation.

[icarusuk 0]

does ofp scripting do differentials??

[Guest]

No, you have to implement a solver. An implementation of for instance Runge-Kutta's method should not be too difficult.

Edit: Forgot to mention: You must then transform the second order diffeq to a system of first order of diffeqs.

s'' = k*(s')^2 - g can be transformed to:

x = y'

x' = k*(x)^2-g

Solve for x with the method of your choice.

[Dschulle 0]

@Runge-Kutta look here:

http://www.ofpec.com/editors/funcref.php?filter_func=25

[icarusuk 0]

Im assuming K is a constant? What might it be?

[Guest]

That's the drag constant. You'll have to determine that one by testing. Start with k = 0 and then narrow down to the value that corresponds best with OFP's.

[icarusuk 0]

Sorry to be a pain, but I need a little more help on this.

I understand I have to use the Kutta method to differentiate, but Im not really sure what I have to differentiate and what Im trying to find as a result.

ds^2/dt^2 = -g + k(ds/dt)^2/m

-g is the acceleration due to gravity yes?

K is the drag coefficient

ds/dt is the change in displacement over time (which is squared)

m is the mass ?

So whats my end result wanting to be? I assume I want some kind of function out of the other end but Im unsure as to what exactly Im trying to get, and how to get it.

[Guest]

The Runge-Kutta method solves numerically a differential equation.

The answer you will get is s(t) which is displacement over time.

ds^2/dt^2 = second derivative of the displacement = acceleration

-g = gravitational acceleration

k = drag coefficient

(ds/dt)^2 = velocity squared

Solving the differential equation you will get the height of your projectile at a given time.

[icarusuk 0]

And why do I care about the height of my projectile at any given time?

[Guest]

Because when it becomes zero it is at the impact point. If you know when then you know how fast, so to say. And if you know how fast then you know where

[bn880 5]

there's a riddle for ya

[icarusuk 0]

Hmmm, Im all confused as to why I need to know the final velocity and position the shell lands.

[Guest]

I thought that was what you wanted to get What is it you want to get?

I thought you were making an artillery script and wanted to predict where the shells would land.

Nevertheless the differential equation gives you the complete trajectory which gives all the information you could possibly want.

[icarusuk 0]

Its possiby what I want, I just dont know why I want it!  Though I dont want to predict the landing position, tell it where to land and thus calculate the velocity.

Anyway.  1st year aerodynamics tells me that;

DragForce = 0.5 * p(rho) * V^2 * S * Cd

(I knew there was a point in going to university)

With some expoerimentation I can get the Coeffecient of drag for a shell (Cd) and the surface area (S).  But I think what I need to do is make some simultaneous equations that will work out it all out.

For example my current script looks like this;

---------------------------------

;Calculate the required velocities using the dynamics equation : u = (s - 0.5 * a * t^2) / t

;This gives a value of V that is perfect for an enviroment with no drag, which is obviously not the enviroment OFP models.

_XvelocityPerfect = _Xdistance/_FlightTime

_ZVelocityPerfect = _Zdistance/_FlightTime

_YvelocityPerfect = (_Ydistance - (0.5*(-9.81)*(_FlightTime^2))) / _FlightTime

;Flashpoint models drag.  Unfortunately dynamically calculating drag is a right pain, so I use a half arsed model and hope

;for the best.

;D= 0.5 * p * V^2 * S * Cd

_SurfaceArea = 0.01

_DragCoeff = 0.4

_density = 1.225

_XDragForce = 0.5 * _density * (_XvelocityPerfect^2) * _SurfaceArea * _DragCoeff

_ZDragForce = 0.5 * _density * (_ZvelocityPerfect^2) * _SurfaceArea * _DragCoeff

_YDragForce = 0.5 * _density * (_YvelocityPerfect^2) * _SurfaceArea * _DragCoeff

;Its not the best way of doing it, but now the velocity is adjusted for drag

_Xvelocity = (_Xdistance - (0.5*(+_Xdragforce)*(_FlightTime^2))) / _FlightTime

_Zvelocity = (_Zdistance - (0.5*(+_Zdragforce)*(_FlightTime^2))) / _FlightTime

_Yvelocity = (_Ydistance - (0.5*(-9.81 + _YDragForce)*(_FlightTime^2))) / _FlightTime

;Create the first mortar round slightly above the launcher.

_shell1 = _shelltype Camcreate [Getpos _Launcher Select 0, Getpos _Launcher Select 1, (Getpos _Launcher Select 2)+1]

;Give the Shell the correct velocity

_Shell1 SetVelocity [ _Xvelocity , _ZVelocity , _YVelocity ];

------------------------------------------

Now hopefully you can see what Ive done there.  I work out what I call a perfect velocity (no drag), then using that drag force I work out the velocities again with the same equation and get a value of V that is adjusted for drag.

I know its a very shoddy way around it, it doesnt adjust the drag force for the variable velocity through the arc, but its the start of somthing.

This way round however, at "runway length" the shells are near perfect, however the shorter the range between target and launcher the more the shell drops short.

Any ideas without huge differential equations I dont understand?

[Guest]

Perhaps if you used a more complex polynomial model instead of a square one. So instead of

DragForce = k*v^2 you can say

DragForce = a*v^2 + b*v + c and decide a,b,c by measurement.

The other really simple alternative is as I suggested earlier using a neural network. As inputs you put the desired impact point and if you want a high or low orbit. As output you set the desired launch velocities. Shoot of a couple of shells at random and register where they land. Shove it all into the neural network and train it. It is a fairly simple problem so it should converge to a solution pretty fast. After it is trained just use it by saying where you want to land and the net will say what velocities you have to fire it at.

[icarusuk 0]

Ok, so where can I get my hands on one of these nural networks then?

You mentioned a CoC one, but said you needed to write more documentation for it, I took that to mean its not released yet.

[Spitfire 0]

*waits for Suma to step in at any moment to reveal the drag coefficient and the correct formula used in the game*

[icarusuk 0]

I was hoping someone would

[Spitfire 0]

It can be a real pain finding out the formula for the force of air resistance (it's not necessarily that real life physics formula you suggested). Even if it was, I remember it being a b*tch to find out the drag coefficients of various objects in physics lab a couple of years ago. Therefore using Denoir's fuzzy logic neural network solution could give better results.

In addition, I would rather put a gun in my mouth than design a solver from the scratch with OFP scripting commands.

While I haven't looked through all artillery/projectile motion scripts I've seen for OFP, I got the impression that no-one has actually made a working drag formula based on real life physics, but used workarounds producing rough approximations instead.

All I can say, good luck with your attempt

[Prospero 1]

I s'pose this may be helpful. Then again, it may not be. I wrote it as a framework for this sort of stuff a while back. Knock yourselves out.

----------

; Assume that we've already created a reference object called "referenceobject" on the map in the Mission Editor to define the start position for the simulation. The following line of code will place this reference object 3 metres off the ground - for purely arbitrary reasons. Also assume that we've already created a trigger called "base" to aid in the determination of absolute height.

referenceobject setpos [(getpos referenceobject select 0), (getpos referenceobject select 1), (getpos referenceobject select 2) + 3]

; This is the global array we'll use to output the required simulation variables.

Output = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

; Physical constants. At least, we'll treat them as constants. Note that rho (it's a letter in the greek alphabet) is the air density. Gravity can be "switched off" by making g equal to zero (very useful for test purposes).

;_g = 9.806650

_g = 0

_rho = 1.225

; We convert initial x, y and z position states from ENU to NED since this flight-model adheres to normal aeronautical convention for body-axes and (flat Earth) Earth-axes coordinate systems, namely NED. To transform ENU coordinates (the OFP world coordinate system is ENU) to NED coordinates (or the reverse process), simply swap over x and y, and change the sign of z. The following four lines read, "Initial x coordinate in NED equals initial y coordinate in ENU... etc." Also note that the NED z value we'll be feeding into the sim is absolute height, rather than height AGL.

_x = getpos referenceobject select 1

_y = getpos referenceobject select 0

base setpos [_y, _x, 0]

_z = -((getpos referenceobject select 2) - (getpos base select 2))

; Initialize various other aircraft states. For convenience, I've used OFP's rad function here so that one can plug in values in degrees per second for the body-axes angular rates and degrees for the Euler angles. The rad function simply converts the values to radians for use in the flight-model.

_u = 0

_v = 0

_w = 0

_p = rad(0)

_q = rad(0)

_r = rad(0)

_phi = rad(0)

_theta = rad(0)

_psi = rad(0)

_oudot = 0

_ovdot = 0

_owdot = 0

_opdot = 0

_oqdot = 0

_ordot = 0

_oxdot = 0

_oydot = 0

_ozdot = 0

_Vt = 0

_alpha = 0

_beta = 0

_lVt = 0

_lalpha = 0

_lbeta = 0

; Build the initial attitude quaternion from the initial Euler angles. The Euler angles are phi (roll), theta (pitch) and psi (yaw). We will also initialize the quaternion parameters' rates. Incidentally, we often talk about and list the Euler angles in this order, but when it comes to applying rotations, we will always perfom yaw first (around the z-axis), pitch second (around the y-axis), and roll third (around the x-axis). This is the so-called 3-2-1 order prescribed by normal aeronautical convention. I should also mention that I will never refer to a quaternion's four parameters by numbered indices, like 0, 1, 2 and 3. Of course, this would be the natural thing to do if the quaternion was stored in an array, but it can quickly lead to confusion because some people put the "angle" parameter first, whilst other people put it last. Therefore, for the sake of clarity, I will ensure that the variable for the "angle" parameter always contains a letter "w", and the variables for the three parameters that make up the Euler "axis vector" always contain either a letter "x", "y" or "z" as appropriate. This way, there should be no confusion.

_q1w = cos(deg(_phi / 2)) * cos(deg(_theta / 2)) * cos(deg(_psi / 2)) + sin(deg(_phi / 2)) * sin(deg(_theta / 2)) * sin(deg(_psi / 2))

_q1x = sin(deg(_phi / 2)) * cos(deg(_theta / 2)) * cos(deg(_psi / 2)) - cos(deg(_phi / 2)) * sin(deg(_theta / 2)) * sin(deg(_psi / 2))

_q1y = sin(deg(_phi / 2)) * cos(deg(_theta / 2)) * sin(deg(_psi / 2)) + cos(deg(_phi / 2)) * sin(deg(_theta / 2)) * cos(deg(_psi / 2))

_q1z = cos(deg(_phi / 2)) * cos(deg(_theta / 2)) * sin(deg(_psi / 2)) - sin(deg(_phi / 2)) * sin(deg(_theta / 2)) * cos(deg(_psi / 2))

_odq1w = 0

_odq1x = 0

_odq1y = 0

_odq1z = 0

; Normalize the initial attitude quaternion. Note that we will always normalize a quaternion before building a rotation matrix from its parameters.

_e = sqrt(_q1w^2 + _q1x^2 + _q1y^2 + _q1z^2)

_q1w = _q1w / _e

_q1x = _q1x / _e

_q1y = _q1y / _e

_q1z = _q1z / _e

; Build the initial Direction Cosine Matrix (also known as the DCM) from the initial attitude quaternion.

_m11 = _q1w^2 + _q1x^2 - _q1y^2 - _q1z^2

_m12 = 2 * (_q1x * _q1y + _q1w * _q1z)

_m13 = 2 * (_q1x * _q1z - _q1w * _q1y)

_m21 = 2 * (_q1x * _q1y - _q1w * _q1z)

_m22 = _q1w^2 - _q1x^2 + _q1y^2 - _q1z^2

_m23 = 2 * (_q1y * _q1z + _q1w * _q1x)

_m31 = 2 * (_q1x * _q1z + _q1w * _q1y)

_m32 = 2 * (_q1y * _q1z - _q1w * _q1x)

_m33 = _q1w^2 - _q1x^2 - _q1y^2 + _q1z^2

; Initialize the control inputs. Conveniently, up to six axes can be handled by my ProSphere Controller software - one for each degree of freedom. ProSphere Controller consists of a set of OFP scripts I wrote which a) read the "main" game movement keys, b) process the keystrokes to provide six simultaneous pseudo-analog axes of real-time control input, and c) load the axes data into a globally accessible array called PS.

_ps = PS

; Initialize variables related to what I always refer to as the "weather-wind". Bear in mind that its xwnd, ywnd and zwnd vector components are specified in NED Earth-axes.

_xwnd = 0

_ywnd = 0

_zwnd = 0

; Overall aircraft mass, mass moments of inertia and ixz product of inertia. We'll treat them all as constants. Here, we're using typical values for a small helicopter UAV such as the Yamaha R50. The R50 weighs about 45 kg and has a rotor diameter of approximately 3 metres.

_m = 45.0

_ix = 2.0

_iy = 6.2

_iz = 6.0

_ixz = 0.0

_id = _ix * _iz - _ixz^2

; Aircraft component data.

; ...blah blah blah...

; Initialize variables related to the timing scheme.

_odt = 0

_ot = Time

; Main loop.

#lp

~0.0001

_nt = Time

;_dt = _nt - _ot

_dt = 0.01

; Get the control inputs from ProSphere Controller.

_ps = PS

; The weather-wind velocity vector's components in NED Earth-axes (in metres per second). For now, we'll specify constant zero wind.

_xwnd = 0

_ywnd = 0

_zwnd = 0

; Transform from Earth-axes to body-axes using the Direction Cosine Matrix.

_uwnd = _xwnd * _m11 + _ywnd * _m12 + _zwnd * _m13

_vwnd = _xwnd * _m21 + _ywnd * _m22 + _zwnd * _m23

_wwnd = _xwnd * _m31 + _ywnd * _m32 + _zwnd * _m33

; Subtract these values from u, v, and w. We'll then be using the resulting relative values (ur, vr and wr) to calculate Vt, alpha and beta, and hence the lift, drag and sideforce forces in wind-axes for each of the aircraft's components.

_ur = _u - _uwnd

_vr = _v - _vwnd

_wr = _w - _wwnd

; Determine aircraft CG's velocity, alpha and beta. I've put in some logic to truncate the value returned by the atan2 function to a range of +/- pi/2 radians for both alpha and beta values. The logic should be removed if you want to use a range of +/- pi radians instead. I've also put in a trap for when the atan2 function is passed two zero-value arguments.

_Vt = sqrt(_ur^2 + _vr^2 + _wr^2)

?(_wr == 0) && (_ur == 0): goto "sp1"

_alpha = rad(_wr atan2 _ur)

?(_alpha > pi / 2) && (_alpha <= pi): _alpha = pi - _alpha

?(_alpha > -pi) && (_alpha <= -pi / 2): _alpha = abs(_alpha) - pi

#sp1

?(_vr == 0) && (_ur == 0): goto "sp2"

_beta = rad(_vr atan2 _ur)

?(_beta > pi / 2) && (_beta <= pi): _beta = pi - _beta

?(_beta > -pi) && (_beta <= -pi / 2): _beta = abs(_beta) - pi

#sp2

; /----------------------------------------------/

; /**********************************************/

; / The aircraft component build-up starts here. /

; /**********************************************/

; /----------------------------------------------/

; Now we determine the local velocity, local alpha and local beta at the position of the structural component's aerodynamic centre (often refered to as its ac). We will call the local velocity components lur, lvr and lwr. The cross product of the aircraft's angular rate vector and the structural component's position vector is used to calculate them.

_lur = _ur + (_q * _htdz -_r * _htdy)

_lvr = _vr + (_r * _htdx - _p * _htdz)

_lwr = _wr + (_p * _htdy -_q * _htdx)

; Local velocity, local alpha and local beta.

_lVt = sqrt(_lur^2 + _lvr^2 + _lwr^2)

?(_lwr == 0) && (_lur == 0): goto "sp3"

_lalpha = rad(_lwr atan2 _lur)

?(_lalpha > pi / 2) && (_lalpha <= pi): _lalpha = pi - _lalpha

?(_lalpha > -pi) && (_lalpha <= -pi / 2): _lalpha = abs(_lalpha) - pi

#sp3

?(_lvr == 0) && (_lur == 0): goto "sp4"

_lbeta = rad(_lvr atan2 _lur)

?(_lbeta > pi / 2) && (_lbeta <= pi): _lbeta = pi - _lbeta

?(_lbeta > -pi) && (_lbeta <= -pi / 2): _lbeta = abs(_lbeta) - pi

#sp4

; ...blah blah blah...

; Relate the signs of these forces to wind-axes convention.

_x1 = -_drag

_y1 = _sideforce

_z1 = -_lift

; We now need to express the wind-axes force vector in aircraft body-axes. I find it helpful to think of the process as follows. Rather than rotate the wind-axes vector, we'll leave it where it is, and instead imagine that we're rotating the wind-axes coordinate frame so that it is aligned with the body-axes frame. Then we'll express the original vector in this new coordinate frame. The wind-axes coordinate system can be aligned with the body-axes coordinate system by two consecutive rotations. First, a rotation of -lbeta (minus lbeta) around the z-axis, followed by a rotation of lalpha around the updated y-axis. The idea then, is to build an "overall" rotation matrix that will accomplish this. In fact, the process turns out to be identical to that which we'd use for rotating a vector - and indeed, it appears that this is precisely what we're doing - but in the final stage of the process we don't apply the rotation matrix to the coordinate frame, but rather we apply its inverse (or transpose) to the vector. What we end up with is an expression for the original vector in the new coordinate system, rather than an expression for the rotated vector in the original coordinate system. For clarity, I have reproduced the process in full here, rather than take any shortcuts.

; Build the primitive quaternion representing a rotation by -lbeta around the z-axis.

_q2wa = cos(deg(-_lbeta/2))

_q2xa = 0

_q2ya = 0

_q2za = sin(deg(-_lbeta/2))

; Build the primitive quaternion representing a rotation by lalpha around the y-axis.

_q2wb = cos(deg(_lalpha/2))

_q2xb = 0

_q2yb = sin(deg(_lalpha/2))

_q2zb = 0

; Combine the two quarternions by quaternion multiplication in the order of z-axis rotation first, y-axis rotation second. The result is a single quaternion describing both rotations in one shot - as if they were "consecutive".

_q2w = -_q2xa * _q2xb - _q2ya * _q2yb - _q2za * _q2zb + _q2wa * _q2wb

_q2x = _q2wa * _q2xb - _q2za * _q2yb + _q2ya * _q2zb + _q2xa * _q2wb

_q2y = _q2za * _q2xb + _q2wa * _q2yb - _q2xa * _q2zb + _q2ya * _q2wb

_q2z = -_q2ya * _q2xb + _q2xa * _q2yb + _q2wa * _q2zb + _q2za * _q2wb

; Normalize the resulting quaternion.

_e = sqrt(_q2w^2 + _q2x^2 + _q2y^2 + _q2z^2)

_q2w = _q2w / _e

_q2x = _q2x / _e

_q2y = _q2y / _e

_q2z = _q2z / _e

; Build the rotation matrix from the quaternion parameters.

_r11 = _q2w^2 + _q2x^2 - _q2y^2 - _q2z^2

_r12 = 2 * (_q2x * _q2y - _q2z * _q2w)

_r13 = 2 * (_q2x * _q2z + _q2y * _q2w)

_r21 = 2 * (_q2x * _q2y + _q2z * _q2w)

_r22 = _q2w^2 - _q2x^2 + _q2y^2 - _q2z^2

_r23 = 2 * (_q2y * _q2z - _q2x * _q2w)

_r31 = 2 * (_q2x * _q2z - _q2y * _q2w)

_r32 = 2 * (_q2y * _q2z + _q2x * _q2w)

_r33 = _q2w^2 - _q2x^2 - _q2y^2 + _q2z^2

; Now, as stated previously, the inverse (or transpose) of this rotation matrix can be used to express any wind-axes vector in body-axes. By the way, the inverse and transpose are exactly the same thing since we are dealing with an orthogonal rotation matrix. Think transpose - it's a lot easier, not to mention faster to calculate - we simply swap the rows and columns of the matrix.

;_x2 = _x1 * _r11 + _y1 * _r21 + _z1 * _r31

;_y2 = _x1 * _r12 + _y1 * _r22 + _z1 * _r32

;_z2 = _x1 * _r13 + _y1 * _r23 + _z1 * _r33

; Alternatively, one can work out aerodynamic forces in body-axes directly. This is a useful shortcut, especially if we just want to approximate aerodynamic forces acting on simple shapes.

_x2 = -1.0 * abs(_lur) * _lur * _rho / 2

_y2 = -1.0 * abs(_lvr) * _lvr * _rho / 2

_z2 = -1.0 * abs(_lwr) * _lwr * _rho / 2

; The forces on, and moments around the aircraft CG position arising from the forces just calculated for the structural component in question. We can use the cross product of the structure's position vector and the structure's force vector to work out the moments around the CG.

_htx = _x2

_hty = _y2

_htz = _z2

_htl = _htdy * _z2 - _htdz * _y2

_htm = _htdz * _x2 - _htdx * _z2

_htn = _htdx * _y2 - _htdy * _x2

; /--------------------------------------------/

; /********************************************/

; / The aircraft component build-up ends here. /

; /********************************************/

; /--------------------------------------------/

; Determine total forces and moments acting on the aircraft. Note that I am also feeding in some control input directly - for test purposes.

_Xtot = (_ps select 4) * 10000 + _htx

_Ytot = (_ps select 3) * 10000 + _hty

_Ztot = -(_ps select 2) * 10000 + _htz

_Ltot = (_ps select 5) * 500 + _htl

_Mtot = -(_ps select 1) * 500 + _htm

_Ntot = (_ps select 0) * 500 + _htn

; Now we'll deal with the rigid body Equations Of Motion (the EOM).

; First, determine the required sines and cosines of the Euler angles calculated in the previous cycle, bearing in mind that the Euler angles, body angular rates and body angular accelerations are expressed in radians within this maths flight-model. Operation Flashpoint trig functions, however, expect to be passed angles expressed in degrees, so we'll take care of this by converting them appropriately as we go (using the deg function). Note, by the way, that C/C++, FORTRAN, Visual BASIC and MATLAB differ in that their trig functions deal in radians, so no conversion would be necessary if one was to implement this code using these languages.

_sph = sin(deg(_phi))

_cph = cos(deg(_phi))

_sth = sin(deg(_theta))

_cth = cos(deg(_theta))

; Determine the body-axes acceleration components, taking gravity into account.

_udot = (_r * _v) - (_q * _w) + (_Xtot / _m) - (_g * _sth)

_vdot = (_p * _w) - (_r * _u) + (_Ytot / _m) + (_g * _sph * _cth)

_wdot = (_q * _u) - (_p * _v) + (_Ztot / _m) + (_g * _cph * _cth)

_pdot = (_iz / _id) * (_Ltot + _ixz * _p * _q - (_iz - _iy) * _q * _r) + (_ixz / _id) * (_Ntot - _ixz * _q * _r - (_iy - _ix) * _p * _q)

_qdot = (_Mtot - (_ix - _iz) * _p * _r - _ixz * (_p^2 - _r^2)) / _iy

_rdot = (_ixz / _id) * (_Ltot + _ixz * _p * _q - (_iz - _iy) * _q * _r) + (_ix / _id) * (_Ntot - _ixz * _q * _r - (_iy - _ix) * _p * _q)

; Integrate trapezoidally to determine body-axes velocity components.

_u = _u + 0.5 * (_udot * _dt + _oudot * _odt)

_v = _v + 0.5 * (_vdot * _dt + _ovdot * _odt)

_w = _w + 0.5 * (_wdot * _dt + _owdot * _odt)

_p = _p + 0.5 * (_pdot * _dt + _opdot * _odt)

_q = _q + 0.5 * (_qdot * _dt + _oqdot * _odt)

_r = _r + 0.5 * (_rdot * _dt + _ordot * _odt)

; Update the attitude quarternion. First, we calculate the rates of change of its parameters based on p, q and r (the body-axes angular velocity components).

_dq1w = -0.5 * (_q1x * _p + _q1y * _q + _q1z * _r)

_dq1x = 0.5 * (_q1w * _p - _q1z * _q + _q1y * _r)

_dq1y = 0.5 * (_q1z * _p + _q1w * _q - _q1x * _r)

_dq1z = 0.5 * (-_q1y * _p + _q1x * _q + _q1w * _r)

; Integrate trapezoidally.

_q1w = _q1w + 0.5 * (_dq1w * _dt + _odq1w * _odt)

_q1x = _q1x + 0.5 * (_dq1x * _dt + _odq1x * _odt)

_q1y = _q1y + 0.5 * (_dq1y * _dt + _odq1y * _odt)

_q1z = _q1z + 0.5 * (_dq1z * _dt + _odq1z * _odt)

; Normalize the updated quaternion. I have not put in a catch for the possibility of a zero divisor since a quaternion with all four parameters equal to zero would not represent a valid rotation in the first place.

_e = sqrt(_q1w^2 + _q1x^2 + _q1y^2 + _q1z^2)

_q1w = _q1w / _e

_q1x = _q1x / _e

_q1y = _q1y / _e

_q1z = _q1z / _e

; We now build the Direction Cosine Matrix from the updated quaternion parameters. Converting the attitude to this form is useful for a number of reasons. First, it will allow us to extract the Euler angles. Second, it will prove useful for expressing the weather-wind vector, originally described in Earth-axes, in body-axes. Third, by using its transpose, we can express the linear velocity vector, originally calculated in body-axes, in Earth-axes. Incidentally, the flight-model wants to deal with Euler angles as input and output, but we perform the attitude calculation with quaternion maths internally to alleviate problems such as numerical drift and gimbal lock.

_m11 = _q1w^2 + _q1x^2 - _q1y^2 - _q1z^2

_m12 = 2 * (_q1x * _q1y + _q1w * _q1z)

_m13 = 2 * (_q1x * _q1z - _q1w * _q1y)

_m21 = 2 * (_q1x * _q1y - _q1w * _q1z)

_m22 = _q1w^2 - _q1x^2 + _q1y^2 - _q1z^2

_m23 = 2 * (_q1y * _q1z + _q1w * _q1x)

_m31 = 2 * (_q1x * _q1z + _q1w * _q1y)

_m32 = 2 * (_q1y * _q1z - _q1w * _q1x)

_m33 = _q1w^2 - _q1x^2 - _q1y^2 + _q1z^2

; Now we extract the Euler angles from the Direction Cosine Matrix, which, according to aeronautical convention, are applied in the order yaw first (around the z-axis), pitch second (around the y-axis), and roll third (around the x-axis). Remember, Operation Flashpoint's arc trig functions return degrees rather than radians, so we convert to radians as we go. I've put in a catch for when the atan2 function is passed two zero-value arguments (which occurs when theta is +90 or -90 degrees) - we just skip the calculations for phi and psi entirely - and have also ensured that the asin function isn't fed with a value of magnitude greater than unity due to rounding error. I realise that the atan2 catch needs to be improved. We should be using Jack D. Crenshaw's way of determining phi and psi. His method supposedly obviates the singularity that normally prevents calculating them when theta is +/- 90 degrees.

?(abs(_m13) > 1): goto "sb1"

_theta = rad(asin(-_m13))

#rn1

?(_m23 == 0) && (_m33 == 0): goto "sp5"

_phi = rad(_m23 atan2 _m33)

#sp5

?(_m12 == 0) && (_m11 == 0): goto "sp6"

_psi = rad(_m12 atan2 _m11)

#sp6

; Transform body-axes linear velocity components to determine Earth-axes linear velocity components in NED. We can use the inverse (or transpose) of the Direction Cosine Matrix to do this.

_xdot = _u * _m11 + _v * _m21 + _w * _m31

_ydot = _u * _m12 + _v * _m22 + _w * _m32

_zdot = _u * _m13 + _v * _m23 + _w * _m33

; Integrate trapezoidally to determine Earth-axes position in NED.

_x = _x + 0.5 * (_xdot * _dt + _oxdot * _odt)

_y = _y + 0.5 * (_ydot * _dt + _oydot * _odt)

_z = _z + 0.5 * (_zdot * _dt + _ozdot * _odt)

; Print useful data such as the state variables to the screen.

;hint format ["X  %1\nY  %2\nZ  %3\nL  %4\nM  %5\nN  %6\n\nudot  %7\nvdot  %8\nwdot  %9\npdot  %10\nqdot  %11\nrdot  %12\n\nu  %13\nv  %14\nw  %15\np  %16\nq  %17\nr  %18\n\nVt  %19\nalpha  %20\nbeta  %21\n\nxdot  %22\nydot  %23\nzdot  %24\n\nx  %25\ny  %26\nz  %27\nphi  %28\ntheta  %29\npsi  %30\n\ndt  %31\n\nlVt  %32\n", _Xtot, _Ytot, _Ztot, _Ltot, _Mtot, _Ntot, _udot, _vdot, _wdot, _pdot, _qdot, _rdot, _u, _v, _w, _p, _q, _r, _Vt, deg(_alpha), deg(_beta), _xdot, _ydot, _zdot, _x, _y, _z, deg(_phi), deg(_theta), deg(_psi), _dt, _lVt]

; Store values which will be used in the trapezoidal integration scheme next cycle.

_oudot = _udot

_ovdot = _vdot

_owdot = _wdot

_opdot = _pdot

_oqdot = _qdot

_ordot = _rdot

_oxdot = _xdot

_oydot = _ydot

_ozdot = _zdot

_odq1w = _dq1w

_odq1x = _dq1x

_odq1y = _dq1y

_odq1z = _dq1z

_odt = _dt

; Bear in mind that the elements of the output array are, variously, in units of metres per second, metres and radians. Also note that they are all in NED. If, for example, one has an external display script that works in ENU (the OFP world coordinate system), one would need to transform the elements from NED Earth-axes to ENU Earth-axes - typically by swapping over the x and y terms, and changing the sign of the z term when applying setPos or setVelocity commands. The Direction Cosine Matrix (again in NED) is also included as part of the output.

Output = [_xdot, _ydot, _zdot, _x, _y, _z, _phi, _theta, _psi, _m11, _m12, _m13, _m21, _m22, _m23, _m31, _m32, _m33]

_ot = _nt

goto "lp"

#sb1

?(_m13 > 0): _m13 = 1; goto "rn1"

_m13 = -1

goto "rn1"

----------

Prospero

[Spitfire 0]

Geez, parts of this would've been helpful to me a couple of months ago when I had to dig out my old math books to use Euler's theorems to rotate vectors. Well, that code is finished now, so thanks for nothing, Prospero

[icarusuk 0]

Ok, so when my brain dribbled out of my ear, that was a bad thing right?

Errr,sodding heck. You know how there are loads of arty scripts around? How do they combat drag?

[Doolittle 0]

They don't.  All of the original artillery scripts created shells above where they were to land, then let gravity do the rest.  It could be dead on.

When 1.85 came out, we could then add velocity to objects, thus you have the "true" artillery where you try to hit something by firing a shell from some other location, giving it initial velocity and then letting the game world handle the rest.

I wrote a script that did that but I couldn't figure out drag.  And after that you still have to figure out different Y heights of target and cannon.

I know gravity affects the shell in the game, but does drag too, as far as up/down motion is concerned?

I've tried simple solutions like make the distance actually some percentage further than it actually is but that doesn't work.

What I think you need to do is after you have initial velocity you find out how long it will be in the air, then you say that every second the forward speed drops some percentage and so you add onto the initial what that lost speed would be.

I tried firing a round straight up into the air and then marking how long it was in the air to get the exact gravity constant.  Also the mass of the object seems to change.  A table or phone flies farther than a shell.

Check out http://www.imgstudio.com/photo-s....ittle01, where I have pics of testing.

Doolittle