Trebuchet simulation&some newbe question

The July issue of "Scientific American" describes the Trebuchet,

an ancient (500 AD) war machine capable of hurling massive

projectiles an astonishing distance.

The magazine cover shows a car flying through the air after being

launched by a modern reconstruction of the machine. The article has

pictures of an upright piano subjected to a similar indignity.

The authors mention that a computer program was used to explore

performance parameters of the machine so I resolved to see what I

could do with the simulation problem.

As might be expected for such an ancient war engine, the basic

idea is quite simple, but I found that there are some subtle

nuances in the design when at last I got my program to run and

produce what I think are reasonable results.

The Trebuchet is an unequal arm lever pivoted on an elevated

bearing so the arm is free to swing in a vertical plane. A heavy

counterweight is fixed to the short end of the lever and a sling

(the David & Goliath type) is attached to the long end. The

mechanism is prepared to launch a projectile by pulling down the

long end of the lever, placing a projectile in the sling pocket and

hooking one of the sling ropes to a hook that is designed to

release when the arm reaches the proper angle. The other sling rope

remains attached to the arm.

My simulation revealed that designing this hook might be a mite

tricky along with fixing the relation between lengths of the sling

and the long arm of the lever.

After pursuing many dead end paths, I realized that the machine is

a compound pendulum where the overall system (arm, counterweight,

and projectile) rotates about the pivot with an angular

acceleration equal to the unbalanced torque divided by the moment

of inertia of the system.

The sling in turn rotates about the long end of the lever arm with

an angular acceleration given by a similar equation but in this

case the moment of inertia is constant; the mass (in slugs) of the

projectile times the square of the length of the sling rope.

To solve the problem first and second degree differential

equations relative to both the arm and sling must be integrated to

obtain the angle of the lever from its acceleration and the angle

of the sling from a force balance on the projectile. I maintained

that simple Euler integration should be adequate because, given

small time increments in difference equations, all the derivatives

would change smoothly without inconvenient discontinuities. A

friend suggested that I use Runge-Kutta integration because the

double integration might introduce excessive error propagation.

My latest version of the simulation makes use of Adams predictor-

corrector integration which, I believe, produces a better result

but at the expense of complexity that I prefer to avoid here.

This brings me to my 'newbe' question about integration using J.

In fact I would appreciate it if someone could show me how to solve

the simulation in J. I have spent a good many hours trying the

iteration function f^:(n) g for integration subroutines, a loop

with label_loop. at the head of a list of J statements and

if. (200 >.N=.N+1) goto_loop.

end.

)

I'm missing something. It seems to me this should work but it

doesn't. I keep getting errors even when I simplify to the point of

just incrementing the counter, N. I'm using J for windows V2.05.

Here is my simulation program:

Please forgive me I'm using 'pigeon J' in order to transmit this

over the net so the APL program won't run until J symbols are

translated to corresponding APL characters. I have provided a

translation table at the end of the nomenclature list for those who

remain confirmed APL users. I have added one symbol :: to represent

the diamond APL symbol used to separate independent statements on

the same line.

Nomenclature list:

W1 = weight of counterweight (lbs)

m1 = mass of counterweight (slugs = W1% 32.2)

W2 = weight of projectile; m2 = W2%32.2

L1 = length short end of lever arm (ft); L2 = length of long end

SL = length of sling (ft); P = X-Y coordinates long end of lever

T = coordinates of projectile relative to P

S = coordinates of projectile (S=P+T )

dT = integration time step

a = (dV/dT); V = (dS/dT) projectile acceleration & velocity

wp (for w-prime) the angular acceleration of the lever arm

w = angular velocity of the lever arm

wp2 = angular acceleration of sling and projectile

w2 = angular velocity of sling and projectile

Translation table:

% = divide * = multiply o. = o (APL circle function)

=. is APL assignment arrow _ is APL negation function

^ = * APL power function $ is rho, the APL shape fuunction

NB. = APL lamp symbol (comment)

Simulation program with typical parameter values:

[0] {del.} ZOT (credits to the anteater in the comic strip BC)

[1] W1 =. 2000 :: W2 =. 20 :: L1=. 6 :: L2=. 30

[2] m1=.W1%32.2 :: m2=.W2%32.2 :: N=.1 :: dT=.0.01 :: w=.w2=.0

[3] A=. o.%6:: P=. L2*(1 _1)*1 2o.A:: T=. _8 _6 :: SL=.(+/T*T)^0.5

[4] A2 =. _1oT[1]%SL :: S=. P+T:: a=.V=.0 0::

[5] C=. m1*L1*L1 :: I2=. m2*SL*SL

[6] TB=.1 17$ 0,P,S,A,w,wp,A2,w2,wp2,V,a,T

[7]

[8] Ln0: I=.(C +m2*(s=.+/S*S)) :: s=. s^0.5 NB. Moment of inertia

[9] wp=. ((W1*L1*1oA) -W2*S[1]) % I NB.System angular acceleration

[10] w =. w +dT*wp NB. angular velocity

[11] A =. A +dT*w NB. new lever arm angle

[12] P =. L2*(1 _1)*1 2 o A NB. end of lever arm

[13] Z =. (w,wp) o.* s * 2 1 o A NB. angular velocity, acceleration

times vector s from origin to S

[14] V =.Z[1;] :: a =. Z[2;] NB.projectile velocity & acceleration

[15] f =. m2*a :: f[2]=. (f[2]>W2)*f[2]-W2 NB.subtract gravity

force after lift-off

[16] wp2=.(-/f*T[2 1]) % I2 NB. vector cross product of force and

normal to vector T gives unbalanced

torque. Divide by sling inertia for

acceleration about point P.

[17] w2 =. w2 +dT*wp2 NB. angular velocity of sling

[18] A2 =. A2 -dT*w2 NB. new sling angle

[19] S =. P + T=. SL*(1 _1)*1 2 o A2 NB. new T & projectile coord.

[20] NB. store new values in result matrix TB

[21] TB =. TB,[1] (N*dT),P,S,A,w,wp,A2,w2,wp2,V,a,T

[22] goto Ln0 if. (160 > N=. N+1) NB. increment counter and repeat

calculations if limit is not

exceeded.

[23] TB[;6 9] =. TB[;6 9]* 180% o.1 NB. I like angles in degrees

I welcome coments and suggestions toward improving ZOT